key: cord-0018871-4kn6d408
authors: Lee, Yunjeong; Lee, Dong Han; Kwon, Hee-Dae; Kim, Changsoo; Lee, Jeehyun
title: Estimation of the reproduction number of influenza A(H1N1)pdm09 in South Korea using heterogeneous models
date: 2021-07-07
journal: BMC Infect Dis
DOI: 10.1186/s12879-021-06121-8
sha: a0fd9195818d70040aee58a3633de1ac9b38ff18
doc_id: 18871
cord_uid: 4kn6d408

BACKGROUND: The reproduction number is one of the most crucial parameters in determining disease dynamics, providing a summary measure of the transmission potential. However, estimating this value is particularly challenging owing to the characteristics of epidemic data, including non-reproducibility and incompleteness. METHODS: In this study, we propose mathematical models with different population structures; each of these models can produce data on the number of cases of the influenza A(H1N1)pdm09 epidemic in South Korea. These structured models incorporating the heterogeneity of age and region are used to estimate the reproduction numbers at various terminal times. Subsequently, the age- and region-specific reproduction numbers are also computed to analyze the differences illustrated in the incidence data. RESULTS: Incorporation of the age-structure or region-structure allows for robust estimation of parameters, while the basic SIR model provides estimated values beyond the reasonable range with severe fluctuation. The estimated duration of infectious period using age-structured model is around 3.8 and the reproduction number was estimated to be 1.6. The estimated duration of infectious period using region-structured model is around 2.1 and the reproduction number was estimated to be 1.4. The estimated age- and region-specific reproduction numbers are consistent with cumulative incidence for corresponding groups. CONCLUSIONS: Numerical results reveal that the introduction of heterogeneity into the population to represent the general characteristics of dynamics is essential for the robust estimation of parameters.

The reproduction number is defined as the average number of secondary cases generated by a typical primary case. It is a measure of the transmission potential associated with the contact rate, duration of infectivity, and probability of transmission per contact. The maximum reproduction number is attained when an infectious person is introduced into a totally susceptible population et al. demonstrated the manner in which sanitary measures reduce the prevalence of an infected population [1] . Estimates of the reproduction number were shown to decrease from 1.4-1.5 initially to 1.1-1.2 later in the summer, which was most likely because of the vacation period and the seasonality of influenza transmission [5] . In addition to capturing temporal dynamics, it is important to consider heterogeneous patterns of the transmission. It is well known that school-age children are disproportionately responsible for influenza transmissions. Estimates of the age-specific reproduction number help with our understanding of the role of each group in the transmission dynamics and with devising effective targeting mitigation strategies.

If all incident cases could be traced back to their index cases, estimating the reproduction number would simply be a matter of counting the number of secondary cases. However, with most epidemics, only the epidemic curve is observed, and there is no available information regarding who infected whom. To appropriately estimate the reproduction number from the influenza outbreak data, it is essential that the selected model capture the underlying dynamics embedded in the data. The objective of this study is to estimate the reproduction numbers based on the incidence data.

When the World Health Organization announced the emergence of influenza A(H1N1)pdm09 (pH1N1) in 2009 [6] , the first probable patient in South Korea was identified on April 28. A total of 763,759 confirmed cases, of which 270 were fatal, were reported by the end of August 2010 [7] . During the initial epidemic phase, the main control measure was containment through quarantine and isolation. Surveillance programs in schools and medical facilities were implemented, and all confirmed cases were investigated. However, when community outbreaks were detected in June, the intervention policy switched from containment to mitigation, including vaccination and antiviral prescription. Vaccination was started on October 27, 2009, and 12.7 million people were vaccinated by the end of August 2010. Before August 20, antiviral agents were prescribed to patients with acute febrile respiratory illness (AFRI) and who had a history of travel abroad or contact with a confirmed patient. However, when the number of community-acquired cases increased, antiviral agents were prescribed to patients with AFRI symptoms.

According to the database, 3,087,788 courses of antiviral agents were prescribed from August 21, 2009 to April 30, 2010. The daily number of incident patients was estimated based on the amount of prescribed antiviral agents (Refer to [8] for details). Figure 1 shows the temporal incidence distribution of pH1N1 in South Korea. The amount of antiviral agents prescribed and the number of incident patients soared from mid-October, reached its peak at the end of October, and started declining in mid-November. Demographic and regional characteristics are illustrated in Table 1 and Fig. 2 . The incidence rate is higher in children and students than in other age group individuals ( Table 1 ). The rate is higher in urban areas than in rural areas and is the highest in the national capital and the south-eastern region (Fig. 2) .

In this paper, two different structured models are proposed to estimate reproduction numbers on the basis of the epidemic curve. We begin by introducing a basic SIR model to describe a single outbreak and build ageand region-structured models by incorporating population heterogeneity. Numerical simulations are conducted to analyze the impact of terminal time and the effect of heterogeneous structures on the estimation of parameters. Finally, the proposed models are applied to the 2009 incidence data of novel pH1N1 in South Korea to compute the age-and region-specific reproduction numbers.

Cumulative incidence by region, which is own figure drawn from the incidence data

We consider the standard SIR model to represent single-outbreak influenza dynamics. The model classifies individuals into three key compartments: susceptible, infected, and recovered. The nonlinear system of differential equations describing the dynamics is given by the following equation.

The state variables S(t), I(t), and R(t) denote the number of individuals who are susceptible, infected, and recovered, respectively, at time t. The number of contact events sufficient for transmitting an infection during unit time per individual is βN based on the mass action incidence. Infective individuals leave the compartment at the recovery rate γ , thereby acquiring immunity to the disease. We can drop the equation for R(t) as it has no effect on the dynamics of others and hence is determined once S(t) and I(t) are known. Based on the SIR model, we have the timedependent net reproduction number R t = βS(t)/γ , which quantifies the level of transmission at time t. Note that R t is the per-infective rate at which new infections occur within the average duration of infection at time t.

An age-structured model is employed to estimate the reproduction number of pH1N1 because the transmission rate is higher in preschool and schoolchildren than in other age group individuals, in general. We consider a subgroup SIR model where the population is divided into n a age groups with different transmission dynamics. We denote the number of susceptible and infected individuals within the i th age group by S i and I i , respectively. Let β ij refer to the transmission from the j th age group to the i th age group, and β=[ β ij ] denote the transmission matrix, also known as Who-Acquires-Infection-From-Whom matrix. Putting these elements together, we have the following system of differential equations.

In a general structured model of the form (1) with n a distinct classes, n 2 a transmission terms are required. However, one transmission term is available at most for each class. The typical way to address this lack of specificity is to constrain the structure of the transmission matrix and/or to use prior knowledge of social mixing behavior. For an age-structured model, we assume that the transmission rates are proportional to the rates of social contact, which can be estimated from contact patterns. A large multi-country population-based survey conducted in Europe as a part of the POLYMOD [9] enables us to implement this approach. The transmission is modeled as the product of the contact rate in the survey and an age-specific proportionality factor to account for characteristics related to susceptibility and infectiousness, which are not captured by contact rates. This leads to

where c ij is the contact rate and q i and σ are proportionality factors. Based on the age-structured SIR model (1), the reproduction number can be calculated by following Driessche and Watmough [10]. It is the spectral radius of the next generation matrix M where

The details are given in Appendix.

The second mechanism incorporates a heterogeneous population based on regions to account for the wave of the pH1N1 pandemic. We denote the number of susceptible and infected individuals within the i th region by S i and I i , respectively. Let β ij refer to the transmission from the j th subgroup to the i th subgroup and β=[ β ij ] denote the transmission matrix. In the same manner as the age-structured model, we have the following system of differential equations.

We assume that transmission rates between distinct regions in the region-structured model can be expressed as the frequency of transportations multiplied by a regionspecific proportionality factor. The transportation information was extracted from the highway portal site and Kakao map for number of buses and highway traffic, respectively [11, 12] . Let the number of buses and highway traffic from region j to region i be denoted by w i,j and W i,j , respectively. Note that w is symmetric because the bus route is circular, although W is not necessarily. The transmission rate can be written as

where q i is the proportionality factor, and σ l are σ g can be chosen such that they balance the weight between different types of transportations.

The same argument as that presented in the agestructured model gives the expression of the effective reproduction number R t , which is the spectral radius of the following next generation matrix

Study subjects were patients who were prescribed antiviral agents from the national stockpile from August 21, 2009 to April 30, 2010. Because of mandatory antiviral agent management program during study period, all patients who were prescribed antiviral agent were included in this study. The data employed to estimate the parameters are the daily number of incident patients in Fig. 1 . It was estimated based on the aggregation of prescribed antiviral agents from deidentified database. This study was approved by the Institutional Review Board (IRB) of Yonsei University Health System. Since this study used retrospective data and the study subjects were anonymized, the IRB waived the requirement for written consent from the patients. Our goal is to estimate the optimal parameters that provide the states that are best fit to the given data. This section briefly reviews the parameter estimation technique of the least squares method. In general, parameter estimation is conducted by minimizing the cost function, which measures the difference between the model prediction and observation. The simplex algorithm proposed by John Nelder and Roger Mead is applied to solve the optimization problem [13] . Let θ be the parameter set and time points t j (j = 1, ..., N) are uniformly distributed with daily time step. The data vector y j (j = 1, ..., N) denotes the number of cases at time t j . It is a scalar for basic SIR model, but it is a vector structured by age and region for age-stratified and region-stratified models, respectively. For example, the y j is a column vector of length 15 for the age-structured model. We recast the mathematical model as

and assume a statistical model for measurement of the form

where f t j ; θ is the model prediction at time t j with parameter θ and the measurement error ε j ∼ N 0, a 2 . The least squares estimator can be obtained by minimizing the following cost function over the given parameter space θ [14] :

The parameter sets to be estimated for the basic SIR model, age-structured model and region-structured model, are θ basic = {β, γ } and initial values of S and I, θ age = q i , σ , γ i fori = 1, 2, · · · , n a and initial values of S i and I i , θ region = q i , σ l , σ g , γ i fori = 1, 2, · · · , n r and initial values of S i and I i where n a and n r denote the number of age groups and regions, respectively.

We illustrate the proposed methodology and investigate its performance by applying it to 2009 incidence data of pH1N1 in South Korea. This section presents the results of the estimation obtained by applying the least squares method to basic SIR, age-and region-structured models. In each experiment, data with different time periods by varying the terminal time are tested to determine the earliest stage of the epidemic sufficient to provide a reasonable estimation. Figure 3 displays the predicted incidence based on the basic SIR model compared with the observed data. Predictions using data only during the initial growth phase cannot effectively exhibit the dynamics and substantially overestimate the spread of the infection. The results of simulation improved after the peak of the epidemic, and the wave is roughly generated at a later stage.

However, the simple SIR model does not provide a reasonable estimation of parameters. The estimated values of γ and R 0 demonstrate a large variation and remains outside of the feasible range for the influenza, regardless of the time period for data in Fig. 4 . The plausible reason for this involves the model assumptions that are too simple to capture the underlying mechanisms.

For the age-structured SIR model (1), the total population is split into 14 subgroups of 5-year age bands and one with 70 years and older (i.e., 0-4, 5-9, ..., 65-69, 70 + ). This incorporates a heterogeneous population into the model in order to reflect different transmission rates in each age group. In Fig. 5 , the outbreaks are simulated using data during various time periods in the same manner as mentioned above. Results of both the models show similar trends as long as the terminal time is earlier than mid-November when the gentle growth begins during the decline stage. Additionally, as the growth begins to decline, the age-structured SIR model fits the incidence data better than the basic SIR model. In Fig. 4 , the reproduction number starts increasing in early October, peaks at 2.5 on October 17, and then decreases to unity at the end of October. Real-time estimation demonstrated that the effective reproduction number rose sharply during mid-October when the number of patients increased dramatically. The reproduction number fell below unity at the end of October and stayed lower than unity indicating that the epidemic starts decreasing, which is consistent with the incidence data. In the age-stratified model, heterogeneity was incorporated by WAIFW matrix where the transmission was assumed to be proportional to the contacts. The effective contacts were measured by POLY-MOD contact survey, which showed a clear evidence for an age-dependency in contact patterns. Taking heterogeneous mixing into the model enabled better description of the dynamics, because the trend in behavior was consistent with the demographic characteristics of cases (as 

shown in Table 1 ). Estimated parameters are possible indicators to determine the feasibility of models. Incorporation of the age structure allows for robust estimation of parameters, while the basic SIR model provides estimated values beyond the reasonable range with severe fluctuation in Fig. 4 . Table 2 summarizes the parameter estimates using three different models. The estimated duration of infectious period using age-structured model is around 3.8. The reproduction number was estimated to be 1.6 which is similar to those obtained in Mexico, the United States, New Zealand, Peru, and Chile [2, [15] [16] [17] [18] .

The general characteristics of regional difference led us to consider a second type of heterogeneity in the model. The nation is split into 252 in the region-structured model (2) , where the transmission rates are implemented based on transportation patterns. Figure 6 compares the predicted cases based on the region-structured SIR model with the observed data over the course of the epidemic. As it was discussed in the previous experiment, it is not earlier than the epidemic peak for estimation to start adjusting to outbreak data. Since this outbreak, the incidence data is well described in the form of the characteristic exponential rise, turnover, and decline pattern predicted by the process model. The estimated duration of infectious period using region-structured model is around 2.1 and the reproduction number was estimated to be 1.4 ( Table 2 ). The time-dependent effective reproduction number is also illustrated in Fig. 4 , which demonstrates a pattern similar to that obtained using the age-structured SIR model. Estimated duration of infectious period and reproduction numbers using three different models are compared in Fig. 4 . Values of the cost function defined by (3) are also provided in Fig. 7 , which shows the goodness-of-fit in the order of region-structured, age-structured and simple SIR model.

It is widely known that the transmission is considerably different among various age groups. We also observe from the pH1N1 epidemic data that the incidence rate is higher in children and students than in other age groups ( Table 1) . Estimates of the age-specific reproduction number help in clarifying the role of each age group in the transmission dynamics and in suggesting guidelines for effective targeting intervention strategies. The estimated age-specific reproduction numbers are displayed in Fig. 8 . The result is closely related to the cumulative incidence for each age group because it is often the contact rate within the same age group is higher than with other groups. The incidence rate is higher in urban areas than in rural areas, and the highest in the national capital and the south-eastern region, as shown in Fig. 2 . We estimated the region-specific reproduction number and observed that it is more than two in some areas and less than one in the others (Fig. 9 ). This is consistent with regions having larger cumulative incidence with a similar argument regarding contact patterns to age-specific cases.

An estimation of reproduction numbers is crucial because it provides a measure of the transmission potential when an infection is spreading throughout a population. The reproductive numbers in the early phase of Influenza A(H1N1)pdm09 have been estimated in several countries with different settings, yielding median 1.46 and range 1.0-3.6 [19] . Many of these studies focused on cases confirmed in the early stage of the pandemic. Because laboratory tests focused on severe cases and there are possible changes in laboratory testing and notification rates, the number of confirmed cases does not necessarily represent the underlying epidemic. It also does not reflect the dynamics during the period of the highest level of activity, which is the winter in temperate climates. Some studies used the number of cases from sentinel surveillance that is much less than the actual number of influenza patients. It is necessary to estimate the reproductive number using the number of all the patients throughout a pandemic, including the period with the highest level of activity. In this study, the reproductive number was estimated based on the national data of incidence deduced from antiviral agent prescription in South Korea during the pandemic.

We discussed parameter estimation methodologies based on deterministic SIR models that included age or spatial structures with the main aim being to estimate the effective reproduction numbers, R 0 and R t . There could be many modelling choices to compute these important epidemiological parameters, including simple SIR model with time varying parameters [20, 21] or renewal equations [22, 23] . We proposed one possible approach to introduce population heterogeneity since we observed demographic and regional characteristics of incidence data. Age-structured and region-structured models were employed to describe the underlying epidemic process, in particular. To avoid exacerbating non-identifiability problem by increasing the complexity of a model, age-and region-specific data were used to estimate parameters for age-and region-structured model, respectively. And Fig. 9 Region-specific reproduction number (left) and cumulative incidence for each region (right) estimation showed that the reproduction number started increasing in early October, peaked at 2.5 on October 17, and then decreased to unity at the end of October. The effective number rose sharply during the mid-October when the number of patients increased dramatically. The reproduction number fell below unity at the end of October and remained lower that unity, indicating that the epidemic starts decreasing, which is consistent with the incidence data. Subsequently, age-specific and region-specific basic reproduction numbers were estimated to account for the differences of incidence. We observe from the pH1N1 epidemic data that the incidence rate is higher in children and students than in other age groups. The estimated agespecific reproduction numbers agree with the cumulative incidence for each age group because the mixing is assortative. The incidence rate is higher in urban areas than in rural areas, highest in the national capital and in the south-eastern region. We estimated the region-specific reproduction number whose trend is similar to the number of cases in each region. Estimates of the age-specific and region-specific reproduction number help to predict the transmission dynamics, and to suggest guidelines for effective targeting intervention strategies.

This study has both limitations and strengths. First, the number of cases is estimated from the amount of prescribed antiviral agents assuming the time lag between symptom onset and antiviral agent prescription, the proportion of prescription and the proportion of pH1N1 confirmation among AFRI patients. Also, vaccination is not considered in the model. However, the effect of vaccination on the transmission of pH1N1 may have been insignificant because the vaccination for general group was initiated in January 2010. The effective contacts were employed from POLYMOD contact survey, which possibly yields discrepancy in mixing pattern of Korea [24] . The use of POLYMOD as well as the potential nonidentifiability of complex models are limitations of this study. We will be able to improve the outcome as we gather more information, because additional knowledge is required to achieve a better result.

On the contrary, the present research has its strengths compared to previous studies. The reproduction number was estimated based on the national level antiviral agent prescription data in South Korea throughout the pandemic including the period of the highest level of activity. The real-time estimation incorporating population structures can be used to predict the disease dynamics, thereby providing guidelines for the optimal implementation of preventive measures, such as school closing and distribution of antiviral agents.

Numerical results reveal that the introduction of heterogeneity into the population and sufficient data to represent general characteristics of dynamics are essential to the robust estimation of parameters. Real-time estimation showed that the reproduction number started increasing in early October, peaked on October 17, and then decreased to fell below unity at the end of October, which is consistent with the incidence data. The estimated age-and region-specific reproduction numbers are also consistent with cumulative incidence for corresponding groups.

The reproduction number for the age-structured SIR model (1), can be calculated following the approach of Driessche and Watmough [10]. Let F i be the new infections and V i be the transitions of i th compartment, then

and

for i = 1, · · · , n a . Subsequently, the derivatives of F =[ S n a σ c n a 1 S n a σ c n a 2 · · · S n a q n a c n a n a respectively, and the next generation matrix is

S 1 q 1 c 11 S 1 σ c 12 · · · S 1 σ c 1n a S 2 σ c 21 S 1 q 2 c 22 · · · S 2 σ c 2n a . . . . . . . . . . . .

S n a σ c n a 1 S n a σ c n a 2 · · · S n a q n a c n a n a ⎤ ⎥ ⎥ ⎥ ⎦ .

Thus, the reproduction number is the spectral radius of FV −1 and the age-specific reproduction number is the column sum of FV −1 corresponding to the age of interest. 

Ahued Ortega A, Villaseñor Ruíz I. Modelling of the influenza A(H1N1)v outbreak in Mexico City

Pandemic potential of a strain of influenza A (H1N1): early findings

Pros and cons of estimating the reproduction number from early epidemic growth rate of influenza A (H1N1) 2009

The transmissibility and control of pandemic influenza A (H1N1) virus

The effective reproduction number of pandemic influenza: prospective estimation

World Health Organization. Influenza-like illness in the United States and Mexico

Analysis of reported pandemic influenza (A/H1N1 2009) virus infections in Korea

Estimation of the reproduction number of pandemic influenza A(H1N1) 2009 from national data of antiviral agent prescription

):e74. 10. van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Highway public data portal

A simplex method for function minimization

Modeling and Inverse Problems in the Presence of Uncertainty

Transmissibility of 2009 pandemic influenza A (H1N1) in New Zealand: effective reproduction number and influence of age, ethnicity and importations

Epidemiological and transmissibility analysis of influenza a (h1n1) v in a southern hemisphere setting: Peru

The transmissibility and control of pandemic influenza A (H1N1) virus

Outbreak of 2009 pandemic influenza A (H1N1)

Estimates of the reproduction number for seasonal, pandemic, and zoonotic influenza: a systematic review of the literature

Accounting for non-stationarity in epidemiology by embedding time-varying parameters in stochastic models

Assessing the impact of non-pharmaceutical interventions on SARS-CoV-2 transmission in Switzerland

A new framework and software to estimate time-varying reproduction numbers during epidemics

Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures

A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Not applicable.Authors' contributions DHL and JL designed the work, and CK contributed to collect data. YL carried out data processing and numerical implementation, and JL, HDK, and YL interpreted the results. JL wrote the manuscript and DHL, CK, HDK and YL assisted with drafting the work. All authors read and approved the final manuscript.

The work of Jeehyun Lee was supported by NRF-2015R1A5A1009350 and HG20C0003030020. The work of Hee-Dae Kwon was supported by NRF-2016R1D1A1B04931897 and the NST grant by the Korean government (MSIP) (No. CRC-16-01-KRICT). The funding source had no role in the study process, including the design, collection, analysis, and interpretation of data and in writing the manuscript.

All data generated or analysed during this study are included in the published article [8] . The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

This study was approved by the Institutional Review Board (IRB) of Yonsei University Health System. The administrative permissions was not required to access the data used in this study and the data was anonymised before its use.

Not applicable.

The authors declare that they have no competing interests.Author details 1 Department of Computational Science and Engineering, Yonsei University, 50, Yonsei-ro, 03722 Seoul, South Korea. 2 Korea Disease Control and