MODELS FOR TENT CATERPILLAR-VIRUS INTERACTIONS by Sarah Jenelle Beukema B.A. Kalamazoo College, 1988 A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in The Faculty of Graduate Studies Department of Zoology and Institute of Applied Math We accept this thesis as conforming to the required standard The University of British Columbia April 1992 © Sarah Jenelle Beukema, 1992 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of ^ The University of British Columbia Vancouver, Canada Date DE-6 (2/88) ^ ABSTRACT Many species of forest Lepidoptera show eight to twelve year population cycles which may involve viral disease. To examine possible interactions between viral disease and population cycles of forest Lepidoptera I explored some models for insect-virus dynamics. All of the models produced population oscillations in their original form. However, after they were modified to conform more closely to the tent caterpillar system, none of the models produced realistic cycles. I then developed a new model specifically for the tent caterpillar system that included important features such as: reduced fecundity of individuals that had been exposed to the disease, transmission of the disease from mother to progeny, a free-living infective stage of the virus, and a horizontal transmission rate that varied with the larval stage, the number of individuals, and amount of virus present. Cyclic dynamics resulted from some simulations. The parameters producing the cycles were similar to actual data. However, unlike natural populations of tent caterpillars, in the simulated population the average fecundity decreased before the population started to decline and survival decreased at approximately the same time as the population. Important further research in the field should include investigation of the distribution and survival of free-living virus and factors that would reduce caterpillar survival at peak populations but not affect fecundity. TABLE OF CONTENTS ABSTRACT ^ ii TABLE OF CONTENTS ^ iii TABLES ^ FIGURES ^ vi ACKNOWLEDGEMENTS ^ vii 1. INTRODUCTION ^ 1.1 Overview ^ 1.2 Organization ^ 1 1 3 2. BIOLOGY OF TENT CATERPILLARS ^ 4 3. MODELS ^ 3.1 What is a model? ^ 3.2 Modelling Background ^ 3.2.1 Exponential and logistic growth ^ 3.2.2 Lotka & Volterra ^ 3.2.3 Kermack McKendrick ^ 3.2.4 Nicholson & Bailey and beyond ^ 3.3 Models that produce cycles ^ 3.3.1 Predator-prey ^ 3.3.2 Weather ^ 3.3.3 Plant quality ^ 3.3.4 Parasites ^ 3.3.5 Diseases: systems with immunity ^ 3.3.6 Diseases: systems with no immunity ^ 3.3.7 Disease and genetics ^ 8 8 9 9 11 14 16 17 18 20 21 23 23 24 25 4. INSECT-VIRUS MODELS ^ 4.1 Introduction ^ 4.2 Anderson and May ^ 4.3 Brown ^ 4.4 Regniere ^ 4.4.1 Original model ^ 4.4.2 Modifications ^ 4.5 Diekmann and Kretzchmar ^ 27 27 28 32 36 36 39 56 5. A NEW MODEL ^ 5.1 Introduction ^ 5.2 Formulation of the model ^ 5.2.1 Model overview ^ 5.2.2 Larval stage ^ 5.2.3 Pupal stage ^ 5.2.4 Adult stage ^ 5.3 Parameter values ^ 60 60 61 61 64 67 68 72 nl 6. SIMULATION RESULTS AND ANALYSIS ^ 77 6.1 Introduction ^ 77 6.2 Simulation results ^ 78 6.3 Analysis ^ 90 6.3.1 Introduction ^ 90 6.3.2 Virus reproduction rate ^ 91 6.3.3 Stability^ 92 6.3.4 Comparison with simulation results ^ 95 6.4 Summary ^ 97 7. CONCLUSIONS ^ 99 8. LITERATURE CITED ^ 106 Al. APPENDIX 1: ANALYSIS OF MODELS ^ 110 A1.1 Introduction ^ 110 A1.2 Lotka-Volterra ^ 110 A1.3 Nicholson-Bailey ^ 111 A2. APPENDIX 2: ANALYSIS OF THE NEW MODEL ^ 113 A2.1 Introduction ^ 113 A2.2 Virus reproductive rate ^ 113 A2.3 Equilibrium analysis ^ 114 A2.4 Stability: numerical analysis ^ 118 iv TABLES Table 1. Effect of increasing parameters in the modified R,egniere models. . . . 43 Table 2. Summary of the different forms of p t used in the modified Regniere model. ^ 49 Table 3. Parameter estimations for b and c based on data^ 54 Table 4. Variables and parameters used in the new model ^ 63 Table 5. Parameter values and effect of increasing them in the new model. .^79 v ^ FIGURES Figure 1. The basic life cycle of the tent caterpillar ^ 5 Figure 2. Effect of increasing r on the behaviour of the discrete logistic equation. 12 Figure 3. Predator-prey cycles from the Lotka-Volterra model. ^ 19 Figure 4. Output from Brown's original model ^ 34 Figure 5. A schematic diagram of Regniere's model ^ 37 Figure 6. Existence and period of cycles in the modified Regniere model. . . . ^ 41 Figure 7. Effect of varying net fecundity of virus-exposed individuals on equations 4.7. ^44 Figure 8. Existence and period of cycles in the modified Regniere model which includes free-living virus. 46 Figure 9. A comparison of the different equations for the infection probability that were used in equations 4.11. ^ 50 Figure 10. A comparison of the different forms between birth terms and for the infection rate 0 59 Figure 11. Lifecycle of the tent caterpillar with the corresponding parameters and variables from the model. ^ 62 Figure 12. Existence and period of cycles in the model as the parameters /3, cr y , and f e are varied.^ L 81 Figure 13. Maximum and minimum clutch sizes produced during each of the cycles indicated in Figure 14. 84 Figure 14. Log number of tents and average clutch size in two actual populations and a simulated population. ^ 85 Figure 15. Early survival rates (from egg to fourth instar) in a simulated and a natural population 87 Figure 16. A comparison between various forms for to t , the amount of virus ingested at each instar. ^ 88 Figure 17. Constant p vs equilibrium number of susceptibles ^ Figure 18. Parameters x and z from the analysis vs constant p^ vi 94 117 ACKNOWLEDGEMENTS I would like to thank Don Ludwig and Judy Myers for always being there to guide me and for their helpful comments at all stages of the work. I also deeply appreciate the encouragement and support from Gordon Davidson and Andrew Trites. vii 1. INTRODUCTION 1.1 Overview Many forest lepidopteran species show eight to twelve year population cycles (Myers, 1988). At peak densities, the larvae often defoliate trees and are considered to be severe pests. Tent caterpillars are a typical cyclic species. Viral disease is usually associated with decline populations of cyclic caterpillars and may be involved in the population dynamics Numerous models for virus-insect interactions have been used to explore possible dynamics between virus and its hosts (Anderson and May, 1980; Brown, 1987; Regnieie, 1984; Diekmann and Kretzchmar, 1991). To investigate the applicability of these models to a realistic field situation, I adapted several of them to details of the tent caterpillar system. All of the models produced population oscillations (usually in the form of limit cycles) in their original form. However, after they were modified to conform more closely to the tent caterpillar system by adding vertical transmission of the virus, lower fecundity of exposed insects and a free-living virus stage, and simulated using realistic parameter values, none of the models produced reasonable cycles. Even though no previous insect-virus model produces reasonable cycles, it is still possible that virus creates the natural population cycles. To examine this, I created a new model specifically for the tent caterpillar system. It included various important features such as: reduced fecundity of individuals that had been exposed to the disease, transmission of the disease from mother to progeny, free-living infective stage of the virus, and a horizontal transmission rate that varied throughout the larval stage and that depended on the number of individuals and amount of virus present. 1 The model was simulated using a wide variety of parameter values and analysis was done on a simplified version of the model. Although the analysis showed that the simplified model would not produce cycles, eight to twelve year oscillations did occur with some parameter values in the full model. However, in all cycles that are produced, the patterns shown in the oscillations are different than those seen in natural populations. Natural populations show that the average clutch size remains high even in the first year of a population decline while survival rates decrease before the population decline. In the model, the clutch size is reduced before the population decreases, and is the cause of the population decline. Survival rates decrease at the same time as the population. Thus, on the basis of current belief about the interactions between virus and insects and about the biology of the virus, it is unlikely that virus is the sole cause of cycles. However, it is possible that virus combines with some other factor to create the cycles. 2 1.2 Organization The thesis is organized into several parts. In the Chapter 2 I describe the biology of tent caterpillars. Chapter 3 contains background information about models including a definition of a model, a description of how models can be used in biology, a brief account of some of the important models in biology and a discussion of various hypotheses and models that produce cycles. In Chapter 4 insect-virus models are examined in more detail, and various modifications to the models are described. Finally, a new model created for tent caterpillar-virus interactions is formulated in Chapter 5 and the results of simulations and analysis are discussed in Chapter 6. All the actual analyses of the various models are in Appendices Al and A2. 3 2. BIOLOGY OF TENT CATERPILLARS At least sixteen species of forest lepidoptera in the North Temperate zone show population cycles. Some well studied examples are the spruce budworm, the tussock moth and the larch budworm, which reach outbreak levels every eight to twelve years and periodically are major defoliators of trees (Myers, 1988). Tent caterpillars, Malacosoma spp, also show cyclic dynamics and the western and forest tent caterpillars have been extensively studied (Daniel, 1990; Myers, 1990). Tent caterpillars are univoltine. The adult female lays her eggs in a single egg mass in July and then dies. Neonate caterpillars overwinter in the eggs and hatch early the following spring. Caterpillars from a single egg mass remain together and form silk tents in trees. The size of the tents is determined by the age and number of the caterpillars in the group. As a consequence, the size and number of tents are good indicators of population size (Myers, 1990). Larvae feed gregariously until the fifth instar when they may start wandering individually and may move more than thirteen meters from their tent before pupation 50 to 80 days after hatching (Wellington et al., 1975). Adult moths mate shortly after emergence and females lay their eggs soon after mating. The basic life cycle is shown in Figure 1. Many agents kill caterpillars, including a variety of parasitiods and predatory insects (Smith and Goyer, 1986). Caterpillars can also be infected by a nuclear polyhedrosis virus, bacterial and fungal diseases. Much work has been done on different nuclear polyhedrosis viruses of Lepidoptera since they are usually species specific and have great potential as biological control agents (Entwistle, 1983a). 4 Eggs overwinter (July - April) reproduce, die Adults Virus ^Larvae (Instars 1-5) (April - June) (July ) Pupae (June - July) Figure 1. The basic life cycle of the tent caterpillar. 5 Nuclear polyhedral virus (NPV) can be transmitted in two ways. Most common is transmission within a generation (horizontal transmission). As early as eight days after infection, the infected larva dies (Clark, 1958). The corpse disintegrates rapidly and the virus propagules (polyherdral inclusion bodies, PIBS) are released onto the surrounding leaf. Other caterpillars become infected by ingesting the virus while feeding (Clark, 1955; Entwistle, 1983a). Caterpillars can also spread the disease by moving over the contaminated leaves (Stairs, 1965). The larger (older) the caterpillar, the more viral propagules are produced when the organism dies. However, larger larvae are also less susceptible to the disease, so that more PIBS must be ingested to initiate infection (Entwistle, 1983b). NPV is also apparently transmitted from a mother to her progeny (vertical transmission) (Clark, 1955, 1958; Myers, 1990; Wellington, 1962), most likely on the coating surrounding the egg. A portion of the caterpillars will become infected at hatching if they ingest viral polyhedra. If eggs laid by females from a virus infested area are brought into a disease-free laborabory or field site, the larvae arising from these eggs will die from NPV infections (Clark, 1955, 1958). Whether a pool of virus exists on the food plants or in the soil is debated (Clark, 1956; Entwistle, 1983a). If such a reserve exists, it could also be a means of infection between generations. Clark (1956) showed that caterpillars arising from eggs from a disease-free area but that had been transferred to a site with a previous population of diseased individuals showed viral infection, indicating that there is a free-living virus pool and that the virus can survive over the winter. 6 Many features of the relationship between the larvae and the virus are not yet well understood. Can a sublethally infected caterpillar (or a caterpillar that is exposed to the virus during the late larval stage but does not die from the virus) act as a carrier for the virus and pass it on to her offspring? If so, what proportion of viral transmission occurs in this way? Laboratory experiments using another lepidopteran, Spodoptera ornithogalli (yellow striped armyworm), have shown that although exposure to the disease (through the larvae eating food infected with the virus PIB's) is not sufficient for vertical transmission, vertical transmission can occur via adult moths that are contaminated with virus. The larvae that emerge from the eggs laid by these infected moths can become infected (Young, 1990). Tent caterpillars from eggs collected from declining populations die from viral and fungal disease when reared in the laboratory (Myers, pers. comm. and 1990; Clark, 1955). Regardless of how the disease is transmitted vertically, some larvae are born infected (or become infected as they are emerging from the egg). This leads to further questions such as if an organism is born infected, will it survive to reproductive age and if it does, will it be able to reproduce? For the yellow striped armyworm mentioned above, the larvae can be exposed to the disease (from eating food with PIB's) and not die. However the adults may suffer reduced fecundity as a consequence (Young, 1990). It is possible the same holds for tent caterpillars. Even though these questions remain to be answered, effects of such assumptions can be explored using models. 7 3. MODELS 3.1 What is a Model? A model is a means of describing something using words, objects and/or equations. In the physical sciences, models are usually considered to be quite reliable since people are usually concerned with relatively deterministic and well understood events. However, in biology events are much more variable, and there tend to be fewer data to work from. Therefore, models in biology are often more speculative, but can be used to clarify ideas or to explore the predictions of a hypothesis (Starfield and Bleloch, 1986). Models can also suggest what data should be collected or what experiments should be performed in the field to test hypotheses. A 'working' model is often considered to be one in which there is agreement between observations and the model output. It is important to realize that in biology, a 'working' model does not necessarily mimic the way nature operates: often there are many different models that give the same results. For example, as will be seen, cycles can arise from a wide variety of different hypotheses. If one model produces a desired result, all that can be concluded is that the hypothesis that the model was based on is still feasible. If a model does not produce output which corresponds to natural occurances, it indicates that the hypothesized interactions are incomplete, the functions to describe these interactions are inappropriate or that an alternate hypothesis should be formulated. Consequently, more is often learned from a model that does not produce an expected outcome. The first section of this chapter describes some of the classic models in biology, namely the exponential and logistic models, the Lotka-Volterra models for predator8 ^ prey interactions, the Kermack-McKendrick epidemiology model and the NicholsonBailey model for host-parasitoid interactions. Each of these models is still used today as a basis for creating more detailed models. In the second section, I discuss existing hypotheses and models for population cycles. 3.2 Modelling Background 3.2.2 Exponential and Logistic Growth. While mathematical models have only recently been widely used in biology, the basic ideas of simple population growth extend back to the eighteenth century when interest in human demography led to calculations of the doubling time and carrying capacity for the human population (Kingsland, 1985; Hutchinson, 1978). Although he was not the first, T.R. Malthus is usually credited with the realization that populations will grow geometrically when unchecked in any way (Hutchinson, 1978; Malthus, 1798). The model for this basic population growth can be formulated by stating the tautology that the rate of change in the size of a population depends on the births minus the deaths: in mathematical terms ^ dN dt = bN — dN (3.1) where b and d are constant birth and death rates. If b and d are independent of N and t, equation 3.1 can be solved to give: ^N(t) = No e (b—d)t .^ (3.2) This exponential growth model is often labelled as the Malthusian model, and the parameter r = b— d is sometimes referred to as the Malthusian parameter (Hutchinson, 9 1978). Although the model was originally formulated as a model for human population growth, it is now generally used only in the study of population growth in bacteria, yeast and other unicellular organisms. No population can continue to grow indefinitely in exponential mode. Eventually some factor, such as lack of space or food, will limit this growth. In 1836, Verhulst first proposed the logistic model which includes these ideas: dN = r(1 - -k-)N dt (3.3) (Hutchinson, 1978). The birthrate is the term r(1-N/k) where k is the carrying capacity for the population (the maximum population size the environment can sustain, due to space or food or some other limiting resource) (Murray, 1989). This model was also derived independently by several researchers (including A.J. Lotka and R. Pearl) in the 1920's, and popularized at that time by Pearl (Kingsland, 1985). Consequently, the logistic model is sometimes referred to as the Pearl-Verhulst equation. Upon examination of the model, it is clear that if the population N is small, the birthrate term will be close to r and equation 3.3 will behave essentially like equation 3.1 and exhibit exponential growth. However, when the population is near k, the birthrate term will be very small. Equation 3.3 may be solved in the form No kert N(t) = k No(ert 1)• (3.4) It can be seen that as t^oo, N -p k. If the population is larger than k, the birthrate will be negative and the population will decline. So for any initial population size and parameter values, the population will always approach k as t^oo. 10 There is also a discrete version of the logistic equation: Nt T )Nt Nt+1 = Nt + r(1 — (3.5) (May, 1973). This model behaves very differently from the continuous model. Figure 2 shows that, as the value of r increases, the model behavior changes. For low values of r, the population will stabilize. If r is increased, the model will show cycles that decay down to some stable level. A larger r produces stable two-point cycles. The cycle period increases as r increases until eventually the population behaviour is chaotic (May, 1973). 3.2.3 Lotka & Volterra The next major breakthrough in modelling came almost simultaneously from a physical chemist, A.J. Lotka, and a mathematical physicist, V. Volterra in the 1920's (Kingsland, 1985). They each proposed a system of equations to model the interactions between two populations. Lotka derived the model as an extension to his work with single species populations. Volterra wanted to explain the oscillations in the catch levels of certain species of fish by using predator-prey interactions (Volterra, 1926). He used basic balance equations involving the rate of change of the prey population N and the predator population P, and made the following assumptions: 1. Each population is made up of organisms that are exactly the same as each other, (for example, in age and behaviour), and this composition of the population does not vary over time. 2. Prey would increase exponentially at a rate r without predators present and predators would decrease exponentially at a rate d with no prey present. 11 I^1 ^I^1^1^1^1^1^1^1 \ , \ \ \^I 1 \ 4 0 B) r=2.1 0 10^20 ^ 30 Time Figure 2. Effect of increasing r on the behaviour the discrete logistic model. The first graph show the population stabilizing. The next graphs show two-point cycles, four-point cycles, and finally, chaos. 12 3. The encounter rate between prey and predators is governed by mass action. Thus, predation reduces prey growth at a rate proportional to the predator and prey populations bN P (a different b from Section 3.2.2, equations 3.1 and 3.2). 4. Predator growth due to prey is also a rate proportional to the size of both populations cNP. 5. Predators and prey population sizes react instantaneously to changes. (Volterra, 1926; Kingsland, 1985; Murray, 1989). Thus the system of equations is: dN = rN — bNP (3.6a) dt dP = cNP — dP. (3.6b) dt These equations are usually called the Lotka-Volterra model. The analysis of this model is shown in Appendix 1.2, where it can be seen that the populations will always cycle. The model is not useful in modelling real populations since the predicted population behaviour (the period and amplitude of the cycles) is entirely dependent on the initial conditions and on the parameter values, none of which can ever be estimated exactly. However, the model is important as a basis for building other, perhaps more realistic, predator-prey models. 13 3.2.4 Kermack Si McKendrick Soon after Volterra proposed his model, two other classic models were created. The first of these was an epidemic model by Kermack and McKendrick (1927). At that time, there was much controversy about what caused the termination of an epidemic. The two theories that were popular were that: a) there were no more susceptible people and b) the virulence of the disease had decreased during the epidemic (Kermack and McKendrick, 1927). Kermack and McKendrick decided to use a simple model to examine the question. They looked at the effects of disease in a closed population (i.e., no births, deaths or other movement in or out of the population), in which the population was divided into three classes: susceptibles S, infecteds I and removed R. Contact with an infected individual is necessary for transmission of the disease to occur, so the susceptibles become infected at a rate which is proportional to the number of infecteds and susceptibles individuals present in the population (0/5). Thus, dS OIS. at (3.7a) The infected individuals are removed from the population (either by immunity or death) at a rate v. So the equations for the final two classes are dI dt OIS — vI dR = vI. di (3.7b) (3.7c) Through analysis of this model, Kermack and McKendrick determined that the epidemic dies out before the susceptible population is exhausted. There is some population 14 threshold level that can be easily derived. First, rewrite equation 3.7b as dI Tit = (s - ')i (3.8) Let V NT=. (3.9) Assume that the initial conditions for the model state that there is a large amount of susceptibles (So ) present, and few infecteds (Is ). Then, according to equation 3.8, the rate of increase of the infecteds will only be positive if S o > NT (Anderson, 1991). The rate of increase of infecteds must be positive for an epidemic to occur. If the population is smaller than NT, the infecteds will eventually be removed, and an epidemic cannot take place. This result is referred to as the "threshold theorem". This theorem and this model are the foundations for much of subsequent epidemiology modelling (Anderson, 1991). 15 3.2.5 Nicholson & Bailey and beyond The entomologist, A.J. Nicholson and physicist, V.A. Bailey proposed the second of the classic models in 1935. Nicholson originally wanted to examine possible hypotheses of population regulation, and started with host-parasite interactions (Kingsland, 1985). Bailey assisted him with the mathematics. Initially, they examined host-parasitoid interactions in which the parasitoid kills the host. In their model, they assumed that only the fraction of unparasitized hosts can reproduce. The number of parasites (P) depends on reproductive rate of the parasitiods (c) in the parasitized hosts (N(1 — f)). Assuming random encounters, the fraction that escapes parasitism (f) can be expressed using a Poisson distribution (Edelstein-Keshet, 1986), ^f = e—aPt ^ (3.10) and the final equations are: Nt+1^AN t a Pt ^ Pt+1 = CAri (1 —^e—aPt ^ (3.11a) (3.11b) This model is analyzed in Appendix 1.3, and is shown to produce growing, unstable, oscillations. This model has been the basis for a large number of other parasite models (see Edelstein-Keshet (1986) for a brief description of some of the model's modifications). In the late 1960's and early 1970's modelling in population ecology became quite popular. There were several books written at that time (Pielou, 1969; May, 1973) which discussed several of the basic models, their properties, applications, and some modifications. Since then, modelling has become increasingly popular, both among 16 mathematicians and biologists. As computers become faster and more powerful, models are becoming more complex and analytically intractable. But the basic purpose remains the same: to use models as an effort to clarify and to understand what is happening in nature. 3.3 Models That Produce Cycles There are many hypotheses about potential causes of population cycles. Some of these include predator-prey interactions, weather patterns, plant quality, parasites and diseases. In all cases, the hypotheses contain the assumption that it is just one major factor that is driving the cycles. Finerty (1980) discusses the existence of population cycles in small mammals and the relevance of many of these hypotheses to the cycles. Models, both theoretical and simulation, have been created for each of these hypotheses, and many show cycles. In this section, I discuss existing models for population cycles of organisms in general, not just those of insects, and mention some of the different hypotheses for population cycles and some of the major models created to support these hypotheses. The relevance of some of these hypotheses to the tent caterpillar system is also discussed. 17 3.3.1 Predator-prey The predator-prey hypothesis states that as the prey numbers increase, the predator numbers will also increase until such a point when the predators actually reduce the growth rate and number of prey. This causes the predator numbers to decrease from lack of food. When the predators become few enough, the prey population can start to grow again. Figure 3 gives an example of predator-prey cycles. The classic Lotka-Volterra model (see equations 3.4) for predator-prey interactions show these cycles. However these are neutral cycles, meaning that the period and amplitude of the cycle depend on the parameter values and the size of the starting population. The model is structurally unstable since small changes in parameters give different results. If a researcher were trying to use the Lotka-Volterra model to predict behaviour of the population, and if the parameter estimates were even slightly different from the true parameter values, the predicted period and amplitude of the cycle would be different. Thus the model is impractical, since parameters such as reproductive and mortality rates vary from year to year and population to population. This model is useful as a starting point for creating other predator-prey models (or plant-herbivore models). However, predation does not seem to be a major factor in the tent caterpillar system and it is more likely that something else is driving the cycles. 18 1 3 i 3 i 20 Generalions Figure 3. Predator-prey cycles 19 I^I^I 40 30^ 3.3.2 Weather Weather variation has been hypothesized as a control of insect outbreaks (Martinat, 1987). Bad weather will tend to keep insect populations at a low level (Andrewartha and Birch, 1954). If there is a period of weather that is beneficial to the insects' survival or fecundity for several generations, then the insects may be able to increase rapidly (Andrewartha and Birch, 1954). Weather can have both direct and indirect effects on insects (Martinat, 1987). Direct effects act on the insect without any intermediary agent and include the effects of weather on the behaviour or physiology of the insect (Martinat, 1987). For example, the development rates of many insects are connected to the temperature. Indirect effects are those in which the weather affects food sources or predators (including parasites) of the insect. Such a case could be that the effect of bad weather may cause a plant to become stressed and to produce fewer defenses or to alter in nutritional value. If the plants had fewer defenses, insects would be able to eat more of them, and the insects' survival would increase (Martinat, 1987). Weather patterns were included as a critical factor in the simulation models of tent caterpillars by Wellington et al. (1975; Thompson et al., 1979). Growth, mortality, disease, and dispersal were all simlated to depend on the weather patterns of the area. They also included other aspects of the biology such as parasitism and the condition of the larvae. The weather patterns that they used initially in the model were all based on available data about the actual weather in the study area. In their model, they were able to produce population trends that corresponded to their data. It seems possible that the weather was driving the behaviour of the system, especially since when they 20 altered the weather patterns in the model, some of the results changed and no longer matched the behaviour of the natural system. In other tent caterpillar systems however, the cycles do not seem to be related to the weather (Daniel, 1990; Myers, 1981). After examining maps and data, Daniel (1990) concluded that larval feeding temperature and overwintering temperature do not explain the dynamics of the forest tent caterpillar in Ontario, although the temperatures may determine the long-term susceptibility of a particular area to population outbreaks. Myers (1981) also noted that since in her populations, rainfall often can be correlated with the increase phase of the cycle but not the decrease phase, and temperature patterns cannot be correlated with any part of the cycle, weather is not the major factor determining the regular population oscillations. 3.3.3 Plant quality Many plants have some form of response, chemical or otherwise, that is induced when they are attacked. When herbivores are at high population densities and the plants are severely stressed by the herbivores' feeding, the plants can induce a response which will increase as the rate of feeding increases (Rhoades, 1983). In some cases, these responses can be detrimental to the herbivore and can cause reduced fecundity, death or generally affect herbivore population growth in some other manner. Thus the herbivore population would decrease and, after some period of low herbivory, the plant quality would increase again, creating a delay in the system which could lead to herbivore population cycles (Haukioja, 1980). Two models for plant-herbivore interactions in which plant quality is important 21 have been formulated. A theoretical model demonstrating that in plant-herbivore interactions the plant's quality could be a governing factor for cycles, was created by Edelstein-Keshet and R,ausher (1989). In the model, cycles occur for a limited set of circumstances. A detailed model incorporating food quality in the larch budmoth Zeiro- phera diniana system produced cycles that corresponded to the data for that system (Baltensweiler and Fischlin, 1988). However, the plant-quality hypothesis is probably not applicable to the tent caterpillar system. In this system the quality of the trees that the caterpillars eat varies, but it is unlikely that it has major impact. Work by Myers (1981; 1988a) has shown that, even when plants are protected from severe feeding stress or are subject to simulated drought conditions (both of which should cause the plants to have a higher quality), caterpillars do not survive better. If plant quality were the governing factor for cycles, the caterpillars should have a higher survival rate on plants with a higher quality. 22 3.3.4 Parasites Many models have explored the interactions between parasites and hosts. These include the Nicholson-Bailey parasite model (see Section 2.5) and its many modifications. Anderson and May (1978; May and Anderson, 1978) created many models exploring different relationships in host-parasite interactions. They found that the relationships that produced cycles were a parasite-induced reduction in the host's reproduction or time delays in the development or transmission of the parasite's infective stage (May and Anderson, 1978). Time-lags will usually generate cycles in a model (May, 1981), and are present in all systems since organisms cannot respond instantaneously to changes in nature. O 3.3.5 Diseases: systems with immunity Considerable theoretical work has been done on the effect of disease on animal and insect populations. The models fall into two main categories, those in which the host can recover and sometimes become immune to the pathogen, and those in which there is no recovery or immunity. The first type of model is used primarily for human diseases but could easily be adapted for any organism with an immune system. Bailey (1975) reviews many of these models, from the earliest epidemiology models of Kermack and McKendrick to those published in the early 1970's. More recent work specializing in nonlinear contact rates includes Hethcote and van den Driessche (1991), and Liu et al., (1986; 1987). In general, these models assume a large, closed population, for which birth and death are ignored and most of the differences between models arises from changing conact rates or assumptions about immunity and recovery. There is also a 23 smaller class of models that incorporates a changing population size (Anderson, 1979; May, 1979). Both types of models (those with and without a varying population size) assume recovery or immunity. Some of these models have been reasonably successful at capturing the cyclic dynamics of human measles and flu outbreaks (Bailey, 1975; Hassell and May, 1989). 3.3.6 Diseases: systems with no immunity Models of systems without immunity from the disease are of particular interest when studying insects, since although populations of insects can increase their resistance to a virus (Briese and Podgwaite, 1985), the insects have no apparent immune system. These models assume open populations and include births and deaths. The disease models give an indication of the many different situations in which the pathogen can regulate the host, what factors are necessary for disease persistence, and when the system is unstable (Hassell and May, 1989; Anderson and May, 1981; Anderson and May, 1980, 1981; review by Onstad and Carruthers, 1990). For cycles to exist, the equilibria must be unstable. Many disease-insect models have shown that imperfect vertical transmission of disease will give population cycles of hosts (Anderson and May, 1980; Brown, 1987; Regniek, 1984; Diekmann and Kretzchmar, 1991; Hochberg, 1989). A more detailed discussion of some of these disease-insect models occurs in Chapter 4. 24 3.3.7 Disease and genetics The effect of disease on the genetics of a population is important in another hypothesis, formulated to explain the cycles in forest lepidoptera, and incorporating the variability in individual resistance to the virus. This hypothesis proposes a trade-off between two different types of organisms, ones that have a high fecundity but are very susceptible to the virus, and ones that are much more resistant to the virus but have a low fecundity. The susceptible individuals with high fecundity are selected for when the population is near low density, so it is primarily these susceptibles that cause the population increase. However, when the population is high, more larvae become infected with virus and it is the resistant type that survives the viral epizootic, and the population declines (Myers, 1988b, 1990). This hypothesis involving the genetics of the organism is very similar to the Chitty Hypothesis for small mammals, except that the Chitty Hypothesis proposes a trade-off between aggression and fertility: at high population densities, there are more organisms that are highly aggressive but have low fecundity so the population declines and those animals with high fecundity that are not as aggressive will be selected for and will eventually cause the population to increase again. This hypothesis appears to be likely to create population cycles. However, the Chitty Hypothesis has been extensively modelled and most models do not produce population cycles. Those that do generate oscillations are often not biologically reasonable or only cycle for a very narrow range of parameter values (Stenseth, 1981). One possible exception is the model by Hunt (1982, 1983). She shows mathematically that cycles are possible. However, it is difficult to interpret her 25 results biologically and to reproduce them through simulation. A model of genetically controlled disease resistance and variation in fecundity for the tent caterpillars is similar to the models for the Chitty Hypothesis, and it does not produce cycles. 26 4. INSECT-VIRUS MODELS 4.1 Introduction There are many interactions that can give cyclic dynamics in models: predator-prey interactions, weather-insect interactions, host-parasite interactions, plant-herbivore interactions and host-disease interactions. All the models referred to in the previous chapter, except the one by Wellington et al. (1975), are very general and do not necessarily consider the specific biology of any particular system. Indeed, as they have been formulated, they could be used to explain the population cycles either of insects or of small mammals (except for differences in immunity to disease). Obviously, the biologies of insects and mammals are different in many ways. It would therefore be interesting to examine one hypothesis, or one set of models, in more detail and to apply it to a system and see if the cycles are robust. Here, disease models for insects will be considered. Four existing models will be examined and adapted to to details of tent caterpillar populations. 27 4.2 Anderson and May Anderson and May (1980) modelled host-disease interactions using differential equations for continuous overlapping generations of classes of infective (I) and susceptible (S) insects, and a free-living infective stage of the pathogen (V). The model was formulated in three equations. The first is for the susceptibles (S). The rate of change of the susceptible population depends on the birthrate (a) of all individuals, both susceptible and infected, minus the natural mortality (d) of the susceptibles minus the newly infected individuals. dS — a(S + I) — dS — 13V S dt (4.1a) The term 13V S is the infection rate of a susceptible individual and is proportional to the rate of encounters (3) between the susceptible and the free-living stage of the pathogen. The change in the infected population (I) results from the addition of the newly infected individuals minus dead infected individuals (natural mortality rate (b) + disease induced mortality (a)). So, dI dt- = i3VS — (a + d)I. (4.1b) The third equation is for the free-living stages of the virus (V). These infective stages are produced from infected hosts at a rate A and die at a rate p,. They are also lost through ingestion into susceptible and infected hosts. This gives, dV . AI — (p, + 0(S + I))V. dt (4.1c) Through analysis of this model the authors show that it is possible to produce limit cycle behaviour from the model. The crucial factor in the model is the mortality rate 28 of the free-living virus (u). If there were no free-living stage of the virus (so susceptible insects became infected directly from other infected insects), the model does not produce population cycles. The authors compared the results from this model to the behavior of the larch budmoth, a lepidopteran that has cyclic dynamics. Given parameter values based on data, the model produced ten-year cycles, the same period as the larch budmoth population shows. However, although the range in the predicted total population was similar to the actual population size (the log of the population size ranged from about 3.5 to 6.5 while the log of the actual population ranged from approximately 2 to 7), the predicted and actual prevalence of infection was quite different (between 30% and 50% in the actual population compared with 80% in the model). Vezina and Peterman (1985) also used the Anderson and May model (equations 4.1) to model the Douglas-fir tussock moth ( Orgyia pseudotsugata), another lepidopteran species that has 7— 10 year population cycles. They simulated this model and several variations of the model that included more of the specific biology of the moth. The first modification was to make the mortality rates density-dependent by changing the minimum natural death rate d to d' where d' = d c(S + I) (c is some measure of the severity of density-dependent natural mortality (Vezina and Peterman, 1985)). The disease-induced mortality a becomes a = m — d' where m is total host mortality. The second modification was to add an incubation period between the time an insect became diseased and when it infected others. A new class of organisms that are infected but not yet infectious (E) was made. Organisms move into that class when they become 29 infected and move out at a rate v. The new system of equations is: dS dt dE dt dI a(S + I) – dS – 131/ S^ (4.2a) ^– (d v)E^ (4.2b) vE – (a + d)I^ (4.2c) dV = AI – (ft + 13(S + .0)V.^ (4.2d) dt The third model Vezina and Peterman (1985) examined included vertical transmission. This means that not all the new larvae are susceptible: a proportion k of them are born infected (and a proportion 1 – k of them are born susceptible). Thus, equations 4.1 become dS dt aS a(1 – k,)I – dS – S dI — = pvs- (a ± d)1" akI dt dV =^– (,u, 13(S + I))V. dt (4.3a) (4.3b) (4.3c) Finally, they looked at a model in which all these components were combined. dS – aS all – k)./- – dS – /3V 5' dt dE dt /.3V S – (d' v)E akE dI dt vE – (a + d')I dt AI – (ft + 13(S + .0)V. (4.4a) (4.4b) (4.4c) (4.4d) Venzina and Peterman (1985) obtained realistic parameter values for the models from field studies reported in the literature and adapted the values to fit into a 30 continuous-time model (so that all values were calculated per year and there would be discrete generations). They found that it was impossible to obtain cyclic dynamics with appropriate periodicity when realistic parameters were used. They concluded that it was not virus alone that generated population cycles so the models were not an adequate representation of the tussock moth system. They hypothesized that a more successful model would include the interactions between the tussuck moth and its enemies and food. One problem with the models of Anderson and May (1980) and Venzina and Peterman (1985) is that they assume that the transmission of the virus is proportional to the initial densities of susceptible and infected individuals or to the densities of susceptibles and free-living virus. This assumes that each instar is as infectious and resistant as all other instars. Dwyer (1991) performed some experiments in the field with the Douglasfir tussock moth to determine the transmission coefficient found that (13) at different instars. He /3 was significantly higher for older instar larvae, indicating that transmis- sion over the entire year is not strictly proportional to the intial numbers of first instar larvae. It is possible that if the models were adapted to account for a changing would occur using realistic parameter values. 31 3, cycles / 4.3 Brown Brown (1987) created a general simulation model for insect-virus populations. The model was quite detailed and included a five day latent period (between becoming infected and dying and infecting others), an environmental pool of free-living stages of the virus, and a random temperature component that varied insect growth rates. Also, the amount of virus infected larvae produce when they die increased with the age of the larvae. In the model, the probability of becoming infected at any time depends on what Brown called the pathogen burden (`pb') and on the density of the hosts (`de'). The pathogen burden is the number of infectious units that are available to infect hosts per host, or, pb th (4.5) where 'pa' is the number of available pathogens and `th' is the total number of hosts. Infection is considered to be a random process so a Poisson distribution is used, and the equation for the probability of infection is p = 1 - e—de*pb (4.6) Notice that if the density (`de') is large, p will be close to one. Biologically, this simulates greater stress and increased susceptiblity at high densities (Brown, 1987). However, p is independent of the age of the larvae, and thus does not account for any increasing resistance to the disease. The model simulates overlapping generations (so that on any given day there are insects of all age classes) and adult insects lay several eggs on each of several days. Both susceptible and exposed adults lay eggs, and they have the same fecundity. Even 32 though this simulation model shows oscillations, the cycles are unstable, with increasing amplitudes (Figure 4). I adapted this model to the tent caterpillars system by making the following changes. 1. Exposed individuals have a lower fecundity than non-infected individuals. The original model assumed that exposed and susceptible adults would have the same fecundity, represented by the variable `ov'. I added a new variable `ove', to represent the fecundity of the exposed individuals, with `ove' < `ov'. 2. The disease is carried by a proportion (q) of the eggs laid by the exposed adults (i.e., the larvae will be in the exposed class when they emerge). The rest of the eggs will be susceptible. Thus the equation governing the reproduction, S1 = ov * (s7 + E. ) ^ (4.7) becomes Segg OV * S7 + ove * Ea * (1 — q)^(4.8a) Eegg = ove * Ea * q.^ (4.8b) where: Si is the number of susceptible eggs in the original model, S7 , Ea are the susceptible and exposed adults, Segg Eegg are the susceptible and exposed eggs in the modified model. 3. Reproduction only occurs towards the end of the season, with the adults laying all of their eggs in a single batch and then dying. The original model had the adults 33 i i 1 1 0 500 1000 1500 Days Figure 4. Output from Brown's original model, showing unstable oscillations. 34 laying a few eggs each day of their adult life. Thus, any time an organism becomes an adult, she would lay a batch of eggs on that day and then die. 4. Only the eggs and the free-living virus survive between seasons. I decided that the maximum total life expectency (from egg to adult) for a tent caterpillar was 100 days. This means that at the end of the season, (at the end of day 100), all insects in all classes except Segg and Eegg , are killed. This modified model no longer produced cycles. The simulated population stabilized or grew exponentially, depending on the birth rates that were used and the reduction in fecundity attributed to the exposed insects. If there only a small reduction in fecundity due to infection, the population would be able to grow exponentially due to high fecundity. When there was a larger sublethal effect of virus on fecundity, the population stabilized. 35 4.4 Regniere 4.4.1 Original model R,egniere (1984) produced a model which most closely represents forest Lepidoptera. He used difference equations which included vertical transmission of the virus (k), no recovery from infection, and lower fecundity of diseased organisms due to sublethal effects of the virus (Figure 5). Difference equations are more appropriate than differential equations since the insects only reproduce once each year. In the model, the number of susceptible individuals at the end of the generation can be calculated by examining three categories of organisms: 1. The number of surviving susceptible individuals that did not contact the disease and gave birth to susceptible offspring ((1 — p t )cSt ). Since the function pi describes the proportion of insects that become infected (see equation 4.10), the term 1 —p t gives the proportion of insects that do not become infected. 2. Those susceptible organisms that contacted the disease but did not pass it on to their offspring ((1 — k)bp t St ). This implies that either these organisms were carriers for the virus but they did not transmit it to their offspring, or that these were organisms that were more resistant to the disease. 3. The surviving infected individuals that did not pass the infection to their young ((1 — k)alt ). Adding items 1, 2, and 3 leads to the equation: St+i^(1 — k)(aIt + bpt St) + c(1 36 — pt)S t •^(4.9a) f s a s (1-p)S (1-kie ae p S (1-k)f.ai I i disease free k^PS feae kf. a. I i^1 diseased Figure 5. A schematic diagram of Ileguiere's model (1984) (based on the diagram in his paper). In the model, a b fe a,, and c La s . The top of the diagram represents generation t while the bottom represents generation t 1. 37 The parameters a, b and c are the net fecundities of infected, exposed and susceptible individuals respectively, and incorporate insect survival, fecundity, and the proportion of eggs that hatch. The net fecundities are related to each other by a < b < c, assuming that infected organisms are adversely affected by the disease more than exposed or susceptible insects. Regniere assumed that a > 0, so that infected indivivals could reproduce. For the population to survive, c must be greater than 1. Vertical transmission is included in the model in the parameter k, the proportion of offspring born to exposed or infected mothers, that are infected (k < 1). Similarly, the number of new infected individuals is calculated by the number of susceptibles that contacted the disease and gave birth to infected offspring (kbptSt), plus the number of infected individuals that also passed the infection to their offspring (ka/t ). = k(a/t bptSt) (4.9b) The proportion of insects becoming infected, Pt , is based on assumptions of random movements of individuals and increases with the number of infecteds present. Pt = (1 — e —h-rt)# (4.10) where h is the rate at which the healthy larvae contact disease propagules (Regniere, 1984). The parameter /3 determines the aggregation of disease propagules near the diseased insects, where 13 > 1. (Low near insects while high /3 implies that the disease is highly concentrated /3 indicates that the disease is more widely dispersed (Regniere, 1984).) The notation is slightly changed from the original paper: Regniere used x and y where I have used I and S for clarity. 38 This model produces cycles with the eight to twelve year period similar to those observed in natural populations. However, there are aspects of the biology of forest Lepidoptera that are not yet incorporated. I modified the model to add some of these aspects. In general, as the model became more detailed, it lost any tendency to cycle. 4.4.2 Modifications Organisms that become infected at hatching probably die before reproducing (since infected larvae usually live for about ten days (Entwistle, 1983a)). Therefore, I modified the original model to include this fact: kbp t St St+i = (1 — Obpt St + c(1 — pt)S (4.11a) (4.11b) As shown in Figure 6, it was still possible to get ten to twelve year cycles in the population, but with a narrow range of parameter values. Specifically, stable cycles (of any length) primarily occur only when /3 is close to one (the minimum value of 0 as indicated by Regniere (1984)), meaning that the disease is highly concentrated near insects and the insects have a greater probability of contacting the disease. Biologically this would be the same as increaseing the amount of virus infected larvae produced when they died, or reducing the amount of virus needed for infection of a larva. When /3 is larger than one, (implying that the disease is more scattered), there is a much narrower range of parameters that will cause cycles. In many instances, the infected organisms quickly die out and only the susceptible ones are left. Because the disease is more widely dispersed, the chance of becoming infected is very low. Consequently, fewer infections occur and the virus disappears. 39 Figure 6 demonstrates the effect of changing three of the parameters. In all cases, c, the net fecundity of susceptible organisms, is held constant, and is equal to the upper value of b, the net fecundity of exposed individuals. If c were to change, the graphs would also change, although the general pattern of the regions with cycles would remain the same. Figure 6 also shows that it is not necessary to have any sublethal effect of the virus (reduced survival or fecundity, characterised by b < c), for oscillations to exist. In Figure 6, a distinction is made between cycles and persistent oscillations (for mathematical, not biological reasons). Cycles are those in which each oscillation is identical with all others. Persistent oscillations are those in which the amplitude of each cycle varies in a consistent, repeatable manner (for example, every third oscillation is identical), and the periodicity is relatively consistant. These persistent oscillations are stable to small perturbations (i.e., if the population is perturbed, it will quickly return to the same periodicity and amplitude as before the perturbation). Either of these types of cycles could be relevant to natural populations, especially since in natural populations each oscillation is not exactly identical to all the others. To get some idea of the effect of various parameters on the period and amplitude of the cycles the model was simulated using one set of parameter values. Then one of the parameter values was increased and the model was simulated again. Table 1 shows a summary of the effects of increasing each of the parameters on the cycle period and amplitude, and any constraints on the values that the parameters can take. One important constraint is the value of c the net fecundity of susceptibles, which must be greater than one for the population to survive. If c is less than one, the population is 40 0.9 B) k unstable cycles or no steady state 0.9 C) k unstable cycles ^r--^or no steady state 0.3 ^ ^ b 3 Stable cycles Persistent oscillations ^ Stable equilibrium Figure 6. An example of the effects on the existence and period of cycles in the modified Regniere model (equations 4.11). as parameters k (vertical transmission), b (net fecundity of exposed individuals). and /1 (degree of clumping of disease) are varied. The numbers show the period of the cycles that are between eight and twelve years. In figure A) # = 1,^B) = 2, and in C) = 3. In all cases, c 3. 41 not reproducing enough to replace itself, even if there was no virus present. Figure 7 gives an example of the oscillations that are seen before and after b, the net fecundity of disease-exposed insects, is increased. When b 1, exposed organisms give birth to exactly enough organisms to replace themselves since few of them survive but those that do survive give birth to a large number of offspring (see Table 3). As can be seen, the cycles become lower and shorter. This occurs because there are many more infected organisms being born, raising the infection rate for the susceptible organisms. Consequently, the total population will not become as dense as when this fecundity rate is lower. The cycles (as graphed on a log scale) look very different than natural population cycles, both in range and in pattern. Next, a free-living virus pool, V, was added to the modified model (equations 4.11). The virus has a survival rate, a, and w propagules are removed per susceptible as the susceptibles become infected. Infected individuals contribute to this pool when they die and release 7 virus propagules. Thus, the equation for the virus is: Vt+1 = av V + 7It — wPt St • (4.12) Since organisms are infected either from the free-living virus or from direct exposure to the dead infected individuals, the probability of contacting the infection should be changed to reflect this. Pt = (1 — e ht — where hi and -h„Vt )0 (4.13) h t, replace the h of equation 4.10, and concern the rate at which a suscep- tible contacts the virus. 42 Parameter b c k h a av 7 w Meaning net fecundity net fecundity vertical transmission contact rate clumping degree virus survival virus production virus removal Period — — + = + + + + Amplitude — + + — + + — + Range <c >1 <1 >1 _ <1 Table 1. Effect of increasing different parameters in the modified Regniere models (equations 4.11 and 4.12). A `+' indicates that the period length or amplitude height increased as the parameter increased, a `—' indicates that it decreased and an ' =' that it did not change. The range indicates any constraints on the parameter, if any. All parameters are > 0. 43 Effect of varying net fecundity of virus-exposed individuals 0 10^20^30 40 50 Generations Figure 7. Effect of varying net fecundity of virus-exposed individuals on equations 4.11, the modified Regniere model. Increasing b gives lower amplitude, shorter period cycles. 44 The new model includes equations 4.11, 4.12, and 4.13. Cycles could occur but their periods and amplitudes were very sensitive to changes in some parameter values. As shown in Table 1, increasing the amount of virus produced by a dead infected insect (7) only decreased the height of the peaks. If the survival rate of the free-living virus (u,,) was decreased, cycles would appear in cases where there had previously been no cycles. As in the previous case (equations 4.11), varying b, the net fecundity of infected insects, changes the period and amplitude of the cycles. In addition, the existence of cycles was very sensitive to the value of 0. Again, cycles with a 8-12 year period would only appear for a narrow range of parameter values if was much different from one. Figure 8 shows the effect of varying k, b, and on cycle existence and period. The parameters hi and hi, only affected the amplitudes of the cycles and not any of the observed dynamics In terms of the biology of the system, these results indicate that if there is more virus present, the cycles should be shorter and the peak populations lower. This is because susceptible organisms would have a greater chance of becoming infected. Infected and disease carrying organisms have lower survival and fecundity rates, so the population size starts to decline earlier (and does not grow as large). This argument fails only with vertical transmission (k). By increasing vertical transmission, the population is able to reach higher levels at the peak of the cycle. Although this seems counter-intuitive, upon examination it can be seen that the low point of the cycle decreases as the vertical transmission increases so the average population level decreases. Since the population cannot increase as quickly, it take longer for the population to reach peak levels and 45 B) 0.9 k unstable cydes or no steady state C) 0.9 k unstable cydes or no steady state ^ ^ 3 b 0.3 Stable cydes Persistent oscillations Stable equilibrium Figure 8. An example of the effects on the existence and period of cycles in the model which includes a free-living stage of the virus (equations 4.11 and 4.12), as parameters k (vertical transmission), b (net fecundity of exposed individuals), and /3 (degree of dumping of disease) are varied. Numbers show the period of the cycles that are between eight and twelve years. In figure A) = 1, in B) # = 2, and in C)/3 = 3. In all cases, c = 3. 46 the cycles become longer. The only exception occurs if transmission is perfect, in which case the population quickly dies out, since every insect that gets exposed to the virus gives it to her offspring, who will then die before reaching reproductive age. Even with changing parameter values in the model that included the free-living stage of the virus (equations 4.11 and 4.12), the cycle periods were usually longer than those of most Lepidoptera (see Figure 8), and, as Figure 8 shows, the range of parameter values that produced oscillations of the right period was narrow. Since equations 4.11 were the only modification that produced cycles of the right length for a fairly wide range of parameters, I concentrated on this version of the Regniere model. Regniere stated that his results were not dependent on the form of p t (the probability of becoming infected), although changing some parameter values such as 0 changes other parameter values for the occurence of cycles. This was tested by comparing three different equations for pt in which P t depended only on the size of the infected population and not the size of the healthy population. The first p t is the original one, seen in equation 4.10: Pt = (1^c hit )13 (4.10) The second is Pt = ( hIt 1+ (4.14) and the third is the negative binomial equation as described by May (1978): hit Pt = 1 — (1 + — ) — '3 (4.15) In all three equations, as the number of infectives becomes large, Pt —> 1, regardless of 47 the size of the healthly population. A graphical comparison of these models is shown in Figure 9 and a comparison of their results is in Table 2. Oscillations only occur in two of the cases although not necessarily for the same ranges of parameter values. For example, if equation 4.10 is used, stable cycles essentially only appear when /3 is close to one. When equation 4.15 is used, persistent oscillations occur if /3 > 2. For both of these equations, the value of /3 that produces cycles in the model makes the value of P t increase rapidly as the number of infected individuals increases. Generating cycles is more complicated when using equation 4.14. Although by varying the parameters h and /3 it was possible to get P t to resemble the curves from equations 4.10 and 4.15 (Figure 9), even after changing the other parameters the model did not produce cycles. So, the existence of cycles is dependent on how the probability of becoming infected (P t ) is modelled. Further modifications of P t were done to include some dependence on the healthy population. Various forms were tried (Table 2). One of these assumes that the likelihood of a susceptible organism contacting an infected larva is related to the total number of encounters with all larvae (N t ): Pt hIt ^ 1 + hNt (4.16) where N is the total population (S I). Thus, the probability of becoming infected depends on the proportion of infected individuals in the population. If there are very few infected organisms (compared with the total population), the probability of a susceptible individual encountering an infected one is small. Using this Pt , no cycles occurred; the susceptible organisms always took over the population. This was because if the number 48 Pt (1 _ c hit r i^hIt^V3 k 1.-1-rt ,/ 1 — (1 + li L) 1 — ^ (1 hIt 1-EhN +^ v hst —v St Equation 4.10 4.14 Cycles? if /3 near 1 existence highly sensitive to parameter values 4.15 oscillations if 0 > 2 4.16 4.17 none narrow range of v Table 2. Summary of the various forms of Pt described in the text and results when simulated in equations 4.11 (the modified Regniere equations). 49 I^ O^10^ I 20 it Inbled Figure 9. A comparison of the different equations for the probability of becoming infected that were used in the modified Regniere model (equations 4.11). In equation 4.10, 0 = 1 and h = .5, in equation 4.14, 0 = 1 and h = 1 and in 4.15, 0 = 2.25 and h = .5. Cyles occur for equation 4.10, persistant oscillations in 4.15 and no cycles in 4.14. 50 of infectives was low compared to the total number of hosts, the chances of contacting an infected individual would be small (Pt small). Few new infections would result, allowing the susceptible population to increase and become a larger portion of the total population. Eventually, the infection would disappear. The negative binomial form for pt that is used in equation 4.15 can also be modified to include some dependence on the number of susceptibles present. Hassell (1980) noted a relationship in his data between the number of hosts (susceptibles in the present situation) and the degree of dumping of the pest (virus here). He therefore decided that the dumping parameter (0) should not be a constant. The new form for p t then became: S hIt s Pt = 1 — (1 +^ry t t (4.17) where v is the slope of the relationship between the degree of clumping and the number of susceptible organisms. This is equivalent to stating that as the number of organisms increases, they become more randomly distributed or that the infection becomes less dumped. Cycles occurred when the model was run with this pt (equation 4.17) but the cycles were highly dependent on the value of v that was used. If v was small (less than .5), the population stabilized. If a larger v was used, the population would show unstable cycles. There was a narrow range in between these two values that produced stable cycles. As can be seen, the manner in which the probablity of becoming infected is modelled has a strong effect on the behaviour of the model. Biologically, this is a problem. For this model to be realistic, it is necessary to decide what form this probability would 51 have. Ideally, this could be determined from data of population sizes and infection rates. However, upon examination of Figure 9, it is easy to see that each of the three functions would fit data equally well (assuming first that these data could be collected), and that it would be virtually impossible to determine which, if any, was the true probability function. Two hypothesis can be made about the model. First the model is a close representation of the tent caterpillar system. In this case, each of the functions for the probability of infection could be tried in the model in turn, and the one that caused the model to produce cycles similar to those that are seen in nature would be the one that represents what is actually happening. The second, more likely, hypothesis that could be proposed is that since the model is so sensitive to how the probability of becoming infected is modelled, the model must not be an accurate representation of the important parts of the Lepidopteran system. I originally ran the models with arbitrary parameter values, usually based on parameter values used in Regniere's paper. To see if the model could produce cycles using realistic parameter values, the parameter values were changed to reflect actual data from tent caterpillars more closely. This was done by simply choosing constant parameter values that seemed to be within the ranges of the means for different years. For example, average fecundity ranged from 140 to 220 eggs per batch, the early survival rate (eggs to fourth instar) varied between .2 and .6, and the late survival rate (fourth instar to adult) ranged from .005 to .2. I assumed that the susceptible larvae that had been exposed to the disease would have lower survival than those that had not been exposed 52 to the disease and if they survived would have lower fecundity. Thus the parameters b and c, the net fecundities (survival x fecundity), could be estimated from the data by choosing lower fecundity and survival rates for the exposed insects than for the susceptible insects. Table 3 lists parameter estimates for fecundity and survival. Parameter b, the net fecundity of exposed insects was calculated by multiplying together all the values in the column marked 'low' while c, net fecundity of susceptible organsims was calculated by multiplying the columns labelled 'high'. Table 3 only reflects how two of the parameters were chosen. Most of the other parameters were the same for all simulations. The parameter h, the rate at which healthy larvae contact disease propagules, has no effect on the period (Table 1) so its value is inconsequential. I chose /3 = 1 since from previous simulations, that was the value that resulted in the highest incidence of cycles (for example, see Figures 6 and 8). For each combination of b and c, I allowed k, the vertical transmission rate to vary between 0 and .9 (i.e., for each set of b and c, I simulated the model ten times, once for each k). The parameter values that I used, and the results from the model are shown in Table 3. When using parameter values that were within the ranges of the means mentioned above, cycles occurred in only one instance and for only one level of vertical transmission (k). I also chose to use some fecundity values that were lower than any of the means. This seemed reasonable since a low mean value allows for the possibility of some of the individual fecundities being quite low. Cycles occurred in two of the three instances mentioned in Table 3. In one case, the cycles persisted for a wide range of vertical 53 Fecundity low high 50 200 50 200 100 200 120 210 140 220 140 250 150 200 Early Surv. low high .3 .5 .3 .5 .5 .3 .3 .5 .46 .34 .2 .5 .3 .6 Late^Surv. low high .2 .3 .2 .04 .2 .04 .01 .01 .2 .04 .02 .01 .2 .3 Parameter b c 3 .6 1.2 .36 1.9 .28 9 30 20 20 1.05 10 2.5 36 Cycles? 8-12 yr no k = .3-.9* k = .9* no k = .8* no no * (8 or 9 yr cycles) Table 3. Parameter estimations for b and c (the net fecundity of exposed and susceptible individuals) based on data, and the results from various simulations of model 4.11 using b, c as given in the table, /3 = 1 and k varying from 0 to .9 (i.e., each set of b and c, was simulated ten times, once for each k). 54 transmissions levels. The fact that oscillations occurred more frequently when a very low fecundity level is assumed indicates that if the cycles are a reasonable representation of the natural cycles, the virus must have a large impact on the fecundity of the exposed organisms. Although cycles can be generated with realistic parameter values, there are two factors about the cycles that make it unlikely that the model is realistically capturing the dynamics seen in the field. The first factor is that the cycles in the modelled population had a very different form from that seen in natural populations. Real population oscillations tend to show a few years at the peak and a few years near the bottom of the cycle (see Figure 14) rather than the sharp peaks and valleys that are demonstrated in Figure 7. Most insect-virus models show cycles similar to those in Figure 7. The second point is that cycles were not the general case when using relatively realistic parameter values (or even unrealistic ones) and the existence of the cycles is sensitive to the parameters that are chosen. Table 3 shows that changing the survival rate at only one stage can alter whether or not cycles are present. For cycles to be considered to be possible in the field, they must be robust to a range of parameters since in the real world, the parameter values differ each year. If virus is important in governing the cycles, there must be some important part of the virus-insect interaction that is missing from this model. 55 4.5 Diekmann and Kretzchmar Recently, Diekmann and Kretzchmar (1991) produced a continuous time disease model that is based on two classes of organisms (not necessarily insects), susceptibles (S) and infectives (I), and incorporates various aspects of birth, death and infections. Although continuous time models are not usually applicable to univoltine organisms such as insects, I adapted the model to fit some of the characteristics of the lepidopteran system. The infection rate 0 is a function which depends on the proportion of infectives in the total population. KI where IC (4.18) represents the contact rate between susceptibles and infectives. There are two types of mortality: natural mortality (it) and mortality caused by being infected (a). An unusual aspect of the Diekmann and Kretzchmar model was that the cyclic dynamics depended on the way in which reproduction was modelled. The authors assumed that both male and female infected adults would contribute to some reduction in fecundity, and that this effect was multiplicative (i.e., if a susceptible and an infective adult mated, the reduction in fertility would be 4" while if two infecteds mated, this reduction would be e) . This leads to the system of equations: ,5/ + e /2 µS — OS dt^S + I^ dI dS 0 S2 + 2 (4.19a) Tit- = os - ( a + p)I (4.19b) The authors show that the behaviour of the model can be summarized using six different phase-plane diagrams (Diekmann and Kretzchmar, 1991). Which diagram is 56 applicable to any case depends on the parameter values that are used, especially the parameter K, (contact rate). The authors give complicated formulas for determining, without simulation, what the behaviour of the model will be based on the parameters, and they describe six possible behaviours. Only two of these six cases contain the possibility of cycles. In both cases, the value of 4 is small, less than .5 (meaning that . the disease has a large effect on the fecundity of infected organisms). One of the cases, the cycles are independent of initial conditions and may be biologically relevant for some system. In the second case in which cycles occur, the behavior is dependent on the initial conditions as well as the parameter ranges. This second case has no clear biological interpretation. Contrary to the model, in the forest Lepidoptera only the females contribute to any reduction in fecundity if infected. However, the existence of cycles is entirely dependent on the form of the birth term (Diekmann and Kretzschmar, 1991). When pairwise formation is not included, or when only the infected mother causes reduction in fertility, the term 0 S2 + 2V - S + 4-2 12 (S + .0 (4.20) from equation 4.20a becomes 8 (S + 4 I). - (4.21) It can be shown that when equation 4.21 is substituted into equation 4.19a, only two of the six behaviours described by Diekmann and Kretzchmar (1991) are possible, neither of which include cycles. In one case, the population stabilizes at an equilibrium while in the other, the population grows exponentially although the proportion of infected 57 individuals in the population stays the same. From various simulations using different forms for birth, I found that in general, if the birth term is written as 13Sf(x) where x a necessary condition for cycles to occur is f(x) < 1. This does not mean much biologically but if equation 4.22 is rewritten in this form as birth = ps f (x)^ (4.22) with ^f (x) = 1 4 x^ . (4.23) it is easy to see that for all positive values of 4, f(x) > 1 in equation 4.23. Removal of pair formation eliminates the major cause of reduced population at high levels of infection. When the number of infectives is high, the birthrate under the pair formation model (using equation 4.20) is quite low while in a modified version, equation 4.21 for example, the birthrate is relatively high. This is demonstrated in Figure 10a. To compensate for this, the infection term (0) can be altered from its form in equation 4.18 to: c + s2 +12 (4.24) in an effort to reduce the population at high levels of infection. A visual comparison of the two forms of 0 can be seen in Figure 10b. As 0, the rate of infection, increases, more susceptibles will become infected thus increasing the death rate, reducing fecundity, and reducing the population level. Even these modifications did not cause the model to produce cycles without pair formation. 58 20^40^60 ^ ^ 100 80 # Infected 0 ^^ ^ ^ 10 5 15 20 # Infected Figure 10. (a) A comparison between the birth term which assumes that both infected males and females contribute to the reduction in fecundity (equation 4.20) and the birth term in which only infected females affect reproductive output (equation 4.21). The second equation increases with the number of infectives much more rapidly than the first. In the graph, e = .3. (b) A comparison of the different forms for the infection rate ek that was used in the model after the birth term was modified. Note that equation 4.24 gives a much higher infection rate than the original, equation 4.18. 59 5. A NEW MODEL 5.1 Introduction Models of disease dynamics in the literature either are not appropriate to virus-lepidopteran interactions, or oscillations no longer occur when the models are modified to resemble the biology of forest Lepidoptera more closely. Therefore, I developed a new model which was based on the biology of the tent caterpillar. If the model were to produce reasonable results with one lepidopteran species, it could be adapted for other species. The formulation of this model is based in part on the idea of Hochberg et al., (1990) to separate the within-season dynamics (infections) from the between-season dynamics (reproduction). This leads to one of the most important differences between this model and any of the models that were discussed in Chapter 4 (except Brown's, (1984)) which is that the probability of infection can vary during the larval growing season as the number of larvae, the susceptibility of the larvae, and the amount of virus changes. This chapter describes the new model. The first section derives the model and gives the basic equations that are used. The second section discusses how the basic parameter values were derived. Results from the simulations and the analysis are in the following chapter. 60 5.2 Formulation of the model 5.2.1 Model overview The tent caterpillar-virus model is concerned with the dynamics of individuals within different family groups (tents). There are three types of larvae: susceptible (S), exposed (E), and infected (I). Susceptible larvae have not come into any contact with the disease, while exposed larvae have actually consumed a full dose of the virus but have not become infected. The exposed larvae could also be called carriers since even though they do not die from the virus, they can pass it on to their offspring. In the model, there is no immigration or emigration of caterpillars or of virus. In nature, infection can occur in two ways: horizontally or vertically. In the model, horizontal transmission takes place in the larval stage as the infection is passed from one larva to another. Vertical transmission occurs from mother to offspring and is represented in the model during reproduction. The model is divided into three periods: larval, pupal and adult. Schematically, the overall model is shown in Figure 11. All parameters and variables used in the model are also described in Table 4. 61 A) Eggs overwinter reproduce, die e Adults ^Virus Larvae 3,E1 (Instars 1-5) V a" v death Pupae (1-a p ) B) death 1- a s /- a e P4 , V , death . death 4 1 death Figure 11. The lifecycle of the tent caterpillars with the corresponding parameters and variables from the model (see Table 4). A) Overall lifecycle. B) Larval stage. 62 Variable S E I V T H Parameter t y as a, cry aP Pt 0 q - - wt 'Yt a dy a .f.9 fe k c Meaning susceptible exposed infected virus tents total hosts (S Units larvae/tent larvae/tent larvae/tent propagules/tent number larvae/tent + E + I) Meaning Units Where Used instar (mini-intervals) time 5.1 generation time 5.8-10 survival rate of susceptibles 5.1 survival rate of exposeds 5.1 survival rate of virus 5.1 survival rate of pupae 5.6 probablility of exposure to virus 5.3 degree of dumping of virus 5.3 proportion of exposures resulting 5.1 in infection virus ingested propagules/larva 5.4 virus produced propagules/larva 5.5 scaling factor 5.5 density-dependence 5.7 threshold at which density5.7 dependence becomes important eggs hatched to susceptibles #/susceptible 5.9 eggs hatched to exposeds #/exposed 5.10 vertical transmission 5.10 proportion of new tents in 5.11 areas with no virus Table 4. Variables and parameters used in the model. 63 5.2.2 Larval Stage The first stage is the larval stage. Since there are five instars, this stage is divided into five mini-time-steps, designated by the subscript t. The interactions between the virus and the caterpillars occur within each such interval, and are pictured in Figure 11b. Since every interaction in the larval stage takes place within a year, the year (or generation) subscript (y) does not change. Each interaction also takes place within a tent or a type of tent (described in Section 5.2.4). Within the course of each instar, the larvae have probability P t of coming in contact with the virus, and thus becoming exposed or infected (a proportion q of larvae that contact virus become infected, the rest (1 — q) become exposed), or a chance 1 — Pt of avoiding the virus and remaining susceptible. The larvae will survive from other mortality sources at a rate u, or (r e , depending on whether the larvae were considered susceptible or exposed at beginning of the instar. Thus, the number of susceptibles at the end of an interval is determined by the proportion of susceptibles that survive the instar (u s ) and that did not come into contact with any virus (1 — pt ). Sy,t+1 (1 P-t)Crs Sy,t (5.1a) The probability of infection, Pt that is used here depends on the total population and on the amount of virus present, and is explicitly defined later. The assumptions are similar for the exposed larvae. Of the survivors, some will not have contacted any virus (1 —po t ), or will have contacted virus but not become infected (p t (1 — q)). In either case, they will remain in the exposed class. There are also some 64 surviving susceptible larvae that contacted virus but did not become infected. E y , t+1 = pt (1 — q)( 0" s Sy ,t Ue Ey ,t) + ( 1 — Pt)C e Ey ,t (5.1b) The number of infected larvae is based on the number of surviving susceptible and exposed larvae that become infected after contacting some virus. The infected larvae all die during the time interval. Thus, the equation for the infecteds is: ly , t +1 = pt g(u8 Sy3 t 0-eEy,t) (5.1c) Virus is removed from the environment by organisms when they contact the virus (since infections and exposures result from ingesting the virus) and through death. Infected larvae release virus into the environment when they die. Consequently, the amount of virus at the end of a period is determined by the surviving virus (Cry Vy,t )5 minus the virus that was ingested by susceptible and exposed larvae (w t ), plus the virus that is released by the newly infected organisms when they die (-y t ). Vy,t-1-1 = Civ Vy,t 'Yt ,t+1 WEPt(0 Sy,t Cre Ey ,t) (5.1 d) The amounts of virus ingested or released by larvae vary with the instar. Their forms are given in equations 5.4 (wt ) and 5.5 (7 t ). Since all infected larvae die within the instar period, and since the amount of virus depends on the new number of infecteds (4, t+1 ), equation 5.1c could be eliminated and substituted into the virus equation (equation 5.1d). The term P t , the probability of exposure to the virus, which governs the transmission of the virus between organisms within a tent (horizontal transmission) varies with 65 time. The expression for pt is based on the distribution of the virus. Suppose there is Vy , t amount of virus present near a tent containing Hy , t larvae. Of this virus, there are Vy , t /w t doses (i.e., there is enough virus to infect at most Vy , t /w t larvae). Assume these doses are distributed near the larvae according to the negative binomial distribution with a mean Vy , t /(w i liy , t ) and clumping parameter 0. A low value of /3 implies that the virus is highly dumped near the larvae, while if p is large, the virus is more randomly distributed. In practice, it is very difficult to estimate 0. The probability of the larvae escaping all the virus is given by: Po = (1 + VY ' t ) -13 . 13Wt Hy ,t (5.2) Therefore, the probability of a larva encountering the virus at any time is: Pt = 1 - (1 + „, VYT' _tT )-# pw-tii y ,t (5.3) where Hy , t = Sy ,t + Ey,t + 4, t , the total population. There are other functional relationships that vary with time. The amount of virus that is ingested by an organism (w t ) increases with time. This is equivalent to increasing a larva's resistance to the virus as the larva grows. The amount of virus produced by a larva when it dies ('yt ) also increases with the size of the larva at death. The actual equations for these relationships is not important to the dynamics as long as both are increasing and the amount of virus produced by a larva is greater than the amount ingested (ry t > w t ). After discovering this, I used the forms: lot = 10t -4^(5.4) 7t = awt^ 66 (5.5) where t is the number of the instar and ranges from 0 to 4 and a is some number greater than one, used to ensure that -y t > w t . The derivation of equations 5.4 and 5.5 and their relationship to data is discussed further in Section 5.3, equations 5.16-5.18. 5.2.3 Pupal Stage The second part of the model involves the pupae. I assumed that any infected larvae die and do not pupate, while the exposed and susceptible larvae pupate. Pupal survival (rnp ) was the same for exposed and susceptible organisms. Sy ,6 = Crp Sy ,5 (5.6a) Ey, 6 = CfpEy 5 (5.6b) 4,6 , = 0 (5.6c) On the right hand side of the equations t = 5, representing fifth instar organisms (there were five intervals in the first part of the model). Adults, denoted by the subscript 6, emerge from the pupae. This part of the model is not essential to the general dynamics, since it only scales the numbers of susceptible and exposed larvae downwards and removes any remaining infected larvae. Therefore, this segment could easily be combined with the third section of the model. 67 5.2.4 Adult Stage The final stage of the model is the adult stage, in which reproduction and overwintering occurs. Overwintering is only implied in the model since the winter mortality of eggs is included in the birth terms, f 8 and f6 . Although it is unrealistic, I assume that the virus has no mortality during the winter. The consequence of changing this assumption is discussed in the next chapter. The number of surviving adults is found for each type: susceptible and exposed. I assume that half of the adults are male and do not produce new tents. This assumption does not affect the overall dynamics. At this stage there is some density-dependence (d y ) which is based on the total number of tents from the previous year as follows: dy a + Ty • (5.7) The term a determines the level at which the effects of density-dependent become important. The term dy will be small if there is a large number of tents and consequently, it accounts for possible effects from crowding such as increased mortality or inability of the females to find viable oviposition sites. This term will affect the number of adults that will lay eggs and hence, the number of new tents. Since a tent is formed by the offspring of a single female, the number of tents each year is directly related to the number of females that lay eggs which hatch. The number of tents can thus be determined by adding the number of adult susceptible and exposed individuals (denoted by the subscript 6) in all tents, dividing by two for the number of adult females, and accounting for density-dependent effects using dy : Ty+1 = —1 (S 6 + Ey,6).^ Y 2 68 ' (5.8) For most simulations of the model two types of tents are used: those that contain only susceptible larvae and those that contain both susceptible and infected individuals. Female susceptibles produce egg masses which, when they hatch, will contain only susceptible larvae. These tents will be in areas where there is virus already. The number of eggs that a female susceptible lays that hatch is denoted by the parameter fs • Type 1 (susceptible, virus present): Sy+1 ,0 = d 2 Vy+1,0^Vy ,6^ (5.9a) (5.9b) Adults exposed as larvae produce offspring that are a combination of susceptible and infected, the proportion determined by a parameter k (vertical transmission). I assumed that all these egg masses were laid in areas that already contained some virus. Exposed adults lay fewer eggs than susceptible adults due to sublethal effects of the virus. This gives the second type of tent: Type 2 (combination, virus present): Sy+1 , 0 = -2-(1 — k)fe Ey , 6 (5.10a) dy ly +1 , 0^ICJ e-C-fy,6 (5.10b) Vy+1,0 Vy,6 + 74+1,0 (5.10c) 2 2 where f e is the number of eggs that an exposed female lays that hatch. The values on the left-hand side of equations 5.9 and 5.10 are the number of larvae in the next generation in each of the tent types. 69 In determining both types of tents, I assumed that they are in areas that already contain some virus. At low population densities that may be an erroneous assumption. Consequently, in a few simulations, I also examined the behaviour of a model in which I had added two other types of tents: ones that contain only susceptibles and ones that contain both susceptibles and infecteds. Neither of these types of tent were in areas that already had virus present. In determining all four sets of tents, I decided that a proportion c of each type of tent would be laid in areas with no virus. The new set of equations are: Type 1 (susceptible, virus present): sy+ 1,0 = 2 (1 — cVsSy , 6 ^ (5.11a) Type 2 (combination, virus present): Sy-f-1,0 dy — (1 — c)(1 — k)fe Ey , 6 2 (5.11b) -4+1,0 c (1 — c)k fe .Ey , 6 A 2 ( 5. 1 1c) Type 3 (susceptible, no virus present): S Y+1 0 = --v-cf s Sy,6 (5.11d) ^2 Type 4 (combination, no virus present): Sy +1,0 = -4+1,0 = 4 41 — k) feE 2 ck fe Ey , 6 y ,6 ^ ^ (5.11e) (5.11f ) The presence of tent type 4 (tents with both susceptible and infected larvae but with no virus present) did not affect the dynamics, and only marginally affected the level 70 at which the population stabilized. However, tent type 3 (tents containing susceptibles only with no virus present) did affect the dynamics. Type 3 tents, since there was no infection present, either in the larvae or in the environment, acted as a refuge from the virus, causing a higher population level and more stable dynamics. Consequently, most simulations were done assuming only two types of tents (equations 5.9 and 5.10). If cycles were to occur, they would be most likely to arise with only the two tent types present. 71 5.3 Parameter values Most of the parameter values that were used in the simulation were originally based on values mentioned in the literature, mainly from Entwistle (1983a) and Myers (pers. comm.). The parameter values that are given below are a basic value used in most simulations, and in all cases, simulations were also done with different values (See Section 6.2, Table 5). The survival rates of the larvae at various stages are loosely based on data from Myers (pers. comm.). In the data, the averages of the early survival rates of the larvae (up to the fourth instar) range from .12/season to .82/season. Since some of the mortality is due to natural causes and some to infection, the values for a , and - Te, the survival rates of the susceptible and exposed larvae respectively, were chosen from the higher end of the scale and were calculated by choosing values for early survival rates of susceptible and exposed larvae as .7/season and .5/season respectively. These values are then raised to the power 1/5 to get the survival rate per instar, assuming that the survival rate is the same at each stage. Thus, a, = .7 1 / 5 = .93/instar (5.12a) cr e = .5 1 / 5 = .87/instar. (5.12b) The survival rate of the virus can be calculated either from its half-life or from the proportion of virus remaining after a certain length of time. Entwistle (1983a) states that NPV virus has a half-life of one day but that there is more than 1% of free-living virus left after twenty days. For both these statements to be true, the death rate of the virus cannot be a strict exponential decrease. In fact, virus mortality usually depends 72 on how exposed to the sun it is. However for convenience, exponential decay will be assumed. Given the statements above, two possibilities arise. Using the hypothesis that after twenty days there is is at least 1% of the virus left, the survival rate of the virus, a, during one instar (ten days) is: o , > .1/ten day period (5.13a) - If the half-life of the virus is one day, the survival rate becomes: cry = .001/period. (5.13b) In most simulations, a value of .01 was used, although the effects of the higher and lower a, values were also explored. In the model, I assumed that the pupal survival rate was the same for exposed and susceptible organisms. From the data, late survival rates (which includes pupal and adult) range from .001/season to .2/season. For most simulations, the pupal survival rate was taken to be o-1, = .1/season. (5.14) At the time of reproduction, I assumed that a proportion k of the offspring from exposed adults are infected at birth. The value of the parameter k cannot be determined from the data so a value of .7 was used for most simulations (i.e., 70% of the offspring from exposed organisms are infected at birth). The number of eggs laid per female is given in the data with a range of 140 to 230 eggs/batch. I assumed that the exposed adults would lay fewer eggs than the susceptible adults (140 vs 230), and that only 60% of the eggs from either type of adult would hatch. The parameters f8 and f, reflect the 73 number of eggs laid and hatched by susceptible and exposed females respectively, and for most simulations were determined to be: ^f, = 140^ fe = 100 (5.15a) (5.15b) As mentioned before, the amount of virus ingested and released by a larva varies with the age and size of the larva. Although experiments for determining the amount of virus needed to infect a larva 50% of the time (LD 50 value) are not always reliable, they give some indication about the levels of virus that may be important. Entwistle (1983a) gives a range of lepidopteran LD 50 values as from 10 propagules when the larvae first hatch to 10 6 propagules for fifth instar larvae. (A propagule is a virus polyhedral inclusion body (PIB) which is the virus surrounded by protein and is visable with a microscope (Entwistle,1983a).) Assuming that the infection rate increases by an order of magnitude at each stage, the relationship for the amount of virus ingested per larva (million propagules/larvae) as given in equation 5.4 is wt 10t -4 where (5.16) t indexess the instar and ranges from 0 to 4. Since w t is the amount of virus that will cause infection 50% of the time, the proportion of contacts resulting in infection (q) was taken to be 50% (or q = .5). Entwistle (1983a) also gives information for the the amount of virus per milligram that lepidopteran larvae can produce and for the range of total virus per larvae that can be produced. Since the model does not calculate the weights of the larvae in the tents, 74 the range was initially used as the basis for determining how much virus was produced by the infected larvae. Entwistle (1983a) gives the range as 10 9 to 1.5 x 10 10 . I assumed that the lower end of this range corresponded to the larvae that died while in the first instar and that the upper value was the amount produced by the fifth instar larvae. I also assumed that the progression from one to the other was linear. Thus, measuring the amount of virus produced by an infected larvae in millions of propagules per larvae, 7t = (1 + 3.5t) x 10 3 (5.17) where t is a measure of age and ranges from 0 to 4. When the model was simulated using this relationship for the amount of virus produced, the population usually died out because of these high production levels and because I assumed that all virus that was released, and that did not die, was available to the larvae. To alleviate this difficulty, I changed the production curve to make it a function of the amount of virus ingested. Thus, as given in equation 5.5, 7t = aw t (5.18) where a = 4.5 (a value determined through simulation). This makes the lower end of the range for the virus produced much lower than given by Entwistle (1983a). Since there were no data for the degree of clumping of the virus (0), the level at which density becomes important (a), and the proportion of organisms born in an environment with no free-living virus present (c), the parameter values used were chosen for other reasons. The parameter /3 must be greater than 0 and I assumed that it must also be relatively small (compared with oo) so I chose 3 = 1 for most simulations. The proportion of larvae in areas with no free-living virus is also likely to be small so I used 75 c = .2 for the simulations in which more than two tent types were possible. I let a = 500 since for most simulations, it caused the population to be near levels that were close to recorded population sizes. 76 6. SIMULATION RESULTS AND ANALYSIS 6.1 Introduction To determine the possible behaviours of the model, I simulated the model on the computer and analyzed a simplified form of it. Although the analysis of the simpler model did not indicate that cycles were possible, oscillations occurred in some simulations of the full model. However, these oscillations had different characteristics than are observed in natural populations. In the first section, the results from various simulations are examined. The second section lists the important results from the analysis of a simplified version of the model and their interpretations and implications in biological terms. The actual analysis is shown in Appendix 2. The applicability of the analysis to the actual model is also discussed. In all cases, as a first step in deciding whether or not the model 'worked', cycles, preferably ones with an eight to twelve year period, were the criterion. If this criterion were fulfilled, I then examined the predicted and the average clutch sizes that the females laid and related them to the ranges given in the data. Finally, I also compared some of the patterns in the changes in fecundity and survival that arise during the model's cycles with those in actual cycles. 77 6.2 Simulation Results Simulations were done with the basic parameter values derived in Section 5.3 and by varying the parameter values individually as shown in Table 5. In all cases the population stabilized at some level, either at zero (no organisms present) or at a higher level. This population size was always stable; if the population was perturbed away from that level, it returned rapidly. Table 5 shows the effect of varying each parameter individually on the level at which the total population H or the amount of virus V stabilized. Persistant oscillations did not occur under any circumstance that was examined when only one parameter was varied within the ranges indicated in Table 5 and all other values were at the base levels shown in the table. Some of the behaviours given in Table 5 may seem counter- intuitive. One such case is the fact that when the virus survival rate (up) is increased (implying that there should be more virus present), the amount of virus present when the population stabilizes actually decreases. This occurs because the caterpillar population has also declined. Thus, fewer larvae become infected, and they become infected at a younger age. Consequently, not as much virus is produced since in the model, the only way the amount of virus can increase is by being eaten by a larvae which dies. A similar argument can be made for the effect of the vertical transmission rate (k). The behaviour of the model was very sensitive to some of the parameters that were used in the model. As shown in Table 5, these were: q, the proportion of exposures to virus that result in infection and a, a parameter used in calculating the amount of virus that is released by an infected larvae when it dies. The sensitivity was such that if the 78 Parameter as a, - vv oP k f8 fe - Basic Value .93 .87 .01 .1 .7 a a 140 100 .5 4.5 500 13 1 q Range Explored .75 — .99 .75 — .99 .001 — .2 .01 —.2 0—1 100 — 400 50 — 200 .3 — .9 1-10 50 — 1100 a — 10 Stable Level Sensitivity H V + — — + + + — + + + — — + — + + — — + — s s s s Cycles Period Amplitude + + + + — + + + + — — v v — + — = + + + s 15-40% change in population level with 10% change in base parameter value. v 60-90% change in population level with 10% change in base parameter value. Table 5. Effect of increasing different parameters on the level at which the population stabilizes and on the period and amplitude of cycles in the new model. The basic value was determined in Section 5.3 and was used for all simulations except where noted. The range column indicates the range through which the parameters were varied. A `+' indicates that response increased, a `—' that it decreased, and an `=', that it remained the same. Base parameter values used to calculate the cycles were: = 4, f, = 320, a , = .1 and all others as in this table. - 79 values were much higher than were given here, the population would die out and if they were much lower, the infected part of the population would disappear. Thus, all results are highly dependent on these parameters. The model was also simulated by varying more than one parameter at a time from the base levels. In some cases, eight to twelve year cycles did occur (Figure 12). In cases when oscillations occurred, I examined the effect of varying different parameters to determine their effect on the period and amplitude of the cycles (Table 5). Figure 12 gives an example of the effects of changing some of the key parameters: a y , the survival rate of the virus, /3, the degree of dumping of the virus, fs and fe , the net fecundities of the susceptible and exposed adults. The period of the cycles, once they had occurred, was not very sensitive to changes in the parameter values (Figure 12). In general, the cycles increased in length as a, and were increased and decreased slightly as f3 and fe were increased. Increasing the susceptibles' birth rate (18 ) allowed the population to recover more quickly from the decline since the susceptibles were the largest proportion of the population when the population was at low levels. The cycles occurred when 13, 18 , and cr, were much higher than the base values given in Table 5. However, this does not necessarily indicate that the parameter values that produce oscillations are unrealistic. For example, the parameter 13, the degree of clumping of the virus near the insects, cannot be determined from the available data and would be difficult to determine empirically. I had originally assumed that since during the season the virus is found primarily near tents (since the larvae do not wander far from the tent), that /3 would be small, close to one. The values of /3 that are shown in 80 fe 50 s ma v ^ 6^12^12 11.5^11^11^11 5 ^12 11.5^11^11^11^11 4^12^11^11^11^11^11 11^11 10.5^10.5^10^10 3^ 2 10^10 10 9.5 ma_v ..05, fe=100 160^200 B) 6 5 beta 4 3 2 1 sigma_v . .1, fe 100 10 10^10 10^9.5 10 10 9.5^9.5 9.5^9.5 160^200 C) 240^280 fs 240^280 fs 320 360 9.5 9.5 9.5 9 9.5 9.5 9.5 9 320 360 160 sigma_v . .2, fe . 100 D) 200 240 280 320 360 fe 150, sigma_v 6 5 beta 4 3 2 1 6 5 beta 4 3 2 1 160 200 240 280 320 360 fs 160 200 240 280 320 360 fs 8-12 year oscillations >12 year oscillations decaying oscillations 7 Figure 12. Occurrence and period of cycles in the model as the parameters /3 (degree of dumping of the virus), o , (virus survival rate), f, and fe (number of larvae born per susceptible and exposed females respectively) are varied. Cycle lengths have been averaged so a period of 9.5 indicates that the cycles alternated between 9 years and 10 years in length. A * indicates that the cycles were extremely irregular but tended to remain between eight and twelve years. - 81 Figure 12 are of the same order of magnitude as one and therefore may still qualify as `small' (especially since the maximum value of /3 is oo). As discussed above, based on information in the literature a v , the survival rate of the virus, might be around .001/period or .1/period. Although the base value was taken to be a, = .01/period, the higher values of o , that produce cycles are closer to - the upper estimate (e.g., o t, = .05, .1 or .2). If u„ < .01 (closer to the lower estimate), - cycles did not occur. The upper calculated value for a, was based on the observation that there was still more than 1% of the virus left after a certain amount of time. If the simulated cycles are a reasonable approximation of natural oscillations, it implies that observations about the presence or absence of virus may be more reliable than calculations of half-life. Actual virus likely does not show exponential decay since the lifetime of the free-living virus is usually determined by the amount of exposure to sun and rain it receives. If the virus is protected, it will live longer than if it is exposed. In the model, I assumed that the virus had no mortality during the winter. The model was very sensitive to this assumption. If I allowed less than 90% of the virus to survive the winter, the cycles disappeared, unless the value for the within season survival rate, cry , was increased. However, increasing a„ to a level at which cycles would again occur, usually leads to a„ values that are much larger than data suggest. Finally, the values for f,, the number of new susceptibles per susceptible female, were also much higher than originally estimated. To see if these values could be reasonable, the maximum and minimum average number of eggs per female was calculated during the course of a cycle for the instances in which Figure 12 shows that cycles 82 occurred (Figure 13). These values were then compared with the range shown in the data. In most cases, the average clutch size was much larger than the data indicated, especially under the assumption that not all the eggs hatched. In fact, almost all the instances in which the dutch size was similar to the data occurred if egg survival was perfect (i.e., if all the eggs that were laid, hatched). One exception to this occurs if the number of new larvae from an exposed female (f e ) is very small (Figure 13E). In natural populations, the egg hatch rate is near 100% (Myers, pers. comm.). Therefore, the birth parameters that produce oscillations are reasonable. Next, the cycles were examined to see if they were similar to natural oscillations in more ways than the period and the range of dutch sizes. It is not necessary to examine the predicted number of tents or population size since this can be easily adjusted by altering a, the parameter that determines the population size at which densitydependence starts to have an effect. This parameter does not affect the period or average clutch sizes of the cycles. Instead, the patterns of population increase and decrease were examined. In natural populations, it has been observed that the average dutch size does not start to decrease until the first year of the population decline (Figure 14A,C). This is not observed in the simulated cycles (Figure 14B,D). The cycles produced by the model show clearly that the average clutch size starts to decrease before the population, and that it is this reduced clutch size that causes the population decline. I then compared the early survival rates (from egg to fourth instar) of the natural populations and the simulated population since the populations in all cases declined in the same generation but the declines were not all due to reduced fecundity. The survival 83 egg survival rate = 1 sigma_v = .05, fe = 100 240 205 6 280 165 175 180 235 270 5 185 beta 195 4 A) 320 185 310 200 290 230 240 280 320 360 fs D) 6 5 beta 4 3 2 C) 6 5 beta 4 3 egg survival rate = 1 sigma_v = .1, fe = 100 175 290 330 130 150 155 170 210 290 325 140 155 130 150 170 205 285 325 140 145 155 160 200 275 315 165 155 170 200 240 280 320 360 fs C) 6 5 beta 4 3 2 egg survival rate = 1 fe = 50, sigma_v = .1 145 185 225 265 305 345 70 75 75 80 80 80 145 185 225 265 305 345 75 75 80 85 90 95 145 185 225 305 265 345 75 80 95 90 90 95 140 180 220 260 300 340 80 85 90 95 100 105 195 270 310 120 135 140 160 200 240 280 320 360 fs egg survival rate = .6 E)^fe = 50, sigma v = .1 245 6 310 375 445 510 580 120 125 130 130 135 135 240 310 375 440 5 510 575 120 130 135 beta 140 145 160 4 •^240'i 305 370 440 505 570 125::. 135 140 155 150 160 3 7167293 436406g 365 145 155 160 170 175 2 320 450 385 515 200 210 225 235 160 200 240 280 320 360 fs . ---- egg survival rate = 1 sigma_v = .2, fe = 100 290 335 105 105 290 335 110 105 240 285 330 115 110 110 240 280 200 320 120 125 120 115 190 270 310 130 130 135 135 240 280 320 360 Is F) 6 5 beta 4 3 - egg survival rate = 1 fe = 150, sigma v = .1 205 245 280 315 185 200 195 205 240 275 315 195 205 210 235 270 310 205 210 215 255 295 230 235 240 280 320 360 Is Figure 13. Maximum and minimum clutch sizes produced during each of the cycles indicated in Figure 12. In A, B, C. D, F, egg survival is 100%. In E, egg survival is 60%. Areas shaded grey roughly correspond with the actual data range of 140-230 eggs /clutch. 84 Mandarte Westham 0 C A N q N -2 0 # tents # eggs/batch 1^I 1^2 a) 1 3 1 4 1 5 1 6 # tents # eggs/batch 0 -g -C C.) CCS .0 Cl) 11111111111111 1^2^3^4^$^6^7^6^9 10 11 12 13 14 1111111111^1 7^5^9^10^11^12^13^14^IS^16^17 Model Model C) 0 0 • \ \^ 1 .....^\ \ \^n^.. ; ^\._.//^ CM 0 \ .^... B ,- . i • • • .^• .^. .^.^• .^•^. .^.^. .^.^. .^.^. .^.^. .^.^. .^. ,^. •• ^.. -2 0 # tents # eggs/batch 0 # tents # eggs/batch -8 1^I^1^1^1^I^I^I^I^1^1^I^1^I 1^23^4^56^7^59^10 11 12 13 14 11111111111111111 123455 7 5 9 10 11 12 13 14 15 16 17 Generations Figure 14. Log number of tents and average clutch size in two actual populations: A) Mandarte, C) Westham and a simulated population (B.D). Parameter values used in the simulation were: /3 = 4, f, = 280, fe = 100, a = 200 and all others at the base value given in Table 5. 85 rates in the natural population declined before those of the simulated population (Figure 15), and it is likely that these lower survival rates contribute to the natural population decline. Therefore, although the population levels are cycling with the same period, the oscillations produced by the model are a different type of cycle than is seen in nature. Some of the assumptions that were used in deriving the basic parameter values may not be valid. Therefore, some of them were changed slightly to see the effect on the behaviour of the model. The first of the changes was made to the survival rates of the larvae. The basic parameter value given in Table 5 is a constant, assuming that each instar has the same mortality rate (due to all causes other than virus) as each of the others. If instead it is assumed that the survival rate is lower for the younger instars and higher for the older instars, and if the overall survival rate is the same (.7 or .5), the dynamics are not greatly affected. I also assumed in that the amount of virus that larvae ingested (that gave them a 50% chance of becoming infected) was a function that increased linearly on a log scale. Various other functions for the amount of virus ingested were also tried, including some that decreased or remained constant as the larvae grew. Some of these are pictured in Figure 16. As long as wt was increasing, the model produced approximately the same behaviour, and this was the only instance in which cycles occurred. This is support for the knowledge that larvae need more virus to become infected as they grow larger. Originally the model was simulated by differentiating between each individual tent. Each tent initially either contained only susceptible organisms or contained a mixture of susceptible and infected organisms and each tent was in an area with differing amounts 86 ^ Mandarte 0 ri . ^ .^. .^. ^.^ :^ . . ^.^ :^ . 1/ .^ / O \ ■■• .^ 0 ..^ ... .^ .;^ ^. O ^..^.^ '....^ . '.. :^ • • - cc! .... ‘t. ): ....^ O ^..'. O : ^• - 1 2 3 4 5 6 7 8 9^11^13^15 17 ci) c a) rn Model 0 -J O (V) B ••■ • \^ • • • /^• • • O co • O • • • O •^ •^• • --• - 1 2 3 4 5 6 7 8 9^11^13^15^17 Generations Figure 15. Early survival rates (from egg to fourth instar) in a simulated and a natural population. Generations are the same as in Figure 14A and B. Parameter values used in the simulation are indicated in Figure 14. 87 ^ ....^ / ... .., ..•^ .... a^ / . . N..^. N^. . N^./ N. ^. ^ N ... -..^ -/ . . / / C /^ ./^ / .)•: . ., ./^\^ , / N. /^ /^ / /^ /^ /^ ...- 7 d z/ / N N. z / / N. N.. ...,^ ^ ./ .-,--- .•-• / / •.. /^ / / ,^ /^ % / ....°^ N N. N..N, . N. --,^ N. 0.•••• 2^ 3^ 4 ^ 5 Instar Figure 16. A comparison between the different forms for w t , the amount of virus ingested at each instar. Lines a, b and c give similar results in the model (including cycles). Line d does not give cycles. Line e usually allows the population to die out. 88 of free-living virus. The model was also simulated by only keeping track of the number of organisms in each of four different groups of tents (see equations 5.11) or in each of two distinct types (see equations 5.9 and 5.10). Finally, the model was also simulated by ignoring all distinctions between tents. The division of the model into individual tents was not critical to the behaviour of the model although the separation into different tent types was important. The model behaved in the same way using individual tents and using four different tent types. However, when the model was simulated with no separations between tents, it behaved quite differently; the dynamics were much less stable, showing more incidences of decaying cycles, and a greater likelihood of the population becoming extinct. Ignoring the differences between tents means that the population is modelled as if there are interactions between all the larvae and as if virus can move between sites. When the population is divided into different tent types, there is still the implication that virus can move between tents that are of the same type, but that there is no interaction between larvae from different types of tents. Those populations that start with both virus and infected larvae present will have more virus present, fewer larvae and therefore higher probabilities of larvae contacting the virus than tents of a different type. The only difference between dividing the population into two types of tents and four types of tents was that, as mentioned in section 5.2, using four groups of tents gave a more stable, higher level population. For example, all the cycles shown in Figure 12 were generated from the model that assumes that all eggs are laid in areas in which there is virus already present. If I assume that some of the tents are in areas in which 89 there is no virus previously (using four types of tents), cycles do not occur. There are two exceptions to this. If the proportion of tents laid in areas with no virus is very small (c < .1) or if the only eggs that are laid in areas with no virus contain infected larvae, then some cases of shorter period, smaller amplitude cycles do occur. The differences between using two or four types of tents occurs because one of the categories of tents contains only susceptible organisms with no free- living virus nearby. This class of disease-free tents acts as a refuge from the virus. The organisms in these tents have relatively high survival and birth rates, since there is only natural mortality and no disease-induced reduction in fertility. Therefore, when this class is present, it has a stabilizing influence on the overall dynamics 6.3 Analysis 6.3.1 Introduction This section examines the possible behaviours of the model determined by analysis performed on simplified equations. The details are given in Appendix 2. The first section discusses the reproductive rate of the virus, and its effects on model dynamics. Next, I calculate the conditions for stability of the simplified model, and discuss the implications for the more complex model. 90 6.3.2 Virus reproductive rate First, I calculated an approximate reproductive rate (R t ) for the virus. This is the rate at which the virus will increase in one within-year time-step. The calculations are shown in Appendix 2.2, and they give the result: (6.1) Vy ,t+1 = RtVy ,t Rt = 0" v + (7tq - wt)(0 .9Sy ,t - + Cr e E y ,t) wtHy,t (6.2) where H y ,t = Sy ,t + Ey , t + 4, 2 . (6.3) It is well known that, for simple difference equations such as the logistic, the reproductive rate, r, of an organism must be at least one for the organism to survive, and that increasing r will lead to cycles (May, 1976; see also Figure 1). However in this model, R t varies with time so the product of the Rt 's must be at least one (i.e. r = Ro R1 R2 R3 R4 R5 > 1). This means that any single R t may be much less than one. Since the calculations were done for very low levels of virus, r is only accurate as an approxiation to the actual reproductive rate of the virus if V is small, less than 10 -4 . Since virus is measured in millions of propagules (i. e., if V = 10 - 4 then there are 100 propagules present), this is not an unreasonable requirement. 91 6.3.3 Stability In order to determine analytically the possible behaviours of the model, the model had to be simplified. Since pt , the probability of contacting the virus, was the most complex term, it was simplified by setting p = a constant. The stability conditions could then easily be determined for a modified version of the model in which p was constant and in which all tent types were combined into a single type. The results of the analysis of the simplified version will not give the exact behaviour of the current model, but it may give some clues as to possible behaviours or key parameters. From Appendix 2.3, the stability condition is given in equations A2.14 and can be summarized to say that for the population to stabilize at some non-zero level, ib > 2 where 7,b fa x + (1 — ic)fe z (6.4) where x and z are combinations of parameters only (see equations A2.19 and A2.22). This can be interpreted easily in terms of the biology of the system. The term lk represents the reproduction of the organisms. In formulating the model, I assumed that only half the organisms give birth. Therefore, there must be at least two offspring per mother for the population to continue. The conditions state that if there are fewer than two offspring per mother, then the population will die out (go to the trivial equilibrium), and if there are more than two, the population will go to the upper equilibrium. If I had assumed a different sex ratio in the model, the stability conditions would reflect the new sex ratio. Basically, for the population to survive, the females must give birth to enough offspring to replace themselves and those adults that do not give birth. 92 The upper equilibrium is determined by d, the term governing the density-dependence in the model. If there is no density-dependence acting at the time of reproduction (i.e., d = 1 always), there is only one equilibrium value, the trivial one, where the population size is zero. This means that the population would either become extinct (if the trivial equilibrium were stable) or would grow without bound (if the equilibrium were unstable). Clearly, due to the way that the model is formulated, the density-dependence is necessary to stabilize the population. Simplifying the model by allowing p to be constant leads to the unreasonable assumption that the amount of virus present at any time does not contribute directly to the stability or instability of the system since the modification causes S and I to be independent of the virus (see system A2.10). However, even with the simplification, the constant value for p that is chosen can reflect the amount of virus present, so a small p implies that there is little virus present, and the population will stabilize at a relatively high level, while if p is dose to one, there is a large amount of virus present, and the population will either become extinct or stabilize at a low level. The relationship between p and the level at which the population stabilizes is shown in Figure 17. Thus, the value of p that is chosen (and indirectly the amount of virus) affects the size of the population at equilibrium, and even whether or not the population will survive. In summary, the analysis of a simplified version of the model indicates that cycles cannot occur. The analysis shows that the population will die out if the females do not lay enough viable eggs to equal the total adult population. If there are more than this number of offspring produced, the population will grow to some level determined by the 93 0.0^ 0.2^0.4^0.6 0.8 Infection probability (p) Figure 17. The relationship between the constant value of p (the probability of becoming exposed to the virus) and the number of susceptibles at equilibrium, according to the analysis. The point pc is where the equilibrium population size is zero. 94 1 .0 effects of density-dependent reproduction. The results from the analysis do not account for all the behaviours that are observed in model simulations so it can be assumed that the fact that Pt , the probability of becoming infected, varies with time is an important component in the observed dynamics. 6.3.4 Comparison with results of simulations In the first part of the analysis, the virus reproductive rate r was calculated and it was hypothesized that increasing r would lead to cycles. In fact, with all parameters at the base values given in Table 5, increasing r just lowers the level at which the population stabilizes until eventually, for high enough r, the population becomes extinct. This occurs because as r increases, there is more virus present, allowing more organisms to become exposed or infected and either to die or to have a reduced fecundity. Thus the population size will be lower. However, as observed in Figure 12, if other parameters are also moved away from the calculated base value, increasing r leads to cycles, and higher r produces longer length cycles. There are two weaknesses with the analysis. The first is that all tent types are combined into one type. As seen from simulations, the behaviour of the combined model is different from the behaviour of the separated model. The model was less stable using a single tent type. The second, more important defect of the analysis is that it assumes that p, the probability of exposure to the infection, remains constant regardless of how much virus is present. In spite of these inadequacies, the results can be extended to some degree to the actual model where pt varies during the season, and between years (equations 5.1 and 95 5.3). From many simulations, it was determined that the geometric mean of the pt for each time i within one season, will give a reasonable estimate about whether or not the population will stabilize (or die out), as long as the population is near equilibrium. For example, given a set of parameter values, the critical value pc can be estimated. The value pc is the p for which the non-trivial equilibrium equals the trivial equilibrium (i.e., the population dies out). This point is indicated in Figure 17. When the model is simulated, the geometric mean, 15 of the pt values for each step can be calculated as p = (ponp2p3p4) 1 / 5 • If p is less than pc , then the population will stabilize, otherwise it will die out. Thus, the analysis can be applied (if the population is already near equilibrium) by using p as the constant p. Cycles are observed from the calculation of p by examining its value over time. If the population is going to stabilize, p will be fairly similar in each time step. If cycles occur, then p will be less than p c for several time steps (allowing the population to move towards the steady state) and then will be greater than pc for several steps (as the population declines). Not all the results from the analysis carry over to the actual model. If the population stabilizes, it has fewer organisms present than were predicted from the analysis. This occurs for two reasons. The first is that po and p i , the pt values for the first two time-steps of the season, are often elevated due to vertical transmission of the virus and because the young instar larvae are more susceptible to the virus. These Pt values cause high levels of early mortality to a large number of organisms. The second reason is that the analysis was done assuming all tents had the same amount of accessible virus. 96 When the model is simulated, some tents have a large amount of virus present while others have little or none. This causes high variation in mortality between tent types, and usually results in having fewer organisms present than predicted from the analysis. 6.4 Summary For this new model, given in equations 5.1 to 5.10, eight to twelve year cycles can occur under certain parameter combinations (as in Figure 12). In most cases, if I assumed that not all eggs hatched, the average clutch sizes were much larger than is seen in natural populations. With the assumption that all eggs hatched, there were more occurrences of clutch sizes with a reasonable range (see Figure 13). However, the range of parameters producing cycles of the right length with reasonable average clutch sizes was very narrow. Also, the cycles showed a different pattern than is seen in natural populations (Figures 14, 15). Simulated cycles showed that the population declined after the average clutch size decreased while the reverse is true in natural populations. Also, in natural populations the survival rates dropped earlier in the cycle than in the simulated cycles. However, from the simulations and analysis, several features of structure and values of parameters are important in creating cycles. 1. Density-dependence affecting the proportion of adults reproducing. Without any sort of density-dependence in the model, the population did not cycle or stabilize at any level. Birth rates are high enough (when using parameters based on data as in Table 5) so that even with virus reducing fecundity, the population is able to grow quickly and the susceptibles become the major proportion of the population. 97 2. Sub-lethal virus affecting fecundity levels. If the virus does not affect fecundity, then f8 = fe . However, cycles never appear under this circumstance, even if there is a large sublethal effect on survival in the larval stage. It is the lowered fecundity that is the primary cause of the population decline in the model. In fact, the larger the difference between f8 and fe the , more likely that cycles will occur (Figure 12). It is not necessary to also have the sublethal effects on survival during the larval stage in order to produce cycles. 3. Survival rate of free-living virus. If there is not enough virus present in the environment (because the survival rate is too low or because it does not survive over the winter), it allows the population to escape and grow to some stable level. Conversely, if the survival rate of the virus is too high, the population could start to show erractic cycles or become extinct. The parameter k controlling vertical transmission has a similar effect in that as it increases, more virus will move into the next generation and keep the population lower. 4. Infection rate. The fact that the analysis and the results disagree implies that the fact that the probability of infection varies with time and density is crucial to the production of cycles. Another critical factor included in the infection rate is the degree of clumping of the virus near the larvae, something that is difficult to measure. 98 7. CONCLUSIONS Models were used to test the hypothesis that the eight to ten year population cycles that are found in tent caterpillars could be produced by a virus. Four insect-virus models from the literature were examined. Each of these models initially showed that population cycles were possible. All of these models lost their tendency to cycle if realistic parameter values were used. The same results were found by Vezina and Peterman (1985) who adapted to Anderson and May (1980) model for biologically resonable details of the douglas-fir tussock moth. Since none of the existing models produced reasonable cycles, I created a new virus-insect model with particular attention to the tent caterpillar system. This model included many characteristics about the biology of tent caterpillars and the interactions between the larvae and the virus. One of these characteristics not included in any previous model examined was that the caterpillars' susceptibility to the virus changes as the larvae become older. Thus, the probability of becoming infected varies during the larval growing season, depending on how much virus is available and the age of the larvae. Heterogeneity was incorporated in my model by varying the amounts of virus available to caterpillars. The model was examined both through simulations and through analysis of a simplified version of the model in which the probability of becoming infected was held constant. This analysis showed that no cycles were possible in the simplified version of the model. However, simulations of the model showed that cycles could occur. This implies that is important that the probability of becoming infected (p) varies with the 99 amount of virus and number of larvae present. Cycles occurred when using relatively realistic parameter values and the period of the cycles produced by the model corresponded well with those of actual cycles. However, the simulated oscillations showed that the average clutch size started to decline before the population size decreased and survival rate decreased at the same time as the population. In natural populations, average clutch size does not go down until the first year of the population decline while the survival rate decreases before the population decline. Thus there is a difference of one to two generations between the simulated and actual cycles of dutch size and survival rates. There are three explanations for this: 1. There are aspects of the biology of the virus or the virus-insect interaction that have not yet been discovered but that are crucial to the character of the cycles. 2. There is some factor or factors, not including virus, that is or are governing or influencing the cycles. 3. Virus works with some other component to produce the patterns seen in the cycles. Each of these arguments will be examined in turn. First, there may be some crucial, as yet unknown biology about the virus or the manner in which the virus interacts with the insects. This conclusion affects how the model is actually constructed, so under this conclusion, it is still possible that virus is the single cause of population cycles and that the only reason that proper oscillations have not yet been produced is that something is missing or inaccurate. This unknown biological factor could involve either interactions that are occurring or the measurement of key parameters. In the first case, it is possible that the virus is 100 transported between tents in some manner (by parasites or wind, for example) or that the virus has some other detrimental effect on exposed larvae besides reducing adult fecundity. Also, little is known about the free-living stage of the virus. As mentioned previously, there is some debate as to how long the virus can live in the environment, and whether the free-living form is available to first instar larvae. Finally, some of the mechanisms in the model, such as the probability of becoming infected, may change as the number of organisms increases. For example, as the population size grows, the larvae will be more likely to encounter larvae or virus from another tent and so increase their likelihood of becoming infected. The probability of becoming infected would also change if the susceptibility of the larvae changed during the cycle. Research in these areas may produce different or additional results than are portrayed in the current model. The second area included in the first conclusion is the measurement of various parameter values that seem to govern the interactions. Many aspects of the biology of the virus are difficult to observe or to measure. For example, often it is not possible to tell if caterpillars have been infected by the virus until after they have died. Therefore, it is difficult to measure transmission and infection rates. As a consequence, the effect of virus on insects that have been exposed to the disease but have not died from it, is not well understood. Since there is always variability among individuals in the field, if different larval mortality rates or adult fecundities are measured, it is difficult to determine why these differences arose, or how much of the difference may be attributed to virus (e.g., are two values different due to differences in environment or due to a sublethal effect of virus?). The survival rate of the virus and the degree of clumping of 101 the free-living virus near larvae are both critical parameters in the production of cycles, but are almost impossible to determine. These are all aspects of the biology that can be studied further and it is possible that when searching for clarity in the determination of parameter values, additional information about interactions between the virus and the larvae will be discovered. The second possible explanation states that there may be some other factor or factors, aside from virus, that is or are acting to create the cycles. Many of these other possible components, such as plant quality, weather, genetics, or parasites, have been modelled before by other researchers and are discussed in Section 4.2. None of the hypotheses seem entirely plausible for the tent caterpillar case unless one of the factors behaves similarly to the virus in its effect on the larvae (such as reducing fecundity). For example, if a mother experienced bad growing conditions, possibly from weather or food, she may have reduced fecundity or her offspring may have reduced survival. It is more probable that it is one or more elements acting together to produce cycles. For example (although this has not been modelled), it may be that weather and parasites combine to produce oscillations. Few models that have been produced examine the effect of multiple factors, excluding virus, on insect populations. The third explanation is essentially a sub case of the second since it suggests that it may be virus acting with at least one other system component to create the cycles. These various combinations could include possibilities such as weather and virus, parasites and virus, plant quality and virus or even weather, plant quality and virus. It also may be that at high densities the larvae become much more susceptible to a variety of mortality 102 factors, including disease and that such mortality factors would act to decrease survival rates as are seen in the natural populations. Several models have examined the possibility of a multiple component interaction producing cycles (Berrryman, pers. comm.; Hochberg et al., 1990; Wellington, et al., 1975). Berryman maintains that the cycles are caused by the action of the parasites and that the virus just increases the magnitude of the cycles. Hochberg et al. (1990), show that cycles are possible when both a parasite and a pathogen are present in the system. As discussed in Section 3.3.2, Wellington et al. (1975) produced a model that combined effects of parasites, virus and weather. That model could be further tested using more recently collected data. These models provide initial support for the hypothesis that more than one factor creates the cycle. Further modifications of the models are necessary and some support from a natural system is still needed to see if the models are valid for the tent caterpillar system. The results from my model also support this third explanation. Using only virus, it was possible to produce cycles with the right period. The model cycles differed from natural cycles in that the population decline was caused by a virus build-up which led to more exposed individuals and lower average fecundities, leading to a population decline. In the natural system, the average fecundities do not decrease until after the population declines and it is lowered survival rates that seem to influence the population decline. Therefore, if my model is a reasonable representation, some other factor (not included in the model) may be working with the virus to reduce population levels before the virus can reduce the average fecundity. In particular, some component in the natural 103 system must be acting to reduce survival without affecting fecundity. The new model, like many others, is a purely deterministic model, and assumes that all parameter values such as survival rates and basic fecundities remain constant from year to year. In nature, this is obviously not the case. There are always stochastic events occurring which may increase or decrease any of these values, or change some of the relationships. A valid deterministic model should not be sensitive to small changes in parameter values. If the stochastically induced changes are small, many of the results that have been discussed remain similar but instead of all the results being smooth, there is variation around the basic lines. However if these changes are large, the deterministic model no longer produces similar results and is not an adequate representation of the real world. Finally, the use of deterministic models makes it easy to explore various hypotheses and to see the effect of various changes or assumptions. It would be possible to add to and to adapt this model to further examine some of the above explanations. There are always more modifications that can be made to any model and more possibilities that can be explored. One obvious change that could be made to my model is to alter the clumping parameter /3. It is likely that in a natural population, does not remain constant. At low populations, virus will be confined to the areas near the tents and would be low. However, at high populations, tents and caterpillars are abundant and there are more interactions between caterpillars of different tents. The virus will then be less likely to be as clumped near particular tents, and the parameter may be higher. Another modification to the model would be to add some other factor that affects 104 the tent caterpillars, such as parasites, to see if the addition of this factor would cause the population to cycle in the same patterns as are seen nature. Ideally, this factor would act to reduce survival levels, but not fecundity, near the peak of the cycle. It is likely that simply adding further detail to the virus-insect interactions would not have this effect. Finally, since this model is developed specifically for tent caterpillars, it would be interesting to adapt it for other Lepidoptera to see if any of the results could be generalized. The results from the model also suggest key parameters that ideally should be measured in the field. Two of these are related to the free-living virus. If the incidence of free-living virus (i.e., where it is found) could be measured, the clumping parameter could be estimated (see Krebs, 1989 pp81-89). Another key parameter was the survival rate of the free-living virus, both during the larval growing season and during the winter. This parameter is also key in Hochberg et al.'s (1990) model for insectpathogen interactions. Other parameters that it would be useful to measure would be the vertical transmission rate (what proportion of a mother's offspring are infected), and the fecundity of females exposed as larvae. 105 7. LITERATURE CITED Anderson, R.M. 1979. The influence of parasitic infection on the dynamics of host population growth. in Population Dynamics, 20th Symposium. Brit. Ecol. Soc. London, 1978 (eds Anderson, Turner, and Taylor), Blackwell. . 1991. Discussion: The Kermack-McKendrick epidemic threshold theorem. Bull. Math. Bio. 53:3-32. Anderson, R.M. and R.M. May. 1978. Regulation and stability of host-parasite population interactions. I. Regulatory processes. J. Anim. Ecol. 47:219-247. . 1979. Population biology of infectious diseases. Part I. Nature. 280:361-367. ^. 1980. Infectious diseases and population cycles of forest insects. Science 210:658661. . 1981. The population dynamics of microparasites and their invertebrate hosts. Philos. Trans. R. Soc. London. Ser. B. 29:451-524. Andrewartha, H.G. and L.C. Birch. 1954. The Distribution and Abundance of Animals. University of Chicago Press, Chicago. Bailey, N.T.J. 1975. The Mathematical Theory of Infectious Diseases. 2nd ed. Macmillan Pub. Co., Inc. New York. Baltensweiler, W and A. Fischlin. 1988. The larch budmoth in the Alps. in Dynamics of Forest Insect Populations. (ed. A. Berryman) Plenum. Pub. Corp. pp 331-351. Briese, D.T. and J.D. Podgwaite. 1985. Development of viral resistance in insect populations. in Viral Insecticides for Biological Control. (ed. K. Maramorosch and K.E. Sherman) Academic Press, Orlando. pp 361-98. Brown, G.C. 1987. Modeling. in Epizootiology of Insect Diseases. (eds J.R. Fuxa and Y. Tanada) John Wiley & Sons. pp 43-68. Clark, E.C. 1955. Observations on the ecology of a polyhedrosis of the Great Basin tent caterpillar Malacosoma fragilis. Ecology 36:373-376. ^. 1956. Survival and transmission of a virus causing polyhedrosis in Malacosoma fragile. Ecology 37:728-732. ^. 1958. Ecology of the polyhedrosis of tent caterpillars. Ecology 39:132-139. Daniel, C. 1990. Climate and outbreaks of the forest tent caterpillar in Ontario. M.Sc. thesis, University of British Columbia. Diekmann, 0. and M. Kretzschmar. 1991. Patterns in the effects of infectious diseases on population growth. J. Math. Bio. 29:539-570. Dwyer, G. 1991. The roles of density, stage, and patchiness in the tranrnission of an insect virus. Ecology 72:559-574. Edelstein-Keshet, L. 1986. Mathematical Biology. Random House Inc. New York. Edelstein-Keshet, L. and M.D. Rausher. 1989. The effects of inducible plant defenses on herbivore populations: 1. Mobile herbivores in continuous time. Am. Nat. 133:787809. 106 Entwistle, P.F. 1983a. Control of insects by virus diseases. Biocontrol News and Information. 4:202-228. Commonwealth Institute of Biological Control. ^. 1983b Viruses for insect pest control. Span. 26. Finerty, J.P. 1980. The Population Ecology of Cycles in Small Mammals. Yale University Press, New Haven. Hassell, M.P. 1980. Foraging strategies, population models and biological control: a case study. J. Anim. Ecol. 49:603-628. Hassell, M.P. and R.M. May. 1989. The population biology of host-parasite, host parasitoid associations. in Perspectives in Ecological Theory. Princeton University Press, Princeton. pp 319-347. Haukioja, E. 1980. On the role of plant defenses in the fluctuation of herbivore populations. Oikos 35:202-213 Hethcote, H.W. and P. van den Driessche. 1991. Some epidemiological models with nonlinear incidence. J. Math. Bio. 29:271-287 Hochberg, M.E. 1989. The potential role of pathogens in biological control. Nature 337:262-265. Hochberg, M.E., M.P. Hassell and R.M. May. 1990. The dynamics of host-parasitoidpathogen interactions. Am. Nat. 135:74-94. Hunt, F. 1982. Regulation of population cycles by genetic feedback: existence of periodic solutions of a mathematical model. J. Math. Bio. 13:271-282. . 1983. A mathematical analysis of the Chitty Hypothesis. in Population Biology. (eds H.I. Freedman and C. Strobeck) Lecture Notes in Biomathematics. SpringerVerlag. 52:33-40. Hutchinson, G.E. 1978. An Introduction to Population Ecology. Yale University Press, New Haven. Kermack, W.O. and A.G. McKendrick. 1927. Contributions to the mathematical theory of epidemics-I. Proc. Roy. Soc. 115A:700-721. Kingsland, S.E. 1985. Modelling Nature. University of Chicago Press, Chicago. Krebs, C.J. 1989. Ecological Modelling. Harper & Row, New York. pp 81-89. Liu, W.M., H.W. Hethcote and S.A. Levin. 1987. Dynamical behavior of epidemiological models with nonlinear incidence rates. J. Math. Bio. 25:359-380. Liu, W.M., S.A. Levin and Y. Iwasa. 1986. Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models. J. Math. Bio. 23:187-204. Malthus, T.R. 1798. An Essay on the Principle of Population. in An Essay on the Principle of Population (ed A. Flew) Penguin, England, 1970. Martinat, P.J. 1987. The role of climatic variation and weather in forest insect outbreaks. in Insect Outbreaks. (eds P.Barbosa and J.C. Schultz) Academic Press, Inc., San Diego May, R.M. 1973. Stability and Complexity in Model Ecosystems. Princeton University Press, Princeton. 107 ^. 1978. Host-parasitoid systems in patchy environments: a phenomenological model. J. Anim. Ecol. 47:833-843. ^. 1981. Models for two interacting populations. in Theoretical Ecology 2nd ed (ed. R.M. May) Blackwell Scientific Pub. Oxford. May, R.M. and R.M. Anderson. 1978. Regulation and stability of host-parasite population interactions. II. Destabilization processes. J. Anim. Ecol. 47:249-267. . 1979. Population biology of infectious diseases. Part II. Nature. 280:455-461. Murray, J.D. 1989. Mathematical Biology Springer-Verlag, Berlin. Myers, J.H. 1981. Interactions between western tent caterpillars and wild rose: A test of some general plant herbivore hypotheses. J. Anim. Ecol. 50:11-25. . The induced defense hypothesis: Does it apply to the poplation dynamics of insects? in Chemical Mediation of Coevolution (ed K.C. Spenser) Academic Press Inc. San Diego. pp 345-366. ^. 1988b. Can a general hypothesis explain population cycles of forest lepidoptera? Advances in Ecological Research 18:179-242. . 1990. Population cycles of Western tent caterpillars: experimental introductions and synchrony of fluctuations. Ecology 71:986-993. Onstad, D.W. and R.I. Carruthers. 1990. Epidemiological models of insect diseases. Ann. Rev. Entomol. 35:399-419. Pielou, E.C. 1969. An Introduction to Mathematical Ecology. John Wiley & Sons, New York. Press, W.H., B.P. Flannery, S.A. Teukolsky and W.T. Vettering. 1988. Numerical Recipes in C. Cambridge Press, Cambridge. pp 387-393. Regniere, J. 1984. Vertical transmission of diseases and population dynamics of insects with discrete generations: a model. J. Theor. Bio. 107:287-301. Rhoades, D.F. 1983. Herbivore population dynamics and plant chemistry. in Variable Plants and Herbivores In Natural and Managed Systems. (eds R.F. Denno and M.S. McClure) Academic Press, New York. pp 155-220. Smith, J.D. and R.A. Goyer. 1986. Population fluctuations and causes of mortality for the forest tent caterpillar, Malacosoma disstria (Lepidoptera: Lasiocampidae), on three different sites in southern Louisiana. Environ. Entomol. 15:1184-1188. Stairs, G.R. 1965. Artificial initiation of virus epizootics in forest tent caterpillars populations. Can. Ent. 97:1059-1062. Starfield, A.M. and A.L. Bleloch. 1986. Building Models for Conservation and Wildlife Management. MacMillan Pub. Co. New York. Stenseth, N.C. 1981. On Chitty's theory for fluctuating populations: the importance of genetic polymorphism in the generation of regular density cycles. J. Theor. Bio. 90:9-36. Thompson, W.A., I.B. Vertinsky and W.G. Wellington. 1979. The dynamics of outbreaks: further simulation experiments with the Western tent caterpillar. Res. Pop. Ecol. 20:188-200 108 Vezina, A. and R.M. Peterman. 1985. Tests of the role of a nuclear polyhedrosis virus in the population dynamics of its host, douglas-fir tussock moth, Orgyia pseudotsugata (Lepidoptera: Lymantriidae). Oecologia 67:260-266. Volterra, V. 1926. Variazioni e fluttuazioni del numero d'individui in specie animali conviventi. Mem. Acad. Lincei. 2:31-113. translated in Chapman, R.N. 1931. Animal Ecology, Appendix. pp 409-448. Wellington, W.G. 1962. Population quatlity and the maintenance of nuclear polyhedrosis between outbreaks of Malacosoma pluviale (Dyar). J. Insect Pathol. 4:285-305 Wellington, W.G., P.J. Cameron, W.A. Thompson, I.B. Vertinsky and A.S. Landsberg. 1975. A stochastic model for assessing the effects of external and internal heterogeneity on an insect population. Res. Popul. Ecol. 17:1-28. Young, S.Y. 1990. Effect of nuclear polyhedrosis virus infection in Spodoptera ornithogalli larvae on post larval stages and dissemination by adults. J. Invert. Path. 55:69-75. 109 ^ ^ Appendix 1: ANALYSIS OF MODELS A1.1 Introduction Many differential models and difference models can be analyzed easily by following the same basic steps in both types of systems: 1. Find the null dines and steady states of the equations. 2. Linearize about each steady state and find the Jacobian and eigenvalues of the Jacobian for each. 3. From the eigenvalues, determine the conditions under which each steady state is st able. In this section, I analyze two of the models discussed in Chapter 2, the Lotka-Volterra equations and the Nicholson-Bailey model. All analyses are based on Edelstein-Keshet (1986). A1.2 Lotka-Volterra The Lotka-Volterra model (equations 3.6) can be analyzed using the steps mentioned above. 1. Find and graph the null dines and steady states of the equations: • Null dines: These are found by setting the left-hand side of each of the equations to zero. So equation 3.6a becomes: 0 rN — bPN (A1.1) N= 0 or P r lb. (A1.2) 0 cPN — dP (A1.3) which implies: Equation 3.6b gives: which implies P= 0 or N dlc^(A1.4) • Steady states: These are found wherever a null dine from each of equations A1.2 and A1.4 cross (but not where two null dines from the same equation cross). For the Lotka-Volterra equations, there are two such equilibria: trivial:^(N,P) = (0,0) non-trivial: (N ,P)^(d1c,r lb) 2. Find the Jacobian and eigenvalues for each of the steady states: • at (0,0): Jacobian: rr^0 1 (A1.5) [0 —d • The eigenvalues are: r and —d. Since one eigenvalue is greater than zero and one is less than zero, (0, 0) corresponds to a saddle point. • at (d1c,r1b): Jacobian: [ 0^—bd/c1 (A1.6) [^_I' cr lb^0 110 ^ Here the eigenvalues are: +i\/rd. Notice that the real part of the eigenvalues is 0, so the equilibrium is neutrally stable with oscillations about it. A1.3 Nicholson-Bailey Next, the Nicholson-Bailey model (equations 3.11) for insect-parasitiod systems will be analyzed following the steps above. 1. Null Clines and Steady States The null dines of difference equations can be found, by definition, when the population is not changing or Nt+i Nt . Thus, solving for the isoclines in equation 3.11a leads to the conditions: N=0 or P ln r a (A1.7) while equation 3.11b gives the isocline: N P c(1 — e — aPt) (A1.8) Thus, there are two steady states: and (N, P) (0,0) (A1.9) — D\^rinr^lnr \ \Iv^)^ac(r — 1)^a ). (A1.10) The first, equation A1.9, is the trivial equilibrium while equation A1.10 is the non-trivial eqilibriu m. 2. The Jacobian The general matrix is e— a15^—aN e — a l ( A 1 .11 ) c(1 — e — a 1-5 ) cal Tre — a15 [ When evaluated at the trivial equilibrium (A1.9), matrix A1.11 becomes: ^ (A1.12) [Or 01 Clearly the eigenvalues of this matrix are r and 0. In difference equations, a steady state is stable if all the eigenvalues of the Jacobioan evaluated at that steady state are between -1 and 1. Thus, in this case, the trivial equilibrium is stable if r < 1. This means that the population will die out if r < 1. The second case is more complicated. When the Jacobian is evaluated at the non-trivial equilibrium (A1.10), matrix A1.11 becomes: r In r c(r —1) In r (r-1) 111 • (A1.13) ^ So, let i3 =1+ ln r ^ (r ^1) (A1.14) and, r In r 7= ^ (r — 1)^•^ (A1.15) Now the stability condition given by Edelstein-Keshet (1986) for difference equations, 2 < 1 + -y < 1 /31, can be checked. It is easy to show (Edelstein-Keshet, p82, 1986) that the first half, 2 < 1 + -y, is true only if r < 1 while the second part, 1+7 < 101, holds only if r > 1. Therefore, the non-trivial equilibrium is never stable (and if the trivial equilibrium is also unstable, the population will show increasing, unstable, oscillations). 112 Appendix 2: ANALYSIS OF THE MODEL A2.1 Introduction In this appendix, I provide the details for the analysis of the new model (equations 5.1). In the first part, I derive an approximate reproductive rate for the virus. In the second part, since the techniques for analyzing systems of equations are already mentioned in Appendix 1, the equations are simplified and are put into a form for which these techniques can be used, and the stability of the system determined. The third section describes some of the techniques that I used to attempt to numerically determine the stability of the system and the possible existence of cycles. A2.2 Virus reproductive rate The first thing that can be calculated is the reproductive rate for the virus. After substituting equation 5.1c into equation 5.1d, the basic equation for the virus becomes: Vy ,t+1 = Cr y Vy ,t^(rt^wt )Pt ( 0 -3 Sy,t^(7e Ey ,t)^(A2.1) where, as before, t^ pt = 1 — (1 + V Ow t Hy , t (5.3) Equation A2.1 is non-linear, due to the form of P t . This makes any analysis difficult. Therefore, to simplify the analysis, pt can be linearized about V = 0 to get: Pt Vt tot Hy , t (A2.2) Substituting A2.2 into A2.1, 1 Vt+1^7v V, ,t^ Wt H y,t Wt)(C r S y,t^aeEy ,t) (A2.3) and rearranging the terms gives: (A2.4) Vy,t+1^RtVy,t e (ytq Wt)(6rsSy,t UeEy,t) R =^ WtHy,t (A2.5) where, Hy,t = Sy ,t^Ey ,t^,t^ (A2.6) Since Pt was linearized about V = 0, A2.4 is valid only when the amount of virus is small. In theory, it would be possible to find other expressions for the reproductive rate of the virus that would be accurate at different amounts of virus. This would be done by linearlizing p about some other steady state (where none of the populations are changing). In practice, this is not possible analytically without making other major simplifications to the model. 113 A2.3 Equilibrium Analysis To determine analytically the behaviour of the system near equilibrium the model must first be simplified. Since it was the term p t , that was creating the difficulties in analysis, I simplified the model by assuming that p remained constant. The analysis can then be easily performed on the adult stage of the model, by combining all tents. Therefore, merging equations 5.8-5.10, the total numbers of organisms are: Sy +i,o = (f8Sy s+ 1 — i)feEy,6)^(A2.7a) , 2 dy 2 ( r 4+1,0 = —Kj_G y , 6^(A2.7b) Ty+1 = Y'6 + Ey ,6^ 2 (A2.7c) and, as before, d is given by dy^1^ a + Ty (5.7) It follows from the assumptions in the model that the number of adults in year y is dependent on the number of young starting the same year, and, due to the simplification, that the dependence is linear. Thus, ^ Sy,6^XSy,0 ^ Ey,6 = ZSy,0 (A2.8a) (A2.8b) where x and z are independent of the state variables. (The exact formulas for x and z are not important for the stability calculations but are shown in equations A2.19 and A2.22 and discussed further at that place.) Substituting equations A2.8 into equations A2.7 and dropping the t subscript since t = 0 on both sides of the equation, gives: dy Sy+1 = 2 (Lx + (1 — — n)fe z)Sy^(A2.9a) dy KfezSy (A2.9b) 2 dy Ty+1 = — (x z)Sy .^ (A2.9c) 2 Using the standard techniques demonstrated in Appendix 1.3 (step 1), the equilibrium values for this system (equations A2.9) can be found at: /y+1 = — (S, I, T) = (0, 0, 0)^ (A2.10) and at, — 2) 2(x z) S= ^ /= T Kfe za(0 2) 2(x + z) — alk 2 a where, 114 (A2.11) f8 x + (1 — ti)fe z.^ (A2.12) The term represents the amount of reproduction of the organisms. The second equilibrium, given by equations A2.11, does not exist if density-dependence was not present. To see this, the density-dependence term dy in equations A2.9 would be replaced with the constant 1. If the equilibrium of A2.9 was calculated as before, the only possible steady state would be the trivial one given in equation A2.10. The Jacobian for the system (equations A2.9) is J cr 2(ce-I-T) atcf, z 2(cr+T) cx(x+z) 2(a+T) cP 2(a+T) 2 aKf2z 2(a-FT) 2 a(x+ z) 0 2(a-FT) 2 3 (A2.13) ^_ Substituting the trivial equilibrium (equation A2.10) into matrix A2.13, and calculating the eigenvalues of the resulting matrix, it can be seen that the eigenvalues are 0, 0, and Ik/2. Similarily, the eigenvalues of the Jacobian evaluated at the non-trivial equilibrium (A2.12) are 0, 0 and 2/0. In Appendix 1.3, it was noted that if Ai is an eigenvalue, the condition for stability of an equilibrium is that 'Ail < 1 for all i. Clearly, 0 < 1, so the stability of each eigenvalue depends on the value of 0. There are two possibilities: < 2 A2.10 stable (A2.14) > 2 A2.11 stable (A2.15) Cycles are not possible in either case. As 7/, approaches 2 from above, the non-trivial equilibrium point appraoches the trivial one, until, when = 2, the non-trivial equilibrium equals the trivial one. For cycles to be possible, either both equilibria must be saddle-nodes, in which case a saddle node connector cycle may be present, or at least two of the eigenvalues of a steady state must be 'Ail = 1, in which case a Hopf bifurcation leading to limit cycles may be present (Guckenheimer and Holmes, 1983). The virus V can also be included in the analysis. Following the same simplification process that was done above to change equations A2.8 into A2.10, the amount of virus from year to year is V +1 = a6Vy coS Y ZSy (A2.16) where w is some combination of parameters only, and is independent of all the variables. (As with x and z, the exact form of w is not important to the stability calculations.) When the virus (in the form of equation A2.16) is included in the analysis, the only change in the results is that the resulting Jacobian has another eigenvalue that equals . Since this is always less than one, the stability conditions given in equations A2.14 and A2.15 do not change. The parameter p was assumed to be constant for the purposes of doing this analysis. However, the value of p affects the results. To see what happens as p is varied, the forms 115 for x and z must be calculated. They are found by solving the difference equations for the within-year dynamics. (5.1a) Syt+i = Q 3 (1 — 13)S y^ becomes Sy , t = (u 3 (1 — p)) t Sy , o . ^ (A2.17) Equation 5.1a could not be solved explicitly if p were dependent on any of the other model variables. Since there were six time-steps for this part of the model, let t = 6. Also, the survival from the pupal stage must be included. So, equation A2.17 becomes Sy,6 = (u s (1 — p)) 6 up Sy ,0 ^ (A2.18) This is now in the form of equation A2.9a so clearly x = (o- 3 (1 — p)) 6 gyp .^ (A2.19) To calculate z, the same procedure is followed, using equation A2.17, and remembering that no organisms are born exposed (Ey , 0 0). Equation 5.1b can be rearranged to give Ey , t+i = p(1 — q)usSy ,t (1 — Pg)creEy ,t • (A2.20) This then becomes Ey ,t^ Cris (1 — py \ 11 pq)u e 13(1 —^ QS (1 — 1)) — (1 — Sy (A2.21) After allowing for pupal mortality, and noticing that equation A2.21 is in the same form as equation A2.9b, it is clear that z =^(1 — p) 6 —(1 — pq) 6 0 ^) '78( 1 — — (1 — pq ) cr e 6 P( 1 0 0 sup. - (A2.22) Upon examination of x and z given in equations A2.19 and A2.22 and shown in Figure 18, it can be seen that x > z for small p (since all parameters involved in equations A2.19 and A2.22 are < 1). As p approaches one (its upper bound), x goes to 0 and z becomes small. Since 7b depends on x and z, it will also be small, less than 2 (unless fe z the birth-rate of exposed mothers is very high). As p becomes small, lk will grow. In general, there is a critical value, pc which corresponds to the case 7b = 2. If p > pc then the population will die out while if p < p c the non-trivial equilibrium is stable (see Figure 17). There is no explicit formula for pc . Thus, the value of p that is chosen affects the level at which the population stabilizes, and even whether or not the population will survive. Obviously, the other parameters involved in will also help determine what value has. Similar analysis could be done for the other parameter values, but since the results from these other analysis would also be dependent on p having a constant value, and would have less meaning in the full model, the analysis will not be presented. 116 0.0^ 0.2^ 0.4^0.6 ^ 0.8 ^ 1.0 Infection probability (p) Figure 18 Parameters x and z versus the infection probability, p. Parameters are: a, = .93, a, .87 and q = .5 (the base values given in Table 5). 117 A2.4 Stability: numerical analysis Since the general analysis performed in section A2.3 is necessarily inaccurate due to the simplification that was made, I also attempted some numerical analysis of the full model. As with the previous analysis, if a steady state can be found, the eigenvalues of the Jacobian of the matrix will describe the behaviour of the populations near the steady state. One way for cycles to occur is to have a Hopf bifurcation. The first criterion for finding such a bifurcation is that for some set of parameters there must be at least one pair of complex conjugate eigenvalues with lA i I = 1 (Guckenheimer and Holmes, 1983). The bifurcation will arise at this point. In order to find the eigenvalues, it was first necessary to find a steady state and to calculate the Jacobian at this point. Initially, I found a steady state by simulating the model until the population stabilized (the population size no longer varied). I then took the partial derivatives of the model at this steady state and created a Jacobian matrix. Next, I used an eigenvalue finding program (Press et. al., 1988, pp 387-393) to find the eigenvalues of the Jacobian. In practice, there was a major problem with computer generated error. The method that I used to find the eigenvalues stated that the error would be on the order of the euclidean norm of the Jacobian. This was quite large (sometimes up to 80%). In order to find a Hopf bifurcation, it is important to know when the eigenvalues change from less than one to greater than one. Having large error makes it impossible to find this point. The was evident from the simulations that I did. For example, if the steady state is stable, the magnitudes of the eigenvalues should all be less than one. The method I initially used to find the steady state (by running the model on a computer until the population size no longer changed) guaranteed that the steady state was a stable one. Simulations perturbing the population away from the steady state confirmed this since the population always returned within three generations. However, the eigenvalues that were found were not all less than one. If the true value of an eigenvalue was near one, even a small amount of error could be enough to make it larger than one, or to create a complex value when the true value was real. This error made this technique useless and no further numerical analysis was performed using this or any other method. 118
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Models for tent caterpillar-virus interactions
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Models for tent caterpillar-virus interactions Beukema, Sarah Jenelle 1992
pdf
Page Metadata
Item Metadata
Title | Models for tent caterpillar-virus interactions |
Creator |
Beukema, Sarah Jenelle |
Date Issued | 1992 |
Description | Many species of forest Lepidoptera show eight to twelve year population cycles which may involve viral disease. To examine possible interactions between viral disease and population cycles of forest Lepidoptera I explored some models for insect-virus dynamics. All of the models produced population oscillations in their original form. However, after they were modified to conform more closely to the tent caterpillar system, none of the models produced realistic cycles. I then developed a new model specifically for the tent caterpillar system that included important features such as: reduced fecundity of individuals that had been exposed to the disease, transmission of the disease from mother to progeny, a free-living infective stage of the virus, and a horizontal transmission rate that varied with the larval stage, the number of individuals, and amount of virus present. Cyclic dynamics resulted from some simulations. The parameters producing the cycles were similar to actual data. However, unlike natural populations of tent caterpillars, in the simulated population the average fecundity decreased before the population started to decline and survival decreased at approximately the same time as the population. Important further research in the field should include investigation of the distribution and survival of free-living virus and factors that would reduce caterpillar survival at peak populations but not affect fecundity. |
Extent | 4620643 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-12-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0086808 |
URI | http://hdl.handle.net/2429/3276 |
Degree |
Master of Science - MSc |
Program |
Zoology |
Affiliation |
Science, Faculty of Zoology, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1992-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1992_spring_beukema_sarah_janelle.pdf [ 4.41MB ]
- Metadata
- JSON: 831-1.0086808.json
- JSON-LD: 831-1.0086808-ld.json
- RDF/XML (Pretty): 831-1.0086808-rdf.xml
- RDF/JSON: 831-1.0086808-rdf.json
- Turtle: 831-1.0086808-turtle.txt
- N-Triples: 831-1.0086808-rdf-ntriples.txt
- Original Record: 831-1.0086808-source.json
- Full Text
- 831-1.0086808-fulltext.txt
- Citation
- 831-1.0086808.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0086808/manifest