Mathematical Biosciences 188 (2004) 17–28
www.elsevier.com/locate/mbs

Evolution and persistence of inﬂuenza A and other diseases
Simon A. Levin *, Jonathan Dushoﬀ, Joshua B. Plotkin
Department of Ecology and Evolutionary Biology, Princeton University, 203 Eno Hall, Princeton, NJ 08544 1003, USA
Received 17 December 2002; received in revised form 13 August 2003; accepted 13 August 2003

Abstract
The evolution of the etiological agents of disease presents one of the greatest challenges for their control,and makes essential complementing standard epidemiological investigations with broader approaches thatallow for evolutionary change. Given the stunning genetic diversity that is possible for many such agents,
such as the inﬂuenza virus, it is impossible to represent all of the diversity manifest at the level of amino acidsequences. We show that drift-variant inﬂuenza strains naturally cluster into groups which are associatedwith functionally important epitopic regions. Dominant clusters typically replace each other every 2–5
years, and this feature is fundamental to the development of vaccination strategies. We furthermore showthat stochastic ﬂuctuations can greatly magnify small interference eﬀects among strains, or even among
subtypes, leading for example to competitive exclusion in situations where such eﬀects would be unexpectedbased on the usual deterministic models. We suggest that this eﬀect might be involved in the explanations ofsome persistent empirical anomalies.
 2003 Elsevier Inc. All rights reserved.
Keywords: Disease evolution; Inﬂuenza; Hemagglutinin; Stochastic models

1. Introduction
The classical theory of epidemics is one of the jewels of mathematical biology, elegant in its
form and providing insights and guidance into the management of infectious diseases. The basicmodel (Fig. 1) divides the host population into classes – typically susceptible (S), infectious (I) andremoved (recovered and immune, R) – according to their disease status, and develops a set of
diﬀerential equations to track changes in abundance (see [2]).

*Corresponding author. Tel.: +1-609 258 6880; fax: +1-609 258 6819.
E-mail address: simon@eno.princeton.edu (S.A. Levin).

0025-5564/$ - see front matter  2003 Elsevier Inc. All rights reserved.
doi:10.1016/j.mbs.2003.08.010

18

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

Fig. 1. A diagram of the S-I-R ÔboxÕ model. Individuals become infected (move from S to I) at a rate proportional to thenumber of infectious individuals in the population, the only non-linear process in the model. Additionally new individuals are born into the susceptible class, and all individuals die at rate l. The basic reproductive number for thismodel is given by R0 ¼ b=ðl þ cÞ. Since the population size is constant, the equation for R is redundant, and is omitted.

In the classical S-I-R model, with constant population size, no sustained oscillations are possible. Depending on the parameters and population size, the disease either dies out, or is sustainedat an endemic level. If the duration of immunity is long compared to the duration of infectiousness (as in most human diseases) disease prevalence shows damped oscillations. But realdiseases show persistent ﬂuctuations, spatially and temporally, on a wide variety of scales. Whatthen is lacking from the basic models, and how can they be adjusted to account for the persistenceof ﬂuctuations?
First of all, it is important to separate the various phenomena – ﬂuctuations on the scale of days
or weeks have very diﬀerent explanations from those on the scale of years or generations. For
inﬂuenza A, for example, there are clear tendencies for seasonal outbreaks in the Fall, although
the actual date of emergence and spatial locus is somewhat unpredictable within the envelope of
possibilities. These patterns apply both in the Northern Hemisphere and the Southern, shifted byapproximately 6 months. The extent to which these patterns relate to the intrinsic properties of thevirus in changing climatic conditions, to the social dynamics of children entering crowded
schoolrooms, or to other factors remains a point of debate, however, as does the extent to whichsuch exogenous factors are ampliﬁed by the oscillatory tendency of the disease system. Furthermore, it is not clear whether these are self-sustaining patterns within each hemisphere, inwhich case a puzzle is to explain how and where the virus survives in between epidemics, orwhether the two hemispheres constantly re-seed each other through migrants. The truth probably
involves both mechanisms.
Over longer time scales, diseases may oscillate due to evolutionary changes, combined withchanges in the epidemiological status of the host population. Indeed, the continual evolution ofpathogens such as the inﬂuenza A virus are crucial to their survival, involving an existential gameof cat and mouse with their hosts. Understanding the role of evolution can thus be a powerful toolin controlling disease. For example, major new subtypes of the inﬂuenza A virus, typically causedby reassortment between an avian inﬂuenza virus and a human strain in a coinfected host, cause
devastating pandemics on generational time scales. The 1918 Spanish Flu, for example, caused
perhaps 40 million deaths [40] and great morbidity, largely among individuals in their twenties
and thirties. The emergence of new subtypes through reassortment is termed ÔshiftÕ, and raises a

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

19large number of theoretical and practical challenges. Fluctuations also occur over shorter timescales, i.e. from year to year, due to mutation-based ÔdriftÕ variants. It is these changes that occupy
most of the attention of those concerned with determining recommended vaccines on an annual
basis. Indeed, despite the fact that the annual epidemics associated with drift variants are muchless pronounced than when pandemics occur, their cumulative impact still greatly outweighs
pandemics [37,38].

2. Oscillations in disease models
In general, there are a wide variety of mechanisms that can lead to oscillations in disease
dynamics [31] These include:• Stochastic factors,• Seasonal forcing,• Age structure, which can introduce implicit delays,• Fluctuating population sizes, perhaps driven in part by disease dynamics,• Non-linear incidence functions, in which the number of new cases does not increase in simple
proportion to the number of susceptibles or number of infectives [27],• Spatio-temporal interactions [5,20],• Interactions among multiple diseases or strains.
Various attempts have been made to look at all of these mechanisms, and their inﬂuence upon
ﬂuctuations. Each has its own signature, and decomposing observed ﬂuctuations into those pieces
attributable to the various possible causes is a daunting challenge, made more diﬃcult by the nonlinear reinforcement of one mechanism by another. For example, Castillo-Chavez et al. [10]showed that introducing age structure alone, or two interacting strains of the inﬂuenza A virus
alone, increased the tendency to oscillate but still resulted in damped oscillations. Coupling together multiple strains with age structure is thus analogous to coupling two damped oscillators,
and was found to lead to sustained oscillations.
In more recent work, Andreasen et al. [3] found that four interacting strains, with appropriate
patterns of cross-reactivity, could sustain oscillations without any other model complexities, suchas age structure; and Lin et al. [26] showed that indeed this is possible even with as few as threestrains, one of which interacts strongly with the other two, which in turn interact weakly with eachother. Even with only three strains, the dynamical system is complicated, with eight levels ofsusceptibility (depending on what combination of strains an individual has experienced in thepast) and 12 kinds of infected individuals (distinguished by current infection and past history).
Thus, even for this simpliﬁed system, there are 20 coupled diﬀerential equations to be explored.
In all of this work, the deﬁnition of a strain is left vague. Mathematically, it refers to a genetic
unit that can be treated as homogeneous in its properties, and has deﬁnable properties of crossreactivity with other genetic units. If this unit is identiﬁed with subtypes, then the partitioning is in
general too coarse – observed patterns of cross-reactivity among cocirculating subtypes appear tobe weak to non-existent, although short-term immunological memory (or even the temporary
isolation of infected individuals) can have surprisingly strong eﬀects upon macroscopic dynamics

20

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28(see Section 4). On the other hand, if each unique RNA sequence is treated as a unique strain, one
has both a hopelessly complicated problem, with literally thousands of equations [3,18,19], andfurthermore with more detail than is appropriate to explain observed patterns. Thus, a relevant
question is whether there is an intermediate scale, involving clusters of types (or ÔquasispeciesÕ)that deﬁne the natural scale for explaining ﬂuctuations on the order of one or several years. We
turn to this question in the next section.

3. Genotype clusters in the inﬂuenza A virus
Theoretical work on the epidemiological dynamics of interacting disease types makes clear that
it is important to clarify what exactly a functional type is. As pointed out in the last section, toomuch detail can create a hopeless morass of equations, while too little detail can miss crucial
interactions. In this section, we turn our attention to determining an appropriate scale of strain
aggregation for dealing with the annual epidemics of inﬂuenza A, and with changes that occur
over a few seasons. In the next section, we look at longer time scales.
Within a given subtype, inﬂuenza A evolves at an astonishing rate, producing a large amount of
RNA sequence diversity. Considerable evidence [13,14,34] suggests that rapidly evolving RNA
viruses eﬀectively experience selection as clusters rather than as individual genotypes. Theseclusters qualify as ÔquasispeciesÕ in the broad sense of a cluster of genotypes that may be seen as aunit of selection [14], and not in the narrower sense (often used in the HIV literature) of a clusterof genotypes within a single infected individual.
The drift variants of inﬂuenza A are produced by gradual substitutions to its genome, especiallyin its hemagglutinin gene. Hemagglutinin, or HA, is the primary target of the antigens producedby the immune system to neutralize an inﬂuenza infection. Recovery from inﬂuenza infection
brings lasting immunity to the infecting genotype, but most people are susceptible to re-infection
by a new drift variant within a few years. Thus antigenic drift requires that vaccines be updated
annually to correspond to the dominant epidemic genotype(s) of HA.
The phenomenon of acquired immunity suggests that the inﬂuenza virus may be under frequency-dependent selection. Novel variants of inﬂuenza may initially enjoy a selective advantage
due to the immunologically naive pool of human hosts; but after a given variant has caused an
epidemic, it is then handicapped by the immunity it has stimulated in the human population.
Phylogenetic reconstructions of HA evolution [16,17] conﬁrm that modiﬁcations to the immunogenic parts of HA accrue at a dramatic rate. Those sites of HA involved in antigen determination also exhibit signiﬁcantly more non-synonymous nucleotide substitutions than synonymous
substitutions [8,21].
Although phylogenetic algorithms can estimate the evolutionary relationships among sequences, inﬂuenza phylogenies, even those using extensive databases, are often plagued by poorbootstrap values and instabilities of tree topology, which have been systematically studied by only
a few authors [7,9,16]. The problem of resolving an evolutionary time series to the level ofindividual sequences is thus diﬃcult, but perhaps unnecessary.
We have attempted to identify swarms of functionally related inﬂuenza genotypes by detecting
a natural length scale at which HA amino acid sequences aggregate into clusters. We have studied

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

21

a database of 560 aligned HA1 sequences, each 987 nucleotides, from inﬂuenza viruses isolated
between 1968 and 2000 [28,29,36].
In order to identify viral clusters we ﬁrst use the Hamming metric to deﬁne the distance between
two HA1 sequences. That is, we count the number of pairwise amino acids that diﬀer between twoaligned sequences. (Alternative metrics that weight the diﬀerences between amino acids according
to their stereochemical properties yield similar results.)
For each threshold distance d, we deﬁne ÔclustersÕ in the sequence database as the connectedcomponents of the graph Gd that results from placing edges between sequences within a distance dof each other. For each threshold distance d, we deﬁne the mean cluster size CðdÞ as the ﬁrst
moment of the (discrete) distribution of connected component sizes of Gd , divided by the numberof sequences. Equivalently, CðdÞ equals the probability that two randomly chosen sequences lie inthe same cluster. As we vary d we obtain a hierarchy of graphs, from the completely disconnected
graph to the fully connected graph.
We detect a natural, non-random scale of aggregration by inspecting the percolation behaviorof the mean cluster size as d varies. According to the fundamental theorem of continuum percolation [32], a large database of random sequences would exhibit a mean cluster-size curve with asingle sharp phase transition from CðdÞ ¼ 0 for small d to CðdÞ ¼ 1 for large d. Therefore, nontrivial plateaus in the cluster-size curve correspond to stable length scales at which the sequences
form non-random clusters.
In order to produce the smooth cluster-size curve in Fig. 2, for each threshold distance, hdi, we
in fact average over many realizations of the clustering (cf. [35]). In each realization, we connect

1
0.9
Mean cluster size

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

2

4

6
8
10
Threshold distance

12

14

Fig. 2. The mean cluster-size curve for 560 sequences of HA1. This curve shows the relationship between the thresholddistance d (at which to connect two sequence into the same cluster) and the resulting normalized mean cluster size CðdÞ.
Plateaus in the cluster-size curve correspond to stable length scales at which the sequences form non-random clusters.
The smooth cluster-size curve results from averaging over 100 probabilistic Gaussian draws for each mean distanceparameter d, with a 5% coeﬃcient of variation (cf. [35]). The HA1 data exhibit two signiﬁcant plateaus corresponding toclusterings at d ¼ 2–3 and d ¼ 4–5. The long tail for d P 6 corresponds to the gradual accumulation of outlier sequences. When d ¼ 2, there are 174 resulting clusters with Cð2Þ ¼ 0:0614; at this scale, the mean cluster size is 34.4
sequences.

22

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

Cluster size

40

30

20

10

1985

1990
1995
Calendar year

2000

Fig. 3. The number of HA1 sequences within each cluster plotted as a function of calendar year of isolation. The
clustering shown here corresponds to d ¼ 2 amino acids (see Fig. 1). Each cluster is indicated by a diﬀerent color, withthe eight largest clusters shown in bold. The dashed line indicates the total number of isolates in the data set each year.
The dominant sequence clusters tend to replace each other every 2–5 years. The dominant cluster in each year accounts
for more than 25% of the sequences isolated that year. (The number of sequences each year does not reﬂect the severityof infections, but rather temporal biases in the sequence data set.)

two sequences by an edge if they lie within a randomly chosen distance d, with mean hdi and 5%
coeﬃcient of variation. (The connection distance d is chosen from a truncated normal distribution.) Thus the cluster-size curve in Fig. 2 is, in essence, the average over many slight perturbations of the deterministic clustering method described above.
Fig. 2 shows that there is a natural, non-random partitioning of the HA sequences into disjoint
clusters at a threshold distance of d ¼ 2 amino acid changes. In other words, d ¼ 2 provides theﬁnest natural scale at which the sequences aggregate. (There is another non-random partition
corresponding to d ¼ 4, at which most sequences fall into three large clusters.) We interpret theexistence of a natural scale of non-random aggregation among HA amino acid sequences asevidence that the inﬂuenza A virus forms viral quasispecies.
The time series of inﬂuenza sequence clusters strengthens our interpretation of viral clusters as
selective units. Each of the 560 HA sequences is associated with a calendar year of viral isolationfrom an infected human. In Fig. 3 we plot the number of sequences in each cluster as a function ofisolation year. Even though we did not use any information about isolation year when clustering
the data, Fig. 3 shows that the resulting viral clusters are localized in time. There is no cluster thathas members spanning more than seven collection years. Instead, dominant clusters of viral sequences tend to replace one another every 2–5 years, in agreement with the time scale of dominantantigenic replacements [11].
Note that the cluster time series seen in Fig. 2 is not periodic within the time-span of two
decades. Once HA evolves away from a given region of sequence space, it does not later revisit
that region. This result, seen here in terms of cluster structure, is consistent with the one-trunk
phylogenetic reconstructions of HA [16].
Fig. 2 demonstrates that the general pattern of inﬂuenza A drift evolution is characterized by
non-random clusters of highly related genotypes that replace each other every 2–5 years. The

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

23

notion of an inﬂuenza ÔstrainÕ, which has hitherto been left unspeciﬁed, may arguably now be
identiﬁed with the observed dominant, non-random sequence clusters. At this scale of precision,
we see that inﬂuenza A virus exhibits signiﬁcant evolutionary change on a 2–5 year time scale, as
opposed to the annual periodicity of epidemics.
Our analysis has identiﬁed viral clusters on the basis of the Hamming distance between aminoacid sequences. The Hamming metric on sequences has been used as an eﬀective model ofantigenic distances or B-cells in general [23] and for inﬂuenza in particular [39]. Nevertheless,comparison of sequence data to direct immunological assays should be used to quantify theexact correlation between amino acid composition and antigenic properties [24]. In other words,
we have good reason to believe that all the sequences within a given cluster have similar patterns of cross-reactivity with outside sequences; but direct immunological conﬁrmation is stillrequired.
Once a natural scale of aggregation has been identiﬁed, the resulting empirical clusters can beinspected in greater detail. Previously, we have analyzed the spatio-temporal spread of sequences
within a cluster, and we have discussed the potential utility of HA clustering for the optimal
design of annual inﬂuenza vaccines [36]. We have also investigated patterns in the variation,cluster by cluster, of the dominant antibody-combining regions of HA, which may hold clues tothe mechanism of inﬂuenzaÕs antigenic plasticity. For the purposes of the present discussion,
however, it suﬃces to summarize that (i) we can detect non-random scales of viral sequence
aggregation; (ii) at the ﬁnest such scale, clusters are localized in time; and (iii) dominant clusters
replace one another every 2–5 years (as seen in Fig. 3).
The work described in this section provides a natural way to cluster strains in order to explore
the evolution and eﬀects of patterns of cross-reactivity in the inﬂuenza A virus. The dominant
clusters that have been identiﬁed can be associated with key epitopic regions, which determine
how the virus attaches to host cells. Given that there are a small number of such regions, it isreasonable to expect them to wax and wane in importance, demonstrating recurrent behavior overdecadal time scales.

4. Demographic stochasticity, seasonality and viral competition
As pointed out in Section 2, one of the factors that can reinforce the damped oscillations of the
S-I-R model, and allow persistent ﬂuctuations, is stochasticity. A natural form of stochasticity toconsider is demographic stochasticity, or the stochasticity that arises purely from considering the
population to be composed of discrete individuals: thus the instantaneous rate of new infections(from the determistic model) becomes an instantaneous probability that a single individual will be
infected.
Demographic stochasticity has surprisingly large eﬀects in S-I-R models, particularly when the
duration of immunity is long compared to the duration of infectiousness. Fig. 4 shows time series
for the number of infectious individuals in a deterministic S-I-R model, and in the corresponding
model with demographic stochasticity. While the deterministic model shows damped oscillations,oscillations in the stochastic model persist. The period of oscillations in the stochastic model ispredicted very well by the period in the deterministic model, while the long-term magnitude ofoscillations (assuming the disease has not gone stochastically extinct) is predicted very well by

24

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28
4000

Stochastic
Deterministic

3500
3000
2500
2000
1500
1000
500
0

0

200

400
600
800
Time (disease generations)1000

Fig. 4. The deterministic S-I-R model and one realization of the stochastic S-I-R model with the same parameters andinitial conditions. Even though the equilibrium number of infectious individuals is on the order of thousands, demographic stochasticity leads to large, persistent ﬂuctuations in the proportion infectious. The stochastic model is realized
by converting each of the arrows in Fig. 1 into an instantaneous rate of birth, death, infection or recovery.

simple diﬀusion approximations [33]. Simple approximations fail, however, to capture the meantime to extinction of stochastic S-I-R systems.
Since inﬂuenza shows annual cycles in prevalence, it is natural to ask about seasonal forcing,
and the extent to which exogenous forcing can amplify an endogenous tendency to oscillate. This
question has been explored in tremendous detail for several childhood diseases, most notably
measles [4,22]. There has been less work in this direction with inﬂuenza, both because there arefewer data available about inﬂuenza prevalence, and because the complication of drift evolution
makes the question much more diﬃcult.
The endogenous period of oscillation is critical to understanding how disease models are affected by seasonal forcing. For diseases with relatively high basic reproductive number R0 , andwhere the ratio of the duration of immunity to the duration of infection is high compared
pﬃﬃﬃﬃﬃﬃ to R0 ,the natural period of oscillation near the equilibrium is known to be approximately 2p ad , whered is the length of the cycle of infection (latency plus infectiousness) and a is the average age of ﬁrst
infection. The age at ﬁrst infection can be further approximated as a ¼ 1=ðR0 xÞ, where x is therate at which people exit the recovered class (through death or loss of immunity). Taking values of
6–10 days for d and 4–12 months for a yields a natural period of between 0.5 and 1 year. The
strong seasonal behavior observed in inﬂuenza is more consistent with values in the upper end of
this range, which show more regular seasonal oscillations, and also stronger magniﬁcation offorcing (see Fig. 5).
Stochastic models also shed light on competitive interactions, including those between drift
variants within a subtype of the inﬂuenza A virus, as discussed above, and between subtypes. At
least in some cases, stochastic ﬂuctuations can interact with other model factors to magnify
greatly the eﬀects of small competitive interactions.
Competing strains can interact if immune responses to one strain aﬀect the ability of others toinvade. Immune interaction can take the form of Ôdirect interferenceÕ, whereby a current infection

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

25

7000
6000
5000
4000
3000
2000
1000
0
0

1

2

3

4

5
6
Time (years)

7

8

9

10

Fig. 5. Time series for the stochastic S-I-R model with annual forcing and diﬀerent delays. In both simulations
N ¼ 500 000 and R0 varies sinusoidally between 9 and 11. The red curve has an infectious cycle of 0.02 year andduration of immunity of 4 years, giving a natural period of around 0.56 year; the blue curve has an infectious cycle of
0.025 year and duration of immunity of 8 years, giving a natural period of around 0.99 year.

3000
2500
2000
1500
1000
500
0
0

50

100

150

200

Time (disease generations)Fig. 6. Time series for a stochastic S-I-R with age structure, showing two strains that compete only by direct interference. Although the eﬀect of these strains on each other is small at the individual level, it is very large at the population level.

with one strain precludes infection by other strains; or of Ôcross-immunityÕ, whereby the long-term
immunity to a strain from which an individual has recovered reduces the ability of other strains to
infect that individual.
Fig. 6 shows a stochastic simulation with two strains that have no cross-immunity, but only
direct interference – so that an individual cannot be simultaneously infected by both strains. Since
each strain infects less than 1% of the population at any given time, simple deterministic
approximations would predict a negligible eﬀect of each strain on the other. When stochasticity is

26

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

combined with age structure (in the form of normally distributed lifetimes), however, the twostrains rapidly begin oscillating out of phase, greatly increasing the magnitude of oscillations, anddriving one strain extinct in a relatively short time, despite the fact that either strain could persist
for quite a long time on its own.
Thus, a very small competitive interaction at the individual level, which would have negligible
eﬀects in a deterministic approximation, has a very large eﬀect at the population level when
individuals are treated discretely. While more study is needed to determine the generality of thisphenomenon, it clearly has potentially important implications to understanding ﬂu ÔshiftsÕ and theriddle of why the old subtype is usually observed to go extinct when a new subtype comes in,despite the fact that the measured cross-immunity between diﬀerent subtypes is very small.

5. Discussion
For many diseases, the key to persistence and the bane of control is their continual and rapid
evolution. Diseases such as inﬂuenza A, which confers lifelong immunity to survivors of a particular strain, nonetheless reemerge year after year, causing tremendous morbidity and mortality,
through an ability to change form and escape host immune systems. Understanding this evolutionary change is fundamental to management, and to reduction in the diseaseÕs toll. Inﬂuenza A
is not unique in this – other diseases such as malaria, AIDS and hepatitis also can be caused by avariety of agents that are variants on a single theme.
There is, thus, a compelling need to extend the conventional modeling framework for infectious
diseases to deal with evolutionary changes in the infectious agent. This is not uncharted territory,
and has been an active area of research for 30 years or more [1,3,6,25,30]. Still, problems remain
regarding how best to represent antigenic diversity. For example, in the case of the hemagglutinin
gene of inﬂuenza A, thousands of diﬀerent drift variants can be identiﬁed from infected individuals in as short a time as a decade. To model the dynamics of such an assemblage, each ofwhich interacts uniquely with each other, is a task obviously beyond computation.
It is crucial, therefore, to ﬁnd ways to partition variants into a small number of clusters, so thatthe variants within each cluster can be treated as functionally equivalent, interacting strongly
among themselves and more weakly with strains in other clusters. Such techniques have long beenused informally for a variety of diseases, for example myxomatosis [15], and indeed form the basis
for vaccine choice [11,12] for inﬂuenza A. Still, as molecular techniques become more powerful,and computational algorithms more available, it is worthwhile to develop formal approaches tothe problem of clustering strains, and to investigating the relationship between genetic sequencesof hemagglutinin and their immunological properties.
We have shown that indeed it is possible to separate the near-inﬁnity of strains into a small
number of clusters, associated with the key epitopic regions of hemagglutinin (see also [36]). This
work then provides a foundation for multitype modeling of infectious dynamics, and for guiding
vaccine choice.
Ultimately, the diﬃculties faced in traditional epidemiological modeling relate to the necessityof lumping individuals into categories. With even 3 virus types, there are 20 categories of individuals; with 4, the number increases to 48. An alternative way to proceed is to develop individual-based models, in which unique life histories can be followed, and stochastic inﬂuences

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28

27

easily incorporated. We show that such models can sustain oscillatory behavior when theirdeterministic equivalents cannot. Most importantly, these individual-based models illuminate therole of cross-reactivity among clusters, or among subtypes. Remarkably, it appears that only
minimal interaction among strains or clusters can translate into strong inﬂuences on oscillatory
behavior, and on patterns of cocirculation, at the population level. These results beg for deeper
investigation, but may shed light on the puzzle of why subtypes, which interact only weakly,
typically replace one another in temporal sequences. Similarly, we ﬁnd that only minimal amountsof seasonal forcing can translate into large oscillations in disease prevalence.

Acknowledgements
The authors are grateful for valuable discussions with Viggo Andreasen, David Earn, Ellis
McKenzie, Tom Reichert, Lone Simonsen, and Martin Weigert, and for the comments of two
anonymous reviewers. We thank Los Alamos National Laboratory and everybody who has
contributed to the Inﬂuenza Sequence Database, including those who have contributed unpublished sequences. This work was supported by a grant from the National Institutes of Health,number 1-R01-GM60729-01 to S.A.L.; J.B.P. acknowledges support from the the National Science Foundation and the Burroughs Wellcome Fund award number 1001782.

References[1] R.M. Anderson, R.M. May, Population biology of infectious diseases 1, Nature 280 (1979) 361.
[2] R.M. Anderson, R.M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford Science Publications,
Oxford, 1991.
[3] V. Andreasen, J. Lin, S.A. Levin, The dynamics of cocirculating inﬂuenza strains conferring partial crossimmunity, J. Math. Biol. 35 (1997) 825.
[4] B. Bolker, B.T. Grenfell, Space, persistence and dynamics of measles epidemics, Philos. Trans. R. Soc. Lond. Ser.
B––Biol. Sci. 348 (1995) 309.
[5] B.M. Bolker, B.T. Grenfell, Impact of vaccination on the spatial correlation and persistence of measles dynamics,
Proc. Nat. Acad. Sci. USA 93 (1996) 12648.
[6] H.J. Bremermann, J. Pickering, A game-theoretical model of parasite virulence, J. Theor. Biol. 100 (1983) 411.
[7] R.M. Bush, Opinion: predicting adaptive evolution, Nat. Rev. Genet. 2 (2001) 387.
[8] R.M. Bush, W.M. Fitch, C.A. Bender, N.J. Cox, Positive selection on the H3 hemagglutinin gene of humaninﬂuenza virus A, Mol. Biol. Evol. 16 (1999) 1457.
[9] R.M. Bush, C.B. Smith, N.J. Cox, W.M. Fitch, Eﬀects of passage history and sampling bias on phylogenetic
reconstruction of human inﬂuenza A evolution, Proc. Nat. Acad. Sci. USA 97 (2000) 6974.
[10] C. Castillo-Chavez, H.W. Hethcote, V. Andreasen, S.A. Levin, W.M. Liu, Epidemiological models with age
structure, proportionate mixing, and cross-immunity, J. Math. Biol. 27 (1989) 233.
[11] N.J. Cox, C.A. Bender, The molecular epidemiology of inﬂuenza viruses, Semin. Virol. 6 (1995) 359.
[12] N.J. Cox, K. Subbarao, Inﬂuenza, Lancet 354 (1999) 1277.
[13] E. Domingo, J.J. Holland, C.K. Biebricher, Quasispecies and RNA Virus Evolution: Principles and Consequences,
R.G. Landes Co., Austin, 2002.
[14] E. Domingo, L. Mededez-Arias, M. Quinones-Mateu, A. Holguin, M. Gutierrez-Rivas, M. Martinez, J. Quer,
I. Novella, J. Holland, Viral quasispecies and the problem of vaccine-escape and drug-resistant mutants, Prog.
Drug Res. 48 (1997) 99.

28

S.A. Levin et al. / Mathematical Biosciences 188 (2004) 17–28[15] F. Fenner, I.D. Marshall, Occurrence of attenuated strains of myxoma virus in Europe, Nature 176 (1955) 782.
[16] W.M. Fitch, R.M. Bush, C.A. Bender, N.J. Cox, Long term trends in the evolution of H(3) HA1 human inﬂuenza
type A, Proc. Nat. Acad. Sci. USA 94 (1997) 7712.
[17] W.M. Fitch, R.M. Bush, C.A. Bender, K. Subbarao, N.J. Cox, Predicting the evolution of human inﬂuenza A,
J. Hered. 91 (2000) 183.
[18] J.R. Gog, B.T. Grenfell, Dynamics and selection of many-strain pathogens, Proc. Nat. Acad. Sci. USA 99 (2002)
17209.
[19] J.R. Gog, J. Swinton, A status-based approach to multiple strain dynamics, J. Math. Biol. 44 (2002) 169.
[20] B.T. Grenfell, O.N. Bjornstad, J. Kappey, Travelling waves and spatial hierarchies in measles epidemics, Nature
414 (2001) 716.
[21] Y. Ina, T. Gojobori, Statistical analysis of nucleotide-sequences of the hemagglutinin gene of human inﬂuenza A
viruses, Proc. Nat. Acad. Sci. USA 91 (1994) 8388.
[22] M.J. Keeling, B.T. Grenfell, Understanding the persistence of measles: Reconciling theory simulation andobservation, Proc. R. Soc. Lond. Ser. B––Biol. Sci. 269 (2002) 335.
[23] V. Kohler, R. Puzone, P.E. Seiden, F. Celada, A systematic approach to vaccine complexity using an automaton
model of the cellular and humoral immune system, Vaccine 19 (2000) 862.
[24] A. Lapedes, R. Farber, The geometry of shape space: application to inﬂuenza, J. Theor. Biol. 212 (2001) 57.
[25] S. Levin, D. Pimentel, Selection of intermediate rates of increase in parasite-host systems, Am. Natural. 117 (1981)
308.
[26] J. Lin, V. Andreasen, S.A. Levin, Dynamics of inﬂuenza A drift: the linear three-strain model, Math. Biosci. 162(1999) 33.
[27] W.M. Liu, H.W. Hethcote, S.A. Levin, Dynamic behavior of epidemiologic models with nonlinear incidence rates,
J. Math. Biol. 25 (1987) 359.
[28] Los Alamos National Labs. Inﬂuenza sequence database. Available from <http://www.ﬂu.lanl.gov/>, World wide
web site.
[29] C. Macken, H. Lu, J. Goodman, L. Boykin, The value of a database in surveillance and vaccine selection, in:
A.D.M.E. Osterhaus, N. Cox, A.W. Hampson (Eds.), Options for the Control of Inﬂuenza IV, Elsevier,
Amsterdam, 2001.
[30] R.M. May, R.M. Anderson, Population biology of infectious diseases 2, Nature 280 (1979) 455.
[31] R.M. May, R.M. Anderson, Periodicity in epidemiological models, in: S.A. Levin, T.G. Hallam, L.J. Gross (Eds.),
Applied Mathematical Ecology, Springer, Berlin, 1989, p. 193.
[32] R. Meester, R. Roy, Continuum Percolation, Cambridge University, Cambridge, 1996.
[33] I. Nasell, On the time to extinction in recurrent epidemics, J. R. Stat. Soc. Ser. B––Stat. Methodol. 61 (1999) 309.
[34] M.A. Nowak, R. May, Virus Dynamics: The Mathematical Foundations of Immunology and Virology, Oxford
University, Oxford, 2000.
[35] J.B. Plotkin, J.M. Chave, P.S. Ashton, Cluster analysis of spatial patterns in Malaysian tree species, Am. Natural.
160 (2002) 629.
[36] J.B. Plotkin, J. Dushoﬀ, S.A. Levin, Hemagglutinin sequence clusters and the antigenic evolution of inﬂuenza A
virus, Proc. Nat. Acad. Sci. USA 99 (2002) 6263.
[37] L. Simonsen, M.J. Clarke, L.B. Schonberger, N.H. Arden, N.J. Cox, K. Fukuda, Pandemic versus epidemic
inﬂuenza mortality: a pattern of changing age distribution, J. Infect. Dis. 178 (1998) 53.
[38] L. Simonsen, K. Fukuda, L.B. Schonberger, N.J. Cox, The impact of inﬂuenza epidemics on hospitalizations,
J. Infect. Dis. 181 (2000) 831.
[39] D.J. Smith, S. Forrest, D.H. Ackley, A.S. Perelson, Variable eﬃcacy of repeated annual inﬂuenza vaccination,
Proc. Nat. Acad. Sci. USA 96 (24) (1999) 14001.
[40] J.K. Taubenberger, A.H. Reid, T.G. Fanning, The 1918 inﬂuenza virus: a killer comes into view, Virology 274(2000) 241.