key: cord-0898434-7pfzgrvm
authors: Matsuyama, Ryota; Akhmetzhanov, Andrei R.; Endo, Akira; Lee, Hyojung; Yamaguchi, Takayuki; Tsuzuki, Shinya; Nishiura, Hiroshi
title: Uncertainty and sensitivity analysis of the basic reproduction number of diphtheria: a case study of a Rohingya refugee camp in Bangladesh, November–December 2017
date: 2018-04-02
journal: PeerJ
DOI: 10.7717/peerj.4583
sha: 13f8fe9cd4b06843fce134699ecc8b89cf6a2353
doc_id: 898434
cord_uid: 7pfzgrvm

BACKGROUND: A Rohingya refugee camp in Cox’s Bazar, Bangladesh experienced a large-scale diphtheria epidemic in 2017. The background information of previously immune fraction among refugees cannot be explicitly estimated, and thus we conducted an uncertainty analysis of the basic reproduction number, R(0). METHODS: A renewal process model was devised to estimate the R(0) and ascertainment rate of cases, and loss of susceptible individuals was modeled as one minus the sum of initially immune fraction and the fraction naturally infected during the epidemic. To account for the uncertainty of initially immune fraction, we employed a Latin Hypercube sampling (LHS) method. RESULTS: R(0) ranged from 4.7 to 14.8 with the median estimate at 7.2. R(0) was positively correlated with ascertainment rates. Sensitivity analysis indicated that R(0) would become smaller with greater variance of the generation time. DISCUSSION: Estimated R(0) was broadly consistent with published estimate from endemic data, indicating that the vaccination coverage of 86% has to be satisfied to prevent the epidemic by means of mass vaccination. LHS was particularly useful in the setting of a refugee camp in which the background health status is poorly quantified.

Diphtheria, a bacterial disease caused by Corynebacterium diphtheriae, is a vaccinepreventable disease. Symptomatic patients initially complain sore throat and fever. Additionally, a gray or white patch causes the "croup," blocking the airway and causing a barking cough. Due to widespread use of diphtheria-tetanus-pertussis (DTP) vaccine globally, the incidence has steadily declined over time, and thus, diphtheria is commonly perceived as a disease of pre-vaccination era. Nevertheless, sporadic cases and even epidemics of the disease have been yet reported especially in politically unstable areas, and many cases have been considered as arising from susceptible pockets of the vulnerable population (Rusmil et al., 2015; Hosseinpoor et al., 2016; Sangal et al., 2017) .

In 2017, multiple diphtheria outbreaks were reported in refugee camps, including those in Yemen and Bangladesh (World Health Organization (WHO), 2017a). Of these, a Rohingya refugee camp in Bangladesh, which is temporarily located in Cox's Bazar, experienced a large-scale diphtheria epidemic. As of December 26, 2017, the cumulative number of 2,526 cases and 27 deaths were reported (World Health Organization (WHO), 2017a). To interrupt chains of transmission, emergency vaccination has been conducted among children since December 12, 2017, achieving the overall coverage greater than 90% by the end of 2017 (World Health Organization (WHO), 2018). Due to vaccination effort and other countermeasures, including contact tracing and hospital admission of cases, the epidemic has been brought under control, with incidence beginning to decline by the end of December 2017 (World Health Organization (WHO), 2017a).

Considering that diphtheria has become a rare disease in industrialized countries, epidemiological information on model parameters that govern the transmission dynamics has become very limited, and thus, it is valuable to assess how transmissible diphtheria would be through the analysis of the recent outbreak data. The basic reproduction number, R 0 , is interpreted as the average number of secondary cases that are produced by a single primary case in a fully susceptible population, acting as the critical measure of the transmissibility. To date, an explicit epidemiological estimate of R 0 for diphtheria has been reported only by Anderson & May (1982) : using a static modeling approach to age-dependent incidence data with an assumption of the endemic equilibrium, R 0 was estimated as 6.6 in Pennsylvania, 1910s and 6.4 in Virginia and New York from 1934 to 1947. Subsequently, a few additional modeling studies of diphtheria took place (Kolibo & Romaniuk, 2001; Sornbundit, Triampo & Modchang, 2017; Torrea, Ortega & Torera, 2017) , but none of these offered an empirical estimate of R 0 .

Here we analyze the epidemiological dataset of diphtheria in Rohingya refugee camp, 2017, aiming to estimate R 0 in this particular epidemic setting. Given that the epidemic occurred among refugees, we explicitly account for uncertainties associated with unknown background information, including the fraction of previously immune individuals and ascertainment rate of cases.

The latest epidemic curve was extracted from the report of the World Health Organization (WHO) Regional Office for South East Asia (SEARO) (World Health Organization (WHO), 2017a). Figure 1 shows the latest available epidemic curve. As of December 26, 2017 (the latest date of observation), a total of 2,526 cases have been reported. Cases consist of (i) confirmed cases: cases reported as positive for C. diphtheriae by multiplex assay, (ii) probable cases: cases with upper respiratory tract illness with laryngitis or nasopharyngitis or tonsillitis AND sore throat or difficulty swallowing and an adherent membrane (pseudomembrane) OR gross cervical lymphadenopathy, and (iii) suspected cases: any case with a clinical suspicion of diphtheria, including cases that are unclassified due to missing values (World Health Organization (WHO), 2017b). Prior to December 11, 2017, cases satisfied only the condition of suspected cases. The definition was improved on and after December 12 (World Health Organization (WHO), 2017c), enforcing to count mainly probable cases, while not restricting the reporting of suspected cases. For this reason, cases reported by December 11 are considered to have been likely over-ascertained compared with cases that were reported later under improved case definition. Mass child vaccination started on December 12, and according to the administrative coverage, greater than 90% vaccination coverage was achieved by December 30 (World Health Organization (WHO), 2018). Vaccine-induced immunity requires at least 7-14 days to become effective, and moreover, the reporting delay was assumed to be about four days (based on additional analyses; results not shown). For these reasons, the dataset from December 23-26 was discarded as the number of cases may be influenced by emergency vaccination and also biased by the reporting delay.

To formulate the epidemiological model, here we mathematically capture the epidemiological process of secondary transmission using R 0 and the serial interval, i.e., the time from illness onset in the primary case to illness onset in the secondary case. We assume that secondary transmission does not take place before illness onset. According to a classical study by Stocks (1930) in the United Kingdom (UK), the time interval from first to second diphtheria cases in the household revealed a bimodal shape. Following Klinkenberg & Nishiura (2011) , the first peak corresponds to an independent infection in the community and the second peak reflects within-household transmission. As one of peak time-lags in the observed time interval was observed on day 8, we assumed that the mean serial interval was eight days, and we imposed an assumption that the coefficient of variation (CV) of the serial interval distribution was 50%, and later varied it from 25% to 75% as part of the sensitivity analysis.

To capture the epidemiological phenomena of reproduction, it has been shown that the renewal process can capture the serial stochastic dependence structure (Nishiura, 2010) . Let i t be the number of new cases on day t. g t represents the distribution of the serial interval. To describe the time-dependent incidence i t on day t, we have

(1)

where s t represents the fraction of susceptible individuals on day t. The renewal process of this type is not original to the present study and has been applied to other settings including real-time epidemic modeling studies (Asai & Nishiura, 2017; Dinh et al., 2016; Ejima & Nishiura, 2018; Endo & Nishiura, 2015; Nishiura et al., 2010 Nishiura et al., , 2016 Tsuzuki et al., 2017) . It should be noted that the incidence i t includes both symptomatic and asymptomatic cases. Let c t be the reported number of cases on calendar day t. Supposing that only the fraction a t among the total number of infections are diagnosed and reported, c t satisfies

where a t is modeled as a function of t. Because the case definition was improved from December 12, 2017 onward, the ascertainment rate likely varied around that time. Namely, we set a t = a 1 for time by December 11 and a 2 on December 12 and later. We assumed that only the ascertainment rate changed as a function of time, and also that R 0 and depletion of susceptible individuals were unaffected by time.

We model the fraction susceptible s t on day t in the following way. Let n represent the previously immunized fraction so that only fraction (1 -n) of the population is susceptible at the beginning of the epidemic. In addition to the previously immune fraction, s t decreases when natural infection takes place. Suppose that the total population size was N, s t is written as

We assume that N is equal to the population size of epidemic area within Rohingya refugee camp as 579,384 persons (Banerji & Ahmed, 2017) . Accordingly, the renewal equation is written as

where t in the right-hand side indicates the time since infection (or the so-called "infection-age"). We assume that c t follows a Poisson distribution. The likelihood to estimate u consisting of the parameters R 0 , n, and a t is derived as

where T is the latest time of observation (i.e., December 22 in our case study) and c T = (c 1 , c 2 , : : : , c T ).

While we specified unknown parameters as R 0 , n, and a t , it is expected that R 0 is correlated with initially immune fraction n and also a t . Thus, it is vital to quantify R 0 while accounting for the uncertainty of other model parameters. In the present study, n was estimated through uncertainty analysis, while a t was estimated as a step function governed by two parameters, a 1 and a 2 . Uncertainty in parameter values can be addressed by randomly sampling the uncertain parameter value from probability distributions (Gilbert et al., 2014) . Here we use the Latin Hypercube sampling (LHS) method (Sanchez & Blower, 1997) in which a symmetric triangular distribution of n was assumed to be in the range from 0.0 to 0.7; the health survey of Rohingya population indicated that overall 30.8% of children had received no vaccinations (Guzek et al., 2017) , and the remaining 69.2% receives at least a single immunization, which is not necessarily against diphtheria, thus, vaccination coverage for diphtheria should be less than 69.2%. Thus, we expect that the actual coverage is nearby the mid-point of the range from 0.0 to 0.7. In addition, we examined the sensitivity of R 0 to variations in the length of the serial interval.

The present study analyzed data that is publicly available. As such, the datasets used in our study were de-identified and fully anonymized in advance, and the analysis of publicly available data without identity information does not require ethical approval. Figure 2 shows univariate distributions of estimated parameters R 0 and a t based on LHS (n = 1,000) of parameter n from 0 to 0.70 with the peak at 0.35. a 1 Reflects ascertainment by December 11, 2017, while a 2 shows the same on and after December 12. R 0 took the minimum and maximum estimates at 4.7 and 14.8, respectively, with the median estimate at 7.2. The distribution was skewed to the right with the mode at 6.7. Excluding lower and upper tails, 950 samples (95%) of R 0 (or what it can be assumed as the 95% tolerance intervals) were in the range of 5.0-12.3. Distributions of a 1 and a 2 were also right skewed. a 1 Ranged from 0.007 to 0.021 with the median 0.010, while a 2 ranged from 0.003 to 0.011 with the median 0.005. Lower and upper 95% tolerance intervals of a 1 and a 2 were (0.007, 0.017) and (0.004, 0.009), respectively. A decline of a t was observed, because the incidence was supposed to have grown more during the early phase (i.e. by December 11) given the estimated epidemiological parameters. In practical sense, the timing coincided with the involvement of contact tracing practice enforced by the Médecins Sans Frontières. Figure 3 shows the distributions of two estimated parameters in two-dimensional spaces and also the comparison between observed and predicted epidemic curve. As can be expected from Eq. (4), R 0 and n were positively correlated given an epidemic curve. Specifically, as n increases R 0 must also increase to achieve an identical epidemic curve. Similarly, R 0 was positively correlated with ascertainment rates, a 1 and a 2 . While the ascertainment might have been lowered due to the stricter case definition, the small estimate of a 2 can also indicate that a substantial fraction of undiagnosed individuals existed and the susceptible fraction was then gradually depleted in the population. The observed and predicted epidemic curves are compared in the right lower panel of Fig. 3 . While the model is kept simple with four unknown parameters, the predicted epidemic curve overall captured the observed pattern.

Sensitivity of R 0 to the variation in the serial interval is shown in Fig. 4 . We varied the CV of the serial interval distribution from 25% to 75%. When the CV was 25%, the median and mode of R 0 from LHS were 9.4 and 8.2, respectively. When the CV was 75%, the median and mode of R 0 were estimated to be 5.7 and 5.2, respectively. That is, as the variance increases, the estimate of R 0 decreases. However, the variation in R 0 induced by the CV was smaller than the uncertainty associated with the initially immune fraction. 

The present study estimated R 0 of diphtheria in the Rohingya refugee camp, explicitly accounting for case ascertainment and previously immune fraction. Since previously immune fraction n of the refugee population was not precisely known, uncertainty analysis of R 0 was conducted with an input parameter assumption for n employing the LHS method. R 0 ranged from 4.7 to 14.8 with the median estimate at 7.2. To our knowledge, the present study is the first to statistically estimate R 0 of diphtheria from epidemic data. For the statistical estimation of R 0 , the renewal process model was employed, which has an advantage to handle the right-censored data during the course of an epidemic, compared with other available methods for completely observed data, e.g., Nishiura (2010) . Estimated median R 0 was broadly consistent with the value ranging from 6 to 7 as indicated by Anderson & May (1982) based on a static model for endemic data that uses the age-dependent incidence in the UK. We have shown that the frequently quoted estimate agrees well with dynamically estimated R 0 from the refugee camp in the present day. Assuming R 0 = 7, to control diphtheria by means of mass vaccination, the coverage greater than 86% must be satisfied. Since our study focused on uncertainty and sensitivity analyses, the exact estimate of R 0 cannot be pointed out. However, despite the uncertainty regarding n in this population, we estimated a distribution of R 0 consistent with previous estimates. While the mode of distribution for R 0 was 6.7, the validity of this value depends on the validity of our prior distribution of n, which was not supported by any published evidence of this refugee population. Nevertheless, a Demographic and Health Survey data of the Rohingya population in Myanmar indicated a close value from 40% to 50% as the coverage of DTP (Ministry of Health and Sports, 2017). Ascertainment rates were jointly estimated only by using the epidemiological case data and the population size. What we have shown in the present study is that when we have an access to not only the initial growth rate of the epidemic but also the incidence data around the time at which peak incidence is observed, R 0 and susceptible fraction can potentially be jointly quantified. Even without explicit estimate of the initially immune fraction, we have shown that an indication of the possible value of R 0 can be obtained through uncertainty analysis. LHS appeared to be particularly useful in the setting of refugee camp in which the background health status is not well quantified (Helton & Davis, 2002; Nishiura et al., 2017) . LHS can offer probabilistic distribution of the outcome measure, R 0 in our case, and this method appeared to be particularly useful when one or more uncertain input information exist (Elderd, Dukic & Dwyer, 2006; Coelho, Codeço & Struchiner, 2008; Samsuzzoha, Singh & Lucy, 2013; Gilbert et al., 2014) . While Bayesian modeling has replaced LHS to some extent of uncertainty analysis as it can also offer posterior distributions of even uncertain parameters (Elderd, Dukic & Dwyer, 2006 ; Coelho, Codeço & Struchiner, 2008) , there could be an issue of identifiability when two or more parameters are evidently correlated, e.g., as anticipated between R 0 and n in our model (4). In such an instance, we cannot be sure if the limited epidemic data with the Bayesian estimation method can offer identifiable distributions for all parameters, and then LHS can remain to act as a useful tool for uncertainty analysis. The estimated small ascertainment rate during the epidemic is considered as reflecting the time-dependent diagnostic practice at a local level. In fact, the small ascertainment rate can lead us to observing the low case fatality risk of diphtheria in the Cox's Bazar that has been estimated to be small (27/2,526 = 1%), although the right censoring would of course matter for the real-time interpretation Mizumoto et al., 2015) . Nevertheless, cases could have died in the community unnoticed with low specificity of the case definition, and this is endorsed as a possible reason for observing the small number of deaths by the World Health Organization (WHO) (2017d). A follow-up study in this regard should be conducted in the future.

Several limitations must be noted. First, our model rested on a homogeneous mixing assumption. No heterogeneous patterns of transmission, including contact patterns and age-dependency were taken into account due to shortage of information. In addition, the time-dependent heterogeneity, including the reduced transmission potential due to contact tracing and rapid hospitalization, was not explicitly taken into account due to insufficiency of the data. If there were any additional indications or datasets that would allow explicit quantification of the effective reproduction number from December 12, that could give additional insights into the success of control measures. Second, for similar reasons, no spatial information was explicitly incorporated into the model. Third, a little more realistic features of refugee population, such as the impact of migration on the epidemic were unfortunately discarded in the present study. Similarly, one could investigate how overcrowding and malnutrition in the deprived population would help enhance the spread of diphtheria, given sufficient data backup from epidemiological investigations.

While these features need to be explicitly quantified in the future, we believe that our study adds an important piece of evidence to the literature on diphtheria. The transmissibility of diphtheria in the refugee population was estimated to be consistent with that in an endemic setting and mass vaccination must satisfy at least the coverage of 86% to halt the major epidemic of diphtheria.

The present study estimated R 0 of diphtheria in the Rohingya refugee camp, explicitly accounting for case ascertainment and previously immune fraction. Since previously immune fraction n of the refugee population was not precisely known, uncertainty analysis of R 0 was conducted with an input parameter assumption for n employing the LHS. R 0 ranged from 4.7 to 14.8 with the median estimate at 7.2. LHS can offer probabilistic distribution of the outcome measure, and this method appeared to be particularly useful in the setting of refugee camp in which the background health status is poorly quantified.

low-income and middle-income countries: a multicountry study of household health surveys

The correlation between infectivity and incubation period of measles, estimated from households with two cases

Mathematical model of the infection process in diphtheria for determining the therapeutic dose of antitoxic anti-diphtheria serum

Nay Pyi Taw: Ministry of Health and SportsRockville: The DHS Program

Estimating the risk of Middle East respiratory syndrome (MERS) death during the course of the outbreak in the Republic of Korea

Correcting the actual reproduction number: A simple method to estimate R 0 from early epidemic growth data

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

Transmission potential of Zika virus infection in the South Pacific

Early epidemiological assessment of the virulence of emerging infectious diseases: a case study of an influenza pandemic

Transmission dynamics of cholera in Yemen, 2017: a real time forecasting

The immunogenicity, safety, and consistency of an Indonesia combined DTP-HB-Hib vaccine in expanded program on immunization schedule

Uncertainty and sensitivity analysis of the basic reproduction number of a vaccinated epidemic model of influenza

Uncertainty and sensitivity analysis of the basic reproductive rate, tuberculosis as an example

Resurgence of diphtheria in North Kerala, India, 2016: laboratory supported case-based surveillance outcomes

Mathematical modeling of diphtheria transmission in Thailand

Infectiousness and immunity in regard to chickenpox, whooping-cough, diphtheria, scarlet fever and measles

A modeling of a diphtheria epidemic in the refugees camps

Dynamics of the pneumonic plague epidemic in Madagascar

Diphtheria outbreak response update. Cox's Bazar, Bangladesh-update #2. Geneva: World Health Organization

Geneva: World Health Organization

World Health Organization (WHO). 2017d. Mortality and morbidity weekly bulletin

Weekly Situation Report #9. Bangladesh: World Health Organization

The present study was financially supported by the Japanese Society for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research (KAKENHI) numbers for Hiroshi Nishiura: 16KT0130, 16K15356, and 17H04701, Ryota Matsuyama with 16H06581, and Shinya Tsuzuki with 17H06487. Hiroshi Nishiura received funding from the Japan Agency for Medical Research and Development, JSPS Program for Advancing Strategic International Networks to Accelerate the Circulation of Talented Researchers, The Telecommunications Advancement Foundation, Inamori Foundation, and the Japan Science and Technology Agency (JST) Core Research for Evolutional Science and Technology (CREST) program (JPMJCR1413). Akira Endo was supported by The Nakajima Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.Guzek J, Siddiqui R, White K, van 

The following grant information was disclosed by the authors: Japanese Society for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research ( 

Hiroshi Nishiura is an Academic Editor for PeerJ.

Ryota Matsuyama conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft. Andrei R. Akhmetzhanov performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, authored or reviewed drafts of the paper, approved the final draft. Akira Endo performed the experiments, analyzed the data, authored or reviewed drafts of the paper, approved the final draft. Hyojung Lee performed the experiments, analyzed the data, contributed reagents/ materials/analysis tools, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft.Takayuki Yamaguchi performed the experiments, authored or reviewed drafts of the paper, approved the final draft. Shinya Tsuzuki performed the experiments, authored or reviewed drafts of the paper, approved the final draft. Hiroshi Nishiura conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft.

Supplemental information for this article can be found online at http://dx.doi.org/ 10.7717/peerj.4583#supplemental-information.