Revealing the True Incidence of Pandemic A(H1N1)pdm09 Influenza in Finland during the First Two Seasons — An Analysis Based on a Dynamic Transmission Model

The threat of the new pandemic influenza A(H1N1)pdm09 imposed a heavy burden on the public health system in Finland in 2009-2010. An extensive vaccination campaign was set up in the middle of the first pandemic season. However, the true number of infected individuals remains uncertain as the surveillance missed a large portion of mild infections. We constructed a transmission model to simulate the spread of influenza in the Finnish population. We used the model to analyse the two first years (2009-2011) of A(H1N1)pdm09 in Finland. Using data from the national surveillance of influenza and data on close person-to-person (social) contacts in the population, we estimated that 6% (90% credible interval 5.1 – 6.7%) of the population was infected with A(H1N1)pdm09 in the first pandemic season (2009/2010) and an additional 3% (2.5 – 3.5%) in the second season (2010/2011). Vaccination had a substantial impact in mitigating the second season. The dynamic approach allowed us to discover how the proportion of detected cases changed over the course of the epidemic. The role of time-varying reproduction number, capturing the effects of weather and changes in behaviour, was important in shaping the epidemic.

, D (hosp) a,t , D (IC) a,t } a∈1...16; t∈0...112 . Above a is the age group index and t is the week index. We solve the computational problem using the following algorithms: • Markov chain Monte Carlo (MCMC). It was used as a main inference tool.
• Iterative component-wise optimization. Auxiliary methods. Used to find the maximum of the posterior. The maximum is then used as the initial point for MCMC. Optimization was also used for the sensitivity analysis.
• Exact approximate MCMC (Andrieu and Roberts, 2009) targeting the tempered marginal posterior distributions p(θ|D) 1 5 and p(θ|D) 1 25 . Auxiliary methods. The tempered posteriors were studied to ensure that the peak area of the target posterior is unimodal and wellbehaving.
• Importance sampler targeting the likelihood p(D|θ). Auxiliary methods, used within Exact approximate MCMC.
The following text gives an overview of the computational methods. All probabilities and importance weights were log-transformed. The model parameters p, s (sev/inf) , s (IC/sev) , d (hosp) were explored in a logit form. Parameter q was explored in a log form. The model parameters w t and d (mild) t were deterministicly transformed parameters ε and σ, and therefore not taken into analysis explicitly. We use the following notation: N (µ, σ 2 ) is a normal distribution with mean µ and variance σ 2 ; N(µ, Σ) is a multivariate normal distribution with mean µ and covariance matrix Σ; U (0, 1) is a number sampled uniformly from the interval (0, 1), '←' is an assignment operator.

MCMC
We used MCMC as a main inference tool targeting the full posterior of the model unknowns p(θ, H|D). The chain was initiated with the maximum estimated with optimization (see below). The parameters ε, σ and the hidden states H were updated in two different ways during the odd and even iterations to improve mixing.
The model parameters p, q, s (sev/inf) , s (IC/sev) were updated using the following Metropolis-Hastings steps: 1. The proposed value for p is sampled in logit form logit(p * ) ∼ N(logit(p), Σ p ), here Σ p is the proposal covariance matrix for the parameter logit(p).

The proposal is accepted
here π() is a prior for logit(p).
3. The parameters q, s (sev/inf) and s (IC/sev) were updated in the same way using steps 1-2. For the parameter q the log transformation is used instead of logit.
The detection probability for a hospitalized case d (hosp) was updated using the following Metropolis-Hastings steps: 1. The proposed values are sampled in logit form here π() is a prior for logit (d (hosp) ).
On the odd iterations, parameters ε, σ were updated using the following Metropolis-Hastings steps: 1. Randomly sample two integers i and j, such that 0 < i < j < T , here T = 113 weeks is length of the studied period. 4. The head and tail parts of the vector ε in the same way using the step 3: 5. The parameter σ is updated in the same way using steps 1-3.
On the even iterations, parameters ε, σ were updated using the following approximation to the Metropolis-Hastings steps.
1. The proposed value is sampled here Σ ε is a proposal covariance matrix for the parameter ε.

The first component of the proposal is accepted
here π() is the prior.

The second component of the proposal is accepted
4. The process is sequentially repeated for all components of ε.
5. The parameter σ us updated in the same way using steps 1-4.
On the odd iterations, hidden states were updated using a particle Gibbs sampler (Andrieu et al., 2010) step with n = 200 particles. The importance function used in the Gibbs step is described in the Importance Sampler section.
1. Initialize the set of n particles Initialize importance weights: and g is the importance function.
The first particle plays the role of a reference. It is not sampled randomly, but is set to be equal to the hidden state H from the current state of the chain: During iterations t ∈ 20, 25, 30 . . . 110, after the particle and weights are propagated, conduct resampling. Resample the set of particles using their normalized weights. Keep the first (reference) particle intact.
Resampling is done using a resampling scheme (Douc et al., 2005) 5. After all weeks are sampled, select one of the particles using the normalized weights to be the proposal. Proposal is always accepted.
On the even iterations, hidden states were updated using the algorithm: 1. Sample a proposal for the hidden states of the first age strata, keeping the rest of the hidden states intact: ..16 , S 2...16 , D, θ); Here I a = {I a,t } t∈0...T , S a = {S a,t } t∈0...T and g is the importance function.
2. The proposal is accepted ( .
3. Each other stratum is updated in the same way using steps 1-2 .
We used 15 independently run chains, each with 55000 iteration, discarding the initial 15000 iterations as warm-up and recording every 20 th iteration. Each iteration included updating p, q, s (sev/inf) , s (IC/sev) , d (hosp) , ε, σ and the hidden states. Proposal variances Σ and σ 2 were adopted during the warm-up phase to obtain an approximate acceptance rate of 1/3 (Brooks et al., 2011).

Optimization
Before starting the MCMC simulation to explore the posterior, the maximum of the posterior first was searched by an ad hoc optimization algorithm. The MCMC chains were then initiated with the found maximum point. In the sensitivity analyses, the results were compared in terms of the posterior maximum points found by the algorithm. The optimization is initiated at a random point. For each parameter x ∈ θ or hidden value x ∈ H the proposal distance step x is set to the arbitrary non-zero value. Parameters and hidden values are updated in a slightly different way because all parameters are continuous while all hidden values are discrete.
The hidden values x ∈ H are updated in the following way: 1. Let x 0 denote the current value of the hidden value x. If x = 0, set step x = 1 2. A new value for x is proposed: 3. The part of the posterior density that depends on x is evaluated for the two values of x: 4. if P 0 > P * , keep the value x = x 0 and set step x ← round(−0.5 step x + 0.6) 5. if P * > P 0 , set the value x ← x * and set The parameters x ∈ θ are updated in the following way: 1. Let x 0 denotes the current value of the hidden value x. Two new values for x are proposed: The part of the posterior depending on x is evaluated for the different values of x: 3. if P 0 > P * , P * * keep the value x = x 0 and set (a) if P 0 = P * * , set step x ← +0.1 step x ; (b) else, set step x ← P * * − P 0 2 max(P * − P * * , P * − P 0 ) step x 5. if P * * > P 0 , P * set the value x ← x * * ; set We run the optimization algorithm for 5000 iterations.

Exact approximate MCMC
We used Exact approximate MCMC (Andrieu and Roberts, 2009) to target the tempered marginal posterior p(θ|D) 1/25 . The chain was initiated with the maximum estimated with the optimization (see above). The algorithm proceeded by sequentially updating the components of θ. The sampler was used twice with the temperature set to temp = 1/25 and temp = 1/5. The model parameters p, q, s (sev/inf) , s (IC/sev) were updated using the following Metropolis-Hastings steps: 1. The proposed value of p is sampled at the logit scale: logit(p * ) ∼ N(logit(p), Σ p ), here Σ p is the proposal covariance matrix for logit(p).
2. The likelihood for the proposal is estimated with importance sampling L * ←p(D|p * , . . . ).
4. The parameters q, s (sev/inf) and s (IC/sev) are updated in the same way using steps 1-3. For the parameter q the log transformation is used instead of logit.
The detection probability for the hospitalized cases d (hosp) was updated using the following Metropolis-Hastings steps: 1. The proposed values are sampled in logit form here σ 2 d (hosp) is the proposal variance. 2. The likelihood for the proposal is estimated with importance sampling L * ←p(D|d (hosp) * , . . . ).

The proposal is accepted
The parameters ε, σ (determining random effects w and d) were updated using the algorithm: 1. Randomly sample two integers i and j, such that 0 < i < j < T , here T = 113 here T = 113 weeks is length of the studied period.

5.
Steps 3-4 are repeated for the head and tail parts of the vector ε: 6. The parameter σ is updated in the same way using steps 1-5.
We drew 6000 samples, discarding the initial 2000 as warm-up. Each iteration includes updating p, q, s (sev/inf) , s (IC/sev) , d (hosp) , ε and σ. At the end of each non warm-up iteration the parameter vector was saved as one sample from the posterior p(θ|D) temp . The proposal variances Σ and σ 2 were adopted during the warm-up phase to get close the optimal acceptance rate 1/3 (Brooks et al., 2011). so work distribution proportion for system with Q accelerators should be W CP U : W ACC Every worker uses OpenMP to expose parallelism of the loop doing Monte Carlo casts. We used offload acceleration model to support Intel Xeon PHI capabilities. The code performing Monte Carlo cast is the same for both microarchitectures (x86 and MIC). Since we run bunch of separate computational threads, we need bunch of PRNG sequences considered to be uncorrelated to get necessary source of randomness (dedicated sequence for each thread). We use Intel MKL Mersenne Twister generator V SL BRN G M T 2203 with 6024 pregenerated sequences designed specifically for massively parallel computations (this implies limit on parallelismsampler can not utilize more than 6024 computational cores).
The rest of the code was written in Python 2.7 with external libraries NumPy and SciPy for mathemathical function and Pillow for visualization. All algorithms were executed on the RSATU server. The main MCMC algorithm took 1 month. The auxiliary MCMC took about 3 days each. The optimization required 2 days. RSATU server is NUMA-machine with two eight-core Intel Xeon E5-2690 CPUs (total 32 threads with Hyper-Threading) and one Intel Xeon PHI 3100 accelerator with 57 cores (total 228 threads -each core supports four hardware threads). The machine had 64Gb DDR3 RAM and the accelerator had 6Gb of on-board GDDR5 RAM. We run 32+228=250 threads to make computations and utilize 250 PRNG sequences accordingly. The proportion of work distribution was ≈ 1.6 : 1.