@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Arts and Sciences, Irving K. Barber School of (Okanagan)"@en, "Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCO"@en ; dcterms:creator "Strohm, Shaun"@en ; dcterms:issued "2013-08-04T14:39:14Z"@en, "2013"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """In this thesis, we use a reaction-diffusion equation with chemotaxis to model the interaction between Mountain Pine Beetle (MPB, Dendroctonius ponderosae), Mountain Pine Beetle pheromones, and susceptible trees. The goal is to understand how the movement and attack of MPB is affected by management activities. We investigate the spatial pattern formation of attack clusters in a system for Mountain Pine Beetle. Mathematical analysis is utilized to discover the spacing between beetle attacks on the susceptible landscape. The model predictions are verified by analyzing aerial detection survey data of Mountain Pine Beetle attack from the Sawtooth National Recreation Area. We find that the distance between Mountain Pine Beetle attack clusters predicted by our model and observed in the Sawtooth National Recreation Area are the same. These results clarify the spatial mechanisms controlling the transition from incipient to epidemic populations and may eventually lead to control measures which protect forests from MPB outbreaks. Our next avenue of investigation is using an experimental study and theoretical work to help understand the effects of habitat fragmentation on the movement of the MPB. The experimental study consists of trap catch data for MPB in different domains of fragmented habitat. We simulate the experimental system using our mathematical model, testing different hypothesis on initial position of MPB emergence and diffusion speed. Our study provides support for the hypothesis that MPB may move faster in harvested landscapes, and that MPB emerge uniformly over the landscape. Finally, we use a multi-year spatially explicit model to test the effectiveness of the management strategies of baiting and tree-removal and prescribed burning. We find that baiting and tree-removal is successful at reducing MPB density and forest impact, as long as MPB emergence densities are not too small. We predict that tree removal without baiting can be more successful than combined baiting and tree removal if the searched area has a large density of MPB. Finally, analysis of our model indicates that prescribed burning can be more effective than clearcutting given certain assumptions about the reproductive output and attractiveness of burned trees."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/44798?expand=metadata"@en ; skos:note "Dispersal of Mountain Pine Beetle andImpacts of ManagementbyShaun StrohmB.Sc. Hons., The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe College of Graduate Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)July 2013c? Shaun Strohm, 2013AbstractIn this thesis, we use a reaction-diffusion equation with chemotaxis tomodel the interaction between Mountain Pine Beetle (MPB, Dendroctoniusponderosae), Mountain Pine Beetle pheromones, and susceptible trees. Thegoal is to understand how the movement and attack of MPB is affected bymanagement activities.We investigate the spatial pattern formation of attack clusters in a sys-tem for Mountain Pine Beetle. Mathematical analysis is utilized to discoverthe spacing between beetle attacks on the susceptible landscape. The modelpredictions are verified by analyzing aerial detection survey data of Moun-tain Pine Beetle Attack from the Sawtooth National Recreation Area. Wefind that the distance between Mountain Pine Beetle attack clusters pre-dicted by our model and observed in the Sawtooth National RecreationArea are the same. These results clarify the spatial mechanisms controllingthe transition from incipient to epidemic populations and may eventuallylead to control measures which protect forests from MPB outbreaks.Our next avenue of investigation is using an experimental study andtheoretical work to help understand the effects of habitat fragmentation onthe movement of the MPB. The experimental study consists of trap catchdata for MPB in different domains of fragmented habitat. We simulatethe experimental system using our mathematical model, testing differenthypothesis on initial position of MPB emergence and diffusion speed. Ourstudy provides support for the hypothesis that MPB may move faster inharvested landscapes, and that MPB emerge uniformly over the landscape.Finally, we use a multi-year spatially explicit model to test the effective-ness of the management strategies of baiting and tree-removal and prescribedburning. We find that baiting and tree-removal is successful at reducingiiAbstractMPB density and forest impact, as long as MPB emergence densities arenot too small. We predict that tree removal without baiting can be moresuccessful than combined baiting and tree removal if the searched area has alarge density of MPB. Finally, analysis of our model indicates that prescribedburning can be more effective than clearcutting given certain assumptionsabout the reproductive output and attractiveness of burned trees.iiiPreface? The mathematical model was constructed by Shaun Strohm. To en-sure biological relevance we consulted with MPB biologist Mary Reid(University of Calgary) and Kurt Trzcinski (Department of Fisheriesand Oceans, Dartmouth, NS).? Chapter 2 is a study of the attack patterns formed by MPB at incipientepidemic levels. Mountain Pine Beetle attack data from the SawtoothNational recreation area were provided by Jim Powell (Utah StateUniversity). More specifically, the data are USDA Forest Service aerialdetection survey; full details are provided in Crabb et al. [CPB12].Analysis of the data and mathematical analysis of the model weredone by Shaun Strohm. This work has been accepted in the Bulletinof Mathematical Biology and is currently in press. Jim Powell, ourco-author, provided valuable insight into this project.? Chapter 3 is an elucidation of Mountain Pine Beetle trap catch datathrough a simulation study of the mathematical model. The data wereprovided by Mary Reid. All simulation work and analysis were doneby Shaun Strohm. Mary Reid was consulted regularly on this project.? Chapter 4 is a study of the effects of different management activitieson Mountain Pine Beetle Movement. All mathematical analysis andsimulation work was done by Shaun Strohm. To maintain relevanceto management and MPB biology we consulted with Mary Reid anda representative from Parks Canada in Banff, Jane Park.? The simulation model was written using a Fortran code base developedoriginally by Rebecca Tyson, and heavily modified by Shaun StrohmivPrefaceto solve the Mountain Pine Beetle model equations. The model usesCLAWpack, a program which is capable of solving very stiff hyperbolicsystems with high accuracy. This program is developed by RandallLeVeque [LT96, LeV97].? All the research subprojects above were completed under the super-vision of Rebecca Tyson. My committee members Dan Coombs andSylvie Desjardins also provided feedback and helpful suggestions.? All thesis text is original, unpublished, and was written exclusively byShaun Strohm.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 1Chapter 2: Pattern Formation in a Model for Mountain PineBeetle Dispersal: Linking Model Predictions toData . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . 72.2.1 Non-dimensionalization . . . . . . . . . . . . . . . . . 102.3 Model Pattern Formation . . . . . . . . . . . . . . . . . . . . 122.3.1 Spatially Uniform Steady State . . . . . . . . . . . . . 132.3.2 Linear Analysis . . . . . . . . . . . . . . . . . . . . . . 142.3.3 Analysis of the Dispersion Relation . . . . . . . . . . . 162.3.4 Numerical Analysis of the Dispersion Relation . . . . 192.4 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.5 Management Implications . . . . . . . . . . . . . . . . . . . . 25viTABLE OF CONTENTS2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Chapter 3: Edge Effects on Mountain Pine Beetle Movement 323.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.2 Experimental . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.3 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . 353.3.1 Non-dimensionalization of Model . . . . . . . . . . . . 393.4 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . 413.5 Numerical Method for Simulations . . . . . . . . . . . . . . . 443.6 Model Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.7 Connections Between Theoretical and Empirical Results . . . 633.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Chapter 4: Impacts of Management on Mountain Pine Bee-tle Movement . . . . . . . . . . . . . . . . . . . . . . 684.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.2 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . 714.2.1 Estimation of parameters . . . . . . . . . . . . . . . . 734.2.2 Non-dimensionalization . . . . . . . . . . . . . . . . . 744.3 Effects of Management . . . . . . . . . . . . . . . . . . . . . . 774.3.1 Baiting and Tree Removal . . . . . . . . . . . . . . . . 774.3.2 Prescribed Burning . . . . . . . . . . . . . . . . . . . . 904.4 Characteristics of MPB spread . . . . . . . . . . . . . . . . . 984.4.1 Interaction Between Source Trees . . . . . . . . . . . . 984.4.2 Clearcut Length . . . . . . . . . . . . . . . . . . . . . 984.4.3 Source Tree Mapping . . . . . . . . . . . . . . . . . . 1014.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 109Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125Appendix A: FFT Analysis and Parameter Sensitivity . . . . . . . 125viiTABLE OF CONTENTSA.1 FFT analysis of MPB data . . . . . . . . . . . . . . . . . . . 125A.2 Parameter Sensitivity . . . . . . . . . . . . . . . . . . . . . . 125Appendix B: Burn Simulations . . . . . . . . . . . . . . . . . . . . 131viiiList of TablesTable 2.1 Parameter Values for (2.1) . . . . . . . . . . . . . . . . 9Table 2.2 Parameter values for (2.4) . . . . . . . . . . . . . . . . . 11Table 3.1 Table of parameter values for the dimensional model (3.1). 38Table 3.2 Table of parameter values for the non-dimensional model(3.3). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41Table 3.3 Order of accuracy for nesting and flying MPB solutionswith decreasing space step size . . . . . . . . . . . . . . 48Table 4.1 Table of parameter (Par) values for the dimensionalmodel (4.1) and (4.2). . . . . . . . . . . . . . . . . . . . 75Table 4.2 Table of parameter values for the non-dimensional model(4.4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Table A.1 Parameter Sensitivity Table . . . . . . . . . . . . . . . . 126ixList of FiguresFigure 2.1 Numerical and Analytical Contour Plots against Wavenum-ber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21Figure 2.2 Average dominant wavenumber over years 1991-2009. 24Figure 2.3 Histogram of dominant wavenumbers . . . . . . . . . 26Figure 2.4 Management Implications Based on Population Density 27Figure 3.1 This plot shows a heterogeneous landscape with 4 baitsites . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 3.2 Convergence of flying and nesting MPB with decreas-ing space step size . . . . . . . . . . . . . . . . . . . . 47Figure 3.3 Log-Log error plots of abs(error) against space step size 47Figure 3.4 Spatial population densities of flying and nesting MPBin a landscape that has a clearcut harvested region anddense forest initial conditions. . . . . . . . . . . . . . . 51Figure 3.5 Spatial population densities of flying and nesting MPBin a landscape that has a thinned harvested region anddense forest initial conditions. . . . . . . . . . . . . . . 52Figure 3.6 Average trap catches in 4 bait trap catches with MPBstarting in dense forest. . . . . . . . . . . . . . . . . . . 54Figure 3.7 Bait pheromone profile at start of simulation . . . . . 55Figure 3.8 Spatial population densities of flying and nesting MPBin a landscape that has a clearcut harvested region anduniform initial conditions. . . . . . . . . . . . . . . . . 58xLIST OF FIGURESFigure 3.9 Spatial population densities of flying and nesting MPBin a landscape that has a thinned harvested region anduniform initial conditions. . . . . . . . . . . . . . . . . 59Figure 3.10 MPB pheromone profile and chemotactic response be-fore significant aggregation . . . . . . . . . . . . . . . . 60Figure 3.11 Average trap catches in 4 bait trap catches with MPBstarting uniformly on the landscape. . . . . . . . . . . 62Figure 4.1 Simple grid representations of area attacked in moni-toring and management zones. . . . . . . . . . . . . . . 79Figure 4.2 Simulation setup of source trees, bait site, and removal 80Figure 4.3 2 year simulations of nesting MPB density with man-agement activities. . . . . . . . . . . . . . . . . . . . . 81Figure 4.4 Length of forest consumed and MPB population sizewith management activities. . . . . . . . . . . . . . . . 82Figure 4.5 MPB population with variable search width and searcheffectiveness. . . . . . . . . . . . . . . . . . . . . . . . . 85Figure 4.6 Length of forest consumed and total MPB populationwith management for 5 tree source tree simulations. . . 87Figure 4.7 Length of forest consumed and total MPB populationwith consecutive management for 5 tree source treesimulations. . . . . . . . . . . . . . . . . . . . . . . . . 89Figure 4.8 Simulation setup for prescribed burn . . . . . . . . . . 92Figure 4.9 Fraction of susceptible trees attacked in regions burnedby a prescribed fire. . . . . . . . . . . . . . . . . . . . . 93Figure 4.10 Fraction of susceptibles consumed with varying dis-tances between source trees. . . . . . . . . . . . . . . . 99Figure 4.11 Clearcut distance needed for different densities of sourcetrees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100Figure 4.12 Number of MPB source trees mapped from one yearto the next. . . . . . . . . . . . . . . . . . . . . . . . . 102xiLIST OF FIGURESFigure A.1 Choice of Region and Example Power Spectral Densityfor 2007 . . . . . . . . . . . . . . . . . . . . . . . . . . 126Figure A.2 Power spectral densities for 1991-1996 . . . . . . . . . 127Figure A.3 Power spectral densities for 1997-2002 . . . . . . . . . 128Figure A.4 Power spectral densities for 2003-2009 . . . . . . . . . 129Figure A.5 Size and Number of Regions from 1991-2009 . . . . . 130Figure B.1 Nesting MPB density in space for prescribed burnsimulations with undamaged emergence area. . . . . . 132Figure B.2 Nesting MPB density in space for prescribed burnsimulations with partially damaged emergence area. . . 133Figure B.3 Nesting MPB density in space for prescribed burnsimulations with scorched emergence area. . . . . . . . 134xiiAcknowledgementsI would like to acknowledge the support of my supervisor, RebeccaTyson. Her amazing guidance and support allowed this work to be possible.Collaborators Mary Reid and Jim Powell provided very helpful discussionsand direction for the work. Other helpful ideas came from discussions withDan Coombs, Sylvie Desjardins, Jane Park, and Kurt Trzcinski. Critical tothe numerical analysis was Clawpack, developed by Randall LeVeque. Com-puting assistance for HPC was provided by Steve Cundy and Wade Klaver.I would like to also acknowledge helpful discussions from Alan Law and themembers of the Tyson Lab group, past and present, including Ben Wilson,Carly Rozins, Meghan Dutot, Garrett Culos, Alex Blaessle, Chris Coutts,and Katrina Williams. Finally, I must acknowledge the amazing support ofmy wife, Elizabeth, as well as my family (especially my parents, Chris andDianne) and friends. I would also like to acknowledge my funding sourcesNSERC, PIMS IGTC, MITACS, and UBC Okanagan.xiiiDedicationI would like to dedicate this thesis to my wife, Elizabeth Strohm, andmy parents, Chris and Dianne Strohm.xivChapter 1IntroductionThe Mountain Pine Beetle (MPB, Dendroctonius ponderosae) is anaggressive bark beetle that has had a major economic impact on the WesternCanadian and United States forestry industries. At endemic (low popula-tion) levels this beetle is a minor pest, killing trees weakened by drought, rootrot, or lightning strikes [PKW+00]. In contrast, at epidemic levels, thou-sands of acres of healthy, vigorous trees can be rapidly eradicated [PKW+00].The MPB preferentially uses lodgepole pine as its host [SW06]. Whatdistinguishes MPB from other tree pests is that it must kill the host tosuccessfully reproduce. The beetles land on the host tree, then burrowinto the bark to reach the nutrient-rich phloegm, where they dig verticalgalleries. The eggs laid in the galleries produce larvae which feed on thephloegm, cutting off the nutrient pathways of the tree. If consumption ofthe phloegm tissue is sufficiently extensive, the tree dies and develops a red-top by the following summer. MPB also have a mutualistic relationship withblue stain fungus [SW06], and can introduce the fungus to attacked trees,which can cause dessication and interrupt transpiration in the host tree.This combination can severely impede the nutrient and water movement(more than just the MPB alone) within the tree and induce mortality in thehost tree.The lodgepole pine has a defense mechanism against beetle attack: Itproduces resin, that fills the holes in the bark produced by the nesting MPBsand crystallizes around the beetles, effectively killing them. Furthermore,the tree can produce oleoresin, which is capable of killing the MPB eggs.The MPB must therefore mass attack in sufficient numbers to overcome hostdefenses and kill the tree.The life cycle of the MPB is generally univoltine (one generation per1Chapter 1. Introductionyear) [SW06] and is very dependent upon ambient temperature. Eggs aregenerally laid from late July to mid-August, and these eggs hatch and de-velop through four larval instars over the winter. The fourth larval instardevelops into a pupa by June and then matures into an adult by late Juneto mid-July. These adults emerge from the tree in mid-July to mid-August,fly around until they find a suitable host, and then attack it and lay theireggs. The flight of these beetles is very interesting since the movement is notgoverned by simple diffusion but also incorporates chemotactic movementdue to semiochemicals produced by the MPBs and the host trees (See Whiteand Powell and references therein [WP97]).Due to the importance of the MPB to the forestry industry, the interac-tion between MPBs and lodgepole pine trees has been extensively studiedand modeled by several authors. The beetle phenology and the impactof temperature on the beetle life cycle has been modeled by Gilbert et al[GPLB04] and Powell et al [PJLB00, PL05, PB09]. Shore and Safranyik[SS92] developed a susceptibility model for stands of trees. Susceptibilitymodels are important in forest management since they can assist in theidentification of high-risk stands of lodgepole pine. A computer model forpopulation dynamics was developed by Raffa et al [RB86]. Discrete modelsfor the spatial movement of MPB were developed by Burnell [Bur77] andGeiszler et al [GGG80]. Burnell [Bur77] incorporated both beetle dispersaland aggregation, while Geiszler et al [GGG80] incorporated the dynamicsof beetle pheromone and modeled variable wind speeds. A discrete modelfor MPB growth and lodgepole pine attack in a 1 ha stand was developedby Safranyik et al [SBTR99]. The purpose of the model was to evaluate theeffectiveness of different management strategies. The authors investigatedstand removal, stand thinning, single-tree application of MSMA (a chemicaldesigned to induce mortality of the MPB), and the use of attractive and re-pulsive pheromone. Since this model was not spatially explicit it neglectedthe spatial variation in the landscape and its effects on the managementactivities. To make the MPB model more usable for forestry management,the model was extended to the stand scale by Riel et al [RFSS04], and thenextended further to the landscape level [RFSS04]. At the stand and land-2Chapter 1. Introductionscape levels, tree removal and green-attack detection were the managementactivities investigated. Riel et al [RFSS04] did not investigate the use ofpheromone baiting at these larger spatial scales.These studies were followed by work using continuous models, such as theordinary differential equation models of Berryman et al [BSW84] and Nelsonet al [NPL+08]. The Berryman et al [BSW84] model uses two coupled ODEsto describe the interaction between the lodgepole pine forest and the MPB.The Nelson et al [NPL+08] model focuses on the interaction of a single treewith the MPB. It incorporates the damage to trees done by MPBs, whichallows the trees to be weakened for attacks by MPBs in subsequent years.The Berryman et al paper [BSW84] focuses on the spatial movement of thebeetle while the work of Nelson et al [NPL+08] focuses on the beetle biology.There have been other models that incorporate both the beetle biologyand the spatial movement of the beetle. Polymenopoulos and Long [PL90]introduced a MPB model with diffusive movement and a simple discretelinear growth term. They made the model more complex through incor-poration of directed movement of beetles toward trees with larger phloem.Later, an integro-difference model developed by Heavilin et al [HP08], whichis continuous in space and discrete in time, was used to model the spreadof the beetle infestation by considering simply the transmission of ?infec-tion? from one tree to another, without regard for the details of the beetledynamics. This model did not explicitly consider the spatial movement ofbeetles due to chemotaxis, which can be done by using a more complex spa-tially explicit partial differential equation model, as developed by Powell etal [PMW98, WP97]. It did, however, incorporate the biology and spatialmovement of the MPB, and the interaction between lodgepole pine trees andbeetles at the stand scale of the forest. An individual-based approximationof this model was developed by Hughes et al [HFSL06] to investigate theeffect of different landscapes on the spatial movement of the MPB.Previous modelling efforts have generally investigated beetle dynamicsover a single season. Few models have looked at beetle populations acrossmultiple years. Of those that have looked at multiple-year dynamics, thereare no models of which we are aware that explicitly model both population3Chapter 1. Introductiondynamics and dispersal, including the interactions between beetles, beetlepheromones, tree kairomones, and the forest landscape. These interactionsare extremely important as they are major determinants of the MPB spreaddynamics over the landscape.In this thesis, we build a multi-year spatially explicit mathematicalmodel that incorporates the interactions between beetles, beetle pheromones,tree kairomones, and the forest landscape. Our model investigates the beetledynamics on the stand scale of 1 ha - 1 km2. We are not aware of this typeof spatially explicit modelling work at this scale investigating the effects ofmanagement.We use our model to address three different questions about MPB spreadand subsequent tree damage. In Chapter 2, we use our mathematical modelto investigate the spacing between clusters of MPB attack on the forest land-scape. To verify the predictions of our mathematical model we utilize aerialdetection survey data from the Sawtooth National Recreation Area (SNRA),and then discuss the management implications of our results. Chapter 3 fo-cuses on the effects of habitat fragmentation on MPB movement. We useour theoretical model to explain results obtained in an experimental MPBstudy. Specifically, we consider a landscape with pheromone traps on eitherside of a boundary between mature and harvested forest. The fragmentedhabitat contains different combinations of dense, thinned, and clearcut forestregions. This work is extremely important as MPB movement in fragmentedhabitats is not well understood [BAL93].These first two chapters provide a basic framework for the more complexstudy and predictions presented in Chapter 4. In this last chapter, weare interested in understanding the impact of various direct and indirectmanagement strategies on limiting the spread of MPB. The managementstrategies we investigate are pheromone baiting, tree-removal, clearcutting,and prescribed burning. We find that the optimal management strategychanges with MPB density and depends on the properties of MPB movementand reproduction.4Chapter 2Pattern Formation in aModel for Mountain PineBeetle Dispersal: LinkingModel Predictions to Data2.1 IntroductionPattern formation is ubiquitous in biology [Mur03]. Nature providesa diverse array of systems with spatial patterns. Examples include staticspatial patterns such as those found on butterfly wings and mammaliancoats, and spatiotemporal patterns such as those exhibited by predator-preypopulations [Mur03]. What is even more interesting is that reaction-diffusionsystems can exhibit all of these various patterns given the correct range ofparameter values. Pattern arises in these systems through diffusion-driveninstability, which was first observed by Turing [Tur52]. An essential elementfor spatial patterning is local activation with long-range inhibition [Mur03].Not all pattern formation is due to diffusion-driven instability. Patternformation can also occur in reaction-diffusion equations with chemotaxis.Examples include models for snake pigmentation patterns and spatial pat-terns formed by colonies of growing bacteria [Mur03]. In particular, Budreneand Berg [BB91, BB95] found very diverse and interesting patterns formedby the bacteria Escherichia Coli and S. typhimurium. Tyson et al [TLM99]were able to reproduce these interesting patterns using a reaction diffusionmodel that incorporated chemotaxis.52.1. IntroductionA very interesting insect to study in regard to pattern formation is theMountain Pine Beetle (Dendroctonus ponderosae Hopkins, MPB). This barkbeetle has had a major economic impact on the Western Canada and UnitedStates forestry industries [SW06]. Many characteristics of the spread andspatial synchrony of MPB have been well researched at both large [HLZ+08,GH08, BS03, PLBW02] and small [RNJ+09, MP91] scales. These studieshave illuminated the factors driving the spatial patterns in beetle spread,such as weather, elevation, and proximity to nearby MPB attacked areas.Additionally, models for MPB movement have been developed to describethe spread and aggregative behaviour of MPB [LWBP98, BPL96, PD11,PL90, SW06, Bur77, GGG80, SBTR99, RFSS04, HP08, PMW98, WP97,HFSL06]. In this paper, we will focus on the spot formation that occursat intermediate spatial scales. Previous cluster analysis focused on eitherlarge scales (kms) or very small scales (100m regions). In this study we areinterested in the pattern formation of clusters at intermediate scales, withdistances between clusters falling within the 0-1000m range. We also focusour search to investigate spot formation in a single year rather than thechange in spot formation across multiple years.Successful aggregation of MPB (in response to a suite of pheromones) iscrucial for reproduction and survival of the species. At incipient epidemicpopulation levels, the pattern of attack of MPB on a landscape is smallisolated spots [SW06]. In contrast there can be high mortality of host treesover thousands of contiguous acres at epidemic population levels. We are in-terested in understanding the pattern formed by MPB during the transitionfrom incipient epidemic to frank outbreaks. To do this work we investigatepattern formation via a spatially explicit model for MPB dispersal [WP97].We then compare the model predictions to 19 years of data from MPB at-tack in the Sawtooth National Recreation Area (SNRA), in Idaho, USA. Wefind that the distance between MPB attack clusters predicted by our modeland observed in the SNRA are the same. This indicates that the biolog-ical behaviours in our model are sufficient to explain the observed attackpattern.We first introduce our mathematical model in Section 2.2. In Section 2.3,62.2. Mathematical Modelwe determine the wavelength in between MPB attack clusters as predicted byour mathematical model using pattern formation analysis. We then calculatethe wavelength between clusters of MPB attack in the data from the SNRAin Idaho, USA and compare it to the model predictions in Section 2.4. Themanagement implications of our results are in Section 2.5. Discussion andfuture work can be found in Section 2.6.2.2 Mathematical ModelWe are interested in the period of emergence, dispersal, and attackof the MPB. These events all occur over the space of one summer: adultsemerge from their host trees, then aggregate on new hosts where they mountan attack and, if successful, lay new eggs. The period between egg-layingand adult emergence occurs during the winter, and is not relevant to themodelling exercise here. That is, we are interested in understanding theattack pattern that results from the dispersal and aggregation stages of theMPB life cycle in a single summer. We thus require a model for MPBmovement through forest habitat that incorporates the interaction betweenthe beetle and its pheromones in a continuous framework over space andtime. The choice of model structure is based on theoretical work by Powellet al [PMW98]. The model equations are?P?t=diffusion? ?? ??p?2P ?chemotaxis? ?? ??[(?ab0 ?Ab0 +A/b1?A)P]?death?????pP ?nesting? ?? ??(x)PP 2P 2 + (kp)2+emergence? ?? ??(x, t) ,(2.1a)?Q?t=death? ?? ???qQ+nesting? ?? ??(x)PP 2P 2 + (kp)2, (2.1b)?A?t=diffusion? ?? ??a?2A+synthesis????a1Q ?degradation?????aA , (2.1c)72.2. Mathematical Modelwhere our three variables are P - the density of flying MPB, Q - the densityof nesting MPB, and A - the concentration of beetle pheromone. The modelwe have chosen allows us to investigate the dynamics of MPB attack in asingle summer season where the emergence rate, ?(x, t), is determined bythe position and severity of MPB attacks in the previous year. That is,given the pattern of MPB attacks in the previous year, we can predict theemergence rate of flying MPB and the resulting pattern of attacks in thecurrent year.The movement of MPB is described by two processes: diffusion andchemotaxis. These are the first two terms in (2.1a). The diffusion compo-nent describes the random movement of flying MPB, while the chemotaxisdescribes the attraction and repulsion of MPB according to the concentra-tion of MPB pheromones. MPB have a biological mechanism whereby thepheromone suite is attractive at low densities and repulsive after the con-centration becomes too high [WP97]. As a result, the density of beetlesattacking a given tree stays below overcrowding levels (though in epidemicsituations, when tree resources are limiting, beetles will attack in higher,suboptimal densities) [RB83].Our model is structurally similar to the model in Powell et al [PMW98],but differs in that the detailed interaction between MPB and lodgepole pinetrees in the original model (holes and resin dynamics) has been replacedby a sigmoidal curve, P 2/(P 2 + kp2), multiplied by a random landing rate,?(x). We discuss both of these terms in some detail here, as they framea novel description of the MPB response to the pheromone and susceptibletree landscape. The random landing rate, ?(x), is spatially dependent basedon the density of susceptible trees on the landscape. A lower susceptible treedensity results in a lower landing rate ?. In this manner, we can includethe effects of spatial heterogeneity on the MPB aggregation behaviour. Weexpect this heterogeneity to affect the spatial distribution of MPB attacks.The type 3 functional response term assumes that the MPB must attack insufficient densities to successfully nest in a lodgepole pine tree [SW06]. Thisfunction is defined such that a low density of attacking beetles has a verylow success rate until the MPB density reaches a population threshold, at82.2. Mathematical ModelTable 2.1: Table of parameter values for the dimensional model (2.1). Notethat the unit ?fh? refers to flight hour of MPB.Parameter Description Units Value?p diffusion of flying MPB hafh 1?a diffusion of beetle pheromones hafh 0.648?a beetle pheromone attractiveness ha2?g?fh 5.7b0 concentration of pheromones at which dis-sipation occurs?gha 5.4b1 concentration at which pheromone is sat-uratedn/a 1? random landing rate of flying MPB hatrees?fh 0.16a1 rate of pheromone increase due to nestingMPB?gfh?mpb 0.02kp flying beetle density required for 50 per-cent success of mass attackmpbha 250000?p death rate of flying MPB fh?1 0.014?q death rate of nesting MPB fh?1 0.001?a degradation of beetle pheromone fh?1 180which point the success of the MPB increases dramatically. This populationdensity threshold is kp, the MPB density required for 50% attack successrate. This parameter was estimated based on the empirical data providedin Raffa et al. [RB83]. All other parameters were based on estimates byBiesinger et al. [BPBL00].We assume that nesting MPB (at density Q) have a small linear deathrate, ?q, since they have successfully penetrated the tree defenses. OnceMPB nest they do not move spatially and therefore (2.1b) contains only re-action terms. The suite of MPB pheromones, A, undergoes three processes.The pheromones undergo movement through diffusion. Furthermore, thepheromones are produced by nesting MPB and degrade naturally at a lin-ear rate, ?a.Previous theoretical and empirical work [BPBL00, RB83] informed theselection of parameter values chosen for this study. These values are dis-played in Table 2.1.92.2. Mathematical Model2.2.1 Non-dimensionalizationTo simplify model analysis we non-dimensionalize the model. Thedimensionless variables are:Q =b0?aa1Q, P =b0?aa1P , A = b0A,t =1?at, (x, y) =??a?a(x, y).(2.2)The choice of non-dimensional scalings can be interpreted biologically.The density of nesting MPB, Q, and the concentration of MPB pheromones,A, were scaled by the density of nesting MPB, 48600 MPB/ha, and theconcentration of pheromone (b0) required for the pheromone to switch fromattractive to repulsive, respectively. The density of flying MPB, P , wasscaled by the same factor as Q to remain consistent. Time was scaled bythe average degradation time of the chemical pheromone. Space was scaledby the average distance that the pheromone will spread before degradation.Inspection of the parameters in Table 2.1 reveals order of magnitudedifferences in the parameter values. We therefore defined a scaling parameterthat identifies parameters as relatively small or large when compared toother parameter values. We chose the scaling parameter describing therelative persistence of a pheromone plume, 1?a , and the life expectancy ofthe dispersing MPB, 1?p . =?p?a.Since ?p (death rate of flying MPB) is very small compared to ?a (degra-dation rate of MPB pheromone), this ratio is a very small quantity. Thisorder parameter allows us to identify parameters that work over fast andslow scales, respectively. With this scaling parameter we define the following102.2. Mathematical ModelTable 2.2: Table of parameter values for the non-dimensional model (2.4).Parameter Value?p 1.54?a 47.5 0.0000778kp 5.14?q 0.0714? 11.4new dimensionless parameters:?p =?p?a, ?a =?ab0?a, kp =kpa1b0?a,? =?a1b0?2a, ?q =?q?a, ? =??a(2.3)Values of the non-dimensional parameters (2.3) are shown in Table 2.2.Substituting the non-dimensional parameters (2.3) and variables (2.2) into(2.1), we arrive at?P?t=?p?2P ? ?a?(1?A1 +A/b1P?A)? P ? ? P3P2+ kp2 + ?, (2.4a)?A?t=?2A+Q?A, (2.4b)?Q?t=? ?qQ+ ?P3P2+ kp2 , (2.4c)For the remainder of the chapter we will drop the bars above the nondi-mensional quantities and assume that we are using the nondimensional vari-ables and parameter values to simplify notation. The model becomes,?P?t=?p?2P ? ?a?(1?A1 +A/b1P?A)? P ? ? P3P 2 + kp2+ ?, (2.5a)?A?t=?2A+Q?A, (2.5b)112.3. Model Pattern Formation?Q?t=? ?qQ+ ?P 3P 2 + kp2. (2.5c)In certain PDE models that include chemotaxis the system of equa-tions can exhibit unrealistic unbounded solutions due to finite-time blowup [HP09]. The proof that our system of equations, (2.5), exhibit globalexistence in two dimensions is outside of the scope of this work. However,we conjecture that the chemotaxis and repulsion of pheromone at high con-centrations is a similar mechanism to the volume-filling equations in Hillenand Painter [HP09], which are known to not exhibit finite-time blow up.2.3 Model Pattern FormationThe type of pattern formation we investigate is diffusion and chemotaxis-driven instability of a spatially uniform steady state [Mur03]. Biologically,this amounts to assuming that early in the season emerging MPB are uni-formly dispersed on a landscape devoid of chemical information and thathot spots of infestation will develop at spatial scales on which the naturalprocesses of dispersal, attack, and pheromone production/dissipation res-onate. Often the distribution of previously attacked trees is clustered andnon-uniform. Early in the season, however, when beetles are just beginningto emerge, there is no pre-existing chemical information. Consequently, therandom dispersal of emerging MPB, coupled with a brief maturation periodduring which beetles are not in search of nesting sites [HFSL06], togethergenerate a largely uniform distribution of beetles. The distribution is clearlynot completely uniform, and the spacing of previously attacked trees shouldhave some effect on the spatial pattern of attack in the following year. Thispre-pattern, however, would most likely amplify and accelarate the patternof MPB attack clusters through a forcing of the inherent spatial resonance.In other words, it is less surprising to see a spatial pattern emerge when itis seeded with a pre-pattern than when it is seeded with no pattern at all.We therefore take the more parsimonious assumption, and set attack andemergence rates to be spatially uniform ?(x) = ? and ?(x, t) = ?(t). We122.3. Model Pattern Formationadditionally set our emergence rate to be constant over time, ?(t) = ?. Toinvestigate potential pattern formation, we find the spatially uniform steadystate in Section 2.3.1, then linearize about this steady state and add spa-tial perturbations to find the dispersion relation in Section 2.3.2 [Mur03].This dispersion relation relates the temporal growth of perturbations to thewavenumber of the pattern. This dispersion relation is studied analyticallyand numerically to determine the dominant wavenumber of the pattern, inSections 2.3.3 and 2.3.4, respectively. This dominant wavenumber predictsthe expected spacing between MPB attacks in a given year.2.3.1 Spatially Uniform Steady StateTo find spatially uniform steady states for (2.5), all spatial and tem-poral derivatives are set to zero. The spatially uniform steady state of themodel, (2.5), is the solution to the system0 =? P ? ? P3P 2 + kp2+ ?, (2.6a)0 =Q?A, (2.6b)0 =? ?qQ+ ?P 3P 2 + kp2, (2.6c)Solving (2.6b) and (2.6c), we find that there is a spatially uniform steadystate given by(P ?, A?, Q?) =(P ?,1?q?(P ?)3(P ?)2 + k2p,1?q?(P ?)3(P ?)2 + k2p). (2.7)P ? is unknown. Rearranging (2.6a), we obtain a cubic equation,P 3 ? ?(1 + ?)P 2 +k2p1 + ?P ??k2p(1 + ?)= 0. (2.8)The real positive root of this cubic equation (2.8) will give us the steadystate density of flying MPB, P ?. The first derivative of (2.6a) is always132.3. Model Pattern Formationnegative and therefore there is only one possible real root. This steadystate is positive by inspection of (2.6a). Emergence rate, ?, is determinedexogenously by the density of MPB attack in the previous year and thetemperature (phenology) [BPL96], which will then determine the uniquedensity of dispersing MPB.This system exhibits three distinct behaviours as the emergence rate,?, is varied. For very low ?, the flying MPB density, P ?, is very low. Atthese population levels we have essentially the trivial steady state and theMPB do not successfully attack any trees in the susceptible landscape. If ?is O(1), the largest root scales like P ? ? 1 . Using this scale one finds thesteady state is approximately(P ?, A?, Q?) =(?(1 + ?),??qP ?,??qP ?). (2.9)At these emergence rates, P ? is at epidemic densities, and therefore the MPBis successful at inducing mortality in any healthy tree. In contrast, when? is small, O(), the roots of P ? are O(1). When ? is at these values, thepopulation of MPB is not large enough to kill susceptible trees easily, andmust successfully aggregate to overcome tree defenses. We will show throughour analysis that for a specific intermediate range of ? values, our uniquereal steady state is an unstable critical point in the presence of diffusionand chemotaxis. Thus, when spatial factors are included, perturbations ofthe uniform steady state lead to the formation of a spatial pattern. Thisscale of aggregation will determine the spacing between MPB attacks on asusceptible landscape.The bifurcation plot of these behaviours around the single steady stateis discussed later in Section 2.5 and is shown in Figure 2.4.2.3.2 Linear AnalysisLinearizing about the spatially uniform steady state, we define f =?a(1?Q?1+Q?/b1)P ? and g = ?(P ?)4+3(P ?)2k2p((P ?)2+k2p)2. Note that since Q? = A?, wereplace A? by Q? in our equations. We add small spatial perturbations by142.3. Model Pattern Formationsubstituting P = P ? + ?P1, Q = Q? + ?Q1, and (? << 1) A = A? + ?A1,into (2.5) to obtain the perturbation equations,?P1?t=?p?2P1 ? f?2A1 ? P1 ? gP1, (2.10a)?Q1?t=? ?qQ1 + gP1, (2.10b)?A1?t=?2A1 +Q1 ?A1. (2.10c)Note that we do not drop terms with  in our linearization, (2.10), andthus the  contained within g has no effect on the linearization or futureanalysis.Method of Annihilators to find dispersion relationOur analysis focusses on one dimensional results, where we do notdifferentiate between the possible differential characteristic distances in be-tween clusters of attack in the x and y directions. We assume that thedomain is large with respect to the insect dispersal distance, and thereforethere is no limitation on possible wavenumbers over the domain. Our anal-ysis in this section would be largely unchanged in two dimensions, otherthan the wavenumber found would consist of both x and y components.For further details see Tyson [Tys96]. Our first-order analysis does not re-quire statement of boundary conditions, although any non-linear analysis(higher-order) would require boundary conditions.Beginning with (2.10), we can rewrite these equations in terms of lineardifferential operators [NSS05]:L1[P1] = ?fL2[A1], (2.11a)L3[Q1] = g[P1], (2.11b)L4[A1] = [Q1], (2.11c)where L1 = (?t ? ?p?xx + ( + g)), L2 = ?xx, L3 = (?t + ?q), and L4 =152.3. Model Pattern Formation(?t ? ?xx + 1).From (2.11) we deduce that L1L3L4[A1] = ?gfL2[A1]. If the linearoperators are then expanded, we have the equation:(?t ? ?p?xx + (+ g))(?t + ?q)(?t ? ?xx + 1)[A1] = ?gf?xx[A1]. (2.12)We assume that the perturbations have an exponential solution of thefollowing form, A1 = c1e?t+i?mx. This substitution in (2.12) produces thedispersion relation that links the temporal growth rate, ?, of patterns totheir spatial wavenumber, ?m,(? + ?p?2m + (+ g))(? + ?q)(? + ?2m + 1) = gf?2m. (2.13)If the polynomial (2.13) is expanded in terms of powers of ?, we have,?3 + ?2(+ ?p?2m + ?2m + g + 1 + ?q) + ?(?2m + ?q + ?p?4m+?p?2m?q + ?p?2m + g?2m + g + g?q + ?2m?q + 2?q + )? gf?2m + g?q + ?p?4m?q + 2?2m?q + g?2m?q + 2?q + ?p?2m?q = 0.(2.14)2.3.3 Analysis of the Dispersion RelationBefore turning to numerical analysis, we first utilize some analyticaltechniques to find the boundary of the region of maximum pattern forma-tion, and to determine the dominant wavenumber. The dominant wavenum-ber is the spatial wavenumber that maximizes the temporal growth rate. Werewrite the dispersion relation (2.14) asp1 = ?3 + a2?2 + a1? + a0, (2.15)wherea2 = + ?p?2m + ?2m + g + 1 + ?q,a1 = ?2m + ?q + ?p?4m + ?p?2m?q + ?p?2m + g?2m + g + g?q + ?2m?q162.3. Model Pattern Formation+ 2?q + ,a0 = ?gf?2m + g?q + ?p?4m?q + 2?2m?q + g?2m?q + 2?q + ?p?2m?q.Thus, p1 is a function of ?.We are interested in situations where pattern formation occurs, thatis, where ?1 = max??C(<(?)|p1(?) = 0)) > 0. Using Descartes? rule ofsigns [SL54], we are able to determine regions in which positive real rootsshould occur. Descartes? rule of signs counts the number of sign changes ofthe coefficients of a polynomial to determine the maximum number of realpositive roots. Furthermore, if there is a maximum of n real positive roots,the number of allowable roots is n, n ? 2, n ? 4, ... because complex rootsmust occur in pairs. Therefore, in order for a real positive root to occur, wemust have a0 < 0 (one sign change). This is equivalent to the condition,gf?2m > g?q + ?p?4m?q + 2?2m?q + g?2m?q + 2?q + ?p?2m?q. (2.16)Technically, Descartes? rule of signs limits us to a single positive realroot, and either zero or two negative real roots. In the case where thereare two negative real roots, we do not need to know anything about thesenegative real roots, as pattern formation occurs if a single root has a positivereal part. When there are not two negative real roots, but a0 < 0, we willargue that the two complex roots have a negative real part.Assume our cubic has a positive real root r1 (r1 > 0), and a pair ofcomplex roots r2 ? r3i. The expanded form of the polynomial is:?3 + ?2(?r1 ? 2r2) + ?(2r1r2 + r22 + r23)? r1(r22 + r23). (2.17)Since a2 > 0 in (2.15), we must have ?r1 ? 2r2 > 0 in (2.17). Thus, sincer1 > 0 by assumption, we have r2 < 0. Therefore, in the case where we havea single positive real root and two complex roots, our complex roots musthave negative real parts.In the case where a0 > 0, Descartes? rule of signs indicates that thereare no positive real roots and that the maximum number of negative real172.3. Model Pattern Formationroots is 3. Therefore, there is either 1 or 3 negative real roots for (2.15).Obviously, in the case of three negative real roots, no pattern formationcan occur. There is the possibility of a single negative root, and a pairof complex roots with positive real parts. This case does not occur in theparameter space we explored in the numerical determination of the roots of(2.15) (in Section 2.3.4).This means that our pattern formation analysis is restricted to thecase where (2.15) has a single positive real root. We focus on this regionwhen trying to determine the maximum region of pattern formation. Pat-terns will first form at wavelengths, ?m, and parameters chosen such that?1 first becomes positive. Therefore, we examine the behaviour of maxi-mum ? with respect to ?m in (2.15). We know that the maximum pat-tern formation region with respect to the wavenumber, ?m, will occur whenp2 =?p1?(?2m)= 0. That is, if we take a derivative of p1 with respect to thesquare of the wavenumber and set it to zero, we can determine the dom-inant wavenumber. The wavenumber at which pattern formation is maxi-mum is called the dominant wavenumber. At this dominant wavenumber,???(?2m)= 0, thus we can reduce our polynomial (2.15) to order 2 by takingthe derivative. In summary, the dominant wavenumber occurs when ?1 > 0,?2 = max??C(<(?)|p2(?) = 0)) > 0, and ?1 = ?2. The last condition mustbe satisfied since both p1 = 0 and p2 = 0 at the dominant wavenumber.Thus we have,p2 =?a2??2m?2 +?a1??2m? +?a0??2m,= d2?2 + d1? + d0.where,d2 = (?p + 1),d1 = (+ 2?p?2m + ?p?q + ?p + g + ?q),d0 = (?gf + 2?p?2m?q + 2?q + g?q + ?p?q),Since g > 0, we have that d2 > 0 and d1 > 0. This means that there182.3. Model Pattern Formationis exactly one positive real root of p2 if d0 < 0. Therefore, ?1 > 0 if andonly if a0 < 0, and ?2 > 0 if and only if d0 < 0. Using these conditionswe can sketch the region of pattern formation (Figure 2.1). Additionally,we numerically calculate (using the method in Section 2.3.4) ?1 and ?2 andfind the squared difference. If the squared difference is zero, this signifies apoint of intersection between the two curves and a maximum value of ? withrespect to ?2m. In short, the dominant wavenumber occurs at the intersectionof ?1 and ?2, which is shown on the contour plot.2.3.4 Numerical Analysis of the Dispersion RelationUsing a root-finding algorithm in Matlab, we calculated ?1 whilevarying ?, the emergence rate (of flying MPB) at steady state, and ?m, thespatial wavenumber. The contour plot produced is shown in Figure 2.1.Additionally, we show the 2-dimensional plot of growth rate with respectto wavenumber. This curve shows the maximum growth rate with respectto emergence rate at each wavenumber. In our calculations we scaled thewavenumber, ?m, so that it would be the reciprocal of wavelength. Thewavelength, wm can be calculated as,wm =2pi?m. (2.18)The maximal eigenvalue of 8.04e-5 in the contour plot is dimensionless,and when redimensionalized becomes 0.0145 fh?1, which is an appropriatetime scale for pattern formation within a single summer (corresponding to80 fh, or 2-4 weeks of the flight season). The maximal eigenvalue occuredat the dominant wavenumber of 2.74 km?1, with an emergence rate of 300MPB/(ha?fh). Assuming an output of approximately 10,000 MPB/tree perflight season [PB09], this corresponds to approximately 2-3 source trees/ha.The resulting steady state density of nesting MPB is Q? ? 20100 MPB/ha.Using the conversion of 800 nesting MPB/tree [PB09], we find the resultingnumber of killed trees due to this attack is approximately 25. Therefore,this pattern of aggregation is important for the transition between incipientepidemic and epidemic densities of MPB. In an incipient epidemic [SW06],192.4. Data Analysisthe MPB population can form small clusters of attack rather than attacksoccuring on large tracts of continuous forest (which will occur at epidemicdensities). These source densities would describe the transition from incipi-ent epidemic to epidemic densities.There is great agreement between the analytical and numerical determi-nations of the dominant wavenumber and the region of pattern formation.The analytical method correctly identifies the region where there is a singlepositive root, and by Figure 2.1, this is exactly the same as the region where?1 > 0 (pattern formation occurs). Additionally, the wavenumber at which?1 = ?2 in Figure 2.1 verifies the numerical calculation of the dominantwavenumber at 2.74 km?1, which corresponds to a wavelength of 364 m.Thus, our analytical and numerical work yields the prediction that attackclusters during the transition between incipient epidemic and epidemic pop-ulation levels will be approximately 364 m apart. This prediction is basedon our model assumptions which include landscape homogeneity, chemotac-tic response of MPB due to pheromone (both attractive/repulsive), and atype three functional response describing the transition between flying andnesting MPB.2.4 Data AnalysisThe second component of this project was to compare spatial dataof MPB attacks to the model predictions. We analyzed MPB attack datafrom the Sawtooth National Recreation Area (SNRA), located in the RockyMountains of central Idaho, to identify any characteristic distances betweenpatches of beetle infestation. Data was provided by USDA Forest Serviceaerial detection survey (ADS) in and around the SNRA. Full details areprovided in Crabb et al. [CPB12]. The data set extends over a period of 19years, 1991-2009. The data are remarkably detailed, taken at a grid-scale of30 m over a region of 275,776 ha. All 19 years include regions where MPBare at incipient epidemic densities; many of these years also have regionswith epidemic densities of MPB. Thus, this data tracks the progression ofan MPB epidemic as captured by dead (red top) trees. Trees develop a red202.4. Data Analysis?1=0max(?1)0 50 100 1500501001502002503003504004505000 50 100 150050100150200250300350400450500?1=?2 ?2=0 ?1=00 50 100 150 2000123456789 x 10?50 1 2 3 4 50123456789 x 10?5bac dWavenumber, ?m (1/km)Wavenumber, ?m (1/km)MaximumRealRoot,?EmergenceRate,?(MPBfh?ha)Figure 2.1: Numerical contour plot of the temporal growth rate, ?1, againstwavenumber and emergence rate of MPB (a). The maximum value of ?1 islabeled with a diamond at (2.74, 299). The surface is unimodal with a singlemaximum. The analytical contour plot of temporal growth rate againstwavenumber and emergence rate of MPB is shown in (b). The contourplot in (b), uses the analysis of the dispersion relation to find the region ofpattern formation (?1 = 0), and the dominant wavenumber (?1 = ?2). Thecontour plot in (a), shows the region of pattern formation and the dominantwavenumber as determined numerically. The plot (c), is a horizontal sliceof the surface in (a) at a fixed emergence rate, ? = 299 MPB/(fh?ha). Azoomed in portion of (c) is shown in (d) where the wavenumber varies from0 to 5. 212.4. Data Analysistop the following summer after begin attacked as a result of beetle-inducedmortality. The initial attacks at incipient epidemic levels resulted in smallclusters of dead trees. During the 19 years, many of these populations hadrisen to epidemic levels and killed a more significant portion of the availablepine trees over the landscape. The later years of this dataset capture theperiod after the epidemic where the MPB population density has decreasedto lower levels. For the purposes of our analysis, the data was definedwith ones given to grid cells (locations) of MPB attack, and zeroes given tolocations with no MPB attack.The data sets chosen for analysis in each year were from areas of incipientepidemic MPB population densities, consistent with the assumptions usedin the linearization of the model. That is, we selected regions with smallspots (? 300 m in diameter) of MPB attack and at least 5 spots per region.Regions with large spots were characteristic of MPB at epidemic densities.For the size of these regions, we picked the largest region such that thesetwo conditions were satisfied. The error in the computation of wavenumberdecreased as the size of the chosen region increased because larger regionsincreased resolution in Fourier space. To calculate the distance betweenspots of MPB attack we used discrete fast Fourier transforms. DiscreteFourier transforms assumes that the data can be decomposed into a finitenumber of sine and cosine functions on a grid. The process returns theamplitude of these sine and cosine functions, the wavenumbers with thelargest amplitude best describe the scale of aggregation of MPB attack in thedata. The particular DFT used was Matlab?s fft2, with which we calculatedthe radial wavenumber, ?r,?r =??2x + ?2y , (2.19)where ?x and ?y are spatial wavenumbers of the data in the x and y direc-tions. Since DFT returns the amplitude of both the sine and cosine com-ponents of the data we need to compute the power, which is the squaredcomplex modulus of the amplitude. This factors both the sine and cosineamplitudes at each wavenumber into a single value, the power. An example222.4. Data Analysisregion and power spectral density is shown in the Appendix (Figure A.1).We are interested mainly in the wavenumber at which the power is max-imum. This is called the dominant wavenumber and is the most influentialwavenumber represented in the data. To have error bounds we computedthe upper (?u) and lower (?l) bounds for the dominant wavenumber, ?d,which were chosen such that:?l = {min(?)|m(?) ? 0.80m(?d)},?u = {max(?)|m(?) ? 0.80m(?d)},where m(?) is the power at the wavenumber, ?.Since multiple regions were chosen in each year, the average dominantwavenumber in a given year was calculated as a weighted average based onthe power:?d = ?i(mi?jmj?i). (2.20)where ?i is the dominant wavenumber for region i and mi is the power at themaximum. Average upper and lower bounds for each year were calculatedsimilarly.The average dominant wavenumber in each year is displayed in Figure 2.2over 1991-2009. The average dominant wavenumber varies between 1.5 and5.5 km?1. The dominant wavenumber appears to be higher in 1991-2000than in the years 2001-2009. The mean dominant wavenumber is calcu-lated to be 2.83 km?1, which is very close to the model predicted dominantwavenumber of 2.74 km?1. For more details and analysis from each year seeAppendix A.1.The frequency of each dominant wavenumber independent of year isgiven in Figure 2.3. The dominant wavenumber from the model, ?d,model, isvalidated by the data, as it is close to the center of the distribution of ?d.In fact, the majority of ?d appearing is enclosed between the lower (?l,model)and upper (?u,model) bounds on the model dominant wavenumber.A weighted histogram was also produced, where each count is scaled bythe relative power at that ?d. This measure is important as it highlights the232.4. Data Analysis1990 1995 2000 2005 20101234567 ?d?l?uMean(?d)?d,modelWavenumber (1/km)YearFigure 2.2: Average dominant wavenumber over years 1991-2009. ?d refersto the average dominant wavenumber, ?l and ?u refer to the lower and upperbounds on ?d, Mean(?d) refers to the average of the dominant wavenumberover all years, and ?d,model refers to the model predicted dominant wavenum-ber. This graph shows the trend in ?d over the years 1991-2009. Points where?u and ?l are near ?d represent years in the data where the power spectraldensity shows a very large sharp peak at the given ?d.242.5. Management Implicationswavenumbers which more strongly represent patterns in MPB attacks inthe data. Similar to the first histogram, ?d,model provided a good estimateof the center of the distribution of ?d (data), and a large proportion thedistribution of ?d was effectively captured within the range between ?l,modeland ?u,model.For the histogram, the upper and lower bounds for the model dominantwavenumber were chosen such that:?l,model = {min(?)|?1(?) ? 0.975?1(?d,model)},?u,model = {max(?)|?1(?) ? 0.975?1(?d,model)},where ?1(?) is the growth rate at the wavenumber, ?.We found an interesting trend when analyzing the data; As time pro-gressed, the spots of MPB attack became larger and farther apart (resultsnot shown). This trend in the spot pattern is intriguing and it would bevaluable to investigate if this trend is characteristic of the progression of anMPB epidemic.2.5 Management ImplicationsHere we consider management implications ensuing from our analy-sis. An important management goal is the disruption of the MPB aggrega-tion, or pattern-formation process. Clearly this can be done by reducing theMPB population density, but what level of reduction is needed? To answerthis question, we look for MPB population densities at which patterns donot occur. The regions of parameter space with and without patterns areshown in Figure 2.4.Patterns do not form at any wavenumber for very large and very smallMPB emergence rates. The very large population densities correspond toepidemic populations, which occur above ?max=423 MPB/(ha?fh). This cor-responds to approximately 4 source trees per hectare. At these densities thespot patterns disappear and are replaced with an area-wide infestation. Weare interested in the low MPB population density below which infestation is252.5. Management Implications0 1 2 3 4 505101520253035400 1 2 3 4 505001000150020002500300035004000Wavenumber (1/km)Wavenumber (1/km)Number of Dominant WavenumbersScaled Number of Dominant Wavenumbers?d,model?l,model ?u,model?d,model?l,model ?u,model abFigure 2.3: The histogram of dominant wavenumbers over all years (a) andthe histogram weighted by the power at each dominant wavenumber over allyears(b). The solid line is the model dominant wavenumber and the dottedlines are upper (5.5069 km?1) and lower (1.3783 km?1) bounds on the modeldominant wavenumber.262.5. Management Implications0 100 200 300 400 5000123456789 x 10?50 50 100 150050100150200250300350400450500Wavenumber(1/km)aMaximum Real RootEmergence Rate (MPB/(ha*fh))Emergence Rate (MPB/(ha*fh))bPattern No Pattern EpidemicEndemicIncipient EpidemicMax EmergenceMin EmergenceFigure 2.4: The contour plot (a) shows the region of pattern formation,and the range of population densities at which the MPB are displayingpatterns consistent with endemic, incipient epidemic, and epidemic densities.The separation in between these densities are ?max=423 MPB/(fh?ha) and?min=82 MPB/(fh?ha). The contour plot in (a) taken at a vertical cross-section at the dominant wavenumber ?d = 2.74, shows the temporal growthof patterns (maximum real root) against the emergence rate of MPB (shownin (b)).272.6. Discussionendemic or absent. This level is ?min= 82 MPB/(ha?fh), which correspondsto approximately 2 source trees per 3 ha. At this density, the pattern forma-tion would not occur, and therefore the MPB would not be able to reproducesuccessfully due to low population densities. Thus, if control measures couldreduce the impact to 2 source trees per 3 ha, our model predicts the MPBpopulations would decrease to endemic levels. The reductions in density tothe MPB populations could be done through various management activitiessuch as thinning, prescribed burning, or tree removal [SW06].2.6 DiscussionModel pattern formation analysis predicts a dominant wavenumberof 2.74, or spots of MPB attack that are 364 m apart. Analysis of SNRAspot data indicates spots are 353 m apart on average, with a wavenumber of2.83, only a 3% difference between the model predictions and field observa-tions. Since our model was parameterized completely independently of theSNRA data, this correspondence between the model and data gives a strongvalidation of our model.The model analysis revealed that clusters of MPB attack could be re-moved if the number of source trees per hectare can be reduced by manage-ment activities to 2 trees per 3 hectares. To accomplish this goal, it mightbe possible to use pheromone baits (repulsive and attractive) to disrupt theaggregation process. Our modelling approach could be used to determinewhether or not judicious placement of pheromone baits could completely orpartially hinder the formation of spot aggregates, and the number and place-ment of baits necessary to prevent the transition from incipient epidemic toepidemic.Past forestry practises and fire suppression have given rise to homoge-neous and even-age stands of lodgepole pine that have led to large outbreaksof MPB [SL00]. Our model could be used to determine what types of distri-butions of susceptible trees (and mixed tree forests) over a landscape wouldprevent pattern formation of MPB attack, and therefore stop a transitionto epidemic densities. For this analysis, we would need to use a multi-year282.6. Discussionmodel, as our current analysis will indicate transition emergence rates, butwill not model the transition from endemic to epidemic densities. Our anal-ysis in Chapter 4 tries to understand the density of source trees needed fortransition from endemic to epidemic densities.Our model assumes the habitat is homogeneous and therefore the wave-length between attacks formed by MPB predicted in our model is drivenby the intrinsic biology of the MPB. This means that at incipient epidemicdensities there can be development of aggregration pattern that is drivenby the MPB movement dynamics and not heterogeneities in the landscape.A previous study by White and Powell [WP97] found that the patterns ob-served at endemic densities were driven by the landscape, while patternsobserved at the epidemic densities are driven by the self-focusing dynamicsof the MPB. Our study adds to this work by finding that the patterns at theincipient epidemic density (between the endemic and epidemic levels studiedby White and Powell) can be explained by the MPB biology.Our results are fairly robust to changes in parameter values. Parametersensitivity analysis showed that our model wavelength prediction is mostsensitive to increases in ?p, the diffusion rate of MPB. This result is expected,intuitively, as diffusion is known to smooth patterns when pattern formationis driven by chemotaxis. All parameters were altered by 10% (while keepingall other parameters constant), and the most sensitive parameter, ?p, onlychanged the wavelength prediction by at most 3.8%. These values can beseen in Appendix A.2. This means that our pattern formation analysisis relatively insensitive to parameter changes. Parameter sensitivity wasmeasured as a ratio of standardized changes in wavenumber to standardizedchanges in parameter values [Hae05]. Future work needs to investigate theeffects of varying multiple parameters simultaneously.An important factor not included in our model is the effect of tempera-ture, which has been shown to have a significant effect on MPB emergenceand spread [HLZ+08, GPLB04, PB09]. If temperature changes, throughglobal warming [PL05], habitats that were previously unsuitable for MPBmay become suitable. Additionally, temperature changes can increase thesynchrony of emergence, which could increase the density of MPB attacks292.6. Discussionand allow for spot formation. An interesting extension of our mathematicalmodel would explore how the predicted wavelength changes as this factor isincluded. This could be done by making the emergence rate, ?, a functionof temperature. This modification would add a new layer to the complexityof the pattern formation analysis, and may require numerical simulations todetermine the expected wavelength between clusters in a given landscape.The functional form of our nesting term is required for pattern formation,as smaller submodels with the same chemotaxis and diffusion terms willresult in pattern formation. This is due to the local activation throughMPB response to pheromone and the long-range inhibition of diffusion. Ourspecific sigmoidal nesting term is required, however, to obtain the correctwavelength and close correspondence with the data.From an analytical standpoint, it would be interesting to complete thesecond-order perturbation analysis of these equations [Mur03, Tys96]. Weonly investigate spot aggregation patterns in the present work, but it maybe that other patterns are possible. An understanding of the possible aggre-gation patterns would provide managers with an additional tool for gaugingthe MPB population level and risk of an epidemic.A final interesting extension of our work would be to determine thetime it takes for the MPB population to reach epidemic levels once thecharacteristic wavelength of pattern formation (364 m) has been established.In order to do this, the model would require a between-season component todescribe the over-winter reproduction and development of MPB. This studyis currently in progress by the Tyson lab.Our modelling approach can be applied to other organisms that ex-hibit patchy spread [SKT95]. Examples include birds such as house finches[LP00], sparrows, and starlings [SKT95]. Additionally, there are insects whoexhibit patchy spread, such as are rice weevils [SKT95], emerald ash borer,leaf-miner moth, pinewood nematode, corn rootworm [CMM+10], and gypsymoth [PM10]. Some plants, such as cheat grass [LP00], also exhibit patchyspread. Note that the patchiness of the MPB spread is an inherent propertyof the chemotactic behaviour of the insects. The patchy spread exhibited bythe other organisms may be due to other factors, such as patchy resources,302.6. Discussionor could be influence by the same inherent mechanisms as the MPB andcomparisons would be both interesting and instructive. Many of the pop-ulations mentioned are invasive species and so an understanding of theirspatial invasion dynamics is vitally important, as invasive species can dev-astate populations of native flora and fauna [MSL+00]. Additionally, oncewe understand how a species may successfully invade we may be use thisknowledge to enhance survival of native species which are threatened byextinction.31Chapter 3Edge Effects on MountainPine Beetle Movement3.1 IntroductionFragmentation of habitat can have diverse effects on ecological com-munities [FCC99], through edge behaviour and alterations in dispersal, mor-tality, and species interactions. Changes in movement behaviour as a re-sult of fragmentation is well-documented [CDO07, FCC99, HC03, HC06,MAD10, Ova04, SW04, SB03] and can have important implications forsurvival of species and the outcomes of competition or predator-prey in-teractions [FCC99, GTLS]. For the purposes of this chapter, we define alandscape as a combination of focal habitat patches (where the individualmay reproduce and/or find food) and matrix habitat (the habitat in whichthe individual cannot reproduce or find food). Matrix habitat is definedas high quality if the individual can disperse throughout this habitat witha low probability of death due to such pressures as predation. We defineedges as boundaries between the focal habitat patch and the matrix habi-tat. The degree of habitat fragmentation effects depends heavily upon thesize of individual focal and matrix habitat patches and the quality of thematrix habitat [HC03, HC06, ST09, ST12]. Fahrig [Fah07] predicted thatpopulations are differentially affected by habitat fragmentation as a result oftheir varying movement rates and the original quality of their matrix habi-tat. For instance, populations with large movement probabilities, that resultfrom high-quality matrix habitat, would be highly vulnerable to focal habi-tat loss or reductions in matrix quality. In contrast, a population with low323.1. Introductionmovement, as a result of low-quality matrix habitat, would be less affectedby focal habitat loss or lower quality matrix habitat. This population would,however, exhibit a strong negative response to lower colonization (success-ful movement to a new focal habitat patch with successful establishment)probability.One important potential change in movement behaviour is the willing-ness of a species to cross habitat boundaries/edges [CC82, Jon77, RH98,SC01, WSM04]. This is often referred to as the permeability of the bound-ary and describes the probability that the species will cross from the focalhabitat patch to the matrix habitat. Intuitively, a boundary with high per-meability has a high probability that the species will cross the edge, whilea low probability of crossing an edge results from an edge with low perme-ability.There are several characteristics of movement behaviour that can be al-tered by increased habitat fragmentation. Chapman et al. [CDO07] foundthat the dispersal paths of leaf beetles in matrix habitats had smaller turn-ing angles and the movement event lasted longer than movement eventsin the focal habitat. The authors also found that if the boundaries wereless permeable, colonization was dramatically reduced. Haynes and Cronin[HC03, HC06] found that populations of planthoppers have different move-ment speeds in brome grass and mud matrix habitats. Additionally, plan-thoppers developed different responses to the edges of the different types ofmatrix habitat. That is, the permeability of the edge was different betweena focal habitat/brome grass matrix edge and a host habitat/mud matrixedge. Studies by Jones [Jon77] and Courtney and Courtney [CC82] foundthat movement and oviposition of butterflies results in clumped egg distri-bution at the edges of focal habitat patches.The Mountain Pine Beetle (MPB) has been extensively studied due tothe economic impact it has on the forestry industries in Western Canadaand the United States [SW06, Hug02]. Nonetheless, the effect of habi-tat fragmentation on the spatial spread of MPB is not well understood[BAL93, Hug02, RWNW]. For MPB, susceptible habitat, where suscepti-ble lodgepole pine trees for MPB attack and reproduction are found, are333.1. Introductiondefined as focal habitat. Theoretical work by Hughes [Hug02] predictedthat habitat fragmentation can slow the spread of MPB if the focal habitatpatch sizes are small enough. In contrast, empirical studies by Robertsonet al. [RWNW] found that the spread of MPB over small scales increasedwith increasing habitat fragmentation. Understanding the spatial spread ofMPB in response to habitat fragmentation is crucial as many current man-agement strategies, such as prescribed burning and removal of susceptiblestands of lodgepole pine, rely on the assumption that MPB spread can beslowed by the introduction of matrix habitat and reduction in the amountof focal habitat. Further studies by Reid [Rei09] found that decreased tor-tuosity in the flight path of MPB can reduce the time required to cross amatrix habitat between two good quality habitats. This result would implythat MPB will move more quickly in matrix habitat, which would violatethe assumption underlying present management in that the introduction ofmatrix habitat will slow the spatial spread of MPB.We want to understand how MPB dispersal in fragmented habitats isreflected in trap catches on the local (stand) scale. In particular, we areinterested if the MPB distribution at the beginning of the attack period ismore uniform over the landscape or localized to the emergence sites. Addi-tionally, we want to understand the rates of diffusion of MPB in fragmentedhabitat. In this chapter we report on an experimental study and use theo-retical work to help understand the effects of habitat fragmentation on themovement of MPB. The experimental study consists of trap catch data forMPB in different domains of fragmented habitat. Specifically, we recordpheromone trap catches on either side of the edge of a fragmented habitatfor different combinations of clearcut, thinned, and dense forest regions inSection 3.2. To help understand the results of this experiment with relationto MPB movement we develop a spatially explicit reaction-diffusion modelwith chemotaxis for MPB in Section 3.3. Simulations of this model are runin Section 3.4, using the numerical scheme outlined in Section 3.5. Section3.6 outlines our simulation results and Section 3.7 ties the empirical andtheoretical results of our study together. Discussion and Future work isreserved for section 3.8.343.2. Experimental3.2 ExperimentalBiological experiments were conducted by Mary Reid in Parson, BritishColumbia over 2007 and 2008. Two fragmented stands were chosen withclearcut harvested regions, and two fragmented stands with thinned har-vested regions were also selected. The dense stands had approximately 2000stems/ha, while thinned stands had only 1000 stems/ha of susceptible lodge-pole pine trees. Statistical analysis of the trap catches across the differentfragmented regions were performed to understand the effect of harvestedregions on the dispersal of MPB. Examining the MPB catches in the baittraps in each landscape, Dr. Reid found the following experimental results: astatistically significant difference was found between the clearcut and denseregion, but no significant difference was found between the thinned anddense region. Comparing across studies Dr. Reid found that the clearcutand thinned bait trap catches were not significantly different. Finally, in agiven habitat (harvested or unharvested) the 20 m and 90 m traps did notdiffer significantly in terms of the MPB trap catches. Further details of theexperiment can be found in Strohm et al. [SRT].3.3 Mathematical ModelTo mathematically analyze MPB movement across habitat edges, weneeded a spatially explicit model which describes the interaction betweenMPB, lodgepole pine trees, MPB pheromones, and pheromone baits. Ourmodel is based on previous theoretical work [WP97] examining MPB disper-sal. This model is a continuous set of Partial Differential Equations (PDEs)describing the period of flight, emergence, and attack. Our five variablesare P - the density of flying MPB, Q- the density of nesting MPB, J - thedensity of pre-adult MPB (beetles not yet emerged from the bole of thetree), A - the concentration of beetle pheromone, and C - the concentrationof tree kairomones. In addition, S - is a measure of the relative density of353.3. Mathematical Modellodgepole pine trees. The model equations are:?P?t=diffusion? ?? ??p?2P ?chemotaxis? ?? ??[(?c?C + ?ab0 ?Ab0 +A/b1?A)P]?death?????pP ?nesting? ?? ?g(P,C)S+emergence?????J ,(3.1a)?Q?t=death? ?? ???qQ+new nesters? ?? ?g(P,C)S , (3.1b)?C?t=diffusion? ?? ??c?2C +synthesis????a2S ?degradation?????cC +baiting????? , (3.1c)?A?t=diffusion? ?? ??a?2A +synthesis????a1Q ?degradation?????aA , (3.1d)?J?t=emergence??????J . (3.1e)where,g(P,C) =landing?????Pmass attack? ?? ?P 2P 2 + (kp)2.All parameters are defined and given specified values in Table 3.1.This model describes MPB emergence and dispersal during a single flightperiod, that is, over a single summer season. Note that we are not trying tomodel the population spread across multiple years.A submodel derived from (3.1) has recently been used to elucidate thespacing observed between MPB attack locations in a single year (Chapter 2).These attack locations appear as spots in the data. The model was shown toreproduce the correct spatial frequency of attacks present in Aireal DetectionSurvey data from the Sawtooth National Recreation Area (Chapter 2). Thiswork provides some validation of the submodel and the parameter valuesused. For this chapter, we build on our earlier work by adding several new363.3. Mathematical Modelcomponents to the submodel, ariving at (3.1).The most important departure from the original submodel is the additionof a spatially-heterogeneous landscape. In our previous work (Chapter 2),we were looking for the pattern of MPB attack on a spatially homogeneouslandscape. Here we are interested in studying the effects of habitat edgeson MPB dispersal, and so we include heterogeneity in the distribution oftrees, S(x). As a result, the interaction between MPB and tree-producedkairomone becomes important necessitating the addition of (3.1c). Thekairomone, C, diffuses, is synthesized by susceptible trees, and degradesover time. We also add a baiting term, ?, that simulates the productionof pheromone by traps. We included this term in the kairomone equationand not in the pheromone equation to differentiate it from the pheromonenaturally produced by the MPB. Kairomone, also appears in (3.1a), as flyingMPB have a chemotactic attraction to kairomone and will be attracted upthe gradient of C [WP97].The indicator of susceptible tree density, S(x), varies between 0 (nosusceptible trees) and 1 (landscape is entirely susceptible). Next, we add thevariable J(x, t), for the density of emerging MPB at any point in time andspace during the summer. Most emergence occurs over a brief time periodduring the flight season, and can extend over longer times under appropriateclimactic conditions [SW06]. We are interested in both scenarios, so we allowthe emergence rate, ?(t), to vary with time. The timing of emergence iscrucial to the successful attack of MPB [GPLB04] and so we expect temporalvariations to be important in the behaviour of the model. For the simulationsshown in this chapter we chose ?(t) such that 99.9% of MPB emerge withinthe first 20 flight hours of the season. We simulated other choices of ?(t)and discuss the results in Section 3.8.The selection of parameter values is informed by previous theoretical andempirical work [BPBL00, RB83]. We chose the parameters a2 and ? basedon the results of numerical simulations of the model in order to simulate thecorrect biological behaviour. A detailed discussion of parameter values forMPB dispersal can be found in Strohm et al. (Chapter 2.3). The valuesused in this study are listed in Table 3.1.373.3. Mathematical ModelTable 3.1: Table of parameter values for the dimensional model (3.1).Parameter Description Units Value?p diffusion of flying MPB hafh 1?c diffusion of host volatiles (kairomones) hafh 0.648?a diffusion of beetle pheromones hafh 0.648?c kairomone attractiveness ha2?g?fh 0.8?a beetle pheromone attractiveness ha2?g?fh 5.7b0 concentration of pheromones at which dis-sipation occurs?gha 5.4b1 concentration at which pheromone is sat-uratedn/a 1? random landing rate of flying MPB hatrees?fh 0.16kp flying beetle density required for 50 per-cent success of mass attackmpbha 250000?p death rate of flying MPB fh?1 0.014?q death rate of nesting MPB fh?1 0.001?c degradation of tree kairomones fh?1 180?a degradation of beetle pheromone fh?1 180?(t) emergence rate of pre-adult beetles fh?1 0.345a1 rate of pheromone increase due to nestingMPB?gfh?mpb 0.02a2 rate of kairomone increase due to suscep-tible trees?gfh?trees 0.02? rate of kairomonne production by traps ?gha?fh 1458383.3. Mathematical Model3.3.1 Non-dimensionalization of ModelWe non-dimensionalize the model for ease of mathematical analysis.We choose the following non-dimensionalizations for our variables:Q =a1b0?aQ, P =a1b0?aP, C =?c?cC,S =??pS, A =1b0A, J =a1b0?aJ,t = ?at, (x, y) =??a?a(x, y).Many of the scalings chosen reflect dynamics of the MPB pheromone. Q andA are scaled such that they represent density of MPB and concentration ofMPB pheromone required to alter the MPB response to pheromone from at-tractive to repulsive. P and J are scaled by the same parameter combinationas Q since they all represent densities of MPB. The tree kairomone, C, isscaled by the ratio of diffusive movement of kairomone to the attractivenessof kairomone to MPB. The susceptible tree density, S, is scaled by the ratioof the average time until natural death of flying MPB divided against theaverage time before flying MPB randomly land on a susceptible tree. Timeand space are scaled by the average time and distance that the pheromonecan spread before degradation.Using these variable scalings we obtain the model equations,?P?t=?p?a?2P ? ?c?a?(P?C)? ?ab0?a?(1?A1 +A/b1P?A)? ?p?aP ? ?p?aP3P2+ (kpa1b0?a )2S +??aJ,(3.2a)?Q?t=? ?q?aQ+?p?aP3P2+ (kpa1b0?a )2S, (3.2b)?C?t=?c?a?2C + ?pa2?c?a??cS ? ?c?aC +?c?a?c?, (3.2c)393.3. Mathematical Model?A?t=?2A+Q?A, (3.2d)?J?t=? ??aJ. (3.2e)Some of the parameter values shown in Table 3.1 differ by orders ofmagnitude, and so processes in the model act over fast and slow time scales[PMW98]. We can explicitly highlight these different timescales by defininga scaling parameter. We choose as our parameter the ratio of average timebefore degradation of pheromone and the average time prior to death flyingMPB, =?p?a.Using this scaling parameter, we define the dimensionless parameters?p =?p?a, ?a =?ab0?a, kp =kpa1b0?a,?c =?c?a, ? =a2?c??c, ? =?c?a?c?,?c =?c?a, ? =??a, ?q =?q?a.The values of the non-dimensional parameters can be found in Table 3.2.Substituting these non-dimensional parameters into (3.2) we arrive at thenon-dimensional model (3.3). equations:?P?t=?p?2P ? ?c?(P?C)? ?a?(1?A1 +A/b1P?A)? P ?  P3P 2 + k2pS + ?J,(3.3a)?Q?t=? ?qQ+ P 3P 2 + k2pS, (3.3b)?C?t=?c?2C + ?S ? ?cC + ?, (3.3c)403.4. Simulation SetupTable 3.2: Table of parameter values for the non-dimensional model (3.3).Parameter Value 0.0000778?p 1.54?c 1?a 47.5? 0.154kp 5.14?q 0.0714?c 1? 24.6? 10?A?t=?2A+Q?A, (3.3d)?J?t=? ?J. (3.3e)Note that we have dropped the overbars on the variables and equations forconvenience.3.4 Simulation SetupThe goal of the model is to help illuminate the results of the baitingexperiment performed in Section 3.2.The simulation was set up (Figure 3.1 according to the landscapes andtrap locations defined in the experimental portion of our study. We simu-late a two-dimensional landscape that is heterogeneous and changes with apiece-wise constant manner. For x < 0, the landscape is unharvested andpopulated with a uniform density of susceptible trees. We define this unhar-vested region as the dense region. For x > 0 the landscape is harvested andpopulated with susceptible trees at a thinned density, or with no trees at all(clearcut). We examined MPB movement behaviour in these two types oflandscapes with pheromone traps placed at -90, -20, 20, and 90 m from theboundary between the dense and harvested regions.413.4. Simulation Setup?90 m ?20 m 20 m 90 my=0Figure 3.1: This plot shows a heterogeneous habitat with 4 bait sites. Thelight region (left) is uniformly dense suitable habitat and the dark region(right) is either a clearcut or a thinned habitat. The bait sites are at?20,?90 and are marked by circles.423.4. Simulation SetupSeveral MPB parameters were varied across the simulations: MPB emer-gence distribution, diffusion rate, and finally MPB population density. Twoemergence distributions were investigated; (1) uniform emergence and (2)emergence only in the dense region. The uniform case arises from the hy-pothesis that flying MPB are first distributed fairly uniformly over the do-main after emergence, before beginning the process of focussing on targettrees (personal communication, Jim Powell). The other initial conditionarises from the alternative hypothesis that beetles remain close to theirsource trees after emergence and so will initially only be found in the denseregion. Different diffusion rates were simulated as there is evidence thatMPB may move more quickly in the harvested region due to the reduceddensity of trees and therefore less tortuous flight paths [Rei09]. Finally, arange of population densities was simulated to understand how the baselineMPB population density affects MPB bait trap catches.At a given time and position, the bait trap catches were calculated byintegrating the flying MPB and nesting MPB density in a 6 m radius aroundthe bait site. We then assumed that 110 of flying MPB are caught in the trap,while all of nesting MPB are caught in the trap. Since we are interested intrap catches over the summer, we averaged the trap catches over time. Thiswill give us an expected qualitative population densities of MPB in each trap.We simulated a very synchronous emergence and bait trap catches weretaken over the period before the MPB population decreased dramaticallydue to death and lack of new emerging MPB. The complex ecology of MPBmakes it difficult to assess the MPB population size in each trap, so ourmethod provides a qualitative understanding of the trap catches.The simulation landscape was chosen to have a size of 100 m in they-direction and 2.82 km in the x-direction (Figure 3.1). Other than thetraps at the center of the vertical landscape, the rest of the landscape isuniform with respect to y. The domain was chosen very large in the x-direction to avoid the effects of boundary conditions on the bait trap catches.The boundary conditions used for this simulation were chosen as zero-flux(reflecting) boundary conditions, so no MPB were lost at the boundary. Weassume that the simulation is embedded in a large domain that is periodic433.5. Numerical Method for Simulationsin y, so we need only simulate a single strip.3.5 Numerical Method for SimulationsThe MPB equations (3.3) are a set of Reaction-Diffusion-Advectionequations. There are optimal schemes for separate Advection, Diffusion andReaction equations, but finding a scheme that can solve all three simul-taneously is difficult [LeV07]. Problems that can arise include finding anappropriate timestep size for stability and grid-scale oscillations in the solu-tions [Tys96]. Therefore, we use the technique of Strang Splitting so that wecan use separate schemes for each type of equation and join them togetherusing this fractional step method. Note that second-order accuracy in solu-tions is obtained, assuming each individual scheme is second-order accurate[LeV07]. We rewrite the system of equations asut = D(u) +A(u) +R(u), (3.4)where D(u) is the diffusion step, A(u) is the advection step, and R(u) isthe reaction step. The splitting technique involves taking a half-step ofadvection, then a half-step of diffusion, followed by a full-step of reaction.After this it proceeds in the the reverse order, with a half-step of diffusion,and then a final half-step of advection. Thus, there is a full-step of diffusion,reaction and advection over the splitting technique. Mathematically, wewriteun+1 = A(?t2)D(?t2)R(?t)D(?t2)A(?t2)un. (3.5)The method of Strang splitting, and the specific order of the diffusion,advection, and reaction steps in (3.5) follows work by Tyson et al [TSL00]and was chosen for its accuracy and stability for PDE systems with chemo-taxis. The advection step is solved using a software package called CLAW-PACK (Conservation LAWs PACKage) [LeV97]. This software uses wavepropagation methods to solve multidimensional hyperbolic systems of partialdifferential equations. The algorithm solves the Riemann problem for the443.5. Numerical Method for Simulationswaves, which requires computation of the flux-difference splitting betweeneach cell. Finally, flux limiters are used to obtain second-order accuracy andmaintain high resolution.The diffusion step is solved using an ADI (Alternating Direction ImplicitMethod, (3.6)) in space and a TR-BDF2 (Trapezoidal, Backwards Differenti-ation Formula of order 2, (3.7)) method in time. The TR-BDF2 method wasused instead of a standard second-order Crank-Nicholson Method since grid-scale oscillations have been found to occur in similar problems [TSL00]. Thereason for the grid-scale oscillations was that the Crank-Nicholson methodis only A-stable, whiled the TR-BDF2 method is both A-stable and L-stable[LeV07]. Finally, the reaction step is solved by using the standard fourth-order Runge-Kutta method, (3.8). The combination of these methods fordiffusion, advection and reaction using Strang Splitting (3.5) has been usedsuccessfully to solve other reaction-diffusion chemotaxis equations [Tys96].The ADI method uses equations [LeV07],U?ij =Unij +k2(D2yUnij +D2xU?ij), (3.6a)Un+1ij =U?ij +k2(D2xU?ij +D2yUn+1ij ). (3.6b)The TR-BDF2 obeys the equations [LeV07],U? =Un +k4(f(Un) + f(U?)), (3.7a)Un+1 =13(4U? ? Un + kf(Un+1)). (3.7b)The Fourth-Order Runge-Kutta Equations are [LeV07],Y1 =Un, (3.8a)Y2 =Un +12kf(Y1, tn), (3.8b)Y3 =Un +12kf(Y2, tn +k2), (3.8c)453.5. Numerical Method for SimulationsY4 =Un + kf(Y3, tn +k2), (3.8d)Un+1 =Un +k6(f(Y1, tn) + 2f(Y2, tn +k2)+ 2f(Y3, tn +k2)+ f(Y4, tn + k)),(3.8e)In (3.6), (3.7), and (3.8), k represents the time step, D2x and D2y, is thesecond derivative in space in the x and y direction, i and j are the x andy coordinates, f is the function such that u?(t) = f(u(t), t), and tn is thecurrent simulation time.Simulations were run in Fortran on the SARAHS cluster at UBC, Okana-gan Campus.To test the method used to solve the reaction step of the PDEs, wesimulated the system using an explicit and implicit method. The implicitmethod, a Backwards Euler Method, and the explicit method, a fourth-orderRunge Kutta routine, converged to the same solution. Therefore, we use theexplicit method to enhance the speed of computations.After consultation with Dr. Randy LeVeque, the author of Clawpack,we decided to check convergence of the solution using grid and time scale-tests on initial stages of pattern formation. It was recommended to usefull problem instead of a subset of the reaction, advection, and diffusionpieces of the equation, because many problems only arise when all portionsof equation are put together. TheThe convergence to the correct solution in nesting and flying MPB fordecreasing space (and time) step sizes is displayed in Figure 3.2. We choseour space step, h to decrease with the time step, k, according to the equationk = 4h [LeV07]. The log-log plot of the absolute value of the error withdecreasing step size is shown in Figure 3.3. It can be seen that the error inthe simulation decreases as step size decreases and that the convergence tothe solution occurs for small space step sizes.To numerically compute the error we run the simulation over a shortperiod of time with the initial pattern formation occuring. Since we do nothave the true solution set of our equations (3.3) we used the method of463.5. Numerical Method for Simulations?20 ?15 ?10 ?5 0 5 10 15 2001000200030004000500060007000?20 ?15 ?10 ?5 0 5 10 15 2000.511.522.533.544.5 x 10?30.031250.06250.1250.250.5124Flying MPB Nesting MPBPopulation DensitySpaceFigure 3.2: The graphs illustrate the convergence of solutions of flying (left)and nesting (right) MPB for decreasing space step sizes, h. The differentspace step sizes are shown in the legend.10?2 10?1 10010010110210310?2 10?1 10010?710?610?510?410?3Log(abs(E))Log(h)Flying MPB Nesting MPBFigure 3.3: Log-log error plots of the absolute value of the error againstspace step size, h, for flying (left) and nesting (right) MPB.473.5. Numerical Method for SimulationsTable 3.3: Order of accuracy for nesting, q, and flying, p, MPB solutionswith decreasing space step size, h. log2(R) is the order of accuracy for agiven space step size.Variable log2(R(4)) log2(R(2)) log2(R(1))q 0.8207 1.5933 1.9301p 0.4353 1.0831 1.2671Variable log2(R(0.5)) log2(R(0.25)) log2(R(0.125))q 1.8407 1.8220 1.8824p 1.2716 1.3210 1.6342Variable mq 1.6482p 1.1687computing errors for numerical solutions (3.9) [LeV07]. That is, for smallstep sizes in space and time we assume thatE(h) = Chm + o(hm), (3.9a)E(h) ? Chm, as h? 0, (3.9b)E(h2)? C(h2)m, (3.9c)R(h) =E(h)E(h/2)= 2m, (3.9d)m ? log2(R(h)). (3.9e)where, E(h) is the error at space step size h, C is the error constant, R isa ratio of errors, and m is the order of accuracy.This allows us to produce Table 3.3, which numerically summarizes theorder of accuracy in our solutions (for flying and nesting MPB) as we de-crease the space step size, h.Based on the order of accuracy from Table 3.3, we find that our numer-ical simulations are above order 1.5 accuracy for nesting MPB density, andgreater than order 1 accuracy for flying MPB density. These two quanti-ties have the non-linear and chemotaxis portions of the code and thus arethe most important variables to calculate error and ensure accuracy of thesolutions.483.6. Model Results3.6 Model ResultsFor our simulations, we used emergence population densities, J0,varying between 31,250 MPB/ha and 250,000 MPB/ha. Assuming thata single source tree can produce on the order of 10,000 MPB per flight sea-son [PB09], these densities correspond to between 3 and 25 source trees perha. Diffusion rates ranging from 1 to 2 ha/fh were considered, consistentwith the estimation in Powell et al. [BPBL00], and will be used to test thehypothesis of faster diffusion in the harvested region.All variables were initially set to zero, except for J , which gives the initialdensity of MPB emerging from trees and becoming flying MPB. We ransimulations of our mathematical model (3.3) with various initial populationdensities and diffusion rates. We tracked the spatial dynamics of flying (P )and nesting (Q) MPB along a linear transect (y = 0, position of bait traps,Figure 3.1) at small and large intial densities of MPB for four cases: clearcutand thinned harvested region with dense forest initial conditions (Figures3.4 and 3.5, respectively), and clearcut and thinned harvested region withuniform initial conditions (Figures 3.8 and 3.9, respectively).When MPB emerged only in the dense forest, the MPB had similar dy-namics in both clearcut (Figure 3.4) and thinned (Figure 3.5) harvested re-gions. Both figures have, in general, higher densities of MPB in the forestedregion than in the harvested region. This difference is due to the fact thatno MPB emerged from this region, they only emerged from the unharvestedregion. This diffusion of MPB into the harvested region creates a decreasingpopulation density curve from the unharvested (left) to harvested (right)region. Therefore, we expect the number of MPB to decrease from left toright (-90 m, -20 m, 20 m, 90 m). In general, this hypothesis is upheld bythe simulation results, with evident spikes in both the nesting and flyingMPB at the location of the bait traps. At high densities, another dynamicemerges, where there is sufficient MPB to begin a successful attack. In thissimulation the MPB attack and nest to the left of the trap at -90 m. Thissuccessful attack pulls other MPB toward it and thus we see the change inthe slope of the flying MPB profile in the plots at 4.8 and 6.4 flight hours493.6. Model Results(fh). The last plot at 8 fh, shows the progression of the successful attackreaching the -90 m trap, and MPB are beginning to attack near this trap ata high density. This pulls MPB not already nesting or caught in traps fromsurrounding regions of -90 m towards that trap, lowering catches in nearbytraps.The major difference between Figures 3.4 and 3.5 is the right hand (x >0) tail in the distribution of nesting MPB. In the clearcut (Figure 3.4), thepopulation of nesting MPB is essentially zero, while it is small but positiveand decaying with x in the thinned region (Figure 3.5). Compare the secondand fourth columns of Figure 3.4 against Figure 3.5. The positive density ofnesting MPB comes from nesting MPB who have diffused into the harvestedregion from the unharvested region.The bait trap catches corresponding to the MPB distributions in Fig-ures 3.4 and 3.5 are shown in Figure 3.6. There is essentially no differencebetween the left and right graphs. This is because the number of MPB whodiffuse from the unharvested region to the harvested region is small enoughthat nesting in the thinned region does not cause a significant increase inthe density of MPB trapped at that point.The bait trap catches in both cases are ordered, with the number oftrapped MPB decreasing from the -90 m to the 90 m trap along the y = 0transect. This makes sense since diffusion pulls MPB away from the sourcesin the unharvested region and into the harvested region. This creates a de-clining density from left to right along the y = 0 transect. It is interestingto note that the attracting region for the -20 m and 20 m bait traps overlap(Figure 3.7) and so it is possible that these two traps could interfere witheach other and disrupt the ordering of the trap catches expected from diffu-sion alone. We do not observe this, however, and so the trap catch patternis dominated by the effects of diffusion.As the density of MPB increases from left to right in each plot in Figure3.6, there are two changes to the graph. First, since a larger density of MPBare emerging, more make it to the trap at 90 m (in the harvested region).Consequently, the difference between the trap catches at 20 m and 90 mdecreases. The second major change is that there is a major aggregation503.6. Model Results?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 100010002000300040005000600070008000?100 ?50 0 50 1004000500060007000800090001000011000120001300014000?100 ?50 0 50 1000.70.80.911.11.21.31.41.51.61.7 x 104?100 ?50 0 50 1000.80.911.11.21.31.41.51.61.71.8 x 104?100 ?50 0 50 1000.911.11.21.31.41.51.61.71.8 x 104?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 10000.050.10.150.20.250.30.35?100 ?50 0 50 100012345678?100 ?50 0 50 1000510152025?100 ?50 0 50 10005101520253035404550?100 ?50 0 50 10001020304050607080?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1000.511.522.533.544.555.5 x 104?100 ?50 0 50 1002345678910 x 104?100 ?50 0 50 1001.41.51.61.71.81.922.12.2 x 104?100 ?50 0 50 1002000300040005000600070008000900010000?100 ?50 0 50 10000.511.522.5 x 106?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1000102030405060708090?100 ?50 0 50 10005001000150020002500?100 ?50 0 50 1000500100015002000250030003500?100 ?50 0 50 1000500100015002000250030003500?100 ?50 0 50 100012345678910 x 104t=0 fht=3.2 fht=4.8 fht=8 fht=1.6 fht=6.4 fhFlying MPB Nesting MPBSpace (m)Flying MPB Nesting MPB High PopulationLow PopulationPopulation Size (MPB/ha)Figure 3.4: Spatial population densities of flying and nesting MPB in alandscape that has a clearcut harvested region and dense forest intial condi-tions. The spatial plots are shown over increasing time simulation in MPBflight hours (fh). The two columns on the left depict the MPB at a lowinitial population density (J = 31, 250 MPB/ha, approximately 3 sourcetrees per ha) and the two right columns illustrate high initial populationdensities (J = 218, 700 MPB/ha, approximately 22 source trees per ha).The simulations shown assume diffusion is 1 ha/fh.513.6. Model Results?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 100010002000300040005000600070008000?100 ?50 0 50 1000.70.80.911.11.21.31.41.51.61.7 x 104?100 ?50 0 50 1004000500060007000800090001000011000120001300014000?100 ?50 0 50 1000.80.911.11.21.31.41.51.61.71.8 x 104?100 ?50 0 50 1000.911.11.21.31.41.51.61.71.8 x 104?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 10000.050.10.150.20.250.30.35?100 ?50 0 50 100012345678?100 ?50 0 50 1000510152025?100 ?50 0 50 10005101520253035404550?100 ?50 0 50 10001020304050607080?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1000.511.522.533.544.555.5 x 104?100 ?50 0 50 1002345678910 x 104?100 ?50 0 50 1001.51.61.71.81.922.12.2 x 104?100 ?50 0 50 1002000300040005000600070008000900010000?100 ?50 0 50 10000.511.522.5 x 106?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1000102030405060708090?100 ?50 0 50 10005001000150020002500?100 ?50 0 50 1000500100015002000250030003500?100 ?50 0 50 1000500100015002000250030003500?100 ?50 0 50 100012345678910 x 104Flying MPB Nesting MPBSpace (m)Flying MPB Nesting MPB High PopulationLow PopulationPopulation Size (MPB/ha)t=0 fht=3.2 fht=6.4 fht=1.6 fht=4.8 fht=8 fhFigure 3.5: Spatial population densities of flying and nesting MPB in alandscape that has a thinned harvested region and dense forest intial condi-tions. The spatial plots are shown over increasing time simulation in MPBflight hours (fh). The two columns on the left depict the MPB at a lowinitial population density (J = 31, 250 MPB/ha, approximately 3 sourcetrees per ha) and the two right columns illustrate high initial populationdensities (J = 218, 700 MPB/ha, approximately 22 source trees per ha).The simulations shown assume diffusion is 1 ha/fh.523.6. Model Resultsevent left of the -90 m trap at high densities of MPB. This is shown in thedistribution of flying and nesting MPB at t = 8 fh for high MPB populationdensities in Figures 3.4 and 3.5. The large density of nesting MPB drawsin nearby flying MPB populations and there is a large increase in the flyingand nesting MPB near the -90 m trap. This major aggregation event is dueto the self-focusing nature of the MPB biology [WP97]. Once a high enoughdensity of flying MPB begin to nest, the nesting MPB produce pheromone,which draws more flying MPB towards that site. Once the site has reacheda critical density, the pheromone becomes repulsive and approaching MPBnest instead in nearby trees. This process is the cause of the major increasein the trap catches at -90 m for densities of MPB above 1.8? 105 per ha.It is unclear whether the MPB in a given landscape are largely from localsources, or if a significant proportion of beetles drift in from elsewhere vialong-distance dispersal. In the first case, we expect MPB to emerge onlyin the unharvested region (dense forest ICs); this was the case simulated inFigures 3.4 and 3.5. In the second case, we expect the initial distributionof MPB to be more uniform, assuming that MPB arrived/emerged in auniform manner in both the dense forested region (unharvested) and theharvested region. The resulting spatial population densities of flying andnesting MPB over time at a high and low population density are displayedin Figures 3.8 and 3.9. Another mechansim which could give rise to a moreuniform emergence density is if there is a sexual maturation time involvedbefore MPB begin to attack susceptible trees and become responsive tochemical signals [HFSL06].In this scenario, instead of the strongly sloped MPB density profile seenin Figures 3.4 and 3.5, the flying MPB have a uniform density, with spikesat the bait locations. If there are trees to nest in at these locations, thesebait traps will also act as the initial positions of MPB aggregation on thesusceptible landscape. The spikes in nesting MPB density are evident inFigures 3.8 and 3.9. At low densities, the nesting MPB change the flyingMPB density as expected. In the clearcut case (Figure 3.8 the nestingMPB produce pheromone which attracts more flying MPB and thus thedensity of flying MPB in the unharvested region is always larger than in533.6. Model Results0 0.5 1 1.5 2 2.5x 1050510152025303540450 0.5 1 1.5 2 2.5x 1050510152025303540450 0.5 1 1.5 2 2.5x 105024681012140 0.5 1 1.5 2 2.5x 10502468101214Average Trap CatchesInitial Population Size (MPB/ha)Regular DiffusionFaster DiffusionClearcut Thinned?90 m?20 m20 m90 mFigure 3.6: Average trap catches in the bait sites at ?20,?90 m with aninitial condition of MPB in dense forested region only. The bait trap catchesare shown for heterogeneous landscapes with the harvested region eitherclearcut or thinned. The bottom two plots show the trap catches expectedif MPB diffuse faster in harvested regions.543.6. Model Results?100 ?50 0 50 10000.020.040.060.080.10.120.14Space(m)Bait Pheromone Concentration (?g/ha)Figure 3.7: This graph represents the bait pheromone profile at the begin-ning of a simulation. It can be seen that the baits are at ? 90 and ? 20 m.The profiles at ? 20 m overlap, therefore sharing an attraction region. Thesimulations shown assume that the diffusion rate is 1 ha/fh.553.6. Model Resultsthe harvested region. In the thinned case (Figure 3.9) the MPB are ableto nest in the harvested region, so the difference in MPB between the tworegions is smaller. The smaller density of susceptible trees in the harvestedregion results in a slower nesting rate, and thus more MPB are found in theunharvested region.At higher densities, the thinned and clearcut results are similar for thefirst 3.2 flight hours. First the nesting occurs exclusively (clearcut) or mainly(thinned) in the unharvested region. Notice that the density of MPB atthe -20 m site is larger than the density of MPB at the -90 m site. Thisbehaviour is due to the nature of the MPB attraction to pheromone. Sincethe landscape changes at 0, the change in MPB pheromone concentrationsharply decreases at this point. Since the MPB respond to the gradientof MPB pheromone (chemotaxis term), and this gradient is largest at theborder between the harvested and unharvested region, the attraction to theunharvested region is maximum near this point. Figure 3.10 shows the MPBpheromone profile (and the resulting chemotactic attraction). Therefore,since the emergence is uniform, MPB will begin to nest in greater numbers atthe -20 m site versus the -90 m site. Once the nesting at the -20 m site begins,the self-focusing behaviour of the MPB strongly affects the spatial dynamicsof the population. At the initial aggregation site, nesting beetles releasepheromone, increasing the gradient there and attracting more flying MPB.This sets up a positive feedback cycle that intensifies the initial aggregation.Once enough MPB have nested and increased the pheromone to a criticallevel the pheromone becomes repulsive and signals for MPB nest in nearbytrees. The nesting density of MPB then broadens to a larger spatial regionof attack.Once repulsion of MPB has broadened the spatial region of attack (highdensity, 4.8 fh) the shape of the nesting MPB curve is determined by thetype of harvesting (thinned or clearcut). If the spatial region is clearcut(Figure 3.8), MPB cannot nest there, so there is a build-up of MPB at theborder between the harvested and unharvested regions. At this point theMPB pheromone gradient is very steep, and thus very attractive, so there isa build-up of MPB at the border of the harvested and unharvested region.563.6. Model ResultsThis leads to a very sharp peak of nesting and flying MPB density near thispoint. To the left of this peak there is a bump in the uniform nesting MPBdensity at the intial site of attraction. This is due to the larger densityof MPB needed to produce the pheromone required before the pheromonesignal became repulsive. For nearby trees, there is already a large concen-tration of pheromone present, so it takes a lower density of nesting MPBto produce the needed pheromone for the pheromone to become repulsive.It can be seen there is also a small increase at -90 m, due to the increasedattractiveness of the bait site from bait pheromone.If the region is thinned (Figure 3.9), the spatial dynamics are a little dif-ferent. Similar to the clearcut region, the flying MPB build up on the edgesof the nesting MPB curve and there are peaks in the nesting MPB densityat the bait sites due to steep pheromone gradients there. The largest peakof nesting MPB occurs at the initial bait site of -20 m (3.2 fh, high density).The behaviour that differentiates the thinned case from the clearcut one isthat the nesting MPB density spreads into the harvested region (since thereare trees to nest in). There is an interesting spike at the border betweenthe harvested and unharvested region (0 m). This is due to the difference inkairomone production and the difference in nesting rate across the boundarybetween the two regions.We calculated bait trap catches based on integration of these curves, theresults are shown in Figure 3.11. Differences between thinned and clearcutare more significant if MPB emergence occurs uniformly (Figure 3.11) in-stead of only in the dense forest (Figure 3.6). The uniform emergence ensuresa higher density of MPB in the harvested regions. This results in a largerproportion to be caught in the traps of the thinned harvested region. Fi-nally, the traps at ? 90 and ? 20 m display more differences with uniforminitial conditions than dense initial conditions.An interesting note is that the trap at 20 m actually has less MPB thanthe trap at 90 m (clearcut). This is due to the MPB pheromone at the -20m trap attracting MPB flying near the 20 m trap (see Figure 3.10). Thisresult is due to the overlap in MPB pheromone from the traps at ?20 m.The same process also explains how the trap in the thinned region, at573.6. Model Results?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1008100820083008400850086008700?100 ?50 0 50 1001.821.841.861.881.91.921.941.96 x 104?100 ?50 0 50 1002.32.352.42.452.52.55 x 104?100 ?50 0 50 1002.52.552.62.652.72.752.82.85 x 104?100 ?50 0 50 1002.552.62.652.72.752.82.852.92.9533.05 x 104?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 10000.050.10.150.20.250.30.350.40.45?100 ?50 0 50 100024681012141618?100 ?50 0 50 10001020304050607080?100 ?50 0 50 100020406080100120140160180?100 ?50 0 50 100050100150200250300?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1005.65.75.85.966.16.26.3 x 104?100 ?50 0 50 10000.511.522.53 x 106?100 ?50 0 50 10000.511.522.533.544.5 x 106?100 ?50 0 50 1000123456 x 106?100 ?50 0 50 10001234567 x 106?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 100050100150?100 ?50 0 50 1000246810121416 x 104?100 ?50 0 50 10000.511.522.53 x 105?100 ?50 0 50 10000.511.522.533.5 x 105?100 ?50 0 50 10000.511.522.533.54 x 105Flying MPB Nesting MPBSpace (m)Flying MPB Nesting MPB High PopulationLow PopulationPopulation Size (MPB/ha)t=0 fht=1.6 fht=3.2 fht=4.8 fht=6.4 fht=8 fhFigure 3.8: Spatial population densities of flying and nesting MPB in alandscape that has a clearcut harvested region and uniform intial condi-tions. The spatial plots are shown over increasing time simulation in MPBflight hours (fh). The two columns on the left depict the MPB at a lowinitial population density (J = 31, 250 MPB/ha, approximately 3 sourcetrees per ha) and the two right columns illustrate high initial populationdensities (J = 218, 700 MPB/ha, approximately 22 source trees per ha).The simulations shown assume diffusion is 1 ha/fh.583.6. Model Results?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1008100820083008400850086008700?100 ?50 0 50 1001.821.841.861.881.91.921.941.96 x 104?100 ?50 0 50 1002.322.342.362.382.42.422.442.462.482.52.52 x 104?100 ?50 0 50 1002.552.62.652.72.752.82.85 x 104?100 ?50 0 50 1002.62.652.72.752.82.852.92.95 x 104?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1000.20.250.30.350.40.450.5?100 ?50 0 50 100681012141618?100 ?50 0 50 1002530354045505560657075?100 ?50 0 50 1006080100120140160180?100 ?50 0 50 100100120140160180200220240260280?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1005.655.75.755.85.855.95.9566.056.16.15 x 104?100 ?50 0 50 10000.511.522.533.5 x 106?100 ?50 0 50 10000.511.522.53 x 106?100 ?50 0 50 10000.511.522.53 x 106?100 ?50 0 50 10000.511.522.5 x 106?100 ?50 0 50 100?1?0.8?0.6?0.4?0.200.20.40.60.81?100 ?50 0 50 1005060708090100110120130140150?100 ?50 0 50 1000246810121416 x 104?100 ?50 0 50 1000246810121416 x 104?100 ?50 0 50 1000246810121416 x 104?100 ?50 0 50 1000246810121416 x 104Flying MPB Nesting MPBSpace (m)Flying MPB Nesting MPB High PopulationLow PopulationPopulation Size (MPB/ha)t=0 fht=1.6 fht=3.2 fht=4.8 fht=6.4 fht=8 fhFigure 3.9: Spatial population densities of flying and nesting MPB in a land-scape that has a thinned harvested region and uniform intial conditions. Thespatial plots are shown over increasing time simulation in MPB flight hours(fh). The two columns on the left depict the MPB at a low initial populationdensity (J = 31, 250 MPB/ha, approximately 3 source trees per ha) and thetwo right columns illustrate high initial population densities (J = 218, 700MPB/ha, approximately 22 source trees per ha). The simulations shownassume diffusion is 1 ha/fh.593.6. Model Results?100 ?50 0 50 10000.0020.0040.0060.0080.010.0120.014?100 ?50 0 50 100?2024681012 x 10?5ClearcutPheromoneChemotaxisSpace (m)Figure 3.10: This graph represents the MPB pheromone profile (left) andresulting chemotactic response (right) near the beginning of the simulationbefore significant nesting occurs. It can be seen that the pheromone profileacross the border between the unharvested and harvested edges has a steepgradient which causes significant attraction (large chemotactic response)to the unharvested region. This specific simulation is at a high density(J = 218, 700 MPB/ha, approximately 22 source trees per ha) for the uni-form initial conditions with a clearcut harvested region. The simulationsshown assume diffusion is 1 ha/fh. Pheromone concentration is in ?g/haand chemotactic response has units 1/fh.603.6. Model Results20 m, has more MPB than the trap at -90 m for larger initial densities ofMPB. The initial MPB aggregation occurs at -20 m and the excess MPBnear this trap may spillover to the 20 m trap. In Figure 3.9 (at 4.8 fh) itis evident that the right peak of flying MPB density (the one progressingthrough the harvested region) is larger than the peak in the unharvestedregion. This is due to the greater attraction of MPB in the harvested regionbecause of a sharper gradient in the MPB pheromone profile. Additionally,the pheromone profile at -90 m is largely smooth. The main reason that theMPB pheromone concentration is lower in the harvested region is the lowerdensity of susceptible trees. This results in a lower density of nesting MPBand thus a lower concentration of pheromone.We tested the effect of faster movement by increasing the diffusion by afactor of 2. This is represented in the plots labeled fast diffusion in Figures3.6 and 3.11. The increase in diffusion will have two main effects. First, ifMPB are initially only in the dense forested region (Figure 3.6), the fasterdiffusion allows more MPB to travel into the harvested region and MPB tomove through the harvested region more quickly. This changes the rate atwhich the MPB approach the pheromone traps. The second effect of thefaster diffusion is that it smooths out the sharp peaks of flying MPB. MPBwill still be attracted to pheromone, but the peaks that form will not beas sharp or large. This is evident in either initial condition, as shown inFigures 3.6 and 3.11.In Figure 3.6, the faster diffusion increased the trap catches in the ? 20m and 90 m traps, and decreased the trap catches at the -90 m trap. Thismakes sense as the higher diffusion rate allowed a larger density of MPB tomove into the harvested region. The trap catch decreasing at -90 m is due toa two-fold effect. First, the diffusion smooths the peak of MPB that formsat this site and second, the faster diffusion of MPB into the harvested regiondecreases the number of MPB present in the unharvested region. Finally,the higher number of MPB in the trap at -20 m is due to the smoothing anddiffusion of MPB away from the trap at -90 m.The effect of faster diffusion with uniform emergence of MPB (Figure3.11), is mainly a decrease in the trap at -20 m. This also changes the613.6. Model Results0 0.5 1 1.5 2 2.5x 1050501001502002500 0.5 1 1.5 2 2.5x 1050204060801001201401600 0.5 1 1.5 2 2.5x 1050204060801001201401600 0.5 1 1.5 2 2.5x 105050100150200250Average Trap CatchesInitial Population Size (MPB/ha)Regular DiffusionFaster DiffusionClearcut Thinned?90 m?20 m20 m90 mFigure 3.11: Average trap catches in the bait sites at ?20,?90 m with aninitial condition of MPB at a uniform density over the domain. The bait trapcatches are shown for heterogeneous landscapes with the harvested regioneither clearcut or thinned. Additionally, the diffusion rate is changed to seethe effect of increased movement rate in harvested regions. Note that thepeaks in trap catches differ between the regular and faster diffusion cases.623.7. Connections Between Theoretical and Empirical Resultstrap catches at 20 m. The smoother pheromone peak at -20 m results ina smaller attraction and consequently a lower density of MPB. In the casewith a thinned harvested region, this means there are fewer MPB which willfind the nearby trap at 20 m. If the harvested region is clearcut, the lowerattraction of pheromone at -20 m (when compared to the regular diffusioncase) results in less MPB moving from the 20 m trap to the -20 m trap.Thus, there is a slight increase in the MPB density at 20 m if the harvestedregion is clearcut. The final consequence of a lower number of MPB at -20m is an increase in the trap catches of MPB at -90 m.3.7 Connections Between Theoretical andEmpirical ResultsWe did not have enough information in the field study to discernthe initial population sizes and whether the initial distribution of MPB wasuniform or not. We therefore discuss the comparison between our simulationresults and data first under the assumption of uniform initial conditions andsecond under the assumption of dense forest initial conditions.Using uniform initial conditions, we find there is a more significant dif-ference between the clearcut and dense region than there is between thethinned and dense region (Figure /reffig:uniformic). This is consistent withour experimental results which show a statistically significant difference be-tween the clearcut and dense region but not between the thinned and denseregion. In general, the difference between the harvested and unharvested re-gions is smaller at lower population sizes but increases as MPB populationsize increases. An exception to this trend is the trap catches at 20 m.Empirical results found no statistically significant difference between thetwo types of harvested regions, thinned and clearcut. This is consistentwith our simulation results for the 90 m trap and for the 20 m trap for lowpopulation densities.The final empirical result is that pairs of traps in the harvested andunharvested region were not significantly different. Our simulation results633.7. Connections Between Theoretical and Empirical Resultsin the harvested region agree with the empirical results. That is, traps inthe harvested region were not significantly different for low population sizes.In contrast to the empirical results, however, only at the lowest populationsize are the simulation predictions of the bait trap catches the same in theunharvested region.Under the assumption that MPB are initially only present in the denseforested region, we find that there are only small differences between alltrap catches. Only at very high densities do the simulation results predictany difference between the harvested and unharvested traps. Thus, thesesimulation results agree with the experimental results that the differencebetween traps is not significant at different locations within the harvestedor unharvested region. Additionally, the simulation results agree that nosignificant difference will be found between the traps between the clearcutand thinned region, as well as no significant difference across the traps in thethinned and dense regions. In contrast to the experimental results, however,our simulation results predicts no significant difference between the clearcutand dense region.Under either initial condition, a faster diffusion of MPB resulted in a de-creased difference between all traps. The simulations shown were done withthe diffusion increased by a factor of 2. Other simulations were performedwith the diffusion increased by a factor of 5 and 10 (results not shown). Asexpected, trap catches became even closer as diffusion increased. Currently,our model simulation results predict that the pairs of traps in the harvestedor unharvested region will only be similar for low population sizes, underuniform initial conditions. If we assume diffusion is increased, our simula-tion results become closer to the experimental results, since the differencesbetween pairs of traps in each region decrease.To test the effect of the pheromone response function chosen [WP97,HFSL06] we ran the simulations with the effects of MPB pheromone ne-glected (results not shown). That is, the pheromone from the bait traps stillworked as usual, but the pheromone produced by the MPB had no attractiveor repulsive effect on the MPB. This led to the prediction of increased MPBnumbers in the harvested region, and smaller differences between pairs of643.8. Discussiontraps in each region when compared to the standard pheromone responseused in this study.If population levels increased beyond the results shown in Section 3.6(results not shown), it was found that the repulsive nature of the pheromoneled to increased densities in the harvested zone. Thus the repulsive nature ofthe pheromone becomes increasingly important at higher densities of MPB.3.8 DiscussionOur experimental results indicate that differences in MPB trap catcheswere not statistically significant between the dense forested region and thinnedforested region as well as the pairs of traps in either the harvested (thinnedor clearcut) or unharvested (dense forest) landscape. Finally, the differencein MPB trap catches between the different types of harvested region wasnot statistically significant.Our simulation results most closely agreed to experimental results whenthe initial distribution of MPB on the landscape was uniform. This resultsupports the hypothesis that MPB were fairly uniformly spread over theharvested and unharvested landscape at the start of the summer. This dis-tribution could be due to the free-flight behaviour of the MPB [Hug02], orcould be due to external sources of MPB blanketing the landscape [dCA12].Furthermore, for the clearcut and the dense forested region to show statis-tically significant differences between trap catches, the population of MPBmust be sufficiently large. That is, the experimental results show a statisti-cally significant difference between the clearcut and dense region, and thus,according to our simulations, the population of MPB must be sufficientlylarge before this difference becomes apparent.In contrast, the simulation results where MPB emerged only from theunharvested region did not agree well with the experimental results. Thedense forest initial conditions lead to smaller trap catches in the harvestedregion then was observed in the experimental trap catches. Additionally, weexpected initial aggregation at the trap sites, which did not occur due tothe diffusive spread of the MPB population in the unharvested region into653.8. Discussionthe harvested one.The experimental results illustrated very little difference between trapsat 20 and 90 m from the edge of the habitat in both the harvested and un-harvested regions. Our simulation results become closer to this experimentalresult as we assume that diffusion increases over the landscape. Faster dif-fusion means smaller catches at each trap, but the trap catches at 20 and 90m from the edge also become more similar. This result supports the hypoth-esis that MPB may move faster in harvested regions due, possibly, to lesstortuous paths [Rei09]. Higher diffusion rates in matrix habitat has beenobserved in studies of other insects such as leaf beetles and planthoppers[CDO07, HC03, HC06].As the density of attacking MPB increases, pheromone interactions be-come more important. In the current chapter we assume that pheromonecan be both attractive and repulsive [WP97], and this combination allows forthe spreading of MPB attack from an initial aggregation. Strictly attractivepheromone responses were suggested by Hughes [Hug02], but in our modelthis would lead to uncharacteristic densities of MPB in a single aggregationat the initial site.In our study, we also made the assumption that MPB had a very syn-chronous emergence. If this is not the case, the density of MPB needs tobe significantly higher to nucleate a successful initial aggregation. In thisscenario, we would expect much lower trap catches than were reported inthe simulation results in Section 3.6 due to lower densities of nesting MPB,and consequently lower densities of flying MPB attracted to the sites withnesting MPB.In the simulation results with uniform initial conditions, we found thatthe pheromone produced by traps was very important in determing theinitial aggregation position of MPB. After this point, however, we found thatthe pheromone produced by MPB had a much larger impact. That is, athigher densities of MPB the pheromone produced by MPB was much moreeffective in determining the behaviour of MPB than the trap pheromone.This supports the conclusion that natural sources of pheromone can largelyaffect the timing and trap catches of MPB, as predicted by Bentz [Ben06].663.8. DiscussionOur model does not include any explicit boundary behaviour at theinterfaces between regions, but several empirical [CC82, Jon77, RH98, SC01,WSM04] and theoretical studies [ML, GTLS] have shown that this may be animportant factor. Another direction for future study is to consider spatiallyexplicit diffusion, changing in each patch. Variations in spread rate havebeen found to be important in other theoretical studies [ML].The difficulties in determining the movement behaviour of MPB is well-known, especially under habitat fragmentation [BAL93] and when using trapcatches to provide information about MPB densities [Ben06]. Our studylinks mathematical and experimental studies to create a powerful approachto elucidate MPB spatial dynamics. Through our combined approach weprovide evidence that MPB move faster in harvested landscapes, and thatMPB emergence and movement behaviour leads to a more uniform densityof MPB prior to attacks on susceptible trees.The results of our study gives new insight into several areas. It has beenhypothesized that MPB movement may be faster in harvested landscapes[Rei09, BPBL00]; our study provides significant support to this hypothesis.Additionally, understanding the dispersal of MPB in fragmented landscapesis crucial to inform management decisions, as MPB moving differently inmatrix habitat will change the optimal strategy for control of the insect. Inaddition, the results of our simulation helps us understand the movementbehaviour of MPB prior to attack. The fairly uniform initial distribution ofMPB is very important, as it allows us to understand the spatial distributionof MPB prior to attack, and will allow further study to target the best wayto control MPB spread.67Chapter 4Impacts of Management onMountain Pine BeetleMovement4.1 IntroductionThe economic and widespread impact of the MPB has led to extensivemanagement efforts, through direct and indirect control [SW06]. Direct con-trol mechanisms have focussed on harvesting newly attacked (green attack)trees to prevent the spread of the MPB before they can emerge from thetree. Additionally, semiochemicals, both attractive and repulsive to the bee-tle, have been used to help prevent spread and to direct MPB movementtowards specific stands. This allows for a more targeted management of thebeetle. A direct control mechanism which has been tested in the past isprescribed burning of beetle infested areas of the landscape. The problemwith this control technique is that the beetle can survive unless the fire isa high-intensity one. Since high-intensity fires are dangerous and very diffi-cult to control, prescribed burning is no longer used as a method of directcontrol. The most common indirect control method has been prescribedburning and harvesting of susceptible lodgepole pine tree stands which havenot yet been attacked by the MPB. This control method is used to decreasethe density of susceptible host trees.In the past two decades, MPB populations rose to epidemic levels andhave had devastating effects on forests in both British Columbia [SW06] andAlberta [TR08]. In particular, Banff National Park (NP) has been affected684.1. Introductionby the MPB since 1997 [TR08]. Management efforts in Banff have focussedon prescribed burning and baiting and felling of green attack trees. Baitingof the trees focusses the attack of nearby MPB, allowing Parks Canadapersonnel to do intensive land surveys around the baited area to fell andburn trees recently attacked by the MPB. Prior to the initiation of MPBcontrol efforts, Banff NP was divided into two areas, one with managementand one without. The resulting data can help elucidate the effectivenessof control efforts. An initial study found that the management activitylimited the long-distance spread of the MPB, but did not reduce the areaof lodgepole pine affected [TR08]. However, this conclusion was confoundedby differences between the two zones: the management zone had a largerarea of dense, susceptible lodgepole pine than the control area. There isneed for a mathematical model to test competing hypotheses to explain theempirical results. This model could also be used to test the effectivenessof management activities in silico, without long-term consequences to thelandscape.The effect of prescribed burning on MPB depends critically on the in-tensity of the fire. High-intensity fires can reduce brood density and re-production [SLSH01], reduce the attractiveness to MPB [SLSH01, Mur77,SG89, ER04], reduce susceptibility of the forest to MPB attack, and reducetraversability of the landscape to MPB [BLB+06]. As mentioned earlier,however, high-intensity fires can be difficult to control and can result inmajor damage to the site (such as soil damage) [SG89]. Suprisingly, lowerintensity fires that only moderately burn trees, may actually increase theattractiveness of these sites to MPB [SLSH01, SG89]. In addition, theremay be increased competition with other insects such as Pine Engravers(personal communication, Jane Park), and possibly increased [MKW03] ordecreased [Rei07] resin production.The pattern of attack of the MPB is intrinsically determined at epidemicand incipient epidemic densities [WP97] (Chapter 2) and is determined bythe susceptible landscape at endemic densities. Similarily, at low densitiesof MPB, stressed, partially burned trees on the edge of a prescribed burn arepreferentially attacked [ER04]. At high densities of MPB, the same pattern694.1. Introductionis not observed, and burned and unburned trees are equally likely to beattacked.Previous modelling efforts have investigated beetle dynamics over a sin-gle season, but few models have looked at beetle populations across mul-tiple years. Of those that have looked at multiple-year dynamics, thereare no models of which we are aware that explicitly model both populationdynamics and dispersal, including the interactions between beetles, beetlepheromones, tree kairomones, and the forest landscape. These interactionsare extremely important as they will determine the spread dynamics of MPBover the landscape.In this chapter, we build a multi-year spatially explicit model (Section4.2) to test the effectiveness of the management strategies of baiting andtree-removal (Section 4.3.1) and prescribed burning (Section 4.3.2). Usingour model, we find that baiting and tree-removal is always more successfulat reducing MPB density and forest impact at high initial MPB populationdensities but may lead to greater forest impact and greater MPB populationgrowth at low initial levels of MPB. We predict that tree removal withoutbaiting can be more successful than combined baiting and tree removal ifthe searched area is optimal (high MPB density). Finally, our analysisindicates that prescribed burning is more effective than unmanaged controls,and can be more effective than clearcutting given assumptions about thereproductive output and attractiveness of burned trees.In Section 4.4 we present indicators of MPB population densities onthe landscape that emerge from simulations of our model in the absenceof control strategies. First, we discuss the number of source trees in acluster needed for MPB to initiate an epidemic (and conversely the thresholdrequired to reduce to endemic population densities). Second, we discussthe distance needed between source trees for MPB to successfully mass-attack the susceptible landscape. Finally, we discuss the width of treesremoved needed for a clearcut to be successful in reducing MPB spread andreproduction.704.2. Mathematical Model4.2 Mathematical ModelOur model consists of a discrete set of difference equations for theoverwintering period (September-June) and a continuous set of reaction-diffusion-chemotaxis equations for the beetle flight, emergence, and attackperiod (July-August). Existing theoretical and empirical work [RB83, PMW98]informed the selection of parameter values, and the choice of structure foreach term in the equations.We assume that the continuous set of Partial Differential Equations(PDEs) describing the period of flight, emergence, and attack incorporatesfive variables. These variables are P - the density of flying MPB, Q- thedensity of nesting MPB, J - the density of pre-adult MPB (beetles notyet emerged from the bole of the tree), A - the concentration of beetlepheromone, and C - the concentration of tree kairomones. The set of dis-crete equations describing the overwintering period incorporates only twovariables, St - the density of susceptible trees and Lt - the density of over-wintering pre-adult beetles. The equations for the summer dispersal are?P?t=diffusion? ?? ??p?2P ?chemotaxis? ?? ??[(?c?C + ?ab0 ?Ab0 +A/b1?A)P]?death?????pP ?nesting? ?? ?g(P,C)St +emergence? ?? ??(t)J ,(4.1a)?Q?t=death? ?? ???qQ+new nesters? ?? ?g(P,C)St , (4.1b)?C?t=diffusion? ?? ??c?2C +synthesis????a2St ?degradation?????cC +baiting????? , (4.1c)?A?t=diffusion? ?? ??a?2A +synthesis????a1Q ?degradation?????aA , (4.1d)?J?t=emergence? ?? ???(t)J , (4.1e)714.2. Mathematical Modeland the equations for the overwintering stage areLt+1 = r1Q(e?r2(r3Q?1r4)r5? e?r6(r3Q?1r4)r7), (4.2a)St+1 = St(1? ?L2t+1L2t+1 + k2j), (4.2b)where,g(P,C) =landing?????Pmass attack? ?? ?P 2P 2 + (kp)2.The functional form of (4.2a) incorporates both competition of larvae forresources and tree defences against MPB. At low densities of nesting MPB,the MPB will not survive as tree defences will overcome the small nestingbrood. At high densities of MPB, competition becomes a factor and survivalof larvae is decreased due to scarce resources. The maximum reproductionand survival of larvae occurs at medium densities of MPB, large enough toovercome tree defences, but not too large as to compete within the host tree.The model we choose, represented by the system of equations in (4.1) and(4.2), is fairly complex. We feel this complexity is necessary to incorporatethe interactions between MPB, their pheromones, and lodgepole pine treesfor a multi-year model. A simpler model may be able to understand partsof this system, but in order to understand the system as a whole, this levelof complexity is required. We assume that the initial conditions at thebeginning of each summer are P (~x, 0) = Q(~x, 0) = C(~x, 0) = A(~x, 0) = 0and J(~x, 0) = Lt(~x). That is, the density of pre-adult beetles in the summeris determined entirely by the density of pre-adult beetles produced from thelarvae that survive the previous winter. The model parameters are discussedbelow.724.2. Mathematical Model4.2.1 Estimation of parametersMany of the parameters for the continuous portion of the model (4.1)were estimated by Powell et al [PLB96, BPBL00] by projecting a spatialpartial diffential equation model onto the local (single-tree) scale and thenutilizing historical MPB data. Most of the parameters for the discrete por-tion of the model (4.1) are estimated from a study by Raffa and Berryman[RB83].The remaining unknown parameters are ?, kp, kj , r3, and r4. We chosethe rate of emergence parameter, ?, such that 99.9% of the emerging MPBleave the host tree after approximately 20 beetle flight hours. This time cor-responds to approximately half the beetle flight period. All other parameterswere chosen based on the data provided by Raffa and Berryman [RB83]. Theparameter kp is the density of flying MPB required for 50% nesting beetlesuccess. Assuming each female beetle makes a single gallery, we estimatethe density of flying MPB required for a 50 percent success of mass attackto be 40 beetles / m2 [RB83]. Since our model uses ha instead of m2 wemust multiply this quantity by a conversion factor:Bh = BmSATA, (4.3)where Bh is the number of beetles per hectare, Bm is the number of beetlesper m2, SA is the surface area of the tree attacked, and TA is the areawithin which a beetle is considered to be ?attacking? a tree. In the discreteequations (4.2), r3 is defined as the area within the tree that the beetle isattacking and r4 is defined as the surface area of the tree attacked. Thus,r3 = TA and r4 = SA.We estimate the surface area as SA = pidh, where the tree diameter, d,can range from 0.1874 m to 0.3456 m and the height h is taken to be 7.5m [RB83]. Therefore, SA can range from 4.42-8.14 m2. We assume thateach tree is centred in a 10 m2 attractive zone. Any beetle flying in thatzone will be attracted to and join the attack on the central tree. Thus,TA = 10m2 = 10?3 ha, which makes kp = 176, 800 ? 325, 600 beetles/ha.We initially assume the intermediate value of kp=250,000.734.2. Mathematical ModelThe parameter kj is the density of surviving MPB, Lt, required for 50%tree mortality to occur. We approximate that the density of pre-adult beetlesrequired for this to occur is 5 beetles/tree. Therefore, using the same esti-mate of TA = 10m2, we find kj is approximately 5000 pre-adult beetles/ha.We choose this small value for kj since a non-zero density of surviving MPB,Lt, is a signal for tree death. In the model, tree resistance to MPB occursin the summer nesting term (4.1) and in the discrete equation for Lt (4.2).The dimensional model (4.1) and (4.2) parameter values can be found inTable 4.1 while the non-dimensional model (4.4) can be found in Table 4.2.4.2.2 Non-dimensionalizationWe non-dimensionalize the model for ease of mathematical analysis.We choose the following non-dimensionalizations for our variables:Q =a1b0?aQ, P =a1b0?aP, C =?c?cC,St =??pSt, A =1b0A, J =a1b0?aJ,t = ?at, (x, y) =??a?a(x, y).By examination of the parameter values in Table 4.1 it can be seen thatmany parameters differ by orders of magnitude. To identify parameterswhich are small (act over slow time scales) and large (act over fast timescales) we choose a scaling parameter, , =?p?a.Using this scaling parameter, we define the dimensionless parameters?p =?p?a, ?a =?ab0?a, kp =kpa1b0?a,?c =?c?a, ? =a2?c??c, ? =?c?a?c?,744.2. Mathematical ModelTable 4.1: Table of parameter (Par) values for the dimensional model (4.1)and (4.2).Par Description Units Value?p diffusion of flying MPB hafh 1?c diffusion of host volatiles (kairomones) hafh 0.648?a diffusion of beetle pheromones hafh 0.648?c kairomone attractiveness ha2?g?fh 0.8?a beetle pheromone attractiveness ha2?g?fh 5.7b0 pheromone dissipation concentration?gha 5.4b1 concentration at which pheromone is satu-ratedn/a 1? random landing rate of flying MPB hatrees?fh 0.16kp flying beetle density required for 50 percentsuccess of mass attackmpbha 250000?p death rate of flying MPB fh?1 0.014?q death rate of nesting MPB fh?1 0.001?c degradation of tree kairomones fh?1 180?a degradation of beetle pheromone fh?1 180?(t) emergence rate of pre-adult beetles fh?1 0.345a1 rate of pheromone increase due to nestingMPB?gfh?mpb 0.02a2 rate of kairomone increase due to susceptibletrees?gfh?trees 0.02r1 production of pre-adult beetles (pupae) fromnesting MPBn/a 26.3r2 density dependent effects of nesting MPB(m2mpb)r50.0015r3 area surrounding tree where MPB are at-tacking treeha 0.001r4 surface area of tree attacked m2 6r5 constant of MPB competition n/a 1.5r6 success of tree defenses(m2mpb)r70.0026r7 cooperative effects of MPB n/a 1.6? maximum success of pre-adult beetles n/a 1kj pre-adult beetle density required for 50 per-cent success in overtaking treempbha 5000? rate of kairomonne production by traps ?gha?fh 1458754.2. Mathematical ModelTable 4.2: Table of parameter values for the non-dimensional model (4.4).Parameter Value 0.0000778?p 1.54?c 1?a 47.5? 0.154kp 5.14?q 0.0714?c 1? 24.6r3 48.6kj 0.103? 10?c =?c?a, ? =??a, ?q =?q?a.Next, my define scaling factors for the terms in my set of discrete equations,St+1 =?p?St+1, Lt+1 =b0?aa1Lt+1,r3 =b0?aa1r3, kj =a1b0?akj .Substituting these values into (3.2) we arrive at the non-dimensionalmodel (4.4). The number of parameters was reduced from 25 to 19 duringthe non-dimensionalization process.Summer Dispersal:?P?t=?p?2P ? ?c?(P?C)? ?a?(1?A1 +A/b1P?A)? P ?  P3P2+ kp2St + ?J,(4.4a)764.3. Effects of Management?Q?t=? ?qQ+ P3P2+ kp2St, (4.4b)?C?t=?c?2C + ?St ? ?cC + ?, (4.4c)?A?t=?2A+Q?A, (4.4d)?J?t=? ?J. (4.4e)Overwintering:Lt+1 =r1Q(e?r2(r3Q?1r4)r5? e?r6(r3Q?1r4)r7), (4.4f)St+1 =St(1? ?Lt+12Lt+12+ kj2 ). (4.4g)Our mathematical model was simulated using Fortran with the numericalscheme as described in Chapter 3.4.3 Effects of ManagementOur study is motivated by questions regarding the impacts of MPBmanagement activities, with specific focus on the management activitiescarried out in Banff National Park in Alberta, Canada. We carried outseveral different studies to investigate the effects of baiting and tree removal(Section 4.3.1), and clearcutting and prescribed burning (Section 4.3.2). Ourstudy of baiting and tree removal focuses on simulations with three (Section4.3.1) and five source trees (Section 4.3.1). Other source tree densities weretested, though the results displayed are sufficient to explain the trends forbaiting and tree removal.4.3.1 Baiting and Tree RemovalAn initial study of the effects of pheromone baiting and tree removalfound unexpected differences in MPB attack between the management and774.3. Effects of Managementmonitoring zones in Banff. It was found that the impact (area of trees con-sumed) of the beetles in the two zones was the same, though managementdid slow the spatial spread of the beetle [TR08]. Figure 4.1 shows an il-lustration of both patterns of attack. In both cases, beetles emerged froma source at the centre of the grid, and the number of grid cells attackedby MPB is the same. The main difference between the two zones is thedistance of beetle spread (distance from the attacked cells to the centre ofthe grid). A confounding factor in the study was that the management areacontained susceptible habitat of higher quality (higher MPB reproductivesuccess) than the monitoring zone. The results also led us to ask if theremight be a situation in which it would be better to only do removal of green-attack trees without the baiting component. Indeed, while Parks Canadaused the baiting and tree removal (in Banff), elsewhere in the provincethe Alberta Environmental and Sustainable Resource Development (ESRD)used solely tree removal with no baiting component. The motivation forusing pheromone baiting is that it concentrates the MPB attack to a knownlocal area, so that the management knows where to search, whereas withoutpheromone baiting the attacked trees are harder to find.Three source tree simulationsTo elucidate the effects of tree baiting and removal we set up a 1-dimensional simulation with three source trees that were attacked in theprevious year. Neighbouring source trees are placed a constant distance ds,apart from the central source tree. This landscape is depicted in Figure 4.2.The simulation included a single bait location, which produced pheromoneduring the first year of study. At the end of the summer dispersal season,once MPB had settled inside new hosts, a search radius, rsearch, around thebait trap was searched. A proportion, p, of MPB and susceptible trees wereremoved from this search radius based on the search success. The removal oftrees and MPB occured only if the MPB density reached a threshold. Belowthis threshold, it was assumed that the tree could successfully defend thelow density of MPB and thus the tree was not attacked at sufficiently high784.3. Effects of ManagementRmonRmanFigure 4.1: Simple grid representations of the monitoring (top) and man-agement (bottom) zones. Areas consumed by the beetles are shown in red,while areas of unattacked forest are shown in green. R is the radius fromthe central emergence region to the farthest attacked stand; with the sub-scripts ?man? referring to the management zone and ?mon? referring to themonitoring zone.794.3. Effects of Managementd dBait Sites sSearch Diameter (100 m)Search Success (80%)Figure 4.2: Simulation initial conditions are shown in the top figure. Thereare three source trees, separated by distance, ds, with a bait site to the rightof the center source tree. Once the summer dispersal season is over, treesthat are attacked are shown in red. Management searches a radius aroundthe bait site (in this example 50 m) and removes attacked trees and MPBwith some success rate (80% in picture).densities and consequently, was not removed.The simulation was run for two years, starting with emergence and dis-persal in the first summer, and with baiting and/or tree removal only oc-curing in the first summer and fall. We need only to simulate two years, asthis is sufficient to illustrate the effects of management. Simulations over 3or more years did not yield novel results which the 2 year process could notexplain. We investigated three separate management scenarios: no man-agement (monitoring/control), baiting only, and baiting with tree removal.Running the control case was important for determining the effects of themanagement activity, and baiting only was investigated to determine theeffects of removal after baiting. These three types of management led to theresults shown in Figures 4.3 and 4.4.The density of nesting MPB at the end of each summer dispersal periodin years 1 and 2 is shown in Figure 4.3. Within each plot, results are shownfor three inter-source tree distances, ds. The control case (no management)804.3. Effects of Management?200 ?100 0 100 20000.511.52 x 105?200 ?100 0 100 20000.511.52 x 105?200 ?100 0 100 20000.511.52 x 105?400 ?200 0 200 40000.511.522.53 x 105?400 ?200 0 200 40000.511.522.53 x 105?400 ?200 0 200 40000.511.522.53 x 105100 m200 m50 mMPB Density (Beetles/ha)Baiting OnlyMonitoring ZoneYear 1Year 2Baiting and Tree RemovalManagement ZoneSpace (m)Figure 4.3: Plots of nesting MPB density over space at the end of each yearfor two years. Simulations are shown in the monitoring zone (no manage-ment, left), and with the management activities of baiting only (middle),and combined baiting and tree removal (right). Three distances, d, betweensource trees are shown; 50, 100, and 200 m. The top panel shows nestingMPB densities at the end of the first summer, while the bottom panel showsnesting MPB densities at the end of the second summer.814.3. Effects of Management50 100 150 200010020030040050060070080050 100 150 20002000400060008000100001200014000No ManagementBaiting OnlyBaiting and Tree RemovalLength of Forest Consumed (m)Distance Between Source Trees (m)Total Beetle PopulationFigure 4.4: Length of forest consumed (left) and MPB population size (right)are shown for different distances between source trees. The three curves rep-resent the impact of MPB with no managment (black), baiting only (blue),and combined baiting with tree removal (red). Simulations were run for twoyears to obtain these densities and length of forest impact.824.3. Effects of Managementis shown in the left panel. As ds is increased, it is harder for MPB toaggregate in sufficient numbers to launch a mass-attack on the trees. Notethat spreading source trees farther apart also simulates lower landscape-levelpopulation densities of MPB. In the case where source trees are too far apart,the MPB will successfully nest in the first year, but will not produce enoughoffspring for successful mass attack and nesting in the second year. Thus,there is an MPB population density below (or ds above) which the beetlepopulation dies out. This situation does not occur however, for baiting only(middle panel, Figure 4.3) for the ds simulated. There is still some minimumpopulation threshold for the baiting only case, but it is not reached for theds simulated. Baiting focusses MPB aggregation on the bait site, increasingthe success of MPB at the initial aggregation stage. Consequently, thespread and impact of a given MPB population is increased in the presenceof baiting. The effect of tree removal combined with baiting are shown inFigure 4.3 (right panel). The combined approach is effective at reducing thepopulation of MPB below the levels seen with monitoring only.A summary of the impact and extent of the MPB infestation after twoyears is shown in Figure 4.4. These plots show the total MPB populationand total susceptible landscape affected by the end of the second year. Notsuprisingly, baiting only is the worst management option, as it improves thesuccess of initial MPB aggregation, and so the length of forest consumedand population size of MPB at the end of the second year is largest. Baitingin combination with tree removal yields the lowest population of nestingMPB remaining. This result gives support to baiting with removal as amanagement tool for reducing the density of MPB. In contrast, the initialsuccess of the MPB aggregation with baiting allows the MPB populationto successfully attack a larger portion of the landscape than it would haveattacked in the control case. The impact of MPB on the landscape is thushigher with baiting and tree removal when the initial MPB source trees arepositioned far apart (akin to a low MPB population density).These simulations were run assuming a constant search radius of 50 mand a search success of 80%. We then repeated the numerical experimentwith a search radius ranging from 25-150 m, and search success ranging834.3. Effects of Managementfrom 0.2-1.0. A search success of 1.0 represents the case where all green-attack trees in the search area are successfully found and removed. Ourresults are displayed in Figure 4.5. As search success or search radius in-creases, the density of nesting MPB after the second year decreases. Theseare both intuitive results. In addition we can use our simulations to de-termine whether it is better to search in a small area with high success orin a large area with low success. According to our results, it appears thatincreasing search radius is more important than increasing search success.In our graphs, when we increased the search radius to 100 m or more, onlythe lowest search success (0.2) had MPB after the second year. We sug-gest that increasing search radius is more important than increasing searchsuccess since increasing search radius will dilute the overall density of MPBemerging the following year. If the MPB population is sufficiently reducedand spread out over an area, it is more difficult for the MPB population tosuccessfully aggregate and mass attack susceptible trees.Five source tree simulationsTo examine the effect of adding more source trees to the landscape,we increased the number of source trees in our simulation from three tofive. For our five tree simulations, we simulated four management activities.Baiting with tree removal, baiting only, tree removal only, and control (nomanagement). Tree removal uses the same conditions (search area, searchsuccess, and removal threshold) as combined baiting and tree removal, buthas no baiting component during the flight period. We had 5 initial sourcetrees equally spaced over the susceptible landscape, and increased the inter-source tree distance, ds. As we increase ds we find that the length of forestconsumed and the total population of nesting MPB in the second year de-creases, which is consistent with our three-tree simulations. We display the5-tree results in Figure 4.6.In these simulations we additionally vary the initial population densityin the source trees to test the effect of initial population size. We find thatthe results are very different as we change the source tree MPB densities.844.3. Effects of Management0.2 0.4 0.6 0.8 1020004000600080001000012000Search EffectivenessTotal Beetle Population 25 m50 m100 m150 mFigure 4.5: MPB population sizes with varying search width and searcheffectiveness. Search widths shown in this simulation are 25, 50, 100, and150 m. Search effectiveness values were 0.2, 0.4, 0.6, 0.8, and 1.0. A searcheffectiveness of 1.0 designates that 100% of the attacked trees in the neigh-bourhood of the bait tree are found and removed854.3. Effects of ManagementAt high densities, the relative impact of each management activity does notchange with variations in ds tested. For both length of forest consumed andtotal MPB population after the second year, the best management activityis tree removal only, followed by baiting and tree removal, then control, andfinally baiting only. For all inter-source tree distances, ds, tested, there werealways MPB remaining after the second year. If we assume however, thatthe source trees had an even smaller initial population density, we find thatpopulations will not survive to the end of the second year for some distancesbetween source trees, and the trends in MPB population size and the lengthof forest consumed are different.From the point of view of MPB population size after the second year, themanagement strategies rank (in ascending population size) as tree removalonly, baiting and tree removal, monitoring only, then baiting only for lowdistances between source trees (Figure 4.6, top panel). As the spatial differ-ence between source trees increases, there is a switch and more MPB mayresult after utilizing baiting and tree removal approach than with monitoringonly. The same trend that is observed in length of forest consumed.The result that monitoring alone will produce a lower number of MPBand a smaller length of forest consumed than baiting and tree removal atlow initial MPB densities (source trees spaced far apart) is due to the en-hanced initial aggregation and increased success of MPB due to the baitingcomponent of management. These results are consistent with the three treesimulations: Baiting and tree removal is a successful management strategyunless MPB population density is low.If the same simulations are run with management efforts occuring inboth years (Figure 4.7), and compared with the single year of managementscenarios (Figure 4.6) the results are intuitively expected. Tree removalonly over both years reduces the impact of MPB even further, regardless ofspacing between source trees, whereas the effects of a second year of baitingand tree removal are dependent on the source tree spacing. If the sourcetrees are close together (MPB density is high), baiting and tree removal overconsecutive years reduces the impact of MPB (length of forest consumedand MPB population survived) below that experienced with a single year of864.3. Effects of Management50 100 150 20005000100001500050 100 150 20011.522.5 x 10450 100 150 200010020030040050060070080050 100 150 2005006007008009001000110012001300Total Beetle PopulationDistance Between Source Trees (m)Length of Forest Consumed (m)Baiting OnlyTree Removal OnlyBaiting and Tree RemovalHigh Initial DensityLow Initial DensityControl (monitoring alone)Figure 4.6: Length of forest consumed (left) and total MPB population(right) with different management types as the distance between source treesis increased. The simulation results are shown after a two year simulationinitiated with 5 source trees. The four management types tested are control(no management), baiting only, tree removal only, and combined baiting andtree removal.874.3. Effects of Managementmanagement. As expected, when the population density is low/source treesare far apart, baiting and tree removal over multiple years increases theimpact of MPB when compared to a single year of management activities.Thus, larger MPB densities and larger forest consumption will result fromconsecutive management activities of baiting and tree removal, if performedat low MPB densities. Comparing the results of baiting only, we find thatconsecutive baiting will increase the impact of MPB when source trees arefarther apart. The difference between single-year and consecutive baitingonly at small source tree distances is not significant.Management implications: Baiting and Tree RemovalAt a first glance our results indicate that tree removal only is a betterstrategy than baiting combined with tree removal. This is conditional onremoving MPB in the correct location. In the simulations we assumed thatthe tree removal only occured in the same spot as the tree removal in thebaiting only case. That is, we removed trees in an area which was rich inMPB and therefore had a large impact on the spread of MPB. If we hadsearched in the wrong area, and removed a smaller proportion of the MPBon the landscape, then the predicted output would be the same as in thecontrol case. Therefore, the success of tree removal only as a strategy mayvary between the curve for tree removal only and the control curve basedon the choice of search location. The advantage of baiting is that the MPBare drawn to a known site, and attacked trees are easier to find.The danger of using baiting is that if the MPB population is low orsource trees are spread far apart, baiting will allow the MPB to aggregateand survive where they would not have naturally been able to do so. In thiscase, baiting will increase the impact of MPB on the susceptible landscape.This effect of baiting was not well understood in past studies.The motivation for this study came fron the differences observed inMPB spread between the management and monitoring zones. The man-agement zone had the same area affected as the monitoring zone, but thelong-distance spread was lower. Based on the results of our simulations, we884.3. Effects of Management50 100 150 200010020030040050060070080050 100 150 20005000100001500050 100 150 200500600700800900100011001200130050 100 150 2000.511.522.5 x 104Distance Between Source Trees (m)Length of Forest Consumed (m)Baiting OnlyTree Removal OnlyBaiting and Tree RemovalHigh Initial DensityLow Initial DensityControl (monitoring alone)Total Beetle PopulationFigure 4.7: Length of forest consumed (left) and total MPB population(right) with different management types as the distance between source treesis increased. The simulation results are shown after a two year simulationwhere management occurs over both years, with 5 source trees initially. Thefour management types tested are control (no management), baiting only,tree removal only, and combined baiting and tree removal.894.3. Effects of Managementbelieve this could have occured because of the interaction of attraction tobait site, removal of trees, and increased MPB aggregation success. First,the baiting in each year concentrates the MPB movement, thus increas-ing the density and success of MPB in attacking susceptible trees on thelandscape. After removal of green-attack trees, any MPB missed would pro-duce the following year?s MPB population. These new MPB would be at alower density than in the monitoring zone (due to removal) and would beattracted to the new pheromone baits, both of which would reduce long-distance spread. The local impact in the management zone, however, wouldbe larger since baiting increases the success of MPB in attacking suscepti-ble trees. Thus, this yearly pattern of bait and removal over the landscapewould result in a more concentrated spread (less long-distance spread) ofMPB, which was exactly the difference observed between the managementand monitoring zones in Banff National Park.4.3.2 Prescribed BurningThe next management activity we tested was prescribed burning.In high-intensity fires, prescribed burning reduces the amount of suscepti-ble lodgepole pine available for the MPB and introduces unsuitable habitatbetween patches of susceptible landscape. Low-intensity fires have beenhypothesized to introduce more competition for reproduction, and trees inthese areas may be more attractive to MPB (Mary Reid, personal communi-cation). More competition through reproduction may occur through otherspecies attacking and reproducing within the same host tree, such as pineengravers. This can occur since partially damaged trees may become suit-able hosts for these other forest pests due to the damage induced by thefire. The reason that trees may be more attractive to MPB is not yet fullyunderstood. In this section we investigate how these two factors change thespread of MPB over the landscape.We simulated four cases: Control (No change between unburned andpartially damaged trees), higher attraction to partially damaged trees (at-tractive), lower reproduction in partially damaged trees (low growth), and904.3. Effects of Managementa combination of lower reproduction and higher attractiveness of partiallydamaged trees (attractive low growth). We chose a symmetric 1-dimensionallandscape that was burned (scorched) in the centre, then partially burned(damaged) outside of that, and finally was surrounded by undamaged sus-ceptible landscape. There were therefore 5 patches, all of the same size,in the order undamaged, partially damaged, scorched, partially damaged,undamaged. Each patch was chosen to be 100 m across, and the simulationsetup is displayed in Figure 4.8.Studies of real landscapes do not yet provide sufficient information forus to select a single initial condition as most realistic. Therefore, we simu-late three initial MPB distributions, with MPB emerging from focal habitat,damaged habitat, or matrix habitat. The MPB will emerge in: (1) The leftundamaged patch (0-100 m), (2) the left partially damaged patch (100-200m), and (3) the scorched patch (200-300 m). We use initial distributionsof MPB on the only one side (left) of the landscape to test scenarios ofMPB invading a landscape, and examining how MPB movement will re-spond to fragemented habitat that has been burned or clearcut. We choseMPB to start in each type of landscape (scorched, partially damaged, andundamaged) to examine the nesting and movement behaviour under ourfour different assumptions. We also test initial densities of MPB between31,250 MPB/ha and 250,000 MPB/ha. This corresponds to the output ofapproximately 3-25 source trees per hectare, assuming 10,000 MPB/tree perflight season [PB09].A summary of the fraction of susceptible forest consumed (m) against thesource tree density (per ha) is shown in Figure 4.9. More detailed outputsof the density of nesting MPB at the end of each year for the three initialconditions and different assumptions on the attraction and reproduction ofpartially damaged trees is found in Appendix B in Figures B.1, B.2, and B.3.A number of patterns emerge from the output of Figure 4.9. The initialMPB density is very important, as the proportion of trees in the landscapekilled by MPB increases at different rates with initial population size. Inall four assumptions on growth and attractiveness, the fastest increase in914.3. Effects of Management0 500100 200 300 400U D S D U0 500100 200 300 400U UU UC or SControlAssumptions on Partially Damaged TreesFigure 4.8: Simulation setup for prescribed burn. ?U? represents un-damaged stands, ?D? represents partially damaged stands, ?C? representsclearcut stands, and ?S? denotes scorched stands. In the control situa-tion (top), no additional assumptions are made upon the partially damagedstands, so they act as essentially the same as undamaged stands. Sincethere is no partially burned trees the centre can be thought of as essentiallya clearcut region in an undamaged forest (or a scorched region with no dif-ference between ?D? and ?U? from the MPB perspective). When additionalassumptions on the growth or attractiveness of partially damaged trees aremade (bottom), the undamaged stands differ from partially damaged stands,thus the stands around the scorched region, ?S?, have been designated aspartially damaged, ?D?. Each type of stand is chosen to be 100 m, with atotal domain size of 500 m.924.3. Effects of Management0 0.5 1 1.5 2 2.5x 10500.20.40.60.810 0.5 1 1.5 2 2.5x 10500.20.40.60.810 0.5 1 1.5 2 2.5x 10500.20.40.60.810 0.5 1 1.5 2 2.5x 10500.20.40.60.81Initial Population Density (MPB/ha)Fraction of Susceptibles Consumed Undamaged ForestPartially Damaged ForestScorched ForestInitial PositionRegular Growth Low GrowthNormal AttractionHigh AttractionFigure 4.9: Fraction of susceptible trees attacked in regions that have beenburned by a prescribed fire. The values vary with the initial source treedensity during emergence in the first year. Simulations were run for threeyears. MPB emergence either occured in the left undamaged forest, in thepartially damaged region, or in the scorched region. 2 different parameterswere varied in the simulation, attractiveness and growth of partially burnedtrees. The simulations were run with normal or low reproduction and normalor high attractiveness in partially burned trees.934.3. Effects of Managementproportion of susceptible landscape killed occured when the MPB emergedfrom the undamaged forest only. This is due to fact that susceptible treesare very easy to locate if MPB begin in the undamaged forest only. Thisrate was reduced if the MPB emerged in the damaged forest only, and wasminimum if the MPB were initially distributed in the scorched forest only.At first glance, this result is suprising since in the control case there is nodifference between the MPB growth or attractiveness in the damaged andundamaged forest, and thus both patches could be thought of as undamagedforest, as displayed in Figure 4.8. The reason that MPB emerging from theundamaged forest lead to a higher proportion of landscape killed is due tothe progression of MPB population spread. In the control case, if the MPBemerge from the undamaged region, they first nest in this region. In thesubsequent year, they nest in the partially damaged region, and then in thethird year they may successfully cross the gap of scorched forest to nest in theright damaged and undamaged landscapes. Note that this trend dependsupon the population size, as low initial populations may not successfullysurvive to year 2 or 3 and high initial population sizes may cross the scorchedforest in year 2. If, on the other hand, the MPB in the control case, emergefrom the damaged region, they first nest in this region. The next year,spread occurs primarily to the left into the undamaged region, instead ofspreading across the scorched region to the damaged and undamaged forestson the right. This leftward spread makes it more difficult for the MPB tospread across the scorched region in the subsequent year, as the distance tocross now includes the damaged region that was consumed in the first year.So the total dispersal distance becomes 200 m instead of 100 m. It is forthis reason that at lower population densities, MPB beginning in the leftundamaged forest will be more successful in killing susceptible trees in theright damaged and undamaged forests than the same population beginningin the damaged forest. The final case of the scorched region initial conditionis intuitively last since the MPB must disperse from the scorched region tothe damaged region before nesting in the first year. This dispersal processresults in additional mortality of MPB in the first year, when compared tothe other initial conditions. This trend in proportion of trees killed is similar944.3. Effects of Managementfor the other assumptions on reproduction and attraction in the damagedforest.The second pattern we observe is a decrease in the proportion of for-est consumed at higher densities of MPB in some cases. The reason forthis decrease is limitation of host trees. As the MPB population spreadsbetween the first and second year the number of susceptible trees availablein the second year is reduced by the length of forest consumed in year 1.At high MPB population densities, the MPB consumes a large portion offorest in the first year and in the second year, has less susceptible landscape(to the left of the scorched region) for reproduction. Since the susceptiblelandscape is smaller, total reproduction of MPB is reduced and there is alower total impact over the three year period. In other words, higher initialdensities of MPB are less successful in crossing the scorched area becausethe reproductive success in the first year creates strong resource limits inthe second year. The property of area limitation is conditional on the sizeof the scorched region and the size of each patch. If the scorched region wassmaller (larger) we would expect smaller (larger) area limitation effects onthe spread of MPB.When considering the case of lower growth we observe the intuitive resultthat the proportion of susceptible trees consumed is smaller than in thecontrol case, with all three initial conditions. If MPB emerge initially inthe undamaged region, the lower reproductive success prior to crossing thescorched area, leads to a lower overall impact on the susceptible landscape.Similarily, if MPB emerge from the scorched region, the spread in the secondyear is heavily affected by the decreased reproductive success in the first year.In contrast, the MPB emerging from the partially damaged region are notas significantly affected as the other two initial conditions. This is becausethe largest limitation on the spread of MPB, with this initial condition, isthe spread in the third year. The reproduction prior to this major dispersalevent (between years 2 and 3) occurs in the undamaged forest, so the effectsof the assumption of lower growth in the partially damaged forest do notaffect the MPB impact with this initial condition as significantly as the MPBimpact with the other two initial conditions.954.3. Effects of ManagementThe assumption of higher attractiveness changes the simulations in twoimportant ways. First, the MPB starting in the undamaged forest are at-tracted to the damaged area, so as the attractiveness of the damaged areaincreases, the initial conditions of starting in the damaged or undamagedareas become essentially the same. The second major difference is thatthe higher attractiveness enhances the success of the MPB in crossing thescorched area. Thus, all initial conditions show increased mortality of sus-ceptible trees on the landscape. That is, MPB populations of smaller sizeare able to persist, and have a larger impact on the susceptible landscape.Additionally, there is no longer a decrease in the proportion of susceptibletrees consumed at high MPB densities. Since the scorched area is not asdifficult to cross, MPB in this simulation cross the scorched region earlier(in the second year), and thus are not subject to a limiting susceptible treedensity. These patterns were not observed in simulations with medium,rather than high, attraction to the partially damaged forest. Therefore, theresults in Figure 4.9 are critically dependent on the degree of attraction ofthe damaged region.Finally, the combination of high attractiveness and lower reproductionin damaged trees leads to curves of proportion of susceptibles consumedthat are intermediate between the curves with high attraction (and nor-mal growth) and low growth (with normal attraction). We predicted thatincreased attractiveness to the partially damaged area (which yield lowergrowth) may reduce the susceptible fraction of trees consumed in compari-son to the control case. This is true when the MPB emerge initially from theundamaged area. For the other two initial conditions, however, the higherattractiveness of the damaged region (causing populations to successfullycross the scorched area) was a larger factor than the lower reproduction inthese regions. Therefore, the fraction of susceptibles consumed when MPBemerge from the scorched or damaged region was larger with high attractive-ness and lower growth than in the control case. In addition, the combinationof lower growth with high attractiveness allowed for smaller area limitationeffects (than in the control case), as the scorched area was easier to cross.Therefore, in terms of the success of the MPB population, there is a964.3. Effects of Managementtrade-off between the damaged trees being more attractive and having alower growth rate. The lower growth rate decreased the MPB population inthe following year, but the higher attractiveness allows the beetles to spreadfurther across patches of scorched timber to reach partially burned stands.The balance point is obviously dependent upon the magnitude of increasedattractiveness and decreased reproduction. Additionally, if the size of allof the patches is increased, the results will differ slightly because the MPBdensity required to cross the gap (scorched area) will also increase.Management Implications: Prescribed BurningOur simulations indicate that lower reproduction of MPB in damagedtrees always results in smaller spread MPB and decreased damage to suscep-tible trees. Additionally, scorched regions of the forest are both unsuitablefor reproduction, and a barrier to the spread of MPB due to dispersal limita-tion and limitation of susceptible host density. Therefore, if MPB do indeedhave a lower reproduction in damaged trees, our simulations would supportthe use of prescribed burning as a management activity.The control simulations assume that there is no difference between thedamaged and undamaged stands. In this case, the scorched area acts thesame as a clearcut introduced in the center of a two-patch system (Figure4.8). As a result, we can compare simulations of clearcutting to prescribedburning by comparing simulations of the control scenario to the other threescenarios.From Figure 4.9, we can see that if damaged stands have lower MPBreproduction, prescribed burning reduces the impact of MPB more thanclearcutting the area equivalent to the scorched area. In contrast, if thedamaged trees are more attractive to beetles, the scorched area becomeseasier to cross. Therefore, if the damaged area is highly attractive, it maybe a better strategy to remove trees through clearcutting instead of using aprescribed burn.In the case where the damaged region is both attractive and has lowerreproduction, the results will be dependent on the magnitude of decreased974.4. Characteristics of MPB spreadreproduction and increased attractiveness. Our results are critically depen-dent on assuming the partially damaged trees are highly attractive. If thepartially damaged trees only have a medium level of attractiveness, no sig-nificant changes from the normal levels of attractiveness are observed. Wetherefore conclude that, if damaged trees are both attractive and have lowerreproductive output, and the level attractiveness is not too high, the bestmanagement strategy will always be to burn a central patch rather thanclearcut the equivalent area to reduce MPB impact.4.4 Characteristics of MPB spread4.4.1 Interaction Between Source TreesMPB from different source trees interact, so we considered how faraway source trees needed to be for MPB from separate source trees to in-fluence one another. We simulated the interaction between source treesusing a periodic environment with three source trees, all equidistant fromeach other. The source tree densities were chosen such that MPB popula-tions would not persist with only a single source tree because of the need toovercome host defences. Thus, the interaction between MPB from differentsource trees was crucial to the survival and reproduction of the population.We found that if the distance between MPB source trees was greater thanor equal to 325 m, the MPB population would not persist (Figure 4.10).4.4.2 Clearcut LengthWe were interested in what width of clearcut would be needed forMPB to be unsuccessful in crossing the region. We inserted a clearcut inthe middle of the domain and observed the subsequent spread of the MPBpopulation. We compared the source tree density (per hectare) againstthe clearcut distance needed to prevent MPB spread across the clearcut inone year, as shown in Figure 4.11. The clearcut distance needed increasesapproximately linearly as the source tree density increases. At 25 sourcetrees per hectare, the minimum clearcut length is 1 km.984.4. Characteristics of MPB spread0 100 200 300 400 50000.20.40.60.81Distance Between Source Trees (m)Fraction of Susceptibles consumedFigure 4.10: Fraction of susceptibles consumed after one year with varyingdistances between source trees. Trees equidistant over the periodic domain,which is equivalent to an infinite domain with source trees spaced at thedistance specified on the x-axis.994.4. Characteristics of MPB spread0 5 10 15 20 2502004006008001000Source Trees per HectareDistance Clearcut Needed (m)Figure 4.11: Clearcut distance needed for MPB populations to not success-fully cross gap. Densities of source trees is varied between 3 and 25 sourcetrees per hectare.1004.5. Discussion4.4.3 Source Tree MappingUsing our multi-year model, we simulated source tree density perhectare needed in a cluster of attack for the MPB to overcome the Alleethreshold. That is, at low densities, MPB are not able to mass-attack trees,and therefore do not have successful reproduction. At higher densities, theMPB population is able to overcome the Allee threshold, and can grow toepidemic levels [HP08].To find the Allee threshold, we considered a single group of source treesand simulated subsequent MPB emergence and attack. We mapped theinitial source tree density from one year to the next. We considered a domainconsisting of uniform susceptible trees and had MPB emerging at the centreof the domain. We ran simulations over a four year period and produceda year-to-year mapping of source tree densities in Figure 4.12. This plotshows that to overcome the Allee effect, there must be approximately 2.3source trees per hectare. Above this threshold, the source trees grow indensity from year to year, while below the threshold the source tree densitydecreases to 0. Thus, we can predict the source tree threshold at whichthe MPB population will transition from endemic to epidemic levels. MPBpopulations may persist at levels below the Allee threshold predicted becauseof stochastic factors which are not captured by our model framework.4.5 DiscussionWith our multi-year spatially explicit mathematical model, we stud-ied the effectiveness of both direct and indirect management strategies. Thedirect management strategies we investigated were the application of treeremoval with and without a trap baiting component. The indirect manage-ment strategies considered were clearcutting and prescribed burning.The direct management strategy of baiting and tree removal was found toalways decrease the MPB population density if the initial MPB populationdensity is sufficiently large. At low initial MPB densities, the baiting compo-nent increases the success of initial aggregation, leading to a higher density1014.5. Discussion0 5 10 15 200510152025source trees per ha in previous yearsource trees per ha in current yearFigure 4.12: Number of source trees per hectare mapped from one year tothe next. The identity line indicates points where the source trees mappedfrom one year to the next are constant. The simulation was run for 4 yearsto obtain these mapping results.1024.5. Discussionof nesting MPB and a larger impact on the susceptible landscape. Depend-ing on the success rate in detection (and removal) of attacked trees, thebaiting component of this management strategy allows populations of MPBwhich would otherwise not survive and reproduce to successfully continuespreading in subsequent years. These results suggest that caution should beused in the application of baiting at very low densities of emerging MPB, orwhen initial source trees are far apart.We propose that the baiting and tree removal process is what led to theunexpected differences between the patterns of MPB attack in the manage-ment and monitoring zones in Banff National Park [TR08]. The manage-ment and monitoring zones had equal areas of susceptible landscape killed,but MPB spread in the management zone was less. The mechanism ofbaiting would have concentrated the MPB attack, thereby reducing long-distance spread, but also leading to higher MPB attack success rates inthe local region of the bait. The tree removal component of this strategylikely reduced local MPB population densities each year enough to reducelong-distance spread, but not enough to reverse the epidemic. Thus, thisdifference in the two zones could be largely attributed to the mechanismsassociated with baiting and tree removal. This observation is consistent withthe results of Gray and Borden [GB89], who found that baiting and tree re-moval does slow the spread of MPB, but has the side-effect of concentratingthe population.When comparing the success of tree removal with and without the baitingcomponent, the results depends critically on the success rate of finding high-densities of MPB-attacked trees. In particular, without baiting, the successrate depends on both the choice of search area and the proportion of attackedtrees detected in that area. If the tree removal occurs in a region with highattack densities, then tree removal without baiting will always lead to lowerMPB population densities and smaller areas of impacted landscape. Green-attack trees can be very difficult to find, however, even if the attack areais generally known (Mary Reid, personal communication). There is thus atrade-off between the increased success rate in locating green-attack treeswhen baiting is used, and the simultaneous increased reproductive success1034.5. Discussionof the MPB. The threshold MPB density at which the benefits of baiting asa management strategy outweigh the drawbacks depends also on our abilityto anticipate the location of MPB attacks when baiting is not used.When management activities are sustained over multiple years the ulti-mate effect depends on the initial density of MPB. If the density is high,baiting and tree removal is a more successful strategy when sustained overmultiple years. If MPB density is low, baiting and tree removal over mul-tiple years ultimately significantly increases the impact of MPB. Baitingonly over multiple years always increases the impact of MPB, whereas treeremoval only over multiple years always decreases the impact of MPB.Direct management strategies employed in the past have met with lim-ited success in supressing the MPB population (See Safranyik and Wilson[SW06] and accompanying references). This was hypothesized to be dueto three problems. First, detection of early MPB attack locations and as-sessment of population growth is difficult. Our simulation results highlightthat management strategies are most effective when early MPB locationsare known, and that assessment of population size is important for choiceof an appropriate management strategy. The second problem is that man-agement efforts were not applied over multiple years. Our simulations showthat significant results can be achieved only if management activities con-tinue over multiple years. The final issue with past management strategiesis that often only a portion of the affected area was targeted. According toour simulations, the success of baiting and tree removal is very dependentupon the search width, and targeting an area with a radius of 100 m ormore of the affected trees is critical. A case study of successful suppressionof MPB occured in Banff in the 1940?s when all three of these shortcomingswere addressed [HM45].Using the indirect method of clearcutting or burning a portion of thelandscape to create a population gap will consistently reduce the spreadand population density of MPB. Prescribed burning can be used to createan intensely scorched region in which MPB cannot nest, with a surroundingregion of partially burned trees. Which method is optimal depends upon theattraction level of MPB and the reproductive success in trees that have been1044.5. Discussionpartially burned. Further work could be done to investigate how introducingrepulsive and attractive pheromones will change these results.If partially damaged trees lead to lower reproduction but have equalattractiveness to MPB when compared to undamaged trees, than creatinga forest gap using prescribed burning is more effective than clearcuttingin reducing the population density and impact of MPB. If partially dam-aged trees are more attractive and have equal reproduction for MPB whencompared to undamaged trees, then the MPB population density can beincreased by burning, compared to the clearcutting case. This result is dueto the increase aggregation successs in partially damaged regions, and theincrease in dispersal success of beetles across gaps of scorched trees. If par-tially damaged trees are both more attractive and have lower reproductionfor MPB, the overall result will depend on the relative magnitudes of attrac-tiveness and reproduction. If attractiveness is sufficiently low, burning theregion will always result in a lower impact (smaller MPB population den-sity and smaller length of forest consumed) than if the equivalent scorchedregion had been clearcut.One counterintuitive result arising from our simulations shows that ifundamaged forest is limiting, it is possible for MPB to become corneredthere. Consequently, a lower initial population of MPB can spread acrossa clearcut or scorched region with higher success than a larger initial pop-ulation of MPB. This pattern results when MPB do not cross the gap andconsume too much of the susceptible habitat on one side of the gap. Thefollowing year, the population is only able to reproduce in the remainingamount of habitat on the one side. This limited habitat leads to a lowerreproduction and thus a smaller overall spread and density. This effect couldbe enhanced if baiting was done in the small region of undamaged forest,which could severely limit the spread to nearby habitat.It would be an interesting avenue of future work if management could finda way to use this area-limiting property in an efficient way. For instance, ifwe know the population level and the position of the MPB population, a burnor clearcut could be made to reduce or stop the MPB spread. Unfortunately,trees do not show the obvious red-top from mortality until after MPB have1054.5. Discussionemerged from the given tree (position unknown) and population levels ofMPB are difficult to determine [SW06]. Thus, this may be very difficultto put into practise. This property, however, may lead to some interestingresults where MPB in management zones at lower densities may lead tohigher impacts than MPB at higher densities.When comparing the indirect results of clearcutting and burning wemade several assumptions which could be investigated further. First, weassumed that the scorched area (very damaged trees) was unusable by MPBfor reproduction and was not attractive to the MPB (Mary Reid, personalcommunication). These assumptions are not currently verified but researchis being done to explore these assumptions further. Second, in interpret-ing the scorched region in the control landscape as equivalent to a clearcutregion, we assumed that movement in a clearcut patch would be identicalto movement in a scorched region. It is possible however, that diffusionin the two patches is different [Rei09] (Chapter 3). A scorched region stillcontains the tree trunks necessitating, possibly, more tortuous beetle flightpaths [Rei09], and thus a different diffusion speed than in a clearcut. Fi-nally, our simulations assumed that dispersal was identical in the scorched,partially damaged, and undamaged stands (apart from the attractive effectsof kairomone and pheromone). In a study by Reid [Rei07], however, it wasfound that movement may increase in burned regions.Our simulation results were based on a very synchronous emergence.Other emergence rates and times were tested, and it was found that loweremergence rates led to higher initial source densities needed for the MPB toreproduce. Thus, the population densities reported in this chapter will varysignificantly based on the emergence profile of the MPB, which can varydramatically based on temperature [GPLB04]. Intuitively, lower emergencerates, and less synchronous emergence leads to a lower success in MPB re-production. Our results can thus be considered worst case scenarios fromthe point of view of emergence. For more accurate predictions of aggrega-tion, growth and spread rates more realistic emergence functions should beconsidered [GPLB04, RPBN12]. For an extremely realistic emergence func-tion, age-dependent differential equations could be used to predict emergence1064.5. Discussiontimes and densities [CP11].When we ran multi-year simulations, we were able to map the den-sity of source trees per hectare from one year to the next. We found thatthe mapping function had an unstable critical point at 2.3 source trees perhectare. Below this threshold, the MPB population decreases to extinction.We interpret ?extinction? in this model as endemic infestation. Above thecritical threshold, the population grows to epidemic levels. This behaviouris consistent with the result found in Chapter 2 predicting endemic popula-tion densities when attack levels are below 0.67 source trees per hectare, andepidemic population densities above 4 source trees per hectare. Incipient epi-demic levels, which transition between endemic and epidemic levels, were inbetween 0.67 and 4 source trees per hectare. A extension of this work couldconsider the Allee threshold present in the discrete portion of our model andinvestigate the effects of the aggregation (due to MPB pheromone) on thisthreshold.Our simulation model predicts that if source trees are spread more than325 m apart, there will be no population of MPB the following year, assum-ing that a single source tree of MPB is not a sufficient density to successfullyreproduce to the following year. This result shows an interesting correspon-dence with previous pattern formation work (Chapter 2). In this work, MPBattack clusters were found to be spaced 353 m on average, which is slightlylarger than the 325 m predicted by our model where source trees would havelittle correspondence. Therefore, this pattern formation wavelength couldresult from the attack clusters being too far apart to overlap and interact.In other words if source trees are too far apart spot formation may arisewith no spread.Our simulations of MPB population spread in the presence of a clearcutshow that a clearcut width of 1 km is sufficient to prevent spread of MPBpopulations as large as 25 source trees/hectare. We caution, however, thatour results are predicated on the assumption that dispersal occurs at thesame speed in clearcut and forested regions. Our result is thus only appli-cable to local MPB dispersal behaviour. It is entirely possible that withincreased wind speeds in large clearcut regions, beetles can be carried across1074.5. Discussionwith relative ease, or are more readily lifted into high air currents that carrythe beetles long distances over the forest canopy [dCA12]. Our model doesnot apply to this mode of dispersal and further studies are needed to de-termine the effectiveness of large clearcuts as barriers to dispersal. Indeedour results in Chapter 3 found that diffusion may be higher in fragmentedhabitats.A further extension of our work would be to consider a mix of directand indirect control strategies [BS03]. We suggest that a successful ap-proach might be to use a prescribed burn to limit MPB spread and growththrough the introduction of gaps, and then use baiting and tree removal toseverely limit the spread of MPB. The introduction of gaps in the susceptiblelandscape will limit the spread of MPB, and the baiting and tree removalconcentrates, but also further limits the ability of the MPB to cross gaps inthe susceptible landscape.From an economic perspective, our models do not factor in the cost ofeach management activity. In addition to MPB dispersal studies, optimalcontrol models factoring in the cost of prescribed burning, clearcutting, andbaiting and tree removal are needed to determine the most cost-effectivecombination of strategies.We note here that burning, clearcutting and targeted tree removal allengender changes to the landscape that are irreversible on a short timescale(decades), and can also be very costly. Simulation approaches like ourscan be used to predic the influence of both direct and indirect managementstrategies. We have been able to successfully explain the non-intuitive MPBspread rates and patterns in response to management activities in Banff NP,and thus provide an informative context for future management decisions.More modelling efforts of this sort are needed, so that we can minimizethe use of destructive management tactics, especially as the MPB epidemiccontinues and spreads into the boreal forest [dCA12].108Chapter 5ConclusionThe ultimate goal of this body of work was to understand the dispersaland attack patterns of the MPB under the influence of management activi-ties. We have succesfully accomplished this goal on the stand spatial scale(1 ha to 1 km2) through the three projects described in this thesis. All ofthese projects use a original spatially explicit mathematical model whichhas been developed based on the existing theoretical and empirical work[WP97, RB83].In Chapter 2, we used our model framework to predict the spacing be-tween MPB attack clusters at incipient epidemic densities. We used MPBattack data from the Sawtooth National Recreation Area to verify our the-oretical predictions, providing strong support for our modelling approach.This gave us possible insight into source tree densities at which MPB tran-sition population levels from endemic to epidemic.A key aspect of MPB dispersal dynamics is movement through frag-mented habitats. We thus considered small-scale regions to develop a de-tailed understanding of movement across heterogeneous regions. This ap-proach led us to the combination of experimental and theoretical work foundin Chapter 3. The dispersal behaviour of the MPB through fragmented habi-tats is largely unknown [BAL93], and therefore our work testing differenthypotheses on MPB movement in various habitats is crucial. We were ableto link theoretical and experimental predictions to select the most likelymechanisms of movement of MPB in fragmented landscapes.The culmination of this work appears in Chapter 4 where we investi-gated the impacts of various direct and indirect management activities onthe spread and impact of the MPB. We determined what is the optimal di-rect management strategy as a function of MPB population density and the109Chapter 5. Conclusionsuccess in detection (and removal) of green-attack trees. We furthermoredetermined that the best choice of indirect management strategy changesbased on the assumptions made regarding the reproduction (of MPB) andattractiveness of partially burned trees. This aspect of MPB biology is asyet unknown, but our work shows that it is a critical consideration in MPBmanagement. Finally, our simulation model was able to provide valuablemanagement predictions for the spacing between source trees, clearcut dis-tance needed, and clustered source tree density required to reduce the MPBpopulation density to endemic levels. Our theoretical framework successfullyexplained the non-intuitive MPB spread rates and patterns observed in re-sponse to management activities in Banff National Park, and thus providesan informative context for future management decisions.The two most important extensions of our model framework are theeffect of temperature on the MPB and the variable resistance of host trees.All stages of the MPB life cycle explicitly depend upon the surroundingtemperature [SW06, GPLB04, CP11]. Recent work by Powell et al. [PB09]used a temperature-dependent model to predict population growth basedon the ambient temperature. This incorporation of temperature into theemergence dynamics of the MPB would lead to a more accurate predictionof MPB population growth and spread in real landscapes.In this work, we restricted our attention to landscapes where trees withina patch were all exactly the same in terms of their attractiveness and re-sistance to MPB. Our only variable quantity was susceptible tree density.However the importance of weakened trees for nucleation of MPB attack iswell-recognized in the literature [WP97, NL08]. The successful growth of aMPB population to epidemic levels can begin at these ?nurse trees?, whichare weakened by such factors as drought or lightning strikes. Incorporationof this factor into the model framework could take place by allowing kp tospatially vary.In summary, this body of work has led us to new insights into the move-ment of MPB in fragmented habitats, the patterns of MPB attack, andthe effects of direct and indirect management activities on MPB impactand spread. We have been able to provide valuable management predic-110Chapter 5. Conclusiontions, without engaging any costly and irreversible changes to the landscape.More extensive use of modelling efforts akin to our study are needed to in-form management decisions in the future, especially as the MPB continueto spread in Alberta [dCA12].111Bibliography[BAL93] B.J. Bentz, G.D. Amman, and J.A. Logan. A critical assessmentof risk classification systems for the mountain pine beetle. ForestEcology and Management, 61(3/4):349?366, 1993. ? pages 4,33, 67, 109[BB91] E.O. Budrene and H.C. Berg. Complex patterns formed bymotile cells of escherichia-coli. Nature, 349(6310):630?633, 1991.? pages 5[BB95] E.O. Budrene and H.C. Berg. Dynamics of formation of symmet-rical patterns by chemotactic bacteria. Nature, 376(6535):49?53,1995. ? pages 5[Ben06] B.J. Bentz. Mountain pine beetle population sampling: infer-ences from lindgren pheromone traps and tree emergence cages.Can. J. For. Res., 36:351?360, 2006. ? pages 66, 67[BLB+06] H.J. Barclay, C. Li, L. Benson, S. Taylor, and T. Shore. Effectsof fire return rates on traversability of lodgepole pine forests formountain pine beetle and the use of patch metrics to estimatetraversability. Natural Resources Canada, 2006. ? pages 69[BPBL00] Z. Biesinger, J. Powell, B. Bentz, and J. Logan. Direct and indi-rect parametrization of a localized model for the mountain pinebeetle-lodgepole pine system. Ecological Modelling, 129:273?296, 2000. ? pages 9, 37, 49, 67, 73[BPL96] B.J. Bentz, J.A. Powell, and J.A. Logan. Localized spatial andtemporal attack dynamics of the mountain pine beetle in lodge-112Bibliographypole pine. Intermountain Research Paper, 494, 1996. ? pages6, 14[BS03] J.E. Brooks and J.E. Stone. Mountain Pine Beetle Symposium :Challenges and Solutions. Natural Resources Canada, CanadianForest Service, Pacific Forestry Centre, Victoria, BC, 2003. ?pages 6, 108[BSW84] A.A. Berryman, N.C. Stenseth, and D.J. Wollkind. Metastabil-ity of forest ecosystems infested by bark beetles. Res. Popul.Ecol., 26:13?29, 1984. ? pages 3[Bur77] D.G. Burnell. A dispersal-aggregation model for mountain pinebeetle in lodgepole pine stands. Res. Popul. Ecol., 19:99?106,1977. ? pages 2, 6[CC82] S.P. Courtney and S. Courtney. The edgeeffect in butterflyoviposition : causality in anthocharis cardarnines and relatedspecies. Ecological Entomology, 7:131?137, 1982. ? pages 33,67[CDO07] D.S. Chapman, C. Dytham, and G.S. Oxford. Landscape andfine-scale movements of a leaf beetle: the importance of bound-ary behaviour. Oecologia, 154(1):55?64, November 2007. ?pages 32, 33, 66[CMM+10] L.R. Carrasco, J.D. Mumford, A. MacLeod, T. Harwood,G. Grabenweger, A.W. Leach, J.D. Knight, and R.H.A. Baker.Unveiling human-assisted dispersal mechanisms in invasive alieninsects: Integration of spatial stochastic simulation and phenol-ogy models. Ecological Modelling, 221:2068?2075, 2010.? pages30[CP11] C.A. Cobbold and J.A. Powell. Evolution stabilises the synchro-nising dynamics of poikilotherm life cycles. Bull. Math. Biol.,73:1052?1081, 2011. ? pages 107, 110113Bibliography[CPB12] B.A. Crabb, J.A. Powell, and B.J. Bentz. Development and as-sessment of 30-m pine density maps for landscape-level modelingof mountain pine beetle dynamics. USDA FS Rocky MountainResearch Station Research Note, 2012. ? pages iv, 20[dCA12] H.-M.C. de la Giroday, A.L. Carroll, and B.H. Aukema. Breachof the northern rocky mountain geoclimatic barrier: initiationof range expansion by the mountain pine beetle. J. Biogeogr.,39:1112?1123, 2012. ? pages 65, 108, 111[ER04] C.M. Elkin and M.L. Reid. Attack and reproductive success ofmountain pine beetles (Coleoptera: Scolytidae) in fire-damagedlodgepole pines. Environ. Entomol., 33(4):1070?1080, 2004. ?pages 69[Fah07] L. Fahrig. Non-optimal landscapes animal movement in human-altered landscapes. Functional Ecology, 21(6):1003?1015, 2007.? pages 32[FCC99] W.F. Fagan, R.S. Cantrell, and C. Cosner. How habitatedges change species interactions. The American Naturalist,153(2):165?182, 1999. ? pages 32[GB89] D.R. Gray and J.H. Borden. Containment and concentrationof mountain pine beetle (Coleoptera: Scolytidae) infestationswith semiochemicals: validation by sampling of baited and sur-rounding zones. Journal of Economic Entomology, 82:1399?1405, 1989. ? pages 103[GGG80] D.R. Geiszler, V.F. Gallucci, and R.I. Gara. Modeling the dy-namics of mountain pine beetle aggregation in a lodgepole pinestand. Oecologia (Berl.), 46:244?253, 1980. ? pages 2, 6[GH08] J.G.P Gamarra and F. He. Spatial scaling of mountain pine bee-tle infestations. Journal of Animal Ecology, 77:796?801, 2008.? pages 6114Bibliography[GPLB04] E. Gilbert, J.A. Powell, J.A. Logan, and B.J. Bentz. Compar-ison of three models predicting developmental milestones givenenvironmental and individual variation. Bulletin of Mathemat-ical Biology, 66:1821?1850, 2004. ? pages 2, 29, 37, 106, 110[GTLS] T. Gauduchon, R.C. Tyson, F. Lutscher, and S. Strohm. Theeffect of habitat fragmentation on populations with cyclic dy-namics and a discontinuous density across patch interfaces. Inpreparation. ? pages 32, 67[Hae05] J.W. Haefner, editor. Modeling biological systems. Springer,New York, NY, 2005. ? pages 29, 126[HC03] K.J. Haynes and J.T. Cronin. Matrix composition affects thespatial ecology of a prairie planthopper. Ecology, 84(11):2856?2866, 2003. ? pages 32, 33, 66[HC06] K.J. Haynes and J.T. Cronin. Interpatch movement and edgeeffects : the role of behavioral responses to the landscape matrix.Oikos, 113:43?54, 2006. ? pages 32, 33, 66[HFSL06] J. Hughes, A. Fall, L. Safranyik, and K. Lertzman. Modeling theeffect of landscape pattern on mountain pine beetles. TechnicalReport BC-X-407, Natural Resources Canada, Canadian For-est Service, Pacific Forestry Centre, Victoria, British Columbia,2006. ? pages 3, 6, 12, 53, 64[HLZ+08] B. H. Aukema, A. L. Carroll, Y. Zheng, J. Zhu, K. F. Raffa,R. D. Moore, K. Stahl, and S. W. Taylor. Movement of outbreakpopulations of mountain pine beetle: influences of spatiotempo-ral patterns and climate. Ecography, 31:348?358, 2008. ? pages6, 29[HM45] G.R. Hopping and W.G. Mathers. Observations on outbreaksand control of the mountain pine beetle in the lodgepole pinestands of western Canada. The Forestry Chronicle, pages 1?11,June 1945. ? pages 104115Bibliography[HP08] J. Heavilin and J. Powell. A novel method of fitting spatio-temporal models to data, with applications to the dynamics ofmountain pine beetles. Natural Resource Modeling, 21(4):489?524, 2008. ? pages 3, 6, 101[HP09] T. Hillen and K. Painter. A user?s guide to pde models forchemotaxis. J. Math. Biol., 58(1):183?217, 2009. ? pages 12[Hug02] J. Hughes. Modeling the effect of landscape pattern on mountainpine beetle. PhD thesis, Simon Fraser University, 2002. ? pages33, 34, 65, 66[Jon77] R.E. Jones. Movement patterns and egg distribution in cabbagebutterflies movement patterns and egg distribution. Journal ofAnimal Ecology, 46(1):195?212, 1977. ? pages 33, 67[LeV97] R.J. LeVeque. Wave propogation algorithms for multidimen-sional hyperbolic systems. Journal of Computation Physics,131:327?353, 1997. ? pages v, 44[LeV07] R.J. LeVeque. Finite Difference Methods for Ordinary and Par-tial Differential Equations. SIAM, Philadelphia, PA, 2007. ?pages 44, 45, 46, 48[LP00] M.A. Lewis and S. Pacala. Modeling and analysis of stochasticinvasion processes. J. Math. Biol., 41:387?429, 2000. ? pages30[LT96] R.J. LeVeque and R.C. Tyson. Clawpack version 2.0, 1996. ?pages v[LWBP98] J.A. Logan, P. White, B.J. Bentz, and J.A. Powell. Model anal-ysis of spatial patterns in mountain pine beetle outbreaks. The-oretical Population Biology, 53:236?255, 1998. ? pages 6[MAD10] S. Matsumura, R. Arlinghaus, and U. Dieckmann. Foragingon spatially distributed resources with sub-optimal movement,116Bibliographyimperfect information, and travelling costs: Departures from theideal free distribution. Oikos, 119:1469?1483, September 2010.? pages 32[MKW03] C.W. McHugh, T.E. Kolb, and J.L. Wilson. Bark beetle attackson ponderosa pine following fire in northern Arizona. Environ.Entomol., 32(3):510?522, 2003. ? pages 69[ML] G.A. Maciel and F. Lutscher. How individual movement re-sponse to habitat edges affects population persistence and spa-tial spread. Am. Nat., in press. ? pages 67[MP91] R.G. Mitchell and H.K. Preisler. Analysis of spatial patterns oflodgepole pine attacked by outbreak populations of the moun-tain pine beetle. Forest Science, 37(5):1390?1408, 1991.? pages6[MSL+00] R.N. Mack, D. Simberloff, W.M. Lonsdale, Evans H., M. Clout,and F. Bazzaz. Biotic invasions: Causes, epidemiology, conse-quences, and global control. Issues in Ecology, 5:1?20, 2000. ?pages 31[Mur77] S.J. Muraro. Preliminary results of using fire to control moun-tain pine beetle. Canadian Forestry Service, Pacific ForestryCentre, 1977. ? pages 69[Mur03] J.D. Murray. Mathematical Biology. Springer-Verlag, Berlin,Heidelberg, 2003. ? pages 5, 12, 13, 30[NL08] W.A. Nelson and M.A. Lewis. Connecting host physiology tohost resistance in the confer-bark beetle system. Theor Ecol,1:163?177, 2008. ? pages 110[NPL+08] W.A. Nelson, A. Potapov, M.A. Lewis, A.E. Hundsdorfer, andF. He. Balancing ecological complexity in predictive models: areassessment of risk models in the mountain pine beetle system.Journal of Applied Ecology, 45:248?257, 2008. ? pages 3117Bibliography[NSS05] R.K. Nagle, E.B. Saff, and A.D. Snider, editors. Fundamentalsof Differential Equations and Boundary Value Problems. Pear-son, Boston, Ma, 2005. ? pages 15[Ova04] O. Ovaskainen. Habitat-specific movement parameters esti-mated using markrecapture data and a diffusion model. Ecology,85(1):242?257, January 2004. ? pages 32[PB09] J. Powell and B. Bentz. Connecting phenological predictionswith population growth rates for mountain pine beetle, an out-break insect. Landscape Ecology, 24:657?672, 2009. ? pages 2,19, 29, 49, 91, 110[PD11] L. Pe?rez and S. Dragic?evic?. ForestSimMPB: A swarming intelli-gence and agent-based modeling approach for mountain pinebeetle outbreaks. Ecological Informatics, 6:62?72, 2011. ?pages 6[PJLB00] J.A. Powell, J.L. Jenkins, J.A. Logan, and B.J. Bentz. Sea-sonal temperature alone can synchronize life cycles. Bulletin ofMathematical Biology, 62:977?998, 2000. ? pages 2[PKW+00] J. Powell, B. Kennedy, P. White, B. Bentz, J. Logan, andD. Roberts. Mathematical elements of attack risk analysis formountain pine beetles. Journal of Theoretical Biology, 204:601?620, 2000. ? pages 1[PL90] A.D. Polymenopoulos and G. Long. Estimation and evaluationmethods for population growth models with spatial diffusion:dynamics of mountain pine beetle. Ecological Modelling, 51:97?121, 1990. ? pages 3, 6[PL05] J.A. Powell and J.A. Logan. Insect seasonality - circle mapanalysis of temperature-driven life cycles. Theoretical Popula-tion Biology, 67:161?179, 2005. ? pages 2, 29118Bibliography[PLB96] J.A. Powell, J.A. Logan, and B.J. Bentz. Local projections fora global model of mountain pine beetle attacks. J Theor Biol,179(3):243?260, 1996. ? pages 73[PLBW02] M. Peltonen, A.M. Liebhold, O.N. Bjornstad, and D.W.Williams. Spatial synchrony in forest insect outbreaks : Roles ofregional stochasticity and dispersal. Ecology, 83(11):3120?3129,2002. ? pages 6[PM10] S Petrovskii and K McKay. Biological invasion and biologicalcontrol: A case study of the gypsy moth spread. Aspects ofApplied Biology, 104:37?48, 2010. ? pages 30[PMW98] J.A. Powell, T. Mcmillen, and P. White. Connecting a chemo-tactic model for mass attack to a rapid integro-difference emula-tion strategy. SIAM Journal of Applied Mathematics, 59(2):547?552, 1998. ? pages 3, 6, 7, 8, 40, 71[RB83] K.F. Raffa and A.A. Berryman. The role of host plant resis-tance in the colonization behaviour and ecology of bark beetles(Coleoptera: Scolytidae). Ecological Monographs, 53(1):27?49,1983. ? pages 8, 9, 37, 71, 73, 109[RB86] K.F. Raffa and A.A. Berryman. A mechanistic computer modelof mountain pine beetle populations interacting with lodgepolepine stands and its implications for forest managers. Forest Sci.,32(3):789?805, 1986. ? pages 2[Rei07] M. Reid. Environmental effects on host selection and dispersalof mountain pine beetle. Natural Resources Canada, 2007. ?pages 69, 106[Rei09] T.G. Reid. Dispersal studies of mountain pine beetles. Master?sthesis, University of Calgary, 2009. ? pages 34, 43, 66, 67, 106[RFSS04] W.G. Riel, A. Fall, T.L. Shore, and L. Safranyik. A spa-tiotemporal simulation of mountain pine beetle impacts on the119Bibliographylandscape. Technical Report BC-X-399, Natural ResourcesCanada, Canadian Forest Service, Pacific Forestry Centre, Vic-toria, British Columbia, 2004. ? pages 2, 3, 6[RH98] L.C. Remer and S.B. Heard. Local movement and edge effectson competition and coexistence in ephemeral-patch models. TheAmerican naturalist, 152(6):896?904, December 1998. ? pages33, 67[RNJ+09] C. Robertson, T.A. Nelson, D.E. Jelinski, M.A. Wulder, andB. Boots. Spatial-temporal analysis of species range expansion:the case of the mountain pine beetle, Dendroctonus ponderosae.Journal of Biogeography, 36(8):1446?1458, 2009. ? pages 6[RPBN12] J. Regniere, J. Powell, B. Bentz, and V. Nealis. Effects of tem-perature on development, survival and reproduction of insects:experimental design, data analysis and modeling. J. InsectPhysiol., 58(5):634?647, 2012. ? pages 106[RWNW] C. Robertson, M.A. Wulder, T.A. Nelson, and J.C. White.The effect of landscape pattern on mountain pine beetlespread. http://rose.geog.mcgill.ca/ski/system/files/fm/2009/Robertson_abstract.pdf, Downloaded on March 10,2013. ? pages 33, 34[SB03] N. Schtickzelle and M. Baguette. Behavioural to habitat patchboundaries responses restrict dispersal and generate emigration-patch area relationships in fragmented landscapes. Journal ofAnimal Ecology, 72(4):533?545, 2003. ? pages 32[SBTR99] L. Safranyik, H. Barclay, A. Thomson, and W.G. Riel. A pop-ulation dynamics model for the mountain pine beetle, Dendroc-tonus Ponderosae Hopk. (Coleoptera : Scolytidae). TechnicalReport BC-X-386, Natural Resources Canada, Canadian For-est Service, Pacific Forestry Centre, Victoria, British Columbia,1999. ? pages 2, 6120Bibliography[SC01] C.B. Shultz and E.E. Crone. Edge-mediated dispersal behaviorin a prairie butterfly. Ecology, 82(7):1879?1892, 2001. ? pages33, 67[SG89] A.J. Stock and R.A. Gorley. Observations on a trial of broadcastburning to control an infestation of the mountain pine beetleDendroctonus ponderosae. The Canadian Entomologist, pages521?523, June 1989. ? pages 69[SKT95] N. Shigesada, K. Kawasaki, and Y. Takeda. Modeling strati-fied diffusion in biological invasions. The American Naturalist,146(2):229?251, 1995. ? pages 30[SL54] D.E. Smith and M.L. Latham. The geometry of Rene Descarteswith a facsimile of the first edition. Dover Publications, NewYork, New York, 1954. ? pages 17[SL00] S. Samman and J. Logan. Assessment and response to barkbeetle outbreaks in the rock mountain area. Technical ReportRMRS-GTR-62, U.S. Department of Agriculture, Forest Ser-vice, Rocky Mountain Research Station, Ogden, UT, U.S., 2000.? pages 28[SLSH01] L. Safranyik, D.A. Linton, T.L. Shore, and B.C. Hawkes. The ef-fects of prescribed burning on mountain pine beetle in lodgepolepine. Technical Report BC-X-391, Natural Resources Canada,Pacific Forestry Centre, Victoria, British Columbia, 2001. ?pages 69[SRT] S. Strohm, M.L. Reid, and R.C. Tyson. Edge effects on moun-tain pine beetle movement. In preparation. ? pages 35[SS92] T.L. Shore and L. Safranyik. Susceptibility and risk rat-ing systems for the mountain pine beetle in lodgepole pinestands. Technical Report BC-X-336, Forestry Canada, Pacificand Yukon Region, Pacific Forestry Centre, Victoria, BritishColumbia, 1992. ? pages 2121Bibliography[ST09] S. Strohm and R.C. Tyson. The effect of habitat fragmentationon cyclic population dynamics: A numerical study. Bulletin ofMathematical Biology, 71:1323?1348, 2009. ? pages 32[ST12] S. Strohm and R.C. Tyson. The effect of habitat fragmentationon cyclic population dynamics: A reduction to ordinary differ-ential equations. Theoretical Ecology, 5(4):495?516, 2012. ?pages 32[SW04] R.L. Schooley and J.A. Wiens. Movements of cactus bugs: Patchtransfers, matrix resistance, and edge permeability. LandscapeEcology, 19:801?810, October 2004. ? pages 32[SW06] L. Safranyik and B. Wilson, editors. The Mountain Pine Bee-tle. Natural Resources Canada, Canadian Forest Service, PacificForestry Center, Victoria, BC, Canada, 2006. ? pages 1, 2, 6,8, 19, 28, 33, 37, 68, 104, 106, 110[TLM99] R. Tyson, S.R. Lubkin, and J.D. Murray. Model and analysisof chemotactic bacterial patterns in a liquid medium. Journalof Mathematical Biology, 38(4):359?375, 1999. ? pages 5[TR08] M.K. Trzcinski and M.L. Reid. Effect of management onthe spatial spread of mountain pine beetle (Dendroctonus pon-derosae) in Banff National Park. Forest Ecology and Manage-ment, 256:1418?1426, 2008. ? pages 68, 69, 78, 103[TSL00] R. Tyson, L.G. Stern, and R.J. LeVeque. Fractional step meth-ods applied to a chemotaxis model. Journal of MathematicalBiology, 41:455?475, 2000. ? pages 44, 45[Tur52] A.M. Turing. The chemical basis of morphogenesis. Phi. Trans.R. Soc. Lond. B, 237:37?72, 1952. ? pages 5[Tys96] R.C. Tyson. Pattern formation by E. coli - Mathematical andnumerical investigation of a biological phenomenon. PhD thesis,University of Washington, 1996. ? pages 15, 30, 44, 45122Bibliography[WP97] P. White and J. Powell. Phase transition from environmental todynamic determinism in mountain pine beetle attack. Bulletinof Mathematical Biology, 59:609?643, 1997. ? pages 2, 3, 6, 8,29, 35, 37, 53, 64, 66, 69, 109, 110[WSM04] J. Whittington, C.C. St. Clair, and G. Mercer. Path tortuos-ity and the permeability of roads and trails to wolf movement.Ecology and Society, 9(1), 2004. ? pages 33, 67123Appendix124Appendix AFFT Analysis and ParameterSensitivityA.1 FFT analysis of MPB dataAn example of the FFT analysis of the data (in 2007) is shown in FigureA.1. The detailed results of the FFT analysis of the data are shown inFigures A.2, A.3, and A.4. Each year the landscape was analyzed for regionsof incipient epidemic densities of MPB. The size and number of regionschosen in each year is displayed in Figure A.5.A.2 Parameter Sensitivity125A.2. Parameter Sensitivity0 5 10 15 20 250123456 x 10?4Distance (km)PowerDistance (km) Wavenumber (1/km)Figure A.1: The SNRA (top left) is plotted for 2007, where white areasdenote regions of MPB attack. The red highlighted region (left bottom)is analyzed using Discrete fast Fourier transforms to determine the powerspectral density (right). In this graph, the dominant wavenumber (?d) isrepresented by a solid red dot, and the upper (?u) and lower (?l) bounds onthe dominant wavenumber are displayed with vertical green and black lines,respectively.Table A.1: Table of parameter sensitivity for the non-dimensional model(2.4). The simple-index, S, computes the ratio of the standardized changein the wavenumber to the standardized change in the parameter values[Hae05]. S+ denotes a 10% increase in the parameter value while S? de-notes a 10% decrease in the parameter value. Negative values correspond tothe wavenumber increasing (decreasing) while the parameter is decreasing(increasing). Values of zero signify no change in the wavenumber with theparameter change.Par ? b1 ?q  kp ?p ?aS? 0.0423 0.0423 0.1273 0.2562 0 -0.33030.1273S+ 0 0 0.1022 0.2538 -0.0517-0.41930.1022126A.2. Parameter Sensitivity0 5 10 15 20 2500.20.40.60.811.21.4 x 10?40 5 10 15 20 25012x 10?40 5 10 15 20 250123456 x 10?30 5 10 15 20 250123456 x 10?30 5 10 15 20 2500.511.522.53 x 10?30 5 10 15 20 2500.511.522.53 x 10?30 5 10 15 20 2500.511.522.533.5 x 10?40 5 10 15 20 25012345678 x 10?40 5 10 15 20 2500.511.522.53 x 10?30 5 10 15 20 2500.20.40.60.811.21.41.61.8 x 10?30 5 10 15 20 2500.20.40.60.811.21.41.61.82 x 10?30 5 10 15 20 2500.511.522.53 x 10?3PowerWavenumber (1/km)199119931994199519961992Sum SumFigure A.2: Plots of power against spatial radial wavenumber, ?r, for years1991-1996 (progressing top to bottom, then left to right) for all data ana-lyzed. Different colours represent the analysis of a different region on thelandscape. The right figure in each year is the sum of the power spectraldensities over all regions in a given year.127A.2. Parameter Sensitivity0 5 10 15 20 2500.20.40.60.811.21.4 x 10?40 5 10 15 20 2501x 10?40 5 10 15 20 250123456 x 10?40 5 10 15 20 2500.20.40.60.811.2 x 10?30 5 10 15 20 2500.20.40.60.811.21.4 x 10?30 5 10 15 20 2500.511.522.533.544.55 x 10?30 5 10 15 20 25012345678 x 10?30 5 10 15 20 2500.0020.0040.0060.0080.010.0120 5 10 15 20 25012x 10?40 5 10 15 20 25012345678 x 10?40 5 10 15 20 2500.0010.0020.0030.0040.0050.0060.0070.0080.0090.010 5 10 15 20 2500.0020.0040.0060.0080.010.0120.0140.0160.0180.02PowerWavenumber (1/km)199719992000200120021998Sum SumFigure A.3: Plots of power against spatial radial wavenumber, ?r, for years1997-2002 (progressing top to bottom, then left to right) for all data ana-lyzed. Different colours represent the analysis of a different region on thelandscape. The right figure in each year is the sum of the power spectraldensities over all regions in a given year.128A.2. Parameter Sensitivity0 5 10 15 20 2500.20.40.60.811.21.4 x 10?30 5 10 15 20 2500.20.40.60.811.21.41.61.82 x 10?30 5 10 15 20 2500.0020.0040.0060.0080.010.0120.0140.0160.0180.020 5 10 15 20 2500.010.020.030.040.050.060.070.080 5 10 15 20 2500.0050.010.0150.020.0250 5 10 15 20 2500.0050.010.0150.020.0250.030.0350.040.0450 5 10 15 20 2500.20.40.60.811.21.41.61.82 x 10?30 5 10 15 20 250123456 x 10?30 5 10 15 20 2500.511.522.533.54 x 10?30 5 10 15 20 2500.0050.010.0150.020.0250.030 5 10 15 20 2500.511.522.533.544.55 x 10?30 5 10 15 20 2500.0050.010.0150.020.0250 5 10 15 20 2500.0020.0040.0060.0080.010.0120 5 10 15 20 2500.010.020.030.040.050.060.070.08Power2003200520082004Wavenumber (1/km)200620092007Sum SumFigure A.4: Plots of power against spatial radial wavenumber, ?r, for years2003-2009 (progressing top to bottom, then left to right) for all data ana-lyzed. Different colours represent the analysis of a different region on thelandscape. The right figure in each year is the sum of the power spectraldensities over all regions in a given year.129A.2. Parameter Sensitivity1990 1995 2000 2005 20100501001502002501990 1995 2000 2005 201005101520253035404550Yeara bNumber of Regions2Region Size (km )Figure A.5: The size (km2) (a) and number of regions (b) over the timeperiod 1991-2009. Multiple regions of the same size in the same year arerepresented by a single dot for clarity.130Appendix BBurn SimulationsWe simulated the burned 1-dimensional domains for three differentinitial conditions, beginning in the left undamaged patch (Figure B.1), be-ginning in the left damaged patch (Figure B.2), and beginning in the centralscorched patch (Figure B.3). The plots depict the nesting MPB density atthe end of every year over the spatial landscape.131Appendix B. Burn Simulations0 100 200 300 400 500024681012 x 1040 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 500024681012 x 1050 100 200 300 400 500024681012 x 1040 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 50000.511.52 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.5 x 1050 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.5 x 1050 100 200 300 400 50000.511.522.53 x 10516 source trees per ha19 source trees per ha22 source trees per ha25 source trees per haNormal AttractionNormal ReproductionNormal AttractionLow ReproductionHigh AttractionNormal ReproductionHigh AttractionLow ReproductionNesting MPB Density (MPB/ha)Year 1 Year 2 Year 3Space (m)Figure B.1: Nesting MPB population density over space for various initialdensity of source trees. The assumptions on the reproduction (low or nor-mal) and reproduction (normal or high attraction) change for each row offigures. This simulation assumes population of MPB emerge from the leftundamaged forest region.132Appendix B. Burn Simulations0 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 50000.511.522.533.5 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.5 x 10516 source trees per ha19 source trees per ha22 source trees per ha25 source trees per haNormal AttractionNormal ReproductionNormal AttractionLow ReproductionHigh AttractionNormal ReproductionHigh AttractionLow ReproductionNesting MPB Density (MPB/ha)Year 1 Year 2 Year 3Space (m)Figure B.2: Nesting MPB population density over space for various initialdensity of source trees. The assumptions on the reproduction (low or nor-mal) and reproduction (normal or high attraction) change for each row offigures. This simulation assumes population of MPB emerge from the leftpartially burned region.133Appendix B. Burn Simulations0 100 200 300 400 50000.511.522.5 x 1050 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.5 x 1050 100 200 300 400 50000.511.522.5 x 1050 100 200 300 400 50002468101214 x 1040 100 200 300 400 50000.511.522.5 x 1050 100 200 300 400 500051015 x 1040 100 200 300 400 50000.511.522.53 x 1050 100 200 300 400 50000.511.522.533.54 x 1050 100 200 300 400 500051015 x 1040 100 200 300 400 50000.511.522.5 x 1050 100 200 300 400 50000.511.522.533.5 x 10516 source trees per ha19 source trees per ha22 source trees per ha25 source trees per haNormal AttractionNormal ReproductionNormal AttractionLow ReproductionHigh AttractionNormal ReproductionHigh AttractionLow ReproductionNesting MPB Density (MPB/ha)Year 1 Year 2 Year 3Space (m)Figure B.3: Nesting MPB population density over space for various ini-tial density of source trees. The assumptions on the reproduction (low ornormal) and reproduction (normal or high attraction) change for each rowof figures. This simulation assumes population of MPB emerge from thecentral scorched forest region.134"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2013-11"@en ; edm:isShownAt "10.14288/1.0074062"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Mathematics"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@en ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Dispersal of Mountain Pine Beetle and impacts of management"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/44798"@en .