In this work, we derive a data-driven integrated chance constraints stochastic optimization approach for finding optimal vaccination strategies under uncertainty to prevent epidemics in a multi-community setting. This new approach involves a stochastic optimization model of disease spread for multi-community settings. The model uses readily available population demographic data together with uncertain disease spread and vaccine efficacy data. This work is motivated by COVID-19 which has continued to spread globally since early 2020 resulting in extremely high mortality, morbidity, and hospitalizations. We now have new variants with uncertain disease transmission characteristics (1), and as can be seen in Table 1, the disease is still active worldwide but no longer at the pandemic scale. Both South Africa and the U.S. have active cases. We believe that it is only a matter of time before another epidemic or pandemic occurs (2). Therefore, it is imperative to be prepared and continue developing new data-driven infectious disease spread models and methods that will aid decision making during future epidemics.

The negative impact of the pandemic worldwide is still felt today, as it has changed our daily lives, not to mention the socio-economic impacts. To fight the pandemic, both pharmaceutical and non-pharmaceutical interventions were implemented and came at high social and economic cost. These costs are estimated to exceed those associated with conventional recessions and modern wars by far (3). It cost more than $16 trillion by 2021 to fight the disease in the U.S. alone (3). In South Africa, for the 2020/21 fiscal year from April 2020 to March 2021, the cost was estimated at R142 billion (about $8.67 billion) of which R21.544 billion (about $1.32 billion) was spent directly on health care services (4). Pharmaceutical interventions included vaccinations and therapies, while non-pharmaceutical interventions included quarantines, lockdowns, and mask mandates.

Clinical trials are usually conducted to estimate the effectiveness of a vaccine or vaccine efficacy (VE), which is the percentage reduction in disease among the vaccinated. There are several ways vaccine efficacy and effectiveness are determined, e.g., vaccine efficacy against infection, severity of symptoms, hospitalization, or mortality. The reported vaccine efficacy values vary significantly, with some vaccines being relatively effective against certain virus variants (5). Vaccine efficacy can also be measured by how it modulates severity of symptoms, reduction in hospitalization as well as the length of hospitalization, and reduction in morbidity. Each of the virus variants has its own spread characteristics and how it infects and affects an individual based on age in terms of their resulting susceptibility and infectivity. Susceptibility is the ability to become infected if in close contact with an infected individual, whereas infectivity is the ability of an infected individual to spread the disease to others. A close contact is when a susceptible individual is within the proximity of an infected individual and may possibly transmit the disease. Vaccine efficacy and an individual's susceptibility, infectivity, and vaccine efficacy are all uncertain. Consequently, this poses challenges in modeling optimal vaccine allocation strategies toward preventing and controlling infectious diseases. This paper makes advances by incorporating uncertainty in close contact within and outside the household, susceptibility, infectivity, and vaccine efficacy into the optimization model.

Non-pharmaceutical interventions during the peak of the pandemic included mask mandates, social distancing, gathering restrictions, quarantining, contact tracing, lockdowns, and border closures. Beginning in early 2020, several countries developed monitoring and alert systems. For example, South Africa enacted five intervention levels, Alert Levels 1 to 5 (6). The lowest alert level (Level 1) allows for most normal activities with precautions and health guidelines followed at all times, while the highest alert level (Level 5) imposes drastic measures such as prohibiting mass gatherings, to contain the spread of the virus and reduce disease burden. Non-pharmaceutical interventions affect the spread of an infectious disease outside the household, which is uncertain at best and is measured in terms of close contact rate.

An important concept in epidemiology to determine whether an epidemic will continue or not is the basic reproduction number R₀, which is the average number of secondary infections produced by a typical case of an infection in a population (7, 8). This means that during an epidemic, R₀>1. However, what is actually measured during an epidemic at a given time t is the effective reproduction number, R(t) (9, 10). During a vaccination campaign, the goal is to have the vaccine induce herd immunity in the population, and in this case, the post-vaccination reproduction number of infected households R_HV, is measured. To prevent epidemics R_HV ≤ 1, and the objective is to optimize vaccination coverage, i.e., to determine the minimum proportion of the population to vaccinate (11). The challenge toward achieving this lies in characterizing R_HV since it is a random variable and involves several parameters (e.g., vaccine efficacy, susceptibility, and infectivity) all of which are uncertain.

We define an optimal vaccination strategy as the fraction (or percentage) of the population in each age group to vaccinate using the minimum number of vaccines to prevent epidemics at a specified reliability level for a given community. Therefore, we take a mathematical programming approach that integrates a disease spread model based on (11) and stochastic programming (12–14). Specifically, we devise an integrated chance constraints stochastic programming (ICC-SP) model that extends a previous model based on chance constraints (15). The idea of using ICC-SP in the epidemiology setting was first introduced by (16) to enable determining data-driven vaccination strategies under a stochastic post vaccination reproduction number. This work continues along the same vein by: (a) incorporating intervention levels into the model as binary decisions; (b) incorporating three types of vaccine efficacy, which do not include mortality: vaccination effectiveness against infection (VEI), vaccination effectiveness against symptomatic cases (VES), and vaccination effectiveness against hospitalization (VEH); and (c) applying the model to historical data for the South Africa setting.

We give a summary of closely related work on optimal vaccination strategies in Table 2. The table gives the citation, the model type (deterministic, simulation, and stochastic), and whether the model includes the following factors: households, age-specific heterogeneity in disease susceptibility and infectivity, vaccine efficacy type (VEI, VES, and VEH) used, and intervention levels. The symbol ✓ in the table means that the model includes the factor, while the symbol X means that the model does not include that factor. The main contributions of this work to the literature on vaccine allocation models include a new data-driven ICC-SP model for determining optimal vaccination strategies for multi-community settings, a computational study illustrating the application of the model using COVID-19 data for the South Africa multi-community setting, and the insights drawn from the study. The ICC-SP model provides an advance over previous related models by incorporating VEI, VES, and VEH, and intervention levels. Public health decision-makers need such data-driven models to aid in determining optimal vaccination strategies to prevent epidemics. The data-driven ICC-SP model is versatile in the sense that it uses readily available census demographic data, available disease spread data, and social interaction data for the communities under study. Thus, the ICC-SP model can be applied to any multi-community setting where the aforementioned data are available. Besides determining optimal vaccination strategies for each community, the model has a unique feature, namely selecting the optimal intervention level to prevent epidemics. This is especially important for public health authorities that need data-driven models to decide the appropriate level of intervention for each community during epidemics.

In this section, we describe the materials and methods used in our study. We start with the description of the model data in the next subsection and then present our data-driven multi-community stochastic model of disease spread in the subsequent subsection. We apply our methodology in a case study that involves census demographic and COVID-19 data for the South Africa setting encompassing five contiguous municipalities in Gauteng province as shown in Figure 1: City of Johannesburg, Ekurhuleni, City of Tshwane, Sedibeng (comprising Emfuleni, Midvaal, and Lesedi), and West Rand (comprising Merafong City, Rand West City, and Mogale City). Since the first case of COVID-19 was reported in South Africa in March 2020, the country experienced several waves of COVID-19 infections, first due to the delta variant, then in November 2021 South Africa was one of the first countries to report the presence of the Omicron variant (17). However, data on the effective reproduction number for each municipality during this period are not available (18).

We implemented the ICC-SP model to investigate the effectiveness of COVID-19 vaccination strategies under vaccination effectiveness against infection (VEI), hospitalization (VEH), and symptomatic (VES) cases using data for Gauteng province described in Section 2.1. The population data for the municipalities used in this study is as follows: Ekurhuleni (4,024,285), Johannesburg (5,079,469), Sedibeng (1,125,281), Tshwane (3,832,516), and West Rand (1,007,756) (20). The ICC-SP model was coded using C++ in Microsoft Visual Basic in the IBM ILOG CPLEX 12.9 Callable Library environment (40). We chose this software environment to enable implementation of a cutting plane decomposition method based on (35) to solve instances of ICC-SP. We created several instances of the model based on a set of predetermined reliability levels α_c to generate optimal vaccination strategies and intervention levels for the five municipalities. All the experiments were conducted on a computer workstation with 4.2 GHz CPU with 16 Cores and 64 GB of RAM. Due to the large-scale nature of the ICC-SP because of having multiple communities, household types, scenarios, and binary decision variables, we used the cutting-plane decomposition method to solve the instances. For the provincial level instances that we created, the average computation (CPU) time to run each ICC-SP problem instance is around 300 seconds with a wall clock time of about one hour. Thus, the ICC-SP model is suitable for up to the provincial level real-time (hourly) optimal decision making and can be used for real-time or near-real-time policy (minutes) requiring good feasible solutions. Scalability to larger regions or national level models would require scenario reduction techniques and other approximations, and further improvements in the decomposition solution approach.

We experimented with the three levels for reliability defined in Section 2.1, namely Low, Medium, and High as illustrated in Figure 1. The reliability levels allowed for assessing the ICC-SP model's optimal solutions across the historical range of values of R_t (effective reproduction number) above one (i.e., z_ωℓc) for each municipality under no vaccines at the lowest intervention ℓ (Level 1). The range of values of z_ωℓc for each municipality across all scenarios after solving the ICC-SP model are as follows: Ekurhuleni [0, 4.16], Johannesburg [0, 4.00], Sedibeng [0, 4.14], Tshwane [0, 4.12], and West Rand [0, 4.16]. The acceptable expected excess α_c value for each municipality c and reliability level was calculated via experimentation as a percentage of 𝔼[z(ω~)]. High reliability was assigned the smallest value, followed by Medium reliability, and then Low reliability as shown in the table. We experimented with limited vaccine availability of up to two-thirds of the total population of the five municipalities. To consider age-related differences in infectivity of COVID-19, we investigated two cases in our experiments. In the first case, we experimented with age-related differences in infectivity such that Group D has lower relative infectivity than Groups C, B, and A, in that order (younger population is more infective). In the second case, we considered Group A with lower relative infectivity than Groups B, C, and D, in that order (older population is more infective). We present detailed results for the first case and then discuss the key findings for both cases.

We first ran the ICC-SP model without fixing intervention levels to determine both the optimal proportion of the population to vaccinate to prevent epidemics and the optimal intervention level ℓ for each community and reliability level α_c. The results are summarized in Table 5. Under VEI as the vaccine efficacy criterion and Low reliability (i.e., relatively high risk in terms of preventing epidemics) the optimal intervention level to prevent epidemics across all municipalities is Level 1 (lowest). However, the optimal proportion of the population to vaccinate is largest for Ekurhuleni and Tshwane (round 61%), followed by Johannesburg and Sedibeng (around 58%), and then West Rand (around 45%). Under Medium reliability (i.e., relatively medium risk in terms of preventing epidemics) the optimal intervention levels to prevent epidemics from highest to lowest are as follows: Tshwane at Level 3 (round 63%), Ekurhuleni at Level 3 (round 62%), Johannesburg at Level 2 (around 72%), Sedibeng at Level 1 (around 78%), and West Rand at Level 1 (around 58%). Under High reliability (i.e., relatively low risk in terms of preventing epidemics) the optimal intervention levels to prevent epidemics from highest to lowest are as follows: Ekurhuleni and Tshwane both at Level 4 (round 54%), Johannesburg at Level 3 (around 76%), Sedibeng at Level 2 (around 93%), and West Rand at Level 1 (around 71%). In terms of general trends across reliability levels, we can see that both Ekurhuleni and Tshwane require the largest proportion to vaccine as well as high intervention levels under both Medium and High reliability levels, followed by Johannesburg, Sedibeng, and then West Rand. We also observe that the proportion to vaccinate under both Medium and High reliability exceed 70% of the population for some of the municipalities. Such a high proportion to vaccinate may not be attainable in practice and would be an indication that the epidemics may not be prevented.

We performed a preliminary empirical analysis toward validating the ICC-SP model against observed outcomes by linking model outputs to the observed effective reproduction number R_t trends. We used the COVID-19 estimated daily R_t trends data for Gauteng province, South Africa, that are available at (24). The data starts from March of 2020 around the onset of the COVID-19 pandemic. In Figure 1, we plot the expected post-vaccination reproduction number R_HV obtained from the ICC-SP model against the observed daily R_t trends in Gauteng. As can be seen from the graph, the ICC-SP model is able to bring down the R_t values very close to one, an indication that the vaccination strategies from the model are able to reduce the spread of the disease as desired. We provide a detailed view of the ICC-SP model expected R_HV values for each of the five municipalities in Gauteng at the bottom of Figure 1. The minimum and maximum R_HV values obtained from the ICC-SP model for each municipality are as follows: Ekurhuleni (1.09, 1.12), Tshwane (1.09, 1.12), Johannesburg (1.04, 1.07), Sedibeng (1.02, 1.07), and West Rand (1.01, 1.03).

When VES is used as the vaccine efficacy criterion, we observe similar trends in terms of having elevated intervention levels and relatively large proportions to vaccinate under Medium and High reliability levels for Ekurhuleni, Johannesburg, and Tshwane. Sedibeng and West Rand both have the lowest intervention level and the proportion to vaccinate; albeit the proportion to vaccine for Sedibeng is around 70% under Medium reliability and just over 80% under High reliability. These percentages may not to be feasible in reality. Under Low reliability, however, all municipalities are at intervention Level 1 and require around 65% or more to vaccinate except for West Rand, which require around 50%. These proportions are more likely to be attainable in practice. Interestingly, when VEH is used as the vaccine efficacy criterion, all municipalities remain at intervention Level 1 and in this case, all the municipalities but one (West Rand) require vaccinating between 60% and 70% of the population under all the three reliability levels. West Rand requires vaccinating between 50% and 60% of its population.

In practice, the same intervention level is usually applied across multiple communities to effectively contain the epidemics. Therefore, we experimented with fixing the intervention level in the ICC-SP model across the five municipalities to determine the optimal proportion to vaccinate to prevent epidemics at a given intervention level. The results for the three reliability levels for VEI, VES, and VEH are plotted in Figures 3–5, respectively. As can be seen in the plots in Figure 3 with VEI as the vaccine efficacy criterion, all five intervention levels are possible under Low reliability level. However, under Medium reliability only intervention Levels 3, 4, and 5 are possible. Under High reliability only Levels 4 and 5 are possible. We observe that, under Medium reliability, the epidemics cannot be controlled with intervention Levels 1 and 2. Similarly, under High reliability the epidemics cannot be prevented with intervention Levels 1, 2, and 3. Regarding the proportion of the population to vaccinate in each municipality to prevent epidemics, we observe from the plots that the proportion decreases with increasing intervention level as expected. We also notice that Ekurhuleni and Tshwane require the largest proportion to vaccinate under all three vaccine efficacy criteria (VEI, VES, and VEH) followed by Johannesburg and Sedibeng. West Rand requires the least proportion to vaccinate across all intervention levels. In terms of reliability levels from Low to High for a fixed intervention level, we can see that larger proportions have to be vaccinated to prevent epidemics. For example, Ekurhuleni under intervention Level 4 requires vaccinating 26.6% of the population to prevent epidemics under Low reliability, 39.4% under Medium reliability, and 53.5% High reliability.

Our ICC-SP model is a heterogeneous model in terms of age and therefore, besides determining the optimal total proportion to vaccine to prevent epidemics in each municipality, the model also computes the corresponding proportion that must be vaccinated in each age group for each household size. Recall that the age groups are A (age ≤ 19years), B (20 ≤ age ≤ 39years), C (40 ≤ age ≤ 64years), and D (age≥65years). In addition, in this case study we assume age-related differences in infectivity such that Group D (older age group) has lower relative infectivity than Groups C, B, and A (younger age group), in that order. We report the computational results using bar graphs in order to visualize the trends in the proportion to vaccinate in each age group. Figures 6, 7 show the plots for VEI as the vaccine efficacy criterion under Medium and High reliability levels, respectively. The graphs provide useful information for public health decision makers, and it is interesting to see the variations in the trends between municipalities at all levels of intervention from Level 1 to 5.

As can be seen in the figures, the proportion to vaccinate for each age group increases with household size in general at each intervention level, with Group D being prominent only in larger households. Overall, we observe that age groups B and C are prioritized with age groups A and D appearing in larger households. Clearly, larger household sizes require vaccinating more of the population in all the age groups while household sizes 1 and 2 (in some cases size 3 as well) do not require vaccinating anyone. We can observe that groups C, B, and A are prioritized in that order starting from smaller to larger households, while Group D is not prioritized in smaller households. As expected, we also see that higher intervention levels require vaccinating relatively smaller proportions of the population to prevent epidemics. Notice that the older age group D is the least to vaccinate at intervention Level 5 in all municipalities. In terms of reliability levels, the results show that the proportion to vaccinate increases with reliability level. Consider, for example, Ekurhuleni under intervention Level 4. To prevent epidemics, all age groups must be vaccinated at 100% only for households of size 10 under Low reliability; sizes 9 and 10 under Medium reliability; and sizes 8, 9, and 10 under High reliability. We also obtained age-related proportions to vaccinate to prevent epidemics using VES and VEH as vaccine efficacy criteria. The results show similar trends to VEI in general, albeit the proportions to vaccinate are relatively lower for VES and VEH, in that order. Groups C and B are prioritized, followed by A and then D. When we consider the case where Group A has lower relative infectivity than Groups B, C, and D, in that order (older population is more infective), the results show that Groups C and B are prioritized, followed by D and then A.

This computational study reveals that of the three vaccine efficacy criteria, VEI is the strictest criterion to determine optimal vaccination strategies, followed by VES and then VEH. VEI considers actual infection of the disease, whether symptomatic or asymptomatic, and typically affects most of the population during an epidemic and is measured by testing. However, VES involves severe symptoms and may only affect a subset (e.g., certain age groups) of the population. Similarly, VEH involves hospitalization and tends to affect a subset of the population of symptomatic cases. However, both VES and VEH can have a significant impact on optimal vaccination strategies, especially for virus variants that cause severe symptoms that lead to hospitalization. In general, the results clearly show that a higher intervention level generally requires vaccinating a relatively smaller proportion of the population to control epidemics. We observe that under High reliability epidemics can only be controlled by imposing a relatively high intervention level in order to achieve vaccination strategies that are practically feasible. This observation confirms the need for intervention control measures coupled with vaccinations to bring epidemics under control.

The results show that vaccination allocations are different based on the population and demographics of each municipality. In particular, the results reveal that Tshwane and Ekurhuleni require the highest proportion to vaccinate under the three vaccine efficacy criteria (VEI, VES, and VEH) followed by Johannesburg, Sedibeng, and West Rand, in that order. West Rand requires the least proportion to vaccinate across all intervention levels. Recall that among the five municipalities, Johannesburg has the largest population of 5,079,469 with a population density of 3,088/km². It is followed by Ekurhuleni with a population of 4,024,285 and a population density of 2,038/km². Next is Tshwane with a population of 3,832,516 and population density of 609/km² followed by Sedibeng with a population of 1,125,282 and a population density of 270/km². West Rand has the smallest population of 1,007,757 and a population density of 247/km². Apart from each municipality's overall population and population density, more granular factors that influence the optimal vaccination strategy include household characteristics such as age heterogeneity and size, other household demographics, geospatial features, and socioeconomic factors.

Prioritizing who to vaccinate in terms of age varies with relative infectivity and susceptibility of each age group. When the younger population is relatively more infective than the older population, Groups C and B are prioritized, followed by D and then A. There is a general trend to prioritize vaccinating the younger and mid-age group population (groups B and C) in smaller families, and the older age group (Group D) in larger households. This is true across the five municipalities, and the reason is that the younger population has higher relative infectivity compared to the older population. When the older population is more infective than the younger population, Groups B and C are prioritized, followed by D and then A. This is due to the older population having higher relative infectivity compared to the younger population. In essence, the combined relative infectivity and susceptibility of the age groups appear to be relatively high for mid-age groups, therefore prioritizing vaccinating age groups B and C to prevent epidemics.

In terms of ethical and equity considerations, prioritizing the younger population (assuming having high relative infectivity) can be effective in reducing and/or preventing epidemics, but can come at the expense of adverse effects (e.g., severity of symptoms, hospitalization, mortality) for the older population, especially if they have high relative susceptibility. On the contrary, prioritizing the older population may minimize those adverse effects, but may result in not containing the disease outbreaks due to the younger population spreading the disease. The ICC-SP model could be extended to include equity constraints in two keyways. First, one can set the relative infectivity and susceptibility parameters in the model for the different age groups appropriately to reflect the needed equity considerations. Second, vaccine constraints can be incorporated into the model to reflect how many vaccines are allocated to each age group to reflect the equity considerations at hand.

Other household demographics that can be critical include family and marital status, gender distribution, and income. How best to deal with a single-family multi-person household from one family may differ significantly from dealing with a multi-family multi-person household. Geospatial features and socioeconomic factors also come into play. Due to the geospatial legacy of South Africa's policies, certain interventions in high-density, low-income township areas may be less effective than in low-density, high-income surburban areas. It may be necessary to vary the nature of the types of interventions not only between municipalities but also within municipalities.

We derive a data-driven integrated chance constraints stochastic programming (ICC-SP) approach for finding optimal vaccination strategies under uncertainty to prevent epidemics in a multi-community setting. This new approach involves a stochastic optimization model of disease spread that includes multi-community population demographic data, uncertainty disease spread data (susceptibility and infectivity), and uncertain vaccine efficacy data. The model determines the optimal vaccination strategies, i.e., that is the proportion of individuals in a given community and household type to vaccinate to bring the reproduction number below one under a given decision-maker's reliability (or risk) level. This data-driven approach was implemented and tested based on COVID-19 data for five neighboring municipalities in Gauteng province, South Africa. The results reveal, among other insights, that to prevent or control epidemics, vaccination strategies should be designed to prioritize vaccinating specific households and age groups with high levels of combined relative susceptibility and infectivity. The results also show that intervention levels play a key role in the containment of epidemics. Specifically, high levels of intervention reduce the proportion of the population that has to be vaccinated. A key advantage of the ICC-SP model is that it is data-driven and its decisions adapt to the underlying stochastic data fed to the model. In addition, the model enables the decision maker to consider different levels of risk, allowing for what-if analyses for public health policy. We believe that to be ready for future epidemics, there is still a strong need for data-driven optimization models to determine optimal vaccination strategies that adapt to new virus variants and uncertain vaccine efficacy. Future work along this line of research includes incorporating vaccine hesitancy and vaccine distribution logistics into the ICC-SP model.