“In Tables 2-3 the 95% credibility interval for SAREMCMC is much wider than for MHMCMC.
We can see that the upper bound values for z (the proportion susceptible before the summer
wave of the all population) are very high (92.1% and 99.7%) for the Leicester and Wigan
populations; the upper bound values of α (the proportion of becoming infective and
symptomatic) are also very high (96.3% and 99.7% for summer wave, 99.1% and 100% for
autumn wave, 99.7% and 99.8% for winter wave) for the Leicester and Wigan populations;
this agrees with the assumption of [14, 25, 26] that with a new virus, the entire population is
susceptible. However, the local search of MHMCMC conflicts with this assumption.”
----------- Global Optimization Search differs from Local Optimization Search very much
“(1) such a complex MSEIRS model is highly over-parameterized with respect to the data
and the estimates are highly correlated; (2) the model should consider the rich data of age
classes ([27]); and (3) the estimates obtained (such as R0) need to be properly compared to
existing estimates in literatures (e.g. [14, 21, 24]) argued.”
“This paper extends M to the prior immune in a very general sense as in [21, 23, 24]. The
prior immune (M) of this paper is not only gained from infants’ mother, but also gained from
seasonal influenza, vaccination, etc ([21, 23, 24] and references therein);”
World Journal of Modelling and SimulationVol. 7 (2011) No. 1, pp. 29-39. ISSN
1 746-7233, England, UK. Submitted 19-NOV-2009, revised on 05-JAN-2010 & 01-MAY-2010,
accepted on 07-JUN-2010, published online 28-JAN-2011:
http://www.wjms.org.uk/wjmsvol07no01paper03.pdf
An Effective Simulated Annealing Refined Replica Exchange Markov
Chain Monte Carlo Algorithm for the Infectious Disease Model of
H1N1 Influenza Pandemic
Jiapu Zhang
Centre for Informatics and Applied Optimization & Graduate School of ITMS, The
University of Ballarat, Mount Helen, Ballarat, VIC 3353, Australia
Mobile: (61) 423 487 360, E-mail: jiapu zhang@hotmail.com, j.zhang@ballarat.edu.au
Abstract: This paper is concerned with a computational algorithm for fitting a deter-
ministic MSEIRS (immune-susceptible-exposed-infectious-recovered-susceptible) epidemic
model for the transmission of influenza (H1N1) to mortality data. The model-fitting is
carried out using a simulated annealing refined replica exchange Markov chain Monte
Carlo algorithm. The effectiveness of the algorithm is illustrated using the triple wave
data from five English towns collected during the 1918-19 influenza pandemic. Numer-
ical results show that the replica exchange (refined by simulated annealing) sampling
technique is superior to other existing sampling techniques such as the Gibbs sampling
technique, the Metropolis-Hastings sampling technique, the Multiple-try Metropolis
technique for the Markov chain Monte Carlo computation. The algorithm presented in
this paper has great promise to be used for carrying out some numerical computations
of the current complex 2009-10 influenza pandemic.
Key words: Replica Exchange; Simulated Annealing; Markov Chain Monte Carlo;
Metropolis-Hastings Sampling; Gibbs Sampling; Multiple-try Metropolis; Extended
MSEIRS Model.
1 Introduction
The world is facing this century’s first influenza pandemic caused by the recent outbreak
of H1N1 swine influenza. In fact H1N1 swine flu pandemic once happened in multiple
waves in 1918-19. A detailed epidemic model for H1N1 swine flu may be demonstrated
as follows ([1, 5, 6, 9, 10, 11, 12, 17, 21, 22, 29, 31, 37]):
M S E1 E2
A
I D
R L
T
Tw λ
γ
(1−α)γ
αγ
ν
ν
µ
ρφ
(1−ρ)φ
φ
1
Extended MSEIRS Model. Within a period of time in resistant state (Tw) the prior immune
individuals (M) become into susceptible individuals (S), then the susceptible individuals are ex-
posed into two sequential latent infections (E1 and E2 with infectious rate γ) as the results of the
exposure to the force of infection (λ). A proportion α then become infective and symptomatic
(I, i.e. clinically ill and infectious) and a proportion 1-α become asymptomatic or unpapered
infective (A). Both I and A pass to the recovered class with immunity (R) with recovery rate
ν; some individuals of I class become in the death class (D) with a proportion µ of infected
individuals who die. From R a proportion ρ develop longer-lasting protection (L), while the re-
mainder pass through a temporary state (T) and eventually become into S again with φ = 2/Tw.
The former achievements on MSEIRS model ([1, 5, 6, 9, 10, 11, 12, 17, 22, 29, 31, 37])
are mainly on the passively immune class M of infants. This paper extends M to the
prior immune in a very general sense as in [21, 23, 24] and is concerned with a practical
computational algorithm for the Extended MSEIRS Model. The prior immune (M)
of this paper is not only gained from infants’ mother, but also gained from seasonal
influenza, vaccination, etc ([21, 23, 24] and references therein); this is a new research
direction of the studies of the spread of H1N1 in the influenza pandemic, though it is
not our interests of this paper.
The aim of this paper is to construct and demonstrate a computational algorithm
for fitting the above extended MSEIRS epidemic model to mortality data. In Section 2
of this paper, a deterministic ordinary differential equation (ODE) mathematical model
is built and a simulated annealing refined replica exchange Markov chain Monte Carlo
(REMCMC) algorithm is designed. The effectiveness of the algorithm is illustrated us-
ing the triple wave data from five English towns collected during the 1918-19 influenza
pandemic ([27]). The MCMC ([2, 4, 5, 15]) model-fitting algorithm is used to esti-
mate parameter distributions of the model. The statistical sampling technique of RE
([13, 18, 19, 32, 33]) is enclosed into the MCMC algorithm, where a simulated annealing
scheme ([3, 16]) is employed to refine the RE. Numerical results in Section 3 of this
paper show that the simulated annealing-refined RE sampling technique is superior to
the Metropolis-Hastings (MH) sampling technique ([4], from which we may know that
Gibbs sampling can be rewritten as a Metropolis algorithm; thus, the RE sampling
is superior to the Gibbs sampling too), the Multiple-try Metropolis (MM) sampling
technique for the Markov chain Monte Carlo computation of the extended MSEIRS
epidemic model. The simulated annealing refined REMCMC algorithm presented in
this paper has great promise to be used for carrying out some numerical computations
of the current complex 2009-10 influenza pandemic.
The organization of this paper is as follows. The ODE mathematical model for the
Extended MSEIRS Model is built in Section 2, and its solver, the simulated annealing
refined replica exchange Markov chain Monte Carlo algorithm, is presented in Section 2.
Section 3, using the triple wave data from five English towns collected during the 1918-
19 influenza pandemic, tests the algorithm for the ODE model. Numerical results are
shown and discussed in Section 3. Compared with other existing sampling algorithms
such as the Metropolis-Hastings algorithm (where the Gibbs sampling is its special
case), the Multiple-try Metropolis algorithm, the algorithm presented in this paper is
2
concluded in Section 4 as a promising algorithm to be used to study the H1N1 influenza
pandemic.
2 Model Building and Solving
R0 denotes the basic reproduction number ([1, 14, 26]). We assume that the infections
of A do not transmit and both A and I occur in proportion, M individuals have prior
immunity including the innate immunity for children, some S individuals who were
exposed in E1 and E2 will develop into A or I and then become into L individuals
and we allow for the possibility that protection may be lost and R individuals become
into S individuals again because immunity waned in individuals in waves or antigenic
drift of the virus. We also assume the populations are homogeneous mixing between
individuals. With these assumptions, we will build the model by the following ODEs
([5, 8, 9, 10, 12, 17, 20, 21, 22, 37]) and then solve it by the simulated annealing refined
REMCMC algorithm.
2.1 An ODE Model
dM
dt
= −
2
φ
MS − µM (1)
dS
dt
= −
R0ν
Nα
IS + φT (2)
dE1
dt
=
R0ν
Nα
IS − γE1 (3)
dE2
dt
= γE1 − γE2 (4)
dA
dt
= (1− α)γE2 − νA (5)
dI
dt
= αγE2 − νI (6)
dD
dt
= µI (7)
dR
dt
= ν(I +A)− φR (8)
dT
dt
= (1 − ρ)φR− φT (9)
dL
dt
= ρφR (10)
M(0) = S(0) = Nz,R(0) = N(1− z), I(0) = I0, . . .
E1(0) = E2(0) = A(0) = T (0) = L(0) = 0. (11)
where N (assumed fixed) is the total number of individuals in the population for the
outbreak in question, z is the proportion initially susceptible, I0 is the given data which
is dependent on the population, Re = zR0 gives the initial effective reproduction num-
ber, and µ = D/AR where AR is the attack rate calculated by AR = αz(1−e−R0AR/α).
3
In order to make the reader know how the above formulas come from, we describe
every formula (according to the chart at the beginning of Section 1) as follows. Equation
(1) describes the variation of the prior immune individuals M, a proportion µ of M were
deaths and the remaining became into the susceptible individuals S within the period
of time in resistant state Tw, where Tw = 2/φ. The force-of-infection λ is given by
λ = R0νNα I. Formulas (2) and (3) can be respectively rewritten as
dS
dt = −λS+φT,
dE1
dt =
λS − γE1. Formula (2) explains the variation of susceptible individuals S: the increase
of S from the recovered individuals with temporary immunity φT and the decrease of
S because of the force of infection λ. Formula (3) describes the variation of individuals
at the first latent state E1 with the latent infectious rate γ. Formula (4) describes the
variation of individuals at the second exposed state E2 of latent infection. Formula
(5) describes the variation of the asymptomatic or unreported infectious individuals
A: a proportion 1 − α of E2 with the latent infectious rate γ become into A and the
recovery rate of A is ν. Similarly as formula (5), formula (6) describes the variation
of the remaining proportion α of E2 which become into the infective and symptomatic
class of individuals I, with the recovery rate ν. Formula (7) describes a proportion µ
of I become deaths. Formula (8) describes the variation of the recovered individuals
R: some I and A individuals recovered with rate ν, but a proportion φ of R become
into the class of susceptible individuals S again. Formula (9) describes (1 − ρ)φ of
R individuals were not susceptible temporarily and a proportion φ of the temporary
class of individuals T become into S soon. Formula (10) describes ρφ of recovered
individuals R become into the class of longer-lasting protected individuals L. Formula
(11) describes the initial values of the variables. The descriptions of formulas (2)-(11)
can be referred to the references such as [5, 8, 9, 10, 12, 17, 20, 21, 22, 37].
2.2 An Simulated Annealing Refined Replica Exchange Markov Chain
Monte Carlo Algorithm
MCMC algorithms are sampling from probability distributions based on constructing
a Markov chain that has the desired distribution as its equilibrium distribution ([37]).
The sampling strategy is very critical for a successful MCMC algorithm. However, in
practice, the MCMC sampling methods such as Gibbs sampling, Metropolis-Hastings
algorithm ([4]), Multiple-try Metropolis algorithm sometimes just randomly walk and
take a long time to explore all the solution space, will often double back and cover
ground already covered, and usually own a slow algorithm convergence. In this paper
a more efficient sampling strategy of simulated annealing-refined RE is enclosed into
the MCMC simulation.
In the Metropolis-Hastings-based MCMC (MHMCMC), a Markov stochastic pro-
cess is built to sample a target probability p(x) = C−1e−f(x) for a current state x,
where f(x) is the objective function (which is the likelihood in this paper) and C
denotes the normalization constant. For simulated annealing, a variable tempera-
ture Temp is introduced into the objective function of the target distribution, i.e.
p(x) = C−1e−f(x)/Temp. For RE, a new state y is generated from the current state
x of the Markov process by drawing y with a transition probability q(x, y), and the
new state y is accepted with the probability min(1, [p(y)q(y, x)]/[p(x)q(x, y)]). In the
4
RE implementation, for a replica j at each iteration step l, the local Markov chain
move from the conformation state xlj to the new state x
l+1
j accepted with the prob-
ability min(1, e−1/Temp(j)(f(x
l+1
j )−f(x
l
j))), and the replica transition from j to j + 1 at
two neighbouring temperatures Temp(j) and Temp(j+1) with the accepted probabil-
ity min(1, e−1/Temp(j)f(x
l+1
j+1
)−1/Temp(j+1)f(x
l+1
j
)+1/Temp(j+1)f(x
l+1
j+1
)+1/Temp(j)f(x
l+1
j
)
). The temperatures
are monotonically decreasing in order for the convenience of simulated annealing in use
and the transition step size of RE is larger for higher temperature and smaller for the
lower temperature. The negative binomial distribution is used to calculate the max-
imum likelihood, where the likelihood is formed by assuming Gaussian errors around
the deterministic model. For the simulated annealing refinement for RE, the neign-
bourhood scheme of [3] is still in use here. The pseudo-code of the simulated annealing
refined REMCMC algorithm is presented as follows:
Define random and constant parameters
Initialization:
offset =0
for (i=1 to replicas) do
Set up MCMC random and constant Parameters for each replica
Temperature scheme to produce Temp(i)
MCMC Initialized at Temp(i) to get MSEIRSinit
endfor
Do MCMC with RE:
for (i set=1 to Total set) do
Set(i) of Total set:
for (j=1 to replicas) do
Input Parameters and MSEIRSinit
Call MCMC to get new Parameters and each likelihood, accepting new MC move with
Metropolis criterion for each replica
Output new Parameters, and the Likelihood last(j), Likelihood last cold(j)
endfor
j = offset +1
while (j+1 ≤ replicas) do
k=j+1
∆ = ( 1Temp(k) −
1
Temp(j) )*(Likelihood last cold(j)−Likelihood last cold(k))
if (∆ ≤ 0) then
Output Likelihood last(j)
swap Parameters and Labels of j and k
else
Generate a random number Rand of [0,1]
if (Rand ≤ e−∆) then
Output Likelihood last(j)
5
swap Parameters and Labels of j and k
endif
endif
j=j+2
endwhile
offset = 1- offset
endfor
Adjust each of Temps by simulated annealing & repeat the algorithm untill reach Temp target.
3 Numerical Results and Discussion
The data used to test the above simulated annealing refined REMCMC (SAREMCMC)
algorithm for fitting the extended MSEIRS Model for the transmission of the H1N1 in-
fluenza is the mortality data from five English towns (Blackburn, Leicester, Manchester,
Newcastle and Wigan) collected during the 1918-19 influenza pandemic ([27]). In the
TABLES 1, 2A-2N, 3-5 of [27] pages 136-143, the details of the triple pandemic wave
data in 1918-19 of English localities are given. This paper just uses the mortality data
(i.e. the 1st, 2nd, 3rd ... 48th week of the deaths of people in each of the five places),
which are listed in Table 1.
10 replicas and 10000 iterations are set for the algorithm, Total set=50. To com-
pare with the SAREMCMC algorithm, the MHMCMC algorithm ([4, 37]) is run for
10000 iterations too. The model curves fitted to the data are shown in Figures 1-5,
from which we can see that the model mortality curve well fits the death data for each
of the five populations by the SAREMCMC algorithm, but the MHMCMC algorithm
can only well fit for the Leicester and Wigan populations. The comparisons of the
successful parameter maximum likelihood estimates and credibility intervals calculated
by the both algorithms of 10000 iterations are listed in Tables 2-3.
In Tables 2-3 the 95% credibility interval for SAREMCMC is much wider than for
MHMCMC. We can see that the upper bound values for z (the proportion susceptible
before the summer wave of the all population) are very high (92.1% and 99.7%) for
the Leicester and Wigan populations; the upper bound values of α (the proportion of
becoming infective and symptomatic) are also very high (96.3% and 99.7% for summer
wave, 99.1% and 100% for autumn wave, 99.7% and 99.8% for winter wave) for the
Leicester and Wigan populations; this agrees with the assumption of [14, 25, 26] that
with a new virus, the entire population is susceptible. However, the local search of
MHMCMC conflicts with this assumption. For the both populations, the effective re-
production number Re is much less than R0 so that the mortality rate is no more than
3.1% as low as reported ([27]). To further confirm the above observations, the MHM-
CMC was run 50000 iterations and the MMMCMC (Multiple-try Metropolis Markov
Chain Monte Carlo) was run for 14999 iterations. In Figures 1-5 we may see that the
performance of MHMCMC of 50000 iterations and of MMMCMC of 14999 iterations
are still not better than that of SAREMCMC of 10000 iterations. This further confirms
6
that REMCMC is much better than MHMCMC and MMMCMC for solving the H1N1
influenza pandemic model.
The detailed analysis for Tables 2-3 is not the interests of this paper because (1)
such a complex MSEIRS model is highly over-parameterized with respect to the data
and the estimates are highly correlated; (2) the model should consider the rich data
of age classes ([27]); and (3) the estimates obtained (such as R0) need to be properly
compared to existing estimates in literatures (e.g. [14, 21, 24]) argued. This will be a
further research direction of the author.
4 Conclusion
This paper is concerned with a computational algorithm for fitting the deterministic
MSEIRS epidemic model for the transmission of H1N1 influenza to mortality data. The
model-fitting is carried out using a simulated annealing refined replica exchange stochas-
tic algorithm. The algorithm is superior to other existing sampling algorithms such as
the Gibbs sampling, the Metropolis-Hastings algorithm, the Multiple-try Metropolis
algorithm through the illustration of using the triple wave mortality data from five En-
glish towns collected during the 1918-19 influenza pandemic. The algorithm presented
in this paper has great promise to be used for carrying out some numerical computa-
tions of the current complex 2009-10 influenza pandemic.
Acknowledgments: The author appreciates the anonymous referees and editors for the nu-
merous insightful comments (05-JAN-2010) on the original manuscript and further comments
(01-MAY-2010) on the revised manuscript, which have improved this paper greatly.
References
[1] Anderson, R., May, R.M. (1992) Infectious diseases of humans: dynamics and
control, Oxford University Press, Oxford.
[2] Andrieu, C., De Freitas, N., Doucet, A., Jordan, M.I. (2003) An Introduction to
MCMC for machine learning. Mach. Learn. 50, 5–43.
[3] Bagirov, A.M., Zhang, J.P. (2003) Comparative analysis of the cutting angle and
simulated annealing methods in global optimization. Opt. 52, 363–378.
[4] Baldi, P., Brunak, S. (2001) Bioinformatics: The Machine Learning Approach
(2nd Edition), MIT Press, Chapter 4.5.
[5] Bootsma, M.C.J., Ferguson, N.M. (2007) The effect of public health measures on
the 1918 influenza pandemic in U.S. cities. P.N.A.S. USA 104 (18), 7588-93.
[6] Bos, M.E.H. et al (2007) Estimating the day of highly pathogenic avian influenza
(H7N7) virus introduction into a poultry flock based on mortality data. Vet. Res.,
38, 493–504.
7
[7] Caley, P., Philp, D.J., McCracken, K. (2008) Quantifying social distancing arising
from pandemic influenza. J. R. Soc. Interface 5, 631–9.
[8] Chen, S.C., Chio, C.P., Jou, L.J., Liao, C.M. (2009) Viral kinetics and exhaled
droplet size affect indoor transmission dynamics of influenza infection. Indoor
Air, 19, 401–13.
[9] Chen, S.C., Liao,