Optimizing a multimodal hub and spoke network for vaccine distribution using regional vulnerability index during healthcare emergencies Open Access

The emergence of SARS-CoV-2 as a novel virus, causing vaccine preventable diseases (VPDs) leading to healthcare emergencies, initially perplexed the global medical community due to its unique contagious behavior. As pharmaceutical approaches to combat the virus were lacking experience, many administrations worldwide initiated non-pharmaceutical approaches as precautions. However, these approaches have proven to be insufficient, as evidenced by the high mortality rates observed during the first and second waves of the pandemic. After observing and analyzing the contagious behavior of the virus, experts suggest that vaccination is the only pharmaceutical solution to prevent infection (Moore et al., 2020). For the implementation of such pharmaceutical approaches as vaccination and its distribution, an analytical estimation of the peak time infection for different regions is necessary (Kröger et al., 2021). Such analytical approaches have been witnessed in the literature for estimating and predicting peak time infection during the coronavirus infection spread. Most such illustrations were presented using the SIR model (Brugnano et al., 2021; Cooper et al., 2020; Kröger et al., 2021; Law et al., 2020). For example, Turkyilmazoglu (2021) has presented an approximate expression to predict the peak time of an epidemic. They also suggest implementing pharmaceutical approaches using the approximation of peak time from the SIR model. However, none of the studies in the literature has used these estimations and predictions for prioritizing the vulnerable regions while distributing vaccines and other medical essentials. This motivates us to quantify the Susceptibility of a region based on the peak time infection for formulating a priority-led optimum distribution model during vaccine scarcity.

Post-vaccine development for the most recent health emergency of COVID-19, several scholars have worked in the area of effective vaccine distribution (Dinleyici et al., 2021; Mandal et al., 2021; Yang et al., 2021). However, there is almost no attention paid to a priority-based vaccine distribution strategy. To address this, we propose a two-phase approach. Phase 1 proposes a metric called Regional Vulnerability Index (RVI) using the static ModSIR model (a modified SIR model) by analyzing the infection pattern in the initial state of infection spread. The index sheds light on the probable epicenters among the different regions in a country. This index is used to grade the vulnerable areas and integrated with Phase 2, which is a P-Center problem for facility locations. In phase 2, we formulated this problem as a Mixed Integer Quadratically Constrained Programming (MIQCP) and proposed three versions of it with different RVI-based weightage schemes, each for a particular scenario. While the first version minimizes the vulnerability-weighted distance to cover susceptible regions by building facility hubs close to vulnerable regions, the second and third versions extend it to minimize the distance weighted with the inverse of RVI and the ratio of population to RVI respectively. The purpose of proposing three variations of the MIQCP with different weighting schemes is to identify the best variation applicable for a fair distribution of vaccines in terms of coverage. However, during a pandemic with insufficient resources, vaccination strategies must prioritize the coverage of more susceptible regions. Thus, we address both population and vulnerability coverage.

Apart from bridging the literature gap in the area of epidemiology-based priority facility location, Phase 2 also addresses the void in multimodal facility-based literature on health emergencies. Although there are studies on multimodal transportation, in most of them, the choice of transportation mode is predefined rather than a decision (Enayati et al., 2023; Li et al., 2023; Ruan et al., 2014, 2016; Zhang and Lam, 2018). Moreover, most emergency-based distribution studies do not consider multimodality, which we believe, is extremely important under such scenarios. This consideration gathers all the resources available under administrative control for transporting medical essentials faster using vehicles suitable for each pair of locations. To address this, the MIQCP model considers multimodality by considering a heterogeneous fleet with vehicles from multiple modes of transportation. The problem setting that the model solves considers a set of locations as candidates for facility locations to serve a set of demand nodes. For transporting vaccines to these demand nodes, the model considers the options of multiple vehicles from multiple modes of transportation based on the available infrastructure. The model decides on the location for opening vaccine facilities, number of vaccine doses to be transported to each node, the vehicle type to be used and the number of trips to be required.

To demonstrate the model, we have collected necessary data (GPS locations, Distance Matrix and Population data) and implemented them in the case of all states and Union Territories (UTs) in India. The results interpret an optimally generated network of hub locations connecting such administrative divisions. We propose a few network metrics namely Degree Centrality, weighted Degree Centrality, Closeness Centrality, average distance and maximum distance from hubs across each set.

The proposed study poses the following four research questions: Average distance of the most vulnerable regions from the facility., how can the classical SIR model be modified to better estimate peak infection during healthcare emergencies? Second, how can insights from epidemiological models be used to prioritize vaccine distribution during health emergencies? Third, how can fairness in population and vulnerability coverage be assessed and validated through network property-based metrics in facility location problems (FLP)? Fourth, in what ways does the introduction of multimodality in the P-Centre FLP enhance the efficiency of vaccine storage and distribution in healthcare emergencies? While answering these questions, the study makes five distinct contributions, which are listed as follows:

• We use a modified SIR model called the ModSIR model for estimating the peak infection rate and peak infection time of a region.

• We propose a Vulnerability index, RVI using the ModSIR model to quantify the susceptibility of a region towards being a potentially adversely affected region under healthcare emergencies.

• We propose three versions of non-linear MIQCP P-Centre formulations for location problems by prioritizing the vulnerable regions using three weighting schemes and discussing their practical application in the healthcare emergency context.

• We introduce multimodality in the P-Centre FLP formulation for vaccine storage and distribution in healthcare emergencies.

• We are the first to validate the fairness in population and vulnerability coverage achieved from FLPs using several network/graph property metrics.

The rest of the paper is organized as follows: The introduction is followed by a literature review in Section 2, where we present the literature on the SIR model and FLPs. Section 3 presents the proposed approach for estimating RVI from the ModSIR model and integrating it into MIQCPs with the mathematical model. Section 4 explains the case study of India where we demonstrate the proposed approach. Section 5 delineates the results and discussion by analyzing the comparison between existing and proposed models and concludes the paper. We present an  Appendix explaining the modified SIR model, ModSIR with an algorithmic approach.

In this section, we discuss the related work in the area of priority-based FLPs in five distinct subsections. In this context, we first review the FLPs in the area of healthcare. Then we discuss the priority-based approach in healthcare emergencies in the pre-COVID and post-COVID era. Moreover, this section also discusses the epidemic models to estimate the peak time of infection and the multimodality in distribution during health emergencies.

The healthcare facility location (HFL) problem has attracted several researchers from the Operations Research community (Ahmadi-Javid et al., 2017). This survey-based study has mentioned that several HFL problems are based on p-center problems. For example, Andersson and Värbrand (2007) propose a model for ambulance allocation by formulating a p-center problem. They have used MILP to solve the same. Kim and Kim (2010) have specifically worked on long-term nursing centers using a p-center problem. de Aguiar and Mota (2012) propose an updated p-center problem for large-scale emergencies. Apart from that, there are also some recent studies addressing vaccine FLPs using p-center problems for uncertain demands (Baron et al., 2019; Lin and Wu, 2022). However, vaccine distribution during pandemics requires some prioritization strategy for its diverse effect on certain groups. Chowdhury et al. (2023) have designed HFL models considering the domino effect associated with disasters like earthquakes where this even is followed by a series of other events like infrastructure collapse, property damage, human injuries, etc. Nevertheless, another humanitarian problem that comes very close to vaccine distribution during a health emergency is the location problem for relief material distribution during emergencies. In this context, Guan et al. (2021) have modeled a location problem for the distribution of relief material during an earthquake. Mollah et al. (2018) have designed a cost-optimization model for relief material distribution during floods. They focused on the solution methodology rather than addressing the priority. Yan et al. (2023) have designed optimization-driven models for emergency relief material. Zhang et al. (2023) have introduced a min-max robust model for identifying locations for relief kit assembly and distribution. Seranilla and Löhndorf (2024) proposed a stochastic-dynamic facility location model that addresses vaccine distribution challenges in developing countries, accounting for facility disruptions due to natural disasters and optimizing cost reductions through approximate dynamic programming techniques. To the best of our knowledge, there are no priority Humanitarian or Healthcare Facility Location problems built around vaccine distribution under health emergencies like pandemics.

In the literature, FLP based on priority covers diverse prioritization criteria. Explored priority FLPs are based on the neediness of rural areas (Datta, 2012), proximity to certain locations, population density for park and ride problems (Rezaei et al., 2022), and prioritizing emerging markets (Mokrini et al., 2019), etc. In the context of HFL, Carlström et al. (2022) have mentioned one priority-based problem focused on vaccination strategies for prioritizing a vulnerable group of the population for influenza vaccine distribution (Huang et al., 2010). Later, Emu et al. (2021) integrated this priority scheme with a distance-based vaccine distribution model. We came across COVID-related vaccine distribution at the last echelon and ultra-cold vaccine distribution (Shukla et al., 2022; Xu et al., 2021). In a recent study, Eksioglu et al. (2024) developed a mixed-integer linear programming model to design drone delivery networks for vaccine supply chains in low-income countries, improving vaccine availability. In an interesting study, Castillo-Neyra et al. (2024) developed optimization algorithms to enhance mass pet vaccination campaigns by minimizing walking distance to vaccination sites, thereby increasing participation and improving spatial coverage. However, most of the prioritized facility location problems are from the pre-COVID period and are based on certain age groups or experts’ opinions rather than any mathematical computation.

Post-Covid, several researchers have proposed different priority schemes supported by mathematical computation for utilization in different sectors (Priya et al., 2022; Qrunfleh et al., 2022). For example, Fu et al. (2021) prioritize different age groups based on transmission rates to model the epidemiology of COVID-19. Babus et al. (2020) have worked on a vaccine allocation problem by comparing age groups and occupation groups. In the context of vaccine facility locations, Bravo et al. (2022) have used the priority-based FLP for targeting the “vaccine deserted” regions in the USA. Recently, Yin et al. (2024) used prioritization in their vaccine allocation model by focusing on regions with higher population density and infection rates. Their model allocates vaccines based on a combination of factors, including the budget level, vaccine costs and vaccine efficacy. Apart from that, we did not come across any priority-based HFL problems in the literature. However, the relevance of these priority schemes under circumstances like pandemics is yet to be studied.

The peak time of infection is one of the most important entities in all the epidemic predictive models found in the literature (Bhardwaj, 2020; Nishiura, 2020; Rodriguez-Palacios et al., 2019). This entity helps the epidemic predictive models to provide an insight into the future infection pattern in a particular region. The epidemic spread models play a vital role in deploying precautionary approaches for containing the pandemic. There are several such models available in the literature such as the SIR model (Kermack and McKendrick, 1927); the Susceptible-Infected-Susceptible (SIS) model (Darabi Sahneh and Scoglio, 2011); the N-intertwined SIS epidemic model (Van Mieghem, 2011), etc. There are other epidemiological models like Susceptible-Infected-Recovered-Susceptible (SIRS), which was first introduced by Anderson and May (1979). Some recent successful implementations of this model have been by Li et al. (2014), Sun et al. (2022), Lan and Yuan (2023), Wang et al. (2023) and many others. Another similar model is the Susceptible-Infected-Recovered-Vaccinated (SIRV), which was also a widely discussed method during the recent pandemic. Some of its recent successful implementations are Liu et al. (2023), Schlickeiser and Kröger (2023), Wen (2022). Such prediction-related studies have also been witnessed where an expression has been presented to forecast the peak time of COVID-19 infections (Kröger et al., 2021). We have come across a study that presents an epidemiological model by integrating Genetic Algorithms and the SIR model to determine the part of the population to be immunized (Rodrigues et al., 2018). The application of statistical machine learning combined with the SIR model to predict the infection progression in India and China is also witnessed in the literature (Das, 2020; Bagal et al., 2020) have used the static SIR model to predict the time for maximum infection in a region. The predictions estimated by the above-mentioned researchers are done to portray a sense of future occurrences. However, using the forecasted estimated values is not further utilized for solving any relief problems during pandemics.

In the context of vaccine distribution, we have encountered literature on the vaccine supply chain using structural equation modeling (Mukherjee et al., 2022). Nevertheless, we have found some literature on vaccine distribution for vaccine-specific distribution models (Goentzel et al., 2022). Cold chains for keeping vaccines cool are also found in the literature (Comes et al., 2018). In recent studies, we have found approaches like system dynamics and regime review (Andiç-Mortan and Gonul Kochan, 2023; Hegde et al., 2023). In the context of other emergency distributions focussed on disaster management, Wan et al. (2023) have solved a multi-period dynamic multi-objective model for emergency material distribution under uncertainty. Peng et al. (2023) have solved a route optimization model for emergency relief material distribution. Wang and Sun (2023) have solved a multiperiod model for sudden large-scale disasters by introducing an equity multiplier. Furthermore, researchers have introduced multimodal modeling of the vaccine supply chain in the healthcare sector in recent studies, which is believed necessary for addressing distribution under emergencies. For example, Enayati et al. (2023) have designed a multimodal vaccine distribution network using drones. Li et al. (2023) have designed and solved a multimodal hub and spoke network for COVID-19 vaccines. Ruan et al. (2016) have modeled multimodal distribution networks for medical essentials in emergencies. Recently, Kar and Jenamani (2024) have proposed a multi-echelon network design for COVID-19 vaccine distribution that considers vehicle selection variables. Azadi et al. (2024) developed a stochastic optimization model for pediatric vaccine distribution in low-income countries, addressing cold chain inefficiencies and optimizing vaccine availability and immunization coverage. To the best of our knowledge, there is almost no litreature on solving the facility location optimization problem addressing the criticality of a region during a pandemic. For health emergencies, we believe all modes of transportation are to be availed for maximum resource utilization. Moreover, we did not come across much litreature that considers the selection of vehicles from multi-modal resources for healthcare supply chain or distribution planning.

In conclusion, we have identified several research gaps in the prioritization techniques used for vaccine distribution during pandemics. The existing HFL problems do not prioritize susceptible regions, and most of them are from the pre-Covid period. Thus, there is a need for prioritization techniques relevant to pandemics with appropriate mathematical foundations. There exist several epidemic predictive models. However, the outcomes of these studies are not used for solving relief problems during pandemics. Additionally, there is a lack of literature that considers multi-modal transportation for HFL problems. These gaps provide an opportunity for our research to address these issues and conduct a study on the optimal distribution of vaccines under the supply of resource scarcity while targeting the susceptible regions.

We follow a two-phase approach to identify the optimal vaccine distribution centers in India. Figure 1 shows the input-process-output model for both phases. Phase 1 uses the modified SIR model, ModSIR, to estimate the RVI to identify vulnerable regions as a proactive measure. Three types of inputs are required for this purpose: infection data, recovery data and population data. Phase 2 uses an aerial/flow path distance matrix, transportation cost, maintenance cost and weights based on estimated RVI, and GPS locations as input for MIQCP formulations. We propose three different weight schemes based on RVI and consequently make three versions of the weighted p-center problem as shown in Figure 1. The results of each variation of MIQCP are sets of hub-and-spoke kinds of networks. To compare the outcome of the MIQCP versions we suggest using the variability of various network properties as described in the subsequent sections. Table 1 presents the lists of the symbols used in the following subsections.

The estimation of the RVI has been done using the outcomes of the static SIR model. The modified SIR model (ModSIR) is a mathematical technique used to describe the spread of infectious diseases as shown in the Algorithm ModSIR shown below. The model assumes that each individual in a population can be classified into one of three categories: susceptible, infected and recovered. The ModSIR model can be used to predict how the infection will spread through a population over time, as well as the expected number of new infections over a given period. The model gets trained using the infection and recovered data for a specific region using the first few days of the infection. Figure 2 compares the working of SIR and ModSIR. The figure shows each of these models has different training data. SIR considers a small period to train itself to estimate the rate of infection (β) and recovery rate (γ). Whereas ModSIR estimates multiple values of these parameters from several periods and considers a combined effect. This approach aligns more with real scenarios. Post-training, predicts the infection pattern for a user-defined period in the future. The model can be used to estimate the minimum time for maximum infection, $Tmin⁡j$ in a region as shown in Figure 2, and the maximum infection percentage $Ijmax⁡$before the peak time of infection. We propose to use these two parameters to estimate the RVI. See the  Appendix for the mathematics of the Modified SIR model. The key difference between the original SIR model and the ModSIR is shown in Table 2.

The entity $Ijmax⁡$indicates the fraction of the population susceptible to infection and $Tmin⁡j$indicates the time. The purpose of RVI is to create a basis for vaccine distribution in a country where the small regions based on administrative boundaries are not truly isolated from each other. Hence, $Ijmax⁡$even if less, can lead to a higher number of infected persons for a region with a high population count. Therefore, we propose a maximum infection rate, $imax⁡j$(equation (1)) to indicate the rapidity with which an infection spreads in a region irrespective of the population size at that location and its normalized value gives RVI as V_j as shown in equation (2). The normalization is done to bring the index value between 0 and 1, where 1 refers to a place that requires immediate attention to vaccine distribution for containing the spread, and 0 refers to a place that requires the least attention:

The Algorithm ModSIR is presented below:

The method to estimate V_j is shown in Algorithm RVI. Following a modified approach presented by Bagal et al. (2020) we make a few assumptions during this estimation. The assumptions are:

• To begin with, there is only a single wave of transmission of an infectious disease.

• Second, one person is infected only once during the wave of infection.

• Finally, the number of infections is assumed as the number of infected people.

Algorithm ModSIR first takes infection and recovered population counts for initializing the model and predicts the spread of the epidemics for each region by giving us two parameters. They are the maximum infection percentage at a location and the minimum time for maximum infection at a location. These parameters are used in Algorithm RVI for estimating the RVI as shown below. This algorithm takes the parameters from the ModSIR model as input and then, it takes each region from the potential facility centers and deploys the SIR model. The outcomes are then individually operated mathematically to get the maximum infection rate. This number shows a wide range in the number line. Therefore, to normalize the number and get the index value between zero to one, we have calculated the numbers as the fraction of the maximum infection rate. The normalized number is the estimated RVI. RVI value zero indicates that the location seeks the least attention from the authority for vaccine distribution, whereas value one indicates a place with immediate attention.

The proposed index provides insights into the vulnerability of different regions. To convert a non-emergency FLP into one that addresses high-vulnerability areas, we have to use these insights. This can be done by integrating the RVI into the objective function of the general FLP by introducing certain RVI-based weights to prioritize regions.

Let us consider an unweighted distance (w₁) as the baseline case for non-prioritized HFLs, over which we want to build a prioritizing FLP. In our literature review, we have come across certain prioritizing FLPs, where techniques where the authors have developed demand-weighted distance to prioritize regions with higher demand (Michael et al., 2013). Let us consider this weight (w₀) as the baseline for the prioritized FLPs. Since the prioritized HFLs are built over the non-prioritized ones, the unweighted distance is named w₁, whereas w₀ is the reference for prioritized ones. By incorporating this weighted distance approach, we aim to enhance the placement of HFLs in a way that better addresses regional vulnerabilities. Building upon w₁, we propose adaptations to these weights as explained in the next section.

To adopt this technique mentioned above (w₀), initially, we conceptualized the weighting scheme, $w2=dijr⋅Vj$ for enhanced coverage of vulnerable regions. This adaptation might work fine while targeting the most vulnerable regions. However, for the problem at hand, this adaptation might deprioritize the less vulnerable regions, which is undesirable for healthcare emergencies. This is because the current less vulnerable regions could be potential hotspots as the timeline progresses. The reason behind this claim is explained in the illustrative example in Illustrative Example 1. Conversely, if we propose another weighting scheme $w3=dijr⋅1/Vj$⁠, which will be an extreme opposite measure of w₂. This is explained in the section Illustrative Example 2. Although these examples are just theories and can be shed light upon by interpreting the results from the computational experiments and numerical results.

For any demand node j, let us consider a weightage scheme represented by $w0=dijr×Dj$⁠. In the context of a minimization model like MIQCP-0, the goal is to establish preferred facilities at locations where the D_j values are maximized, effectively reducing the product, $dijr×Dj$ to zero. However, this approach results in demand nodes with lower D_j values being assigned facilities at distant locations. The reason is that when the coefficients of the objective function are the product of the pair $dijr$ and D_j, the model tends to select a facility location i that minimizes the distance from demand node j. This strategy mitigates the impact of high D_j values. Consequently, facility locations tend to be established around nodes with the most significant D_j values. While this approach is suitable for business or non-emergency applications, it has limitations in healthcare-based scenarios. Based on this example we hypothesize that the weighting scheme $w0=dijr×Dj$ will cause demand nodes with higher D_j values to be prioritized, resulting in lower D_j value nodes being assigned to facilities at distant locations. This is undesirable for healthcare scenarios. For validation of this example, please refer to Section 4.5.

Similarly as stated in Theoretical Illustrative Example 1, in healthcare, using the weighting scheme $w2=dijr×Vj$⁠, the model prioritizes regions based on vulnerability (V_j). However, this deprioritization property becomes problematic when applied to potentially vulnerable nodes, which are susceptible to future virus impacts. In such cases, it is crucial to avoid locating facilities far away from these nodes, ensuring timely access to healthcare services. Based on this example we hypothesize that the weighting scheme $w2=dijr×Vj$ will cause the model to prioritize highly vulnerable regions, potentially deprioritizing less vulnerable nodes that might become critical in the future, which is problematic for healthcare emergencies. For validation of this example, please refer to Section 4.5.

Due to such extreme measures obtained from the use of RVIs in both numerator and denominator in weights w₂ and w₃ respectively, we propose a weight $w4=dijr⋅DjVj$⁠. This nullifies the extreme effects of ignoring the most vulnerable regions. Moreover, as RVI indices are largely driven (not entirely; see equation (1)) by the population of a region, this weight also does not ignore the less vulnerable regions if the population is also used as a weight. This claim can be shed more light upon by observing the numerical results. Although w₄ is an adaptation of w₀ and w₁, w₂ and w₃ are stepping stones for w₄ as explained in the previous two subsections. Thus, to verify these theories, we should compare the results obtained from all the versions of FLPs with different weights. Naturally, the w₀ is the best possible network for non-emergency cases, we need to identify how much we have to compensate to address the vulnerable regions compared to the baseline weights.

In this section, we describe the problem that pertains to finding prioritized facility locations during health emergencies. To set the problem, we have considered a set of potential candidates for facility locations called i ∈ PF = {1, 2,…,|PF|}. Among the |PF| nodes, P nodes are optimally selected, and the associated maintenance cost c_mant_i| i ∈ PF. Moreover, we consider a set of demand nodes j ∈ n = {1, 2,…,|n|}. We assume that total demand is centered at each demand node j ∈ n, and the total demand is denoted as D_j (number of doses). Moreover, for transporting vaccines, we have considered a set of Cold vehicles r ∈ R = {1, 2,…,|R|} where r^th vehicle can be from any mode of transportation and each having a capacity of cap_r liters and a transportation cost of _trans_r $/km. For mathematically connecting the demand and the vehicle capacity, we consider a vaccine dosage of v. liters. That means each vaccine dose is v. liters. Moreover, we have also considered a parameter, $ωijr$⁠, which means if there is necessary infrastructure available for r^th vehicle between i and j. We present a typical non-emergency P-Centre problem as follows:

MIQCP-1

Subject to constraints:

The objective function (3) minimizes the total cost of transportation and facility maintenance. The first part is a nonlinear function explaining the total transportation cost. The second part is linear and explains the facility maintenance cost. This objective does not include w₁, as this weight is a base for introducing weightage. Constraint (4) ensures that each city is assigned to only one facility center. Constraint (5) bounds the number of selected facility centers to P, the number of desired facility centers. Constraint (6) restricts each city from getting allocated to a city other than the optimally located facility centers. Constraint (7) ensures the vehicle selection only if the facility for the vehicle is available between two locations. Constraint (8) explains the demand to be fulfilled. Constraint (9) is used to allocate vaccines only between optimally selected origins and destinations. Constraint (10) ensures a single type of vehicle between two locations. Constraint (11) limits the allocated number of vaccines according to the capacity of the selected type of vehicle. Later the vulnerability index is introduced in the objective function to prioritize the vulnerable regions as described in Table 3. Constraints (12), (13) and (15) force integrality. Constraint (14) shows the domain of vaccine to be transported, and is considered a real number. Since vaccine numbers deal in millions, we can ignore the fractional values. This will reduce the complexity associated with a pure Integer Programming Problem. From this point onwards we replace the weightage scheme w₁ in the objective function (3) with w₀, w₂, w₃ and w₄ will be called MIQCP-0, MIQCP-2, MIQCP-3 and MIQCP-4 respectively. The weight containing the term V_j will convert the problem into a emergency based problem. It should be noted that the constraints in each MIQCP are the same and indicate the following. Therefore, we can demonstrate three distinct versions of MIQCP other than MIQCP-1 by incorporating various weightage schemes.

To evaluate the peaks obtained from the ModSIR algorithm, we compare its results with some baseline models built over SIR. The performance is evaluated based on how close the peaks of these models are to the actual peaks. Thus, we have introduced a few key performance indicators (KPIs) to validate the results. These KPIs include the closeness of peak to actual and percentage improvement of ModSIR peaks over the closeness to actual peaks. The symbols and the interpretations of these KPIs are explained in Table 5.

The results obtained from the proposed formulations show P sets of discrete “hub and spoke” graphs for each MIQCP version. To evaluate their outcomes, we used a few social network metrics and analyzed the variation of the metrics for all the graphs in each MIQCP. The variance is calculated to analyze the equality of the properties of distributed hubs across the country. Variance in the centralities of a set of graphs measures how far their numerical values are from their mean. This is done by calculating the difference between the centrality of each graph and the mean of the set, squaring it, and finally adding up all of the results. The variance will be high if there are extreme values in the centralities of the graph set. Moreover, the distances of the spokes emitting from the hubs of the graphs are also used as performance indicators. The metrics are explained as follows:

Degree Centrality identifies a hub’s importance in a social network (Daskin, 1997). We have applied the statistical variance of this metric to distinguish each hub’s importance in discrete networks from various MIQCP versions. Degree Centrality assigns equal weights to all graph nodes as formulated by Daskin (1997) in Table 4, and its variance interpretation is given in Table 5. To evaluate fairness in coverage, we have introduced weightage vectors comprising the population to be covered by the region. The interpretation of the variance of Population-weighted Degree Centrality is shown in Table 5. Vulnerability-weighted degree centrality, a similar metric, weights hubs by RVI coverage, with interpretation in Table 4. Closeness Centrality gauges node proximity in a social network (Daskin, 1997). For MIQCP, we use hub closeness to assigned locations, with variance indicating proximity variation. Our paper uses average and maximum distances from regions to respective facility centers as performance metrics (Table 5). Along with all the network-based metrics, we propose d_avg_vul≥thrs, which means the average distance of the regions having vulnerability more than a threshold value from their respective facility centers. This metric is extremely important to validate the networks addressing the vulnerable regions.

India has introduced several vaccine distribution drives since its first vaccine in 1804. The government introduced immunization programs like the Universal Immunization Program (UIP) in the year 1885, the Pulse Polio Immunization Programme, Mission Indradhanush, etc. (Chandra and Kumar, 2020; Lahariya, 2014). The latest vaccination program in India was revised in 2014 showing minimal improvement in the coverage, which is not yet updated even under pandemics. Thus, we propose a vaccine distribution strategy focused on pandemics where the distribution criteria are different from normal conditions. However, a large country like India, where different administrative bodies run different geographical areas, requires a systematic approach to vaccine distribution. During a pandemic, a limited number of vaccines must reach the needy regions with high susceptibility rates for limiting infections. The case study is conducted to identify the susceptible regions across the country for vaccine distribution.

This paper has used several data sets to achieve the purpose of solving the FLP. For implementing the ModSIR model, the data sets are collected from a government website (data.gov.in). There are separate data sets for all the states and UTs of India. The data set consists of daily infection, recovery and death data till January 31, 2021. The number of data points varies for each state depending on the infection pattern in a region. However, for a few states and UTs, COVID-19 infection data sets were unavailable. To bridge this gap, we have taken states with similar population densities to interpolate the infection-recovery data. The Integer programming problems require a flow-path distance matrix, an aerial distance matrix and GPS locations. The flow-path distance matrix is calculated using the API Key issued by Bing Maps. Nevertheless, for a few places in India, especially islands, there is no road connection. Thus, the aerial distances are used to link those regions. Moreover, to address the use of all available resources in terms of transportation including air mode of transportation, an aerial distance matrix is calculated. The above-mentioned data sets are enough to achieve the solution to the program described.

The infection data for all the regions had common features in the data set. It had five features namely., date, daily confirmed cases, daily recovered cases, daily deceased, total recovered cases and total confirmed cases. However, since the ModSIR model assumes that the population will remain constant, the daily deceased feature is merged with the recovered column in the pre-processing task. Similarly, features like the “total confirmed” and “total recovered cases” were also removed, as our focus is on daily cases rather than “total confirmed cases”. The distance matrix collected from the Bing Maps produces a square matrix with integer values of distance. We considered 37 potential facility centers for calculating the flow path distance among all the cities forming a 37 × 37 matrix. However, some parts of India are islands that are not connected to mainland India via road. Thus, the distance calculated between the islands of India was the aerial distance between other regions. Among the 36 cities, the model has to choose P number of hubs, which is a user-defined number. In this paper, in the absence of any published literature from GoI, the number P is decided by basing the Indian PDS strategy on Food Corporation of India (FCI). FCI has segregated mainland India into five distinct zones. To match this, we have taken the number P to be five. The data for the parameter related to transportation for the case study of India is given in Table 6.

During the pandemic due to COVID-19, scholars around the world have used the SIR model to predict infection patterns for different countries. The implementation of the SIR model on the infection data sheds light on its behavior when fed with different population sizes and infection data from different regions. However, the scholar community did not use the predicted entities further in any precautionary initiative. With insufficient resources, identifying and addressing the highly vulnerable regions is the challenge. Therefore, we propose an index for identifying the needy regions. This is done by implementing the ModSIR model on the infection data. These predictions are made using the initial infection data for 30 days. The prediction of infection patterns for the next 500 days by analyzing the infection data of the four different locations is shown in Figure 3. In Figure 3(a), the infection pattern for State 1 is explained where the population is 112.4 Million. Within 59 days from the zeroth day of infection, the state records the maximum number of infections with a maximum infection rate of 71%. Whereas, for a smaller and geographically remote state with a population count of just 3.7 Million (State 3), it takes 69 days for the maximum number of infections as explained in Figure 3(c). The maximum infection rate is 82%. Figure 3(b) and (d) explain the infection scenario for State 2 and State 4. Both of these states have similar populations from the northern southern halves of India but differ in the estimated RVI from the SIR model.

To validate the efficacy of the ModSIR model in determining the peak infection time, we compare it with some classical epidemiological models like SIR, SIRS and SIRV models as baseline. The SIRS model, first introduced by Anderson and May (1979), allows for the possibility of re-infection after recovery, adding a parameter for the rate of loss of immunity (δ) in addition to the infection rate (β) and recovery rate (γ) used in the ModSIR model. Recent implementations of the SIRS model include works by Li et al. (2014), Sun et al. (2022), Lan and Yuan (2023) and Wang et al. (2023). The SIRV model is also a widely discussed method in the field of epidemiology. Some of the recent successful implementations of this model are Liu et al. (2023), Schlickeiser and Kröger (2023), Wen (2022). The model requires an additional parameter called the rate of vaccination (ν).

We analyzed the infection patterns for various regions in India, deploying all the baseline models, with data sets ranging from April 2020 to August 2021 representing the first wave of COVID-19 infections. Readers should note that all these epidemiological models have minimal training data, mainly in the early days, so we can predict peak time and take proactive measures. With such small training data, achieving prediction accuracy of day-to-day (vertical axis values) infection is impractical and impossible. Moreover, peak infection time is a very important component of the RVI, and evaluating and comparing them enhances the credibility of the proposed ModSIR.

In the results, we evaluate the performance, defined as how close the predicted peak infection time is to the actual ones, for ModSIR and other models. In “percentage improvement in closeness to actual cases”, ModSIR gives better performance than SIR and SIRV models. For illustration Refer to Table 7. However, the results varied when we extended the ModSIR model to other states. In some states, the ModSIR model performed better in predicting peaks than the SIRS model as detailed in Table 7. The green text indicates, “ModSIR is superior”, while the red text indicates the opposite. Although the SIRS model shows slight improvement in most cases, it requires extensive and rare data sets, which limits its practicality for real-time applications. On the other hand, the ModSIR model, despite being marginally behind SIRS in a few cases, achieves substantial accuracy over SIR and SIRV, with readily available data, demonstrating its robustness and practical applicability. This observation answers the first research question posed in Section 1.

As discussed in Section 3.2, the outcomes of the SIR model are used to calculate the RVI. The estimated RVI for each state and UT is presented in the column named “V_j” in Table 7. The heat map to show the distribution of RVI across India is shown in Figure 4. The numerical values of the RV indices lie between zero and one. As evident from the results, state 28 requires immediate attention for vaccine distribution as its RVI is one indicating the highest infection rate, whereas UT6 requires the least attention as the RVI is estimated to be 0.006 and the infection rate is minimum for this region. Matching the actual vaccination status as of 1.00 PM CET, 3^rd May 2022, state 28 has registered a maximum of 315.5 million vaccines, and UT6 has registered a minimum of 118.5k vaccines (Bing, 2022).

The outcomes obtained from the ModSIR model, after converting to RVI, are used as parameters in the three proposed formulations. The proposed formulations give us a set of P number of hub and spoke graphs for each MIQCP model. In section 3.2, we explain that w₀ and w₁ are the baseline non-emergency-based weights. However, to covert these to emergency-based weights, two weights, namely., w₂ and w₃ were introduced, over which w₄ is built. Thus, for a fair evaluation, MIQCP-4, resulting from w₄, should be compared with all other versions for the following two reasons: First, although w₄ is an adaptation of w₀, for converting a non-prioritized FLP to a prioritized one, w₂ and w₃ are the stepping stones to reach w₄ from the baseline weights. w₄ is built over these two based on some theoretical claims, the proof of which can be interpreted from the results section. Second, although w₀ and w₁ are the baseline weights for prioritized and non-prioritized FLPs, respectively, these weights do not consider any term that reflects the health emergency scenarios. w₂ and w₃ are the emergency-based weights. Thus to evaluate a balanced (between w₂ and w₃) weight w₄, its comparison with emergency-based weights is necessary.

The discrete graphs with P facility centers obtained from each MIQCP are shown in Figure 5. For reference, we have used a heat map representing the population distribution in India in Figure 5(a). As claimed in illustrative Example 1 & illustrative Example 2, we can see that the locations identified from MIQCP-0, are mostly in the regions with the highest population. It is to be mentioned that the MIQCP-1, which is a standard P-Centre formulation, is used as the baseline for comparison. The network generated from this formulation is shown in Figure 5(b). Figure 5(d) depicts that the hubs in the map connect all the vulnerable regions but the facility locations are located far from the allocated regions. Figure 5(e) depicts that the facility locations are located away from the mainland of India. Figure 5(f) depicts that the facility centers are assigned to the cities uniformly in terms of degree centrality as well as high population and vulnerability areas of the country. To explain the effect of the ModSIR parameters (β and γ) on the optimally generated networks, we have presented sensitivity analysis in a supplementary file along with this manuscript. This analysis displays the effect of Phase 1 results on Phase 2. This section summarizes the answer to the second research question posed in Section 1.

For evaluating networks in each set of graphs, we have considered our objective to be the fairness of the result concerning population coverage, vulnerability coverage and closeness to the locations and distances of the locations from hubs. MIQCP-0 is the baseline formulation that minimizes the population-weighted distance and is expected to cover high-population regions with priority (as claimed in Illustrative Example 1). Both Figure 5 and Table 8 prove this remark by positioning its facility locations in highly populated areas and most fairly covering population with a variance in the population coverage of 0.0364 × 10⁵. Since RVI is directly, (not entirely) proportional to the population (see equation (1)), this w₀ corresponding to MIQCP-0 also performs better in vulnerability-weighted degree centrality-based variance. The interpretation of Table 8 summarizes the answer to the third research question posed in Section 1.

MIQCP-1 is the baseline formulation for non-prioritised FLP with an unweighted distance (w₁) is expected to minimize the distance between the demand nodes and their corresponding facility locations. As evident from Table 8, this version excels in average distance, maximum distance from the facility location and total transportation cost. MIQCP-2 is the first version built over the baseline using the technique of Michael et al. (2013) is a vulnerability-weighted distance, expected to prioritize the most vulnerable regions. MIQCP-3 is the exact inverse of the MIQCP-2. From Table 8 it may appear that both these versions of MIQCP are not at par with the other versions when compared with the network-based metrics. However, these versions perform very well when compared based ond_avg_{vul ≥ thrs}.

For the comparison based on d_avg≥thrs, we considered several classes with the top 5, 10, 15, 20, 25, 33 and 36 regions ranked according to their RVI. For each class, the vulnerability of the last region in the group is used as the threshold value (thrs) for vulnerability. The results are summarized in Table 9. For the five, ten and fifteen most vulnerable regions, the average distance from their assigned facility location is the smallest in the network generated by MIQCP-2. For the next three classes, which include the top 20, 25 and 30 regions – meaning less vulnerable regions are also considered – the average distance is minimized in the networks produced by MIQCP-4. Beyond the thirty most vulnerable regions, MIQCP-3 shows the smallest average distance, as these groups include both the most and least vulnerable regions.

MIQCP-4, the balanced version between MIQCP-2 and 3, is expected to prioritize the most vulnerable regions keeping the least ones within a safe reach. To evaluate networks from this version, we can see the network-based KPIs in Table 8. The network from this version excels in the variance of degree centrality, ensuring that the assignment of demand nodes to each facility has been done fairly among all the other versions. MIQCP-4 also excels in closeness centrality. Moreover, the result suggests that, using this network we can address the critical regions with minimum compensation in transportation cost (14.18%), average distance from facility (11%) and maximum distance from the facility (7.78%).

The results can be further interpreted from Figure 6, which is a radar plot. The radius of the plot represents the average distance from the assigned facility location, while the circumference ranks all regions in descending order of vulnerability. Starting from the topmost point and moving clockwise, regions are ranked based on their RVI.

In the case of MIQCP-0, the average distance trajectory starts from zero and increases as we move clockwise, reaching higher radius values (average distance) by the last point. This indicates that this version deprioritizes low-vulnerability regions, as also shown in Table 9. In contrast, MIQCP-1, which does not account for vulnerability, maintains an even average distance from all demand nodes, aiming to minimize overall distances. However, since it does not consider RVI, it places Uttar Pradesh, the most vulnerable state, farthest from the facility location – nearly 600 km away. This also happens with other regions, which could be undesirable during a healthcare emergency.

MIQCP-2 and MIQCP-3 are designed to prioritize the most and least vulnerable regions, respectively. As shown in the figure, MIQCP-2 keeps the most vulnerable regions within an average distance of about 250 km from their respective facility locations. However, as we move towards less vulnerable regions along the circumference, the radius increases to around 500 km. This indicates that similar to MIQCP-0, this version also deprioritizes the lower vulnerable regions, as demonstrated in Illustrative Example 1. Conversely, the distance trajectory in the figure starts by locating high vulnerable regions away from the facility locations (600 km approx.). As we move toward the less vulnerable states along the circumference, the distance trajectory converges toward the center and stabilizes around 300 km.

Among all these versions, MIQCP-4 has been observed to keep a balanced distance compared to all other versions. From start to end, the trajectory keeps all the demand nodes within an average distance of 450 km from the facility locations. Thus, we achieve our objective of prioritizing the high vulnerable regions while keeping the less vulnerable regions within a safe distance from the facility locations.

Validation of Illustrative Example 1:

The results show that MIQCP-0, which uses w₀, positions facilities in highly populated areas and covers the population with a variance of 0.0364, as evident from Figure 6 and Table 8. This model minimizes the population-weighted distance, confirming that facilities are established around nodes with significant D_j values. This supports Illustrative Example 1, as the results indicate that the lower D_j value nodes are assigned to more distant facilities, which is not ideal for healthcare scenarios.

Validation of Illustrative Example 2:

Table 9 and the radar plot in Figure 6 illustrate that MIQCP-2, which uses w₂, prioritizes the most vulnerable regions initially but shows increased average distances for lower vulnerable regions as the vulnerability decreases. This means the model tends to deprioritize less vulnerable regions, validating Illustrative Example 2. The average distance for the most vulnerable regions is the lowest for MIQCP-2, confirming its effectiveness in prioritizing these regions.

In conclusion, while all the versions of MIQCP excel in different KPIs explained above, MIQCP-4 performs in all the areas extremely necessary during a healthcare emergency. For example, in the network-based KPIs, this version performs well in degree centrality-based variance and closeness centrality-based variance. This ensures that the facility center is fairly distributed among all the demand nodes and the nodes are close too, which is extremely necessary during Healthcare Emergencies. Moreover, with minimum compensation, in terms of cost and distance, over a traditional “minimum-cost” formulation (MIQCP-1), we can address all the vulnerable regions. Moreover, as claimed in Illustrative Examples 1 and 2, all other versions with priority-based weights prioritize either high vulnerable regions or low vulnerable regions. Only MIQCP-4 addresses both these classes of regions, which we believe works best in the health emergency. Other evaluations of the networks are explained in the next section.

Although there is no published literature on vaccine distribution strategy acquired by the Government of India, we justify the efficacy of our proposed strategy in two ways. First, by comparing the geographical spread comparing it with the Indian Public Distribution System (PDS); second, by comparing the distribution time and cost based on the data available from various sources. PDS is the only information available for dispatching materials to a large segment of the population maintaining fairness in the distribution strategy. The number of facility centers is estimated from the (PDS) map published by the FCI as shown in Figure 7. GoI has divided mainland India into five different zones geographically, whereas the MIQCP formulations proposed by us are dividing India into different zones optimally. It is to be noted here that, although we compare the segregation of the zones with the PDS (a non-emergency distribution system) strategy acquired by the Indian Government, the objective is not to match the performance of PDS, rather this comparison is conducted to validate that such geographical segregations of India will help achieve fairness. To further analyse these segregations we have individually analyzed the geographical sections obtained from the results of MIQCPs. The separate zones of India are shown in Figure 8 for analysis of the population and vulnerability coverage of each section. Later we interpret the zones and combine them to compare with the zones published by the FCI. The different zones interpreted from the results are explained in Figure 8. The subfigures a–e of Figure 8 explain the zonal segregation based on results MIQCP-0, MIQCP-1, MIQCP-2, MIQCP-3 and MIQCP-4. Geographically, the zonal divisions that come closest to the PDS map (Figure 8) by FCI are MIQCP-4 [Figure 8(d)]. Figure 8 explains the fairness in segregating India based on networks generated from different MIQCPs as compared to the Indian PDS strategy. Figure 9 shows the zone-wise vulnerability and population coverage with all the versions of formulations.

Additionally, we have gathered the information from news bulletins to figure out the actual network followed by the government. Thus, to compare the model performance, we compare the total transportation cost and the total transportation time along the actual network and the prioritized network. In Table 10, we compare the transportation cost and the transportation time for all the zones in India. To calculate the transportation time and the transportation cost we have considered the vehicle types used for the actual vaccine distribution network and the optimized selection of vehicles for the prioritized network. From this comparison, it is evident that the proposed model outperforms the actual strategy obtained by the GoI. The implementation of our two-phase approach can save up to almost 70% of cost and around 76% of transportation time. In Figure 10 we intend to show how simplified the optimized distribution network generated from MIQCP-4 is compared to the actual COVID-19 vaccine distribution network in India. This proves the network efficiency in terms of transportation time, actually answers the fourth research question posed in Section 1.

The results obtained from the ModSIR model demonstrate its efficiency in predicting the peak infection times for various regions, outperforming classical epidemiological models such as SIR and SIRV, especially when considering limited data sets. The ModSIR model shows an improved percentage closeness to actual peak infection times across many states, proving more practical than the SIRS model, which requires rare data sets for higher accuracy. The advantage of ModSIR lies in its ability to work with readily available data, making it robust and well-suited for real-time applications, especially in regions with insufficient resources. Furthermore, the study evaluated the proposed Multi-Objective Integer Quadratically Constrained Program (MIQCP) models designed to optimize the distribution of vaccines based on the RVI. The results showed that MIQCP-4, which balances the prioritization of high- and low-vulnerability regions, provides the most efficient and fair network in terms of population and vulnerability coverage. This model successfully reduces both transportation costs and times compared to the actual network followed during the COVID-19 vaccine distribution in India, with up to 70% savings in cost and 76% savings in time. This suggests that the proposed network model is not only effective in ensuring equitable vaccine distribution but also in optimizing logistical efficiency during healthcare emergencies.

These findings underscore the significance of the ModSIR model in the pandemic response and highlight the practical advantages of optimized vaccine distribution networks in managing public health crises.

This study provides five distinct theoretical contributions by integrating the fields of Epidemiology and Operations Research to address healthcare emergency FLPs. First, we introduce a modified version of the SIR model, referred to as the ModSIR model, to estimate the peak infection rate and peak infection time of regions. These estimations are crucial for prioritizing facility locations during pandemics. Second, we propose a RVI based on the ModSIR model, which quantifies the susceptibility of a region to being severely affected during healthcare emergencies. Third, we develop a non-linear Mixed-Integer Quadratically Constrained Programming (MIQCP) formulation for the P-Centre problem, which incorporates three different weighting schemes. These schemes prioritize vulnerable regions and demonstrate practical applicability in optimizing facility locations during healthcare emergencies. Fourth, we introduce multimodality in the P-Centre FLP to account for vaccine storage and distribution logistics, addressing the complexities of healthcare emergency responses. Fifth, we are the first to validate the fairness in population and vulnerability coverage in FLPs by using various social network and graph-based metrics. These metrics help assess the effectiveness of facility placement, particularly in locating hubs for vaccine distribution.

Though we demonstrate the proposed approach considering the administrative regional divisions in the Indian context in health emergencies, it is equally applicable in other countries when selecting facility locations for vaccine storage and distribution. However, before replicating this study, governing bodies must ensure the availability of data to compute the RVI, distance matrices and transportation costs. Without accurate and comprehensive data, the effectiveness of the model may be compromised, leading to suboptimal facility location decisions. Post data collection; a few key considerations must be addressed when replicating the approach. When implementing Phase 1 of the proposed approach, decision-makers must carefully select an epidemiological model for predicting the infection peak. Unlike models such as SIR, SIRV and SIRS, which assume fixed infection and recovery rates, ModSIR adjusts these rates over multiple periods, capturing the evolving nature of disease transmission. Based on numerical analysis of our use case, ModSIR consistently demonstrates improvements, forecasting infection peaks closer to actual outcomes compared to classical models like SIR and SIRV. However, the performance of the SIRS model is equivalent to ModSIR in many problem instances; it does so, under the assumption that reinfection data is available in every time period. Our experience from infection data collection shows that in reality, the segregation of re-infected populations is indeed difficult. Thus, ModSIR offers an appropriate balance of accuracy and simplicity, while outperforming SIR and SIRV, and providing similar precision to SIRS with fewer data constraints. Thus, it has the potential to become a robust tool for decision-making during health emergencies. Finally, the approach to developing RVI may take different forms based on data availability. Users are advised to normalize these values as presented in this paper and implement the model accordingly.

In Phase 2, several critical points should be addressed. First, the choice of specific weighting schemes should align with the strategic perspective on emergency and non-emergency scenarios. For instance, if the administration aims to establish containment zones during healthcare emergencies, the weighting schemes would differ from those applicable to medical essential distribution. Sensitivity analysis in our study indicates that varying the weights for different parameters (e.g. population density, and infection rates) can significantly affect the selected facility locations. Practitioners are encouraged to conduct similar analyses tailored to their specific requirements. Second, the study provides various metrics to compare the new strategy with existing ones before adoption. This creates a logical basis for replacing the current strategy with the most effective option. Metrics such as cost, degree centrality, closeness to assigned facilities and percentage compensation in cost should be evaluated to ensure that the selected strategy is superior. Third, as the results highlight, the substantial savings of up to 70% in cost and 76% in transportation time compared to traditional strategies. Managers are advised to make sure that the approach performs better than the existing strategy adopted by the governing body.

In addition to these insights, we offer Rule-of-Thumb guidelines for practitioners. First, practitioners should prioritize facility locations in areas with high RVI, while still accounting for less critical regions. This approach ensures that vulnerable areas receive attention while maintaining readiness for outbreaks in lower RVI regions. Instead of solely focusing on cost minimization, managers can aim to build a hub-and-spoke network that effectively serves the most vulnerable areas while minimizing the compensation in cost. Second, practitioners should ensure that the model is adaptable to emergencies beyond epidemiological situations, where different criticality may necessitate estimating other suitable models. Third, practitioners should adjust the model’s weights based on strategic needs during health emergency applications as mentioned in Section 3.2. For example, in crises like fires or earthquakes where only the affected regions matter, vulnerability-weighted distance (Please refer w₂ in Table 3) is ideal. Similarly, during a containment to prevent the virus from gaining a foothold, distance weighted by inverse vulnerability (w₃ in Table 3) works best. Finally, the population-to-vulnerability ratio (w₄ in Table 3) balances coverage and responsiveness, making it well suited for addressing both highly vulnerable regions and less vulnerable areas when necessary while solving the vaccine distribution problem during health emergencies. By incorporating these Rule-of-Thumb suggestions, based on our sensitivity analyses and numerical experiments, practitioners can effectively optimize vaccine storage and distribution in their specific contexts.

This paper presents a two-phase approach to address the four research questions posed at the outset, focusing on healthcare emergencies, epidemiological modeling and facility location optimization. The first phase modifies the classical SIR model to develop the RVI, which significantly improves the estimation of peak infection times across different regions. This modification enhances prediction accuracy, especially in areas with limited data sets, making it suitable for real-time decision-making during healthcare emergencies. In response to the first research question, the classical SIR model was successfully modified to improve peak infection estimation through the creation of the RVI. The second research question was addressed by embedding the RVI into the facility location optimization model, ensuring that vaccine distribution could prioritize the most vulnerable regions. For the third research question, fairness in population and vulnerability coverage was validated using network property-based metrics and KPIs, proving the effectiveness of the proposed distribution strategy. Finally, the fourth research question was answered by introducing multimodality in the P-Centre FLP, enhancing the efficiency of vaccine storage and distribution. The network property-based metrics demonstrated that the ratio of population to RVI effectively balances vulnerability coverage and cost, ensuring fairness in the distribution process. The study showed that introducing multimodality in the P-Centre FLP not only enhances vaccine distribution efficiency but also ensures equitable access across regions. By adapting the model with various weighting schemes, the proposed solution significantly reduces transportation costs and times while maintaining coverage for both high- and low-vulnerability regions. Overall, the study provides a robust and practical solution for managing healthcare emergencies, demonstrating the effectiveness of the modified SIR model in combination with optimized facility location strategies. The findings are relevant for real-world applications, particularly in improving vaccine distribution logistics, as evidenced by the substantial cost and savings in transportation time in the Indian context.

For future work, it is to be noted that the RVI is very specific to a particular wave of infections. In the subsequent wave of infections, the infection rates may be different which will lead to a different set of Vulnerability indices. Thus, we consider the identification of the facility locations from the proposed approach only are temporary locations. However, while the research in the epidemiological field evolves to predict the actual number of waves with a more accurate scenario, the proposed approach can be applied once to address all the vulnerable regions. In the future, researchers can take up the issue, consider the factors affecting accurate predictions and produce a better fit of predicted values to real values. In the future, researchers can connect the subsequent waves to predict the multi-wave scenario of a pandemic. Moreover, if the relevant data is available for all the regions at a national level, this work can be extended using models like SIRS, SIRV and other such epidemiological methods. Regarding the network, optimally located facility locations can be introduced into the vaccine supply chain for vaccine distribution from manufacturers to primary health centers. However, the model in this paper considers the state capitals as potential locations whereas considering all the big cities of the country could have increased the responsiveness to vulnerable regions. Researchers working in the area of the vaccine supply chain can identify those elements and quantify their effects on the criticality of a region to strategize the distribution of vaccines. Furthermore, the result section lacks comparison in terms of facility location maintenance cost, as there is no publicly available data. In case of the availability of this data, the comparison of the model outcome with the maintenance cost of the benchmark network will provide better insights. Apart from emergencies, organizations can use the proposed model to target conventional distribution or service problems by proper identification of the factors affecting their objectives. Moreover, Government organizations can also adopt this technique for public distribution systems or large-scale emergencies by recognizing the suitable criterion.

Initially, this work was funded by iHub Drishti, The Technology Innovation Hub (TIH) under a project named “Rakshak” (IITJ/RAKSHAK/2020–21/01, Dt. 17–09-2020). The authors would also like to take the opportunity to thank the Ministry of Education, Government of India for granting the Prime Minister’s Research Fellowship to one of the authors and for supporting the work.