An Agent-Based Model to study the epidemiological and evolutionary dynamics of Influenza viruses
© Roche et al; licensee BioMed Central Ltd. 2011
Received: 14 December 2010
Accepted: 30 March 2011
Published: 30 March 2011
Influenza A viruses exhibit complex epidemiological patterns in a number of mammalian and avian hosts. Understanding transmission of these viruses necessitates taking into account their evolution, which represents a challenge for developing mathematical models. This is because the phrasing of multi-strain systems in terms of traditional compartmental ODE models either requires simplifying assumptions to be made that overlook important evolutionary processes, or leads to complex dynamical systems that are too cumbersome to analyse.
Here, we develop an Individual-Based Model (IBM) in order to address simultaneously the ecology, epidemiology and evolution of strain-polymorphic pathogens, using Influenza A viruses as an illustrative example.
We carry out careful validation of our IBM against comparable mathematical models to demonstrate the robustness of our algorithm and the sound basis for this novel framework. We discuss how this new approach can give critical insights in the study of influenza evolution.
The ecology and evolution of influenza viruses, especially in human and avian populations, have received considerable attention over the last decade [1–3]. It is increasingly recognised that because of the genetic diversity of these RNA viruses [4, 5] and their high mutation rate , transmission dynamics and evolution need to be considered simultaneously. The evolutionary dynamics of influenza viruses are complicated the frequent occurrence of reassortment events [6, 7], whereby two virus strains co-infecting the same individual exchange genetic material, resulting in a third strain. Despite its demonstrated frequent occurrence [1, 6], reassortment remains largely absent from theoretical studies of influenza evolution. In this respect, it is important to develop an appropriate theoretical framework to explore its evolutionary impact.
The most commonly used approach to analyse influenza dynamics has been mathematical models where state space is expanded to take into account multiple virus strains [8, 2, 9, 3, 10]. The overlap between epidemiological and evolutionary time scales , the observed immune drift in human populations [2, 3], and the high mutation rate  make it essential to integrate a large strain space into models. Consequently, the classic modelling approaches, based on the familiar SIR framework [11–13], are not readily amenable for this purpose because the resulting state space increases exponentially with the number of strains and, consequently, becomes too cumbersome for meaningful analysis. One possibility, recently suggested [14, 3], is not to track co-infections, yielding a model whose number of state variables increase linearly with the number of strains. However, these models by necessity overlook the contribution of reassortment events (but see ).
It is possible to fill this void by developing an Individual-Based Model (IBM, see [16, 17]). The rationale of this kind of modelling is to explore individual heterogeneity, allowing us to track infection history. Extensively used in ecology  and more recently in epidemiology , the main goal of these models have been to focus on very detailed and specific scenarios, with the goal of making quantitative predictions especially about alternative infection control decisions [19, 20, 9]. IBMs have also been shown to be useful in theoretical studies [21–23].
Here, we develop a stochastic IBM to simultaneously address the ecology and evolution of avian influenza viruses. We first describe the current paradigms of influenza mathematical modelling, which drives the structure of our computational framework. After detailing our IBM, we demonstrate that under analogous conditions, our model precisely recaptures the output of classic models. This validation step will allow the future use of our IBM to rely on the background of mathematical epidemiology [11–13]. Finally, we conclude that our model improves the current mathematical models and can be exploited to study evolutionary processes previously intractable, giving the opportunity to address all the components of influenza evolution.
1.1 Transmission process and mean field theory
The goal our of IBM is to extend the capabilities of the mathematical models of Influenza viruses. Here, we describe the current paradigms about the Influenza transmission and how they are modelled using the classic SIR framework. These models underlie the IBM design presented during the next section.
In humans populations, a direct transmission is assumed between individuals as usual when airborne propagation is documented . The situation is more complex in birds. Generally, direct transmission is assumed to encapsulate both fecal/oral and airborne transmission because they are occur on the same time scale and both rely on the proximity of susceptible and infectious individuals . There is increasing evidence for indirect influenza transmission in avian species through an environmental reservoir [26–29]. Experimental studies have shown the long-term persistence of these viruses in aquatic environments [4, 30] and theoretical studies have underlined its importance in epidemiological and evolutionary outcomes [29, 25, 31, 32].
Epidemiological models of influenza viruses generally consider only a limited number of strains, perhaps as few as two [33, 31]. The host population is then classified according to their infectious status with regards to each strain. For a two strain model, for example, individuals are born susceptible to both strains (N SS ) and become infected (N IS or N SI ) at rate λ i . While infected with strain i, they may be exposed to strain j resulting in co-infection (N II ), at a rate (1 - σ) λ j where σ captures the extent of cross-protection. Individuals recover (N RS or N SR ) at rate γ i and are assumed to be invulnerable to that strain.
where σ ij represents the cross-immunity network between strains, m is the mutation rate and i and j are the strain "identities" of the pathogens. Here, we assume that cross-immunity decreases exponentially  with the distance between the infection history and the challenging strain (equation 19) as it has been described for influenza viruses in humans  and horses . This exponential decrease is shaped by the parameter d .
where R j is the rate of event j and RAND 1 is a uniform random deviate in (0, 1). Then, a new uniform random number in (0,1) is generated and multiplied by ∑ n R j , denoted by P . The next event is then determined by ∑j-1< P < ∑ j .
The recent development in mathematical models of influenza viruses fails to follow co-infections and consequently to incorporate reassortment. We choose the framework of IBM (Individual-Based Model, see ) to develop a model fully considering epidemiology and evolution of Influenza viruses.
2.1 Model structure
Parameters of the model
Class in IBM
Value used in this study
Variables of the model
Array with every strains in Infectious state
Array with every strains just enter in Infectious state
Array with infection history
Array containing every strains present in the environment
Viral load of every stains present in the lake
Array with all "Host" objects
2.2 Dynamical Model behaviour
For the sake of simplicity, here we do not explain every functions in detail, but focus instead on the most important ones (all algorithms are detailed in the section S1 in additional file 1). Specifically, we explain the functions where infection, mutation or recovery events are applied. The probability of these events are calculated by applying an exponential distribution on rates, as usual in mathematical epidemiology (P = 1 - exp(- rate × δt) where δt is the model integration time step and rate is the individual rate expressed in the mean field approximation ).
2.2.1 Program scheduler and host status
In IBM, two steps are generally considered at each iteration . The first one determines the status of each individual at the end of this time step. Then, when all individuals have computed their next status, they are updated through a second function. This methodology is applied to avoid synchronization problems which can occur where an individual calculates its status based on other individuals status who have one more time step.
In our case, at each iteration, a function "Step" is called for every individual. If a host becomes infected to one specific strain during this function, a "Host" attribute ("nextPathogens") will be filled with a memory reference to this "Pathogen" object. After all individuals have called the "Step" function, the "Update" function moves this reference to another variable ("currentPathogens"). During the following time step, this reference can move to another variable ("oldPathogens") with a probability P = 1 - exp(- γδ t ).
2.2.2 Infection process
Within the "Step" function, transmission process (direct and/or environmental) takes place. For each strain present in the system at this time step, a probability of infection is computed for both transmission routes. In the case of direct transmission, this probability is equal to P i = (1 - σ)(1 - exp(- βI i δ t )). This probability is added to which describes the environmental transmission (parameters have the meaning than in equations 1-13).
The mathematical formulation assumes that co-infections occur sequentially and not simultaneously. In fact, co-infections can be approximately simultaneous because of the very small time step involved in Gillespie's method. Time step duration is not infinitesimal in the IBM. To avoid an under-estimation of co-infections, we have to consider that co-infections occur sequentially, but potentially during the same time step. Then, we apply a method already used in . We rank randomly infection probabilities for each individual and apply the transition if selected (if RAND 1() < P). Ranking randomly the infections events removes the potential issue that pathogens at the top of the pathogen list can have a greater chance to infect the host.
Upon infection, mutations occurs at a fixed rate. At the end of the host update, a function "mutation" is called for every "Pathogen" objects found in the list "currentPathogens". The probability that a mutation occurs is given by P = 1 - exp(- mδt). When this event is selected, the "identity" will increase or decrease by 1 with a probability of 50% each. Hence, we assume a linear strain space. Other pathogen parameters stay identical (i.e. no evolution of life-history traits).
When all individuals have completed their "Update" function, host demography is applied. Hosts have an expected period before producing one o spring as well as an average lifespan. Each host produces one off spring (susceptible) with a probability P = 1 - exp(- μδt) and dies with the same probability (assuming a constant population size).
2.2.5 Viral demography within environment
2.3 Input and Output Files
The attributes which have to filled by the user are split between three different input files in order to distinguish their class destination. The "main" input file contains the global parameters of the simulation, i.e. the length of a time step, the lake volume or the filenames of other input files. The "species" file describes the ecology of the population considered. It integrates its birth, death and drinking rate, its population size and which transmission route has to be included. Finally, the "pathogen" file depicts the life-history traits of the introduced pathogens, i.e. its mutation and direct-transmission rate, the viral load required to yield 50% of infections, as well as its initial conditions within the host population (identity, infectious population size and environmental concentration). Finally, we have implemented three different outputs which track the pathogen dynamics at the environmental (dynamics of strains concentration), individual (all the infection events) and populational level (dynamics of strain infectious population size).
3 Results and discussion
3.1 Methodology for IBM validation
We compare IBM outputs in specific cases which can be addressed by the stochastic version of the two mathematical frameworks previously exposed (equations 1-13, see also table S1 in additional file 1). The goal of this validation step is to show that our algorithms reproduce correctly the expected behaviour of SIR framework.
Within the existing epidemiological studies using an IBM, a time step of one day is generally used and claimed to be small enough. To assess this, we analyse a first simple epidemiological case where only one strain and direct transmission are involved. These results determine the optimal time step duration to produce similar disease dynamics than in SIR. It is worth to point out that the period between two time steps has to be as small as possible to get the correct dynamics and large enough to reduce significantly the simulation time.
Then, we study a situation where two identical strains and environmental transmission are introduced. We compare the IBM outputs with the stochastic version of the SIR system depicted in equations 1-13. As usual when stochastic outputs are compared , we analyse a bunch of pathogen dynamics characteristics. For both strains, we analyse (i) the epidemics peak (in terms of infectious population size), (ii) the time at epidemics peak (time needed to reach the epidemics peak), (iii) the epidemic duration (time between the start and the end of the epidemics) and (iv) the sampled epidemics size (sum of all new infections during an epidemics). We study also the global properties of the epidemics: (v) the global extinction proportion (no epidemics are observed for both pathogens) and (vi) the ratio between dominance of strain 1 and dominance of strain 2 (epidemics occurs only for strain 1 over epidemics occurs only for strain 2).
Finally, we compare the evolutionary dynamics of our model with the model assuming polarity immunity (described in equations 17-18). We study this dynamics without environmental transmission since the analysis of its influence on evolutionary dynamics is beyond the scope of this paper.
3.2 Dynamical behaviour
With only direct transmission and one strain (β2 = ω2 = 0), we explore here the influence of the IBM time step duration in order to find its optimal value
3.3 Comparison of epidemiological signatures
We now compare IBM and the stochastic SIR outputs where two identical strains and environmental transmission are involved. Here, our goal is to show that our model yields epidemiological dynamics indistinguishable from the SIR model depicted in equations 1-13.
3.4 Evolutionary dynamics
We now focus on the adequacy between evolutionary dynamics produced by our IBM and by stochastic SIR model when only density-dependent transmission is involved (described in equations 17-18, implemented as in ). We assume the same cross-immunity pattern than in equation 19.
4 Conclusions & Discussions
In this paper, we have introduced a stochastic Individual-Based Model (IBM) to study epidemiological and evolutionary dynamics of avian influenza viruses. We have shown that we can accurately reproduce the solutions of classic SIR model as a special case of our model. However, the more general set of conditions that may be represented by the IBM will enable investigation of a much larger class of evolutionary and ecological dynamics.
Additional epidemiological signatures could be analysed. Nevertheless, our goal was to give the proof that our algorithms have been implemented correctly in order to can rely on SIR literature for the future improvements of our IBM. Since six different epidemiological signatures, tackling different aspects of the disease dynamics (coexistence, extinction, epidemics peaks, etc...), have been used (Figure 3), we believe that the overall picture allows us to claim that our IBM mimics SIR models with a computational formulation at an individual scale.
The main difference between IBM and SIR models is the possibility to acquire sequentially several infections within a same time step. We had to make this assumption since the period between two time steps is necessarily bigger in the IBM than in the mathematical framework where the time step is infinitesimal. Results reported here show that, despite this difference, the two formulations are dynamically identical.
We suggest that the discrepancy between the evolutionary dynamics produced by SIR model and by our IBM (Figure 4) is due to the known underestimation of the infectious population size consecutive to the assumption of polarity immunity . Even if it is impossible to quantify this under-estimation, it is worth to point out that this dichotomy cannot be due to a time step too coarse in our IBM. Indeed, a too large IBM time step tends to underestimate infectious dynamics (as shown in Figure 3), instead of the opposite.
The possibility to track co-infections has been central in the design of this model, but it is valuable to underline that our IBM improves the current mean field theories in an additional way. Our model can cope with a strain space that is virtually infinite (constrained by the maximal number of a "double" variable, usually more than 10300). Hence, it is not necessary to invoke pragmatic approximations in order to simplify the model; the full epidemiological and evolutionary dynamics can be explored in our framework.
This IBM can be extended in numerous ways. Excepting all the possible computational additions (space, network interactions, etc..), complex ecological and evolutionary processes can also be included. So far, we have considered a constant host population size, but the modification of the demography function can integrate all kind of dynamics, even data from the field. Similarly, the reassortment process, which has never been theoretically studied as we said before, has not been exposed here because its analysis is beyond the scope of this paper. Nevertheless, it can be included easily through the mutation function. The strain identity should be replaced by a set of characters, representing amino acids for instance. Then, this string can be also divided into different parts to mimic the possibility of different genes. This modification of pathogen genome makes possible the creation of a third "Pathogen" object produced by the exchange of amino acids between strains co-infecting a given individual.
Here, we propose an Individual-Based Model improving the current modelling paradigms on influenza viruses. Its validation against a SIR model allows the future uses to rely on the SIR background. Its ability to be extended opens many new areas of influenza research which were previously constrained by the limitations of the mathematical formulation. To conclude, we shown that modelling at an individual scale allows the study of mathematically inaccessible situations despite an identical behavior. We believe this model offers an unique opportunity to fully address the evolution of influenza viruses.
5 Availability and requirements
Project name: InfluenzaIbm
Project home page: https://sites.google.com/site/rocheben/influenzaibm
Operating system(s): Platform independent
Programming language: C++
Other requirements: None
License: GNU GPL
This work was supported by the Centers for Disease Control and Prevention (5U19Cl000401), the James S. McDonnell Foundation and the National Science Foundation (DEB-0917853). PR was also supported by the RAPIDD program of the Science & Technology Directorate, Department of Homeland Security, and the Fogarty International Center, National Institutes of Health.
- Webster RG, Bean WJ, Gorman OT, Chambers TM, Kawaoka Y: Evolution and ecology of influenza A viruses. Microbiol Rev 1992, 56: 152–179.PubMed CentralPubMedGoogle Scholar
- Ferguson N, Galvani A, Bush R: Ecological and immunological determinants of influenza evolution. Nature 2003, 422(6930):428–33. 10.1038/nature01509View ArticlePubMedGoogle Scholar
- Koelle K, Cobey S, Grenfell BT, Pascual M: Epochal evolution shapes the phylodynamics of interpandemic influenza A (H3N2) in humans. Science 2006, 314: 1898–1903. 10.1126/science.1132745View ArticlePubMedGoogle Scholar
- Webster RG, Yakhno M, Hinshaw VS, Bean WJ, Murti KG: Intestinal influenza: replication and characterization of influenza viruses in ducks. Virology 1978, 84(2):268–278. 10.1016/0042-6822(78)90247-7View ArticlePubMedGoogle Scholar
- Olsen B, Munster VJ, Wallensten A, Waldenström J, Osterhaus ADME, Fouchier RAM: Global Patterns of influenza A Virus in Wild Birds. Science 2006, 312(5772):384–388. 10.1126/science.1122438View ArticlePubMedGoogle Scholar
- Holmes EC, Ghedin E, Miller N, Taylor J, Bao Y, George KS, Grenfell BT, Salzberg SL, Fraser CM, Lipman DJ, Taubenberger JK: Whole-Genome Analysis of Human influenza A Virus Reveals Multiple Persistent Lineages and Reassortment among Recent H3N2 Viruses. PLoS Biol 2005, 3(9):e300. 10.1371/journal.pbio.0030300PubMed CentralView ArticlePubMedGoogle Scholar
- Dugan VG, Chen R, Spiro DJ, Sengamalay N, Zaborsky J, Ghedin E, Nolting J, Swayne DE, Runstadler JA, Happ GM, Senne DA, Wang R, Slemons RD, Holmes EC, Taubenberger JK: The evolutionary genetics and emergence of avian influenza viruses in wild birds. PLoS Pathog 2008, 4(5):e1000076. 10.1371/journal.ppat.1000076PubMed CentralView ArticlePubMedGoogle Scholar
- Andreasen V, Lin J, Levin SA: The dynamics of cocirculating influenza strains conferring partial cross-immunity. J Math Biol 1997, 35(7):825–842. 10.1007/s002850050079View ArticlePubMedGoogle Scholar
- Ferguson NM, Cummings DAT, Cauchemez S, Fraser C, Riley S, Meeyai A, iamsirithaworn S, Burke DS: Strategies for containing an emerging influenza pandemic in Southeast Asia. Nature 2005, 437(7056):209–214. 10.1038/nature04017View ArticlePubMedGoogle Scholar
- Minayev P, Ferguson N: Improving the realism of deterministic multi-strain models: implications for modelling influenza A. J R Soc Interface 2009, 6(35):509–518.PubMed CentralView ArticlePubMedGoogle Scholar
- Kermak W, McKendrik A: A contribution to the mathematical theory of epidemics. Proc Roy Soc Lond 1927, 115: 700–726. 10.1098/rspa.1927.0118View ArticleGoogle Scholar
- Anderson RM, May RM: Infectious diseases of humans: Dynamics and control. Oxford Science Publications; 1991.Google Scholar
- Keeling MJ, Rohani P: Modeling infectious diseases in humans and animals. Princeton University Press; 2008.Google Scholar
- Gog JR, Grenfell BT: Dynamics and selection of many-strain pathogens. Proc Natl Acad Sci USA 2002, 99(26):17209–17214. 10.1073/pnas.252512799PubMed CentralView ArticlePubMedGoogle Scholar
- Day T, André JB, Park A: The evolutionary emergence of pandemic influenza. Proc Biol Sci 2006.Google Scholar
- Bonabeau E: Agent-based modeling: methods and techniques for simulating human systems. Proc Natl Acad Sci USA 2002, 99(Suppl 3):7280–7287. 10.1073/pnas.082080899PubMed CentralView ArticlePubMedGoogle Scholar
- Ferber J: Multi-agent Systems: Introduction to Distributed Artificial Intelligence. Addison-Wiley; 1999.Google Scholar
- Grimm V: Ten years of indivdal-based modelling in ecology: what have we learned and what we could learn in the future? Ecol Mod 1999, 115: 129–148. 10.1016/S0304-3800(98)00188-4View ArticleGoogle Scholar
- Segovia-Juarez J, Ganguli S, Kirschner D: Identifying control mechanisms of granuloma formation during M. tuberculosis infection using an agent-based model. J Theor Biol 2004, 231(3):357–76. 10.1016/j.jtbi.2004.06.031View ArticlePubMedGoogle Scholar
- Muller G, Grébaut P, Gouteux J: An agent-based model of sleeping sickness: simulation trials of a forest focus in southern Cameroon. C R Biol 2004, 327: 1–11. 10.1016/j.crvi.2003.12.002View ArticlePubMedGoogle Scholar
- Grimm V, Revilla E, Berger U, Jeltsch F, Mooij WM, Railsback SF, Thulke HH, Weiner J, Wiegand T, Deangelis DL: Pattern-oriented modeling of agent-based complex systems: lessons from ecology. Science 2005, 310(5750):987–991. 10.1126/science.1116681View ArticlePubMedGoogle Scholar
- Grimm V, Railsback SF: Individual-based modeling and ecology. Princeton University Press; 2005.View ArticleGoogle Scholar
- Roche B, Guégan J, Bousquet F: Multi agent systems and vector-borne disease epidemiology. Tinain model: First step toward real integration. BMC BioInformatics 2009, 9: 435. 10.1186/1471-2105-9-435View ArticleGoogle Scholar
- McCallum H, Barlow N, Hone J: How should pathogen transmission be modelled? Trends in Ecology and Evolution 2001, 16(6):295–300. 10.1016/S0169-5347(01)02144-9View ArticlePubMedGoogle Scholar
- Rohani P, Breban R, Stallknecht DE, Drake JM: Environmental Transmission of Avian influenza Viruses and its Implications for Disease Control. P Natl Acad Sci USA 2009, 106(25):10365–10369. 10.1073/pnas.0809026106View ArticleGoogle Scholar
- Markwell DD, Shortridge KF: Possible waterborne transmission and maintenance of influenza viruses in domestic ducks. Appl Environ Microbiol 1982, 43: 110–115.PubMed CentralPubMedGoogle Scholar
- Laudert E, Sivanandan V, Halvorson D, Shaw D, Webster RG: Biological and molecular characterization of H13N2 influenza type A viruses isolated from turkeys and surface water. Avian Dis 1993, 37(3):793–799. 10.2307/1592031View ArticlePubMedGoogle Scholar
- Sivanandan V, Halvorson DA, Laudert E, Senne DA, Kumar MC: Isolation of H13N2 influenza A virus from turkeys and surface water. Avian Dis 1991, 35(4):974–977. 10.2307/1591638View ArticlePubMedGoogle Scholar
- Roche B, Lebarbenchon C, Gauthier-Clerc M, Chang CM, Thomas F, Renaud F, van der Werf S, Guegan JF: Water-borne transmission drives avian influenza dynamics in wild birds: the case of the 2005–2006 epidemics in the Camargue area. Infect Genet Evol 2009, 9: 800–805. 10.1016/j.meegid.2009.04.009View ArticlePubMedGoogle Scholar
- Brown JD, Goekjian G, Poulson R, Valeika S, Stallknecht DE: Avian influenza virus in water: infectivity is dependent on pH, salinity and temperature. Vet Microbiol 2009, 136(1–2):20–26. 10.1016/j.vetmic.2008.10.027View ArticlePubMedGoogle Scholar
- Roche B, Rohani P: Environmental transmission scrambles coexistence patterns for avian influenza. Epidemics 2010, 2(2):92–98. 10.1016/j.epidem.2010.03.002View ArticlePubMedGoogle Scholar
- Breban R, Drake JM, Stallknecht DE, Rohani P: The role of environmental tranmission in reccurent avian influenza dynamics. PLoS Comput Biol 2009, 5: e1000346. 10.1371/journal.pcbi.1000346PubMed CentralView ArticlePubMedGoogle Scholar
- Breban R, Drake JM, Rohani P: A general multi-strain model with environmental transmission: Invasion conditions for the disease-free and endemic states. JTB 2010, 264: 729–736. 10.1016/j.jtbi.2010.03.005View ArticleGoogle Scholar
- Rohani P, Wearing H, Huang Y: Infectious Disease Ecology. Princeton University Press 2008 chap. Understanding Host-Multi-Pathogen Systens: The Interaction Between Ecology and Immunology;Google Scholar
- Park AW, Daly JM, Lewis NS, Smith DJ, Wood JLN, Grenfell BT: Quantifying the impact of immune escape on transmission dynamics of influenza. Science 2009, 326(5953):726–728. 10.1126/science.1175980PubMed CentralView ArticlePubMedGoogle Scholar
- Gillespie DT: Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry 1977, 81: 2340–2361. 10.1021/j100540a008View ArticleGoogle Scholar
- Gillespie DT: Approximate accelerated stochastic simulation of chemically reacting systems. Journal of Chemical Physics 2001, 115: 1716–1733. 10.1063/1.1378322View ArticleGoogle Scholar
- Stroustrup B: C++ Programming Language. 3rd edition. Addison-Wesley Professional; 1997.Google Scholar
- Treuil JP, Drogoul A, Zucker JD: Modélisation et simulation à base d'agents. Exemples commentés, outils informatiques et questions théoriques. Dunod. 2008.Google Scholar
- Keeling MJ: The ecology and evolution of spatial host-parasite systems. PhD thesis. University of Warwick; 1995.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.