Rubella Vaccine Introduction in the South African Public Vaccination Schedule: Mathematical Modelling for Decision Making

Background: age structured mathematical models have been used to evaluate the impact of rubella-containing vaccine (RCV) introduction into existing measles vaccination programs in several countries. South Africa has a well-established measles vaccination program and is considering RCV introduction. This study aimed to provide a comparison of different scenarios and their relative costs within the context of congenital rubella syndrome (CRS) reduction or elimination. Methods: we used a previously published age-structured deterministic discrete time rubella transmission model. We obtained estimates of vaccine costs from the South African medicines price registry and the World Health Organization. We simulated RCV introduction and extracted estimates of rubella incidence, CRS incidence and effective reproductive number over 30 years. Results: compared to scenarios without mass campaigns, scenarios including mass campaigns resulted in more rapid elimination of rubella and congenital rubella syndrome (CRS). Routine vaccination at 12 months of age coupled with vaccination of nine-year-old children was associated with the lowest RCV cost per CRS case averted for a similar percentage CRS reduction. Conclusion: At 80% RCV coverage, all vaccine introduction scenarios would achieve rubella and CRS elimination in South Africa. Any RCV introduction strategy should consider a combination of routine vaccination in the primary immunization series and additional vaccination of older children.


Supplemental Material Supplement 1: Model Description
We developed a discrete-time stochastic age-structured compartmental rubella transmission model for South Africa, building from previous work describing rubella dynamics [16,20]. The key feature of the model is a matrix that at every time-step defines transitions from every combination of epidemiological stage (maternally immune 'M', susceptible 'S', infected 'I', recovered 'R', and vaccinated 'V', taken to indicate the effectively vaccinated) and age group (1 month age groups up to 20 years old, then 1 year age groups up to 100 years old; 321 total age groups) to every other possible combination of epidemiological stage and age group. The discrete time-step was set to about 16 days (i.e., 24 time steps in a year), the approximate generation time of rubella. We simulated a deterministic run for each of the vaccination scenarios from year 1995 to 2050. Figure S1 displays the epidemiological transitions of the transmission model. The model is agestructured so that each epidemiological transition is age-specific, and depending on the parameter also time-specific. Here, da is the probability of losing maternal immunity by age class a, is the probability an individual in age class a becomes infected, r is the recovery rate, and va,t is the probability an individual in age class a and time-step t is successfully vaccinated.

Epidemiological Parameters
The duration of protection by rubella maternal antibodies ranges between 3 and 9 months; accordingly, we modelled the probability of remaining in the maternally immune epidemiological stage over age (1-da) as an exponential decay function with a constant rate of 0.95 per month [28].
The probability of infection by age, (also called the age-specific force of infection, FOI) is a function of n(t), a vector describing the population at time t, defined as, ( ) = (M , , S , , I , , R , , V , , M , , … V , ) according to where z is the total number of age classes (here z = 321), , , is the rate of transmission between individuals in age class a and j at time-step t, also known as the Who-Acquires-Infection-From-Whom (WAIFW) matrix, and , is the number of infected individuals in age class j and time-step t, while captures the non-modeled heterogeneities in age mixing and the effects of discretization of the underlying continuous process. We fix at 0.97 reflecting values obtained for measles in England and Wales [27] , because discrete-time models that do not incorporate this exponent result in unrealistically unstable dynamics prone to frequent extinction. Given that rubella transmission is frequency dependent, we divide the number of infected individuals in each age class by the total population size at time-step t (∑ ( )).
Transmission to individuals in age group a from individuals in age group j for each time-step t is defined by , , = , , (1 + (2 )), where , is mean transmission from individuals in age group j to age group a, and is a parameter controlling the magnitude of seasonal fluctuations. Previous validation of this model has shown that model results for the burden of CRS were robust to the magnitude of seasonal fluctuations [16]; we set to 0.35 and held it constant over time [16]. Mean transmission from individuals in age class j to age class a, , , was estimated by rescaling populationadjusted age-contact rates (per POLYMOD study based on diary entries [29]) to reflect the assumed basic reproductive number (R0) of rubella. The value of R0 used in this analysis of 7.9 and was obtained from a previously published modelling study estimating R0 for 40 African countries [21]. We proceeded to run simulations with different estimates for R0 in a sensitivity analysis. The highest estimate used was an R0 of 12 which was estimated in Ethiopia [22] and the lowest estimate estimated in Burkina Faso was 3.3 [21].
The recovery rate, r, is equal to 1, such that by the next time-step (or rubella generation) all infected individual will immediately move into the recovered class.
The probability an individual in age class a and time-step t is successfully vaccinated, va,t, depends on the assumed vaccination coverage rate assumed over time and vaccine effectiveness over age. The vaccination coverage rate ranges from 0 to 1 and is vaccination scenario specific ( Table 1 in the main text). Vaccination effectiveness rate over the first 11 months of life was empirically estimated from data extracted from Boulianne et al. 1995 [29] forcing saturating at 97% and staying constant at 97% for all ages 12 months and older.

Demographic Parameters
Demographic parameters (population size, crude birth rates, and age-specific death rates) were country-specific and extracted from the United Nations World Population Prospects 2015 (cran package wpp2015).
The number of births per time-step t were estimated by multiplying the crude birth rate per time-step t (i.e., annual crude birth rate divided by 24 generations in a year) by the total population at time-step t (∑ ( )). Age-specific death rates as of 1995, extracted at five year age intervals, were estimated for all 321 age classes using smoothing splines and held constant over time. We assumed a constant rate of aging into the next age class (i.e., 1 divided by the length of age class a in years multiplied by 24).
To simulate rubella dynamics, we first needed country-specific rubella endemic populations (n(1)). We began with fully susceptible populations based on country-specific population and age structure estimated for 1995. The one year age interval population estimates were stratified into 321 age classes using smoothing splines. In order to move beyond the transient non-seasonal outbreaks to populations representing endemic rubella, we seeded infected individuals into the population and iteratively simulated rubella dynamics for four 20-year increments assuming constant births and deaths. At the end of each 20 year cycle, we rescaled the mean transmission ( , ) by the assumed R0 and the population by the 1995 population and age structure, and then simulated again, four times total. The result was a country-specific population representing endemic rubella in 1995 (n(1)). In 2015, we rescaled the population size (n(t)) based on population total estimates for the respective year to correct for small population size differences that accumulate over time in our model.

Model Outcomes of Interest
Our model outputs the number of individuals in each age class and epidemiological stage at every time-step, allowing us to directly extract the number of rubella cases (i.e., the number of individuals in the 'I' infected epidemiological class) per age and time-step.
Age-and time-specific CRS cases were estimated by multiplying the age-specific number of susceptible individuals and probability of becoming infected over 16 week period (based on model output from each vaccination scenario), the sex ratio of the population and age-specific fertility rate (extracted from the United Nations World Population Prospects 2015), and finally the probability of CRS following rubella infection during the first 16 weeks of pregnancy (estimated 0.65 [14]).
The effective reproduction number (RE) was estimated from the model output using the next generation method [40].
Model diagram of age-structured model. Figure S1. Relationship between data and the age-structured model. Solid lines ending in arrows indicate either data or elements inferred from data (i.e., R0, the appropriate structure of the WAIFW) that directly enter the model. Individuals in the maternal immunity (M), susceptible (S), infected (I), recovered (R) and vaccinated (V) compartments are represented with arrows representing movement between compartments: ɗ is the probability of losing maternal immunity, ϕ is the probability of becoming infected, ᴦ is the recovery rate and v is the probability of being vaccinated.

Supplement 2: CRS incidence over time for scenarios 2 to 6 compared to scenario 1 with extreme values of R0
Prior to RCV introduction, the incidence of CRS when R0 =3.3 was about three-fold that of R0=12. A lower value of R0 implies the rate of infection is lower. As a result, individuals are becoming exposed to the pathogen later in life. The risk of infection is therefore higher in adulthood compared to the case when R0 is higher and this leads to an older age distribution of infected individuals and therefore a higher CRS incidence. Figure S2. CRS incidence over time at 80% RCV coverage for scenarios 2 to 6 compared to scenario 1 with R0 values of 3.3 and 12. The lines for scenarios 3 and 4 overlap with that of scenario 5 so only this line is visible. The vertical dotted line represents the year of RCV introduction. Table S3. Number of CRS cases averted and corresponding number of undiscounted DALYs averted for each scenario involving RCV introduction (2 to 6) compared to scenario 1. Estimates are shown for a range of routine vaccine coverage levels (60% through 95%) and for three time horizons: 10, 20 and 30 years. Figure S4a. Change in RE over time for scenario 2 compared to scenario 1. While RE never drops to values below one for 60% vaccine coverage, it takes between 11 and 14 years for RE to drop below one with vaccine coverage levels of 65% to 95%. The slow drop in RE can be explained by the time it takes for successive vaccinated cohorts to age and achieve sufficient reduction in rubella incidence. . Figure S4b. Change in RE over time for scenario 3 compared to scenario 1. Following RCV introduction, RE immediately drops to values way below one due to the wide age range of vaccinated individuals during the initial mass campaign. There is then a progressive rise in RE corresponding to accumulation of susceptible individuals missed during routine vaccination and the initial SIA, with this rebound being less prominent with increasing routine vaccine coverage. RE eventually goes above one only for 60% RCV coverage. Figure S4c. Change in RE over time for scenario 4 compared to scenario 1. Following RCV introduction, RE immediately drops to values way below one due to the wide age range of vaccinated individuals during the initial mass campaign. There is then a progressive rise in RE corresponding to accumulation of susceptible individuals missed during routine vaccination and the initial SIA, with this rebound being less prominent with increasing routine vaccine coverage. Following the second mass campaign 5 years after RCV introduction, RE drops again but resumes an upward trend, eventually going above one only for 60% RCV coverage. Figure S4d. Change in RE over time for scenario 5 compared to scenario 1. Following RCV introduction, RE immediately drops to values way below one due to the wide age range of vaccinated individuals during the initial mass campaign. There is then a progressive rise in RE corresponding to accumulation of susceptible individuals missed during routine vaccination and the initial SIA, with this rebound being less prominent with increasing routine vaccine coverage. Following subsequent mass campaigns every 5 years, RE drops again but resumes an upward trend. In this scenario, RE never goes above one irrespective of RCV coverage. Figure S4e. Change in RE over time for scenario 6 compared to scenario 1. It takes between 4 and 6 years for RE to drop below one for all vaccine coverage levels with higher vaccine coverages associated with quicker decrease in RE. The slow drop in RE can be explained by the time it takes for successive vaccinated cohorts to age and achieve sufficient reduction in rubella incidence.