流感模型MSEIRS(免疫-易感-潜伏-传染-康复-易感)

“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,