UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Avian influenza epidemic recurrence and approximate stochastic models Mata, May Anne Estrera 2017

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2017_may_mata_mayanne.pdf [ 3.38MB ]
JSON: 24-1.0347223.json
JSON-LD: 24-1.0347223-ld.json
RDF/XML (Pretty): 24-1.0347223-rdf.xml
RDF/JSON: 24-1.0347223-rdf.json
Turtle: 24-1.0347223-turtle.txt
N-Triples: 24-1.0347223-rdf-ntriples.txt
Original Record: 24-1.0347223-source.json
Full Text

Full Text

Avian Influenza Epidemic Recurrenceand Approximate Stochastic ModelsbyMay Anne Estrera MataB.Sc., University of the Philippines Mindanao, 2007M.Sc., University of Washington, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE COLLEGE OF GRADUATE STUDIES(Interdisciplinary Studies)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)April 2017c© May Anne Estrera Mata, 2017Thesis CommitteeThe undersigned certify that they have read, and recommend to theCollege of Graduate Studies for acceptance, a thesis entitled: Avian In-fluenza Epidemic Recurrence and Approximate Stochastic Mod-els submitted by May Anne Estrera Mata in partial fulfilment of therequirements of the degree of Doctor of PhilosophyRebecca C. Tyson, Unit 5 (Mathematics), The University of British ColumbiaSupervisor, Associate ProfessorPriscilla E. Greenwood, Department of Mathematics, The University of British ColumbiaCo-supervisor, Professor EmeritusRobert G. Lalonde, Unit 2 (Biology), The University of British ColumbiaSupervisory Committee Member, Associate ProfessorDavid B. Jack, Unit 3 (Chemistry), The University of British ColumbiaUniversity Examiner, Associate ProfessorSebastian J. Schreiber, Department of Evolution and Ecology, University of CaliforniaExternal Examiner, ProfessorApril 24, 2017(Date Submitted to Grad Studies)iiAbstractThis thesis is mainly concerned with avian flu epidemic recurrence, itscurrent paradigm, and further mathematical research. Generally, this thesisaims to characterise the recurrent pattern of epidemics simulated by stochas-tic avian flu models using mathematical techniques. Of particular interesthere are the stochastic fluctuations observed in recurrent epidemics. Thisthesis has two main parts.The first part presents a thorough analysis of a simple stochastic avianflu model to provide insight into the role of different transmission routes inits recurrent dynamics. Recent modelling work on avian influenza in wildbird population takes into account demographic stochasticity and highlightsthe importance of environmental transmission in determining the outbreakperiodicity, but only for a weak between-host transmission rates. A newanalytic approach is used here to determine the relative contribution ofenvironmental and direct transmission routes to the features of recurrentepidemics. Using an approximation method to describe noise-sustained os-cillations, the recurrent epidemics simulated by the stochastic model is iden-tified to be governed by the product of rotation and a slow-varying standardmean-reverting stochastic process, in a limiting sense. By analytically com-puting the intrinsic frequency and theoretical power spectral density, it canbe shown that the outbreak periodicity can be explained by both types oftransmission, and even by either one in the absence of the other.The final part of this thesis presents a novel approach to understand-ing the role of parametric (e.g. seasonal) forcing and stochasticity in thestochastic fluctuations around a cyclic solution. An approximate descrip-tion about these stochastic fluctuations is developed, which paves the wayfor a new mathematical tool to be used for analysing oscillations gener-ated from the interactions of non-linear terms and stochasticity. The theorydeveloped here is used to explore a stochastic avian flu model with season-ally forced environmental transmission which may be applicable to otherstochastic system with seasonal forcing.This thesis highlights the importance of approximation theory to analysecomplex stochastic systems such as avian flu epidemic recurrence.iiiPrefacei The author of this thesis owns certain copyright or related rights in it(the Copyright) and she has given The University of British ColumbiaOkanagan certain rights to use such Copyright, including for adminis-trative purposes.ii The main content of this thesis manuscript is an original, unpublished,and independent work of the author.iii Chapter 2 is written with the help of different review papers on avianinfluenza that the author of this thesis collected, synthesised, and sum-marised.iv Article versions of Chapter 3 and 4 are prepared for submission to ISI-indexed (or Scopus-indexed) journals.ivTable of ContentsThesis Committee . . . . . . . . . . . . . . . . . . . . . . . . . . iiAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Epidemic recurrence . . . . . . . . . . . . . . . . . . . . . . . 21.2 Modelling epidemic recurrence . . . . . . . . . . . . . . . . . 41.2.1 Damped oscillations . . . . . . . . . . . . . . . . . . . 41.2.2 Noise-sustained oscillations . . . . . . . . . . . . . . . 51.2.3 Parametric forcing: Forced limit cycle . . . . . . . . . 61.2.4 Parametric forcing with noise . . . . . . . . . . . . . . 71.3 Mathematical theories and methodology used in this thesis . 91.4 Thesis structure and contribution . . . . . . . . . . . . . . . . 10Chapter 2: The Epidemiology, Ecology, And Epidemic Re-currence Of Avian Influenza:An Overview . . . . . . . . . . . . . . . . . . . . . . 122.1 Avian influenza (AI) . . . . . . . . . . . . . . . . . . . . . . . 122.2 Avian influenza virus (AIV) . . . . . . . . . . . . . . . . . . . 12vTABLE OF CONTENTS2.3 AIV ecology in wild birds . . . . . . . . . . . . . . . . . . . . 132.4 Transmission of AIVs to and between birds . . . . . . . . . . 142.5 AI epidemic recurrence . . . . . . . . . . . . . . . . . . . . . . 152.5.1 LPAI outbreaks . . . . . . . . . . . . . . . . . . . . . . 152.5.2 HPAI outbreaks . . . . . . . . . . . . . . . . . . . . . 172.5.3 Modelling efforts . . . . . . . . . . . . . . . . . . . . . 182.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19Chapter 3: The Role Of Direct and Environmental Transmis-sion In Stochastic Avian Flu Epidemic Recurrence 213.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 The stochastic avian influenza model . . . . . . . . . . . . . . 233.2.1 Model description . . . . . . . . . . . . . . . . . . . . 233.2.2 Parameter values . . . . . . . . . . . . . . . . . . . . . 283.2.3 Preliminary analysis of the model . . . . . . . . . . . . 303.3 Analytic methods . . . . . . . . . . . . . . . . . . . . . . . . . 333.3.1 Approximate solutions . . . . . . . . . . . . . . . . . . 333.3.2 Power spectral density . . . . . . . . . . . . . . . . . . 343.3.3 The intrinsic frequency and the decay rate of the de-terministic dynamics . . . . . . . . . . . . . . . . . . . 343.3.4 Numerical tools and functions . . . . . . . . . . . . . . 353.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.4.1 An approximate avian flu epidemic process . . . . . . 353.4.2 The relative contribution of the different transmissionroutes . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Chapter 4: Sustained Oscillations in Stochastic ModelsWith Periodic Parametric Forcing . . . . . . . . . 494.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.2 A periodically parametrically forced (PPF) Markov model . . 524.2.1 Description . . . . . . . . . . . . . . . . . . . . . . . . 524.2.2 Stochastic fluctuations around the limit cycle solution 544.3 The approximation . . . . . . . . . . . . . . . . . . . . . . . . 564.3.1 Transformations . . . . . . . . . . . . . . . . . . . . . 564.3.2 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.4 Example: Driven harmonic oscillator with 1D noise . . . . . . 644.4.1 Power spectral density (PSD) . . . . . . . . . . . . . . 654.4.2 Amplitude and phase as stochastic processes . . . . . 674.5 Application: an avian flu model . . . . . . . . . . . . . . . . . 70viTABLE OF CONTENTS4.5.1 Stochastic avian flu model with seasonal forcing . . . 704.5.2 Deterministic dynamics . . . . . . . . . . . . . . . . . 724.5.3 Stochastic fluctuation around the limit cycle and itsapproximation . . . . . . . . . . . . . . . . . . . . . . 754.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Chapter 5: Conclusions and Future Work . . . . . . . . . . . . 815.1 Factors influencing epidemic recurrence . . . . . . . . . . . . 825.2 The stochastic fluctuations around a deterministic skeleton . 835.3 Avian flu in a stochastic and seasonally forced environment . 845.4 Open questions . . . . . . . . . . . . . . . . . . . . . . . . . . 86Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110Appendices A: Supplementary materials for Chapter 3 . . . . . . . 111A.1 Derivation of the avian flu SDE system . . . . . . . . . . . . . 111A.2 Stochastic linearization . . . . . . . . . . . . . . . . . . . . . . 113A.3 Approximate solution for linear diffusion equations in threedimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115A.4 Additional insight from the approximation on the interactionof disease components . . . . . . . . . . . . . . . . . . . . . . 118A.5 The subspace where the cycling takes place . . . . . . . . . . 121A.6 Derivation of the explicit form of the mean-field eigenvalues . 121Appendices B: Supplementary materials for Chapter 4 . . . . . . . 125B.1 Floquet Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 125B.2 Itoˆ formula for multi-dimensional process . . . . . . . . . . . 126viiList of TablesTable 3.1 Descriptions of the parameters in the avian flu modeland the values used in stochastic simulations. Therange of values for β in bold-faced is the set of valuesthat were used in [1]. . . . . . . . . . . . . . . . . . . . 29Table 4.1 Descriptions of the parameters in the avian flu modeland the values used in stochastic simulations of (4.78). 74viiiList of FiguresFigure 1.1 Observed recurrence in epidemics for infectious dis-eases:(a) Reported incidence of measles in Copen-hagen from 1920-1990 [2], (b) monthly cholera mor-tality from 1893-1940 in India [3, 4], and (c) annualprevalence of avian influenza from 1975-2000 in NorthAmerican wild birds (in green square) and shorebirds(red circle) [1]. Figures are adopted from cited sources. 3Figure 1.2 Whooping cough cases from actual data and SEIRmodels with seasonal forcing. The plots are adaptedfrom [5]. . . . . . . . . . . . . . . . . . . . . . . . . . 8Figure 3.1 A schematic diagram of the host-pathogen model foravian influenza adapted from Wang et al. [1]. Thehost population has three compartments correspond-ing to susceptible, infected, and removed individu-als. The solid arrows represent movement of individ-uals from one compartment to another as a result ofbirth, death, infection, and recovery processes. Thedashed arrows represent the increase or decrease inthe number of infected individuals due to the interac-tion (squiggly arrow) of the susceptible host with theavian influenza virus present in the environment. . . . 25ixLIST OF FIGURESFigure 3.2 The real (solid curves) and imaginary (dash-dot curves)parts of the eigenvalues of A0 in (3.10) associatedwith the stability of the endemic steady-state (3.8)plotted against R0 in the case where β = 0.05. Thereare three real eigenvalues when R0 < 1, of which twoare negative (thick solid lines). For R0 > 1, we havea complex-conjugate pair of eigenvalues with nega-tive real parts (thick solid lines) and another nega-tive eigenvalue (thin solid line). The imaginary partsof the complex eigenvalues are shown by the dottedlines. Other parameter values are shown in Table 3.1. 31Figure 3.3 Simulation of the stochastic model (3.2) and its cor-responding deterministic solution for β = 0.05 with(a) R0 = 1.2 and (b) R0 = 2.5, and for R0 = 2.5with (c) β = 0 (ρ = 0.4205) and (d) ρ = 0 (β = 14.5). 32Figure 3.4 Comparisons between the theoretical PSD of the ex-act process ξ(t) (solid line) satisfying (3.21) and theapproximate process ξapp(t)(dashed line) given by (3.25),for the fluctuations of the susceptible, infectious, andthe virus populations. Default parameter values arein Table 3.1 with β = 0.05 and ρ = 0.4. . . . . . . . . 37Figure 3.5 A stochastic realization of the population fluctuationsby simulating (3.21) (solid line) and (3.25) (dashedline) and their corresponding stationary standard de-viations, i.e. typical amplitudes (in horizontal lines)computed using (3.28). . . . . . . . . . . . . . . . . . 38Figure 3.6 (a,b) Plot of the ratio λ/ω as a function of β andρ. The triangular white region in the lower left iswhere no recurrence is observed. Panel (a) considers0 ≤ β ≤ 50 and 0 ≤ ρ ≤ 1 while Panel (b) consid-ers 0 ≤ β ≤ 300 and 0 ≤ ρ ≤ 3. The grey regionin both panels is where λ/ω ≤ 0.35. In Panel (b),the black curve corresponds to λ/ω = 1. (c) Plot ofthe basic reproduction number R0 as a function of βand ρ. The darker (blue online) shade indicates thatno epidemic can be observed in the model for thisparameter region. Parameter values and ranges fromthe literature are in Table 3.1. . . . . . . . . . . . . . 40xLIST OF FIGURESFigure 3.7 Plots of (a) the intrinsic frequency ω and (b) thesteady-state proportions of the susceptible and infec-tious populations as functions of the direct transmis-sion rate β (left) and the environmental transmissionrate ρ (right). The dark region in (a) is where the 2to 8 year recurrence period is observed. Parametervalues are in Table 3.1. . . . . . . . . . . . . . . . . . 42Figure 3.8 Top panels: Theoretical PSDs of the infectious popu-lation fluctuations ξI(t) for (a) ρ = 0.22, 1.1, 2.2 whenβ = 0 and (b) β = 7.5, 37.5, 75 when ρ = 0. Thelinewidth of the PSD curves increases with the trans-mission parameter values. Bottom panels: Approxi-mate dominant outbreak period, 2pi/ω with formula(A.42), as a function of (c) ρ when β = 0 and (d)β when ρ = 0. The square markers in (c) and (d)are located at ρ = 0.22, 1.1, 2.2 and β = 7.5, 37.5, 7.5,respectively. The black dots represent the exact dom-inant outbreak period obtained using the theoreticalPSDs in the top panels. The black dots and squaremarkers are indistinguishable from each other. Thehorizontal lines in (c) and (d) indicate the 2 to 8 yearperiods observed in actual prevalence data (See [1]).The corresponding R0 values are also shown. Defaultvalues for all other parameters are given in Table 3.1. 44Figure 3.9 Typical amplitude (intensity) of the fluctuations inthe proportion of (a) susceptible, (b) infectious, and(c) virus populations, measured by their stationarystandard deviations (3.29) as a function of the envi-ronmental transmission rate ρ and β. All other pa-rameters are at the default values in Table 3.1. . . . 46Figure 4.1 Numerical PSDs of the approximate (solid curves)and exact (dashed curves) solutions of (4.62), given(a,b) damped and (c,d) undamped regimes as definedin Section 4.4.1, obtained using forcing amplitudes(a,b) b = 0.01 and (c,d) b = 0.5 for noise levels σ0 = 1(gray curves) and σ0 = 5 (black curves). Other pa-rameters are held fixed, namely a = 1, b0 = 2,k0 = 0.1for underdamped case, and k0 = 0 for the undampedcase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66xiLIST OF FIGURESFigure 4.2 Deterministic dynamics when (a) b = 0 and (b) b =1. The initial condition is indicated by a black ar-row. The initial conditions used are: (φ1, φ2, ψ) =(0.3295, 0.0298, 0.6935). . . . . . . . . . . . . . . . . . 73Figure 4.3 Limit cycle generated from the solution of (4.82) fordifferent values of the forcing amplitude b. . . . . . . 73Figure 4.4 Computed values of (a) A11(t) and (b) C11(t) in (4.83)-(4.85). The parameter values used are in Table 4.1.Each function has a periodic behaviour with periodequal to 1, the period of the limit cycle and the forc-ing frequency. . . . . . . . . . . . . . . . . . . . . . . 76Figure A.1 (Colour online) A sample path of the approximatefluctuations given by (3.25) when the first term is (a)not set to zero and (b) set to zero. The grey regionis the plane (A.31) that lies in the subspace spannedby the eigenvectors of −0.3091± 0.8377i. . . . . . . . 119Figure A.2 (Colour online) Plot of ζ (left panel) and λ (rightpanel) as functions of β and ρ. The white region iswhere R0 < 1, i.e. noise-sustained oscillations cannotbe observed here. Default parameter values are inTable 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . 120xiiAcknowledgementsThere are so many people to acknowledge for their moral, emotional,and academic support in my doctoral study. Since they are too many tomention, I hereby only recognise the most important contributors of thisgreat work. I wish to thank the following people whom I have close contactwith during my PhD study here in Canada.First of all, I would like to express my sincerest gratitude to my wonder-ful supervisors, Rebecca Tyson and Cindy Greenwood. The positive energyof my academic mentor Rebecca is contagious and I will forever treasure hermotherly care for me in this PhD journey. All pieces of advise and encour-agements she shared are valuable to me and her brilliant ideas especiallywith regards to writing this thesis were inspiring. On the other hand, thisthesis was conceived from a small summer project that was suggested byCindy before I left UBC Vancouver. I am very grateful that Cindy chal-lenged my thoughts and trained me to become a probabilist. Her mentoringstyle taught me to be stronger than I used to be and I will always rememberher words of wisdom. I also would like to thank Leah Keshet for choosingme to be part of her research group when I first applied for PhD at UBCVancouver. Leah helped me see my research potential and the support shehas given me ushered the fulfilment of this dream. Thanks to Bob Lalondefor the insightful conversation and constructive remarks on the epidemiologyand ecology of avian influenza part of my thesis. I thank my thesis exam-iners, David Jack and Sebastian Schreiber, for their time spent to evaluatemy research work.I thank my academic brothers and sisters in the TyLab group: Jessa,Joe, Jimit, Sarah, Maria, Geoff, Charles, and Michael, for the stimulatingdiscussions, questions, and anything-under-the-sun chitchats. To Joey andAli for the good laughs, the late night rides, and inspiring conversationson research writing. To my colleagues in the Computational Ecology Re-search Group (CERG), especially to Bruno, Fabian, Matt, and Rodolphe,for the support and constructive criticisms. More than helping me with pa-per works, providing supplies, and scheduling my teaching duties, I thankthe Unit5 Staff, Hilda and Ashlee, for their moral support.xiiiAcknowledgementsTo my folks at my home university, the University of the PhilippinesMindanao (UPMin), I am very grateful for their continued support and forbelieving me in my pursuit to become a biomathematician. I thank theUPMin officials for granting me this special privilege of leaving the countryto study and getting paid at the same time. I am grateful that they held onto my intention of returning service after my study.To those people who showed me that there is more than a scientificlife, I thank them all a ton for their support especially during the lowestpoint of my study. To my very supportive and amazing friends – Ellaine,Aleza, Stephanie, Aaron, and Brent – for welcoming me like family. To myspiritual mentors, Tita Liza, Tito Odi, Ate Arleen, Kuya Daniel – who hadhelped me walk in faith. To my Christian brothers and sisters all over theworld, especially the Buhangin Community Church, the Willowpark Churchand all the small groups I have attended with, for praying with me. ToSimon, Jenny, Tian, and Isabell for the life-giving stories, encouragementsand exhortations. To my ever-dearest friend Jen M, for the sisterly piecesof advice. To my former roommates Janine, Jocelyn, and Chalise, for thefun adventures and “Canadianizing” me. To Mr. and Mrs. Diopita for theirgenerosity and providing me a comfortable place to stay for my last fewmonths here in Canada.Last but not least, thanks to my biological family especially to mymother, Mimay Luz, who stayed strong for me despite all the crises wefaced and to my only sibling, Janelle, who cheered me up and helped melearn not to take life so seriously. To Kuya Telo and other online friendsfor the good company and for boosting my confidence. And of course, toPatrick, for his love that continually reminds me of my good qualities andinspires me to finish this degree.Finally, I would like to thank the heroes from the Good Book who taughtme these seven valuable lessons during my PhD study:To fight the good fight of the faith, be a good steward of my intellectualgifts, conquer my fears, draw strength from the real Source, be wise as aserpent, be humble as a dove, and believe that all things work together forgood for those who are called according to the LORD’s purpose.Hakavod lashem!xivDedicationDespite all stochastic perturbations,His steadfast love remains as my deterministic skeleton.Still in awe of His faithfulness!Truly, He never left me behind.Hence I offer this great workfor the glory of the Almighty, not mine.Beloved Mayang (2017)xvChapter 1IntroductionThe ever-present outbreaks of avian influenza and its pandemic threatsemphasize the importance of understanding and predicting how the diseasespreads. Mathematical models are important tools in the study of diseasedynamics. Theoretical epidemiological models and methods may not alwaysbe directly applicable to real-life epidemics but they provide qualitative in-sight into the dynamics of the spread of infection. Most of these modelsare mathematically intractable or difficult to analyse and so numerical sim-ulations of these models have been used for predictions. Hence, analysis ofuseful disease models remains as a challenge [6]. Moreover, real-life applica-tions of disease models demand for more advanced mathematical methodsto be developed that can also be useful in other fields of study [7, 8].The methods we developed here can be applied to study infectious dis-eases, such as measles, that display a recurrent pattern. All of the models weconsider in this thesis are formulated as continuous-time Markov processeswhose temporal evolution is described by a system of stochastic differentialequations (SDEs). In the physical literature, the SDE system is typicallyderived through the master equation and system-size expansion [9]. Here,however, we use an alternative approach known as strong approximationtheorems for Markov chains [10] to obtain the full stochastic dynamics. Aswith the system-size expansion, the full stochastic dynamics are decomposedinto a macroscopic part (deterministic dynamics) plus a stochastic fluctuat-ing part scaled by the square-root of the population size.In this thesis, we are concerned with applying and developing mathemat-ical methods to analyse disease models that simulate stochastic populationfluctuations. Our study focusses primarily on understanding recurrent avianflu epidemics. While stochastic simulations can capture the qualitative epi-demic pattern of recurrent diseases, there remains considerable discussionregarding the characteristics of this pattern. The novel aspect of this thesisis the development of an approximate description for the recurrent dynamicsof avian influenza and the subsequent analysis of potential mechanisms thatdrive disease recurrence.11.1. Epidemic recurrenceIn this chapter, we introduce the technical terms and theoretical back-ground which are necessary to support the content of this thesis. Specifically,we describe epidemic recurrence and its models. Moreover, we briefly dis-cuss here the mathematical theories used in this study. Finally, we note indetail the structure of this thesis, its aims, and contributions.1.1 Epidemic recurrenceA disease is said to display recurrent epidemics if its incidence (number ofcases per unit time) fluctuates dramatically over a period of time. Recurrentepidemics have been observed in childhood diseases such as measles, whoop-ing cough, and chicken pox. It is also observed in non-childhood diseasessuch as avian influenza and cholera (see Figure 1.1).A recurrent disease can exhibit various patterns of recurrence in differentlocations and periods of time [11]. Measles, for instance, displays recurrentoutbreaks that have been observed to transition from an annual to a bien-nial pattern in some locations [2], or exhibit an irregular pattern in otherlocations [12]. Patterns of recurrent epidemics also vary from one disease toanother. In contrast with measles data from Copenhagen [2], wavelet anal-ysis of data for avian influenza among wild birds revealed multi-year majoroutbreak periodicity [1, 13].The pattern of recurrent epidemics is characterized by large-amplitudefluctuations, which are said to arise from a resonance phenomenon. Under-standing the underlying mechanism of this ubiquitous feature of recurrentdisease has received attention from mathematical modellers with an aimto provide additional insights into the disease transmission process. Whengiven an epidemiological or ecological time-series with a recurrent pattern,typical questions raised are:1. What are the periodicities?2. Can we identify the sources of these periodicities?3. Why does the amplitude size of the epidemics vary?These questions have been addressed in the context of childhood diseasesusing traditional modelling approaches. One useful tool for distinguishingdifferent periodicities in a time-series is the power spectrum, which showsthe distribution of frequencies present in the time-series data [14, 15]. Thepower spectrum has been used to demonstrate that childhood diseases, e.g.measles and whooping cough, have a strong seasonal component. This tool21.1. Epidemic recurrence(a)(b) (c)Figure 1.1: Observed recurrence in epidemics for infectious diseases:(a) Re-ported incidence of measles in Copenhagen from 1920-1990 [2], (b) monthlycholera mortality from 1893-1940 in India [3, 4], and (c) annual prevalenceof avian influenza from 1975-2000 in North American wild birds (in greensquare) and shorebirds (red circle) [1]. Figures are adopted from citedsources.31.2. Modelling epidemic recurrencehas also been used to illustrate the influence of one disease parameter on theperiodicity of disease outbreaks. In this thesis, we use the power spectrumto study the relationship between the approximate and exact stochasticprocesses describing epidemic recurrence.1.2 Modelling epidemic recurrenceModels in the form of differential equations are widely used in the studyof biological dynamics [16–19]. Epidemic recurrence, in particular, is an-other area in which mathematical modelling has played important role. Itis noteworthy to mention that in modelling epidemic recurrence, extensivedatasets are often available and can be used to validate and inform themodelling exercise. Here we discuss the existing paradigms for modellingepidemic recurrence.1.2.1 Damped oscillationsThe model developed by Kermack and McKendrick [20] set the foun-dation for epidemic modelling studies. Their model describes the spreadof an infectious disease wherein recovered individuals gain immunity to re-infection. Fundamentally, the model considers a population that can bedivided into three classes or compartments and keeps track of how the pop-ulation of each class changes over time. The dynamic variables are the num-ber of susceptible (S(t)), infectious (I(t)), and recovered (R(t)) individualsat time t. The susceptibles are those individuals who have not been infectedby the disease while infectious individuals are those who have the diseaseand the ability to spread it to susceptibles. Recovered individuals, on theother hand, are those who no longer have the disease and are immune tothe disease after infection. The SIR (susceptible-infected-recovered) modelassumes a (homogeneous) well-mixed population and is formulated in termsof three ordinary differential equations. It was learned from the classicalSIR model that an epidemic threshold exists, which is defined in terms ofthe basic reproductive number R0, i.e. the average number of secondaryinfections caused by one infected individual in a pool of susceptibles [20]. Adisease can spread and is said to be endemic in the population when R0 > 1,otherwise, it dies out.A modification of the SIR model serves as a starting point for modellingepidemic recurrence. An SIR model that gives rise to damped oscillatorybehaviour, i.e. decaying oscillations, is a typical starting point for recurrentbehaviour. These models has been used to investigate measles dynamics41.2. Modelling epidemic recurrence[21]. The key point here is that the disease model exhibits an endemicequilibrium wherein infectives tend to oscillate around a mean value overtime, i.e. exhibit damped oscillations, as a new supply of susceptibles enterinto the system. Damped oscillations can also be observed when anotherpopulation class is incorporated to the SIR model. For example, one canadd an exposed class, a population of individuals that are exposed to thedisease but not yet infectious. Such a model is known as an SEIR model,and is biologically meaningful for many childhood diseases [21]. Anotherexample is a simple host-pathogen model such as the SIR-V model, where Vrepresents the virus population. This model is used to study avian influenza(see [1] and Chapter 3).Recurrent behaviour, however, is not entirely captured in models pro-ducing damped oscillations because, in the long-run, the infectives tend to aconstant level. Sustained oscillations obtained by introduction of age struc-ture in the population, delays and non-linear infection terms [22–25], onthe other hand, are not able to predict the regularity in patterns of mostepidemics. From the literature, it is apparent that there must be two mainingredients necessary to model epidemic recurrence accurately: (1) stochas-ticity, because the population is made up of individuals and it changes bymeans of random processes [26]; and (2) parametric forcing, because someparameters vary periodically over time, e.g. seasonally forced transmis-sion rate for measles [3, 12, 27]. These two factors are well understoodindependently but there is a considerable question as to how they interact[5, 11, 28, 29].1.2.2 Noise-sustained oscillationsInteracting populations, in general, behave as stochastic processes. Thediscrete nature of populations means that there is inherent randomnessin the birth and death terms, thus introducing demographic stochasticity.In modelling epidemics, the role of stochasticity has long been recognisedas a key element in the appearance of a recurrent pattern of epidemics.Bartlett [30], in his first stochastic simulations, recognised that the addi-tion of stochastic noise to a simple deterministic model can drive recurrentdynamics. The idea here is that noise maintains oscillations by constantlyperturbing the system away from its steady state [31]. A number of sys-tems exhibit oscillations sustained by noise (or noise-sustained oscillations)[26, 32–34]. One example is the amplified biochemical oscillations studiedby McKane et al. [34] where they show that noise alone can give rise to dra-matic fluctuations and explain that intracellular processes may have taken51.2. Modelling epidemic recurrenceadvantage of inherent stochasticity to amplify the oscillations. From a theo-retical standpoint, noise-sustained oscillations are generated by a stochasticmodel whose deterministic version shows damped oscillatory behaviour andare quantified as coherence resonance [35].According to Coulson et al. [5], the role of stochasticity can either becategorised as active or passive. Noise has an active role when the systemproduces new patterns, due to the interaction of non-linear dynamics andstochasticity, that the deterministic model alone cannot capture [28, 29].Noise has a passive role when perturbation of the ’deterministic skeleton’is sufficient to explain the stochastic dynamics [36]. In the context of dis-ease dynamics, the passive noise interpretation has gained weight becauseit seems to explain the behaviour of time-series and simulations for a broadrange of diseases [5]. In fact, for systems exhibiting noise-sustained oscilla-tions, demographic stochasticity generally has a passive role.In epidemic modelling, the use of noise-sustained oscillations to describedisease recurrence has been observed. For example, Wang et al. [1] addedstochastic noise to an SIR-V model for avian influenza and showed thatthe dominant outbreak period varies with the environmental transmissionrate, i.e. the rate at which the virus enters the host from the environment.Modelling epidemic recurrence with this technique can be challenging, es-pecially when one is interested in identifying the parameters that influencethe periodicity and/or intensity of outbreaks in a system with more thantwo population classes. One of the aims of this thesis is to demonstrate theuse of an approximate stochastic process to study recurrent behaviour in apopulation system composed of four populations.1.2.3 Parametric forcing: Forced limit cycleUp to this point, we have discussed models of epidemic recurrence with-out external parametric forcing. Here we introduce the forced limit cycle asanother paradigm for a recurrent pattern in epidemics. A common source ofparametric forcing is time-dependence in a model parameter. In this thesis,we define a forced limit cycle as a stable limit cycle generated by a determin-istic system with time-periodic parameter. In the context of epidemiologicalmodels, forced limit cycles occur in seasonally forced systems [33, 37, 38].The most popular approach for introducing seasonal forcing in an epi-demic model is by making the transmission rate a periodic function of time.For example, in studying measles recurrence, the transmission rate is sim-ply represented as a step function so that transmission is high during schoolterms and low otherwise [22, 28]. This deterministic forced epidemic model61.2. Modelling epidemic recurrencecan display a wide range of dynamic behaviour such as limit cycles and at-tractors with sub-harmonic frequencies depending on the magnitude of theseasonal forcing [39–41].The seasonality of avian influenza outbreaks points to a seasonally vary-ing transmission rate as a potential driver. This idea has been investigatedby [42] who showed that a simple epidemic model with periodic transmissionrate and disease-induced mortality best explained the periodicity of H5N1avian influenza data from Food and Agriculture Organization of the UnitedNations [43]. As with the measles model mentioned above, this simple forcedmodel for avian flu exhibits complex dynamics, such as period-doubling andchaos, for various forcing magnitudes.With respect to the fluctuating avian flu epidemics seen in wild birds,this thesis explores the possibility of modelling the recurrent pattern as aforced limit cycle. One potential driver for recurrent behaviour that weinvestigate in this study is seasonality in the environmental transmissionrate.1.2.4 Parametric forcing with noiseIt is likely that both seasonal forcing and stochasticity are importantin capturing realistic epidemic dynamics. To illustrate this view, considerthe different plots of whooping cough cases in Figure 1.2. Rohani et al.[29] demonstrated that although the deterministic SEIR model with sea-sonal forcing provides the appropriate ’deterministic skeleton’ for the actualwhooping cough cases reported in Birmingham, adding stochastic noise im-proves the model’s ability to capture the observed qualitative dynamics.Despite all known studies involving discussion of the role of stochasticityand seasonal forcing [37, 44], little is known about the stochastic dynamics ofmodels exhibiting both limit cycle and stochasticity. This point motivatesone of the main questions this thesis aims to answer: how do seasonalityand stochasticity interact within the fluctuations of the system when thedeterministic skeleton of the system is a robust oscillation, i.e. stable limitcycle? Can we mathematically distinguish between the fluctuations aroundan endemic fixed point from the fluctuations around a cyclic solution? Thelatter part of this thesis is devoted to answering these questions by meansof analytic methods for stochastic equations.71.2. Modelling epidemic recurrenceFigure 1.2: Whooping cough cases from actual data and SEIR models withseasonal forcing. The plots are adapted from [5].81.3. Mathematical theories and methodology used in this thesis1.3 Mathematical theories and methodology usedin this thesisBiological systems are characterised by changes in natural variables andso are often described in terms of differential equations. These equationsare also referred to as dynamical systems. Given a dynamical system, one isinterested in obtaining its solution to accurately determine the ’fate’ of thesystem over a period of time. Most differential equations with real-life appli-cations have solutions that are difficult or even impossible to obtain analyt-ically. Qualitative methodologies for studying dynamical systems, such asstability and bifurcation analyses [45], are straightforward to apply in mostdeterministic equations. In the case when random processes are involved,the behaviour of the system is predicted by means of probability distribu-tion and generating functions [46]. This is due to the fact that multiple runsof a stochastic model generate different realisations that show variations.Furthermore, when both randomness and external forcing are incorporatedin the system, the analysis of its dynamical behaviour demands tractablemathematical techniques.In this thesis, we focus our attention on stochastic differential equation(SDE) models to describe the epidemic recurrence observed in avian fludynamics. These stochastic epidemic models are approximations to theMarkov chain model for population growth in the epidemiological system.The SDE models are developed based on the theorem of Kurtz [10], whichshows that the (density-dependent) Markov jump process, normalized bythe total population size N , and the approximating diffusion process canbe constructed on the same probability space. Kurtz [10] proved that thediffusion process deviates from the jump process by O (logN/N) as N →∞. Application of this method to epidemic modelling is found mostly inthe literature on stochastic processes [47–49]. Note that in the physicalliterature, the method for constructing SDE models is based on the system-size expansion [9]. This thesis, however, uses the former approach (i.e.Kurtz’s method) in formulating the SDEs for avian flu epidemics, not onlybecause it is straightforward to implement but to also introduce this methodto the present body of knowledge on stochastic modelling of epidemics.The recurrent dynamics displayed by the SDE models in this thesis areanalysed by studying the equations mathematically. Here we apply thetheorem of Baxendale and Greenwood [50] (the BG theorem) to study thestochastic system without external forcing by constructing an approximateprocess that allows us characterise the simulated avian flu recurrent epi-91.4. Thesis structure and contributiondemics. The power spectral density (PSD) [15] is used to check the proximityof the approximate process to the actual process, and stochastic simulationsvia the Euler-Maruyama scheme [51] are used to generate sample paths.In investigating the stochatic system with external forcing, i.e. seasonalforcing, Floquet theory [52] is applied, which provides the stability informa-tion of the deterministic perturbation around a cyclic solution (e.g. limitcycle). The use of this theory is common in the physical literature that takesinto consideration linear differential equation with time-periodic coefficients.Another novel aspect of this thesis is the extension of the BG theorem todevelop an approximate process for fluctuations in the time course of epi-demics resulting from the interaction of nonlinear dynamics and stochas-ticity. Moreover, we extend BG theorem to noise around the forced limitcycle. Finally, fundamental probability concepts and methodologies, e.g.Itoˆ’s formula [53], Ornstein-Uhlenbeck process [46], etc., are essential forthe development of the theory.1.4 Thesis structure and contributionIn this thesis, we study the recurrent pattern of avian flu epidemics bymeans of analytic techniques. We begin by gathering information about theepidemiology, ecology, and outbreak recurrence of avian influenza (Chapter2).In Chapter 3, we first develop an understanding about the current paradigmfor modelling avian flu recurrent epidemics. This is done by analysing math-ematically, i.e. using approximation technique, a host-pathogen stochasticmodel for avian influenza recurrent epidemics with two types of transmissionmodes. It is found with the aid of analysis that the periodicity of recurrenceis determined by the transmission modes while the noise has a substantialcontribution to the intensity of epidemics.In Chapter 4, we develop a new approximation tool that can be usedto understand epidemic recurrence modelled as stochastic perturbationsaround a limit cycle. Our tool is an extension of the approximate modelfor stochastic fluctuations around an equilibrium, i.e. the result of Baxen-dale and Greenwood [50]. Our approach can be applied to a general classof SDE systems - SDEs with time-dependent drift and diffusion coefficientsdue to periodic parametric forcing. Such systems naturally occur when oneof the parameters is a periodic function of time.Finally, the conclusion of this thesis is presented in Chapter 6. The mainfinding points to the deterministic mechanism playing a key role in influ-101.4. Thesis structure and contributionencing avian flu epidemic recurrence described by stochastic perturbationsaround the deterministic dynamics (e.g. endemic equilibrium or forced limitcycle solution).11Chapter 2The Epidemiology, Ecology,And Epidemic Recurrence OfAvian Influenza:An Overview2.1 Avian influenza (AI)Avian influenza (AI) also known as “bird flu” is a viral disease that com-monly infects poultry and wild birds. In poultry, the disease is characterisedby various symptoms such as nasal discharge, lack of energy, and diarrhea,resulting to “fowl plague” in which mortality of at least 75% has been ob-served [54]. In wild birds, the disease do not usually manifest symptoms [55].More than a decade ago, AI was a disease of limited significance [56, 57] butthis perspective has changed since the emergence of various virus strains(e.g. H1N1 in 1976, 1986, and 1988; H3N2 in 1993; H7N7 in 1996; H9N2 in1998, 1999; H5N1 in 1997, 2003, and 2013; H7N9 in 2014) that are deadlyto humans [58–60]. Since the disease poses threat to human health, severalefforts have been made to control its spread such as culling large number ofpoultry and restrictions on trade of poultry products [61].2.2 Avian influenza virus (AIV)The causative agent of avian influenza is an influenza A virus [62] consist-ing of different combinations of haemaggluttinin (HA) and neuraminidase(NA) surface glycoproteins [63, 64]. A total of 16 HA subtypes (H1-H16)and nine NA subtypes (N1-N9) have been recognised [65], and each virushas one HA and one NA antigen.The avian influenza viruses (AIVs) can be divided into two distinctgroups according to their virulence in poultry, namely the low pathogenicavian influenza virus (LPAIV) and the highly pathogenic avian influenza122.3. AIV ecology in wild birdsvirus (HPAIV). The LPAIVs are commonly isolated from many species ofwild birds [66] and cause little harm in the wild [55, 65]. On the other hand,HPAIVs are not maintained by wild birds [67] but have been isolated fromwild population in the event when the disease has spread to domesticatedbirds [68, 69]. Significantly, certain LPAIV subtypes (e.g. H5 and H7) maybecome highly pathogenic if they are transmitted to poultry [70–72]. Theability of LPAIVs to mutate into HPAIVs [73] in poultry and the diversityof these viruses circulating in bird populations in the wild potentially makeswild birds the primary source of introduction of influenza into humans [74–76].2.3 AIV ecology in wild birdsThe AIVs found in wild birds constitute the historic source of humaninfluenza viruses having a rich pool of diversity that may likely lead to cross-species transmission [55, 77]. A most recent example is the transmissionof H7N9 AIV to humans in China (see [78] and references therein). TheH7N9 virus is associated with a significant risk of mortality to humans. Inthe year 2014, 212 deaths were reported to the World Health Organizationdue to infection with H7N9 alone. Thus, it is important to have a betterunderstanding of the ecology of AIVs in wild birds.Influenza viruses infect a great variety of birds [66, 79, 80]. Various sub-types of AIVs are maintained in and isolated from wild bird populations froma variety of major Families. In 1961, an influenza virus was first isolatedfrom common terns (Sterna hirundo) in South Africa [81]. Further moni-toring suggested that ducks and geese are important natural reservoirs ofAIV [82]. Sharp et al. [83] however suggested that waterfowl do not act as areservoir for all AIVs. It is possible that AIV that ultimately gets transmit-ted to humans is maintained in shorebirds and gulls, where the AIV isolatespredominantly have different subtypes than those isolated from ducks [84].According to Stallknecht and Shane [82], AIVs infect 90 bird species from22 different Families and 12 Orders. On the other hand, Olsen et al. [85]lists at least 105 wild bird species from 26 families. With the large numberof recent AI outbreaks, it seems that the actual number of fully susceptiblespecies is likely to be much greater.The complexity of the ecology of AIV is linked to the biology of thenatural host and environment. For example, it seems incidental that ducksliving in freshwater habitats have relatively high prevalence rates of AIVinfection [86] in contrast to shorebirds and gulls, which are salt water birds.132.4. Transmission of AIVs to and between birdsHowever, as demonstrated by Stallknecht et al. [87] some AIV isolates aremore sensitive to salinity. Moreover, persistence and prevalence of AIV asoutlined by Stallknecht et al. [88] can depend on three factors: (i) concen-tration of AIV shed at an adequate time duration, (ii) the stability of AIVin the environment, and (iii) the AIV concentration required for productiveinfection of the next host (infective dose).Human intervention could play a role in the ecology of AIV. For exam-ple, the movement of infectious agents among numerous domestic avian andmammalian species due to global wildlife trade provides transmission mech-anism into wild bird population when these infected species are released intothe wild [89]. However, as reviewed by Reed et al. [90], human effects onbird movement patterns may have greater significance. For example, climatealterations around urban areas may distort the disease epidemiology as it af-fects the pattern of migration, population densities, and species interactionof many wild bird species.2.4 Transmission of AIVs to and between birdsA potential cause of initial infection for AI outbreaks in poultry couldbe direct or indirect contact with infected waterfowl populations [91–93].However, the knowledge about the transmission mechanism of AIVs betweenbirds is insufficient. It is suggested that AIVs could be transmitted amongbirds via oral-oral (respiratory) route [94] or faecal-oral/faecal-cloacal route[95].In oral-oral route, i.e. direct transmission, AIVs are transmitted throughdirect contact between infected and susceptible birds. Past experimentalassessments on the transmissibility of AIVs in poultry suggested that directtransmission is a complex process as it depends on the virus strain, the birdspecies and many environmental factors [91, 96–99]. For example, HPAIVstend to show much poorer transmission from infected to susceptible hostsin poultry than LPAIVs in both natural and experimental settings. Onthe other hand, low LPAIV samples from mallards may suggest that therespiratory tract plays a limited role in the replication, transmission andecology of LPAIV in dabbling ducks [100], yet oral-oral route may be relevantfor bird species in which faecal-oral transmission would prove difficult.The fecal-oral route is identified to be the main transmission mechanismamong aquatic birds. LPAIVs are thought to be transmitted via the fecal-oral route in wild bird population when infected birds shed virions in theirfaeces and susceptible birds ingest the virions from the contaminated aquatic142.5. AI epidemic recurrencehabitat [55, 87, 101–103]. The filtering or drinking of contaminated watercould be a potential explanation for the high AIV prevalence among birdspecies aggregating in wetlands. Moreover, the prolonged infectious periodsof the virus could potentially allow temporal and spatial connectivity amongbird populations with their respective virus populations. In freshwater lakes,AIVs remain infectious for up to 4 days at 22 degrees Celsius and up to 6months at the freezing point [93]. The virus stays infectious for even longertime durations in ice or frozen ground [104–107]. During the breeding season,infected birds in the Arctic region shed AI virions into the environmentthrough faeces, which then persist in cold water, in ice or frozen groundthroughout the winter. The thawed contaminated ground ice or frozen lakebecomes, in the next spring, a source for re-infection among birds returningto the breeding ground [93, 108] suggesting that bird migration is associatedwith the persistence of AIV in wild bird populations.Virus transmissibility may be related to the amount of virus releasedorally or from excretions. Several factors such as virus strain, bird species,and immune response can affect the virus shedding rate [71]. For example,it is possible that a relatively small quantity of virus is excreted duringthe course of HPAIVs infections, as they cause rapid deaths in poultry.The Asian H5N1 HPAIV, on the other hand, are shed at much greaterconcentrations and for longer durations from the respiratory tract than bythe cloacal route [109, 110]. Such differences may potentially explain theefficiency of virus transmission within bird population or to domestic poultry[88].2.5 AI epidemic recurrenceSurvey-based reports suggest that AI epidemics are recurrent and aretied to the persistence and prevalence of AIV. Here we give an overview ofLPAI and HPAI outbreaks that have been observed in the past as reviewedby several authors [66, 111–113]. We also summarize and highlight theexisting modelling efforts that have been done to understand AI outbreaks,e.g. the role of transmission on its multi-year periodicity.2.5.1 LPAI outbreaksThe prevalence of LPAIV in its natural reservoir, i.e. wild birds, isbelieved to be the source of LPAI infection in poultry. For certain subtypesof LPAIV, when the AI infection moves to poultry it can escalate to thehighly pathogenic type. In this review of LPAI outbreaks in poultry, it is152.5. AI epidemic recurrenceimportant to distinguish between outbreak events caused by non-mutatingand mutating subtypes as the latter type is linked to HPAI outbreaks.H5/H7 subtypesPrior to the discovery that certain subtypes of LPAIV, e.g. H5 or H7, canmutate to HPAIV, little attention was given on identifying LPAI outbreaksin wild bird population. In fact, the relationship between wild birds andLPAIV is assumed to be of a commensal nature [114]. LPAI outbreaks,however, occur frequently in poultry where the association is less benign.Since 2006, H5/H7 LPAI outbreaks were observed in several parts ofEurope [115–117]. Between 2002 and 2006, there were sixty LPAIV strainsisolated from poultry and other captive birds in Europe, Asia, Africa andAustralia [118]. Since many outbreaks may remain unnoticed, the actualnumber of LPAI outbreaks in poultry may be even higher [112].A number of outbreaks can spread over large regions as a result of LPAIVcirculating among poultry. From 1997 to 2001, Italy, for example, hadseveral subtypes of LPAIV circulating in areas populated with a high densityof poultry farms. In 1999 to 2000, an LPAIV (H7N1) that was endemicin poultry gave rise to the emergence of a highly pathogenic variant [119].Although efforts have been done to eradicate H7N1 HPAIV,viruses related toH7 LPAIV appear regularly in Italy and less commonly in other Europeancountries [115]. Similarly, in 1994-1995, another LPAIV (H5N2) causingthe emergence of a HPAIV in Mexico still circulates throughout CentralAmerica.An H7N2 LPAIV that had become endemic in poultry circulated amonglive bird markets and backyard smallholder flocks for 13 years in the north-east region of United States. Fortunately, the virus was finally eliminatedin 2006 without evolving into HPAIV [120]. Several LPAIV subtypes havebeen introduced recently in poultry farms in Germany, the Netherlands,Denmark, and Belgium. For most of these outbreaks, it is assumed thatLPAIVs of H7 subtypes are acquired by direct introduction from their nat-ural reservoir, i.e. wild birds. However molecular analysis suggests that theLPAIV in poultry has become endemic [116, 121].Non-mutating subtypesLPAIVs that do not mutate to HPAIVs may not be a serious problem,however these viruses can still cause severe illness and so should not bedeemed unimportant. Non-mutating LPAIVs continue to spread in various162.5. AI epidemic recurrenceparts of the world: the H3N2 in Italy [122] and China [123], H6N1 in HongKong [124], H6N2 in USA [125], and H9N2 in the Middle East [126]. Inparticular, the subtype H9N2 became very widespread and maintained inpoultry throughout large parts of Asia [127]. According to many case re-ports based on the review by [112], H9N2 has caused high mortality ratesin poultry that is comparable with HPAI outbreaks. Under experimen-tal conditions, H9N2 do not produce severe symptoms in affected poultrybut secondary infections, accompanied by bronchitis virus and Mycoplasmagallsepticum are always involved [128]. Since non-mutating LPAIVs poseharm to poultry and their aggravating symptoms are worrisome, such as thecase of H9N2, eradication (depopulation) or control measures (vaccination)are applied to such outbreaks [116, 129].2.5.2 HPAI outbreaksSeveral descriptions of disease outbreaks or animal plagues are analogousto the present HPAIV outbreaks [130]. An HPAI outbreak was first reportedby an Italian scientist Perroncito [131] who described that the infection ofpoultry around Torino during the Fall and Winter season of 1877-1878. Al-though the infected poultry initially did not show serious illness, the diseasebecame epidemic throughout a larger area causing high mortality amongpoultry [113]. After few decades, the outbreak was found to likely involvethe mutation of LPAIV into an HPAIV [132] as witnessed in the UnitedStates of America [133], Mexico [134], Italy [135], Chile [136] and Canada[137]. Several HPAI outbreaks have been studied from then on.Major HPAI outbreaks occurred around the world in the early 20thcentury, in which poor outbreak detection and management promoted thespread and distribution of the virus. Nevertheless, it seemed that effortswere made to control (e.g. live poultry shipping restrictions, depopulation,quarantine and disinfection) some of these outbreaks [113]. In 1959, thefirst HPAI outbreak took place in Scotland and thereafter other outbreaksoccurred in various places. These HPAI outbreaks are enlisted by Alexan-der and Brown [132], Koch and Elbers [138], Lupiani and Reddy [113], andWorld Animal Health Organization (OIE) [117].Since about 19 years ago, the most worrisome HPAI outbreak, the AsianH5N1 HPAI, emerged and is still severely infecting poultry all around theworld. It has been found that this notorious virus has managed to spill-overfrom poultry to the wild bird reservoir and is maintained in natural hostsdespite its tendency to cause sickness and mortality in these wild birds [138–141].172.5. AI epidemic recurrenceOriginating from Hong Kong and other parts of southern China since2004, it seemed that the Asian H5N1 HPAIV became endemic in birdsthroughout south-east Asia [117]. The virus was said to spread across theHimalayas to the south and westward to Europe and Africa through mi-gratory routes [142]. In 2006, the H5N1 virus prevalence peaked and hasspread to as many as 63 countries across Asia, Europe and Africa. There isevidence that the virus continues to spread, but at a slower rate. In 2008,H5N1 HPAIVs were last detected in Germany. Two years later, the epi-demic ended in Romania, Bulgaria and Russia [117]. However, Asian H5N1HPAIV remains endemic in some Asian countries and Egypt as of 2014 [116].2.5.3 Modelling effortsMathematical modelling has been shown to play an important role in un-derstanding disease dynamics [143, 144]. For AI, models have been used todescribe the disease transmission, evaluate parameters of interest and assessintervention strategies (see [145] and references therein). Other modellingwork emphasizes the importance of combining modelling with experimen-tal or survey data to generate quantitative information that is essential fordesigning control strategies for epidemics [146]. The probability of trans-missions and incidence rates, for example, can be used to assess the role ofimproving bio-security [147]. Other works showing the use of models anddata from other epidemics are described by Garske et al. [148], Mannelliet al. [149], McQuiston et al. [150], Sharkey et al. [151], Nishiguchi et al.[152]. As outlined by Dorjee et al. [145], AI is studied using determinis-tic and stochastic compartmental models [153–158], and network models[151, 159].Past AI models focussed on investigating how the spread of AI to humanscould potentially lead to pandemia [160–162]. A few models were developedinvolving both animals and humans linked by a pathogen. An example isthe simplest model that captures the transmission pathway of H5N1 frompoultry to humans [163]. However, in light of the recurrent outbreaks ofLPAI and HPAI, recent models focus on capturing the pattern of AI out-breaks and identifying the drivers of the disease dynamics. The temporalpattern of AIV prevalence have been used to shed light on the outbreak pe-riodicity of LPAI [1, 164]. These studies support the result of Rohani et al.[165] which emphasises the importance of environmental transmission, i.e.faecal-oral/faecal-cloacal route, in recovering the pattern of periodicity andpersistence of AIV in waterfowls. On the other hand, seasonal patterns ofinfection have been detected where prevalence is highest during late Autumn182.6. Discussionand Winter [3, 42, 76, 166–168]. The patterns are said to be associated withimmune response, bird migration, and/or environmental variability. For in-stance, the high virus prevalence in juvenile mallards for their low immuneresponse and the virus persistence in the environment under colder temper-ature contribute to a spike (or spikes) of influenza activity during wintermonths [169].2.6 DiscussionAvian influenza is indeed a very complex disease in which understandingof its epidemiology, ecology, and epidemic recurrence is limited by monitor-ing methods, genetic sequencing technology, and other data gathering tech-niques. Although it is established that wild birds are the natural reservoir ofAIVs, no accurate tracking survey has been done to show that LPAI-infectedwild birds cause HPAI outbreak in poultry. The recent discovery of HPAIVinfection in wild bird populations poses a challenge to determine whether itis likely that virus in wild birds will become endemic in the case when thevirus ceases to spread in poultry. The diversity and evolution of AIV areother important research subjects. The answers to the research questions onthe biology of each virus subtype and the complex nature of its replicationand interaction with the host cell are useful for vaccine development andother methods to control disease spread.Mathematical models are effective tools to mimic the complex dynamicsof AI epidemics. The models can assist by evaluating current measures andmechanisms to explain outbreak pattern, investigating hypothetical scenar-ios, and making predictions based on empirical data. Some modelling workhas addressed AI systems but more needs to be done. For example, whereassome AI outbreak recurrence models have been used to provide a persis-tence mechanism within host populations, these models were based on thetemporal pattern of AIV prevalence data and not on the incidence of in-fectious cases. It may be true that recurrent prevalence produces recurrentoutbreaks, however such an idea remains poorly tested. We suggest thatcareful attention must be given to using epidemiological terminologies whendescribing models. Further modelling work is needed in order to evaluatethe relative contribution of other possible persistence mechanisms, as thereare other factors that come into play such as spatial and age structures, hostimmunity and evolutionary processes of virus strain. Also other recurrencemodels incorporate seasonality due to environmental factors, but a majordrawback is the lack of consideration for spatial features in the transmis-192.6. Discussionsion dynamics of AI. Nevertheless, these final steps serve as a starting pointto acquire further insights about AI epidemiology, ecology, and epidemicrecurrence.20Chapter 3The Role Of Direct andEnvironmental TransmissionIn Stochastic Avian FluEpidemic Recurrence3.1 BackgroundAvian influenza is an infectious disease present in poultry and wildbirds, and is known to pose threats to humans [170]. The disease pathogenis the avian influenza virus (AIV) whose natural hosts include aquatic birds[83, 171] with wild ducks as its main reservoir [172]. AIV strains can beeither highly pathogenic (HP) or lowly pathogenic (LP) according to theirability to infect hosts. The HPAI viruses are the most virulent and areresponsible for ‘fowl plague’ causing mortality as high as 100% in poultry[66]. The LPAI viruses are endemic in wild bird populations [55, 85] but caneasily be transmitted to domestic stock and then mutate to HPAI type [173].Regardless of the type of virus strain, the prevalence data for avian influenzadisplays recurrent epidemics over time. The underlying mechanism behindthis outbreak pattern is the subject of active investigation [13, 165, 174–176],and is a key consideration in the development of effective control strategiesfor disease mitigation.The virus spreads to healthy individuals either (i) by contact with aninfected host i.e. through inter-host (direct) transmission, or (ii) by hostsacquiring the virus from the environment through drinking or filtering waterwhile feeding, i.e., through environmental (indirect) transmission [13, 175].Recently, a few authors have highlighted the importance of environmentaltransmission as a driver of AIV epidemics. Rohani et al. [165] demonstratedhow neglecting environmental transmission could lead to underestimates forthe explosiveness and duration of AIV epidemics. Breban et al. [164] de-veloped a new host-pathogen model combining within-season transmission213.1. Backgrounddynamics, between-season migration and reproduction, and environmentalvariation, and showed that environmental transmission offers an explanationfor the 2 to 4 year periodicity of AIV epidemics. Wang et al. [1] formulated asimple stochastic model to show that increasing environmental transmissioncan make the outbreak period shorter. Their model predicted an outbreakperiod of 2 to 8 years. Wang et al. [1] found this result consistent withthe observed outbreak period obtained from wavelet analysis of empiricaldata [171]. Together, these papers and others point to environmental trans-mission as the key mechanism behind the approximate periodicity of AIVepidemics. They are based on the assumption that direct transmission isweak between wild birds. Recently, however, a simple SI model with only adirect transmission route, and without stochasticity, was found to providethe best fit to poultry outbreak data [42] suggesting that direct transmissionmay be stronger than originally thought. Motivated by these findings, were-examine, mathematically, the contributions of the different transmissionroutes to the multi-year periodicity of avian flu epidemics.The approximate multi-year periodicity of outbreaks is thought to bedue to the random nature of contagion and recovery processes, togetherwith demographic stochasticity (i.e., uncertainty in birth and death times),and the deterministic dynamics of the system [26]. It has been shown thatdemographic noise can sustain population oscillations that would otherwisedamp to a stable equilibrium. These oscillations are commonly called noisesustained oscillations [31, 177, 178]. Other treatments of this phenomenonin various contexts call it coherence resonance [35], stochastic amplification[32–34], or stochastic resonance [179–181]. Based on known AIV biology,a plausible model [1] suggests that AIV dynamics may arise from noisesustained oscillations.A system exhibiting noise sustained oscillations can be analyzed using arecently developed approximation method. Baxendale and Greenwood [50]showed that a stochastic process of two-dimensional noise sustained oscilla-tions is, in distribution, approximately a rotation whose radius is modulatedby a slowly varying bivariate standard Ornstein-Uhlenbeck (OU) process, awell-studied mean-reverting stochastic process [182].In this thesis, we apply the approximation of Baxendale and Greenwood[50] to a three-dimensional stochastic host-pathogen model, and assess thecontributions of the direct and environmental transmission rates to the re-current epidemics it produces. First, we show that the avian flu epidemicprocess can be approximated by the sum of a scaled univariate OU pro-cess and the product of a rotation and a bivariate slowly varying standardOU processes. Using the approximate stochastic process, we show that the223.2. The stochastic avian influenza modeloutbreak period of the epidemics has a distribution centred at the intrin-sic frequency of the associated deterministic part of the process, i.e., thedeterministic analogue. We obtain the intrinsic frequency as a function ofthe two transmission rates and identify the relationship of each transmissionrate with the dominant outbreak period. Furthermore, we determine howthe outbreak intensity varies over a wide range of direct and environmentaltransmission rates.We organize this chapter as follows. In Section 3.2, we discuss thestochastic avian flu model formulated in terms of stochastic differential equa-tions using the result of Kurtz [10]. In the same section, we also suggestan appropriate range of values to use for the transmission rate parameters.We review previous analytic work. Section 3.3 contains our own analysis,which includes determining the approximate process, the theoretical powerspectral density (PSD), and the formula for the intrinsic frequency. Wepresent in Section 3.4 the interpretation of our analysis, that the disease re-currence observed in stochastic simulations is approximately governed by arotation matrix multiplied by a standard Ornstein-Uhlenbeck (OU) process.We then describe the influence of each transmission route on the dominantperiod of outbreaks (based on the intrinsic frequency) as well as the intensityof outbreaks (based on the stationary standard deviation of the approximateprocess). Finally, in Section 3.5, we discuss the implications of our resultsfor characterizing avian flu epidemics and providing insights into the depen-dence of the frequency and variation of recurrent avian influenza epidemicson each transmission route.3.2 The stochastic avian influenza model3.2.1 Model descriptionWang et al. [1] formulated a stochastic model for avian influenzausing the susceptible-infected-recovered (SIR) framework with an added en-vironmental transmission rate. They found that their model was sufficientto explain the multi-year periodicity of flu outbreaks. In this host-pathogenmodel, there is a susceptible duck population of size S, an infected pop-ulation of size I, a population of recovered individuals of size R, and anenvironmental virus concentration V . Note that S, I, and R are integers,and all populations are functions of time. We also assume that a new sus-ceptible duck is born once a host dies so that the total duck population isconstant, i.e. S+ I +R = N . Thus, we can solve for R in terms of S and I.By considering a short time interval [t, t+ ∆t] and denoting by T (σ′|σ) the233.2. The stochastic avian influenza modeltransition rate from state σ = (S, I, V ) to σ′ = σ ± ν, where νi ∈ {0, 1},the different system events and their corresponding transition rates are:1. Infection T (S − 1, I + 1, V |S, I, V ) = βS IN+ ρSVNV(3.1a)2. Birth and Death T (S + 1, I, V |S, I, V ) = µ(N − S − I), (3.1b)T (S, I + 1, V |S, I, V ) = µI, (3.1c)T (S, I, V + 1|S, I, V ) = τI + δV, (3.1d)T (S, I, V − 1|S, I, V ) = ηV. (3.1e)3. Recovery T (S, I − 1, V |S, I, V ) = γI. (3.1f)We schematically present these processes in Figure 3.1.The first event (3.1a) describes the infection of a susceptible individ-ual. Infection happens when a susceptible host is in close contact with aninfected host or when it acquires the virus directly from the environment.The rate of transmission of the disease is given by the likelihood of contactbetween a susceptible individual and either infected individuals or virionsin the environment, multiplied by the rate at which virions are acquired bysusceptible individuals in each case. Note that the likelihood of contact de-pends on the fraction of infected individuals (I/N) and the concentration ofvirus in the environment normalised by a reference concentration (V/NV ),i.e., the likelihood of contact in a frequency-dependent process.The second category of events (3.1b)-(3.1e) encompasses the birth anddeath processes of the host and the virus. For simplicity, it is assumed herethat the per capita host birth and death rates have the same value µ. Weassume that virus is introduced into the environment at a constant rateδ from alternative hosts. The virus concentration in the environment alsogrows when infected ducks shed virions; this shedding occurs at rate τ . Theclearance rate of the virus in the environment is η.Finally, the third category of events (3.1f) is the recovery of infectedducks at per capita rate γ.The parameters β, µ, τ , δ, η and γ are stochastic rates.243.2. The stochastic avian influenza modelFigure3.1:Aschematicdiagramofthehost-pathogenmodelforavianinfluenzaadaptedfromWangetal.[1].Thehostpopulationhasthreecompartmentscorrespondingtosusceptible,infected,andremovedindividuals.Thesolidarrowsrepresentmovementofindividualsfromonecompartmenttoanotherasaresultofbirth,death,infection,andrecoveryprocesses.Thedashedarrowsrepresenttheincreaseordecreaseinthenumberofinfectedindividualsduetotheinteraction(squigglyarrow)ofthesusceptiblehostwiththeavianinfluenzaviruspresentintheenvironment.253.2. The stochastic avian influenza modelThe resulting stochastic avian flu host-pathogen model (see Appendix A.1)is approximated for large N by the following system of stochastic differentialequations (SDEs):ds = (−βsi− ρsv + µ(1− s)) dt+ 1√N(−G1dW1 +G2dW2 +G3dW3) ,di = (βsi+ ρsv − (µ+ γ)i) dt+ 1√N(G1dW1 −G3dW3 −G4dW4) ,dv = (kτi+ δv − ηv) dt+ 1√NV(G5dW5 −G6dW6) ,(3.2)where,G1 =√βsi+ ρsv, G2 =√µ(1− s− i), G3 =√µi,G4 =√γi, G5 =√kτi+ δv, and G6 =√ηv.(3.3)Here, s = S/N, i = I/N, v = V/NV , and k = N/NV . The SDE for theproportion of recovered ducks r(t) is not necessary because we eliminate r(t)using r(t) = 1− s(t)− i(t), which follows from S + I +R = N . In (3.2), thesecond term vanishes as N,NV →∞ which leads us to the deterministic, orso-called mean-field, dynamics as found in Wang et al. [1]:φ˙1 = −βφ1φ2 − ρφ1ψ + µ(1− φ1),φ˙2 = βφ1φ2 + ρφ1ψ − (µ+ γ)φ2,ψ˙ = κτφ2 + δψ − ηψ.(3.4)The deterministic variables φ1, φ2, ψ represent fractions of the susceptiblehosts, infected hosts, and virus in the environment, respectively, while κ =limN,NV→∞NNV. The disease is epidemic in model (3.4) if the number ofinfected ducks increases, i.e. φ˙2 > 0, which implies thatβφ1µ+ γ+ρφ1ψ(µ+ γ)φ2> 1. (3.5)Note that φ˙1 = 0 and ψ˙ = 0 imply φ2 =µ(1− φ1)µ+ γand ψ =φ2κτη − δ ,respectively. We then write the inequality (3.5) in terms of φ1 and obtainβφ1µ+ γ+ρκτφ1(µ+ γ)(η − δ) > 1. (3.6)263.2. The stochastic avian influenza modelAt the outset of an epidemic, nearly all hosts are susceptible (φ1 ≈ 1).By substituting φ1 = 1 into (3.6), we arrive at the condition for the diseaseto persist,R0 ≡ βµ+ γ+κρτ(η − δ)(µ+ γ) > 1, (3.7)whereR0 is the basic reproduction number, the average number of secondaryinfections resulting from one infected individual in a susceptible population.The deterministic system has a stable endemic equilibrium, which canbe written in terms of R0,(φ∗1, φ∗2, ψ∗) =(1/R0, µµ+ γ(1− 1/R0), κµτ(η − δ)(µ+ γ)(1− 1/R0)).(3.8)When the basic reproductive number R0 > 1, the disease is epidemic. Re-sults (3.4)-(3.8) are also found in Wang et al. [1].We are interested in characterizing the fluctuations of the stochasticmodel around the steady-state solution. Hence, we linearize (3.2) around(φ∗1, φ∗2, ψ∗) and obtain the linear diffusion equation (see Appendix A.2 forthe detailed derivation),dξ = A0ξ dt+ C0 dW, ξ(t),W(t) ∈ R3,A0,C0 ∈ R3×3. (3.9)whereA0 =−βφ∗2 − ρψ∗ − µ −βφ∗1 −ρφ∗1−βφ∗2 βφ∗1 − µ− γ ρφ∗10 κτ δ − η , (3.10)andC0 =C11 C12 0C21 C22 00 0 C331/2 ,C11 = βφ∗1φ∗2 + ρφ∗1ψ∗ + µ(1− φ∗1),C12 = C21 = −βφ∗1φ∗2 − ρφ∗1ψ∗ − µφ∗2,C22 = βφ∗1φ∗2 + ρφ∗1ψ∗ + (µ+ γ)φ∗2,C33 = κτφ∗2 + δψ∗ + ηψ∗.(3.11)Equations (3.9)-(3.11) are the Langevin equations derived by Wang et al.[1]. The drift coefficient matrix A0 is the Jacobian of the deterministic model(3.4) evaluated at the endemic equilibrium (φ∗1, φ∗2, ψ∗). The diffusion matrixC0 is formed using the coefficients of the independent Wiener processes in273.2. The stochastic avian influenza model(3.2). This matrix is the square-root of the covariance matrix B found byWang et al. [1]. This chapter focusses on the study of the linear system(3.9).3.2.2 Parameter valuesAll simulations produced in this work use as default parameters the val-ues in Table 3.1 largely taken from Wang et al. [1]. The values of parametersµ, η, β, and γ are based on empirical studies in the literature (see captionof Table 3.1). No data are available for ρ and δ and so they are variedwithin the range used by Wang et al. [1]. Note that the unit of the sheddingrate τ is virion/mL/duck/year rather than virion/duck/day as erroneouslyreported in Wang et al. [1].The values for the transmission parameters β and ρ can vary widely withseasonal climate changes and from one geographic area to another. We thusgive a range of values for β and ρ and study the system’s response to differentlevels and modes of transmission. As displayed in Table 3.1, we use a widerrange of β values than was used by Wang et al. [1]. Observe that the infectionterm in (3.4) is given by βφ2 = βI/N (rather than simply βI). In ecologicalterms, this means that β is the transmission rate for a frequency-dependentprocess rather than a density-dependent process. Roche and Lebarbenchon[175] find that β ranges from 0.00005 to 1500. Here, we plot our results forβ values ranging from 0 to 300.Avian influenza epidemic models exhibit noise sustained oscillations witha nonzero dominant frequency for a small range of β values, with some fixedvalue of ρ. Knowing only that β falls within a wide range of values, re-investigation of the relative contribution of each of the transmission modesis necessary. In particular, a larger β poses the possibility that direct trans-mission (β) alone, along with stochasticity, can drive avian flu epidemicswith multi-annual periodicity.283.2. The stochastic avian influenza modelTable3.1:Descriptionsoftheparametersintheavianflumodelandthevaluesusedinstochasticsimulations.Therangeofvaluesforβinbold-facedisthesetofvaluesthatwereusedin[1].ParameterValue/RangeDescriptionUnitN103DuckpopulationsizeduckNV105Virusconcentrationintheenvironmentvirion/mLµ0.3Duckbirthanddeathrateyear−1δ0.1Virusreplicationrateyear−1τ104VirussheddingratevirionmL−1/duck/yearη3Virusclearancerateyear−1β0.00005Directtransmissionrateduck−1year−1to1500(0-0.05)ρ0−3Environmentaltransmissionrateyear−1γ5.5Duckrecoveryrateyear−1293.2. The stochastic avian influenza modelAssessing the relative contribution of each of the transmission rates todisease recurrence requires a quantitative comparison of results. From (3.7),we see that if bothβ < µ+ γ and ρ <(η − δ)(µ+ γ)κτ, (3.12)then R0 < 1. Using the values in Table 3.1, we find that when both β < 5.8and ρ < 0.16 we have R0 < Preliminary analysis of the modelIn this section we highlight key results from the analysis of the modelof Wang et al. [1], and extend these results to include the role of R0 in thedynamics of avian flu. It is an established fact that the oscillations of adamped system can be sustained by stochasticity. Wang et al. [1] showedthat the deterministic version of model (3.4) has a stable endemic steady-state (3.8) when R0 > 1. Here, we take the calculation one step further bydetermining the parameter ranges where (3.8) is a stable sink or a stablefocus. We do this by plotting the eigenvalues of the Jacobian in (3.10)against R0 from (3.7). Figure 3.2 displays this plot.In Figure 3.2, we see the eigenvalues of A0 plotted as functions of R0.We use the parameter values in Table 3.1 with β = 0.05. As ρ increases,so does R0. We observe that for R0 < 1, all of the eigenvalues are realand one has positive sign. As expected from Wang et al. [1], the endemicsteady-state is unstable for this case and the system evolves to the disease-free equilibrium. A transition in the signs of the eigenvalues happens atR0 = 1. When R0 > 1, the system gives rise to a complex-conjugate pair ofeigenvalues with negative real parts, and a negative real eigenvalue. Thus,the steady-state (φ∗1, φ∗2, ψ∗) is a stable focus, i.e., the deterministic systemexhibits damped population cycles. Note that a similar result is obtained ifwe fix ρ and let R0 vary with β.Suppose that we represent the complex eigenvalues as −λ ± iω whereλ and ω are magnitudes of the real and imaginary parts, respectively. ForR0 > 1 in Figure 3.2, λ is clearly smaller than ω, because the upper solidcurve is above the lower dashed curve with both curves drawn below thex-axis. Moreover, observe that ω increases faster than λ. Hence we deducethat the ratio λ/ω decreases as R0 increases.303.2. The stochastic avian influenza model0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−10−8−6−4−202EigenvaluesBasic reproductive number (R0)The real (solid) and imaginary (dotted) parts of the eigenvaluesFigure 3.2: The real (solid curves) and imaginary (dash-dot curves) parts ofthe eigenvalues of A0 in (3.10) associated with the stability of the endemicsteady-state (3.8) plotted against R0 in the case where β = 0.05. There arethree real eigenvalues when R0 < 1, of which two are negative (thick solidlines). For R0 > 1, we have a complex-conjugate pair of eigenvalues withnegative real parts (thick solid lines) and another negative eigenvalue (thinsolid line). The imaginary parts of the complex eigenvalues are shown bythe dotted lines. Other parameter values are shown in Table 3.1.In Figure 3.3, we plot typical stochastic realizations of the avian flumodel (3.2) for R0 > 1 using various combinations of parameters. Thesestochastic paths display oscillations sustained by noise.We display the plots of stochastic realizations for the case when β = 0.05but with different ρ values in Figures 3.3(a) and 3.3(b). Here we noticethat higher amplitude and higher frequency of epidemics are observed inFigure 3.3(b) as compared to the fluctuation plotted in Figure 3.3(a) whereρ is twice as much. This observation is consistent with Figure 3.2.313.2. The stochastic avian influenza model10 20 30 40time (year)00.050.1Infectives (proportion)β=0.05, ρ=0.2, R0=1.2, λ=0.17, ω=0.29deterministicstochastic(a)10 20 30 40time (year)00.050.1Infectives (proportion)β=0.05, ρ=0.4, R0=2.39, λ=0.31, ω=0.84deterministicstochastic(b)10 20 30 40time (year)00.050.1Infectives (proportion)β=0, ρ=0.42, R0=2.5, λ=0.32, ω=0.87deterministicstochastic(c)10 20 30 40time (year)00.050.1Infectives (proportion)β=14.5, ρ=0, R0=2.5, λ=0.38, ω=1.57deterministicstochastic(d)Figure 3.3: Simulation of the stochastic model (3.2) and its correspondingdeterministic solution for β = 0.05 with (a) R0 = 1.2 and (b) R0 = 2.5, andfor R0 = 2.5 with (c) β = 0 (ρ = 0.4205) and (d) ρ = 0 (β = 14.5).In Figures 3.3(c) and 3.3(d), we plot stochastic paths for the case wheneither β = 0 or ρ = 0 but with roughly the same R0. In this comparison,R0 ≈ 2.5, the upper bound of the confidence interval for the estimate of R0for avian influenza in wild birds [183]. Comparing Figures 3.3(c) and 3.3(d),we observe that the periodicity and intensity of outbreaks in the two cases isdifferent, even though the basic reproduction number is the same in the twocases. Current theory points toward the transmission mode as a determiningfactor for understanding the periodicity and intensity of avian flu outbreaks.Here, we develop a method that allows us to determine mathematically theeffect of each transmission mode on the recurrence (periodicity and intensity)323.3. Analytic methodsof avian influenza.3.3 Analytic methodsIn this section we develop the analytic tools we need to understand thecontribution of each transmission route to the dynamic behaviour of themodel. In particular, our goal is to develop a mathematical description ofthe noise sustained oscillations that are observed in stochastic simulations,such as those shown in Figure Approximate solutionsWe begin by defining an approximation for the process of oscillationsproduced by the stochastic avian flu model. Our starting point is the lineardiffusion equation given by (3.9), which describes the process ξ = (ξS , ξI , ξV )near the endemic steady-state for large time. Under certain conditions,an approximate solution of (3.9) can be obtained using an extension of aresult from Baxendale and Greenwood [50]. Specifically, if the eigenvaluesof the drift coefficient matrix A0 (see (A.41) below) are of the form −ζ, and−λ ± iω with λ/ω small, for ζ, λ, ω > 0, an approximate solution for (3.9)(see Appendix A.3 for details) is given byξapp(t) = y1(t)Q•1 +σ˜√λ[Q•2,Q•3]R−ωtSλt, (3.13)where Q is the canonical form of the eigenvector matrix associated withthe three eigenvalues. The vector Q•j denotes the jth column of Q. Thestochastic process y1(t) is an Ornstein-Uhlenbeck (OU) process [182] withmean zero, decay rate ζ, and diffusion coefficient σ1 obtained using thefirst row of the matrix Q−1C0. The matrix R−ωt is a rotation matrix. Itdescribes the circular motion of the process with frequency ω. The vectorprocess Sλt is a bi-variate OU process with independent components. Thescalar σ˜ is determined by the last two rows of Q−1C0 (the full expressionis given by (A.22) in Appendix A.3). We express the large-time stationarysolution of (3.2) as:s(t) = φ∗1 +1√NξS(t) ≈ φ∗1 +1√NξappS (t),i(t) = φ∗2 +1√NξI(t) ≈ φ∗2 +1√NξappI (t),v(t) = ψ∗ +1√NVξV (t) ≈ ψ∗ + 1√NVξappV (t).(3.14)333.3. Analytic methodsThe formulation (3.14) shows that the solution of (3.2) near the endemicsteady-state behaves approximately like the product of a rotation and anOU process (3.13).3.3.2 Power spectral densityWe use the theoretical power spectral density (PSD) to determine thedistribution of frequency components within the stochastic process producedby (3.9) with (3.10) and (3.11), and by their approximate form (3.13). Alinear diffusion process is described as a general multivariate OU process,i.e.,dx = −Ax(t) dt+ B dW(t). (3.15)According to Gardiner [46], the PSD of an n-dimensional process is obtainedfrom the main diagonal of the matrix given by:S(f) =12pi(A + if)−1BBT (AT − if)−1. (3.16)From (3.16), it follows that the PSD of ξ(t) satisfying (3.9) is obtained fromthe main diagonal of the matrixS(f) =12pi(−A0 + if)−1C0CT0 (−AT0 − if)−1. (3.17)On the other hand, the theoretical PSD of the approximate solution ξapp(t)in (3.13) is determined using the coefficients of (see Appendix A.3 for deriva-tion)dξapp(t) = A0ξapp dt+ Q · diag(σ1, σ˜, σ˜) dW(t).Thus, the PSD of ξapp is obtained fromS(f) =12pi· diag(σ21, σ˜2, σ˜2) · (−A0 + if)−1QQT (−AT0 − if)−1. (3.18)3.3.3 The intrinsic frequency and the decay rate of thedeterministic dynamicsThe approximation (3.13) depends explicitly on the intrinsic frequencyω and the decay rate λ of the damped oscillations predicted by the de-terministic system (3.4). We can derive these quantities by obtaining theeigenvalues of the Jacobian matrix A0, which are the solutions of the char-acteristic polynomial given byν3 − aν2 − bν − c = 0, (3.19)343.4. Resultswhere:a = (δ − η)− µR0 − γ − µ+ β/R0,b = −µ(η − δ + γ + µ)R0 + µβ/R0,c = −µ(η − δ)(γ + µ)(R0 − 1).(3.20)In Appendix A.6, we derive λ and ω as functions of β and ρ, using theroots of (3.19) and the avian flu parameters in Table 3.1. The results areshown in Figures 3.6 and 3.7(a) where we see λ/ω and ω are functions of βand ρ.3.3.4 Numerical tools and functionsAll numerical computations were done using MATLAB [184]. We com-puted the solution of the deterministic avian flu model (3.4) with the built-infunction ode45(), an ordinary differential equation solver that is based onthe Runge-Kutta (4,5) formula. Stochastic simulations of (3.2) and (3.9)were done using the Euler scheme (or the Euler-Maruyama scheme) [185],a first-order discretization method for stochastic differential equations. Atime step of 0.01 was used for simulations of (3.2), (3.4), and (3.9).3.4 ResultsWe present our results in two parts. First, we substitute the parametervalues from Table 3.1 into the approximate solution (3.13) of (3.9) with(3.10) and (3.11). Second, we use the approximation to understand how thedifferent transmission routes affect the dominant periodicity and the typicalintensity of epidemics (see Section 3.4.2).3.4.1 An approximate avian flu epidemic processIn this section, we use the approximation (3.13) to describe explicitlythe noise sustained oscillations we saw in Figure 3.3. For consistency withprevious work [1], we have chosen β = 0.05 and ρ = 0.4. This choiceof parameters results in the endemic steady-state φ∗1 = 0.419, φ∗2 = 0.03,ψ∗ = 1.04 with R0 ≈ 2.39. The deterministic process defined by (3.4)353.4. Resultspersists. The diffusion equation (3.9) with (3.10) and (3.11) becomesdξ = A0ξ dt+ C0 dW, whereA0 =−0.716 −0.021 −0.1680.410 −5.78 0.1680 100 −2.9 , andC0 = 0.568 −0.161 0−0.161 0.568 00 0 2.494(3.21)The eigenvalues of the Jacobian A0 in (3.21) are −ζ = −8.78 and −λ± iω =−0.309± 0.838i. The ratio λ/ω = 0.3090.838 ≈ 0.369 is sufficiently small that wecan approximate the process ξ(t) by (3.13), as evidenced by the comparisonof PSDs in Figure 3.4.Using the eigenvalues above, we arrive at the canonical matrix of eigen-vectors, i.e., the columns are eigenvectors of A0, given byQ = 0.021 −0.078 0.16−0.059 0.026 0.0080.998 0.984 0 . (3.22)The matrix C0 in the diffusion term of (3.21) gives usΣ = Q−1C0 = 2.13 −6.43 0.834−2.16 6.52 1.692.23 3 0.715 . (3.23)We compute σ1 (see Appendix A.3) by taking the norm of the first row ofthe matrix Σ in (3.23),σ1 = ||(2.13,−6.43, 0.834)|| ≈ 6.82.We form the matrix C˜ = (Σ˜Σ˜ᵀ)1/2 where,Σ˜ =[−2.16 6.52 1.692.23 3 0.715], (3.24)to obtain σ˜2 = Tr(C˜C˜ᵀ)/2 and so σ˜ ≈ 5.68.Using (3.13), the solution to (3.21) near the endemic steady-state isapproximately363.4. ResultsξappS (t)ξappI (t)ξappV (t) = y1(t) 0.021−0.0590.998+ 10.2−0.078 0.160.026 0.0080.983 0R−0.838tS0.309t,(3.25)wheredy1 = −8.78y1 dt+ 6.82 dW1. (3.26)One way to see if the process (3.25) is a reasonable approximate solutionfor (3.21) is to compute the theoretical PSDs of the exact solution ξ(t) of(3.21) and the approximation ξapp(t) in (3.25) using formulae (3.17) and(3.18), respectively. In Figure 3.4, the PSDs of fluctuations ξ(t) and ξapp(t)agree fairly well. The dominant frequency is close to the intrinsic frequencyω = 0.8377.Figure 3.4: Comparisons between the theoretical PSD of the exact processξ(t) (solid line) satisfying (3.21) and the approximate process ξapp(t)(dashedline) given by (3.25), for the fluctuations of the susceptible, infectious, andthe virus populations. Default parameter values are in Table 3.1 with β =0.05 and ρ = 0.4.The process ξapp(t) depends on the OU process y1(t) whose dynamicsare described by (3.26), an OU process with asymptotic mean zero, decayrate 8.78, and diffusion coefficient 6.82. From the first term of (3.25), wefind that the contribution of the process y1(t) to host fluctuation processes,ξS(t) and ξI(t), is very small as compared to that of the virus fluctuationprocess ξV (t). We use the stationary standard deviation to measure the373.4. Resultstypical amplitude of the population fluctuations. Note that the stationaryvariance of the process y1(t) is 46.51/(2 × 8.78) ≈ 2.65, i.e. stationarystandard deviation√2.65 ≈ 1.63. The stationary standard deviation, i.e.typical amplitude, of the stochastic path for the process ξapp(t) (see (3.29)or (A.30) in Appendix A.3), is given bySSDi ≈√2.65q2i1 + 10.22r2i . (3.27)The computed typical amplitudes of the susceptible, infectious, and viruspopulation fluctuations are 1.81, 0.29, and 10.18, respectively. In Figure 3.5,the stochastic realizations for the different population fluctuations given by(3.21) and (3.13) are plotted along with their stationary standard deviationsindicated by horizontal lines to indicate that the stochastic oscillations gen-erally lie within one stationary standard deviation.10 20 30 40time (year)-1.810.001.81S (fluctuation)Susceptibleexactapp10 20 30 40time (year)- (fluctuation)Infectiousexactapp10 20 30 40time (year)-10.180.0010.18V (fluctuation)VirusexactappFigure 3.5: A stochastic realization of the population fluctuations by sim-ulating (3.21) (solid line) and (3.25) (dashed line) and their correspondingstationary standard deviations, i.e. typical amplitudes (in horizontal lines)computed using (3.28).383.4. ResultsWriting (3.25) in polar form (see Appendix A.3), we obtainξS(t) ≈ 0.021y1(t) + 1.82|S0.309t| cos (ϕ0.309t − 0.838t− 2.02),ξI(t) ≈ −0.059y1(t) + 0.275|S0.309t| cos (ϕ0.309t − 0.838t− 0.313),ξV (t) ≈ 0.998y1(t) + 10|S0.309t| cos (ϕ0.309t − 0.838t).(3.28)We know from the approximate form (3.28) that the phase differences be-tween the population fluctuation processes ξi(t) are constants. For in-stance, the susceptible duck population fluctuation is out of phase withother populations in the system. With respect to the virus population,the susceptible duck population exhibits noise-sustained oscillations witha phase advance, with respect to the virus population, of approximately2.02/ω = 2.02/0.838 ≈ 2.4 years. The infectious duck population, on theother hand, oscillates with a phase advance of 0.313/0.838 ≈ 0.37 years fromthe virus population.3.4.2 The relative contribution of the different transmissionroutesOne of the key factors leading to outbreaks is the rate of transmission.As outlined previously, there are two transmission routes for avian flu: directtransmission of virus from one duck to another, and indirect transmission ofviruses via the environment. We use the approximation ξapp(t) to shed lighton the relative importance of direct and indirect transmission to outbreakoccurrence.As a first step, we identify the constraints for β and ρ under which oneexpects to observe recurrent epidemics and where the approximation is valid.We use the parameter ranges from Table 3.1 and display in the top panelsof Figure 3.6 a plot of λ/ω as a function of β and ρ.393.4. Results(a) (b)(c)Figure 3.6: (a,b) Plot of the ratio λ/ω as a function of β and ρ. Thetriangular white region in the lower left is where no recurrence is observed.Panel (a) considers 0 ≤ β ≤ 50 and 0 ≤ ρ ≤ 1 while Panel (b) considers0 ≤ β ≤ 300 and 0 ≤ ρ ≤ 3. The grey region in both panels is whereλ/ω ≤ 0.35. In Panel (b), the black curve corresponds to λ/ω = 1. (c) Plotof the basic reproduction number R0 as a function of β and ρ. The darker(blue online) shade indicates that no epidemic can be observed in the modelfor this parameter region. Parameter values and ranges from the literatureare in Table 3.1.In Figure 3.6(a), we see a triangular region surrounding β = ρ = 0 wherethe approximation is no longer relevant because no recurrent epidemics canbe observed there. In the triangular parameter region, R0 ≤ 1 as seen inFigure 3.6(c) and there are no recurrent epidemics. Within a large portionof the remaining parameter region defined by 0 ≤ β ≤ 15 and 0 ≤ ρ ≤ 1 (see403.4. ResultsFigure 3.6(a)), λ/ω is sufficiently small that the approximation is valid.We further extend the range of transmission parameters to 0 ≤ β ≤ 300and 0 ≤ ρ ≤ 3 as displayed in Figure 3.6(b) to show the behaviour of λ/ωfor larger values of β and ρ. We observe that when β and ρ are both verylarge, as depicted by the white region in the upper portion of Figure 3.6(b),λ/ω > 1 and so the approximation is not necessarily valid. In Figure 3.6(b),a black curve is drawn, which corresponds to λ/ω = 1. The approximationis valid for β and ρ values in the region well below this curve. For instance,it is guaranteed that the approximation can be used to precisely describethe system in our parameter space when 0 < λ/ω ≤ 0.35 (grey region).Since there is a sufficient portion of the region defined by 0 ≤ β ≤ 300 and0 ≤ ρ ≤ 3 where the approximation is applicable, the following analysisfocusses on the region well below the black curve to study how the differentroutes of transmission influence the periodicity and intensity of recurrentepidemics reflected by the model.Dominant outbreak periodThe period of any oscillating function is the inverse of the frequency.We analyze the dominant outbreak period of the simulated epidemic from(3.9) by considering the intrinsic frequency of the deterministic system. Fora parameter range where the approximation (3.13) is close to the exactprocess satisfying (3.9) (see Figure 3.4), the dominant frequencies predictedby the two are very close. In addition, this frequency is close to that fromthe deterministic system when λ is small [186]. Hence we use ω given byformula (A.44) (see Appendix A.6).Using formula (A.44) in Appendix A.6, we compute the intrinsic fre-quency ω = ω(ρ, β) over 0 ≤ ρ ≤ 3 and 0 ≤ β ≤ 300 for the parameters inTable 3.1 to obtain Figure 3.7(a).We are interested in the parameter domain when R0 > 1, that is, wherethe disease has recurring epidemics (see Figure 3.6(c)). In Figure 3.7(a),we find a region where the intrinsic frequency is relatively high, namelyfor low environmental transmission (ρ ≤ 0.5) and high direct transmission(β ≥ 100). When the β and ρ values are both very low or very high, theintrinsic frequency is near zero, and so all populations fluctuate around theirendemic equilibrium very slowly.The steady-state infectious proportion is approximately 0.05 and is ac-companied by a very low steady-state proportion of susceptibles (see Fig-ure 3.7(b)). Thus in this parameter region, the host population is mainlycomposed of individuals that have recovered from the disease.413.4. ResultsWe also find that, for (β, ρ) values where the approximation is valid, thedominant outbreak period is 2 to 8 years, as seen by comparing the greyregion of Figure 3.6(b) with Figure 3.7(a).(a)(b)Figure 3.7: Plots of (a) the intrinsic frequency ω and (b) the steady-stateproportions of the susceptible and infectious populations as functions of thedirect transmission rate β (left) and the environmental transmission rate ρ(right). The dark region in (a) is where the 2 to 8 year recurrence period isobserved. Parameter values are in Table 3.1.423.4. ResultsFor the default parameter values with β = 224 and ρ from 0 to 0.05, thelargest 5.6 radians per year, corresponding to an avian flu outbreak periodof 2pi/5.6 ≈ 1.12 years. For the default parameter values with β = 224, thehighest intrinsic frequency reached is 5.6 radians per year corresponding to0 < ρ < 0.05. This maximum intrinsic frequency corresponds to an avianflu outbreak period of 2pi/5.6 ≈ 1.12 years (Figure 3.7(a)).We determine the outbreak periodicity for various values of either β orρ alone by computing the theoretical PSD of the linear system (3.9) with(3.10) and (3.11) (Figure 3.8). When there is no direct transmission, thedominant frequency of the simulated outbreaks increases with the environ-mental transmission rate (Figure 3.8(a)). Similarly, the dominant frequencyincreases with the direct transmission rate in the absence of environmentaltransmission (Figure 3.8(b)).In parallel with the increasing dominant frequency, the width of thePSD also increases with increasing environmental transmission rate. Thismeans that, at low transmission rates, the outbreak pattern is more regular,and becomes more irregular as transmission rate increases. Consequently,when transmission rates are high, the timing of outbreaks is more difficultto predict.Note here that we have chosen values of ρ (when β = 0) and β (whenρ = 0) that give roughly the same basic reproduction number R0, allowingfor comparison between the theoretical PSDs. Consider the PSDs associ-ated with (ρ, β) = (0.22, 0) and (ρ, β) = (0, 7.5) in Figure 3.8(b) and 3.8(a),respectively. Both parameter values correspond to R0 ≈ 1.3, a value usedby [1] to define a boundary between disease persistence and stochastic ex-tinction for their model. If we compare the PSDs between Figures 3.8(a)and 3.8(b), we find that, with the values in Table 3.1, the model predictslarger-valued PSDs in the case when environmental transmission is absentthan when direct transmission is absent.The corresponding dominant outbreak period is shown in Figure 3.8(c)and 3.8(d). Both dominant periods of the approximate and exact process areshown. The results are indistinguishable. We also observe in Figure 3.8(c)that, for β = 0, an increase in ρ from 0.22 to 1.1 results in a decline ofapproximately 12 years in the outbreak period. However, for ρ = 0, anincrease in β from 7.5 to 37.5 only results in a decrease of approximately7 years in the outbreak period (see Figure 3.8(d)). In other words, whenthe disease is epidemic but R0 is close to 1, increasing the environmentaltransmission can cause a greater drop in the outbreak period than increasingdirect transmission. However, for larger transmission rates we find that thedominant outbreak period declines slowly as transmission rate increases and433.4. Results0 0.2 0.4 0.6 0.8 1Frequency (per year)10-3100Log PSD (infectious)Power spectral density when β = 0ρ=0.22ρ=1.1ρ=2.2(a)0 0.2 0.4 0.6 0.8 1Frequency (per year)10-3100Log PSD (infectious)Power spectral density when ρ = 0β=7.5β=37.5β=75(b)0 1 2 3Environmental transmission rate ρ (per year)2468101214161820Dominant outbreak period (year) No direct transmission (β=0)R0=1.3R0=6.5 R0=13(c)0 20 40 60 80 100Direct transmission rate β (per year)2468101214161820Dominant outbreak period (year) No environmental transmission (ρ=0)R0=1.29R0=6.47 R0=12.9(d)Figure 3.8: Top panels: Theoretical PSDs of the infectious population fluc-tuations ξI(t) for (a) ρ = 0.22, 1.1, 2.2 when β = 0 and (b) β = 7.5, 37.5, 75when ρ = 0. The linewidth of the PSD curves increases with the trans-mission parameter values. Bottom panels: Approximate dominant outbreakperiod, 2pi/ω with formula (A.42), as a function of (c) ρ when β = 0 and(d) β when ρ = 0. The square markers in (c) and (d) are located atρ = 0.22, 1.1, 2.2 and β = 7.5, 37.5, 7.5, respectively. The black dots rep-resent the exact dominant outbreak period obtained using the theoreticalPSDs in the top panels. The black dots and square markers are indistin-guishable from each other. The horizontal lines in (c) and (d) indicate the2 to 8 year periods observed in actual prevalence data (See [1]). The corre-sponding R0 values are also shown. Default values for all other parametersare given in Table 3.1.443.4. Resultsso the effects of changes in the individual transmission rates may be difficultto distinguish.Typical intensity of outbreaksFor each i = S, I, V , the stationary standard deviation of ξappi (t) isSSDi =√q2i1σ212ζ+ r2iσ˜2λ(3.29)for r2i = q2i2 + q2i3 where qij are entries of the matrix Q, e.g. (3.22), todetermine the typical intensity of outbreaks. Here we show how the typicaloutbreak intensity is influenced by each type of transmission. In Figure 3.9,we display the plots of SSDi’s of each process ξappi (t) as functions of β andρ.In Figure 3.9, we have focused on the (ρ, β) region that correspondsto R0 > 1 and λ/ω small, i.e. where the approximation method is valid.Typical amplitudes in simulated epidemics when 0 ≤ ρ ≤ 1 and 0 ≤ β ≤ 100are depicted in Figure 3.9.In Figure 3.9(b), we see that in general, fluctuation amplitudes are higherfor larger values of β, the direct transmission rate. On the other hand, inFigure 3.9(a), the fluctuation amplitudes for the susceptibles are lower forlarger β, and there is high sensitivity to β. Figure 3.9(c) shows an optimalregion in β for the fluctuation amplitude of virus especially for low ρ.453.5. Discussion(a) (b)(c)Figure 3.9: Typical amplitude (intensity) of the fluctuations in the propor-tion of (a) susceptible, (b) infectious, and (c) virus populations, measuredby their stationary standard deviations (3.29) as a function of the environ-mental transmission rate ρ and β. All other parameters are at the defaultvalues in Table DiscussionWe have written the Wang et al. [1] stochastic model for avian influenzaincluding direct and environmental transmission routes as stochastic dif-ferential equations using the method of Kurtz [10], following Greenwoodand Gordillo [47]. Under large host and virus populations, the stochasticmodel approaches the deterministic system wherein, for R0 > 1, the en-demic steady state is a stable focus. Our discussion of the avian flu model463.5. Discussionin Section 3.2 suggests that the disease can persist (R0 > 1) if either one ofthe two transmission routes is sufficiently strong. We have also shown viastochastic simulations that the model gives rise to noise sustained oscillationin the presence of either transmission route.Our analysis allows us to conclude that the temporal pattern of epi-demic recurrence is the sum of two processes: (1) a scaled OU process withlong-term mean zero, and (2) the product of a rotation and a slowly varyingstandard OU process in two dimensions. This structure holds for any diseasewhere the decay rate in the amplitude of successive epidemics is sufficientlyslower than the frequency of recurrent epidemics. That is, as the determin-istic process decays, there will be several noticeable epidemics before thesystem decays to near the steady state.After linearisation, we study the stochastic path in three-dimensionalspace. The sample path behaves as an OU process that travels along theaxis pointing in the direction of an eigenvector associated with the negativereal eigenvalue, and cycles on the subspace spanned by the eigenvectorsassociated with the complex eigenvalues (See Appendix A.4).We have shown that there is good agreement between the theoreticalPSDs of the exact and approximate processes for each population type con-sidered in the system. Although we observe small differences in the PSDs ofthe approximate and exact processes, we find that the two PSDs have closelymatching dominant frequencies. The dominant frequency of the simulatedepidemics is close to the intrinsic frequency ω. We know from e.g. Green-wood et al. [186] that the difference between the dominant and intrinsicfrequencies is O(λ), where λ, which is small here, represents the decay rateof the deterministic solution. Based on the polar form of the approximateprocess (3.28), each population cycles at a frequency corresponding to theintrinsic frequency ω perturbed by a stochastic phase process ϕλt.We also notice that the PSDs of the exact process appear flat for lowfrequency and have a downward slope for high frequency, which are featuresof the PSD of a stationary OU process [46]. This observation is in agreementwith the approximate process we derived (3.25), which is a sum of an OUprocess (3.26), and the product of a rotation and a bi-variate standard OUprocess.Previous studies [1, 13, 165] of the recurrence of avian flu epidemics havefocussed on identifying mechanisms that explain the multi-annual periodic-ity of the disease. Wang et al. [1] claimed that environmental transmissibilityis an important ingredient in explaining the 2 to 8 year period of avian flu.Using theoretical PSDs, they showed that for a fixed direct transmissionrate, the outbreak period decreases as the environmental transmission rate473.5. Discussionis increased. Their results are based on the assumption that direct trans-mission is weak. However, according to Roche and Lebarbenchon [175], thedirect transmission rate β can have a wide range of values. Our results showthat the 2-8 year outbreak period can also be obtained chiefly as a resultof direct transmission, and even in the absence of environmental transmis-sion. Thus we conclude that both transmission rates are important factorsin understanding the multi-year periodicity of disease outbreaks.Our approach also allows us to obtain approximate values for the typicalamplitude of the population fluctuations, because the stationary behaviourof the OU process is known. For the given avian flu parameter values, we findthat the virus population fluctuations peak after the infectious populationfluctuations peak. Typically, the phase lag is on the order of 0.3127/ω =0.31/0.8377 ≈ 0.37 of a year, i.e. approximately 4 months. This lag isreasonable as we have evidence that virus can be excreted by an infectedbird for many days after infection [97].High amplitude epidemics arise when the direct transmission rate is high.In this case, the stochastic perturbation phase process changes slowly andso the dominant frequency of epidemics is comparatively regular.Our analysis emphasizes that the interaction of stochasticity and trans-mission routes indeed plays an important role in determining outbreak pe-riodicity and intensity. We suggest here that the recurrent pattern of avianflu outbreaks in itself is a result of noise amplification wherein its peri-odicity and amplitudes are influenced by either or both of the modes oftransmission. The approach we introduced here could be used to perform asystematic study of other recurrent diseases.48Chapter 4Sustained Oscillations inStochastic ModelsWith Periodic ParametricForcing4.1 MotivationModellers have long been interested in identifying the mechanisms be-hind cyclic population fluctuations in natural systems such as those observedin infectious diseases [3, 187, 188]. In one theory, the fluctuations are at-tributed to the stochastic nature of the system. For example, in populationbiology, the number of individuals over time could fluctuate simply due tothe presence of demographic stochasticity (i.e. uncertainty in the birth anddeath times) together with stochasticity in infection and recovery rates. Inthe absence of any stochasticity, we have damped oscillations over a largeparameter range [32]. Sometimes the cyclic dynamics are modelled as adeterministic system with time-periodic parameter, i.e. parametric forc-ing, such as seasonal forcing driving recurrent epidemics [28, 39]. A rathermore general model of irregular cyclic dynamics has been introduced andsimulated quite recently. The model incorporates both randomness and pe-riodic forcing to describe the fluctuating patterns observed in systems withlarge populations [11, 33, 44] and expressed in terms of Itoˆ stochastic dif-ferential equations (SDEs). These SDEs can be constructed starting froman ODE model derived from compartmental analysis and adding Gaussianwhite noise with covariance matrix that depends on the terms of the ODEs.Periodic parametric forcing is introduced in the model by assigning at leastone of the parameters to be a periodic function of time. Here I coin thephrase stochastic periodically parametrically forced (stochastic PPF) mod-els to mean stochastic models with time-periodic parameters. Such modelsare SDE systems with time-dependent periodic drift and/or diffusion coef-494.1. Motivationficients.Stochastic PPF models is appear in studies of disease dynamics [33, 37,164] that are mainly concerned with characterizing the sustained oscilla-tions exhibited in stochastic simulations of the system. In stochastic PPFmodels, the stochastic fluctuations around a limit cycle are identified by“linearising” the stochastic PPF model around the solution of its determin-istic analogue. Though the nature of stochastic fluctuations around a limitpoint is well understood [26, 35, 50], limited knowledge is available aboutstochastic fluctuations around a limit cycle.Typically, solutions of linear SDE models are studied using numerical ap-proaches [51]. However, it is sometimes possible for approximate solutionsto be written in closed form. In the work of Baxendale and Greenwood[50], for example, an approximate solution was obtained for n-dimensional(n ≥ 2) linear SDE systems whose drift and diffusion coefficient matricesare constant matrices, and where the eigenvalues of the drift coefficient ma-trix include a complex-conjugate pair, −λ ± iω, with λ, ω > 0, and othernegative real eigenvalues. They showed that as long as λ  ω, such aSDE system has an approximate solution in closed form that describes thenature of noise-sustained oscillations as a constant times a circular motionmultiplied by a bivariate standard Ornstein-Uhlenbeck (OU) process [182]time-scaled by λ. They compared the power spectral densities (PSD) of thenumerical and approximate solutions of the SDE system and found agree-ment. The closed form approximation allows the identification of the sourceof the dominant frequency of the fluctuations and the parameters that con-trol amplification. Our aim here is to obtain such a closed form approximatesolution for stochastic PPF models linearized around a limit cycle.From numerical evidence found in the literature [11, 44], we know thatthe stochastic fluctuations around the limit cycle, produced by stochasticPPF models, have three governing frequencies: the frequency the around thelimit cycle in response to the stochastic perturbation, the sum and differenceof the response and limit cycle frequencies. The response frequency is thedominant one among the three and is identified as the imaginary part ofthe complex-conjugate pair of Floquet exponents [189]. Though it has beensuggested that these frequencies resulted from stochastic amplifications, ithas not been understood how they arise from the solutions of the model.We begin with a general formulation of discrete state-space two-dimensionalstochastic PPF models and their macroscopic limit, using the method ofKurtz [10], an approach for developing SDE models which parallels the VanKampen [9] system-size expansion. Assuming that a limit cycle solutionexists in the macroscopic limit, we show that the stochastic fluctuations can504.1. Motivationbe written as a linear diffusion equation with time-dependent drift and dif-fusion coefficients. We use the approach of Baxendale and Greenwood [50]to obtain an approximate solution of the linear stochastic model. First, were-formulate the problem using a series of transformations similar to thoseof Baxendale and Greenwood [50] but for time-varying coefficients. In oneof these transformations, we apply Floquet theory [52].We find that when the Floquet exponents are a complex-conjugate pairwith negative real part and the perturbation rotates more quickly than it de-cays, the stochastic PPF model can approximated by a standard OU processmultiplied by periodic matrices carrying the response frequency predictedby Floquet theory and the limit cycle frequency. Moreover, the process isamplified by the ratio of a time-averaged standard deviation from the time-varying co-variance matrix and the square-root of the magnitude of the realpart of the complex-valued Floquet exponent.We illustrate our approximation result using a linear harmonic oscillatormodel with one-dimensional noise. We compare the power spectral densities(PSDs) of the exact and approximate solutions computed numerically, forvarious forcing amplitudes and noise levels. By considering a simple caseof the model and using a fundamental matrix solution that contains thewell-studied Mathieu functions [190], we provide further insight into theapproximate solution. For instance, in Section 4.4, we compute stochasticdifferential equations for the amplitude and phase processes. In Section4.5, we apply to a stochastic epidemic model with seasonal forcing, in thecontext of a host-pathogen model for avian influenza with time-dependentenvironmental transmission.This chapter is organized as follows. In Section 2, we present the modeland its assumptions, and derive an expression for the linearisation aroundits limit cycle. This fluctuation process is a stochastic PPF model. Wederive in Section 3 a useful approximation to the stochastic PPF model ofSection 2. This explicit form of the solution of the system that allows us totease apart the effects of stochasticity and periodic forcing. In Section 4, weillustrate our approximation technique with a parametrically forced linearoscillator model with noise. In Section 5 is an example of the applicabilityof the approximation to an avian flu model with demographic stochasticityand seasonal forcing. Finally, further analysis and discussion of theoreti-cal results and the significance of the approximation are in Section 6 withemphasis on how the approach can shed light on fluctuations observed inepidemic models.514.2. A periodically parametrically forced (PPF) Markov model4.2 A periodically parametrically forced (PPF)Markov model4.2.1 DescriptionWe consider a continuous-time Z2-valued Markov process X(t) =(X1(t), X2(t))T with transition rates, or jump rates, f (r)(X(t), t) where r ∈R2 ranges over a set of possible jumps between pairs of integers. These ratesexplicitly depend on time t. Following Kurtz [10], we write the process X(t)asX(t) = X(0) +∑rrN (r)(t), (4.1)where N (r)(t) is the number of jumps of type r jumps within the timeinterval [0, t]. Here N (r)(t) is a counting process with rate f (r)(X(t), t) andN (r)(t)− ∫ t0 f (r)(X(s), s) ds is a martingale and this implies thatN (r)(t) = Ψ(r)(∫ t0f (r)(X(s), s) ds), (4.2)where Ψ(r)(u) are independent, unit Poisson processes with u =∫ t0 f(r)(X(s), s) ds =E[N (r)(t)] (Section 10.4 of Kurtz [53]).Suppose that X1(t) + X2(t) = N where N is constant for all t ≥ 0.Combining (4.2) to (4.1), we write the scaled process XN (t) ≡ X(t)N asXN (t) = XN (0) +∑rrNΨ(r)(∫ t0Nf (r)(XN (s), s) ds). (4.3)The integrand in (4.3) represents stochastic rates in density-dependentprocesses. In this thesis, the separate dependence on the time variable inthe jump rate will be a periodic parametric forcing.Using essentially the proof of Kurtz [10], we see that the diffusion ap-proximation of (4.3) can be written similarly asX˜N (t) = X˜N (0) +∑rrNB(r)(∫ t0Nf (r)(X˜N (s), s) ds), (4.4)where B(r)(t) = W (r)(t)+t and W (r)(t) are independent, standard Brownianmotion (Wiener process). Hence, (4.4) can be written as the Itoˆ equation524.2. A periodically parametrically forced (PPF) Markov modelX˜N (t) = X˜N (0) +∫ t0F(X˜N (s), s) ds+∑rr√N(∫ t0√f (r)(X˜N (s), s) dW (r)(s)),(4.5)where F(·) = ∑r rf (r)(·).The process (4.5), in more compact form, isdX˜N (t) = F(X˜N (t), t) dt+1√NC(X˜N (t), t) dW(t). (4.6)Here C(X˜N (t), t) is a 2×2 matrix that is obtained by taking the square-root of the covariance matrix function defined by [Vij ] =∑r rirjf(r)(X˜N (t), t).The vector F ∈ R2 represents the deterministic dynamics of the macroscopiclimit for (4.6) and W(t) is a two-dimensional Brownian motion.As N → ∞, the process X˜N (t) approaches a macroscopic limit φ(t)satisfying the two-dimensional ODEφ˙ = F(φ(t), t) ≡ F(φ(t); θ(t)). (4.7)In (4.7), the notation θ(t) appears in the re-written form of F(φ(t), t) toindicate that the vector-field F varies with t explicitly because a parameterthat periodically varies with time is introduced. Hence, we refer to (4.7) asperiodically parametrically forced (PPF) model.Using the same argument as Kurtz [10], by the Law of Large Numbersfor Ψ(r) under the assumption that F(·) is Lipschitz and as a consequenceof the central limit theorem, it follows that:XN (t) = φ(t) +O(1√N). (4.8)Since XN (t) ≈ X˜N (t), we haveX˜N (t) ≈ φ(t) +O(1√N)≡ φ(t) + 1√Nξ(t). (4.9)The process ξ(t) defined in (4.9) represents the random fluctuationsaround the macroscopic limit φ(t). Analysis of this process when φ(t) has astable fixed point is the common theme in studies on sustained oscillationsaround this fixed point observed in stochastic simulations of various popu-lation models [32, 33, 189]. In this chapter, we study the behaviour of ξ(t)when φ(t) has a stable limit cycle.534.2. A periodically parametrically forced (PPF) Markov model4.2.2 Stochastic fluctuations around the limit cycle solutionExisting work suggests that models of the form (4.7) may have a stablelimit cycle [33, 191]. In this section, we derive a stochastic differential equa-tion (SDE) for the process that describes the stochastic fluctuations arounda stable limit cycle where it exists. Let us assume that the system (4.7) hasa stable limit cycle solution L(t) with period TLC .For large t, keeping in mind that N is also large and using (4.9), we canreplace φ(t) by L(t) and obtain:X˜N (t) ≈ L(t) + 1√Nξ(t). (4.10)We substitute (4.10) into F and C in (4.6) and linearise them around thelimit cycle L(t) to obtainF(X˜N (t), t) ≈ F(L(t) +1√Nξ(t), t)≈ F(L(t), t) + 1√NDF(L(t), t) ξ(t),and C(X˜N (t), t) ≈ C(L(t), t),(4.11)whereDF(L(t), t) =[∂F1∂φ1∂F1∂φ2∂F2∂φ1∂F2∂φ2]∣∣∣∣∣φ(t)=L(t)≡ A(t), (4.12)i.e. the Jacobian of F with respect to the first variable evaluated at L(t). Thematrix (4.12) carries the forcing frequency embedded in the time-periodicforcing parameter θ(t) and the frequency of the limit cycle L(t).We integrate (4.6) and substitute (4.10) and (4.11) into X˜N (t), F andC, respectively to obtainL(t) +1√Nξ(t) ≈ L(0) + 1√Nξ(0) +∫ t0F(L(s), s) ds+1√N∫ t0DF(L(s), s) ξ(s) ds+1√N∫ t0C(L(s), s) dW(s).(4.13)544.2. A periodically parametrically forced (PPF) Markov modelSince L(t) satisfies (4.7),L(t)− L(0) =∫ t0F(L(s), s) ds, (4.14)and. writing C(L(s), s) = C(s), (4.13) simplifies to,ξ(t) = ξ(0) +∫ t0A(s) ξ(s) ds+∫ t0C(s) dW(s). (4.15)In differential form,dξ(t) = A(t)ξ(t) dt+ C(t) dW(t). (4.16)We now refer to (4.16) as the “linearisation” of the stochastic PPF model(4.6). It is an equation of motion for the fluctuations ξ(t), about the de-terministic trajectory L(t). In physics language, such an equation is knownas a linear Langevin equation, which can be obtained by means of a VanKampen [9] system-size expansion as in Black and McKane [191].Consider (4.16) in the absence of noise, i.e. when C(t) = 0. The remain-ing term describes the behaviour of the system when perturbed away fromthe limit cycle L(t). We call this behaviour the response of the system to adeterministic perturbation. We expect this response to exhibit damped os-cillations with a frequency that is not necessarily the same as the frequencyof the limit cycle. We call this new frequency the response frequency.Simulations of systems of the form (4.16) with A(t) given by (4.12)applied to disease and population dynamics produce sustained oscillations[11, 33, 44, 189]. The nature of the oscillations described by ξ(t) has beenstudied using the power spectral density (PSD), i.e. the distribution offrequencies in their stochastic simulations. The PSD of ξ(t) in specific ex-amples [189, 191] consistently shows a dominant peak at the imaginary partof complex Floquet exponents associated to the limit cycle and two smallerpeaks at frequencies equal to the sum and difference of the imaginary partof these Floquet exponents and the limit cycle frequency.The nature of the fluctuation in (4.16) has been analyzed using Floquettheory [52], which is useful when establishing the stability of a limit cycle.For example, Boland et al. [189] showed that stochastic noise gives rise tolarge-amplitude, coherent sustained oscillations for systems with stable limitcycles having complex Floquet multipliers. On the other hand, Black [11]developed an analytic tool based on Floquet quantities to predict the pe-riodicity in oscillations produced by a stochastic SIR model with seasonalforcing. However, it remains to show how the interaction of the stochastic554.3. The approximationnoise and parametric forcing produce large-amplitude oscillations that aresustained over time. Moreover, the differences between the characteristicsof the fluctuations arising from a stochastic PPF model and an unforcedstochastic model remains unclear. In the next section, we derive an approx-imate form of ξ(t) which makes apparent the source of the PSD frequenciesfound by simulation in Boland et al. [189], Black and McKane [191].4.3 The approximationThe solution ξ(t) of the SDE (4.16) fluctuates at a rate that is in-fluenced by parametric forcing and by the limit cycle embedded in A(t),and also varies according to Brownian noise with diffusion matrix C(t) that,again, depends on the forcing parameter θ(t).Clearly, it is not possible to separate the roles of stochasticity and ofparametric forcing by simply examining the numerical simulations of (4.16).Here we develop an approximation that allows us to separate the roles ofstochasticity and forcing in a stochastic PPF model whose deterministicanalogue has a limit cycle solution.4.3.1 TransformationsOur starting point is the system of SDEs defined by (4.16) whereφ(t) in (4.7) has a limit cycle. We apply a series of transformations. First,we re-write (4.16), for brevity, as:dξt = Atξt dt+ Ct dWt, (4.17)The subscript means that each variable is a function of t. Let T be theperiod of A(t).We can see from (4.17) that when Ct is the zero matrix, we obtain alinear ODE system with T -periodic coefficient matrix A(t). Let us considerthe systemx˙ = A(t)x, x(t) ∈ R2. (4.18)One can think of (4.18) as resulting from linearising (4.7) about thelimit cycle solution L(t). In other words, x(t) is a perturbation aroundL(t) whose dynamics is governed by (4.18). Now suppose that x1(t) andx2(t) are linearly independent vector solutions of (4.18), then the matrixΦ(t) = [x1 x2] is called a fundamental matrix solution of564.3. The approximationΦ˙ = A(t)Φ. (4.19)According to Floquet theory [52], if Φ(t) is a fundamental matrix so-lution for (4.18) so is Φ(t + T ) and there exists a constant matrix B thatsatisfies the relation Φ(t+T ) = Φ(t)B. The eigenvalues of B, which we de-note here as ρj for j = 1, 2, are called the Floquet multipliers. The Floquetexponents are µj = ln ρj/T . If their real parts are both negative, then thelimit cycle L(t) is stable. Since B = Φ(t)−1Φ(t+T ) is constant in t, a simpleway of expressing this matrix is by putting t = 0 so that B = Φ(0)−1Φ(T ).One can choose the principal fundamental matrix solution, i.e. Φ(0) = I, sothat the Floquet multipliers are the eigenvalues of B = Φ(T ).Given ρj and µj , according to Floquet theory, there exists a T -periodicfunction p(t) ∈ R2 such that x(t) = diag[eµ1t, eµ2t]·p(t). Thus, the principalfundamental matrix Φ0(t) for (4.18) has the formΦ0(t) = P0(t)Y0(t), (4.20)where the nonsingular 2×2 matrix P0(t) satisfies P0(t+T ) = P0(t), ∀t andY0(t) = diag[eµ1t eµ2t]. Moreover, one sees that Y0(t) satisfies the matrixdifferential equation with constant coefficients of the formY˙0 = D0Y0,with D0 = diag[µ1, µ2].Since the limit cycle L(t) is stable, then the solution of (4.18) has dampedoscillations, which implies that we have complex-valued Floquet exponentswith Re{µj} < 0. As a result, the solutions Φ0(t) are complex-valued. Weused these complex-valued Φ0(t) to construct linearly independent, real-valued solutions for (4.19).As in Grimshaw [52], if we suppose that µ and µ¯ are complex- conjugatepair, then µ = −λ+ iω (for λ, ω > 0) where |ρ| = e−λT and arg ρ = ωT , andp(t) = q(t) + ir(t), where q(t), r(t) ∈ R2 are T -periodic, then the real andimaginary parts of eµtp(t),Re{eµtp(t)} = e−λtq(t) cosωt− e−λtr(t) sinωt,Im{eµtp(t))} = e−λtq(t) sinωt+ e−λtr(t) cosωt, (4.21)form the real-valued fundamental matrix solution for (4.19). Hence, we canrewrite Φ0(t) as a real-valued solution574.3. The approximationΦ0(t) = [Re{eµtp(t)}, Im{eµtp(t)}],= [q(t) r(t)]︸ ︷︷ ︸Q(t)e−λt[cosωt sinωt− sinωt cosωt]︸ ︷︷ ︸Y(t). (4.22)We see from (4.22) that the matrix P0(t) in (4.20) can be written as theinvertible periodic matrix Q(t) = [q(t) r(t)] andY(t) ≡ e−λtR−ωt, (4.23)where Rθ is the rotation matrix. Note that Q(t) has period T . On the otherhand, Y(t) comes from the ODE (4.18) whose period is 2pi/ω not necessarilyequal to T . Additionally, Y(t) satisfies the equationY˙ = ΛY where Λ =[−λ ω−ω −λ].With this information about Φ0, the fundamental matrix solution of thedeterministic equation (4.19), we go back to the stochastic equation (4.17).Let us write ξ(t) = Q(t)y(t) in (4.17). Since y(t) is a stochastic process, wedescribe its dynamics using Itoˆ calculus. By Itoˆ’s formula, we havedy = Q˙−1Qy dt+ Q−1AQy dt+ Q−1C dW. (4.24)We wish to simplify the SDE (4.24) so that it is of the formdy = Λy dt + Q−1C dW. (4.25)Denote by [Q(t)]ij and [Q−1(t)]jk (for i, j, k = 1, 2) the components ofQ(t) and Q−1(t). For every t,2∑j=1[Q]ij [Q−1]jk = δik, (4.26)where δik is the Kronecker delta. The time derivative of (4.26) gives us2∑j=1[[Q]ij(ddt[Q−1]jk)+(ddt[Q]ij)[Q−1]jk]= 0, (4.27)which implies that584.3. The approximation2∑j=1[Q]ij(ddt[Q−1]jk)= −2∑j=1(ddt[Q]ij)[Q−1]jk. (4.28)A shorthand form of (4.28) is:[Q · Q˙−1]ik = −[Q˙ ·Q−1]ik. (4.29)From (4.29), we conclude that Q · Q˙−1 = −Q˙ ·Q−1 and so,Q˙−1 = −Q−1Q˙Q−1. (4.30)Substituting (4.30) into (4.24), we havedy = Q−1(AQ− Q˙)y dt+ Q−1C dW. (4.31)In (4.22), we see thatQ(t) = Φ0(t)Y−1(t) = Φ0(t)eλtRωt. (4.32)It follows thatQ˙ = Φ˙0eλtRωt + λΦ0eλtRωt + eλtΦ0R˙ωt. (4.33)Define the notation R˜ωt byR˙ωt = ω[− sinωt − cosωtcosωt − sinωt]≡ ωR˜ωt. (4.34)From (4.33), (4.34), (4.22), (4.19), (4.23), and simplifying, we obtainQ˙ = (A + λI)Q + ωeλtΦ0R˜. (4.35)Knowing from (4.32) thatΦ0(t) = e−λtQ(t)R−ωt (4.36)and writing ωRωtR˜ =[0 −ωω 0], (4.35) now readsQ˙ = AQ−QΛ. (4.37)Therefore, using (4.37), the SDE system (4.31) simplifies to (4.25) to whichthe proof of Baxendale and Greenwood [50] applies with some minor changes.Following Baxendale and Greenwood [50], we apply three changes ofvariables to (4.25). First, we write y(t) = R−ωtz(t) to simplify the drift594.3. The approximationterm. Again, by Itoˆ’s formula, the two-dimensional SDE for z(t) is found tobe:dz = −λz dt+ RωtQ−1C dW. (4.38)This means thatz(t) = z(0)− λ∫ t0z(s) ds+∫ t0RωsQ−1(s)C(s) dW(s). (4.39)Rescaling time by replacing t by t/λ, (4.39) then reads:z(t/λ) = z(0)− λ∫ t/λ0z(s) ds+∫ t/λ0RωsQ−1(s)C(s) dW(s). (4.40)We apply the substitution s = u/λ =⇒ ds = du/λ to (4.40) and obtainz(t/λ)− z(0) = −∫ t0z(u/λ) du+∫ t0Rωu/λQ−1(u/λ)C(u/λ) dW(u/λ).(4.41)Let C˜t = Q−1(t)C(t). Defineσ2(t) =12Tr(C˜tC˜∗t)(4.42)andσ¯2 =1T∫ T0σ2(t) dt. (4.43)By defining U(t) =√λz(t/λ)/σ¯, Equation (4.41) becomesU(t)−U(0) = −∫ t0U(u) du+∫ t0√λσ¯Rωu/λQ−1(u/λ)C(u/λ) dW(u/λ).(4.44)Using the scaling property of Brownian motion, i.e. W˜(t) =√λW(t/λ),and the notation D(t/λ) = 1σ¯Q−1(t/λ)C(t/λ), (4.44) takes the differentialformdU(t) = −U(t) dt+ Rωt/λD(u/λ) dW˜(t). (4.45)We note that (4.45) takes a similar form as with the SDE found in Eq. (15)in Baxendale and Greenwood [50], which was compared to the bivariatestandard Ornstein-Uhlenbeck process.Using the Martingale Problem method, Baxendale and Greenwood [50]proved that, on any fixed bounded time interval [0, τ ], the distribution ofU˜(t) satisfying the equationdU˜(t) = −U˜(t) dt+ Rωt/λD dW˜(t), where Tr (DD∗) = 2 (4.46)604.3. The approximationconverges to the distribution of the bivariate standard Ornstein-Uhlenbeckprocess S(t) generated by the SDEdS(t) = −S(t) dt+ dW(t) (4.47)as λ/ω → 0. In their context, D is a constant matrix.We check that in our context where D is a function of t, limλ→0 Tr(DD∗) =2 holds for (4.45). Observe thatTr(Dt/λD∗t/λ) =1σ¯2Tr(C˜t/λC˜∗t/λ)=2σ2(t/λ)σ¯2. (4.48)Since σ2(t/λ) is periodic, its frequency increases as λ → 0. This impliesthat in the limit, σ2(t/λ) can be approximated by its time-average over theinterval [0, T ], i.e. σ¯2. Hence, the Tr(Dt/λD∗t/λ) → 2 as λ → 0, and thetheorem of Baxendale and Greenwood [50] applies to (4.45). We concludethat the process U(t) converges weakly to the process S(t).4.3.2 ResultThe approximate processWe conclude that the solution ξ(t) of (4.17) when λ ω can be approx-imated as follows:ξ(t) = Q(t)y(t) = Q(t)R−ωtz(t) =σ¯√λQ(t)R−ωtUλt≈ σ¯√λQ(t)R−ωtSλt ≡ ξapp(t).(4.49)We see from (4.49) that the stochastic fluctuations ξ(t) around the limitcycle oscillate at frequencies which are combinations of the frequency ofthe drift coefficient A(t) and the frequency ω, which arise already from thedeterministic system (4.18). Note that the frequency ω is the imaginarypart of the Floquet exponents associated with the deterministic portion of(4.17). The typical amplitude of the fluctuation ξ(t) depends on the factorσ¯/√λ, where σ¯ is square-root of the average variance σ2(t) over the timeinterval [0, T ] and λ is the magnitude of the real part of the complex Floquetexponents.614.3. The approximationThe phase and amplitude processesWe write R−ωtSλt in polar form to see the role of the two-dimensionalstandard Ornstein-Uhlenbeck process Sλt appearing in (4.49) and obtainξapp(t) =σ¯√λQ(t)Zλt[cos(ωt− ϕλt)sin(ωt− ϕλt)], (4.50)where Zλt = |Sλt| and ϕλt = arg[S1(λt) + iS2(λt)]. In (4.50), we have awell-defined phase process ϕ(t) and a radial process Z(t) which are knownto satisfy the following stochastic diffential equations [46]:dZ(t) =[12Z(t)− Z(t)]dt+ dWZ(t),dϕ(t) =1Z(t)dWϕ(t),(4.51)where WZ(t) and Wϕ(t) are independent Wiener processes: WZ(t) is a stan-dard Brownian motion and Wϕ(t) is a Brownian motion run on a unit circle.Note here that in the context of (4.50), Z(t) and ϕ(t) in (4.51) should berun at a rate λ. Hence, the amplitude of ξapp(t) does not depend only onthe scalar factor and the magnitude of the T -periodic matrix Q0(t) but isalso modulated by the two-dimensional slowly varying standard Ornstein-Uhlenbeck process S(λt). Furthermore, the process Wϕ(t) is responsible forthe random phase slips (instantaneous change in the phase of an oscilla-tor) affecting the deterministic phase of the cycles produced by Q(t)R−ωt.From (4.51), we observe that when the radial process Z(t) is large for arandom time interval, the phase process moves at a slow rate and so im-plies that the phase process of ξ(t) predominantly follows Q(t)R−ωt in anapproximate sense. In this case, the radial process behaves much like anOrnstein-Uhlenbeck process. On the other hand, the increments of ϕ(t)become prominent when Z(t) is small.The form of Q(t) and its consequenceAt this point, we only know that Q(t) ∈ R2×2 is T -periodic, i.e. Q(t)has the same period as A(t). When Q(t) have function entries that belongto a class of smooth functions whose Fourier series coefficients decay rapidlywith increasing frequency index [192], then we can writeQ(t) ∼ Q0 + Q1 cos νt+ Q2 sin νt, (4.52)624.3. The approximationwhere Q0,Q1, and Q2 are constant matrices of Fourier coefficients with[Q0]ij > [Q1]ij > [Q2]ij for i, j = 1, 2 and ν = 2pi/T . Thus,Q(t)R−ωt ∼ Q0R−ωt + Q1R−ωt cos νt+ Q2R−ωt sin νt. (4.53)Using product-to-sum trigonometric identities, we haveR−ωt cos νt =12[cosψ+t + cosψ−t sinψ+t − sinψ−tsinψ−t − sinψ+t cosψ+t + cosψ−t](4.54)andR−ωt sin νt =12[sinψ+t + sinψ−t cosψ−t − cosψ+tcosψ+t − cosψ−t sinψ+t + sinψ−t], (4.55)where ψ±t = (ν±ω)t. Therefore, the entry of Q(t)R−ωt is a sum of sine andcosine functions with frequencies ω and ν ± ω.Given the form of Q(t) in (4.52), we letQk =[Q(1,1)k Q(1,2)kQ(2,1)k Q(2,2)k], k = 0, 1, 2 (4.56)and rewrite (4.50) asξapp(t) =σ¯√λZλt[r(1)0 cos Ψ01 +12r(1)1(cos Ψ+11 + cos Ψ−11)+ 12r(1)2(sin Ψ+21 + sin Ψ−21)r(2)0 cos Ψ02 +12r(2)1(cos Ψ+12 + cos Ψ−12)+ 12r(2)2(sin Ψ+22 + sin Ψ−22)] ,(4.57)wherer(j)k =√(Q(j,1)k)2+(Q(j,2)k)2, (4.58)θ(j)k = arctan(Q(j,2)kQ(j,1)k), (4.59)Ψ0j = ωt− ϕλt − θj0, (4.60)andΨ±ij = νt± ωt∓ ϕλt ∓ θ(j)i . (4.61)In (4.57), with (4.58)-(4.61), we find that the approximate process ξapp(t)representing the stochastic perturbation around the limit cycle is a linearcombination of sinusoidal functions whose central frequencies are ω and ν±634.4. Example: Driven harmonic oscillator with 1D noiseω, which are determined from the deterministic matrix function Q(t)R−ωt.Based on (4.60) and (4.61), each central frequency is perturbed by the slowlyvarying stochastic phase process ϕλt and has θ(j)k , given by (4.59), as phaseshifts. Moreover, the amplitude of ξapp(t) in (4.57) is modulated by theslowly varying amplitude process Zλt (as in (4.50)) and depends on theconstants r(j)k , which is given by (4.58). These constants are associated tosinusoidal waves that determine the rotational dynamics of the process ξ(t).From (4.58), we find that r(j)k is obtained from a sequence of decaying Fouriercoefficients and so we expect that the process ξ(t) predominantly rotates ata central frequency ω, in an approximate sense.In the following section, we will illustrate that ξ(t) ≈ ξapp(t) in an ex-ample by verifying that their numerical power spectral densities (averagedover many realizations) are in good agreement.4.4 Example: Driven harmonic oscillator with1D noiseConsider the stochastic PPF model:dξ(t) = A(t)ξ(t) dt+ C(t) dW(t), where (4.62)A(t) =[0 1−k0 −b0(1 + b cos 2piat)]and C(t) =[0 00 σ0(1 + b cos 2piat)].(4.63)Model (4.62) is linear in ξ(t) = [ξ1, ξ2]T and if we let x = ξ1, we can writeit asx¨+ k0x˙+ b0(1 + b cos 2piat)x = σ0(1 + b cos 2piat)W˙ . (4.64)Equation (4.64) describes a harmonic oscillator which is driven by peri-odic variation of one of its parameters at frequency a and is further affectedby a white noise whose variance depends on time and other parameters inthe model. This white noise acts as an external force applied to the oscil-lator. In the model, the parameters affecting the frequency of the oscillatorare the forcing frequency (a) and amplitude (b).In the absence of parametric forcing and noise, i.e. b = 0, σ0 = 0, (4.64)is the well-known damped harmonic oscillator model. For the harmonicoscillator system, the coefficient k0 represents the damping rate and ω0 =644.4. Example: Driven harmonic oscillator with 1D noise√b0 is the oscillator’s angular frequency when it is not damped. A nonzero k0implies that the harmonic oscillator is damped i.e., the oscillator’s amplitudegradually goes to zero. In this case, the frequency of the oscillator does notonly depend on ω0 but also to the damping ratio ζ = k0/2√b0. In particular,the frequency of a damped harmonic oscillator is given by ω1 = ω0√1− ζ2.When k0 is zero, we recover the simple harmonic oscillator whose frequencyis ω0.4.4.1 Power spectral density (PSD)One aim of this work is to determine that the approximation is a gooddescription for the stochastic oscillations produced by (4.62) as representedby the PSD for various values of the forcing amplitude b and noise levelσ0. Two main parameter regimes are considered here, namely when thedeterministic system is damped (k0 6= 0) or undamped (k0 = 0). We choosethese regimes because both can give rise to noisy but sustained oscillations.We obtained 100 realizations of stochastic simulations and computed theaverage power spectral density (PSD) numerically. We compare the PSDoutputs for different regimes and various levels of b and σ0 (See Figure 4.1).We choose the forcing frequency to be a = 1. We compare these with thePSD of the stochastic simulations (on average), based on the approxima-tion (4.49). In Figure 4.1, there is a good agreement in the PSDs of theapproximate and exact solutions with respect to capturing the governingfrequencies of the process. The approximate process tends to give higherPSD values than the exact process.In the damped regime, the deterministic system displays decaying os-cillations approaching to zero solution while in the undamped regime, thesystem oscillates regularly at a frequency ω0 without damping. In bothregimes, the process ξ(t) have large PSD peaks at the system’s intrinsic fre-quency ω0 indicating that the process’ dominant frequency is amplified bystochastic and parametric forcing. The undamped regime however gives riseto higher PSD values than the PSD values in the damped regime. Otherobserved PSD structure is the broadening of PSD peaks which results fromthe stochastic phase perturbations explained by the form of our approxima-tion (4.49). The flattening of the PSDs in the right end indicates that theprocess is comparable to the Ornstein-Uhlenbeck process.654.4. Example: Driven harmonic oscillator with 1D noise0 ω/2pi 1−ω/2pi 1 1+ω/2pi100105Frequency (time−1)Log PowerPSD of ξ1: Underdamped, Small Forcing  exact,σ0=1exact,σ0=5approx, σ0=1approx, σ0=5(a)0 ω/2pi 1−ω/2pi 1 1+ω/2pi100105Frequency (time−1)Log PowerPSD of ξ1: Underdamped, Large Forcing  exact,σ0=1exact,σ0=5approx, σ0=1approx, σ0=5(b)0 ω/2pi 1−ω/2pi 1 1+ω/2pi100105Frequency (time−1)Log PowerPSD of ξ1: Undamped, Small Forcing  exact,σ0=1exact,σ0=5approx, σ0=1approx, σ0=5(c)0 ω/2pi 1−ω/2pi 1 1+ω/2pi100105Frequency (time−1)Log PowerPSD of ξ1: Undamped, Large Forcing  exact,σ0=1exact,σ0=5approx, σ0=1approx, σ0=5(d)Figure 4.1: Numerical PSDs of the approximate (solid curves) and exact(dashed curves) solutions of (4.62), given (a,b) damped and (c,d) undampedregimes as defined in Section 4.4.1, obtained using forcing amplitudes (a,b)b = 0.01 and (c,d) b = 0.5 for noise levels σ0 = 1 (gray curves) and σ0 = 5(black curves). Other parameters are held fixed, namely a = 1, b0 = 2,k0 =0.1 for underdamped case, and k0 = 0 for the undamped case.664.4. Example: Driven harmonic oscillator with 1D noise4.4.2 Amplitude and phase as stochastic processesWe now apply the approximation to mathematically describe the ampli-tude and phase of the stochastic fluctuations produced by (4.62),(4.63) interms of the model parameters. We use the fundamental matrix solutionΦ0(t) of the systemx˙(t) =[0 1−k0 −b0(1 + b cos 2piat)]x(t) (4.65)and write the approximation asξapp(t) =σ¯√λeλtΦ0(t)Sλt. (4.66)Here (4.66) is obtained by re-writing ξapp(t) in (4.49) using (4.36). Recallthat Φ0(t) is the principal fundamental matrix solution that satisfies (4.19).We compute the fundamental matrix solution of (4.65) using DSolvecommand in Mathematica [193] and obtainΦ(t) = e−12k0t[C (α, q, x) S (α, q, x)−k02 C (α, q, x) + apiC ′(α, q, x) −k02 S (α, q, x) + apiS ′(α, q, x)],(4.67)where α =4b0 − k204a2pi2, q =bb02a2pi2, and x = apit. The functions C and S areknown as Mathieu cosine and sine, respectively, which are unique solutionsto the Mathieu equation [194, 195] given byd2ydx2+ [a− 2q cos (2x)]y = 0.The functions C ′ and S ′ are derivatives with respect to x.For simplicity, we consider the case when the forcing amplitude is small,i.e. b ≈ 0, so that q is small. In theory, when q is small, we have C (α, q, x) ≈cos (√αx) and S (α, q, x) ≈ sin (√αx) /√α [190]. Hence, assuming that theforcing amplitude b is small with all other parameters fixed, it follows thatq  1,C (α, q, x) ≈ cos(12t√4b0 − k20), andS (α, q, x) ≈ 2api√4b0 − k20sin(12t√4b0 − k20).674.4. Example: Driven harmonic oscillator with 1D noiseDefine θ = 12√4b0 − k20, then C ′ = −θ sin(θt)/api and S ′ = cos(θt). There-fore, for small forcing amplitude, the fundamental matrix solution of (4.65)readsΦ(t) = e−12k0t[cos(θt) apiθ sin(θt)−k02 cos(θt)− θ sin(θt) −apik02θ sin(θt) + api cos(θt)].(4.68)We hence construct the principal fundamental matrix solution of (4.65)using Φ0(t) = Φ(t)Φ−1(0) and obtainΦ0(t) = e− 12k0t[cos(θt) + k02θ sin(θt)1θ sin(θt)(−θ − k204θ)sin(θt) −k02θ sin(θt) + cos(θt)]. (4.69)The real and imaginary part of the complex Floquet exponents −λ ± iω,computed by taking the eigenvalues of Φ0(T ), are λ = − cos θT and ω =| sin θT |, respectively. To guarantee that λ > 0 and λ  ω, we must havethat pi/2T < θ < 3pi/4T or 5pi/4T < θ ≤ 3pi/2T .Therefore, for appropriate θ and by (4.66), the approximate solution for(4.62) when b 1 is given byξ(t) ≈ ξapp(t)=σ¯√λe(λ−k0/2)t ×[cos(θt) + k02θ sin(θt)1θ sin(θt)(−θ − k204θ)sin(θt) −k02θ sin(θt) + cos(θt)] [S1(λt)S2(λt)].(4.70)In (4.70), we have ξapp(t) expressed as an explicit SDE system whichhas a unique solution. We write (4.70) in terms of phase-shifted sine andcosine to characterise the amplitude and phase of the approximate process.Let S1(λt) = Zˆλt cos ϕˆλt and S2(λt) = θZˆλt sin ϕˆλt where Zˆ(λt) = |S1(λt) +i1θS2(λt)| and ϕˆ(λt) = arg{S1(λt) + i1θS2(λt)}. By expanding (4.70) andusing formulas for products of sine and cosine functions, we haveξapp(t) =σ¯√λe(λ− k02)tZˆλt ×[cos Θ−t +k04θ sin Θ+t +k04θ sin Θ−t−k04 cos Θ−t + k04 cos Θ+t −k208θ sin Θ+t + (θ − k208θ ) sin Θ−],(4.71)where Θ±t = θt± ϕˆλt.684.4. Example: Driven harmonic oscillator with 1D noiseThe amplitude of ξapp(t) varies stochastically in a slow manner due tothe radial process Zˆλt. On the other hand, the stochastic process rotates ata dominant angular frequency θ and has a phase that varies according tothe stochastic perturbation given by ϕˆλt.We can take a further step in our analysis by obtaining the stochasticdynamics of the radial process Zˆ(t) and phase process ϕˆ(t). To do this, weconsiderlog(S1(t) + i1θS2(t))= log Zˆ(t) + iϕˆ(t). (4.72)We then apply the Itoˆ formula to obtaind(log Zˆ + iϕˆ)=d(S1 + iS2/θ)S1 + iS2/θ− 12[d(S1 + iS2/θ)]2(S1 + iS2/θ)2. (4.73)By expanding the dSi terms in (4.73) using (4.47) and noting that dW1(t)dW2(t) =0, dW 2i (t) = dt, and dWi(t)dt → 0 faster than dW 2i (t) as dt → 0, we findthat (4.73) simplifies tod(log Zˆ + iϕˆ)= −(S1 + iS2/θ)S1 + iS2/θdt+dW1 + idW2/θS1 + iS2/θ− (1− 1/θ)2(S1 + iS2/θ)2dt.(4.74)We write S1 + iS2/θ = Zˆeiϕˆ and group together appropriate terms in (4.74)so that we arrive atd(log Zˆ + iϕˆ)=(−1−(θ − 12θ)cos ϕˆ+ i sin ϕˆZˆ2)dt+(dW1 + idW2/θ)(cos ϕˆ+ i sin ϕˆ)Zˆ.(4.75)We now take the real part of (4.75) and apply the Itoˆ formula to ˆZ(t) =elog Zˆ(t) and obtain the SDE for Zˆ(t), which isdZˆ(t) =(1− θ2θZˆ(t)cos ϕˆ(t)− Zˆ(t))dt+ dWZˆ(t)+(12Zˆ(t)cos2 ϕˆ(t) +12Zˆ(t)θ2sin2 ϕˆ(t))dt.(4.76)The imaginary part of (4.75) gives the SDE for the phase process ϕˆ(t), whichisdϕˆ(t) =1− θ2θZˆ2(t)sin ϕˆ dt+1Zˆ(t)dWϕˆ(t). (4.77)694.5. Application: an avian flu modelThe scalar Wiener processes dWZˆ and dWϕˆ are independent of eachother. As a check, we can set θ = 1 and find that the radial process Zˆ(t) in(4.76) and phase process ϕˆ(t) in (4.77) satisfy (4.51). In this example, wecan see from (4.76) and (4.77) that the angular frequency θ, which dependson the damping rate k0 and the parameter b0, plays an important role indescribing the stochastic dynamics of Zˆ(t) and ϕˆ(t). For instance whenθ  1, we observe that in a random time interval, the phase ϕˆ(t) tends tobe more influenced by its drift term whenever Zˆ(t) is small. This behaviouris not observed in the phase process described in (4.51).4.5 Application: an avian flu modelStochastic epidemic models have been used to capture the qualitativepattern of recurrent epidemics. Some of these stochastic models incorporateboth demographic stochasticity and periodic parametric forcing [33, 37, 189].These stochastic PPF model account for the seasonal but stochastic natureof the epidemics. In this section, we describe a stochastic model for avianflu that takes into account seasonal environmental transmission and showhow our analytic tool can be applied.4.5.1 Stochastic avian flu model with seasonal forcingA simple stochastic avian flu model was recently developed by Wanget al. [1] and further by Mata et al. [196]. There are three stochastic pro-cesses representing the susceptible and infected host populations, and thevirus concentration in the environment. In the model, the disease is as-sumed to be transmitted via direct contact by infected host to a suscepti-ble individual (direct transmission) and/or via ingestion of virus particlesfrom a contaminated aquatic habitat by a susceptible host (environmentaltransmission). The model reveals that the dominant outbreak period of thedisease varies with environmental transmission. Motivated by this resultand given that there is biological evidence supporting seasonal variation inenvironmental transmission, e.g. virus persistence in cold water [93, 169], weexplore here the dynamics of the stochastic avian flu model in the presenceof seasonal forcing. The system of SDEs for the unforced stochastic avian flumodel was originally derived using a Van Kampen [9] system-size expansion[1] and has been developed alternatively using Kurtz [10] method [196]. Weobtain a seasonally forced stochastic model for avian influenza by replacingthe constant environmental transmission rate parameter ρ, in the unforced704.5. Application: an avian flu modelsystem, with a time-periodic sinusoidal function (i.e. ρ(t) = ρ0(1+b cos 2pit))given byds = (−βsi− ρ(t)sv + µ(1− s)) dt+ 1√N(−G1 dW1 +G2 dW2 +G3 dW3) ,di = (βsi+ ρ(t)sv − (µ+ γ)i) dt+ 1√N(G1 dW1 −G3 dW3 −G4 dW4) ,dv = (kτi+ δv − ηv) dt+ 1√NV(G5 dW5 −G6 dW6) ,(4.78)whereG1 =√βsi+ ρ(t)sv, G2 =√µ(1− s− i), G3 =√µi, G4 =√γi,G5 =√kτi+ δv, and G6 =√ηv.(4.79)The stochastic processes s(t), i(t), and v(t) represent proportions of suscepti-ble, infected, and virus populations, respectively. Infection occurs by directtransmission with a constant rate β and by environmental transmission witha time-periodic rate ρ(t). The environmental transmission rate varies arounda reference rate ρ0 in a sinusoidal manner whose intensity is represented bya forcing amplitude b. Note that the variation is expressed in terms of acosine function with one-year period. Viruses are more persistent and soare likely to be transmitted during the winter season, i.e. the beginning andend of the year. For simplicity, a host is assumed to be born and to die atthe same constant rate µ. Infected individuals recover at the rate γ. Thevirus population in the environment increases with the shedding rate τ ofinfected ducks and produced from an alternate population of infected birdsat a rate δ. The population of the virus decreases with the clearance rate η.The host population is also assumed to be constant with size N while theviral reference concentration is represented by a constant parameter NV . Inthe model (4.78), the scaling parameter k = N/NV resulted from the deriva-tion of the SDE, and the Wiener process Wi(t) whose standard deviation isGj . The Wiener processes are independent of each other. Note that G1 istime-periodic since it depends on ρ(t). Using equivalences among Brownianmotions [48], (4.78) can be written asdx(t) = F(x(t), t) dt+ DC(x(t)), t) dW(t), (4.80)714.5. Application: an avian flu modelwhere x(t) = (s(t), i(t), v(t)), D = diag(1√N, 1√N, 1√NV), W(t) ∈ R3 repre-sent independent Wiener processes,F(x, t) =−βsi− ρ(t)sv + µ(1− s)βsi+ ρ(t)sv − (µ+ γ)ikτi+ δv − ηvandC(x, t) =C11 C12 0C21 C22 00 0 C331/2 , (4.81)for C11 = βsi+ ρ(t)sv + µ(1− s), C12 = C21 = −βsi− ρ(t)sv − µi,C22 = βsi+ ρsv + (µ+ γ)i, and C33 = kτi+ δv + ηv.4.5.2 Deterministic dynamicsAs N,NV →∞, (4.78) gives rise to the macroscopic limitφ˙1 = −βφ1φ2 − ρ(t)φ1ψ + µ(1− φ1),φ˙2 = βφ1φ2 + ρ(t)φ1ψ − (µ+ γ)φ2,ψ˙ = κτφ2 + δψ − ηψ.(4.82)The avian flu SIR-V deterministic system (4.82) exhibits dynamics thatis similar to an SIR model with time-dependent forcing studied by Black andMcKane [191]. In the absence of forcing, i.e. b = 0, the avian flu system,using appropriate avian flu parameter values (Table 4.1), exhibits dampedoscillations. On the other hand, when forcing is present (e.g. b = 1), astable limit cycle is observed (See Figure 4.2).The existence of the stable limit cycle generated by (4.82), induced byseasonal forcing, implies that the deterministic system (4.82) displays regularoscillations for large t. In fact, the amplitude of these oscillations increasesas the intensity of forcing is increased (Figure 4.3). The frequency of thelimit cycle is, in this example, equal to the forcing frequency a = 1.724.5. Application: an avian flu model0.60.030.650.70.33ψ0.0250.75φ20.80.32φ10.850.02 0.310.015 0.3(a) (b)Figure 4.2: Deterministic dynamics when (a) b = 0 and (b) b = 1. Theinitial condition is indicated by a black arrow. The initial conditions usedare: (φ1, φ2, ψ) = (0.3295, 0.0298, 0.6935).Figure 4.3: Limit cycle generated from the solution of (4.82) for differentvalues of the forcing amplitude b.734.5. Application: an avian flu modelTable4.1:Descriptionsoftheparametersintheavianflumodelandthevaluesusedinstochasticsimulationsof(4.78).ParameterValue/RangeDescriptionUnitN103DuckpopulationsizeduckNV105Virusconcentrationintheenvironmentvirion/mLµ0.3Duckbirthanddeathrateyear−1δ0.1Virusreplicationrateyear−1τ104VirussheddingratevirionmL−1/duck/yearη3Virusclearancerateyear−1β15Directtransmissionrateduck−1year−1ρ00.5Environmentaltransmissionrateyear−1γ10Duckrecoveryrateyear−1a1Forcingfrequencyyear−1b1Forcingamplitude744.5. Application: an avian flu model4.5.3 Stochastic fluctuation around the limit cycle and itsapproximationWe are interested in the stochastic fluctuation around the limit cyclesolution associated with the deterministic analogue of (4.78). Suppose that(φ¯1, φ¯2, ψ¯) denotes the limit cycle solution of (4.78). The stochastic fluctu-ations around (φ¯1, φ¯2, ψ¯) can be described as in Section 4.2 by the systemof SDEs given bydξ(t) = A(t)ξ(t) dt+ C(t) dW(t), (4.83)whereA(t) =−βφ¯2 − ρ(t)ψ¯ − µ −βφ¯1 −ρ(t)φ¯1−βφ¯2 βφ¯1 − µ− γ ρ(t)φ¯10 κτ δ − η , (4.84)andC(t) =V11 V12 0V21 V22 00 0 V331/2 , whereV11 = βφ¯1φ¯2 + ρ(t)φ¯1ψ¯ + µ(1− φ¯1),V12 = V21 = −βφ¯1φ¯2 − ρ(t)φ¯1ψ¯ − µφ¯2,V22 = βφ¯1φ¯2 + ρ(t)φ¯1ψ¯ + (µ+ γ)φ¯2, andV33 = κτφ¯2 + δψ¯ + ηψ¯.(4.85)Here ξ and W take values in R3, and A,C take values in R3×3. Equations(4.83)-(4.85) are obtained by linearising (4.80) around xLC(t) = (φ¯1, φ¯2, ψ¯)using the substitutionx(t) = xLC(t) + Dξ(t),where D = diag(N−1/2, N−1/2, N−1/2V). Using avian flu parameter valuesdisplayed in Table 4.1, we find that the drift and diffusion matrices areperiodic with period T = 1. We illustrate this point by plotting the valuesof A11(t) and C11(t) in Figure 4.4.754.5. Application: an avian flu model0 0.5 1 1.5 2 2.5 3time (year)-1.4-1.3-1.2-1.1-1-0.9-0.8-0.7-0.6-0.5-0.4valueA11(t)period = 1(a)0 0.5 1 1.5 2 2.5 3time (year)0.450.50.550.60.650.70.75valueC11(t)period=1(b)Figure 4.4: Computed values of (a) A11(t) and (b) C11(t) in (4.83)-(4.85).The parameter values used are in Table 4.1. Each function has a periodicbehaviour with period equal to 1, the period of the limit cycle and the forcingfrequency.Since A(t) in (4.83) is periodic for avian flu parameter values, we useFloquet theory to determine the quantities needed to confirm that (4.83)can be approximated using (4.49).The Floquet matrix associated to the deterministic analogue of (4.83) isB =−0.0787 −1.3454 −0.08840.1208 0.0619 0.00432.6447 3.8381 0.2582 . (4.86)We compute the ratio of the natural logarithm of the eigenvalues of theFloquet matrix B and the period of A(t) to obtain the Floquet exponentswhich are −0.4967 ± 1.3714i and −8.1455. These quantities imply that adeterministic perturbation in three-dimensional space damps toward a limitcycle on a plane. The subspace where the limit cycle lies is spanned bythe eigenvectors of the Floquet matrix B corresponding to the complexeigenvalues.Now, we observe that the ratio of the decay rate and cycle frequency is0.4969/1.3715 = 0.3622, which is significantly less than 1. Hence, one canuse (4.49) to make an approximate description of the stochastic process ξ(t)satisfying (4.83).Though we have limited tools to obtain an explicit form of the matrixQ(t) in terms of avian flu parameters, it is still possible to describe the764.6. Discussionapproximate dynamics of the stochastic fluctuations ξ(t) for the stochasticavian flu model with seasonal forcing. In Mata et al. [196], we know thatwithout seasonal forcing, the stochastic fluctuation is drawn toward a planespanned by the eigenvectors associated by the complex eigenvalues of theconstant drift coefficient matrix. The fluctuations chiefly lie on this planeand its dynamics follow a cyclic path modulated by a bi-variate standardOU process. We expect the stochastic fluctuations in the seasonally forcedcase, i.e. (4.83), to behave in a similar manner with slight modifications.First, note that since the fluctuations are stochastic perturbations aroundthe limit cycle, the subspace where the fluctuations will be attracted tois equivalent to a plane where the limit cycle lies. Second, the governingfrequencies of the stochastic fluctuations are determined by the periodicityof Q(t) and the angular frequency identified using Floquet exponent (e.g.1.3715). Finally, the amplitude of the fluctuations around the limit cycledepends on the amplitude of the entries of Q(t), the decay rate of the deter-ministic perturbation (e.g. 0.4969), and a time-averaged standard deviation.Thus, based on approximation (4.49) and the previous analysis, the stochas-tic population fluctuations in our avian flu system (4.78) satisfying (4.83)must be a sum of a scaled one-dimensional OU process, and a product ofdeterministic and a stochastic (i.e. slowly-varying bi-variate OU) processes.The deterministic processes involve the interactions of (i) periodic trans-formation (Q(t)), (ii) rotation (R−ωt), and ratio of time-averaged standarddeviation (σ¯) and square-root of the decay rate (√λ) of the deterministicperturbation.4.6 DiscussionStochastic models with periodic parametric forcing are seen to produce alimit cycle with noise in stochastic simulations [11, 26, 33]. We have shownthat the stochastic fluctuation around the limit cycle can be approximatedby a product of a scalar σ¯/√λ, a periodic matrix Q(t), a rotation matrixR−ωt, and a slowly varying bivariate standard Ornstein-Uhlenbeck (OU)process Sλt. The parameters λ, ω, and the matrix Q(t) arose from Floquettheory applied to the solution of the deterministic analogue of the linearSDEs with time-periodic coefficients. In particular, λ and ω are magnitudesof the real and imaginary part of the Floquet exponents while Q(t) dependson the forcing frequency a, which determines the periodicity of the deter-ministic solution. In numerical PSDs generated by a model examples, i.e.,the driven harmonic oscillator model with one-dimensional noise and the774.6. Discussionavian flu model (not shown), we see that the approximate process describesthe exact process well in terms of the peak frequencies observed. The PSDoutputs show that stochasticity and parametric forcing result in a processwhose dominant frequency is close to the intrinsic Floquet frequency of thedeterministic analogue.Linear SDE models as in (4.17) arise from systems that incorporateparametric forcing and stochastic noise. Though in our example we onlyassume one forcing parameter that is time-periodic, the approximation resultwe derived here also applies to the case where there is more than one time-periodic forcing parameter. The resulting approximate process in such acase would have a complicated structure for the periodic matrix Q(t).Our work differs from that of Baxendale and Greenwood [50] in that weconsider a diffusion process described by a linear SDE system with time-periodic drift and diffusion coefficients. The mathematical structure of ourapproximate process (4.49) is very similar to the approximate process inBaxendale and Greenwood [50] but with a different interpretation of its fac-tors. In this study, the approximate process is a perturbation process, dueto stochastic noise and parametric forcing, around a stable limit cycle. InBaxendale and Greenwood [50], however, the approximate process is a per-turbation process around a stable limit point. As a result, our approximateprocess contains frequencies that depend on the limit cycle frequency, asidefrom the intrinsic frequency of the deterministic analogue. Note that forboth types of approximate process, the description of sustained oscillationscontains a slowly varying bivariate standard OU process.In our examples, we found that the interaction of stochasticity and para-metric forcing can make hidden frequencies visible. When the intensity offorcing is high and the stochastic noise level is sufficient, the process de-scribed by (4.17) will have three main frequencies: (1) an intrinsic frequencyω/2pi, (2) the sum of the intrinsic frequency and the limit cycle frequency(e.g. 1 + ω/2pi), and (3) the difference of the limit cycle frequency andthe intrinsic frequency (e.g. 1− ω/2pi). According to the approximate form(4.49), we see that these frequencies all come from the product of Q(t)R−ωt,which is related to the fundamental matrix solution associated with the de-terministic analogue of the stochastic PPF model.Given the approximate description of the sustained oscillations producedby the simple stochastic PPF model, such as the driven harmonic oscillatorwith noise and the stochastic avian flu model with seasonal forcing, wearrive at the following conclusions. First, the prominent frequencies of thesustained oscillations in the linear stochastic PPF model can be derived fromits deterministic analogue. For example, if we want to identify the dominant784.6. Discussionfrequency of the stochastic perturbations in a stochastic PPF model, wemust determine the response frequency of the parametrically forced linearsystem (i.e. without noise).Second, our description of the stochastically amplified frequencies en-hances our understanding of processes generated by stochastic models withparametric forcing. Take, for example, the PSD of the oscillations observedfrom the seasonally forced stochastic SIR model of Black and McKane [191].Note that in Black and McKane [191], it was shown by using a Van Kampen[9] system-size expansion that a process in the seasonally forced stochasticepidemic model can be expressed as the sum of the deterministic periodicsolution and the stochastic fluctuation scaled by the reciprocal of the square-root of the population size. The evolution of the stochastic fluctuation isderived by linearising the model around the limit cycle solution. This resultsin a linear stochastic PPF model of the form (4.17). Black and McKane [191]observed that the sharp but narrow peak in their PSD corresponds to thelimit cycle frequency while the broader (stochastic) peaks in the PSD aredue to the stochastic amplification of the transients.With the analysis presented here, we can now see that all of the frequen-cies predicted by the aforementioned stochastic peaks in the PSD arise fromthe deterministic aspect of the model represented by (4.18). The fundamen-tal matrix solution of the deterministic system (4.18) has the form (4.22),which contains frequencies ω from R−ωt and, 2pi/T ≡ ν from Q(t). Thesmoothness condition we have for Q(t), as discussed in Section 4.3.2, allowsus to decompose Q(t) into a constant matrix and some constant matricesmultiplied by sinusoidal functions with frequency ν. Hence, using product-to-sum trigonometric identities, the product Q(t)R−ωt have frequencies ωand ν ± ω. From these general principles, we therefore expect that theproduct in (4.53) leads to the emergence of the central frequencies ω andν ± ω observed in the PSD of ξ(t) where the largest power corresponds tofrequency ω. The form of the approximate process in (4.57) also allows usto conclude that the broad peaks observed in the PSD of ξ(t) are due tothe slowly varying stochastic phase perturbation process around the centralfrequencies, which are predicted by the deterministic aspect of the system(4.17).Third, for a linear stochastic PPF model, the noise intensity level con-tributes only to the overall amplification of the system. Thus, for sufficientpopulation size, the stochastic fluctuation may become dominant over thelimit cycle. The closed form approximation may be a good tool to mathe-matically analyze the stochastic PPF disease model to predict the period-icity of outbreaks in the presence of stochasticity and seasonal forcing. For794.6. Discussionstochastic epidemic modelling, the approximation can be used to comparethe importance of parameters (or drivers) that contribute to fluctuationamplification or emergence of frequencies. The question on whether season-ality matters in the presence of demographic stochasticity may be furtherexplored mathematically now that an approximate form of the fluctuationprocess for stochastic PPF models has been derived. The theory of Mathieuequations may enter into any further progress in the analysis of the nonlinearinteractions occurring in the stochastic fluctuation.Finally, it may be that the approximation method is limited by theunknown aspects of the fundamental matrix solution (4.22). This limitationprovides motivation for modellers to develop mathematical tools to computethe fundamental matrix solutions of particular differential equation modelsfor complex systems.80Chapter 5Conclusions and FutureWorkThe main theme of this thesis is the mathematical analysis of recurrentepidemics, in general, and in avian flu in particular. From our overview ofthe existing literature on the epidemiology, ecology, and epidemic recurrenceof avian influenza in Chapter 2, we see that the recurrent behaviour of thedisease is a complex process and further investigation using models andmathematical analysis is needed. SDE models [197] are a promising tool inthis study, but existing analysis has been restricted to computer simulations.The pattern of epidemic recurrence can be described by the period be-tween epidemics, the amplitude or intensity of epidemics, length of inter-epidemic (endemic) time interval, and the amount of variability. This thesisis focused on understanding recurrent pattern of epidemics with respect tothe period and amplitude of epidemics.Epidemic recurrence can be studied by formulating the system as acontinuous-time Markov jump processes, and then by modelling the dy-namics via a master equation. Unfortunately, the only way to obtain exactsolutions to the master equation is by stochastic simulation. Another ap-proach to model the stochastic dynamics is by approximating the Markovjump process by a system of SDEs either via system-size expansion [9], awell-known method in the physical literature, or equivalently via the methodof Kurtz[198]. The stochastic path of the SDEs is then investigated by com-puting the theoretical power spectra which can be used to study how fluc-tuations vary with the parameter of interest. The major drawback with thisapproach is that it cannot fully explain the mechanisms behind observedfluctuations. Analysing how the dominant period and typical amplitudes offluctuating epidemics vary with the parameters is computationally cumber-some.In this thesis, we emphasise the importance of analytic tools to betterunderstand epidemic models exhibiting recurrent behaviour as they shedlight on the underlying mechanisms that were not realised before. Thischapter is devoted to conclusions drawn from previous results. We broadly815.1. Factors influencing epidemic recurrencesummarise them in the following sections. We also discuss potential researchdirections for future considerations.5.1 Factors influencing epidemic recurrenceGenerally speaking, the factors influencing epidemic recurrence of avianinfluenza can be determined from the epidemiology and ecology of the dis-ease. The differences in immune responses of host individuals, host interac-tions and movement, variations in virus strains, and variability of environ-mental factors all contribute to the recurrence pattern. As an example, wewould expect the recurrence period to be shorter in a population with highrate of exposure to the novel avian influenza viruses, or in a population withlow exposure of virus in the environment but higher movement, i.e. highcontact rate [61]. The growing body of knowledge regarding the epidemi-ology and ecology of avian flu offers an opportunity for the formulation ofnew hypotheses to explain epidemic recurrence. One must be cautious informulating the hypothesis as existing literatures pose problem in basic ter-minologies. For example, “outbreaks” and “epidemics” may have the samemeaning since several authors seemed to be using these terms interchange-ably. Apart from word usage, epidemic recurrence is a complex dynamicalphenomenon that mathematical approaches are ideally suited to untangle[11].Our theoretical approach allows us to study, within a stochastic frame-work, factors that are thought to influence epidemic recurrence. As ourstarting point, we considered the stochastic model developed by Wang et al.[1]. We chose this model because it is simple yet complex enough to mimicthe overall dynamics of avian flu recurrence data, e.g. the dominant out-break period. Using an approach different from that in Wang et al. [1],we re-derive the stochastic avian flu model and so deduce that the elementof stochasticity comes from centred Poisson increments. These incrementsrepresent the uncertainty in our knowledge of various aspects of the dis-ease. This inherent uncertainty makes a stochastic model appropriate fordescribing the avian flu system.The presence of demographic stochasticity makes it difficult to identifyexactly how individual factors influence epidemic recurrence. In modellingterms, this task involves determining the parameters in the model that cangive rise to recurrent epidemics and finding meaningful parameter rangesthat are consistent with biological data. In Chapter 3, we overcame thischallenge by extending the method of Baxendale and Greenwood [50] to the825.2. The stochastic fluctuations around a deterministic skeletonthree-dimensional case. We then applied the new result to the stochasticmodel, and found that avian flu epidemic recurrence can be representedas a sum of two processes: (1) a scaled one-dimensional OU process and(2) a scaled rotation multiplied by a slowly-varying bi-variate standard OUprocess. We conclude from this formulation that the stochastic elementsof avian flu recurrence act as perturbations of the phase and amplitude ofthe transient deterministic dynamics. The stochastic perturbations of OUprocesses follow a distribution that has been studied in the literature [46].By identifying the stochasticity of avian flu as arising from an OU process,we are able to understand the role of stochasticity in modulating the patternof epidemic recurrence for avian flu.One factor that is identified to have a key role in driving avian fluepidemic recurrence, by previous modellers, is environmental transmission.However, our study shows that having a sufficient amount of transmission,from either direct or environmental sources, is an important ingredient forrecurrence. Our analytic approach has shed light on this fact and found thatwhen one transmission rate is weak, the other has a critical role in alteringthe features of epidemic recurrence, e.g. the typical amplitude and dominantperiod. We conclude that in a stochastic setting, either direct or indirecttransmission is an important factor that influences epidemic recurrence andunderstanding their contributions may be necessary for the development ofappropriate policies for controlling the disease.5.2 The stochastic fluctuations around adeterministic skeletonWe defined stochastic fluctuations as the ”fuzziness” around the deter-ministic dynamics. This description suggests that stochastic noise in thisstudy has a passive role.In this thesis, we considered two models of avian flu epidemic recurrence:(1) the stochastic host-pathogen model (the SIR-V model) developed byWang et al. [1] and (2) an extension of the SIR-V model incorporating aseasonally forced environmental transmission rate. The main motivation forconsidering the latter model is the fact that environmental transmission hasan important role in explaining multi-year periodicity in avian flu epidemics[1, 13, 165] and is well-known to show a seasonal pattern. We note here thatthough the second model is more complex than the first one, we can stillmathematically contrast the stochastic fluctuations that they generate.Our novel approach which is discussed in Chapter 4, shows that the pro-835.3. Avian flu in a stochastic and seasonally forced environmentcess describing stochastic fluctuations around a periodic solution induced by(seasonal) forcing can be approximated by another periodic process mod-ulated by a standard Ornstein-Uhlenbeck process. We conclude that, forboth types of stochastic fluctuations (with or without seasonal forcing),stochasticity perturbs the deterministic dynamics in a similar way. Stochas-tic fluctuations around an endemic equilibrium are noise-sustained oscilla-tions whose dominant frequency is close to the intrinsic frequency of thedeterministic version. On the other hand, fluctuations around a forced limitcycle are also noise-sustained oscillations but their dominant frequency isclose to the intrinsic frequency of the deterministic perturbation of the limitcycle, not necessarily the intrinsic frequency of the deterministic version ofthe model. This point implies that the key to characterising these stochasticfluctuations is the understanding of the deterministic perturbation aroundthe limit cycle, which gives rise to the product Q(t)R−ωt in (4.49).Rozhnova and Nunes [44] showed how the population size affects the fullstochastic dynamics of a forced stochastic epidemic system. Since we havean approximate process that describes the stochastic fluctuation aroundthe limit cycle, we can use the approximate process to investigate how theresonance produced by the interaction of a deterministic periodic solutionand the stochastic fluctuation vary with population size. This way, we cananswer the question: To what extent does population size affect the fullstochastic dynamics? The answer to this question relies on the magnitude ofthe stochastic fluctuations which is approximately determined by the scalarfactor σ¯/√λ in (4.49).5.3 Avian flu in a stochastic and seasonallyforced environmentA typical approach for analysing avian flu epidemic recurrence is byincorporating demographic stochasticity. In this formulation, we do nothave to make any assumptions about the form or strength of the noise thesystem is subjected because the variance is defined in the formulation of thestochastic model. In this thesis, we have explored the possibility of modellingavian flu epidemic recurrence by combining stochastic and seasonal effects.The way we include seasonal effects is by representing the environmental(virus-host) transmission rate as a time-periodic function. In this setting,the dynamics of avian flu epidemics exhibit a limit cycle (i.e. forced limitcycle) as the intensity of forcing increases. The frequency of the epidemics,however, does not change with the forcing intensity but depends on the845.3. Avian flu in a stochastic and seasonally forced environmentperiodicity of the environmental transmission rate. Moreover, it has beenshown previously that for a low direct transmission rate, the period of avianflu epidemics can be influenced by the environmental transmission rate thatmeasures how frequently the virus is transmitted to the host from the envi-ronment (See Chapter 3).Our theoretical approach describes avian flu epidemic recurrence, pro-duced from a seasonally forced and stochastic environment, as oscillationsthat can be expressed as a sum of the forced limit cycle and stochastic fluctu-ations scaled by the square-root of the total population size. The stochasticfluctuations have a dominant frequency that largely depends on the intrinsicfrequency of the deterministic perturbation of the limit cycle although otherfrequencies, e.g. the forcing frequency and the limit cycle frequency, are alsoinvolved in driving the fluctuations. Now in the case where the dominantfrequency of the stochastic fluctuations is larger than the frequency of thedeterministic dynamics, we mathematically identified possible mechanismswhereby stochastic effects could dominate over seasonal effects, based on theform of the approximate process we formulated in Chapter 4:1. When the population size is sufficiently small. For a density-dependentprocess, it is know that its stochastic dynamics can be regarded asthe sum of the macroscopic dynamics and the stochastic fluctuationsdivided by the square-root of the population size N . Therefore, if Nis sufficiently small, the stochastic fluctuations term would be greaterthan the macroscopic term hence, will have substantial effects to thefull dynamics.2. The time-average of noise variance is sufficiently large. This pointis due to the fact that the amplitude of the approximate process forstochastic fluctuations tends to increase as the time-average of thenoise variance increases.3. The decay rate of the deterministic perturbation around the periodicsolution is small enough. This point is related to the form of theapproximate process, which is inversely proportional to the square-root of the decay rate. The decay rate is determined by the real partof the Floquet exponent associated to the deterministic perturbationof the limit cycle.For future work, it would be interesting to identify the components thatmake up the time-average noise variance and the decay rate of the deter-ministic perturbation. Furthermore, a detailed numerical or mathematical855.4. Open questionsstudy on the conclusions drawn above would be useful for validation andexplore the active role of stochastic noise in the full dynamics. Our analysiscould pave a way to develop new tools in quantifying the differences betweeneffects of stochasticity and parametric forcing for other SDE models fromdifferent disciplines.5.4 Open questionsThe results and analysis presented in this thesis lead us to ask a varietyof open questions left for the reader to explore. We end this thesis bypresenting the following open questions:1. Are there other biological factors that contribute to periodicity andamplitude of recurrence pattern in the avian flu system?2. It is known that, in a stochastic setting, either direct or environmen-tal transmission has a major influence in the dynamics of avian fluepidemics. How does this result change in the presence of other mech-anism such as seasonality?3. Given that the nature of stochastic fluctuations describing epidemicrecurrence can be determined via an approximate solution, how canthese results be used to understand other diseases displaying epidemicrecurrence? How can we quantify the distinction between two diseases,e.g. measles and whooping cough, displaying epidemic recurrence ac-cording to the nature of their stochastic oscillations?4. How will the results and analysis change if environmental stochasticityis involved?5. How can we use the analytic methods here to approximately determineother features of recurrence pattern, e.g. the length of inter-epidemictime interval, amount of variability?6. How does our mathematical results be used for fitting avian flu datadisplaying recurrent pattern?This thesis has provided new tools for studying stochastic disease dy-namics, or any dynamics with noise-sustained oscillations or noisy limit cy-cles. One can use these tools as starting point to answer the open questionsabove. For instance, the first three questions above can be answered by in-vestigating the effects of other parameters and applying the methods to other865.4. Open questionsdiseases, i.e. as done in Chapter 3 of this thesis. On the other hand, thegeneral derivation of the stochastic system with parametric forcing shownin Chapter 4 may be a useful framework for model involving environmentalstochasticity as mentioned in question 4. The approach for answering thelast two questions are unclear at this point but we think that a careful anal-ysis of the phase and amplitude processes of the OU process, as well as itsassociated curve fitting techniques, is necessary to address them.87Bibliography[1] R H Wang, Z Jin, Q X Liu, J van de Koppel, and D Alonso. Asimple stochastic model with environmental transmission explainsmulti-year periodicity in outbreaks of avian flu. PloS one, 7(2):e28873,January 2012. ISSN 1932-6203. doi: 10.1371/journal.pone.0028873.URL http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3281819&tool=pmcentrez&rendertype=abstracthttp://dx.plos.org/10.1371/journal.pone.0028873.g003. → pagesviii, ix, xi, 2, 3, 5, 6, 18, 22, 23, 25, 26, 27, 28, 29, 30, 35, 43, 44, 46,47, 70, 82, 83, 113, 114, 121, 122, 123[2] Lars Folke Olsen, Gregory L Truty, and William Morris Schaffer. Os-cillations and chaos in epidemics: a nonlinear dynamic study of sixchildhood diseases in copenhagen, denmark. Theoretical populationbiology, 33(3):344–370, 1988. → pages ix, 2, 3[3] S Altizer, A Dobson, P Hosseini, P Hudson, M Pascual, and P Ro-hani. Seasonality and the dynamics of infectious diseases, 2006. ISSN1461023X. → pages ix, 3, 5, 19, 49[4] Mercedes Pascual, Menno J Bouma, and Andrew P Dobson. Choleraand climate: revisiting the quantitative evidence. Microbes and Infec-tion, 4(2):237–245, 2002. → pages ix, 3[5] T Coulson, P Rohani, and M Pascual. Skeletons, noise and populationgrowth: The end of an old debate? Trends in Ecology and Evolution,19(7):359–364, 2004. ISSN 01695347. doi: 10.1016/j.tree.2004.05.008.→ pages ix, 5, 6, 8[6] A J Black and A J McKane. Stochasticity in staged models ofepidemics: quantifying the dynamics of whooping cough. Jour-nal of the Royal Society, Interface / the Royal Society, 7(49):1219–27, August 2010. ISSN 1742-5662. doi: 10.1098/rsif.2009.0514. URL http://www.pubmedcentral.nih.gov/articlerender.88Bibliographyfcgi?artid=2894870&tool=pmcentrez&rendertype=abstract. →pages 1[7] Matthew O Jackson and Leeat Yariv. Diffusion on social networks.Economie publique/Public economics, (16), 2006. → pages 1[8] Noam Berger, Christian Borgs, Jennifer T Chayes, and Amin Saberi.On the spread of viruses on the internet. In Proceedings of the sixteenthannual ACM-SIAM symposium on Discrete algorithms, pages 301–310. Society for Industrial and Applied Mathematics, 2005. → pages1[9] N G Van Kampen. Stochastic processes in physics and chemistry,volume 11. North- Holland personal library, Amsterdam, 2nd edition,1992. ISBN 0444893490. doi: 10.2307/2984076. URL http://books.google.com/books?hl=en&lr=&id=3e7XbMoJzmoC&pgis=1. → pages1, 9, 50, 55, 70, 79, 81, 113[10] T G Kurtz. Strong approximation theorems for density dependentMarkov chains, 1978. ISSN 03044149. → pages 1, 9, 23, 46, 50, 52,53, 70, 113[11] A J Black. The stochastic dynamics of epidemic models. PhD thesis,2010. → pages 2, 5, 49, 50, 55, 77, 82[12] Wayne P London and James A Yorke. Recurrent outbreaks of measles,chickenpox and mumps i. seasonal variation in contact rates. Americanjournal of epidemiology, 98(6):453–468, 1973. → pages 2, 5[13] R Breban and J M Drake. The role of environmen-tal transmission in recurrent avian influenza epidemics.PLoS Computational Biology, 5(4):e1000346, April 2009.ISSN 1553-7358. doi: 10.1371/journal.pcbi.1000346. URLhttp://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2660440&tool=pmcentrez&rendertype=abstracthttp://dx.plos.org/10.1371/journal.pcbi.1000346.g007. → pages 2,21, 47, 83[14] Chris T Bauch and David JD Earn. Transients and attractors inepidemics. Proceedings of the Royal Society of London B: BiologicalSciences, 270(1524):1573–1578, 2003. → pages 289Bibliography[15] Petre Stoica and Randolph L Moses. Spectral analysis of signals, vol-ume 452. Pearson Prentice Hall Upper Saddle River, NJ, 2005. →pages 2, 10[16] May Anne Mata, Meghan Dutot, Leah Edelstein-Keshet, andWilliam R Holmes. A model for intracellular actin waves exploredby nonlinear local perturbation analysis. Journal of theoretical biol-ogy, 334:149–161, 2013. → pages 4[17] Susanne Ditlevsen and Priscilla Greenwood. The morris–lecar neuronmodel embeds a leaky integrate-and-fire model. Journal of mathemat-ical biology, 67(2):239–259, 2013. → pages[18] James Dickson Murray. Mathematical biology [electronic resource].:An introduction. Springer, 2002. → pages[19] James D Murray. Mathematical Biology. II Spatial Models andBiomedical Applications {Interdisciplinary Applied Mathematics V.18}. Springer-Verlag New York Incorporated, 2001. → pages 4[20] William O Kermack and Anderson G McKendrick. Contributions tothe mathematical theory of epidemicsi. Bulletin of mathematical biol-ogy, 53(1):33–55, 1991. → pages 4[21] Roy M Anderson, Robert M May, and B Anderson. Infectious diseasesof humans: dynamics and control, volume 28. Wiley Online Library,1992. → pages 5[22] Dieter Schenzle. An age-structured model of pre-and post-vaccinationmeasles transmission. Mathematical Medicine and Biology, 1(2):169–191, 1984. → pages 5, 6[23] Herbert W Hethcote. An age-structured model for pertussis transmis-sion. Mathematical biosciences, 145(2):89–136, 1997. → pages[24] Herbert W Hethcote and P Van Den Driessche. Two sis epidemiologicmodels with delays. Journal of mathematical biology, 40(1):3–26, 2000.→ pages[25] Herbert W Hethcote and P Van den Driessche. Some epidemiologicalmodels with nonlinear incidence. Journal of Mathematical Biology, 29(3):271–287, 1991. → pages 590Bibliography[26] R M Nisbet and W Gurney. Modelling fluctuating populations. JohnWiley and Sons Limited, 1982. → pages 5, 22, 50, 77[27] Nicholas C Grassly and Christophe Fraser. Seasonal infectious diseaseepidemiology. Proceedings of the Royal Society of London B: BiologicalSciences, 273(1600):2541–2550, 2006. → pages 5[28] MJ Keeling, Pejman Rohani, and BT Grenfell. Season-ally forced disease dynamics explored as switching between at-tractors. Physica D: Nonlinear Phenomena, 148(3-4):317–335,2001. URL http://www.sciencedirect.com/science/article/pii/S0167278900001871. → pages 5, 6, 49[29] Pejman Rohani, Matthew J Keeling, and Bryan T Grenfell. The in-terplay between determinism and stochasticity in childhood diseases.The American Naturalist, 159(5):469–481, 2002. → pages 5, 6, 7[30] MS Bartlett. Deterministic and stochastic models for recurrent epi-demics. In Proceedings of the third Berkeley symposium on mathe-matical statistics and probability, volume 4, page 109, 1956. → pages5[31] J P Aparicio and H G Solari. Sustained oscillations instochastic systems. Mathematical biosciences, 169(1):15–25,2001. URL http://www.sciencedirect.com/science/article/pii/S002555640000050X. → pages 5, 22[32] A J McKane and T J Newman. Predator-prey cycles from reso-nant amplification of demographic stochasticity. Physical review let-ters, 94(21):4, January 2005. doi: 10.1103/PhysRevLett.94.218102.URL http://arxiv.org/abs/q-bio/0501023http://prl.aps.org/abstract/PRL/v94/i21/e218102. → pages 5, 22, 49, 53[33] D Alonso, A J McKane, and M Pascual. Stochastic amplifica-tion in epidemics. Journal of the Royal Society Interface, 4(14):575–582, June 2007. ISSN 1742-5689. doi: 10.1098/rsif.2006.0192.URL http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2373404&tool=pmcentrez&rendertype=abstracthttp:// → pages 6, 49, 50,53, 54, 55, 70, 7791Bibliography[34] A J McKane, J D Nagy, T J Newman, and M O Stefanini. Ampli-fied biochemical oscillations in cellular systems. Journal of Statisti-cal Physics, 128(1-2):165–191, 2007. ISSN 00224715. doi: 10.1007/s10955-006-9221-9. → pages 5, 22[35] R Kuske, L F Gordillo, and P Greenwood. Sustained oscilla-tions via coherence resonance in SIR. Journal of TheoreticalBiology, 245(3):459–469, April 2007. ISSN 0022-5193. doi:10.1016/j.jtbi.2006.10.029. URL http://www.ncbi.nlm.nih.gov/pubmed/17173935http://www.sciencedirect.com/science/article/pii/S0022519306004930. → pages 6, 22, 50[36] Jon Greenman, Masashi Kamo, and Mike Boots. External forcing ofecological and epidemiological systems: a resonance approach. PhysicaD: Nonlinear Phenomena, 190(1):136–151, 2004. → pages 6[37] A J Black and A J McKane. Stochastic formulation of ecological mod-els and their applications. Trends in ecology & evolution, 27(6):337–345, June 2012. ISSN 0169-5347. doi: 10.1016/j.tree.2012.01.014. URLhttp://www.ncbi.nlm.nih.gov/pubmed/22406194http://www.sciencedirect.com/science/article/pii/S016953471200033X.→ pages 6, 7, 50, 70[38] Ira B Schwartz. Multiple stable recurrent outbreaks and predictabilityin seasonally forced nonlinear epidemic models. Journal of mathemat-ical biology, 21(3):347–361, 1985. → pages 6[39] J L Aron and I B Schwartz. Seasonality and period-doubling bifurca-tions in an epidemic model. Journal of theoretical biology, 110:665–679,1984. ISSN 00225193. doi: 10.1016/S0022-5193(84)80150-2. → pages7, 49[40] Yu A Kuznetsov and C Piccardi. Bifurcation analysis of periodic seirand sir epidemic models. Journal of mathematical biology, 32(2):109–121, 1994. → pages[41] David JD Earn, Pejman Rohani, Benjamin M Bolker, and Bryan TGrenfell. A simple model for complex dynamical transitions in epi-demics. Science, 287(5453):667–670, 2000. → pages 7[42] N Tuncer and M Martcheva. Modeling seasonality in avian in-fluenza h5n1. pages 1–27, 2013. ISSN 02183390. doi: 10.1142/S0218339013400044. → pages 7, 19, 2292Bibliography[43] Food and Agriculture Organization of the United Nations. Animalproduction and health division. http://www.fao.org/avianflu/en/overview.htm, 2013. → pages 7[44] G. Rozhnova and a. Nunes. Stochastic effects in a seasonally forcedepidemic model. Physical Review E, 82(4):041906, oct 2010. ISSN1539-3755. doi: 10.1103/PhysRevE.82.041906. URL http://link.aps.org/doi/10.1103/PhysRevE.82.041906. → pages 7, 49, 50, 55,84[45] Steven H Strogatz. Nonlinear dynamics and chaos: with applicationsto physics, biology, chemistry, and engineering. Westview press, 2014.→ pages 9[46] C Gardiner. Handbook of stochastic methods for physics, chemistryand the natural sciences, 1986. → pages 9, 10, 34, 47, 62, 83[47] P E Greenwood and L F Gordillo. Stochastic epidemic modeling.Mathematical and Statistical Estimation Approaches in Epidemiology,pages 1–23, 2009. URL http://link.springer.com/chapter/10.1007/978-90-481-2313-1_2. → pages 9, 46, 111[48] E J Allen, L J Allen, A Arciniega, and P E Greenwood. Constructionof Equivalent Stochastic Differential Equation Models. Stoch. Anal.Appl., 26:274–297, 2008. → pages 71, 113, 116[49] Linda JS Allen. An introduction to stochastic processes with applica-tions to biology. CRC Press, 2010. → pages 9[50] P H Baxendale and P E Greenwood. Sustained oscil-lations for density dependent Markov processes. Jour-nal of Mathematical Biology, 63(3):433–457, September2011. ISSN 1432-1416. doi: 10.1007/s00285-010-0376-2.URL http://www.ncbi.nlm.nih.gov/pubmed/21076832http://link.springer.com/article/10.1007/s00285-010-0376-2. →pages 9, 10, 22, 33, 50, 51, 59, 60, 61, 78, 82, 113, 115, 116[51] G.N. Milstein. Numerical integration of stochastic differential equa-tions. 1995. ISBN 9789048144877. doi: 10.1007/BF01015322. →pages 10, 50[52] Roger Grimshaw. Nonlinear ordinary differential equations, volume 2.CRC Press, 1991. → pages 10, 51, 55, 57, 12593Bibliography[53] Thomas G Kurtz. Lectures on stochastic analysis. Department ofMathematics and Statistics, University of Wisconsin, Madison, WI,pages 53706–1388, 2001. → pages 10, 52[54] MJ Stear. Oie manual of diagnostic tests and vaccines for terrestrialanimals (mammals, birds and bees) 5th edn. volumes 1 and 2. worldorganization for animal health 2004. isbn 92 9044 622 6. 140., 2005.→ pages 12[55] R G Webster and W J Bean. Evolution and ecology of influenza Aviruses. Microbiological Reviews, 56(1):152–179, 1992. URL http://mmbr.asm.org/content/56/1/152.short. → pages 12, 13, 15, 21[56] Ilaria Capua and Dennis J Alexander. Animal and human healthimplications of avian influenza infections. Bioscience Reports, 27(6):359–372, 2007. → pages 12[57] Ilaria Capua and Stefano Marangon. Control and prevention of avianinfluenza in an evolving scenario. Vaccine, 25(30):5645–5652, 2007. →pages 12[58] Centre for Disease Control and Prevention avian influenza a (h7n9)virus. http://www.cdc.gov/flu/avianflu/h7n9-virus.htm, 2016.Accessed: 2016-07-07. → pages 12[59] World Health Organization avian influenza in humans.http://www.who.int/influenza/human_animal_interface/avian_influenza/en/, 2016. Accessed: 2016-07-07. → pages[60] M E E Zowalaty and S A Bustin. Avian influenza: virology, diagnosand surveillance. Future . . . , pages 1209–1227, 2013. URL http://www.futuremedicine.com/doi/abs/10.2217/fmb.13.81. → pages12[61] M Martcheva. Avian Flu: Modeling and Implications for Control.Journal of Biological Systems, 22(01):151–175, March 2014. ISSN0218-3390. doi: 10.1142/S0218339014500090. URL http://www.worldscientific.com/doi/abs/10.1142/S0218339014500090. →pages 12, 82[62] Nicole M Bouvier and Peter Palese. The biology of influenza viruses.Vaccine, 26:D49–D53, 2008. → pages 1294Bibliography[63] David E Swayne. Avian influenza. John Wiley & Sons, 2009. → pages12[64] Chang-Won Lee and Yehia M Saif. Avian influenza virus. Compara-tive immunology, microbiology and infectious diseases, 32(4):301–310,2009. → pages 12[65] Ron AM Fouchier, Vincent Munster, Anders Wallensten, Theo MBestebroer, Sander Herfst, Derek Smith, Guus F Rimmelzwaan, Bjo¨rnOlsen, and Albert DME Osterhaus. Characterization of a novel in-fluenza a virus hemagglutinin subtype (h16) obtained from black-headed gulls. Journal of virology, 79(5):2814–2822, 2005. → pages12, 13[66] D J Alexander. A review of avian influenza in different birdspecies. Veterinary microbiology, 74, 2000. URL http://www.sciencedirect.com/science/article/pii/S0378113500001607.→ pages 13, 15, 21[67] DE Swayne and DL Suarez. Highly pathogenic avian influenza. Re-vue Scientifique et Technique-office International des Epizooties, 19(2):463–475, 2000. → pages 13[68] Christopher S Jennelle, Michelle Carstensen, Erik C Hildebrand, LouisCornicelli, P Wolf, DA Grear, Hon S Ip, Kaci K Vandalen, andLarissa A Minicucci. Surveillance for highly pathogenic avian influenzavirus in wild birds during outbreaks in domestic poultry, minnesota,2015. Emerging infectious diseases, 22(7), 2016. → pages 13[69] Bethany J Hoye, Vincent J Munster, Hiroshi Nishiura, MarcelKlaassen, and Ron AM Fouchier. Surveillance of wild birds for avianinfluenza virus. Emerging infectious diseases, 16(12):1827–1834, 2010.→ pages 13[70] Kennedy F Shortridge, Nan Nan Zhou, Yi Guan, Peng Gao, ToshihiroIto, Yoshihiro Kawaoka, Shantha Kodihalli, Scott Krauss, DeborahMarkwell, K Gopal Murti, et al. Characterization of avian h5n1 in-fluenza viruses from poultry in hong kong. Virology, 252(2):331–342,1998. → pages 13[71] Erica Spackman, Dennis A Senne, Sherrill Davison, and David LSuarez. Sequence analysis of recent h7 avian influenza viruses as-sociated with three different outbreaks in commercial poultry in the95Bibliographyunited states. Journal of virology, 77(24):13399–13402, 2003. → pages15[72] Ruben O Donis. Evolution of h5n1 avian influenza viruses in asia.Emerging infectious diseases, 11(10), 2005. → pages 13[73] M Perdue, J Crawford, M Garcia, J Latimer, and D Swayne. Occur-rence and possible mechanisms of cleavage-site insertions in the avianinfluenza hemagglutinin gene. Avian Diseases, pages 182–193, 2003.→ pages 13[74] JL Heeney. Zoonotic viral diseases and the frontier of early diagnosis,control and prevention. Journal of internal medicine, 260(5):399–408,2006. → pages 13[75] Susan J Baigent and John W McCauley. Influenza type a in humans,mammals and birds: determinants of virus virulence, host-range andinterspecies transmission. Bioessays, 25(7):657–671, 2003. → pages[76] Rudolf Deibel, Diane E Emord, Ward Dukelow, Virginia S Hinshaw,and John M Wood. Influenza viruses and paramyxoviruses in ducksin the atlantic flyway, 1977-1983, including an h5n2 isolate related tothe virulent chicken virus. Avian Diseases, pages 970–985, 1985. →pages 13, 19[77] Linda Widjaja, Scott L Krauss, Richard J Webby, Tao Xie, andRobert G Webster. Matrix gene of influenza a viruses isolated fromwild aquatic birds: ecology and emergence of influenza a viruses. Jour-nal of virology, 78(16):8771–8779, 2004. → pages 13[78] World Health Organization avian influenza a (h7n9) virus.http://www.who.int/influenza/human_animal_interface/influenza_h7n9/en/, 2016. Accessed: 2016-07-07. → pages13[79] VS Hinshaw, RG Webster, and RJ Rodriguez. Influenza a viruses:combinations of hemagglutinin and neuraminidase subtypes isolatedfrom animals and other sources. Archives of virology, 67(3):191–201,1981. → pages 13[80] Dennis J Alexander. Ecology of avian influenza in domestic birds.Emergence and control of zoonotic ortho-and paramyxovirus diseases,pages 25–33, 2001. → pages 1396Bibliography[81] WB Becker. The isolation and classification of tern virus: influenzavirus a/tern/south africa/1961. Journal of Hygiene, 64(03):309–320,1966. → pages 13[82] DE Stallknecht and SM Shane. Host range of avian influenza virusin free-living birds. Veterinary research communications, 12(2-3):125–141, 1988. → pages 13[83] GB Sharp, Y Kawaoka, SM Wright, B Turner, V Hinshaw, andRG Webster. Wild ducks are the reservoir for only a limited numberof influenza a subtypes. Epidemiology and Infection, 110(01):161–176,1993. → pages 13, 21[84] Yoshihiro Kawaoka, Thomas M Chambers, William L Sladen, andRobert Gwebster. Is the gene pool of influenza viruses in shorebirdsand gulls different from that in wild ducks? Virology, 163(1):247–250,1988. → pages 13[85] B Olsen, V J Munster, and A Wallensten. Global Patterns of Influenzaa Virus in Wild Birds. science, 312(5772):384–388, 2006. URL http://www.sciencemag.org/content/312/5772/384.short. → pages 13,21[86] Jeong-Ki Kim, Nicholas J Negovetich, Heather L Forrest, andRobert G Webster. Ducks: the trojan horses of h5n1 influenza. In-fluenza and other respiratory viruses, 3(4):121–128, 2009. → pages13[87] DE Stallknecht, MT Kearney, SM Shane, and PJ Zwank. Effects ofph, temperature, and salinity on persistence of avian influenza virusesin water. Avian diseases, pages 412–418, 1990. → pages 14, 15[88] David E Stallknecht, Justin D Brown, and DE Swayne. Ecology ofavian influenza in wild birds. Avian influenza, 1:43–58, 2008. → pages14, 15[89] William B Karesh, Robert A Cook, Elizabeth L Bennett, and JamesNewcomb. Wildlife trade and global disease emergence. Emerg InfectDis, 11(7):1000–1002, 2005. → pages 14[90] Kurt D Reed, Jennifer K Meece, James S Henkel, and Sanjay KShukla. Birds, migration and emerging zoonoses: West nile virus,lyme disease, influenza a and enteropathogens. Clinical medicine &research, 1(1):5–12, 2003. → pages 1497Bibliography[91] HA Westbury, AJ Turner, and L Kovesdy. The pathogenicity of threeaustralian fowl plague viruses for chickens, turkeys and ducks. Veteri-nary Microbiology, 4(3):223–234, 1979. → pages 14[92] VS Hinshaw, RG Webster, and B Turner. The perpetuation of or-thomyxoviruses and paramyxoviruses in canadian waterfowl. Cana-dian Journal of Microbiology, 26(5):622–629, 1980. → pages[93] DE Stallknecht, SM Shane, MT Kearney, and PJ Zwank. Persistenceof avian influenza viruses in water. Avian diseases, pages 406–411,1990. → pages 14, 15, 70[94] KM Sturm-Ramirez, DJ Hulse-Post, EA Govorkova, J Humberd,P Seiler, P Puthavathana, C Buranathai, TD Nguyen, A Chaisingh,HT Long, et al. Are ducks contributing to the endemicity of highlypathogenic h5n1 influenza virus in asia? Journal of virology, 79(17):11269–11279, 2005. → pages 14[95] Robert G Webster, M Yakhno, V S Hinshaw, W J Bean, and K CMurti. Intestinal influenza: replication and characterization of in-fluenza viruses in ducks. Virology, 84(2):268–278, 1978. → pages 14[96] DJ Alexander, WH Allan, DG Parsons, and G Parsons. Thepathogenicity of four avian influenza viruses for fowls, turkeys andducks. Research in veterinary science, 24(2):242–247, 1978. → pages14[97] D J Alexander, G Parsons, and R J Manvell. Experimental assessmentof the pathogenicity of eight avian influenza A viruses of H5 subtypefor chickens, turkeys, ducks and quail. Avian pathology : journal ofthe W.V.P.A, 15(4):647–662, 1986. ISSN 0307-9457. doi: 10.1080/03079458608436328. → pages 48[98] O Narayan, G Lang, and BT Rouse. A new influenza a virus infectionin turkeys. Archiv fu¨r die gesamte Virusforschung, 26(1-2):149–165,1969. → pages[99] HA Westbury, AJ Turner, and L Amon. Transmissibility of two avianinfluenza a viruses (h7 n6) between chickens. Avian Pathology, 10(4):481–487, 1981. → pages 14[100] Patrik Ellstro¨m, Neus Latorre-Margalef, Petra Griekspoor, JonasWaldenstro¨m, Jenny Olofsson, John Wahlgren, and Bjo¨rn Olsen. Sam-pling for low-pathogenic avian influenza a virus in wild mallard ducks:98Bibliographyoropharyngeal versus cloacal swabbing. Vaccine, 26(35):4414–4416,2008. → pages 14[101] Justin D Brown, David E Stallknecht, Joan R Beck, David L Suarez,and David E Swayne. Susceptibility of north american ducks and gullsto h5n1 highly pathogenic avian influenza viruses. Emerg Infect Dis,12(11):1663–1670, 2006. → pages 15[102] David E Stallknecht and Justin D Brown. Wild birds and the epi-demiology of avian influenza. Journal of wildlife diseases, 43(3), 2007.→ pages[103] Virginia S Hinshaw, RG Webster, and B Turner. Water-borne trans-mission of influenza a viruses? Intervirology, 11(1):66–68, 1979. →pages 15[104] Scott O Rogers, William T Starmer, and John D Castello. Recyclingof pathogenic microbes through survival in ice. Medical hypotheses, 63(5):773–777, 2004. → pages 15[105] Dany Shoham. Biotic-abiotic mechanisms for long-term preservationand reemergence of influenza type a virus genes. Progress in medicalvirology, 40:178–192, 1993. → pages[106] Alvin W Smith, Douglas E Skilling, John D Castello, and Scott ORogers. Ice as a reservoir for pathogenic human viruses: specifically,caliciviruses, influenza viruses, and enteroviruses. Medical Hypotheses,63(4):560–566, 2004. → pages[107] Gang Zhang, Dany Shoham, David Gilichinsky, Sergei Davydov,John D Castello, and Scott O Rogers. Evidence of influenza a virusrna in siberian lake ice. Journal of virology, 80(24):12229–12235, 2006.→ pages 15[108] T Ito, K Okazaki, Y Kawaoka, A Takada, RG Webster, and H Kida.Perpetuation of influenza a viruses in alaskan waterfowl reservoirs.Archives of virology, 140(7):1163–1172, 1995. → pages 15[109] Terrence M Tumpey, David L Suarez, Laura EL Perkins, Dennis ASenne, Jae-gil Lee, Youn-Jeong Lee, In-Pil Mo, Haan-Woo Sung, andDavid E Swayne. Characterization of a highly pathogenic h5n1 avianinfluenza a virus isolated from duck meat. Journal of virology, 76(12):6344–6355, 2002. → pages 1599Bibliography[110] N Kishida, Y Sakoda, N Isoda, K Matsuda, M Eto, Y Sunaga,T Umemura, and H Kida. Pathogenicity of h5 influenza viruses forducks. Archives of virology, 150(7):1383–1392, 2005. → pages 15[111] Ilaria Capua and Dennis J Alexander. Avian influenza: recent devel-opments. Avian Pathology, 33(4):393–404, 2004. → pages 15[112] Dennis J Alexander. An overview of the epidemiology of avian in-fluenza. Vaccine, 25(30):5637–5644, 2007. → pages 16, 17[113] Blanca Lupiani and Sanjay M Reddy. The history of avian influenza.Comparative Immunology, Microbiology and Infectious Diseases, 32(4):311–323, 2009. → pages 15, 17[114] RAM Fouchier, VJ Munster, et al. Epidemiology of low pathogenicavian influenza viruses in wild birds. Revue scientifique et technique,28(1):49, 2009. → pages 16[115] Ian H Brown. Summary of avian influenza activity in europe, asia,and africa, 2006-2009. Avian diseases, 54(s1):187–193, 2010. → pages16[116] Gerwin Claes, Sylvie Marche´, Jeroen Dewulf, Thierry van den Berg,and Be´ne´dicte Lambrecht. An experimental model to analyse the riskof introduction of a duck-originated h5 low-pathogenic avian influenzavirus in poultry through close contact and contaminative transmission.Epidemiology and infection, 142(09):1836–1847, 2014. → pages 16, 17,18[117] World animal health information database detailed country(ies) dis-ease incidence. http://www.oie.int/wahis_2/public/wahid.php/Diseaseinformation/statusdetail, 2016. Accessed: 2016-07-07. →pages 16, 17, 18[118] Dennis J Alexander. Summary of avian influenza activity in europe,asia, africa, and australasia, 2002-2006. Avian diseases, 51(s1):161–166, 2007. → pages 16[119] I Capua, S Marangon, M Dalla Pozza, C Terregino, and G Cattoli.Avian influenza in italy 1997-2001. Avian diseases, 47(s3):839–843,2003. → pages 16100Bibliography[120] Susan C Trock and John P Huntley. Surveillance and control of avianinfluenza in the new york live bird markets. Avian diseases, 54(s1):340–344, 2010. → pages 16[121] C Probst, JM Gethmann, HJ Petermann, J Neudecker, K Jacobsen,and FJ Conraths. Low pathogenic avian influenza h7n7 in domesticpoultry in germany in 2011. Veterinary Record, pages vetrec–2012,2012. → pages 16[122] Laura Campitelli, Concetta Fabiani, Simona Puzelli, AlessandroFioretti, Emanuela Foni, Alessandra De Marco, Scott Krauss,Robert G Webster, and Isabella Donatelli. H3n2 influenza virusesfrom domestic chickens in italy: an increasing role for chickens in theecology of influenza? Journal of General Virology, 83(2):413–420,2002. → pages 17[123] Yi Peng, Zhi-xun Xie, Jia-bo Liu, Yao-shan Pang, Xian-wen Deng,Zhi-qin Xie, Li-ji Xie, Qing Fan, and Si-si Luo. Epidemiological surveil-lance of low pathogenic avian influenza virus (lpaiv) from poultry inguangxi province, southern china. PLoS One, 8(10):e77132, 2013. →pages 17[124] PS Chin, E Hoffmann, R Webby, RG Webster, Y Guan, M Peiris, andKF Shortridge. Molecular evolution of h6 influenza viruses from poul-try in southeastern china: prevalence of h6n1 influenza viruses possess-ing seven a/hong kong/156/97 (h5n1)-like genes in poultry. Journalof virology, 76(2):507–516, 2002. → pages 17[125] PR Woolcock, DL Suarez, and D Kuney. Low-pathogenicity avian in-fluenza virus (h6n2) in chickens in california, 2000-02. Avian diseases,47(s3):872–881, 2003. → pages 17[126] Saad Gharaibeh. Pathogenicity of an avian influenza virus serotypeh9n2 in chickens. Avian diseases, 52(1):106–110, 2008. → pages 17[127] DJ Alexander. Report on avian influenza in the eastern hemisphereduring 1997–2002. 2009. → pages 17[128] Hassan Nili and Keramat Asasi. Natural cases and an experimentalstudy of h9n2 avian influenza in commercial broiler chickens of iran.Avian Pathology, 31(3):247–252, 2002. → pages 17[129] Dennis A Senne. Avian influenza in north and south america, 2002-2005. Avian diseases, 51(s1):167–173, 2007. → pages 17101Bibliography[130] George Fleming. Animal plagues: their history, nature, and preven-tion, volume 1. London, Chapman and Hall, 1871. → pages 17[131] Edoardo Perroncito. Epizoozia tifoide nei gallinacei. 1878. → pages17[132] DJ Alexander and IH Brown. History of highly pathogenic avian in-fluenza. Revue scientifique et technique (International Office of Epi-zootics), 28(1):19–38, 2009. → pages 17[133] WJ Bean, Y Kawaoka, JM Wood, JE Pearson, and RG Webster. Char-acterization of virulent and avirulent a/chicken/pennsylvania/83 in-fluenza a viruses: potential role of defective interfering rnas in nature.Journal of Virology, 54(1):151–160, 1985. → pages 17[134] T Horimoto, Eduardo Rivera, J Pearson, D Senne, S Krauss,Y Kawaoka, and RG Webster. Origin and molecular changes asso-ciated with emergence of a highly pathogenic h5n2 influenza virus inmexico. Virology, 213(1):223–230, 1995. → pages 17[135] Ilaria Capua and Stefano Marangon. The avian influenza epidemic initaly, 19992000: A review. Avian Pathology, 29(4):289–294, 2000. →pages 17[136] David L Suarez, Dennis A Senne, Jill Banks, Ian H Brown, Steve CEssen, Chang-Won Lee, Ruth J Manvell, Christian Mathieu-Benson,Valentina Moreno, Janice C Pedersen, et al. Recombination resultingin virulence shift in avian influenza outbreak, chile. Emerg Infect Dis,10(4):693–699, 2004. → pages 17[137] John Pasick, Katherine Handel, John Robinson, John Copps, DeidreRidd, Kevin Hills, Helen Kehler, Colleen Cottam-Birt, James Neufeld,Yohannes Berhane, et al. Intersegmental recombination between thehaemagglutinin and matrix genes was responsible for the emergenceof a highly pathogenic h7n3 avian influenza virus in british columbia.Journal of General Virology, 86(3):727–731, 2005. → pages 17[138] Guus Koch and Armin RW Elbers. Outdoor ranging of poultry:a major risk factor for the introduction and development of high-pathogenicity avian influenza. NJAS-Wageningen Journal of Life Sci-ences, 54(2):179–194, 2006. → pages 17102Bibliography[139] Kennedy F Shortridge. Poultry and the influenza h5n1 outbreak inhong kong, 1997: abridged chronology and virus isolation. Vaccine,17:S26–S29, 1999. → pages[140] LD Sims, TM Ellis, KK Liu, K Dyrting, H Wong, M Peiris, Y Guan,and KF Shortridge. Avian influenza in hong kong 1997-2002. Aviandiseases, 47(s3):832–838, 2003. → pages[141] LD Sims, J Domenech, C Benigno, S Kahn, A Kamata, J Lubroth,V Martin, and P Roeder. Origin and evolution of highly pathogenich5n1 avian influenza in asia. Veterinary Record, 157(6):159, 2005. →pages 17[142] Honglin Chen, GJD Smith, SY Zhang, K Qin, J Wang, KS Li,RG Webster, JSM Peiris, and Y Guan. Avian flu: H5n1 virus out-break in migratory waterfowl. Nature, 436(7048):191–192, 2005. →pages 18[143] F Brauer, P van der Driessche, and J Wu. Mathematical Epidemiol-ogy: Lecture Notes in Mathematics. Springer-Verlag Berlin Heidelberg,2008. ISBN 978-3-540-78910-9. doi: 10.1136/bmj.1.5082.1287-a. →pages 18[144] JAP Heesterbeek. Mathematical epidemiology of infectious diseases:model building, analysis and interpretation, volume 5. John Wiley &Sons, 2000. → pages 18[145] S Dorjee, Z Poljak, CW Revie, J Bridgland, B McNab, E Leger, andJ Sanchez. A review of simulation modelling approaches used for thespread of zoonotic influenza viruses in animal and human populations.Zoonoses and public health, 60(6):383–411, 2013. → pages 18[146] MC de Jong and TJ Hagenaars. Modelling control of avian influenzain poultry: the link with data. Revue scientifique et technique (Inter-national Office of Epizootics), 28(1):371–377, 2009. → pages 18[147] A Ssematimba. Mechanisms of Avian Influenza virus transmis-sion between farms: combining data collection and mathemat-ical modelling. Doctoral degree, Wageningen University, 2013.URL http://www.researchgate.net/publication/235641881_Mechanisms_of_Avian_Influenza_virus_transmission_between_farms_combining_data_collection_and_mathematical_modelling/file/32bfe51238408a38d1.pdf. → pages 18103Bibliography[148] Tini Garske, Paul Clarke, and Azra C Ghani. The transmissibility ofhighly pathogenic avian influenza in commercial poultry in industri-alised countries. PLoS One, 2(4):e349, 2007. → pages 18[149] A Mannelli, N Ferre, and S Marangon. Analysis of the 1999–2000highly pathogenic avian influenza (h7n1) epidemic in the main poultry-production area in northern italy. Preventive veterinary medicine, 73(4):273–285, 2006. → pages 18[150] Jennifer H McQuiston, Lindsey P Garber, Barbara A Porter-Spalding,John W Hahn, F William Pierson, Sherrilyn H Wainwright, Dennis ASenne, Thomas J Brignole, Bruce L Akey, and Thomas J Holt. Eval-uation of risk factors for the spread of low pathogenicity h7n2 avianinfluenza virus among commercial poultry farms. Journal of the Amer-ican Veterinary Medical Association, 226(5):767–772, 2005. → pages18[151] Kieran J Sharkey, Roger G Bowers, Kenton L Morgan, Susan E Robin-son, and Robert M Christley. Epidemiological consequences of anincursion of highly pathogenic h5n1 avian influenza into the britishpoultry flock. Proceedings of the Royal Society of London B: Biologi-cal Sciences, 275(1630):19–28, 2008. → pages 18[152] A Nishiguchi, S Kobayashi, T Yamamoto, Y Ouchi, T Sugizaki, andT Tsutsui. Risk factors for the introduction of avian influenza virusinto commercial layer chicken farms during the outbreaks caused bya low-pathogenic h5n2 virus in japan in 2005. Zoonoses and publichealth, 54(9-10):337–343, 2007. → pages 18[153] Marian EH Bos, Michiel Van Boven, Mirjam Nielen, AnnemarieBouma, Armin RW Elbers, Gonnie Nodelijk, Guus Koch, Arjan Stege-man, and Mart CM De Jong. Estimating the day of highly pathogenicavian influenza (h7n7) virus introduction into a poultry flock basedon mortality data. Veterinary research, 38(3):493–504, 2007. → pages18[154] V Guberti, M Scremin, L Busani, L Bonfanti, and C Terregino. Asimulation model for low-pathogenicity avian influenza viruses in dab-bling ducks in europe. Avian diseases, 51(s1):275–278, 2007. → pages[155] Levan Elbakidze. Modeling of avian influenza mitigation policieswithin the backyard segment of the poultry sector. Journal of Agri-cultural and Resource Economics, pages 195–211, 2008. → pages104Bibliography[156] Shingo Iwami, Takafumi Suzuki, and Yasuhiro Takeuchi. Paradox ofvaccination: is vaccination really effective against avian flu epidemics?PLoS One, 4(3):e4915, 2009. → pages[157] JA Van der Goot, G Koch, MCM De Jong, and M Van Boven. Quan-tification of the effect of vaccination on transmission of avian influenza(h7n7) in chickens. Proceedings of the National Academy of Sciencesof the United States of America, 102(50):18141–18146, 2005. → pages[158] V Bavinck, A Bouma, M Van Boven, MEH Bos, E Stassen, andJA Stegeman. The role of backyard poultry flocks in the epidemicof highly pathogenic avian influenza virus (h7n7) in the netherlandsin 2003. Preventive veterinary medicine, 88(4):247–254, 2009.→ pages18[159] Juan Pablo Aparicio and Mercedes Pascual. Building epidemiologicalmodels from r0: an implicit treatment of transmission in networks.Proceedings of the Royal Society of London B: Biological Sciences, 274(1609):505–512, 2007. → pages 18[160] N M Ferguson, D A T Cummings, S Cauchemez, C Fraser, S Riley,A Meeyai, S Iamsirithaworn, and D S Burke. Strategies for containingan emerging influenza pandemic in Southeast Asia. Nature, 437(7056):209–214, 2005. ISSN 0028-0836. doi: 10.1038/nature04017. → pages18[161] Timothy C Germann, Kai Kadau, Ira M Longini, and Catherine AMacken. Mitigation strategies for pandemic influenza in the unitedstates. Proceedings of the National Academy of Sciences, 103(15):5935–5940, 2006. → pages[162] Rebecca F Grais, J Hugh Ellis, and Gregory E Glass. Assessing the im-pact of airline travel on the geographic spread of pandemic influenza.European journal of epidemiology, 18(11):1065–1072, 2003. → pages18[163] Shingo Iwami, Yasuhiro Takeuchi, and Xianning Liu. Avian–humaninfluenza epidemic model. Mathematical biosciences, 207(1):1–25,2007. → pages 18[164] R Breban, J M Drake, D E Stallknecht, and P Rohani. The roleof environmental transmission in recurrent avian influenza epidemics105BibliographySupporting Information. Journal of Experimental Biology, pages 1–9,2009. → pages 18, 21, 50[165] Pejman Rohani, Romulus Breban, David E Stallknecht, and John MDrake. Environmental transmission of low pathogenicity avian in-fluenza viruses and its implications for pathogen invasion. Proceedingsof the National Academy of Sciences, 106(25):10365–10369, 2009. →pages 18, 21, 47, 83[166] Scott F Dowell. Seasonal variation in host susceptibility and cyclesof certain infectious diseases. Emerging infectious diseases, 7(3):369,2001. → pages 19[167] JC Downie, V Hinshaw, and WG Laver. The ecology of influenza.isolation of type’a’influenza viruses from australian pelagic birds. TheAustralian journal of experimental biology and medical science, 55(6):635, 1977. → pages[168] Andrew W Park and Kathryn Glass. Dynamic patterns of avian andhuman influenza in east and southeast asia. The Lancet infectiousdiseases, 7(8):543–548, 2007. → pages 19[169] D Causey and S V Edwards. Ecology of avian influenza virusin birds. Journal of Infectious Diseases, 197 Suppl:S29–S33,February 2008. ISSN 0022-1899. doi: 10.1086/524991. URLhttp://www.ncbi.nlm.nih.gov/pubmed/18269325http://jid.oxfordjournals.org/content/197/Supplement_1/S29.short. →pages 19, 70[170] World Health Organization. Avian influenza in humans.http://www.who.int/influenza/human_animal_interface/avian_influenza/en/, 2015. Accessed: 2015-10-26. → pages21[171] S Krauss, D Walker, S P Pryor, L Niles, L Chenghong, V S Hinshaw,and R G Webster. Influenza A viruses of migrating wild aquatic birdsin North America. Vector borne and zoonotic diseases (Larchmont,N.Y.), 4:177–189, 2004. ISSN 1530-3667. doi: 10.1089/vbz.2004.4.177.→ pages 21, 22[172] J K Kim and N J Negovetich. Ducks: the Trojan horsesof H5N1 influenza. Influenza and other respiratory viruses,3(4):121–128, 2009. doi: 10.1111/j.1750-2659.2009.00084.x.106BibliographyDucks. URL http://onlinelibrary.wiley.com/doi/10.1111/j.1750-2659.2009.00084.x/full. → pages 21[173] L Z Garamszegi and A P Møller. Prevalence of avian influenza andhost ecology. Proceedings of the Royal Society B, 274(1621):2003–2012, August 2007. ISSN 0962-8452. doi: 10.1098/rspb.2007.0124.URL http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2275171&tool=pmcentrez&rendertype=abstracthttp://rspb.royalsocietypublishing.org/content/274/1621/2003.short. → pages 21[174] K A Herrick, F Huettmann, and M A Lindgren. A globalmodel of avian influenza prediction in wild birds: the impor-tance of northern regions. Veterinary Research, 44:42, 2013.URL http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3687566&tool=pmcentrez&rendertype=abstracthttp://www.biomedcentral.com/content/pdf/1297-9716-44-42.pdf.→ pages 21[175] B Roche and C Lebarbenchon. Water-borne transmissiondrives avian influenza dynamics in wild birds: the caseof the 20052006 epidemics in the Camargue area. Infec-tion, Genetics and Evolution, 9(5):800–805, September 2009.ISSN 1567-7257. doi: 10.1016/j.meegid.2009.04.009. URLhttp://www.ncbi.nlm.nih.gov/pubmed/19379841http://www.sciencedirect.com/science/article/pii/S156713480900077X.→ pages 21, 28, 48[176] C F Clancy, M J A O’Callaghan, and T C Kelly. A multi-scale problem arising in a model of avian flu virus in a seabirdcolony. Journal of Physics: Conference Series, 55:45–54, De-cember 2006. ISSN 1742-6588. doi: 10.1088/1742-6596/55/1/004. URL http://stacks.iop.org/1742-6596/55/i=1/a=004?key=crossref.6679b89e39a8072315c8e403dfc665d6. → pages 21[177] S Danø, P G Sø rensen, and F Hynne. Sustained oscillations in livingcells. Nature, 402(6759):320–322, 1999. ISSN 0028-0836. doi: 10.1038/46329. → pages 22[178] T Tome´ and M J de Oliveira. Role of noise in pop-ulation dynamics cycles. Physical Review E, 79(6):061128,June 2009. ISSN 1539-3755. doi: 10.1103/PhysRevE.79.107Bibliography061128. URL http://link.aps.org/doi/10.1103/PhysRevE.79.061128http://pre.aps.org/abstract/PRE/v79/i6/e061128. →pages 22[179] B Mcnamara. Theory of stochastic resonance. Physical review A, 39(9), 1989. → pages 22[180] H Gang. Physical Review Letters, (6). → pages[181] M I Dykman and P V E McClintock. What can stochastic reso-nance do? Nature, 391(6665):344, 1998. ISSN 0028-0836. doi: 10.1038/34812. URL http://www.nature.com/nature/journal/v391/n6665/abs/391344a0.html. → pages 22[182] G E Uhlenbeck and L S Ornstein. On the theory of Brownian motion.Physical review, 36(1905), 1930. URL http://journals.aps.org/pr/abstract/10.1103/PhysRev.36.823. → pages 22, 33, 50, 116[183] J a van der Goot, M C M de Jong, G Koch, and M Van Boven.Comparison of the transmission characteristics of low and highpathogenicity avian influenza A virus (H5N2). Epidemiology andinfection, 131(2):1003–1013, 2003. ISSN 09502688. doi: 10.1017/S0950268803001067. → pages 32[184] MATLAB. version 7.10.0 (R2010a). The MathWorks Inc., Natick,Massachusetts, 2010. → pages 35[185] B Øksendal. Stochastic differential equations. 2003.URL http://link.springer.com/content/pdf/10.1007/978-3-642-14394-6_5.pdf. → pages 35[186] Priscilla E Greenwood, Mark D McDonnell, and Lawrence M Ward.Dynamics of Gamma Bursts in Local Field Potentials. Neural Com-putation, pages 1–30, 2014. ISSN 0899-7667. doi: 10.1162/NECOa 00688. URL http://dx.doi.org/10.1162/NECO{_}a{_}00688. →pages 41, 47[187] MS Bartlett. Measles Periodicity and Community Size. Journal ofthe Royal Statistical Society. Series A ( . . . , 120(1):48–70, 1957. URLhttp://www.jstor.org/stable/2342553. → pages 49[188] Robert M. May and Roy M. Anderson. Population biology of infectiousdiseases: Part II, 1979. ISSN 0028-0836. → pages 49108Bibliography[189] RP Boland, Tobias Galla, and AJ McKane. Limit cycles, complexFloquet multipliers, and intrinsic noise. Physical Review E, 79(5):15, mar 2009. URL http://arxiv.org/abs/0903.5248http://pre.aps.org/abstract/PRE/v79/i5/e051131. → pages 50, 53, 55, 56, 70[190] Eric W. Weisstein. Mathieu function. From MathWorld—A Wol-fram Web Resource, 2016. URL http://mathworld.wolfram.com/MathieuFunction.html. Last visited on 13/4/2012. → pages 51, 67[191] AJ Black and AJ McKane. Stochastic amplification in an epidemicmodel with seasonal forcing. Journal of theoretical biology, 267(1):85–94, nov 2010. ISSN 1095-8541. doi: 10.1016/j.jtbi.2010.08.014. URLhttp://www.ncbi.nlm.nih.gov/pubmed/20723547http://www.sciencedirect.com/science/article/pii/S0022519310004236.→ pages 54, 55, 56, 72, 79[192] Elias M Stein and Rami Shakarchi. Fourier analysis: an introduction,volume 1. Princeton University Press, 2011. → pages 62[193] Mathematica. Mathematica 8.0, 2010. → pages 67[194] DLMF. Mathieu Functions and Hill’s Equation. From NIST DigitalLibrary of Mathematical Functions. http://dlmf.nist.gov/, Release1.0.10 of 2015-08-07. URL http://dlmf.nist.gov/. → pages 67[195] Milton Abramowitz, Irene A Stegun, et al. Handbook of mathematicalfunctions. Applied mathematics series, 55:62, 1966. → pages 67[196] May Anne Mata, Priscilla Greenwood, and Rebecca Tyson. The roleof direct and environmental transmission in stochastic avian flu recur-rence. Paper in submission, forthcoming. → pages 70, 77[197] Edward Allen. Modeling with Itoˆ stochastic differential equations, vol-ume 22. Springer Science & Business Media, 2007. → pages 81[198] T G Kurtz. Strong approximation theorems for density dependentMarkov chains, 1978. ISSN 03044149. → pages 81[199] M B Connor. A historical survey of methods of solving cubic equations.1956. URL http://scholarship.richmond.edu/masters-theses/114/. → pages 122109Appendices110Appendices ASupplementary materialsfor Chapter 3A.1 Derivation of the avian flu SDE systemFirst, we define the probability of a jump or an increment ∆ ~X asP (∆ ~X = σt+∆t − σt) = T (σ′|σ)∆t.Note that the increments of stochastic processes St, It and Vt are ∆S =St+∆t − St, ∆I = It+∆t − It, and ∆V = Vt+∆t − Vt, respectively. Then,the expected values of the increments given the transition probabilities inSection 3.2.1 areE[∆S] = −(βSNI + ρSVNV)∆t+ µ(N − S − I)∆t+ µI∆t,E[∆I] =(βSNI + ρSVNV)∆t− µI∆t− γI∆t,E[∆V ] = τI∆t+ δV∆t− ηV∆t.(A.1)Now, each increment can be expressed as the expected value of the incre-ment plus a sum of centred increments [47]. Hence, we write the incrementsas:∆S =(−β SNI − ρS VNV+ µ(N − S − I) + µI)∆t−∆Z1 + ∆Z2 + ∆Z3,∆I =(βSNI + ρSVNV− µI − γI)∆t+ ∆Z1 −∆Z3 −∆Z4,∆V = (τI + δV − ηV ) ∆t+ ∆Z5 −∆Z6.(A.2)Here the quantities ∆Zi are conditionally centred Poisson incrementswith mean zero with conditional variances that are related to the tran-sition rates. The Poisson increment ∆Z1 corresponding to infection of asusceptible individual has a conditional variance(β SN I + ρSVNV)∆t. The111A.1. Derivation of the avian flu SDE systemincrements ∆Z2, and ∆Z3 corresponding to births in susceptible class (ordeaths in recovered and infected class) respectively have conditional vari-ances µ(N−S−I)∆t, and µI∆t. On the other hand, the Poisson increment∆Z4 that corresponds to recovery of an infected individual has a conditionalvariance equal to γI∆t. Finally, the two increments corresponding to thereplication and decay of viruses ∆Z5 and ∆Z6 must have conditional vari-ances (τI + δV ) ∆t and ηV∆t, respectively. Divide (A.2) by N and NVappropriately and take ∆t→ 0 to obtaindS =(−β SNI − ρS VNV+ µ(N − S − I) + µI)dt− dZ1 + dZ2 + dZ3,dI =(βSNI + ρSVNV− µI − γI)dt+ dZ1 − dZ3 − dZ4,dV = (τI + δV − ηV ) dt+ dZ5 − dZ6.(A.3)Suppose we replace the Poisson increments in (A.2) by multiples ofWiener increments, i.e. ∆Zi → gi∆Wi, with same standard deviations asthe Poisson increments they replace. By doing the same limiting process∆t→ 0, we obtain the stochastic differential equations (SDE):dS =(−β SNI − ρS VNV+ µ(N − S − I) + µI)dt− g1dW1 + g2dW2 + g3dW3,dI =(βSNI + ρSVNV− µI − γI)dt+ g1dW1 − g3dW3 − g4dW4,dV = (τI + δV − ηV ) dt+ g5dW5 − g6dW6,(A.4)whereg1 =√βSNI + ρSVNV, g2 =√µ(N − S − I), g3 =√µI,g4 =√γI, g5 =√τI + δV , and g6 =√ηV .(A.5)Furthermore, we can re-write (A.4) by expressing the host and viruspopulations as proportions rather than absolute numbers, i.e.s =SN, i =IN, v =VNV, and k =NNV.112A.2. Stochastic linearizationThe corresponding SDEs for the proportions of ducks and virus are thengiven by (3.2).The approximation (A.4) is an example of a result of Kurtz [10]. Analternate approach is to use a van Kampen [9] system-size expansion of theKolmogorov (Master) equation, see e.g. in Baxendale and Greenwood [50].A.2 Stochastic linearizationIn matrix notation, (3.2) can be written as:dx = F(x(t)) dt+ DG(x(t)) dW (A.6)where D = diag( 1√N, 1√N, 1√NV), dW(t) = (dW1, dW2, dW3, dW4, dW5, dW6)T ,and G(x(t)) =−G1 G2 G3 0 0 0G1 0 −G3 −G4 0 00 0 0 0 G5 −G6 .Note that x = (s, i, v) which depends on N and NV and limN,NV→∞F(x)is a vector whose components are the right-hand side of (3.4). It has beenpointed out by Allen et al. [48] that one can construct a stochastic system,which is the same in distribution such that all matrices in the diffusion termof (A.6) are square matrices whose sizes are equal to the dimension of x, i.e.in this case, a matrix C ∈ R3×3 such that (A.6) would be equivalent in lawto the stochastic systemdx˜ = F(x˜(t)) dt+ DC(x˜(t)) dW˜. (A.7)The Wiener processes W˜ ∈ R3×1 and W ∈ R6×1 both have independentterms. Moreover, the stochastic processes x˜ in (A.7) are different from theoriginally defined stochastic processes found in (A.6) but it can be shownthat their stochastic paths are the same. Thus, x˜ can be replaced by theR3-valued stochastic process x that is considered originally. Matrices G andC are related through the 3× 3 matrix V, where V = GGᵀ and C = V1/2.An explicit computation of V confirms that it is the general form for thenoise covariance matrix B that was described in Wang et al. [1]. In otherwords,V(x, t) =βsi+ ρsv + µ(1− s) −βsi− ρsv − µi 0−βsi− ρsv − µi βsi+ ρsv + (µ+ γ)i 00 0 kτi+ δv + ηv .(A.8)113A.2. Stochastic linearizationLetting N,NV → ∞ so that s → φ1, i → φ2, and v → ψ and t → ∞ wehave φ1 → φ∗1, φ2 → φ∗2, and ψ → ψ∗ implies that limN,NV ,t→∞V = Bwhich is a constant matrix whose entries are displayed as follows wherexeq ≡ (φ1 = φ∗1, φ2 = φ∗2, ψ = ψ∗), the equilibrium state of the deterministicsystem:B =B11 B12 0B21 B22 00 0 B33 , whereB11 = βφ∗1φ∗2 + ρφ∗1ψ∗ + µ(1− φ∗1),B12 = B21 = −βφ∗1φ∗2 − ρφ∗1ψ∗ − µφ∗2,B22 = βφ∗1φ∗2 + ρφ∗1ψ∗ + (µ+ γ)φ∗2, andB33 = κτφ∗2 + δψ∗ + ηψ∗.(A.9)It remains to show that the set of Langevin equations obtained by Wanget al. [1] can be constructed from the linear stochastic differential equations(the tilde in (A.7) is dropped for brevity)dx = F(x(t)) dt+ DC(x(t)) dW. (A.10)Recall that the diagonal matrix D is given in (A.6) and C = V1/2 wherethe entries of V is described in (A.8). The system (A.10) with the stochasticterm generates the deterministic process (3.4), since (A.10) becomes (3.4)as N,NV →∞ which means that this term describes the average dynamicsof the processes. On the other hand, the second term is referred to as thediffusion term. It represents the variation from the average dynamics, theO(N−1/2) fluctuations of x(t) away from the deterministic process. Thediffusion term prevents a damped system from settling to an equilibriumstate.We linearize (A.10) using the substitution x(t) = xeq + Dξ(t) and obtainDdξ = F(xeq) dt+ DJ(xeq)ξ dt+ DC(xeq) dW. (A.11)The Jacobian of F(x) evaluated at xeq is denoted by J(xeq). Now, F(xeq) =0 and so simplifies (A.11), after pre-multiplying by D−1, todξ = J(xeq)ξ dt+ C(xeq) dW. (A.12)Eq. (A.12) is the Langevin (i.e. stochastic) equation in [1](See Eq.6) writtenin slightly different form. In particular, the two equations would be equiv-alent if we divide (A.12) by dt and denote A = J(xeq) and represent the114A.3. Approximate solution for linear diffusion equations in three dimensionsdiffusion term as ζ(t), i.e. Gaussian white noise with correlation function〈ζ(t), ζ(t′)T 〉 = Bδ(t−t′). In (3.9), we have A0 = J(xeq) and C0 = C(xeq).A.3 Approximate solution for linear diffusionequations in three dimensionsWe follow Baxendale and Greenwood [50] to derive the approximatesolution for our example where the diffusion processes have values in R3.Consider the stochastic systemdξ = A0ξ dt+ C0 dW, ξ(t),W(t) ∈ R3,A0,C0 ∈ R3×3. (A.13)where,A0 =−βφ∗2 − ρψ∗ − µ −βφ∗1 −ρφ∗1−βφ∗2 βφ∗1 − µ− γ ρφ∗10 κτ δ − η , (A.14)andC0 =βφ∗1φ∗2 + ρφ∗1ψ∗ + µ(1− φ∗1) −βφ∗1φ∗2 − ρφ∗1ψ∗ − µφ∗2 0−βφ∗1φ∗2 − ρφ∗1ψ∗ − µφ∗2 βφ∗1φ∗2 + ρφ∗1ψ∗ + (µ+ γ)φ∗2 00 0 κτφ∗2 + δψ∗ + ηψ∗1/2 .(A.15)Here W(t) contains independent Wiener processes (or Brownian motion).Suppose that A0 has eigenvalues −ζ and −λ± iω for ζ, λ, ω ∈ R+. Onecan find a matrix Q ∈ R3×3 such thatQ−1A0Q = Λ ≡−ζ 0 00 −λ ω0 −ω −λ . (A.16)The matrix Λ is called the real block diagonal form of the eigenvalue matrixof A0 so it follows that Q is the real block diagonal form of the associatedmatrix of eigenvectors. By pre-multiplying (A.13) with Q−1 and using thesubstitution y(t) = Q−1ξ(t),we havedy = Λy dt+ Q−1C0 dW. (A.17)Let Σ = Q−1C0 and denote Σ•j and Σi• as its jth column vector andith row vector, respectively. With y = [y1, y2, y3]ᵀ, we write (A.17) asdy1 = −ζy1 dt+ Σ1• dW, (A.18a)dy˜ = Λ˜y˜ dt+ Σ˜ dW, (A.18b)115A.3. Approximate solution for linear diffusion equations in three dimensionswhere y˜ = [y2, y3]ᵀ, Λ˜ =[−λ ω−ω −λ], and Σ˜ = [Σ2•,Σ3•]ᵀ .Now, using the result of [48], we find that the SDE (A.18a) is equivalent tody1 = −ζy1 dt+ σ1 dW1, (A.19)where σ21 = Σ1•Σᵀ1• is the variance of the stationary distribution of y1(t)and W1(t) is a one-dimensional Wiener process. It is apparent that (A.19)describes an Ornstein-Uhlenbeck process [182] in one dimension with a sta-tionary variance σ21/2ζ. The square-root of this stationary variance corre-sponds to the standard deviation typically observed in the process y1(t).On the other hand, (A.18b) is equivalent to:dy˜ = Λ˜y˜ dt+ C˜ dW˜, (A.20)where C˜ = (Σ˜Σ˜ᵀ)1/2 and W˜(t) is a two-dimensional Wiener process. Theapproximate solution of (A.20) is related to a two-dimensional Ornstein-Uhlenbeck process as proven by Baxendale and Greenwood [50]. The ap-proximation is reasonable under the assumption that λ  ω. Thus, if itis assumed that λ  ω then by the theorem of Baxendale and Greenwood[50], the approximate solution for y˜ is:y˜ ≈ y˜app = σ˜√λR−ωtSλt, (A.21)whereσ˜2 =12Tr(C˜C˜ᵀ). (A.22)Thus,ξ ≈ ξapp ≡ Qyapp = y1Q•1 + yapp2 Q•2 + yapp3 Q•3= y1Q•1 + [Q•2,Q•3]y˜app.(A.23)More precisely,ξapp(t) = y1(t)Q•1 +σ˜√λ[Q•2,Q•3]R−ωtSλt. (A.24)Now, we know that in polar coordinatesR−ωtSλt =[cosωt sinωt− sinωt cosωt] [S1(λt)S2(λt)]=[S1(λt) cosωt+ S2(λt) sinωt−S1(λt) sinωt+ S2(λt) cosωt].(A.25)116A.3. Approximate solution for linear diffusion equations in three dimensionsUsing the formula x cos θ+y sin θ = z cos(θ − ϕ) where z = |x+ iy| and ϕ =arg(x+ iy) for x, y ∈ R and i = √−1, we then write S1(λt) = z(λt) cosϕ(λt)and S2(λt) = z(λt) sinϕ(λt) with z(λt) =√S21 + S22 ≡ |S(λt)| and ϕ(λt) =tan−1(S2/S1) to obtainR−ωtSλt = z(λt)[cos (ϕ(λt)− ωt)sin (ϕ(λt)− ωt)]≡ |S(λt)|[cos (ϕλt − ωt)sin (ϕλt − ωt)]. (A.26)Applying (A.26) to the second term of (A.24) yieldsξapp1 (t)ξapp2 (t)ξapp3 (t) = y1(t)q11q21q31+ σ˜√λ|S(λt)|q12 cos (ϕλt − ωt) + q13 sin (ϕλt − ωt)q22 cos (ϕλt − ωt) + q23 sin (ϕλt − ωt)q32 cos (ϕλt − ωt) + q33 sin (ϕλt − ωt) ,(A.27)where Q = [qij ].We define qi2 = ri cos θi and qi3 = ri sin θi where ri =√q2i2 + q2i3 andθi = tan−1(q13/q12) so that the approximate fluctuation of each componenttakes the form:ξappi (t) = qi1y1(t) +σ˜√λ|S(λt)|ri cos(ϕλt − ωt− θi). (A.28)The polar form of the approximation reveals that each model componentfluctuates according to a combination of a univariate and bi-variate OUprocesses. The first term of the approximation contains a one-dimensionalOU process weighted by a scalar determined from the transformation matrixQ while the second term contains the two-dimensional OU process thatvaries slowly and λt is a quantity that influences the radius and phase of thecircular path. The stationary variance of ξappi (t) is the sum of the stationaryvariance of each term in (A.28). This means that the long-term variance ofa fluctuation isq2i1σ212ζ+ r2iσ˜2λ(A.29)Hence, the typical magnitude of ξappi (t), i.e. stationary standard deviation,isSSDi =√q2i1σ212ζ+ r2iσ˜2λ. (A.30)Note that the fluctuation of each component i has a constant phase shift θi,which is useful in computing phase differences between disease components.117A.4. Additional insight from the approximation on the interaction of disease componentsA.4 Additional insight from the approximationon the interaction of disease componentsIn Section 3.4.1, we showed that for the given set of avian flu parametervalues, the system exhibits noise-sustained oscillations which can be viewedas a sum of two processes given by (3.25): (i) a process proportional to theone-dimensional OU process and (ii) a process proportional to the productof a rotation matrix and a standard OU process.We describe here, using (3.13), the behaviour of the sample path inthree-dimensional space. The first term of (3.13) means that sample pathbehaves as an Ornstein-Uhlenbeck process y1(t) that travels along the axisthat points to the direction of Q•1, i.e. the eigenvector associated to −ζ. Inaddition, the second term of (3.13) implies that the sample path cycles onthe subspace spanned by the last two column vectors of the transformationmatrix Q, i.e. eigenvectors of the eigenvalues −λ ± iω. This subspacecontains a plane whose equation (see (A.35) in Appendix A.5 for generalformulation), for our chosen set of parameters, is given by:ξS − 19.376ξI + 0.5814ξV = 0. (A.31)Figure A.1 shows the plane (A.31) and a realization of the stochasticsimulation of (3.25). In Figure A.1(a), we observe that the sample path lieschiefly on or near the plane (A.31). However, if we neglect the first termof (3.25), the dynamics of the fluctuations lie entirely on this plane (seeFigure A.1(b)). Thus, the portion of the sample path that departs from theplane is clearly due to the one-dimensional OU process whereas the secondterm constrains the sample path to move within the plane.From (3.25), we know that the stationary standard deviation of y1(t) is1.63 which is small compared to the constant σ˜/√λ = 10.21 that appearsin the second term of the equation. Therefore, we can neglect the first termof (3.25) and show that avian flu epidemics cycle on the plane (A.31). Thiscan be achieved mathematically when the magnitude of the real eigenvalueζ is larger than λ, which means that the approximate process approachesthe hyperplane in fast manner.118A.4. Additional insight from the approximation on the interaction of disease components−303−101−20020ξ1ξ2ξ 3(a)−303−101−20020ξ1ξ2ξ 3(b)Figure A.1: (Colour online) A sample path of the approximate fluctuationsgiven by (3.25) when the first term is (a) not set to zero and (b) set to zero.The grey region is the plane (A.31) that lies in the subspace spanned by theeigenvectors of −0.3091± 0.8377i.In Figure A.2, we compute the magnitude ζ over combinations of β and ρvalues and found that the epidemic cycles occur primarily in the hyperplanewhen β is below 100 and ρ is high, i.e. where ζ is larger than λ.119A.4. Additional insight from the approximation on the interaction of disease componentsFigure A.2: (Colour online) Plot of ζ (left panel) and λ (right panel) asfunctions of β and ρ. The white region is where R0 < 1, i.e. noise-sustainedoscillations cannot be observed here. Default parameter values are in Ta-ble 3.1.The fact that avian flu dynamics could primarily occur in the planesuggest that under certain conditions for each transmission route, one canproject the avian flu system (3.9) onto the plane (A.31) and so simplify theanalysis. For instance, using the equation of the plane, we can write onecomponent in terms of the other and convert the three-dimensional linearavian flu SDE system (3.9) into a two-dimensional one. The possibility ofmodelling recurrent avian flu epidemics as stochastic system in two dimen-sions must therefore be explored.120A.5. The subspace where the cycling takes placeA.5 The subspace where the cycling takes placeFor the case when the stationary standard deviation (s.s.d.) of the secondterm is very large compared to the s.s.d. of the first term in our approxi-mation, the first term of (A.24) is negligible and we expect the stochasticpath to primarily lie in a plane, i.e. a subspace of R3, spanned by the lasttwo column vectors of Q (Q•2 and Q•3). Here we show a general way forcomputing the equation of this plane.The sample path defined by the fluctuations ξi(t) is centred at (0, 0, 0) andso the equation of the plane should take the form:a1ξ1 + a2ξ2 + a3ξ3 = 0 (A.32)We know that the vectors Q•2 and Q•3 span the plane, hence must satisfy(A.32).Therefore,[a1, a2, a3] · [Q•2,Q•3] = 0. (A.33)By Gaussian elimination or by manipulating the explicit form of the linearsystem, we can eliminate a3 in (A.33) and a little algebra turns (A.33) intoa simpler equation,det(M1)a1+det(M2)a2 = 0 where M1 =[q12 q13q32 q33]and M2 =[q22 q23q32 q33].(A.34)Now we require det(M2) 6= 0 so that a2 = −det(M1)det(M2)a1, which givesa3 =a1q32(det(M1)det(M2)q22 − q12). Therefore, assuming that a1 6= 0, the de-sired equation of the plane isξ1 − det(M1)det(M2)ξ2 +1q32(det(M1)det(M2)q22 − q12)ξ3 = 0. (A.35)A.6 Derivation of the explicit form of themean-field eigenvaluesOur starting point is the Jacobian evaluated at the stable endemic equi-librium point [1], i.e.,J(xeq) =−βφ∗2 − ρψ∗ − µ −βφ∗1 −ρφ∗1−βφ∗2 βφ∗1 − µ− γ ρφ∗10 κτ δ − η , (A.36)121A.6. Derivation of the explicit form of the mean-field eigenvalueswhere φ∗1 =1R0 , φ∗2 =µµ+ γ(1− 1R0), and ψ∗ =κµτ(η − δ)(µ+ γ)(1− 1R0)for the basic reproduction numberR0 = βµ+ γ+κρτ(η − δ)(µ+ γ) .The condition R0 > 1 must be satisfied for the non-trivial steady-statexeq = (φ∗1, φ∗2, ψ∗) to exist.The eigenvalues of J(xeq) determine the local dynamics of the deter-ministic SIR-V system close to the non-trivial equilibrium point xeq. Now,denote the eigenvalues of J(xeq) as ν. It follows that |J(xeq)−νI| = 0 givesrise to a cubic polynomial of the formν3 − aν2 − bν − c = 0, (A.37)where:a = (δ − η)− µR0 − γ − µ+ β/R0b = −µ(η − δ + γ + µ)R0 + µβ/R0c = −µ(η − δ)(γ + µ)(R0 − 1).(A.38)Equations (A.37) and (A.38) in fact appeared in the Appendix sectionof [1], where it was proven that the endemic equilibrium is stable. Now,substitute ν = y +a3to yield the normal form transformation,y3+py+q = 0 where p =13(−a2−3b) and q = 127(−2a3−9ab−27c).(A.39)This method is also known as Vieta’s subsitution [199]. Equation (A.39)has been well-studied and has known solutions in general form:y1 = Y+ + Y−,y2 = −12(Y+ + Y−) + i√32(Y+ − Y−),y3 = −12(Y+ + Y−)− i√32(Y+ − Y−),(A.40)where: Y± =(−q2±√q24+p327)1/3and i =√−1.We are interested in the case when all three roots exist with two of thembeing complex conjugates. This is satisfied by assuming thatq24+p327> 0.122A.6. Derivation of the explicit form of the mean-field eigenvaluesBy back substitution, the solutions of (A.37) are:ν1 = Y+ + Y− +a3,ν2 = −12(Y+ + Y−) +a3+ i√32(Y+ − Y−),ν3 = −12(Y+ + Y−) +a3− i√32(Y+ − Y−).(A.41)The eigenvalues ν2 and ν3 are conjugate pairs whose real part is negativeas confirmed by Wang et al. [1]. The magnitude of the real and imaginarypart corresponds to the decay rate λ and the intrinsic frequency ω, respec-tively, of the deterministic system linearized at the endemic equilibriumstate. Therefore,λ =∣∣∣∣−12(Y+ + Y−) + a3∣∣∣∣ ,ω =∣∣∣∣∣√32(Y+ − Y−)∣∣∣∣∣ .(A.42)Using avian flu parameters in Table 3.1, we can write λ and ω in termsof β and ρ, as follows:λ(ρ, β) ≈∣∣∣− 2.9− 0.0172β − 0.5945ρ+ β0.5172β + 17.84ρ− 0.5(3√F(ρ, β) + G(ρβ) + 3√F(ρ, β)− G(ρ, β)) ∣∣∣, (A.43)ω(ρ, β) ≈√32∣∣∣− 3√F(ρ, β) + G(ρ, β) + 3√F(ρ, β)− G(ρ, β)∣∣∣, (A.44)whereF(ρ, β) =−0.5P1(ρ, β)(0.1724β + 5.945ρ)3,G(ρ, β) = 0.56√81P1(ρ, β)2(0.1724β + 5.945ρ)6+12P2(ρ, β)3(0.1724β + 5.945ρ)6,P1(ρ, β) =6∑i=06∑j=0Mi+1,j+1βiρj ,P2(ρ, β) =6∑i=06∑j=0Ni+1,j+1βiρj .(A.45)123A.6. Derivation of the explicit form of the mean-field eigenvaluesHereM =0 0 0 9190 3152 −646 88.30 0 236 311 −119 15.4 00 1.39 13.5 −8.33 1.11 0 0−0.008 0.306 −0.284 0.043 0 0 00.003 −0.005 0.001 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0(A.46)andN =0 0 −892 183 −37.5 0 00 −19 23 −4.35 0 0 0−0.135 0.871 −0.189 0 0 0 00.01 −0.004 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0. (A.47)Additionally, the first eigenvalue ν1 < 0 (eigenvalue with largest negativereal part) and so the decay rate of the OU process y1(t) in (3.13) is ζ = |ν1|.124Appendices BSupplementary materialsfor Chapter 4B.1 Floquet TheoryWe present here the formal statement of the theory preceded by def-initions of essential terminologies [52]. One can refer to the notes of Dr.Michael J. Ward for examples and detailed discussion of the theory (see http://www.emba.uvm.edu/~jxyang/teaching/Floquet_theory_Ward.pdf).Definition B.1 (Principal fundamental matrix). Let x1(t), ...,xn(t) be nlinearly independent solutions ofx˙(t) = A(t)x(t). (B.1)The n×n matrix X(t) = [x1(t) · · ·xn(t)] is a fundamental matrix solution of(B.1) if X˙ = AX for a given X(t0). Moreover, if X(t0) = I then X(t) ≡ X0is the principal fundamental matrix.Definition B.2 (Floquet multipliers and exponents). For t ∈ [0, T ], de-fine B = X−1(0)X(T ), e.g. B = X0(T ). The eigenvalues ρ1, . . . , ρn of Bare called the Floquet multipliers of X˙(t) = A(t)X(t) or equivalently, of(B.1). The Floquet exponents are µ1, . . . , µn satisfying ρj = eµjT . Notethat ρj , µj ∈ C.Theorem B.3 (Floquet theory). Let ρ be a Floquet multiplier for (B.1)and let µ be the corresponding Floquet exponent so that ρ = eµT . Thenthere exists a solution x(t) of (B.1) such that x(t + T ) = ρx(t) for all t.Further, there exists an n-dimensional T -periodic function p(t) such thatx(t) = eµtp(t).125B.2. Itoˆ formula for multi-dimensional processAnalogous to the interpretation of the eigenvalues in continuous- ordiscrete-time systems with A(t) = A0 a constant matrix, the long-timebehaviour of the solutions to (B.1) can be predicted using both the Floquetexponents µj and multipliers ρj . Below we enumerate the possible solutionbehaviours:1. If Re(µj) < 0, or equivalently |ρj | < 1, for all j, then x(t) → 0 ast→∞.2. If Re(µj) > 0, or equivalently |ρj | > 1, for some j, then x(t) → ∞ ast→∞.3. If Re(µj) = 0, or equivalently |ρj | = 1, for some j, then we have apseudo-periodic solution.In application, Floquet exponents are used to describe the growth rate ofvarious perturbations over a cycle on the averageCorollary B.4 (Lyapunov transformation). Under the transformation x(t) =P (t)y(t), which is invertible and periodic, the periodic system given by (B.1)yields a time-invariant system.B.2 Itoˆ formula for multi-dimensional processTheorem B.5. Suppose a multi-dimensional process X satisfiesdX = u dt+ V dW ,where X ∈ Rn,u ∈ Rn, V ∈ Rn×n,W ∈ Rn. Let Y (t) = (g1(t,X(t)), . . . , gp(t,X)(t)).Then the process Y (t) satisfiesdYk(t) =∂gk∂t(t,X) dt+∑i∂gk∂xi(t,X) dXi +12∑i,j∂2gk∂xi∂xj(t,X) d〈Xi, Xj〉,where 〈Xi, Xj〉 denotes the bracket process. When the SDEs dX = u dt +V dW is linear in X, the last term in the Itoˆ formula above vanishes since∂2gk∂xi∂xj= 0.126


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items