F. R. Adler and J. M. C. Pearce-Duvet and M. D. Dearing.

17 Pages·2008·0.45 MB·English

BulletinofMathematicalBiology(2008)70:236–252
DOI10.1007/s11538-007-9251-8
ORIGINAL ARTICLE
How Host Population Dynamics Translate into Time-Lagged
Prevalence: An Investigation of Sin Nombre Virus
in Deer Mice
FrederickR.Adlera,∗,JessicaM.C.Pearce-Duvetb,M.DeniseDearingb
a
DepartmentofMathematicsandDepartmentofBiology,155South1400East,
UniversityofUtah,SaltLakeCity,UT84112,USA
b
DepartmentofBiology,257South1400East,UniversityofUtah,SaltLakeCity,
UT84112,USA
Received:5January2007/Accepted:22June2007/Publishedonline:16August2007
©SocietyforMathematicalBiology2007
Abstract HumancasesofhantaviruspulmonarysyndromecausedbySinNombrevirus
aretheendpointofcomplexecologicalcascadefromweatherconditions,populationdy-
namics of deer mice, to prevalence of SNV in deer mice. Using population trajectories
from the literature and mathematical modeling, we analyze the time lag between deer
mousepopulationpeaksandpeaksinSNVantibodyprevalenceindeermice.Becausethe
virus is not transmitted vertically, rapid population growth can lead initially to reduced
prevalence,buttheresultinghigherpopulationsizemaylaterincreasecontactratesand
generateincreasedprevalence.Incorporatingthesefactors,thepredictedtimelagranges
from0to18months,andtakesonlargervalueswhenhostpopulationsizevarieswitha
longerperiodorhigheramplitude,whenmeanprevalenceislowandwhentransmission
is frequency-dependent. Population size variation due to variation in birth rates rather
thandeathratesalsoincreasesthelag.Predictingfuturehumanoutbreaksofhantavirus
pulmonarysyndromemayrequiretakingtheseeffectsintoaccount.
Keywords Hantavirus·SinNombrevirus·Timelags
1. Introduction
The 1993–1994 outbreak of hantavirus pulmonary syndrome (HPS) in humans in the
American southwest was caused by Sin Nombre virus (SNV), and was quickly cor-
related with the high rainfall associated with the large El Nino (ENSO) event of the
previous year that led to the increased populations of deer mice (Peromyscus manic-
ulatus), the primary host of the virus (Schmaljohn and Hjelle, 1997; Engelthaler et
al., 1999). However, further study has shown that the link between weather and HPS
∗
Correspondingauthor.
E-mailaddress:[email protected](FrederickR.Adler).
TimeLagsbetweenHostPopulationDynamicsandDisease 237
cases is quite complex (Glass et al., 2000), involving at best a noisy trophic cascade
from weather, plants, deer mouse population size, and deer mouse infection preva-
lence before humans are finally affected (Engelthaler et al., 1999; Mills et al., 1999;
Yates et al., 2002; Kolivras and Comrie, 2004). The magnitude and timing of the re-
sponse at each step depends on the detailed behavior of individual species: the many
species of plants, the feeding and reproduction of deer mice, the transmission dy-
namics of the virus between deer mice, and the transmission to humans. The ecolog-
ical steps, from rain to rodents, have proven to be quite unpredictable (Brown and
Ernest, 2002). Indeed, there is debate about whether the 1997 ENSO event led to
an increase in HPS cases beginning in the spring of 1998 (Hjelle and Glass, 2000;
KolivrasandComrie,2004).
Although frequently fatal to humans, the effects of SNV on deer mice are nonexis-
tent or weak (Helle and Yates, 2001), although there is evidence of higher mortality in
subadults(Douglassetal.,2001).Deer miceapparently remain seropositivethroughout
life,andtheinfectioncanbemodeledasifthereisnorecovery(Millsetal.,1999;Kuenzi
etal.,2005),althoughantibodylossmaybepossible(Booneetal.,1998).Transmission
isprimarilythroughbiting(Bottenetal.,2002),althoughthereisevidenceoforal/fecal
transmission in burrows (Mills et al., 1999), with no vertical transmission from parents
tooffspring.Prevalenceofinfectionisgenerallyhigherinmalesandinolderindividuals,
and large older males are particularly likely to carry the infection (Boone et al., 1998;
Abbottetal.,1999).
Indiseaseswithoutverticaltransmission,populationdynamicscanhavetwoconflict-
ingeffectsonprevalence(Millsetal.,1999).Duringperiodsofrapidpopulationgrowth,
the lack of vertical transmission leads to dilution and potentially reduced prevalence.
However,theresultinghigherpopulationsizeislikelytoincreasecontactratesbetween
deermicethatcouldeventuallyleadtoanincreaseinprevalence.
Theseconflictingeffectshaveappearedintheliteratureasinconsistentresultsregard-
ingthelinkbetweenpopulationsizeandSNVseroprevalence(Table1).Onlyafewstud-
ieshavefoundapositiverelationshipbetweencurrentpopulationsizeandseroprevalence,
andothershavefoundapositiverelationshiponlyinthepresenceofalag(Davisetal.,
2005).Inthesystemunderconsiderationinthispaper,severalstudieshavefoundnolinear
relationshipbetweencurrentpopulationsizeandprevalence(GrahamandChomel,1997;
Millsetal.,1997;Booneetal.,1998),andourreanalysisofdatapresentedintwootherpa-
persshowsthesameresult(Calisheretal.,1999;Rootetal.,1999).Thereisevidencefora
thresholdeffectinthatsufficientlysmallpopulationsshowedaprevalenceofzero(Boone
etal.,1998).NegativecorrelationswerefoundintwooutoftentimeseriesofSNVsero-
prevalenceversuspopulationsizeinMontana(Douglassetal.,2001).Usingmultivariate
analysis,onestudyfoundapositiveeffectofdeermousedensityonprevalence(Biggset
al.,2000).
Studiesofotherhantaviruseshavefoundsimilarlyinconsistentresults.Intheclosely
relatedLimestoneCanyonvirus(Sanchezetal.,2001),higherprevalencewasobserved
during periods of low population density (Abbott et al., 1999), consistent with a nega-
tivecorrelation.Arecentstudyfoundanegativecorrelationbetweendensityandpreva-
lence in El Moro Canyon virus in western harvest mice (Calisher et al., 2005). In bank
voles,Puumalavirusseroprevalencewaswellpredicted inpartbythephaseofpopula-
tion,beinghighestafterthepeak,implyingbothdensity-dependenceanddelayeddensity-
dependence(Olssonetal.,2002,2003).
238 Adleretal.
Table1 Correlationbetweencurrentdensityandprevalence
Meanprevalence(%) Relationshipa Reference
9.5 0 (Calisheretal.,1999)
11.0 0 (Millsetal.,1997)
11.9 − (Calisheretal.,2005)
14.3 0 (Rootetal.,1999)
15.8 − (Douglassetal.,2001)
16.5 + (Biggsetal.,2000)
17.0 0 (GrahamandChomel,1997)
17.8 − (Abbottetal.,1999)
18.5 0 (Booneetal.,1998)
a+ indicates a significant positive relationship, − a significant negative
relationship,and0thelackofasignificantrelationship
In plague, another disease of rodents, one study found that a model with a two year
delay performed well in predicting presence or absence of disease as a function of the
proportionofburrowsoccupiedbygreatgerbils(Davisetal.,2004).Frequencyofhuman
plague cases increased with a one or two year lag as a function of winter or summer
precipitation(Enscoreetal.,2002),althoughthiscouldbeduetoadirectnumberseffect
ratherthanincreasesinseroprevalence.
The interplay between dilution due to births and increased prevalence due to higher
contact rates depends sensitively on how contact rates vary with host density. Models
range from a constant contact or frequency-dependent model, where the contact rate is
independentofpopulationsizeandtheforceofinfectionisproportionaltothefrequency
ofinfectedhosts,toamassactionordensity-dependentmodel,wherethecontactrateis
proportional to population size and the force of infection is proportional to the number
of infected hosts (McCallum et al., 2001). Intermediate models, where the contact rate
saturatesashostdensitybecomeslarge,canresultfrombehavioralchangesatlowden-
sity (Ramsey et al., 2002) or from host heterogeneity (Barlow, 2000). Such models can
produce substantially different predictions about disease dynamics and persistence, and
may be distinguishable with field data (Caley and Ramsey, 2001). The contact rate de-
pends on the territorial behavior of hosts (Wolff, 1985). As density increases from low
levels, deer mice decrease home range size, thus potentially leaving contact rates un-
changed.Asdensityincreasesfurther,deermiceincreaseaggressionandperhapseffective
contactrateswithoutfurtherdecreaseinhomerangesize.
Mathematical models of the temporal dynamics of hantaviruses have yet to explic-
itly consider lags. The models have focused instead on spatial spread and the role of
refugia for the virus when population sizes are small (Abramson and Kenkre, 2002).
Furtheranalysisshowedthatseasonalitycanpromotediseasepersistenceandcauseout-
breaks(Bucetaetal.,2004).Inparticular,themaximuminfectedpopulationoccursnear
theendoftheharshseason,andthussubstantiallydelayedfromtheperiodofoptimalre-
production.Arecentpapershowsthathumanoutbreakscanbeinducedbypopulationsize
increasesinrodentsthatquicklygeneratemanyhighlyinfectiousnewinfections(Sauvage
etal.,2007).
Thispaperidentifiesasimpleepidemiologicalmodelthatcancapturethemajortrends
inprevalence,andusesthatmodeltostudytherelationshipbetweenpopulationsizeand
TimeLagsbetweenHostPopulationDynamicsandDisease 239
prevalence.Thekeyparametersarebirth,death,andtransmissionrates,whichwehypoth-
esizecreatecharacteristictimelagsbetweenpopulationsizeandprevalence.Webeginby
fittingpublisheddata(Yatesetal.,2002)withasimplemodel,andillustratethesensitiv-
ityoftheresultstothemodelparameters,particularlythosedescribingtransmission.We
thenusethissimplemodeltoshowthatthelaghasacomplexrelationshipwithpopula-
tiondynamicswhichdependsonthedetailsofthetransmissionprocessinadditiontohost
demographicparameters.
2. Modelandparameterestimates
Our model is an SI model of susceptible and infected individuals (Anderson and May,
1992)
dS
=−cg(N)pS+bN−δS,
dt
dI
=cg(N)pS−δI,
dt
whereS,I andN representthedensityofsusceptible,infected,andtotalhosts,respec-
tively,andp= I istheprevalence.Theparameterc describestheeffectivecontactrate
N
per day, and b and δ give the birth and death rates, respectively, and can be functions
of time. Table 2 gives a complete description of the parameters in the model and their
baselinevalues.
The function g(N) describes the dependence of contact rate on density. We chose a
variantofasymptotictransmission(McCallumetal.,2001)withtheform
(cid:2) (cid:3)(cid:2) (cid:3)
N +K N
g(N)= d . (1)
N N+K
d
HereK givesthedensitywherecontactis50%ofmaximum,andN isascalingfactor
d
equal to the density where g(N)=1. If K =0, g(N)=1 for all N, corresponding to
amodelwithconstantcontactratesandfrequency-dependenttransmission.AsK→∞,
Table2 Variablesandparametersinthemodel
Symbol Description Baselinevalue
N Densityofdeermice
I Densityofinfecteddeermice
S Densityofsusceptibledeermice
p Prevalenceofinfection
δ Mortalityrateperday 0.01
b Birthrateperday(afternatalmortality) 0.02
c Effectivecontactrateperday 0.0105withvariationinbirths
0.022withvariationindeaths
(cid:4) (cid:5)(cid:4) (cid:5)
g(N) Dependenceofcontactrateondensity g(N)= NdN+dK NN+K
K Densitywherecontactis50%ofmaximum 50mice/ha
Nd Densitywhereg(N)=1 7mice/ha
240 Adleretal.
g(N)→N/N ,correspondingtoamassactionmodelanddensity-dependenttransmis-
d
sion.
TheequationscanberewrittenintermsoftheprevalencepandtotaldensityN as
dp
=cg(N)p(1−p)−bp, (2)
dt
dN
=bN−δN. (3)
dt
Mortality does not differ between susceptible and infected deer mice and thus plays no
directroleintheequationforprevalence.Ahigherbirthratebreducesprevalencebecause
thereisnoverticaltransmission.Ifcontactratesincreasewithincreasingpopulationsize,
ahigherbirthratewillleadtoincreasedpopulationsizeandadelayedincreaseinpreva-
lence.Thediseasewillpersistinanisolatedpopulationwithconstantcoefficientsifthe
averageeffectivecontactratecg(N)exceedstheaverageofthebirthrateb.
Ourparameterestimatesforthebaselinebirthanddeathratescomefromstudiesthat
deducedeathratesfrompopulationchangeduringthenonbreedingwinterseason(roughly
δ ≈0.008) and birth rates from the growth during the breeding season (r =b−δ ≈
0.007)(PetticrewandSadleir,1974;Fairbairn,1977).ThebaselinevalueofKwaschosen
basedontheassumptionthatmousebehaviorchangesatdensitiesofover30deermice
perhectare(Wolff,1985).ThescalingfactorN isthemeandensityfoundinYatesetal.
d
(2002).Thetransmissionratescwerefoundbymatchingtheobservedmeanprevalence
in Yates et al. (2002). With a large value of K, our model approximates “pseudo mass
actionmodel”inastudyofcowpoxinbankvoles(Begonetal.,1998),whichfound,in
theunitsofthecurrentwork,cintherangefrom0.006to0.008,slightlylowerthanour
estimates.
3. Comparingthemodelwithdata
WeextracteddataonthedensityofdeermiceandofseropositivedeermicefromYates
etal.(2002),andshowtherawdataandsmootheddataoveran87monthperiod(created
with the supsmu function in R (R Development Core Team, 2005)) in Fig. 1a. We use
smoothed data in the following analyses to average over factors creating noise in the
originaldataandtointerpolatevaluesformissingdatapoints.
We examined the correlations between prevalence and lagged population size. We
found the strongest correlation with population size 15 months earlier, although there
wassignificantcorrelationwithalagof12monthsasfoundbyYatesandothers(Fig.1b).
Usingthetimeseriesforthedensityofallmice,wecomputedthegrowthrater=b−δ
asafunctionoftimethatexactlymatchesthesmootheddensitytrajectory.Wecreatedtwo
versionsofourmodel:
• Variablebirthrates:Setthedeathrateδ(t)tobeconstantandsolveforthebirthrate
b(t)asafunctionoftime.
• Variabledeathrates:Setthebirthrateb(t)tobeconstantandsolveforthedeathrate
δ(t)asafunctionoftime.
TimeLagsbetweenHostPopulationDynamicsandDisease 241
Fig.1 (a)ComparisonofrawdataandsmootheddatafromYatesetal.(2002)showingthedensityof
deermiceandofseropositivedeermiceoveran87monthperiod.(b)Plotofprevalenceagainstdensity
withtheoptimallagof15months,chosenasthelagthatmaximizesr2.
Thisfitdoesnotuseanyinformationaboutseroprevalence.Withourbaselineparameter
values, the model predicts the dynamics of the infection reasonably well (Fig. 2). The
modelwithvaryingdeathratesrequiresalargervalueofthecontactrate c becauseour
baselinebirthrateishigherthanthebaselinedeathrate(Fig.2b).
242 Adleretal.
Fig.2 Modeleddensityofseropositivemicewith(a)variablebirthrates,and(b)variabledeathrates.The
modeleddensityofallmiceissettoexactlymatchthesmootheddata.ParametervaluesasinTable2.
The two cases make slightly different predictions about the correlations with lagged
populationsizeandgrowthrate.Themodelwithvariablebirthratesmatchestheobserved
lag somewhat better (Fig. 3) but has lower r2 values. It generates a significant positive
correlationwhenthelagisatleast6months,andasignificantnegativecorrelationwhen
TimeLagsbetweenHostPopulationDynamicsandDisease 243
Fig.3 Thecoefficientofdeterminationr2ofprevalenceasafunctionoflaggeddensityforarangeoflags
fortheactualdata(opencircles),andthetwosimulationsshowninFig.2(b’swithvariablebirthrates,d’s
withvariabledeathrates).
thelagis0or1month.Themodelwithvariabledeathratesgeneratesasignificantpositive
correlationwhenthelagisatleast2months.
The model results are extremely sensitive to the choice of the contact rate c. Small
increasesordecreasesleadnotonlytolargechangesinthenumberofinfectedindividuals,
butalsotochangesinthelagbetweenpeaksofpopulationdensityandofseropositivity,
withhighercontactratesgeneratingasubstantiallyshorterlag(Fig.4).
4. Theory:lagsinducedbyoscillatingpopulationsize
What aspects of the data produce the observed lags? In this section, we use linear per-
turbationtheorytocalculatethelagsinducedbysinusoidalpopulationoscillationswith
different periods, and then use Fourier analysis to decompose the observed population
sizedataintoasumofsinusoidaloscillationstoidentifykeyaspectsofthedata.
4.1. Linearperturbationtheory
AsdescribedindetailintheAppendix,weassumethatthepopulationsizeN(t)variesin
asmallamplitudesinusoidaloscillationaroundameanvalueofN ,or
0
N(t)=N +(cid:3)αcos(2πωt).
0
Here,(cid:3) isasmallparameter,α scalesthemagnitudeoftheoscillation,andω isthefre-
quencyoftheoscillation.Usingtheequationforpopulationsize(Eq.(3)),wecancompute
thesmallamplitudesinusoidaloscillationdescribingthebirthordeathrate,dependingon
thecaseunderconsideration.
244 Adleretal.
Fig.4 Effectsofchangingthecontactrateconthedynamicsofseroprevalencewith(a)variablebirth
rates,and(b)variabledeathrates.ParametervaluesasinTable2.
Substitutingtheseoscillationsintothebasicequationforprevalence(Eq.(2)),wecan
write
(cid:4) (cid:5)
p(t)=p∗+(cid:3)Zcos 2πω(t−φ) ,
TimeLagsbetweenHostPopulationDynamicsandDisease 245
Fig.5 (a)Comparisonofcomputedlagsfromthelineartheorywiththestatisticallyestimatedlagsfrom
simulations.Solidsquaresshowsimulationswithα=1,alowamplitudeoscillation,andopensquares
simulationswithα=4,anamplitudecomparabletothatintherealdata.Thestatisticallyestimatedlags
arethosegivingthehighest r2 value.Parameterswerevariedseparatelyaroundthebaselinevaluesin
Table2,withcrunningfrom0.01to0.03inthecasewithvariablebirthratesandfrom0.02to0.06in
thecasewithvariablebirthratesand K runningfrom0to1024.Theunderlyingperiodis87months
throughout.Allsimulationswererunfor240months.(b)Theeffectofmeanprevalenceontheobserved
lag,usingtheparametervaluesinTable2withα=4andanunderlyingperiodof87monthsthroughout.
Opencirclesshowthecasewithvariablebirthrates,andsolidcirclesshowthecasewithvariabledeath
rates.(c)Theeffectofperiodontheobservedlag,usingperiodsinmultiplesof12monthsfrom12to
84.OtherparametersasinTable2.(d)TheeffectofK onlagsusingtheparametervaluesinTable2
andanunderlyingoscillationwithamplitude4andperiod87months,butadjustingctomatchthemean
prevalenceobservedinthedata.
wherep∗isthemeanprevalence,Zisthescaledamplitude,andφisthephaseshiftorlag
inducedbythedynamics.Wecanusemethodsfromperturbationtheorytosolveexplicitly
forZandφintermsoftheparametervalues.
Figure5ashowsthattheexactlagfromthelineartheorycompareswellwiththestatis-
ticalfitstothesimulation,andthathigheramplitudeoscillations(ofasimilarmagnitude
to those observed) produce longer lags than small amplitude oscillations. Much of the
variationinlagsforafixedperiodofthepopulationsizeoscillationisgeneratedbytwo
factors(Fig.5b).Highermeanprevalenceproducessubstantiallyshorterlagsasdoesat-
tributingthevariationinpopulationsizetovariationinthedeathrate.Theobservedlag
does not vary linearly with the period of the population size oscillation, but eventually
saturatestoafixedvaluethatdependsontheotherparameters(Fig.5c).Ifwecontrolfor
themeanprevalencebyadjustingc,largervaluesofK leadtoshorterlagswithsubstan-