1
Supporting Information
2
3
Timing and severity of immunizing diseases in rabbits is
4
controlled by seasonal matching of host and pathogen dynamics
5
Konstans Wells1,*, Barry W. Brook1, Robert C. Lacy2, Greg J. Mutze3, David E. Peacock3, Ron G.
6
Sinclair3, Nina Schwensow1, Phillip Cassey1, Robert B. O’Hara4, Damien A. Fordham1
7
8
1 The Environment Institute and School of Earth and Environmental Sciences, The University of
9
Adelaide, SA 5005, Australia
10
2 Chicago Zoological Society, Brookfield, Illinois, United States of America
11
3 Department of Primary Industries and Regions, Biosecurity SA, Adelaide, SA 5001, Australia
12
4 Biodiversity and Climate Research Centre (BIK-F), Senckenberganlage 25, 60325 Frankfurt am
13
Main, Germany
14
*
Author for correspondence: konstans.wells@adelaide.edu.au
15
16
17
Demographic and epidemiological model
18
An individual-based model of seasonal rabbit demography and two co-circulating diseases
19
in a meta-model framework – overview and specification
20
Overview
21
To build our individual-based model of rabbit demography and coupled epidemiology of rabbit
22
haemorrhagic disease (RHD) and myxomatosis we used the freely available software packages Vortex
23
10.0 and Outbreak 2.1, which were linked through the software Meta-Model Manager 1.0 [1, 2]. The
24
software and user manuals can be downloaded at http://www.Vortex10.org. In brief, a demographic
25
model in Vortex simulates the fate and reproduction of individuals over discrete time steps with
26
various deterministic and stochastic forces [3, 4]. The software was initially conceptualized to model
27
population viability over years as discrete time steps [3]. However, ‘days per year’ can be specified
28
for modelling shorter time steps, allowing parameters such as demographic rates to be varied over
29
shorter time steps. Linking epidemiological models in Outbreak through Meta Model Manager,
Supporting Information
30
allows the disease state and individual fate (i.e. death through disease) of organisms to be modified
31
for discrete time steps encapsulated within the time steps defined in Vortex (typically, disease
32
transmission dynamics are modelled with daily intervals). The different models are linked through
33
defined state variables describing age, disease state, and other characteristics of individuals, and
34
Meta-Model Manager facilitates the matching of events in time and space.
35
36
We developed a seasonal model in which we varied reproductive efforts over weeks. We
37
defined ‘years’ in Vortex as 7-day time steps and allowed reproductive efforts to vary in different time
38
steps. We achieved this by modelling seasonal parameters as functions. For example, if reproduction
39
is assumed to vary among calendar weeks in Vortex, the parameter ‘PercentBreed’ can be specified as
“= ((Y%52=0)*X1 *(K-N)/K) + ((Y%52=1)*X2 *(K-N)/K) +((Y%52=2)*X3, …, +
40
41
((Y%52=51)*X52 *(K-N)/K)”,
42
with Xw being the weekly parameter value in week w and (K-N)/K) accounts for density-dependent
43
regulation of reproduction.
44
45
We modelled the increasing susceptibility of infants (1-90 days old) to rabbit haemorrhagic
46
disease virus (RHDV)[5] and the decreasing susceptibility to myxoma viruses (MYXV)[6] with a
47
logistic model as a function of age. For this, we assumed that the recovery rate from RHD after the
48
insusceptibility period declined from a maximum of 100 % (i.e. a fixed intercept μJuv,RHD = 6 in the
49
logit-link logistic regression model) to the recovery rates to those of adults (V). The regression
50
coefficient (Juv-RHD, ‘recovery adjustment factor’ for juveniles) for infant age, which measure the rate
51
of change in infant recovery with increasing age was sampled during simulations (Figure A.1) and
52
translates into the dynamic change in infant recovery rate with Juv,RHD:
logit(Juv,RHD) = μRHD + Juv,RHD Ageadjusted .
53
54
Here, the variable ‘Ageadjusted’ accounting for the insusceptibility period, i.e. the onset for an
55
increasing RHD susceptibility is at an assumed age of 22 days for RHD (insusceptibility period for 21
56
days, [7]).
57
In Outbreak, the dynamic model for recovery rates can be expressed as a function for ‘probRecovery’
58
as
59
60
61
“= ((Age<90)* RHD + (1 - RHD)* Juv,RHD) + ((Age>90) *RHD”
whereby
Juv,RHD = exp(μJuv,RHD + Juv,RHD *A - 21)/(1+exp(μ Juv,RHD + Juv,RHD *A - 21).
62
63
An inverse relationship was assumed for myxomatosis (the negative value of the logistic
64
model) to account for decreasing susceptibility of infants with increasing age. Here, we sampled the
65
regression intercept (μJuv,Myxo) rather than the slope as an unknown parameter; the model assumes
2
Supporting Information
66
maximum susceptibility to myxomatosis for juveniles after the insusceptibility period, which is
67
determined by μJuv,Myxo and thus assume to be the critical factor for disease dynamics (see Figure A.1).
68
69
70
Figure A.1. Illustration of the dynamic modelling of decreasing recovery rates of infants for RHD and
71
increasing recovery rates for myxomatosis. The rate of change in infant recovery towards the values of adult
72
recovery rats is models with a logistic function after the assumed time of insusceptibility (grey bars) to the age
73
of 90 days, when recovery rates are assumed to approach those of adults (V). For this, regression slopes
74
(RHD,Juv) are sampled for RHD and regression intercepts (μJuv,Myxo) are sampled for myxomatosis disease
75
dynamics. Dashed lines indicate of possible behaviour of the dynamic model across the sampled parameter
76
values (Juv,RHD, μJuv,Myxo).
77
78
79
We modelled waning maternal antibodies against RHDV and MYXV and the resulting
80
decrease in recovery rates from disease of infants (1-90 days of age) as an additional component for
81
‘probRecovery’, assuming a linear and constant decline.
82
For the implementation in Outbreak, we assumed a state variable ImAB of value ImAB = 100 if infants
83
received antibodies from resistant does and ImAB = 0 otherwise. A is a function variable of individual
84
age (see Vortex/Outbreak software manual). Since the effect of maternal antibodies is additive to the
85
changes in infant recovery rates due to increasing/decreasing susceptibility the model for
86
‘probRecovery’ changes to
87
“= ((Age<90)* RHD + (1 - RHD)* Juv,RHD + (1 - RHD,Juv)* 0.01* (MAX(0;(ImAB - A)))) +
88
((Age>90) *RHD)”,
89
whereby
90
RHD,Juv =RHD + (1 - RHDt)* Juv,RHD.
3
Supporting Information
91
92
Note that the actual effect of maternal antibodies in each simulation are also determined by
93
the simulated time period of full protection (tMRHD), during which no infections are possible. For the
94
sake of model parsimony and lack of detailed information, we assumed transmission probability βRHD
95
to be constant and not reduced by maternal antibodies after the time period of full protection.
96
Age of individuals in Vortex are counted as time steps (i.e. weeks in our model), while the age
97
of individuals in Outbreak are defined in ‘days’, so that the time periods of disease states such as
98
maternal antibodies can be flexibly accounted for.
99
Meta-Model Manager communicated all relevant information between the different software
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
components to fully take the changes in individual state variables over time into account.
4
Supporting Information
117
118
Figure A.2. Flow diagram for an individual-based epidemiological model, used for examining the effect of
119
seasonal timing of host reproduction and virus exposure (RHD and myxomatosis) on disease dynamics and
120
population persistence. Panel (a) shows the possible individual progression through different disease states (P:
121
Pre-/Insusceptible, S: Susceptible, E: Exposed, I: Infected, R: Resistant) for an individual without maternally
122
acquired antibodies, (b) the respective dynamics for individuals with maternally acquired antibodies (M:
123
Maternal antibodies full protection from infection and disease in juveniles) is illustrated in panel (b); The
124
subscript juv indicates that an individual is juvenile (< 12 weeks old), for which changes in susceptibility to
125
diseases is a function of age and the parameters and (i.e. increase in RHD, decrease in Myxo). Disease
126
exposure of susceptible individuals depends on the transmission rate β(I) and/or the probability [t] that a virus
127
is introduced at time t, capturing the seasonality of virus activity. For infected juveniles, the recovery rate is a
128
function of the dynamical state of juvenile susceptibility and adult recovery rate . For individuals with maternal
129
antibodies, changes in juvenile susceptibility and recovery rates are also dependent on the parameter , which
130
describes the linear decline of maternal antibodies of aging juveniles. The parameters P and I describe the time
131
length of transition steps. Note that individual fate in our model depends also on demographic dynamics and the
132
disease-induced mortality from the co-circulating virus. In particular, the interaction with seasonal birth rates
133
determines the phenological matching between demographic and epidemiological dynamics.
134
135
136
137
Using R for sampling and simulation in Meta Model Manager
5
6
Supporting Information
138
The utilised software can be both operated from ready-to-use user interfaces or, alternatively, via
139
command-line operation. We ran Meta-Model Manager from the freely available software R [8]. In
140
brief, we defined relevant parameters as either single numerical constants (i.e. fixed parameters) or
141
vectors with variable values for each simulation sampled from a latin hypercube with the R package
142
LHS [9]. We established template files for each submodel (1 Vortex, 2 Outbreak, and 1 Meta-Model
143
Manager for the study). From the R environment, we then loaded these template files, replaced the
144
values for key parameters for each simulation and saved the edited files with the simulation number as
145
an index. The full set of simulations was then operated with system commands, i.e. a command that
146
runs Meta-Model Manager and calls for each simulation. By scanning the output files with R, we
147
assembled output values such as the number of individuals in different disease states each day of the
148
different simulations into arrays that allowed a large range of subsequent analysis and graphical
149
display. The code for repeating our study is available in the Vortex library at
150
www.vortex10.org/Downloads/OzRabbitDzFiles.zip.
151
152
153
154
155
156
Model specification
157
Table A.1. First-order independent model parameters and their ranges (minimum/ maximum) used for
158
simulations. For each simulation, key parameters of most interest were drawn from a hypercube
159
(‘Sample’; abbreviation given as used in the main article). Some parameters are further varied over
160
time steps by sampling from uniform distributions over the defined ranges (indicated with ‘uniform’).
Parameter name
Range / Unit
Sampling/
Description
uncertainty
Justification/
Reference
Rabbit demography
Survival probability
juv: 0.42
subad: 0.55
ad: 0.75
Sample
Survival probability
(annual rate)
depending on age
with juveniles < 12
weeks, subadults 1224 weeks, adults > 24
weeks.
[10]
7
Supporting Information
Annual peak in reproduction
(RepPeak)
20 – 50
Sample
[calendar week]
Calendar week of
Field data, [11, 12]
maximum
reproduction (used as
mean in normal
distribution to fit
relative reproductive
efforts over weeks.
Variation in reproduction
(RepVar)
2 – 10
Sample
[calendar week]
Variation in relative
Field data, [11, 12]
weekly reproductive
efforts in different
simulations (used as
SD in normal
distribution).
Annual reproductive effort
(RepEff)
Litter size
50 – 200 [%
Sample
Total reproductive
female
effort [%] in relation
reproducing]
to populations size N.
1-5
uniform
Number of offspring
[12, 13]
Expert guess
per litter.
Reproductive age
26 [weeks]
Age of first
Field data
reproduction.
Maximum age
416 [weeks]
Maximum lifetime of
[14]
rabbits in calendar
weeks.
Carrying capacity
500
Number of
[individuals]
Field data, [10]
individuals defining
carrying capacity; we
assumed density
dependent
reproduction, scaled
by (K-Nt)/K.
Environmental stochasticity
(EV)
0.01 – 0.5 [SD
Sample
Explorative
of vital rates,
weekly]
Disease epidemiology
Transmission probability (βV)
0.3 – 0.9
Sample
Probability that a
Explorative
virus is transmitted
from an infected
individual.
Maternal antibodies (tMV)
1 – 50
RHD: Sample
Time period of
Explorative
8
Supporting Information
[ days]
Myxo: none
maternal antibody
protection against
infection and disease.
Insusceptibility (newborns)
RHD: 21
Time period before
Myxo: 1
newborns may
[days]
Exposure period
RHD: 1
[6, 7]
become susceptible.
uniform
[6]
uniform
[6]
Myxo: 1 – 4
[days]
Infection period
RHD: 1-2
Myxo: 8 – 12
[days]
Recovery rate (V)
0.2 – 0.9
Sample
Probability of
Explorative
survival of infected
individuals > 90 days
old.
Juvenile susceptibility factor
(Juv,RHD, μJuv,Myxo)
Juv,RHD: - 0.1 –
Sample
Decreasing (RHD) /
Explorative,
-1
increasing (Myxo)
similar dynamics
μJuv,Myxo: -6 - 6
recovery of infants
described in [5]
aged 1-90 days,
modelled as a
coefficient in logitlink logistic model.
Virus introduction rate (pIntrV)
0 – 0.1 [%]
Sample
Introduction of
Explorative
viruses from external
sources such as
arthropod vectors or
carcasses.
First calendar week of virus
RHD: 1 – 52
introduction (wkIntroV)
Myxo: 1 – 52
Sample
First calendar week
Explorative
each year in which
[calendar week]
viruses are
introduced/ host are
exposed to the virus.
Last calendar week of virus
RHD: 1 – 52
introduction
Myxo: 1 – 52
Sample
Last calendar week a
Explorative
virus may be (value
[calendar week]
first week of
introduction).
Exposure time
RHD: 1
Myxo: 1 - 4
Infection time
RHD: 1 - 2
uniform
[15]
9
Supporting Information
Myxo: 8 12
[days]
Model initialisation
Initial N
200
Initial population size
arbitrary (first 5
for all simulations.
years of simulation
not considered in
analysis)
Initial infection rate
5%
Proportion of
arbitrary (first 5
individuals of N0
years of simulation
infected.
not considered in
analysis)
161
162
163
164
165
166
167
168
References
169
[1] Lacy, R.C., Miller, P.S., Nyhus, P.J., Pollak, J.P., Raboy, B.E. & Zeigler, S.L. 2013
170
Metamodels for transdisciplinary analysis of wildlife population dynamics. PLoS ONE 8,
171
e84211. (doi:10.1371/journal.pone.0084211).
172
[2] Lacy, R.C. & Pollak, J.P. 2014 VORTEX: A stochastic simulation of the extinction
173
process. Version 10.0. (Brookfield, Illinois, USA, Chicago Zoological Society.
174
[3] Lacy, R. 2000 Structure of the VORTEX simulation model for population viability
175
analysis. Ecological Bulletins 48, 191-203.
176
[4] Bradshaw, C.J.A., McMahon, C.R., Miller, P.S., Lacy, R.C., Watts, M.J., Verant, M.L.,
177
Pollak, J.P., Fordham, D.A., Prowse, T.A.A. & Brook, B.W. 2012 Novel coupling of
178
individual-based epidemiological and demographic models predicts realistic dynamics of
179
tuberculosis in alien buffalo. Journal of Applied Ecology 49, 268-277. (doi:10.1111/j.1365-
180
2664.2011.02081.x).
Supporting Information
181
[5] Robinson, A.J., So, P.T.M., Muller, W.J., Cooke, B.D. & Capucci, L. 2002 Statistical
182
models for the effect of age and maternal antibodies on the development of rabbit
183
haemorrhagic disease in Australian wild rabbits. Wildlife Research 29, 663-671.
184
(doi:10.1071/wr00119).
185
[6] Kerr, P.J. 2012 Myxomatosis in Australia and Europe: a model for emerging infectious
186
diseases. Antiviral Research 93, 387-415. (doi:10.1016/j.antiviral.2012.01.009).
187
[7] Abrantes, J., van der Loo, W., Le Pendu, J. & Esteves, P.J. 2012 Rabbit haemorrhagic
188
disease (RHD) and rabbit haemorrhagic disease virus (RHDV): a review. Veterinary
189
Research 43. (doi:10.1186/1297-9716-43-12).
190
[8] R Development Core Team. 2013 R: A language and environment for statistical
191
computing. Vienna, Austria, R Foundation for Statistical Computing. http://www.
192
http://cran.r-project.org/
193
[9] Carnell, R. 2012 lhs: Latin hypercube samples. R package version 0.10.
194
[10] Fordham, D.A., Sinclair, R.G., Peacock, D.E., Mutze, G.J., Kovaliski, J., Cassey, P.,
195
Capucci, L. & Brook, B.W. 2012 European rabbit survival and recruitment are linked to
196
epidemiological and environmental conditions in their exotic range. Austral Ecology 37, 945-
197
957. (doi:10.1111/j.1442-9993.2011.02354.x).
198
[11] Gilbert, N., Myers, K., Cooke, B.D., Dunsmore, J.D., Fullagar, P.J., Gibb, J.A., King,
199
D.R., Parer, I., Wheeler, S.H. & Wood, D.H. 1987 Comparative dynamics of Australasian
200
rabbit populations. Australian Wildlife Research 14, 491-503.
201
[12] Mutze, G., Bird, P., Kovaliski, J., Peacock, D., Jennings, S. & Cooke, B. 2002 Emerging
202
epidemiological patterns in rabbit haemorrhagic disease, its interaction with myxomatosis,
203
and their effects on rabbit populations in South Australia. Wildlife Research 29, 577-590.
204
(doi:10.1071/wr00100).
205
[13] Mutze, G.J., Sinclair, R.G., Peacock, D.E., Capucci, L. & Kovaliski, J. 2014 Is increased
206
juvenile infection the key to recovery of wild rabbit populations from the impact of rabbit
207
haemorrhagic disease? European Journal of Wildlife Research 60, 489-499.
208
(doi:10.1007/s10344-014-0811-6).
209
[14] Peacock, D.E. & Sinclair, R.G. 2009 Longevity record for a wild European rabbit
210
(Oryctolagus cuniculus) from South Australia. Australian Mammalogy 31, 65-66.
211
[15] McPhee, S.R., Butler, K.L., Kovaliski, J., Mutze, G., Capucci, L. & Cooke, B.D. 2009
212
Antibody status and survival of Australian wild rabbits challenged with rabbit haemorrhagic
213
disease virus. Wildlife Research 36, 447-456. (doi:10.1071/wr08137).
10
Supporting Information
214
View publication stats
11