Characterizing measles transmission in India: a dynamic modeling study using verbal autopsy data.

BACKGROUND: Decreasing trends in measles mortality have been reported in recent years. However, such estimates of measles mortality have depended heavily on assumed regional measles case fatality risks (CFRs) and made little use of mortality data from low- and middle-income countries in general and India, the country with the highest measles burden globally, in particular.
METHODS: We constructed a dynamic model of measles transmission in India with parameters that were empirically inferred using spectral analysis from a time series of measles mortality extracted from the Million Death Study, an ongoing longitudinal study recording deaths across 2.4 million Indian households and attributing causes of death using verbal autopsy. The model was then used to estimate the measles CFR, the number of measles deaths, and the impact of vaccination in 2000-2015 among under-five children in India and in the states of Bihar and Uttar Pradesh (UP), two states with large populations and the highest numbers of measles deaths in India.
RESULTS: We obtained the following estimated CFRs among under-five children for the year 2005: 0.63% (95% confidence interval (CI): 0.40-1.00%) for India as a whole, 0.62% (0.38-1.00%) for Bihar, and 1.19% (0.80-1.75%) for UP. During 2000-2015, we estimated that 607,000 (95% CI: 383,000-958,000) under-five deaths attributed to measles occurred in India as a whole. If no routine vaccination or supplemental immunization activities had occurred from 2000 to 2015, an additional 1.6 (1.0-2.6) million deaths for under-five children would have occurred across India.
CONCLUSIONS: We developed a data- and model-driven estimation of the historical measles dynamics, CFR, and vaccination impact in India, extracting the periodicity of epidemics using spectral and coherence analysis, which allowed us to infer key parameters driving measles transmission dynamics and mortality.

Background

The World Health Assembly officially endorsed in 2005 a goal to reduce measles-related deaths by 90% globally between 2000 and 2010 [1]. The Measles & Rubella Initiative (www.measlesrubellainitiative.org), a consortium of leading global health agencies launched in 2001, has been supporting the World Health Organization (WHO) strategies required to reduce measles mortality [2], including routine immunization of children, and supplementing it by a second dose opportunity for measles vaccine [3, 4]. In countries with routine immunization systems that have achieved good population coverage, the second dose of measles vaccine is usually included in the routine vaccination schedule [4]. Conversely, in countries where routine immunization services have not yet met these targets, the second dose of measles vaccine tends to be delivered through supplemental immunization activities (SIAs) [4].

Researchers at the WHO have reported decreasing trends in measles mortality [5]. Simons and colleagues [6] estimated that global measles-related deaths decreased from 535,000 (95% confidence interval (CI): 347,000–976,000) in 2000 to 139,000 (71,000–448,000) in 2010. Measles mortality was estimated to have reduced by more than 75% in all regions but the WHO Southeast Asia region [6]. In this region, India accounted for about 50% of global measles mortality in 2010 [6]. Yet, the WHO estimates depended heavily on measles case fatality risks (CFRs) reviewed from the literature spanning more than two decades [7]. Such reviewed CFRs were drawn from a variety of outbreak situations, surveillance methods, reporting conditions, and locations (e.g., community vs. in hospitals), leading to great heterogeneity and vast confidence intervals, and were subsequently pooled into regional estimates [7] for use in WHO’s mortality estimation [5, 6]. These estimates made little use of measles mortality data from low- and middle-income countries in general, particularly regarding India, the country with the highest measles burden globally. For instance, available measles data for India include not only national case notifications but also mortality data such as those from the Million Death Study [8-10]. The Million Death Study is an ongoing longitudinal study recording deaths across 2.4 million Indian households (14 million individuals) and attributing causes of death using verbal autopsy. Notably, the observed periodicity in measles deaths, similar to the periodicity in measles cases, can help us identify some of the intrinsic features of measles transmission in India. The periodicity in measles transmission may be difficult to extract because of multiple frequencies and noise contained in the measles mortality time series. However, the use of mathematical techniques such as the Fourier transform [11], which we implement in this study, is ideal for extracting such periodic features and therefore represents a novel aspect of the methodological approach we develop here and subsequently apply to India.

In this paper, we constructed a dynamic model of measles transmission in India that was empirically inferred from a time series of measles mortality extracted from the Million Death Study. The model was then used to estimate among under-five children: (1) the measles CFR, (2) the number of measles-related deaths, and (3) the impact that routine and SIA immunization has had since 2000. Model development and estimation was also pursued for Bihar and Uttar Pradesh (UP), two states with large populations, high under-five mortality rates [12], relatively low vaccine coverage [13], and the highest numbers of measles-related deaths in the country [10].

Methods

Measles deaths data

We used measles deaths data from the Million Death Study (MDS), an ongoing longitudinal study recording deaths occurring in 2.4 million households across India and attributing causes of death using verbal autopsy [8-10]. Deaths attributed to measles were available over a period of 3 years (2000–2003), and we analyzed weekly occurrence of measles deaths over this 3-year period in India and in the states of Bihar and UP since the day of the first measles death recorded (Fig. 1). We used anonymous secondary data (Additional file 1: Web appendix I provides MDS research ethics details), and the full protocol of the MDS can be consulted online (http://www.cghr.org/project.htm).

Fig. 1

Weekly measles deaths recorded over the period 2000–2003, for India (a), Bihar (b), and Uttar Pradesh (c), since the day of the first measles death (source Million Death Study), with their associated frequency spectra with strength of harmonic on y-axis and period up to 2 years on x-axis. The black lines correspond to the raw data and corresponding spectra, while the red line denotes the data and spectra post-detrending with the Baxter-King filter

Dynamic model of measles transmission

We adapted the Dynamic Measles Immunization Calculation Engine (DynaMICE), an existing dynamic compartmental discrete-time model of measles transmission and vaccination stratified by age previously used to examine measles transmission in high burden countries [14] (see Additional file 1: Web appendix II). As discussed in [14], the model population occupied states representing being susceptible to measles, being infected with measles, and having recovered from measles with lifelong immunity thereafter. The rate at which infection occurred depended on the existing proportion of the population who were already infected, the effective contact rate between different age groups, and a periodic forcing term (using a sine function) representing seasonal changes in measles transmissibility. Effective contact rates and mixing between individuals in different age groups were based on the British arm of the POLYMOD social contact survey [15], as with our previous work [14]. We also simulated measles case importation into the population at a rate of one per week, to mitigate against unpredictable effects on infection dynamics due to tiny fractions of infected and numerical solver accuracy.

Although individuals in the population were categorized according to discrete age classes from 0 to 80 years, a proportion of each age class had their age incremented by one at every time step. In keeping with the model described in great detail in [14], vaccinated individuals were assumed to have a reduced risk of measles infection. Vaccine efficacy was assumed to be 85% for the first dose when vaccinating before age one, 95% after age one, and 98% for two doses, as suggested by a meta-analysis [16]. Vaccines were assumed to be “all or nothing”: individuals receiving the vaccine were either completely protected or not at all. We assumed that vaccination gave lifetime protection if it successfully elicited an immune response, and that vaccinating already infected individuals did not increase the rate of infection clearance. The model was programmed using the software R version 3.2.2 (www.r-project.org).

As in our previous study [14], the population age distribution was obtained from government sources [17] along with crude birth rates and death rates to inform initial age-specific population sizes and mortality rates. These demographic parameters were assumed to be static over time (constant population growth) to avoid difficulties around long-term population projection. The average infectious period of measles was assumed to be 14 days [18]. Routine measles vaccination was assumed to be delivered before age one, with the SIA being given to all children between one and ten years. The probability of receiving a routine and SIA dose was assumed to be uncorrelated. Routine measles vaccine coverage was obtained from government sources [13]. The timing of the rollout and coverage estimates of the SIA doses were based on that achieved by the 2010–2013 measles SIAs [19, 20]. The SIA started in 2010 and increased in coverage year on year over three phases to reach a cumulative coverage estimate of 91% across states receiving the SIA. Different states had different timings; for example, Bihar received the SIA earlier (phases 2 and 3) than UP (phase 3, which started in 2012) [20].

Model calibration

Model dynamics depend primarily on parameters (amplitude of seasonal variation in measles a 0, basic reproduction number R 0) which cannot be directly observed, so we instead infer these parameters by fitting our model to the measles deaths data. First, we used the measles mortality time series (Fig. 1) to infer the periodicity in the oscillations in the number of measles deaths over time in India, Bihar, and UP. To do so, we first detrended the data using a Baxter-King filter [21] and then conducted a spectral analysis of each of the three time series applying a Fourier transform [11]. Subsequently, we obtained a spectrum composed of specific frequencies (or periods) and corresponding amplitudes. In our analysis, we used the five longest periods (1.50, 1.00, 0.75, 0.60, and 0.50 years).

Second, for India, Bihar, and UP, the dynamic model was run for 100 years to allow it to reach a pre-vaccination equilibrium, then run for a further 100 years with routine vaccination to reach a post-vaccination equilibrium, as implemented in [14]. We assumed routine first dose coverage only (at 56%), as there were no SIAs in the years for which measles mortality data were available (2000–2003). Subsequently, we obtained a periodic time series of measles cases for distinct values of the basic reproduction number (R 0) and the amplitude of the forcing term (a 0). R 0, which was defined using the dominant eigenvalue of the next generation matrix, was varied between 10 and 25 (based on ranges given in [18]), while a 0 was varied between 0 and 0.5 (Additional file 1: Web appendix III, Figure S1). R 0 was set by adjusting the probability of an infected individual transmitting measles to a susceptible individual following an effective contact. The spectrum of the time series of measles cases (once post-vaccination equilibrium was reached) was then extracted for each combination of values of R 0 and a 0. Assuming that periodicity in measles mortality would be the same as periodicity in measles incidence (i.e., a constant CFR), we then compared the model spectrum with the mortality data spectra for India, Bihar, and UP.

We used the coherence function [22] as a measure of similarity between the two time series (Additional file 1: Web appendix III, section 2). Each (R 0; a 0) combination was simulated and run until it reached equilibrium. For each period of interest from the mortality data spectrum, we estimated the coherence between the time series estimated from the (R 0; a 0) combination and the measles mortality time series. By combining the five coherence estimates for each period (1.50, 1.00, 0.75, 0.60, and 0.50 years), we produced a pooled (or combined) coherence [23] for each (R 0; a 0) combination (see the explanation in Additional file 1: Web appendix III), from which we could subsequently calculate a CFR estimate. We also validated our procedure to ensure that a 3-year time series provided sufficient information to recover values of R 0 and a 0 from the underlying infection process (see Additional file 1: Web appendix III, section 2).

Model analysis

First, based on the simulated models obtained, we estimated the total number of measles cases, per given age group, for India, Bihar, and UP for each R 0 and seasonal variation a 0 combination. Subsequently, to estimate CFR, we compared the total numbers of measles cases estimated among under-five children with the numbers of measles deaths among under-five children as extrapolated from the same measles deaths data at the national and state levels for 2005 by Morris and colleagues [10]. If the pooled coherence was significant, we used the associated coherence level for each simulated (R 0; a 0) combination as weight to sample from the estimated cases (i.e., case estimates from an (R 0; a 0) combination that had a coherence value of 0.8 would have twice the probability of being sampled than those with a coherence of 0.4). This method produced a weighted distribution of cases that we could then combine with the sampled distributions of estimated deaths during the same time period [10] (Additional file 1: Web appendix III, Table S1). This yielded a CFR distribution for measles among under-five children for India as a whole, Bihar, and UP, for 2005, provided with 95% CIs embedded in the coherence estimation procedure (Additional file 1: Web appendix III, section 2). Second, using the R 0 and a 0 parameters that produced the greatest pooled coherence across the five chosen periods and the estimated CFRs, we derived the number of measles deaths which would have occurred among under-five children over the period 2000–2015 in India, Bihar, and UP. Finally, we estimated the number of measles deaths for under-five children that would occur over 2000–2015 had (1) the SIAs not happened, and (2) routine vaccination not occurred.

Results

The measles frequency spectrum for India showed strong components for periods of 0.5, 1.0, and 1.5 years. The Bihar frequency spectrum had a principal component of 1.0 year, while the UP spectrum had strong components at periods of 0.5 and 1.0 year (Fig. 1). Subsequently, we ran the coherence estimations. Based on the estimations for the pooled coherence (Fig. 2), we extracted the following combinations for R 0 and a 0: R 0 = 24 and a 0 = 0.18 for India as a whole, R 0 = 20 and a 0 = 0.10 for Bihar, and R 0 = 14 and a 0 = 0.16 for UP.

Fig. 2

Pooled coherence estimate across five periods (0.5, 0.6, 0.75, 1.0, 1.5 years) for India (a), Bihar (b), and Uttar Pradesh (c), plotted against estimated basic reproduction number (R 0) and amplitude of the forcing term (Amp or a 0). Pink indicates a better match to periodicity in the mortality data. The location of the maximum coherence is identified with an asterisk (*), and the 95% significance contour (coherence = 0.11) is marked in dark blue

Furthermore, we plotted the pooled coherence against CFR for India, Bihar, and UP (Fig. 3), to display the associated CFR values to the coherence values. This is done because the varying of R 0 and a 0 sets both the periodicity and the number of cases in the model, which determine the coherence and CFR, respectively. With the pooled coherence estimates, we obtained the following estimates of CFRs among under-five children calibrated to the year 2005: 0.63% (95% CI: 0.40–1.00%) for India as a whole; 0.62% (0.38–1.00%) for Bihar; 1.19% (0.80–1.75%) for UP.

Fig. 3

Pooled coherence across five periods (0.5, 0.6, 0.75, 1.0, and 1.5 years) plotted against estimated case fatality risk (CFR) calculated from corresponding R 0 and a 0 for India (a), Bihar (b), and Uttar Pradesh (c). The left and right green lines delimit the range of the CFR where coherence is above the 95% significance level. The middle green line indicates the CFR for which the highest pooled coherence was estimated. The bottom red line indicates the upper 95% significance limit for the coherence function

Over 2000–2015, assuming CFRs remained constant over time, we estimated that a total of 607,000 (95% CI: 383,000–958,000) under-five deaths attributed to measles occurred in India as a whole. Likewise, 100,000 (61,000–161,000) under-five deaths caused by measles occurred in Bihar, which represented about 17% of all Indian under-five measles deaths; 313,000 (211,000–461,000) under-five deaths caused by measles occurred in UP, which represented around 52% of all Indian under-five measles deaths.

If the SIAs of 2010–2013 had not occurred in India, 66,000 (95% CI: 41,000–108,000) excess under-five deaths would have occurred; among which 7000 (4100–11,300) would have occurred in Bihar and 40,000 (27,000–58,000) in UP. Lastly, if no routine vaccination coverage or SIA occurred from 2000 to 2015, an additional 1.6 million (95% CI: 1.0–2.6 million) under-five deaths would have occurred across India with 110,000 (70,000–169,000) in Bihar and 491,000 (331,000–724,000) in UP.

Discussion

We present in this paper a data- and model-driven estimation of the historical measles dynamics, CFR, and vaccination impact in India as a whole as well as in two key states (Bihar and UP) for measles mortality. We used measles mortality data from verbal autopsy studies from a nationally representative longitudinal cohort of 2.4 million households and extracted the periodicity of measles epidemics, which allowed us to infer key parameters driving the dynamics of measles transmission. To do this we used spectral and coherence analysis together with statistical inference using an age-stratified dynamic compartmental model of measles transmission. This represents a first-of-a-kind modeling approach for measles transmission in the developing world.

Spectral analysis of measles mortality data suggested that the time series consisted of the superposition of several cycles of measles epidemics with different inter-epidemic periods. We were able to estimate very high basic reproduction numbers (R 0 > 20) for the country as a whole and Bihar. The case of UP with its lower estimated basic reproduction number (R 0 = 14) nevertheless points to the great complexity and potentially significant heterogeneity of measles transmission and mixing contact patterns within the state. In addition, our modeling approach allowed us to estimate a measles CFR between 0.40% and 1.80% in India, which offers a defensible range of CFRs compared with those previously reported in the literature [7]. It also enabled us to quantify the burden of measles mortality in India, for which we derived estimates in the range of other studies [24, 25]. Lastly, our analysis confirmed the likely high impact [26-28] of measles SIAs in India which were rolled out for about 110 million Indian children aged 9 months to 10 years [20] and could have potentially averted 66,000 deaths [19]. This is also consistent with the significant child survival impact of SIA recently estimated for 25 sub-Saharan African countries [29].

Our modeling approach gives policymakers tools to estimate the impact of measles routine immunization and SIAs in India, Bihar, and UP. India is highly diverse in terms of geography, health systems, and local epidemiology. Due to a lack of robust mortality data at state and district levels, this analysis modeled an increase in overall vaccine coverage across India to the average level seen in the recent SIA, rather than the details of the SIA itself. As additional data become available, our model will be able to capture finer scale disease dynamics and more accurately assess intervention impact.

Our study presents a number of strengths and limitations. Our model is one of the few dynamic models [6, 30] of measles transmission calibrated to measles data from low- and middle-income countries; it is to the best of our knowledge the only such model calibrated to measles mortality data taking advantage of the spectral features (periodicity) of measles infection [31]. As such, it enables us to directly estimate the measles CFR for India, setting itself free from the vast uncertainty in measles CFRs pooled from previous studies [7]. In general, most attempts at estimating measles mortality and vaccine impact in low- and middle-income countries, and especially in India, have focused on assuming a measles CFR taken from the literature [5, 6, 27, 28, 32] and on good reporting of measles case notifications during outbreaks, which is unlikely. Hence, our paper proposes a new approach that could strengthen measles mortality estimates and also be applied to other settings where historical time series of measles cases or deaths are available. It also builds on a rich mortality dataset [8] rather than on under-reported case notifications [33]. Yet, it has a number of limitations. First, the time series used was short, with only 3 years of measles deaths data, which increased the uncertainty of the spectral estimation. Hence, the longest cycles in the spectrum estimated (with period 1.5 years) may actually be even longer (e.g., 2 years). This is supported by the fact that most measles epidemics observed in the real world have integer-valued periods. Other datasets may present a longer time period; however, they often have a much coarser nature (e.g., WHO measles case notification data [33]), such as annual data rather than weekly data as with the MDS. This can thus be equally restrictive in the periods calculable with Fourier analysis (Additional file 1: Web appendix III, section 4). Second, it is possible that some measles deaths were not recognized as such in the MDS study, and as a result our CFR estimation would correspond to a lower bound because of such under-ascertainment. However, the MDS dataset has been well validated and scrutinized to reduce biases and misclassifications in the assessment of under-five deaths [9, 10], so we would anticipate this under-ascertainment to be minimal. Third, small number issues prevented us from examining additional age groups (above age 5) and additional states besides the highly populated states of Bihar and UP. Fourth, our mortality data dated back from 2000 through 2003, which suggests that our extrapolation into the future should be interpreted with caution. For example, the CFR could well decrease over time as under-five mortality decreases and health services improve. Fifth, our modeling assumption of equilibrium behavior is a simplification, because when vaccine coverage and birth rates are changing, transitions in measles epidemic cycles can also occur [31, 34]. However, this assumption may not be highly inaccurate, because measles coverage in India has been relatively stable in 2000–2003, and the crude birth rate has been only gradually decreasing from 1990 to 2005 [35, 36].

Conclusions

India is a key country for measles elimination objectives, and between 2010 and 2013 India carried out measles SIA campaigns targeting states where access to health services is often limited. Our modeling approach enables us to project measles incidence and mortality into the future, allowing us to assess whether India could reach measles elimination goals through both the scale-up of routine immunization and the implementation of SIA mass campaigns. A similar approach (e.g., Fourier analysis) could be applied in any setting with a time series of measles deaths or notifications indicating the periodicity of outbreaks, regardless of under-reporting in the data. Hence, one could examine the progress of other measles priority countries toward achievement of measles control and elimination.