Design of Targeted Community-Based Resource Allocation in the Presence of Vaccine Hesitancy via a Data-Driven Compartmental Stochastic Optimization Model

Epidemiological models are key to understanding and predicting the course of infectious disease outbreaks. With the emergence of COVID-19, the field of mathematical epidemiology quickly attracted broader interest. Taking advantage of existing approaches, researchers adapted and applied various methods specifically to study COVID-19, including compartmental models, structured metapopulations, agent-based networks, deep learning, and complex networks (Wu et al., 2020; Adiga et al., 2020; Rodríguez et al., 2021; Zhou et al., 2020a; Chang et al., 2021a; Kerr et al., 2021; Chang et al., 2021b). Among these, compartmental models are essential in infectious disease research due to their effective representation of the dynamics of disease transmission. In these models, the population is divided into several subgroups or compartments, each representing a specific stage in the progression of the disease. Consider, for example, the SIR model that divides the population into three compartments: Susceptible, Infectious, and Recovered. This model provides a basic framework for tracking disease spread by monitoring changes in population proportions within these compartments. Building on the SIR model, the SEIR model adds an Exposed compartment, highlighting a critical phase where individuals have contracted the disease but are not yet infectious. This addition is particularly important for diseases that have an incubation period, such as COVID-19. Several researchers extended the SEIR model to capture the unique transmission dynamics of the COVID-19 virus, as well as the impacts of vaccination and public health mandates on the spread of the disease. Kong et al. (2022) provide a detailed overview of these extensions, which consider various stages of disease progression, intervention strategies, and demographic factors (Kimathi et al., 2021; Tuite et al., 2020; Foy et al., 2021).

Classic epidemiological models assume that model parameters, such as the transmission rate, the recovery rate, and the contact rate, are constant throughout the course of the epidemic being modeled. Most recently, epidemiological models have been updated to incorporate the variability and uncertainty of these parameters. Models that take into account the randomness that exists in the system can lead to better decision-making. For example, (Lekone and Finkenstädt, 2006; Gray et al., 2011; Faranda and Alberti, 2020; Ganyani et al., 2021; Maltsev and Stern, 2021; Gatto et al., 2020; Giordano et al., 2020; Kretzschmar et al., 2020; Kucharski et al., 2020) propose stochastic epidemiological models. These models use historical data to estimate how parameters change over time, and thus, provide a more realistic representation of the dynamics of infectious diseases.

VH has consistently presented a significant challenge to the healthcare system, a concern that has been recognized from the World Health Organization (WHO) as a global health threat ((WHO, 2019). The advent of COVID-19 vaccines once again highlighted this issue. Several compartmental models capture the impact of vaccine availability and hesitancy on transmission rate. Buonomo et al. (2022) demonstrate that voluntary vaccination alone cannot stop the spread of the disease. Choi and Shim (2020) develop a game-theoretic epidemiological model to explore vaccination and social distancing as strategies to control the spread of the disease. They use a group behavioral model to determine optimal strategies for individuals. Awad et al. (2021) propose a compartmental model that uses geospatial and VH information to develop strategies to distribute vaccines. Jing et al. (2023) examine the role of VH in the emergence of new SARS-CoV-2 variants. Bui et al. (2024) present a compartmental model that captures the impact of VH on disease outcomes, such as the expected number of hospitalizations and deaths. The model runs for different values of VH (which are taken from a certain distribution) to estimate the impacts of uncertain VH on model outcomes. The authors performed a sensitivity analysis to assess the impact of the availability of healthcare resources and vaccines on the expected number of deaths. Despite the variety of compartmental models available, none examines the impacts of the temporal and spatial variations of VH on demand for healthcare resources. This gap underscores the need for models that use predictions of temporal and spatial disease spread to determine the effective allocation of critical resources.

Based on our review of the literature, epidemiological models are typically used to support decision making via “What-if" analysis. This approach explores how changes in model parameters or assumptions affect the spread of the disease and the need for resources. Some researchers integrate epidemiological models into optimization models to support optimal decision making during disease outbreaks. For example, Abdin et al. (2023) propose a non-linear programming formulation of a compartmental model to present the complex dynamics of infectious diseases, offering a nuanced perspective on disease transmission. Yin and Buyuktahtakin (2021a) propose an MSP compartmental model to support resource allocation during the Ebola outbreak in West Africa. The MSP uses information on stochastic disease progression to determine resource allocation to control the infectious disease outbreak. The model determines the fair allocation of critical resources.

Our study investigates the impact of vaccination and VH on disease transmission dynamics. The model we propose enables decision makers to assess the impact that the timing, scale, and scope of resource allocation have on the spread of the disease and disease outcomes. While previous research has addressed aspects of resource management during disease outbreaks, none has presented a comprehensive, data-driven model, as we propose. Our model highlights the advantages of merging the predictive capabilities of compartmental models with optimization techniques, providing valuable insights into effective disease control strategies.

MSP models are used to support sequential decision making under uncertainty. The simplest MSP model, the two-stage stochastic program (2-SP), is commonly used to minimize expected costs (maximize expected benefits) in decisions made today and tomorrow in the face of uncertainties faced tomorrow. MSPs are an extension of 2-SPs, which involves making decisions sequentially based on new information that emerges over time. For an overview of SP and MSP models and their applications, we refer the reader to Pereira and Pinto (1991); Goh et al. (2007); Zhang and Cardin (2017); Li and Grossmann (2021).

MSP models pose significant computational challenges that stem from the inherent uncertainty in the data and the intricate nested structure characteristic of multistage decision-making problems. A common approach to handling uncertainty is to approximate the underlying stochastic process using a scenario tree. A scenario tree branches at each stage to consider a finite number of possible scenarios. When the number of scenarios is small, the deterministic equivalent of the MSP problem can be solved using commercial optimization solvers, such as CPLEX or Gurobi. Lagrangian and Benders decomposition methods are often used to solve problems with a large number of scenarios. These methods break down the problem into smaller subproblems that are easier to solve (Laporte and Louveaux, 1993; Guignard, 2003). Rolling-horizon and heuristic approaches are also used to solve larger-scale MSP models (Christian and Cremaschi, 2015).

MSP models are often used for resource allocation in the healthcare sector, as the size and timing of the demand for these resources are often uncertain, and the availability of resources is limited. For example, Yin and Buyuktahtakin (2021b) develop an MSP model to effectively manage resource allocation during the 2018-2020 Ebola outbreak. In response to the COVID-19 pandemic, Yin et al. (2023) propose an MSP model for allocating ventilators across counties in New York and New Jersey. Their MSP incorporates a compartmental model to account for disease dynamics under varying levels of uncertainty. Hosseini-Motlagh et al. (2023) introduce a multi-stage fuzzy SP approach for patient allocation across healthcare facilities during COVID-19. Tanner et al. (2008) and Yarmand et al. (2014) propose SP models to determine optimal vaccination strategies. Mehrotra et al. (2020) propose an SP model for ventilator allocation during COVID-19, focusing on cost minimization. These publications highlight the relevance of MSPs to support sequential decisions related to resource allocation under uncertainty.

Integrating compartmental models within MSP resource allocation models provides a robust framework for planning responses to disease outbreaks. This approach enables high-quality decision making over time, taking into account the unpredictable nature of disease spread and resource availability. While some researchers are investigating this modeling approach, there is more to be done. Our proposed data-driven compartmental MSP model addresses VH and fairness in resource allocation within this modeling framework. Our proposed model supports effective vaccination strategies and resource allocations that ensure greater immunity and mitigate the impact of outbreaks.

A WHO report highlights the importance of following ethical principles in the allocation of healthcare resources by analyzing the utility, efficiency, and fairness approaches taken by decision makers (WHO, 2004). Utilitarian decision makers implement strategies that yield the highest overall benefit or the lowest overall cost. Efficient decision makers advocate using the least amount of resources to achieve the maximum possible impact. Fair decision makers aim to ensure equitable allocation and reduce disparities between groups and individuals. Aligning resource allocation strategies with these principles can be challenging as each strategy, while meeting one principle, falls short of meeting others. The work of Lane et al. (2017) highlights the complexity of ensuring equity in allocating healthcare resources, emphasizing the need for clear definitions that reflect societal values to ensure equal access to healthcare, address patient needs effectively, and distribute potential benefits fairly. Equitable strategies can build public trust and encourage cooperation, which are essential during health crises.

Several OR models for resource allocation take a utilitarian approach to optimize logistics and operations to control disease spread (Sun et al., 2014; Liu et al., 2015; Liu and Zhang, 2016; Zaric and Brandeau, 2001; Day et al., 2020). Other studies introduce fairness metrics and propose models that optimize fairness in resource allocation. For example, Orgut et al. (2016) proposes a model that allocates and distributes food donations proportionally to demand, aiming for an equitable distribution while minimizing waste. Similarly, Davis et al. (2015) propose a multi-step optimization approach to improve fairness in the allocation of organs (kidneys) across the U.S., thereby improving equity of organ transplant. Dönmez et al. (2022) presents a multi-objective, multi-period, nonlinear model for equitable personal protective equipment (PPE) allocation to health centers in pandemics, with the aim of minimizing infections and ensuring fair distribution amidst resource scarcity. The research by Yin and Buyuktahtakin (2021a) introduces an MSP model for a fair allocation of resources, such as beds and treatment centers, during the West African Ebola outbreak, balancing efficiency and equity in epidemic responses.

Similarly to the literature, our proposed model allows decision makers to evaluate the trade-offs between fairness and efficiency in allocating critical healthcare resources. Unique to our approach is the incorporation of uncertainty in VH, which introduces additional complexity to resource allocation strategies. We analyze the ‘price of fairness’, the efficiency loss incurred when striving for equitable resource allocation over time and across different geographic regions. The MSP compartmental model adapts to the evolving landscape of healthcare needs, ensuring that resource allocation remains responsive and equitable even under unpredictable changes in public sentiment towards vaccination.

Motivation: Our proposed model considers unccertainty in VH. This is motivated by observations made from the COVID-19 Trends and Impact Survey (CTIS) conducted by the Delphi Group at Carnegie Mellon University, in collaboration with Facebook (Salomon et al., 2021), and from data about vaccination rates collected by the Centers for Disease Control and Prevention CDC (2021).

Figure 1(a) presents VH estimates of every county in the states of Missouri (MO) and Rhode Island (RI), providing a snapshot of the willingness of the public to vaccinate. The VH estimate is derived from responses to the following question in the CTIS survey: “If a vaccine to prevent COVID-19 was offered to you today, would you choose to get vaccinated?”. Participants were asked to select one of the following answers “Yes, definitely,” “Yes, probably,” “No, probably not,” or “No, definitely not.” The last three responses are considered indicative of VH. Each thin line represents VH progression of an individual county. The thicker line depicts the average VH trend for the respective state.

Figure 1(b) shows the vaccination rates in counties of Missouri and Rhode Island based on data from CDC (2021). The data show fully vaccinated individuals who have received the complete dose - either two doses of a two-dose vaccine or one dose of a single-dose vaccine. These data highlight a general decline in VH over time and underscore the varied vaccination sentiments in counties. Figure 1(c) illustrates a negative correlation between average vaccination rates and average VH estimates with correlation coefficients of $\rho=-0.85$ for Rhode Island and $\rho=-0.93$ for Missouri. The high correlation coefficients in both states indicate that higher vaccination rates are associated with lower VH.

Within a state, varying attitudes towards vaccination and differing vaccination rates significantly impact disease transmission and the resulting need for healthcare resources. We capture variations in VH over time and space via a scenario tree that we describe next.

Scenario tree: Scenario trees provide a systematic approach to represent uncertainty that is suitable for quantitative models. They provide likely scenarios of the future with associated probabilities.

Each scenario $\omega$ ($\omega\in\Omega$) represents the changes of VH for a specific population group, over time. We derive temporal changes in VH from the historical data provided by the CTIS. We fit a distribution to the rate of change of VH. Numerical analysis shows that the (monthly) rate of change of VH is normally distributed. The corresponding mean ($\mu_{t}$) and standard deviation ($\sigma_{t}$) change from one period (month) to the next (see Figure 12). To construct a decision tree of manageable size, we approximate the normal distribution with a discrete triangular distribution that assumes the values $\mu_{t},\mu_{t}-\sigma_{t},\mu_{t}+\sigma_{t}$ with probabilities of 0.158, 0.684, and 0.158, respectively.

Figure 2 presents the scenario tree for a problem with $T=3$ stages. The resulting scenario tree has $3^{2}=9$ unique scenarios. Let $p_{\omega}$ denote the probability associated with each of these scenarios. The tree showcases a baseline (green path), an optimistic (blue path), and a pessimistic (orange path) scenario. In the optimistic scenario, the VH changes by $\mu_{t}-\sigma_{t}$ in each stage. In the pessimistic scenario, the VH changes by $\mu_{t}+\sigma_{t}$ at each stage. In the baseline scenario, the VH changes by $\mu_{t}$ at each stage. Consequently, the probability associated with the pessimistic and optimistic scenarios is $0.158^{2}=0.024$. The probability associated with the baseline scenarios is $0.684^{2}=0.468$.

The SEIR model divides the population into four groups (compartments): Susceptible ($\boldsymbol{S}$), Exposed ($\boldsymbol{E}$), Infectious ($\boldsymbol{I}$), and Recovered ($\boldsymbol{R}$). The model captures the transition of individuals through these compartments over time. It also accounts for the latency period before an individual becomes infectious, which is important in modeling the spread of diseases such as COVID-19. In a previous study, our team extended this model to address the specific complexities of COVID-19, such as vaccine availability, VH, and the scarcity of critical healthcare resources in the spread of the disease (Bui et al., 2024). For the convenience of the reader, we summarize the model in the following paragraphs.

A schematic representation of the proposed susceptible-vaccinated-exposed-infectious-hospitalized-recovered (SVEIHR) model is shown in Figure 3. In this model, susceptible individuals can be exposed or opt to vaccinate, moving forward to the vaccinated compartment (V). The size of compartment V is dictated by $\rho$, representing the maximum vaccination rate achievable in the absence of VH. The willingness of people to get vaccinated is represented by $(1-\tilde{h})$, where, $\tilde{h}$ is a stochastic parameter that represents the proportion of the population reluctant to vaccinate. Furthermore, the model considers that vaccination does not confer absolute immunity; therefore, vaccinated individuals can progress to an Exposed but Vaccinated compartment ($\boldsymbol{EV}$), with a probability determined by vaccine efficacy, $\epsilon$.

The model captures different postexposure health trajectories based on incubation rates $\alpha$ and the severity of the disease. Individuals progress to Mild ($\boldsymbol{I^{\text{m}}}$), Severe ($\boldsymbol{I^{\text{s}}}$), or Critical Condition ($\boldsymbol{H^{\text{c}}}$) compartments with probabilities $p^{\text{m}},p^{\text{s}},$ and $p^{\text{c}}$, respectively. Vaccinated individuals are less likely to experience severe/critical health outcomes; thus, several infected individuals are likely to recover and move to the compartment ($\boldsymbol{R}$).

The model factors the availability of critical healthcare resources, such as hospital beds (Bed) and ventilators (Vent), in determining patients’ health outcomes. Refer to the example presented in Appendix A for a demonstration of the impact of the number of beds and ventilators on the model outcomes. Recovery ($\gamma$) and mortality ($\mu$) rates vary depending on the severity of the disease and the availability of healthcare resources. The model accounts for the transition of critically ill patients to the compartment ($\boldsymbol{D}$) when ventilators are unavailable, emphasizing the dire consequences of resource limitations during a pandemic.

This section outlines the formulation of the proposed multi-stage stochastic program (MSP) model. The model determines an optimal sequence of decisions related to ventilator allocation across various regions over time. The MSP minimizes the total number of deaths by dynamically adjusting the distribution of healthcare resources in response to changing VH rates and state transitions within the SVEIHR model. Resource allocation decisions take into account not just immediate health outcomes, but also the need for equitable access to care. By integrating the compartmental model with the MSP model, we now better understand how VH and disease dynamics affects healthcare systems, and use this information to optimize the allocation of ventilators. We refer to this proposed model as the compartmental MSP model. For a detailed description of the model notations used, please see Appendix C.

Model assumptions: Our model operates on three primary sets of assumptions that are critical to its structure and function:

Assumption 1: Uncertainty in vaccine hesitancy (VH): The model considers three potential outcomes of VH in each decision stage. Let ${h}_{t-1}$ represent the proportion of the population reluctant to vaccinate in period $t-1$. Then, in period $t$, the proportion of the population reluctant to vaccinate for the optimistic, expected, and pessimistic scenarios are $h_{t-1}(1+\mu_{t}-\sigma_{t});h_{t-1}(1+\mu_{t})$; and $h_{t-1}(1+\mu_{t}+\sigma_{t})$. Each of these outcomes is associated with a fixed probability that does not change over the course of the study. These probabilities remain constant and unaffected by changing decisions. We note that decision makers can customize these outcomes (number of outcomes and values) based on available data.

Assumption 2: Disease dynamics: Since SVEIHR is a compartmental epidemiological model, the following assumptions of SEIR model apply. ($i$) Every individual in a population is equally likely to interact with others. As a result, disease progression and transitions between compartments are consistent throughout the population. However, it should be noted that in a real-world setting, elements such as social networks and age demographics can significantly influence disease transmission patterns. ($ii$) The size of the total population does not change due to natural births, non-COVID-19-related deaths, or migration. Since the time frame considered in our study spans several months, recovered individuals are assumed to have full and lasting immunity to the disease. ($iii$) Infection, recovery, and exposure rates are assumed to be deterministic. However, we recognize that in practice, they can fluctuate for many reasons, from policy interventions and demographic differences to viral mutations.

Assumption 3: Healthcare system: The available healthcare resources (i.e., beds, ventilators) are only used to treat COVID-19 patients. Furthermore, our proposed MSP model does not account for potential reallocations of healthcare resources.

Mathematical model formulation: Classical SIR are continuous-time models. However, several researchers propose discrete-time approximations of these models that are easier to solve (Wacker and Schlüter, 2020). These studies show that many of the desired properties of the time-continuous SIR models remain valid in the time-discrete approximations. We use a discrete-time model to mimic the dynamic changes in the size of compartments of the SVEIHR model over time. Let $\mathcal{T}$ represent the set of time periods in the planning horizon. Let $\mathcal{J}$ represent the set of decision periods (stages). Note that, between every two consecutive decision periods, there are several time periods, thus $\mathcal{J}\subset\mathcal{T}$. We denote by $\mathcal{R}$ the set of regions studied. We distinguish between a decision stage $j$ and a time period $t$ to keep the discrete approximation of the SEIR model effective, ensuring short time intervals between transitions, while avoiding too many decision stages for the MSP model.

In a decision period (stage), we determine the allocation of ventilators. Our decision variables are $x_{\omega,j,r}$, which represent the number of ventilators assigned to the region $r\in\mathcal{R}$ at the beginning of period $j\in\mathcal{J}$ in scenario $\omega\in\Omega$. Let $\mathcal{X}_{\omega,t,r}$ represent the cumulative number of ventilators allocated to region $r$ by time $t$.

The decisions made affect the state variables, which are the number of susceptible individuals who are unvaccinated ($S_{\omega,t,r}$) and vaccinated ($V_{\omega,t,r}$); the number of exposed individuals who are unvaccinated ($E_{\omega,t,r}$) and vaccinated ($EV_{\omega,t,r}$); the number of infectious individuals with mild ($I^{\text{m}}_{\omega,t,r}$) and and severe ($I^{\text{s}}_{\omega,t,r}$) symptoms; the number of hospitalized individuals in severe ($H^{\text{s}}_{\omega,t,r}$) and in critical ($H^{\text{c}}_{\omega,t,r}$) condition; the number of individuals admitted to a hospital and in severe ($\mathcal{A}^{\text{s}}_{\omega,t,r}$) or critical ($\mathcal{A}^{\text{c}}_{\omega,t,r}$ ) condition; the number of individuals in severe ( $\mathcal{K}^{\text{s}}_{\omega,t,r}$) and critical ($\mathcal{K}^{\text{c}}_{\omega,t,r}$) condition not admitted to the hospital due to capacity constraints; the number of recovered individuals ($R_{\omega,t,r}$); and the number of deaths ($D_{\omega,t,r}$). Let $\mathcal{L}=\{S,V,E,EV,I^{\text{m}},I^{\text{s}},H^{\text{s}},H^{\text{c}},R,D\}$ represent the set of our state variables.

Next, we provide detailed descriptions of the objective function and constraints of the proposed compartmental MSP.

Objective function: The objective is to minimize the total expected number of deaths across regions during the planning horizon.

Constraints: We organize our constraints into five groups (G1, G2,…, G5), each tailored to fulfill a specific function. This structure enhances the model’s clarity and highlights the essential contribution of each constraint category to the model’s objective.

(G1:) SVEIHR compartmental model: This group of constraints captures the dynamic changes of state variables over time. We begin by initializing the state variables $\ell\in\mathcal{L}$ at $t=0$. Let $\boldsymbol{\pi}_{\ell,r}$ represent the corresponding initial values.

Constraints (3)-(12) describe the transition of the population between different compartments of the epidemiological model over time (as illustrated in Figure  3). Constraint (3) updates the size of the susceptible population in compartment S. To determine the size of S at the end of period $t+1$, we update the size of S at the end of period $t$ by subtracting the number of individuals exposed, vaccinated and subtracting (adding) those who immigrated (emigrated) to the region from (to) other regions considered in the analysis. In Appendix F, we provide descriptions of the effects of migration in our model. Constraint (4) updates the size of the vaccinated population, compartment V. To determine the size of V at the end of the period $t+1$, we update the size of V at the end of period $t$ by subtracting the number of individuals exposed and adding the number of individuals vaccinated. The efficacy rate of the vaccine, $\epsilon$, affects the size of this compartment. Constraint (5) updates the size of the unvaccinated population who were exposed, compartment E. To determine the size of E at the end of period $t+1$, we update the size of E at the end of period $t$ subtracting the number of individuals exposed who became infected, and by adding the number of individuals exposed during period $t+1$. Constraint (6) updates the size of the vaccinated population who are exposed, compartment EV. This constraint updates the size of EV by taking a similar approach to  (5). The difference between constraints (5) and (6) is that vaccinated individuals can recover after exposure and do not become critically ill. Constraints (7) and (8) update the size of compartments $I^{\text{m}}$ and $I^{\text{c}}$, respectively. To determine the size of compartment $I^{\text{m}}$ at the end of $t+1$, constraint (7) considers new infections, and recoveries. To determine the size of compartment $I^{\text{s}}$ at the end of $t+1$, constraint (8) considers new infections, recoveries, deaths and transitions to the hospitalization compartment. Constraints (9) and (10) update the size of compartments $H^{\text{s}}$ and $H^{\text{c}}$. To determine the number of severely hospitalized individuals (the size of compartment $H^{\text{s}}$) at the end of period $t+1$, constraint (9) considers new admissions, recoveries, and deaths. To determine the number of critically hospitalized individuals (the size of compartment $H^{\text{c}}$) at the end of period $t+1$, constraint (10) also considers new admissions, recoveries, and deaths. Constraint (11) determines the number of recovered individuals (the size of the compartment R) at the end of the period $t+1$ by aggregating the number of those recovered from various compartments (of infected and hospitalized individuals), each with specific recovery rates. Lastly, the constraint (12) determines the size of the compartment D. It compiles the death toll within the region, including deaths from the inability to hospitalize individuals in severe and critical conditions, and deaths among individuals in critical conditions, adjusting to their respective mortality rates.

Constraints (3)-(12) are for all $t\in\{1,\ldots,|\mathcal{T}|-1\},r\in\mathcal{R},\omega\in\Omega$:

(G2) Resource capacity constraints: Constraints (13) to (15) determine the values of the decision variables. The initial number of ventilators in hospitals is limited. Consider that $\Delta_{j}$ is the number of new ventilators available for allocation in period (stage) $j\in\mathcal{J}$. Constraints (13) determine the allocation of these ventilators to hospitals in different regions. Constraints (14) determine the total number of ventilators available in region $r$ at period $j$. Since the number of time periods between every two consecutive decision periods is greater than one, we ensure (via constraints (15)) that the total number of available ventilators does not change until the next decision period.

Constraint (16) connects the decision variables with the state variables. This constraint determines the number of critically ill patients admitted to a hospital, $\mathcal{A}^{\text{c}}_{\omega,t,r}$. This number depends on the demand for critical care, the number of ventilators available, and the number of beds available.

Constraint (17) determines the number of severely ill patients admitted to a hospital, $\mathcal{A}^{\text{s}}_{\omega,t,r}$. This number depends on the demand for care and beds available after the admission of patients in critical condition. In Appendix D we detail the linearization of “min" functions in constraints (16) and (17). The restrictions (18) and (19) determine the number of patients with severe and critical conditions who were not admitted to a hospital due to limited resources. Constraints (16) to (19) are for all $t\in\mathcal{T},r\in\mathcal{R},\omega\in\Omega$.

(G3:) Resource allocation constraints: Equity in healthcare resource distribution, especially during emergencies, is vital for public health policy. Defining and measuring equity in public health decisions is a complex task with no widely accepted standard. We present several approaches used by the literature to ensure fair allocation of healthcare resources (see our discussion of (20)-(22)). Each constraint, when applied individually, offers a distinct approach to resource allocation. Constraint (20), adapted from (Yin and Buyuktahtakin, 2021a), uses information about the proportion of the expected critically ill individuals in a region (the first term in this constraint) and the corresponding population proportion of the region (the second term in this constraint) to allocate ventilators. It uses an equity threshold parameter, denoted by $k\in(0,1]$, to gauge need-based equity in resource allocation. As $k$ approaches zero, equity in resource allocation is enforced. As the value of $k$ approaches one, the equity requirements become less stringent. In Appendix D, you will find our approach to linearizing this constraint.

The restrictions (21) present a variation in the population-proportional method for the allocation of resources (Dangerfield et al., 2019; World Health Organization, 2020). We include an equity threshold parameter $\zeta$. This constraint ensures that each region receives ventilators in proportion to its share of the total population, with $\zeta$ providing adjustable scaling to fine-tune how strictly the allocation of resources follows population size. For $\zeta=1$, the additional ventilators available for use during period $j$, denoted as $\Delta_{j}$, are distributed according to the population size of each region. As $\zeta$ decreases, this constraint relaxes, and the equity requirements become less stringent.

Both constraints (20) and (21), focus on ensuring an equitable distribution of critical resources. One of the constraints uses a need-based approach, and the other population-based approach to resource allocation.

Constraints (22) ensure that ventilators are distributed equally, ignoring differences in population size or the specific needs of each region. As a result, the same number of ventilators is allocated to every region.

Notice that we solve the compartmental MSP model with either constraints (20), or (21), or (22) to enforce equity or equality in resource allocation. We also solve the compartmental MSP model without constraints  (20),  (21), and (22) to demonstrate the effects of a utilitarian approach that aims to minimize the overall cost of resource allocation.

(G4) Non-anticipativity constraints (NACs): These constraints ensure that decisions made at a given time (decision stage) are based solely on the information available up to that point. This prevents the use of future, yet unknown, information in decision making. Specifically, constraints (23) ensure that decisions that share an identical scenario path up to time $j$ must be consistent and identical across shared paths. In (23), $\mathcal{S}_{\omega,j,r}$ represents the set of scenarios that are indistinguishable from $\omega$ at time $t$. In Appendix E, we provide an example that demonstrates the implementation of these constraints.

(G5) Other constraints: Constraints (24) restrict the decision variables to be non-negative integers. Constraints (25) and (26) restrict the state variables to be nonnegative. The following constraints are also for all $\omega\in\Omega,r\in\mathcal{R}$.

Data collection: We collect data about county-level vaccination uptake from January to May 2021 from the CDC (2021), and data about \AcVH from the CTIS dataset (Salomon et al., 2021). Data on healthcare resources come from the CovidCareMap (2020) project. Data are organized by facility, county, referral region, and state. Table 1 summarizes the data for Arkansas.

Figure 4(a) presents the county-level distribution of hospitals and licensed beds in Arkansas at the beginning of the COVID-19 pandemic. Figure 4(b) illustrates the relationship between population size and the number of licensed beds in each county. We cluster counties of Arkansas into four regions based on their VH. Below, we provide details of our clustering approach. Figure 4(c) summarizes the total number of licensed beds and intensive care unit (ICU) beds in each region. We assume that each ICU bed has a ventilator. The results of these figures emphasize the variations in healthcare capacity across the state.

Clustering of the data: In an effort to reduce the size of the decision tree of the MSP model, we cluster counties of Arkansas into four regions. The counties are clustered according to their VH similarity. We measure similarity using the Euclidean distance between county-level VH time series data. We use the SciPy Python dendrogram tool to visualize the clustering of counties that belong to the same region (Figure 5(a)) (Virtanen et al., 2020). We use as a threshold a distance of 0.4 and identify four distinct regions. Figures 5(b)-(e) illustrate VH trends within each region. Each thin line represents the VH of a county, and the dotted line represents the average VH of that region. In particular, Region 3, marked red, starts with the highest VH and, despite the declining trend, remains the highest throughout the study period. In contrast, Region 2, marked green, begins with a VH of 0.39, which then drops dramatically, ending as the region with the lowest VH. Figure 5(f) maps the regions geographically, highlighting that Arkansas’s major cities, Fayetteville, Fort Smith, and Little Rock, are in Region 2, which has the lowest VH.

Other data processing: Our model employs the parameter $\Delta_{j}$ to indicate the number of ventilators available for distribution among different regions at the start of each decision period, $j$. The COVID-19 pandemic significantly increased the demand for ventilators, a critical medical device that saved many lives. In response, under Defense Production Act (DPA), various manufacturers, including auto manufacturers, quickly began to increase the production of ventilators or repurpose their facilities for this task (Albergotti and Siddiqui, 2020). For example, General Motors (GM) was contracted to supply 30,000 medical ventilators for the national stockpile (Massey and Cole, 2020). In April 2020, AdvaMed reported that its member companies were ramping up production from an average of 2,000 to 3,000 ventilators per week to an expected 5,000 to 7,000 per week. This was a notable increase from the 700 ventilators they produced weekly for the U.S. market in 2019 AdvaMed (2020). Due to the lack of real-world data related to the additional ventilators allocated to Arkansas after the enactment of DPA, we use the following approach to generate these data.

We use equations (27) to (29) to determine $\Delta_{j}$. Let $\underline{p}$ represent the initial stockpile of ventilators in Arkansas prior to enactment of DPA, and let $\overline{p}$ represent an upper bound on the total number of available ventilators. Let $p_{j}$ represent the additional ventilators made available in period $j$.

We begin with an initial stockpile ($\underline{p}$) of 100 ventilators. In each subsequent time period, we increase the number of ventilators ($p_{j}$) by 50. As a result, the number of ventilators available for distribution in Arkansas follows this pattern: $\Delta_{1}=100$, $\Delta_{5}=150$, $\Delta_{9}=200$, $\Delta_{13}=250$, and $\Delta_{17}=300$.

Please refer to tables in Appendix G for a complete list of the parameter values used in our MSP compartmental model. The table also presents the tuned parameters, which are outcomes of the model calibration process detailed in Section 4.2.

The calibration of SVEIHR focuses on fine-tuning the model parameters to ensure that it accurately represents the dynamics of disease progression in each region of study. We collected data from COVID-19 Open Data (Wahltinez et al., 2020). This dataset provides county-level data. For the purpose of this study, we aggregate county-level data for Arkansas into the four aforementioned regions.

We use data from two critical time periods in the U.S. The first spans 20 weeks, from October 25, 2020 to March 7, 2021. This period includes the peak of the wave of coronavirus infections during the Fall of 2020. We use these data to evaluate SVEIHR’s power in predicting new infections. The second spans 20 weeks, from January 10 to May 23, 2021, and corresponds to the initial stages of vaccine rollout. We use these data to evaluate SVEIHR’s power in predicting vaccinations. These study periods do not coincide with the start of the pandemic. Therefore, we adjust the size of the initial susceptible population to account for individuals who were infected, recovered, or died prior to the start of the study period.

Figure 6 summarizes the results of our model calibration of the baseline scenario. Each plot in this figure presents the actual data and SVEIHR model predictions of the weekly number of infections and vaccination intakes for each region of study. New infections at time $t$ are calculated using this equation: $\alpha(p^{\text{m}}+p^{\text{s}})E_{t,r}+\alpha(p^{\text{m},v}+p^{\text{s},v})EV_{t,r}$. The cumulative number of vaccinations up to week $t$ is calculated using $\sum_{t}(1-h_{t,r})\rho_{r}S_{t,r}$. We employ the Optuna framework to optimize parameter tuning (Akiba et al., 2019). Table 2 summarizes the different performance metrics that we use to evaluate the performance of SVEIHR model of the baseline scenario, including root-mean-square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). The low value of MAPE highlights the accuracy of SVEIHR in tracking new infections and vaccination trends in each region.

The value of stochastic solution (VSS) quantifies the benefit of incorporating uncertainty into decision-making processes by comparing the results of the compartmental MSP with those of its deterministic counterpart. To compute the VSS for each decision stage $j\in\mathcal{J}$, we adopt the methodology proposed by Escudero et al. (2007). In this method, the expected value (EV) model, which operates on the expected values of uncertain parameters, is used as the deterministic benchmark. Let $\mathit{EEV}_{j}$ denote the optimal value of the proposed compartmental MSP model when the solution up to the decision period immediately preceding $j$ is fixed to the solution from the deterministic model. Let $\mathcal{Z}^{*}$ denote the optimal objective function value of the compartmental MSP model. Therefore, VSS at stage $j$, represented as $\mathit{VSS}_{j}=\mathit{EEV}_{j}-\mathcal{Z}^{*}$, captures the value derived from accommodating uncertainty in the model up to decision stage $j$.

The evolution of the VSS across the decision stages is illustrated in Figure 7(a). For the model used in this illustration, there are 5 decision stages that correspond to weeks 1, 5, 9, 13, and 17. The results of this figure highlight the benefits derived from the proposed compartmental MSP model. The VSS starts at zero, which is expected since $\mathit{EEV}_{1}=\mathcal{Z}^{*}$. As the model progresses through the stages, the VSS trends upward, reflecting the benefits of incorporating uncertainty into strategic decision making. We notice a decrease in the rate of change in VSS from week 13 to 17. This is mainly because the mean change in VH is lower during these weeks. The results in Figure 12 show that the mean change in VH from week 13 to week 17 ranged from 0.03 to 0.05, which is lower than 0.06-0.09 seen from weeks 9 to 13.

Note that $VSS_{j}=87$ in the illustrative example in Figure 7(a) appears to be a rather small number compared to the optimal solution $\mathcal{Z}^{*}=43,900.$ In this example, a critically ill patient uses a ventilator for three weeks, the recovery rate is 51% and $\overline{p}=300$. This leads to a total of 1,000 ventilators available for allocation and use over five months, and a maximum of 729 critically ill patients to be saved using these ventilators. Thus, a $VSS_{j}=87$ represents 12.21% of this total, which is pretty high and highlights the value of using the compartmental MSP.