Fair-split distribution of multi-dose vaccines with prioritized age groups and dynamic demand: The case study of COVID-19

The emergence of the SARS-CoV-2 virus and new viral variations with higher transmission and mortality rates have highlighted the urgency to accelerate vaccination to mitigate the morbidity and mortality of the COVID-19 pandemic. For this purpose, this paper formulates a new multi-vaccine, multi-depot location-inventory-routing problem for vaccine distribution. The proposed model addresses a wide variety of vaccination concerns: prioritizing age groups, fair distribution, multi-dose injection, dynamic demand, etc. To solve large-size instances of the model, we employ a Benders decomposition algorithm with a number of acceleration techniques. To monitor the dynamic demand of vaccines, we propose a new adjusted susceptible-infectious-recovered (SIR) epidemiological model, where infected individuals are tested and quarantined. The solution to the optimal control problem dynamically allocates the vaccine demand to reach the endemic equilibrium point. Finally, to illustrate the applicability and performance of the proposed model and solution approach, the paper reports extensive numerical experiments on a real case study of the vaccination campaign in France. The computational results show that the proposed Benders decomposition algorithm is 12 times faster, and its solutions are, on average, 16% better in terms of quality than the Gurobi solver under a limited CPU time. In terms of vaccination strategies, our results suggest that delaying the recommended time interval between doses of injection by a factor of 1.5 reduces the unmet demand up to 50%. Furthermore, we observed that the mortality is a convex function of fairness and an appropriate level of fairness should be adapted through the vaccination.


Introduction
As of 20 March 2022, the ongoing SARS-CoV-2 virus (COVID-19, hereafter) pandemic has engendered approximately 427 million confirmed cases and 6.09 million fatalities worldwide ( Johns Hopkins University & Medicine, 2021 ). Since the beginning of the pandemic, almost all countries around the world have adopted a variety of measures to restrict the propagation of COVID-19. These measures include individual precautions ranging from wearing face masks to social distancing, non-pharmaceutical measures against COVID-19 to confine its propagation. This challenge becomes even more critical knowing that the number of companies producing COVID-19 vaccines is minimal, and almost all countries have to meet their needs by purchasing vaccines from them. However, countries that can produce vaccines prioritize meeting their own needs first ( Alam et al., 2021 ). This issue forces other countries to be more accurate in estimating their demand when purchasing vaccines. Nonetheless, the precision of such estimation in the COVID-19 pandemic could be jeopardized due to many factors, including the emergence of new variants of the virus, noncompliance with health protocols by some individuals, monopolies in the distribution and consumption of vaccines by developed countries, and insufficient infrastructure and distribution capacity of vaccines ( Alam et al., 2021;Burgos et al., 2021;Tyagi et al., 2021;Wouters et al., 2021 ). Such factors lead to important fluctuations and dynamicity in vaccine demand ( Besiou & Van Wassenhove, 2021 ). Hence, providing a precise estimation of the progression of the pandemic allows more fruitful management of the logistics challenges. In this regard, in order to provide an accurate estimation of disease propagation, a broad range of efficient approaches have been employed in the literature, including artificial intelligence approaches such as time series ( Betcheva et al., 2021;Qi et al., 2020 ), machine learning ( Kushwaha et al., 2020 ), and deep learning ( Yang et al., 2020 ), agent-based approaches ( Kai et al., 2020 ), metrological and meta-population models ( Ma et al., 2021 ), and compartmental epidemiological models ( Arfan et al., 2021 ). From the technical point of view, apart from compartmental epidemiological models, the other approaches hang on data accessibility. Hence, these approaches almost become inapplicable since data is unavailable and unreliable in the early phases of a pandemic ( Gnanvi et al., 2021 ). Also, historical data cannot be exploited due to the occurrence of different virus variants. Moreover, these approaches cannot consider various government strategies, such as quarantine and raising people's awareness ( Biswas & Alfandari, 2022;Nikolopoulos et al., 2021 ). Finally, they cannot optimize the measures related to implementing the policies established by the governments, such as reducing the number of deaths or minimizing system costs to control the propagation of a pandemic. In this regard, the synergy between compartmental epidemiological models and the optimal control problem (OCP) as a powerful dynamic optimization approach not only tackles the barriers mentioned above ( Castilho, 2006 ) but also offers an endemic equilibrium point of the COVID-19 contagious influence ( Gaff & Schaefer, 2009 ). However, due to the complexity of implementing these methods, researchers have rarely used them.
Speaking of which, a careful examination of the COVID-19 pandemic reveals that additional factors should also be addressed in new pandemics when planning vaccine distribution. This paper accounts for four additional factors to design an efficient vaccination network. These factors are classified into prioritizing the population for vaccination , fair distribution of vaccine , multi-dose vaccination planning , and split delivery .
First, since the supply of COVID-19 vaccines has initially been restricted, the world realized that to pilot a vaccination campaign efficiently, governments all over the world should make a tough decision on who to vaccinate first. In fact, health policy-makers all around the world came to a global consensus on vaccination strategies to prioritize the vaccination of front-line workers (e.g., healthcare workers) against the disease, then vulnerable populations such as older people and people with weak immune systems. Accordingly, such prioritization decisions should be taken into account when designing a vaccine distribution network ( Mohammadi et al., 2022 ).
Second, the lack of COVID-19 vaccines raises the issue of a fair distribution of the produced vaccines among or within countries. In fact, vaccine distribution presents a slew of complicated and contentious matters, including public health, finance, public persuasion, and diplomacy ( DeRoo et al., 2020;Forni & Mantovani, 2021 ). Various global organizations, national administrators, and vaccine manufacturers realized that ethics is vital in making decisions. However, inadequate progress has been achieved in defining the fair distribution of vaccines globally and nationally. An "equitable vaccine distribution" is often praised without specifying a structure or offering any suggestions. This paper investigates the fairness in the distribution of vaccines among different vaccination centers, each being responsible to vaccinate the population of their corresponding region.
Third, another factor affecting the operation of vaccine distribution centers is multi-dose vials of vaccines recommended/imposed by the manufacturers to achieve a sufficient level of immunity against COVID-19. Such recommendation/necessity for multi-dose vaccination encounters the distribution center with important challenges. As a matter of fact, multi-dose vaccination imposes a commitment to the distribution network to imperatively fulfill the demand for the second dose, if the first dose has previously been fulfilled. Another point to note is that individuals must receive the second dose of a vaccine after a given time interval, depending on the type of their first-dose vaccine ( Silva-Cayetano et al., 2021 ). In this regard, a number of recent studies have articulated that injecting a second dose from another type of vaccine is possible if it is compatible with the first dose, although the majority of evidence has recommended that it is better if two doses have the same kind ( Sampath et al., 2021 ).
Fourth, the delivery of multi-dose vaccinations puts far more pressure on the distribution network. This becomes even more problematic when we know that COVID-19 vaccines must be transported in refrigerators or freezers, limiting the number of accessible vehicles. One of the most efficient strategies to cope with this obstacle is to use the split delivery method ( Haddad et al., 2018;Mohammadi et al., 2020;Veysmoradi et al., 2018 ).
According to the above-mentioned description, we address the two problems as a comprehensive two-stage framework. The first stage concentrates on estimating the dynamic demand for vaccines and assesses the number of vaccines that need to be provided to control the pandemic. The second stage focuses on offering a dedicated multi-vaccine, multi-depot location-inventory-routing problem (MVMDLIRP) to distribute vaccines among vaccination centers, wherein a broad range of characterful concerns of the COVID-19 vaccine distribution is reflected. The first stage offers an adjusted SIR model with the optimal control problem for controlling the dynamicity of demand for vaccines, and the resulting model provides the number of vaccines needed during the vaccination campaign. The second stage formulates a dedicated MVMDLIRP, where a broad range of decisions are considered, including supply, locationallocation, shipment, inventory, and routing concerning various concerns, including perishability, split delivery, capacity, population priorities, and multi-dose delivery. To the best of our knowledge, the MVMDLIRP has yet to be investigated, let alone the other reflected considerations in the current study that make it more realistic and sophisticated. Hence, the contributions of this study regarding the investigated literature in the next section are as follows: • Introducing a comprehensive and exclusive optimization framework to distribute vaccines to control the pandemic situation, • Proposing an adjusted SIR model with the optimal control problem to overcome the demand dynamicity and determine the number of vaccines that need to be supplied to control the pandemic, and • Formulating the MVMDLIRP regarding a broad range of applicable concerns, including prioritizing age groups, fair distribution, multi-dose injection, split delivery of vaccines, vaccine fraction during transit and vial-opening, vaccine deterioration, and dynamic demand.
The remainder of the paper is organized as follows. Section 2 scrutinizes related literature. Section 3.1 provides the description of the SIR epidemiological model with the optimal control problem. Section 4 presents the MVMDLIRP formulation. Section 5 offers the details of the accelerated Benders decomposition algorithm. Section 6 describes the case study and provides numerical results. Next, Section 7 provides a set of sensitivity analyses and managerial insights. Finally, this study is concluded, and future research directions are proposed in Section 8 .

Literature review
The current research focuses on two main phases: estimating the demand for vaccines through disease progression modeling via epidemiological models and the logistics planning of vaccine distribution. Hence, this section scrutinized the related literature of the two mentioned subjects to illustrate research gaps and our contributions. Asano et al. (2008) suggested a modified Susceptible-Infected-Recovered (SIR) model with the optimal control problem for managing vaccine bait distribution to control the spread of rabies in raccoons, where the rate of vaccine bait distribution was considered as a control variable and determined the number of vaccines needed to control the disease. Nguyen & Carlson (2016) extended the conventional SIR model as a semi-Markovian process to reflect the real-time circumstances of an outbreak for controlling vaccine allocation. Büyüktahtakın et al. (2018) proposed an epidemics-logistics optimization model to control the Ebola virus. This model's decisions include allocating resources regarding capacity and budget restrictions. The objective function minimized the total number of infected people and deaths of infected individuals, where a modified SIR model was proposed to model the Ebola disease transmission. Duijzer et al. (2018b) considered the conventional SIR model to investigate two vaccination strategies, including early and later vaccines. They suggested that a hybrid strategy can diminish the number of infected individuals by more than 50% compared to the best absolute strategy. Duijzer et al. (2018a) employed the conventional SIR model to introduce an analytical model for allocating the necessitated number of vaccines in order to control an epidemic, where the objective function maximized the health benefit, which was defined as the total number of people who escaped from infection with or without vaccination. Enayati & Özaltın (2020) proposed a modified susceptible-exposedinfected-recovered (SEIR) model to manage influenza vaccine distribution by minimizing the number of vaccines used to immunize people and control such an outbreak, where the Gini coefficient reflected the equity concern. Westerink- Duijzer et al. (2020) applied the conventional SIR model to formulate an influenza disease transmission as part of an analytical model to explore sharing vaccines as a redistribution problem, where health agencies might cooperate based on cooperative game theory. Gamchi et al. (2021) proposed a bi-objective mathematical model for influenza vaccine distribution, where the conventional SIR model with the optimal control problem, adapted from Alcaraz & Vargas-De-León (2012) , was employed to prioritize regions and individuals for vaccination. Then a classical vehicle routing problem was used to distribute vaccines. The first objective function minimized the social costs of infected individuals before and after vaccination, and the second objective function minimized the fixed costs of vehicles used in the distribution process. Chen et al. (2014) proposed a mathematical model to formulate an inventory-transshipment problem for vaccine distribution in developing countries. The objective function maximizes the number of immunized children and additional doses of vaccines. The authors consider two types of vaccines (frozen and refrigerated), vaccine fractions during transit and vial-opening, and deterioration of vaccine inventory have been the main contributions. Lim et al. (2022) proposed a multi-modal location-transshipment model to design a distribution network of vaccines, where the objective function minimizes the total costs of opening and operating facilities and transportation. The main feature is the consideration of a replenishment frequency, where two discrete options, including monthly and quarterly, were reckoned. Abbasi et al. (2020) proposed a transshipment-assignment problem to allocate COVID-19 vaccine to people, where the objective function minimized the infection risk of people based on the susceptibility rating of target groups, transshipment time, and the number of over-supplied vaccines.

Phase II-Vaccine distribution problem
They also considered the possibility of transmitting vaccines between vaccination centers, and different target groups to receive vaccines. Larissa et al. (2021) proposed a mathematical model to formulate a classical vehicle routing and scheduling problem, where the objective function minimized the transit time and vehicle and road penalties, which represented unreliable conditions that can jeopardize vaccine distribution. However, none of the specific considerations of vaccine distribution were considered in the suggested model. Yang et al. (2021) presented a multi-period location-routing model to distribute vaccines.
The considered problem accounts for vehicle and facility capacities, and the objective function minimized the costs of opening facilities and transportation. The authors applied the replenishment frequency suggested by Lim et al. (2022) . Rastegar et al. (2021) proposed a multi-period inventory-location model to distribute the influenza vaccine during the COVID-19 pandemic. The model accounts for fairness concerns, where the objective function maximized the minimum delivery-to-demand ratio as an equity measurement. They considered multiple target groups for vaccination without priorities, but a minimum coverage rate for each group was considered as another equitable feature. Tavana et al. (2021) made a similar study to Rastegar et al. (2021) , with the exception that they formulated the COVID-19 vaccine distribution. The main feature was the consideration of various vaccines, which required very cold and ultra-cold refrigeration for equipping distribution centers and their storage. Chen et al. (2021) proposed an agent-based simulation model to allocate COVID-19 vaccines to people with the help of social contact network. They concluded that the proposed simulation model was much more effective than other simulations whose dialectics were based on the number of infected, hospitalized, or dead people.  proposed a bi-objective location-inventory model for a vaccination planning problem, where the first objective minimized the total operational costs and the second one minimized the travel distance of people to vaccination centers. Also, they offered a genetic algorithm and dynamic programming method to solve the proposed model. Gilani & Sahebi (2022) considered a sustainable production-inventory-transshipment vaccine problem as a vaccine supply chain during the COVID-19 pandemic, where the decisions of manufacturers' capacity and capacity expansion, the locations of the packing centers of different vaccines and their distribution centers, and transshipment of vaccines were considered. They proposed a multi-objective mathematical model under uncertain supply to formulate the addressed problem. The first objective function minimized the total costs, the second one minimized the total environmental effects, and the third maximized the number of created jobs. Moreover, they offered a robust data-driven approach to tackle the challenges of the uncertain supply of vaccines. Finally, Table 1 summarizes the significant concerns of the investigated literature to demonstrate research gaps and our contributions.

Epidemiological model
In order to formulate the propagation of contagious disease and calculate the number of susceptible, infected, and recovered individuals in a population, the standard SIR model with the possibility of vaccination, as the well-known epidemiological model, is depicted in Fig. 1 (a). As can be seen, susceptible individuals become either infected once exposed to the disease or recover if vaccinated, where the population size is fixed and predetermined. Notably, the standard model cannot consider the unique concerns of spreading COVID-19 because one of the crucial factors in model-ing the spread of this virus is considering the test results on infected people. So, according to the test results, people must be divided into two groups: quarantined and unquarantined ( Anand et al., 2020 ). Consequently, the adjusted version of the SIR model with the possibility of vaccination is offered to model the spread of COVID-19. The configuration of this model is depicted in Fig. 1 (b), and its notations are presented in Table 2 .
As can be seen, the susceptible individuals ( S) become either infected ( I) once exposed to the disease via the term (β 1 I + β 2 Q + β 3 U ) S, or they get recovered ( R ) if vaccinated via term uS. Unlike the classical SIR model, the adjusted SIR model considers that infected individuals become either quarantined ( Q) or unquarantined ( U) via rates kτ and k (1 − τ ) , respectively. In fact, the Q compartment consists of the infected individuals whose test has been positive. These individuals are quarantined/isolated from the rest of the population, and they no longer jeopardize the health of other susceptible individuals. On the other hand, the U compartment represents a group of individuals who are infected but have not still been tested. Contrary to the Q individuals, the U individuals pose the risk of infection to the susceptible individuals. More importantly, these two compartments are correlated with the testing factor τ . Indeed, the higher the level of this factor, the higher the number of individuals who become quarantined. Hence, there are fewer unquarantined people and less risk of infection for susceptible individuals. Finally, the quarantined and unquarantined individuals are differently removed from the system and become recovered ( R ). It is indeed assumed that the vaccinated individuals will never get infected again. The justified SIR model of Fig. 1 (b) is mathematically formulated as the system of Eq. (1) . Indeed, the system of Eq. (1) represents the balance of input and output time-dependent flows at each state and models its fluctuation and dynamicity.
The basic reproduction ratio R 0 ( Wallinga et al., 2010 ) is defined as the number of new infections caused by a single infectious individual in a completely susceptible population ( Duijzer et al., 2018b ). Accordingly, in the case of R 0 ≥ 1 , the pandemic is anticipated to progress. Next, the method proposed by ( In the following, the equivalent linearized matrices of H (x ) and v (x ) at disease-free equilibrium E 0 = (S 0 , 0 , 0 , 0 , 0) can be constructed as (5) and (6) , respectively.
Finally, the basic reproduction number is equal to the spectral The endemic equilibrium point is then obtained by solving the system of Eq. (7) .

The optimal control problem (OCP)
At this stage, we investigate how intervention measures can mitigate the population's disease burden. For this aim, the system of Eq. (1) is revised to reflect the impact of the control variable u (t) over time. In this regard, the aim is to minimize the cost incurred by infected and quarantined individuals and the immunization program in order to maximize the number of recovered individuals using the possible minimal control variable u (t) . The optimal value of the control variable u (t) (i.e., u * (t) ) can be found through the optimal control problem (8) and (9) . Indeed, this value makes the system follow a trajectory state variable that minimizes the performance measure J, where positive parameters w 1 , w 2 , and w 3 are weight constants balancing the units of the integrand.
The objective function (8) denotes the total incurred cost, where the first term calculates the cost when infected individuals consult medical professionals concerning symptoms and intend to take a test. The second term also represents the cost imposed on the government for quarantined individuals (e.g., social and home-care aids). The third term computes the cost of the immunization program. Noteworthy, the square of the control variable u (t) reflects the severity of vaccination side effects ( Laarbi et al., 2013 ). So, the region for the control intervention u (t) ∈ [0 , 1] is given as final time up to which the control policy is executed, and also u (t) is a measurable and bounded function. The optimal control intervention u * (t) exists in that minimizes the cost function J. Theorem 1. The optimal control intervention u * in of the optimal control problem (8) and (9) In what follows, in order to obtain the optimal control variable u * , the Hamiltonian function is formulated by introducing adjoint variables λ = (λ 1 , λ 2 , . . . , λ 5 ) ∈ R 5 , and minimized by applying Pontryagin's Maximum Principle ( Pontryagin, 1987 ) (please see more explanations in Appendix B ).

Theorem 2.
If u * is the optimal control variable and S * , I * , Q * , U * , R * are optimal state variables of the optimal control problem (8) and (9) , there exist then adjoint variables λ = (λ 1 , λ 2 , . . . , λ 5 ) ∈ R 5 that satisfy the canonical system of Eq. (10) : With the transversality conditions λ i (T f ) = 0 , i = 1 , . . . , 5 , the optimal control variable u * is given as: Finally, we can establish the optimal system by replacing the obtained optimal control variable (11) at the optimal control problem (8) and (9) with optimal state variables. Noteworthy, as far as the pandemic is not over, uS(t ) helps estimate the demand for vaccines for the population that should be provided. More importantly, at the endemic equilibrium point (i.e., R 0 < 1 ), the demand for the vaccines can be estimated as u * S * resulting from the optimal system. Accordingly, the government should increase the value of u until it reaches its optimal value u * to enter the endemic equilibrium point.

Mathematical formulation
This section formulates the MVMDLIRP to distribute purchased vaccines from a central storage hub to a set of vaccination centers through a set of intermediate distribution centers/depots. In this network, the location of the central storage hub and vaccination center is fixed, while the locations of the distribution centers need to be determined. Both distribution and vaccination centers are able to hold an inventory of vaccines in different periods. Moreover, two groups of vaccines are available that differ in their physicochemical properties, including mRNA vaccines (e.g., Pfizer-BioNTech and Moderna) and Adenovirus vector vaccines (e.g., As-traZeneca and Sputnik V). These vaccines require different storage modes while the former should be stored in ultra-cold freezers between [ −80 • C, −15 • C], the latter should be stored in refrigerators between [2 • C, 8 • C]. Therefore, different inventory holding costs are applied to each type of vaccine. It is also considered that a fraction of the vaccines might perish/deteriorate from one period to the next and be lost. To make the problem more realistic, the population is divided into different classes based on their age, and each class has a particular priority to be vaccinated.
According to the description, in order to provide an applicable plan for vaccine distribution, three main categories of decisions should be made simultaneously: design , distribution , and inventory .
The design decision focuses on locating distribution centers from a set of potential locations. Notably, these potential places have specialized storage facilities and usually only need to make small changes in the deployment of equipment. Therefore, when the pandemic spread is controlled, we need to close these centers or reduce their numbers. So, this can be performed with minimal energy, time, and cost, and these facilities can return to their previous activities. Therefore, this level of decision should be reviewed periodically according to the volume of demand for vaccination. Therefore, in the presented model, we have periodically reviewed the decisions related to this level, confirming that the location decision in the current problem does not belong to strategic level decisions. Next, decisions related to the distribution of vaccines should be made. In this regard, the vaccine distribution should be performed from opened distribution centers as intermediate nodes between the central storage hub and vaccination centers. In fact, the inventory levels in opened distribution centers and vaccination centers are the decision criteria for supplying vaccines from the central storage hub and sending them to vaccination centers.
Additionally, the delivery of vaccines from opened distribution centers to vaccination centers is performed via a limited number of refrigerated trucks, and each truck has a limited capacity. At each distribution center, the trucks are loaded with appropriate vaccines, and each truck is responsible for delivering a number of vaccines to different vaccination centers. Consequently, the shipping of vaccines from distribution centers to vaccination centers is addressed as a vehicle routing problem. What is more, due to the dynamic nature of demand, resource limitations, and capacities of trucks and vaccination centers, in most cases, it is impossible to meet the vaccine demand of vaccination centers in one visit. Therefore, in order to overcome this problem, split delivery is intended for the distribution of vaccines. Thus, several trucks from different distribution centers can meet the demand for a vaccination center.
Moreover, since the amount of vaccines produced is much lower than the current demand, a fair distribution system is considered between vaccination centers. For the sake of fairness, a service gap level is defined, which guarantees that the difference between the ratio of satisfied demand to the demand of two different vaccination centers cannot exceed the specified service gap. Therefore, we can guarantee that all vaccination centers can supply the vaccine demand of the people assigned to them within an acceptable boundary. Also, a predetermined percentage of each vaccination center's total demand is guaranteed to be met. In this regard, the shortage of vaccines is also considered by unmet demand that happens when the demand for vaccines is more significant than the supply power. Additionally, the cost of unmet demand is correspondent by morbidity and mortality costs that an unvaccinated individual imposes on the health system. Notably, we have considered this cost to be the same for all target age groups because, on the one hand, this issue is of equal importance to all residents in a community, and on the other hand, it expresses a kind of fairness in the distribution of vaccines.
More importantly, COVID-19 vaccines should be injected in two doses, which should be done in a predefined time interval. Although this time has a lower bound according to the assessment of the World Health Organization and vaccine-producing companies, this lag is not necessarily fixed. In fact, the inventory of purchased vaccines, the severity of the disease's spread, and the government's policies in the vaccination process can change the time interval. In addition, as previously discussed, the second dose should be the same as the first dose, so each vaccination center must respond to two types of demand in each period. The first type, which we call new demand, is related to people who want to inject their first dose. The second type of demand is related to people who want to inject their second dose, and the type of vaccine required for the second injection must be the same as the first dose. Therefore, each period's cumulative demand must be calculated for each distribution center. The main assumptions of the described problem are as follows: • The locations of the central storage hub and vaccination centers are predefined; • The capacities of the central storage hub, distribution centers, vaccination centers, and vehicles are limited; • The shortage of vaccines is allowed, which is considered by unmet demand; • COVID-19 vaccines must be injected twice with a predefined time interval; • The second dose of vaccines must be injected the same as the first dose for each individual; • People are classified into a number of classes based on their ages; • The satisfaction of demand for vaccines is based on the priority of each class of people; • Each vaccine center can be visited by vehicles more than once, which is implied by the split delivery; • The demand for vaccines is uncertain and has a dynamic pattern, which is estimated by the proposed adjusted SIR model with OCP; • Vaccines can become unusable due to fractions during shipping or opening and deterioration; • The unit cost of unmet demand is an estimation of the cost that an unvaccinated individual imposes on the health system (e.g., test, treatment, hospitalization, and death costs); • The government obliges the people to receive the total dose of vaccines in order to receive social services. Table 3 lists the notations to formulate the MVMDLIRP.

The MVMDLIRP formulation
This section formulates the MVMDLIRP in terms of a mixedinteger non-linear programming model using the notations described in Section 4.1 , where the demand for vaccines is estimated from the proposed modified SIR model (9) . The non-linear terms are linearized afterward.

Objective function
The objective function of the MVMDLIRP formulation is presented as (12) which minimizes the total costs of the system.
The objective function (12) includes seven different terms, including 1) fixed cost of opening distribution centers, 2) variable transportation cost of shipping vaccines from the central storage hub to distribution centers, 3) variable transportation cost of shipping vaccines from distribution centers to vaccination centers, 4) variable inventory holding cost of vaccines in distribution centers, 5) variable inventory holding cost of vaccines in vaccination centers, 6) variable transportation cost (per distance) of vehicles, and 7) unmet demand cost. The last term also guarantees the prioritization of different age groups for vaccination. In fact, this term seeks to minimize the total remaining fraction of unmet demand, wherein the fraction of unmet demand of a given vaccination center is weighted by its related priority score.

Constraints
In the following, the constraint body of the proposed MVMDLIRP formulation is explained. The first set of constraints, i.e., (13) - (17) , determine different quantities in the model from distributed vaccines to unmet demands.
Constraint (13) guarantees that the number of vaccines distributed from each distribution center should be less than the number of arrival vaccines to that center, incorporating the shipping lost rate. Similarly, constraint (14) holds the same condition on arrival and distributed vaccines at each vaccination center. Constraints (15) and (16) determine the amount of demand that needs to be fulfilled for the first and the second dose of each vaccine, respectively, such that the fulfilled demand of a vaccine for its first dose in period t should be re-fulfilled in period t + L . Finally, constraint (17) monitors the unmet demands at each period based on the actual demand and the amount of demand that each vaccination center fulfills.
By involving inventory and shipping lost rates, the set of constraints (18) and (19) preserve the inventory level of vaccines in distribution and vaccination centers, respectively.
The set of constraints (20) and (24) ensure the capacity restrictions such that (20) represents the supply capacity for the number of vaccines transferred from the central storage hub; constraints (21) limit the number of vaccines delivered to each vaccination center; constraints (22) ensure the capacity of vehicles when shipping vaccines; and finally, constraints (23) and (24) limit the level of inventory at distribution and vaccination centers, respectively. Notations of the proposed multi-period PISP model.

Sets & Indices
The subtour elimination variable In order to address the fairness in distributing the vaccines among different vaccination centers, constraints (25) and (26) are introduced such that the former determines the service gap, i.e., the maximal ratio between fulfilled demand proportions to the whole demand for all vaccines among all population classes. The latter also signifies that all vaccination centers that are served by the same distribution center will receive the same proportions of their demand.
The remaining set of constraints (27) -(42) are to design the distribution network. Constraint (27) ensures that only opened/established distribution centers can distribute vaccines to vaccination centers. Constraints (28) and (29) , respectively, guarantee that each vaccination center is assigned to at least one distribution center at each period as well as more than a single vaccination center is assigned to each distribution center at each period. Constraints (30) and (31) ensure that shipping vehicles serve only established distribution centers and only the established links between them and vaccination centers. Constraint (32) forces that each vehicle can travel through the link from node i to node j if and only if it is decided that node j be on the path of the vehicle v . Constraint (33) guarantees that each vehicle is assigned to at most a single distribution center. Constraint (34) indicates that each node should be on the path of at least one vehicle.
Constraint (35) states that a vaccination center can be assigned to a distribution center if both centers are on the same route. Constraint (36) allows split delivery and ensures that each vaccination center is visited at least once. That is to say, the number of vaccines required by each vaccination center can be fulfilled not as a whole but through smaller deliveries by different vehicles from distribution centers. Constraint (37) ensures that each vehicle should be dispatched from at most a single distribution center. Constraint (38) forces that each distribution center dispatches at least one vehicle. In addition, constraint (39) ensures that at most one vehicle could be assigned for each route. The connectivity condition for vehicles is guaranteed by constraint (40) , and constraint (41) indicates that distribution centers can deliver vaccines only to the assigned vaccination centers. Constraint (42) The proposed formulation of this section is a non-linear programming model due to the multiplication T p jit Y jit in constraints (13) and (14) as well as the multiplication T p jit ZK jv t in constraint (22) . The linearization has been provided in Appendix D .
The employment of SIR models in a vaccine distribution network is twofold: Apriori or Interactive . In an apriori way ( Gamchi et al., 2021 ), the SIR model is used to estimate the dynamic demand with a given vaccination strategy (i.e., fixed rate of vaccination (i.e., u in 1 )). In this approach, once the demand is estimated for the whole planning horizon, the demand is given as an input parameter to the optimization model. Furthermore, the decisions in one period do not affect the estimated demand of further periods, except the unmet demand accumulated over periods. In an interactive way ( Bertsimas et al., 2022 ), the SIR model interacts with an optimization model in each period. In this way, the vaccination rate u is a decision variable, and the optimization model determines it according to the technical constraints of the distribution network. In this way, the decisions in each period directly affect the estimated demand for further periods, and no unmet demand is considered. Accordingly, the SIR model is updated at the end of each period based on the output of the optimization model in that period, and it estimates the demand for the next period. The new demand is then given to the optimization model for further vaccination decisions. These optimization models, which are executed for each period, are called myopic models. In this paper, we have employed the apriori approach with a given vaccination rate u , where the SIR model is executed once throughout the whole planning horizon.

Benders decomposition -BD
This section employs a Benders decomposition algorithm ( Benders, 1962 ) to solve the proposed model for the MVMDLIRP. Benders (1962) proposed the BD algorithm to deal with problems with complicated variables, where it decentralizes the structure of a problem to provide an easier formation of the problem to reduce the computational burden ( Abou-Ismail, 2020; Rahmaniani et al., 2018 ). From then until now, it has been gaining prestige as one of the most powerful exact algorithms to solve a broad range of NP-hard problems, such as vehicle routing ( Corréa et al., 2007 ), facility location ( Boland et al., 2016 ), logistics network design ( Cordeau et al., 2006 ), transportation ( Coelho & Laporte, 2014 ), and inventory vehicle routing ( Alkaabneh et al., 2020 ). What is more, many studies in the literature have used the BD algorithm to solve the vehicle routing problem and its derivatives, such as location-routing ( Çalık et al., 2021 ), location-inventory-routing ( Zheng et al., 2019 ), and production-inventory-routing ( Cordeau et al., 2015 ) problems. As a result, several efficient valid inequalities have been presented for these problems that could limit the solution space and improve the convergence speed of the BD. Hence, since the proposed model in the current research is one of the derivatives of the vehicle routing problem, the BD algorithm is employed to solve the proposed model.

BD in general
The idea behind the BD algorithm is to divide the original problem into a master problem and a set of sub-problems, with the hope that the decomposed problem be easier to be solved ( Abou-Ismail, 2020; Rahmaniani et al., 2018 ). The master problem incorporates a part of the original model with only integer variables. On the other hand, sub-problems are formed by applying programming duality over the rest of the original model, knowing that the value of the integer variables is given. Along with the integer variables in the master problem, an artificial variable is considered that describes a lower bound (upper bound) on the sub-problems' objective function for a minimization (maximization) problem. Through the BD algorithm (see Fig. 2 ), master and sub-problems are solved iteratively, such that the master problem is solved first and the values of the integer variables are determined. Next, the sub-problems are solved for the given value of the integer variables from the master problem. Finally, a feasibility cut (based on the values of the sub-problems' variables) is added to the master problem if some sub-problems are infeasible or unbounded; else, an optimality cut is added. If sub-problems are feasible, an upper bound can be obtained, and if the optimal solution is obtained by solving the master problem, a lower bound can be obtained. This procedure is repeated for further iterations until a stopping criterion is met (i.e., the maximum number of iterations or a threshold for the gap between the obtained lower and upper bounds). Afterward, a set of accelerators, in terms of inequalities , are employed to enhance BD's performance, the so-called accelerated BD (ABD) algorithm. In the following subsection, a set of valid inequalities are introduced.

Valid inequalities
The first set of valid inequalities (44) -(46) have been widely used in the literature for inventory and location routing problems ( Coelho & Laporte, 2014;Darvish et al., 2019 ). These inequalities are employed to empower the linear programming relaxation of the MVMDLIRP formulation.
The idea behind the valid inequality (44) relates to whether the inventory held at each vaccination center at each period is sufficient to fulfill future demands. More precisely, if the inventory held in period t 1 by vaccination center i is sufficient to fulfill its demand for periods [ t 1 , t 2 ] , then no delivery of vaccines to vaccination center i is required. Indeed, if p∈ P I VC On the other hand, if the inventory is not sufficient to fulfill future demands, then a delivery must take place. Moreover, the valid inequality (45) eliminates links between a node and itself. Finally, the valid inequality (46) forces that one route can start and end at a distribution center in each period if the distribution center is opened/established. Another type of valid inequality, called knapsack inequality ( Santoso et al., 2005 ), is used in the proposed ABD algorithm. The knapsack inequality is added to the master problem with the aim of accelerating the branch-and-bound process of the solver. Indeed, the convergence speed of Benders decomposition can be improved by adding the knapsack inequality, which assists progressive solvers like Gurobi in deriving a range of valid inequalities.
where LB n represents the best lower bound found so far (i.e., the objective function value (the best lower bound found by the solver) in the previous iteration in case of feasibility (infeasibility)), and ζ is an additional variable representing the total cost except fixed opening cost distribution centers as well as the variable shipping cost of the vehicles.
The Benders reformulation of the MVMDLIRP is provided in Appendix E .

Computational experiments
This section presents the computational results of the designed vaccine distribution network under demand dynamicity. In the following, Section 6.1 designs a set of experiments based on a real case study, the vaccination campaign during the COVID-19 pandemic in France. Next, Section 6.2 presents a comparative analysis between the proposed ABD algorithm and the Gorubi solver (GRB). Finally, a set of comprehensive sensitivity analyses are conducted in Section 7 to investigate the behavior of the proposed MVMDLIRP model to any changes of its input parameters. These sensitivity analyses help to draw valuable insights for health policymakers ( Fig. 3 ).

Experimental design and case study
This section describes a real case study (see Fig. 3 ) on the vaccination campaign in France. France encompasses a total number of 18 regions (i.e., 13 metropolitan and 5 overseas regions) and is divided into 101 departments (i.e., 96 metropolitan and 5 overseas departments). Among all, a number of 12 metropolitan regions (see Fig. 3 (a) with the total eligible-for-vaccination population of each region in parentheses) are considered to be the potential location for establishing distribution centers (i.e., | J| = 12 ), and a number of 80 metropolitan departments are considered as the vaccination centers (i.e., | I| = 80 ) or the areas that have a particular demand for vaccines. The metropolitan regions in France possess different properties (e.g., demographic, educational, social, economic, etc.), and these disparities among regions make vaccination a complex work for France's health policymakers. A summary of these properties has been provided per region in Table F.1 in Appendix F . In this case study, the population in each area is divided into four age-based classes (i.e., | K| = 4 ; ages: "18-49", "50-64", "65-74", "> 75 "), and each area includes a different number of individuals per each age-based class (see Fig. 3 (b)). The central storage hub of the purchased vaccines is located in the capital (i.e., region de-France). A total number of 24 refrigerated trucks are considered (i.e., | V | = 24 ) to ship the vaccines from distribution centers to vaccination centers.
The dynamic demand for vaccines in each populated area is obtained by running the proposed system dynamics model (9) for that area, where the birth and death rates as well as the population are different for each area. Detailed information on all other parameters of the adjusted SIR model and the MVMDLIRP model is provided in Table 4 . Except for the birth and death rates as well as the population in each set of populated areas, other parameters of the adjusted SIR model are considered similar for the whole of France ( Angeli et al., 2021;OPECST, 2021 ).
In Table 4 , the unit shipment costs c and cs are calculated based on the rent price of refrigerated trucks divided by their holding capacity (per pallet of 10k vaccines). Depending on the type and capacity of the truck, this cost varied from 10 to 40 euros. In addition, the variable (per distance) shipment cost π represents the vehicles' traveling cost (e.g., fuel). The unit holding costs h DC/VC are also proportional to the shipment cost of vaccines, c, since the shipment cost is indeed the cost of utilizing mobile refrigerators to hold vaccines. Furthermore, the unit unmet demand cost is calculated as the expected cost that an unvaccinated individual imposes on the health system (e.g., treatments, medicines, hospitalization, test costs, etc.). Figure 4 illustrates the outcomes of the adjusted SIR model in determining the state variables for two distinct populated areas, the demand points of vaccines. As can be observed, the dynamicity of the state variables in the two areas are significantly different, and they do not follow a specific distribution. These issues demonstrate the necessity of employing epidemiological-based system dynamics models (i.e., the adjusted SIR model) to determine the demand for vaccines. Remarkably, since the pandemic is not over at the time of this research, the state variable uS ( Fig. 4 f) is used to estimate the demand for vaccines in each populated area.
The performance of the proposed ABD algorithm is tested through a set of 36 test problems with a different number of vaccination centers and time periods. All test problems except the last one are indeed smaller parts of the case study with fewer vaccination centers or fewer time periods. Since the adjusted SIR model is run for each populated area for a time horizon of 365 days, and it is also impossible and illogical to plan the MVMDLIRP daily, each period of the MVMDLIRP model represents a week and the demand for each period is the cumulative demand of the whole week.
Both the MVMDLIRP model and the ABD algorithm were coded in Python 3.8 using the Gurobi 9.1.2 library, and all experiments were done on a server containing four Intel XEON processors with 5 gigabytes of RAM memory running at 2.3 gigahertz. Furthermore, two stopping criteria were considered, when executing the ABD algorithm and GRB, as 1) a gap of 1% and 2) a CPU time of 7200 seconds. For the ABD algorithm, the first criterion is the gap (%) between the obtained lower-bound and upper-bound at each algorithmic iteration; and for GRB, it is the gap (%) between the bestfound solution and the current obtained solution.

Numerical results
This section conducts a comparative analysis of the performance of the proposed ABD algorithm with valid inequalities (ABD-w-VI) and without valid inequalities (ABD-wo-VI) compared to GRB. Table 5 shows the results of this comparison for 36 test problems with a different number of vaccination centers (i.e., column | I| ) and different periods (i.e., column "| T | ") in terms of objective function values (i.e., columns "Obj. Values") and computational time (i.e., columns "Time (s)"). In columns "Obj. Values", the values are proportional to the objective value of the test problem #1. The goal of this comparison is to evaluate the performance of the proposed ABD algorithm compared to GRB and evaluate the benefits of valid inequalities for the proposed ABD algorithm.
As can be seen in Table 5 , both ABD-w-VI and ABD-wo-VI algorithms have been able to obtain solutions with a gap of less than 1% for all test problems before reaching the maximum allowable CPU time of 7200 seconds, even for larger test problems. However, GRB has been unable to find the optimal solution for test prob- (70 , 14 − 18) , (80 , 8 − 12) } under a computational effort of 7200 seconds. Furthermore, GRB has not even been able to find a feasible solution for larger test problems (80 , 14 − 24) } under the limited computational time. Moreover, columns "Gap" compare the performance of the ABD-wo-VI algorithm and GRB with the ABD-w-VI algorithm in terms of objective function gap (%) and the computational time ratio. Regarding the objective function gap, we only compare the performance of GRB, since the ABD-wo-VI has reached the same solutions as those of the ABD-w-VI algorithm. In terms of objective function gap, the ABD-w-VI algorithm has obtained much bet- ter solutions than GRB with a mean gap of 16% over test problems for which GRB has reached at least a feasible solution. In terms of computational time, a paired comparative analysis has only been done over the test problems for which both corresponding methods have obtained the optimal solution with a gap of less than 1% within the computational time of 7200 seconds. It can be observed that the ABD-w-VI algorithm is faster than both the ABD-wo-VI algorithm and GRB for all the corresponding test problems with mean ratios of 1.75 and 12, indicating that the ABD-w-VI algorithm is, on average, 1.752 and 12 times faster than the ABD-wo-VI algorithm and GRB for solving the MVMDLIRP model. To better show the difference between algorithms, Fig. 5 illustrates the comparison between three solution techniques (i.e., ABD-w/wo-VI and GRB) in terms of both objective function values ( Fig. 5 (a)) and the computational time ( Fig. 5 (b)).
In order to explore the detail of results, Table 6 reports different information for all test problems, including the number of opened distribution centers ("OCD" in %), the unmet demand ("UMD" in %) for each age-based class, and the contribution of each vaccine to the vaccination of each age-based class ("CVG" in %).
Looking at column "ODC (%)", we observe that the larger the size of the problem (i.e., the number of vaccination centers and planning periods), the higher the number of distribution centers to fulfill the demand. Overall, about 82% of the total number of distribution centers have been used to fulfill the demand in the distribution network. The reason for the value "0" for instances with | T | = 2 (i.e., two weeks) is that the initial inventory of vaccines in the vaccination centers is sufficient to fulfill the demand for such a short planning horizon. When solving the original case study with | I| = 80 for planning horizons of | T | ≥ 8 (i.e., greater than 8 weeks), the whole capacity of the distribution network is utilized by opening 100% of the distribution centers.
The results of columns "UMD (%)" indicate that, in overall, 6.4%, 9.4%, 13.8%, and 23.1% of the demand have remained unmet for age-based classes "> 75 ", "65-74", "50-64", and "18-49", respectively. Comparing age-based classes shows that a lower percentage of the population remains unvaccinated, from elders to younger individuals. This comes from the prioritization of the age-based classes in the objective function. More importantly, it can be seen that the unmet demand of elders remains fixed when increasing the number of vaccination centers (or indirectly increasing the potential demand for vaccination); however, the unmet demand of younger individuals increases. The reason is that no matter how much the population is, the system put its all efforts into vac-   cinating the elders as many as possible. The rest of the results in columns "CVG (%)" show that Pfizer (Pf) is always higher than Moderna (MO) and AstraZeneca (AZ) vaccines. The main reason is the higher quantity of Pfizer vaccines purchased and stored by the French government. These results of Table 6 are helpful to extract some rules for the government when designing a vaccine distribution network. For instance, when planning to vaccinate eligible individuals through a set of | I| = 80 vaccination centers for a horizon of | T | = 24 weeks, the optimal vaccination strategy would be opening the whole distribution centers and allocating, on average, 74% (average over 67. 3, 73.1, 77.0, and 78.4), 17.5% (average over 13.7, 23.5, 18.0, and 14.8), and 8.5% (average over 19.0, 3.5, 5.0, and 6.7) of the available capacity (purchased or stored) of vaccines Pfizer, Moderna, and AstraZeneca vaccines to vaccinate the population.

Sensitivity analyses
This section provides a comprehensive sensitivity analysis of three main criteria, including the total cost of the system, the total inventory of different vaccines, and the total unmet demand of each age-based class with respect to changes of certain input parameters, including the minimum percentage of demand to be met ( β), the maximum supply of vaccines at vaccination centers ( Q VC ), inventory holding capacity at vaccination centers ( Q VC ), inventory holding cost at vaccination centers ( h VC ), the demand of vaccines ( d), the time interval between two doses injection ( L ), fairness service level gap ( δ), and unmet cost ( cu ). Figures 6-13 illustrate these sensitivities, wherein the value "1.00" on x -axis represents the original level of the corresponding parameter in the case study and other values on x -axis are proportional to the original value. Figure 6 illustrates how altering parameter β affects the three mentioned criteria. Precisely, Fig. 6 (a) shows that increasing β forces the system to vaccinate more individuals, and consequently, the corresponding costs also increase linearly. Figure 6 (b) illustrates that the level of inventory for all three vaccines decreases since a larger demand should be fulfilled, and consequently, a larger portion of the inventory should be spent. Since a higher amount of inventory belongs to the Pfizer vaccine, its decrease happens with a higher slope. Finally, Fig. 6 (c) illustrates that imposing the system to fulfill more demands causes more unmet demand for all age-based classes. The reason goes back to the limited availability of vaccines in the distribution system. Indeed, the system's capacity to fulfill extra demand is limited, and the extra demand remains unmet. Figure 7 investigates whether changing the supply capacity of vaccines in vaccination centers impacts the concerned three criteria. Figure 7 (a) shows that the total cost of the system decreases as more vaccines are allowed to be supplied in the vaccination center. The potential reason for this decrease is twofold, the reduction in the total inventory holding cost or the reduction in the unmet demand cost. As Fig. 7 (b) illustrates, the level of inventory increases, while Fig. 7 (c) shows that the unmet demand decreases.  What can be inferred is that the reduction of the unmet demand cost outweighs the increase of the inventory holding cost. Consequently, the total cost of the system decreases. The trends in Fig. 7 (b) and (c) were expected since any increase in the supply capacity of the vaccines permits to store of more vaccines (i.e., an increase of the inventory level), and much larger storage helps to vaccinate more individuals (i.e., less unmet demand). What can be extracted more from Fig. 7 (c) is that the unmet demand of age-based classes with lower priority (i.e., "18-49") always remains higher than the unmet demand of high-priority classes (i.e., "> 75 ").
Moreover, Fig. 8 investigates whether any change in the storage capacity of vaccination centers affects the performance of the distribution system. As can be seen, by increasing the storage capacity, all three concerned criteria show a decreasing trend.
Importantly, significant decreases (e.g., 80% for age-based class "18-49") happen in the unmet demand for all age-based classes when the capacity of the vaccination centers increases up to 40% (from 0.6 to 1.0 on x -axis of Fig. 8 (c)). Indeed, the increase of the capacity could be either increasing the capacity of the current vaccination centers by adding more vaccination lines or even opening new vaccination centers. It can also be observed that increasing the current capacity of the vaccination campaign in France up to 80% (point 1.8 on x -axis) will cause important reductions in the unmet demand in all age-based classes; however, any increase higher than 80% will not affect the unmet demands. Similar to Fig. 7 (c), the unmet demand of age-based class "18-49" always remains higher than other classes; however, a significant decrease happens on the unmet of this class when increasing the capacity.
Since the vaccines require special modes of storage (i.e., ultracold freezers/refrigerators), the inventory holding cost plays a vital role in the distribution centers' performance. In this regard, Fig. 9 illustrates how the system reacts to the increase in the inventory holding cost. As can be seen in Fig. 9 (a), the total cost   9. Impact of inventory holding cost on total cost, inventory, and unmet demand.
of the system increases with the increase of the holding cost. When this parameter increases, the reaction of the system is to absorb such an increase by decreasing the inventory level ( Fig. 9 (b)), with the hope to reduce the total inventory holding cost. On the other hand, any reduction in the inventory level of vaccines signifies less capability of the system to vaccinate the population. Accordingly, the unmet demand increases as illustrated in Fig. 9 (c). In such a situation, the increase of the total unmet demand cost prevails to the reduction of the total holding inventory cost; hence, the system's total cost increases. In addition, Fig. 9 (b) reveals that by increasing the inventory holding cost, the system puts its effort into reducing the share of Pfizer and Moderna since they require ultra-cold freezers and possess a higher holding cost. However, the inventory level of the AstraZeneca vaccine remains stable since it costs less in terms of inventory holding. Figure 10 depicts the impact of demand alteration on the performance of the system and illustrates how the increase in demand affects the cost of the system as well as the levels of inventory and unmet demand. Figure 10 (a) represents the reaction of the objective function to the increase in demand. As can be seen, this reaction consists of two phases: a decreasing trend followed by an increasing trend in the total cost of the system. In the decreasing phase, the demand for vaccines is small, and it is mostly fulfilled via the initial inventory of the vaccines in both distribution and vaccination centers. During the time that the system consumes the initial inventory in this phase, no extra vaccines are neither distributed nor stored in the network. Hence, both inventory holding and transportation costs are reduced, leading to a reduction in the system's total cost. On the other hand, once the system runs out of the initial inventory, the system starts distributing and storing new vaccines. This phenomenon then increases the total cost of the sys-  tem. In Fig. 10 (b), the decrease in the inventory level of Pfizer and Moderna is higher than AstraZeneca, since a more significant initial inventory has been considered for the two former vaccines. Finally, Fig. 10 (c) shows an expected trend for the unmet demand, as any increase in the demand for vaccines increases the potential unmet demand due to the limited availability of the vaccines in the system.
A more interesting result has been obtained when investigating the impact of the time interval between two doses of injection of the vaccines on the performance of the system, as depicted in Fig. 11 . Figure 11 (a) illustrates that if the time interval is reduced to half of their recommended delay, the unmet demand increases by a factor of 2.5 as shown in Fig. 11 (c). Indeed, when the system focuses on a complete vaccination with a short time interval, more individuals are completely vaccinated, but a considerable part of the population remains un-vaccinated (i.e., unmet demand).
Furthermore, if the recommended time interval increases by a factor of 1.5, the unmet demand decreases up to 50%. Increasing the time interval signifies that more individuals receive at least one dose of vaccines. On the other hand, this increase will not always reduce the unmet demand since, in further periods, the accumulated demand for the second dose of vaccines in early periods is summed up with new demands for the first dose; hence, augmented demand for vaccines is imposed on the system. Therefore, the unmet demand starts to be increasing when the system is disabled to absorb such an augmented increase in the demand. Accordingly, Fig. 11 (a) and (b), collectively, show that the optimal performance of the system happens when delaying the time interval between two doses of injection of vaccines with a factor of 1.5. Figure 12 depicts the impact of the fairness service level on the unmet demand. Lower values of this service level guarantee fair access for the population among different regions to vaccines. On the contrary, higher values of this service level allow a higher dispersion of vaccine allocation among different regions. The former situation should provide a higher equity and higher satisfaction among the whole population; however, the latter puts more effort into al-  locating vaccines to high-populated regions with a higher risk of infection due to higher social contacts. Figure 12 illustrates that increasing fairness (i.e., decreasing δ) is not always the best strategy in vaccination campaigns to stop pandemics. In fact, an optimum level of fairness leads to the minimum level of unmet demand. More explanations have been provided through managerial insights in Section 7.2 .
In another sensitivity analysis, we investigate the impact of unmet demand cost cu on the objective function (12) ( Fig. 13 (a)), the total distribution cost (i.e., the objective function (12) , except the last term), and the total amount of unmet demand (i.e., the last part of the objective function (12) ). As can be seen from Fig. 13 (a), increasing the unit of unmet demand cost increases the objective function, since the network aims to vaccinate as many individuals as possible to alleviate the impact of unmet demand. Figure 13 (b) shows, in detail, how increasing the unit unmet demand cost increases (decreases) the total distribution cost (amount of unmet demand).

Discussion and managerial insights
In this section, some of the more interesting results are discussed, and a sort of managerial insights are provided.
We observed that if the recommended time interval increases by a specific factor, the unmet demand decreases significantly. It means that vaccinating a higher number of individuals, even with a single dose, remarkably decreases the death rate. However, it implies vaccines with a high level of efficacy for their first-dose injection. Therefore, the availability of vaccines also plays a vital role in piloting a vaccination campaign. In reality, depending on the availability of vaccines, several countries have recommended an increase in time intervals between two injections ( ECD, 2021 ). In this regard, by June 2021, 16 EU countries increased the time interval to provide more individuals with their first vaccination dose. Romania applied for such extension only in particular cases Dascalu et al. (2021) ; however, seven other EU countries (Iceland, Latvia, Lithuania, Malta, Slovakia, Slovenia, and Spain) did not apply such a strategy. Shortly after, seven countries removed this extension in September 2021 (Austria, Belgium, Czechia, Finland, Ireland, Portugal, and Spain). The possibility of such an increase in the interval between two injections introduced the concept of fractional dosing . Actually, due to the limited availability of vaccines, especially at the beginning of the pandemic, fractional dosing could be a crucial decision to make a trade-off between reducing the number of potential deaths (by vaccinating individuals with complete doses) or decelerating the pandemic (by providing a higher number of individuals with the first dose).
Scientists have always debated the impact of equity and fairness in pandemics ( Enayati & Özaltın, 2020;Mohammadi et al., 2022 ). Although a higher level of fairness in vaccine allocation provides more social satisfaction among the population, the fair distribution of vaccines is not always the best barrier strategy to stop pandemics. In this paper, we observed that the unmet demand is a convex function of fairness with a global minimum. Therefore, low and high levels of fairness do not represent the optimum allocation strategy. In fact, a low level of fairness puts more focus on the vaccination of populated areas with a higher risk of infection while ignoring less populated areas where the virus may transmit without sufficient barriers (e.g., vaccination). On the contrary, a high level of fairness guarantees fair access to the population to vaccines. However, there would be an imbalance between high and low-populated areas regarding barrier actions. Therefore, there should always be a trade-off between effort s f or fair access to vaccination and attempts to control/stop the pandemic in highrisk (populated) areas. The focus on populated areas reveals the importance of the level of social contact between individuals during the progress of the pandemic. Therefore, from an optimization viewpoint, the effort of health policymakers should not only be vaccinating high-risk individuals to decrease deaths but also vaccinating the population with higher social contacts (i.e., populated areas). Although the former may directly decrease the number of deaths among the population, the latter directly breaks/decreases the circulation of the virus in society and indirectly reduces the number of deaths among the population. In this regard, one way to control social contact could be imposing quarantine restrictions in some regions while vaccinating others.
An important practicality of the proposed model for health policymakers can be highlighted from the results shown in Fig. 13 (b). In fact, parameter cu takes the role of an interplay between two viewpoints, cost-oriented or mortality-oriented vaccination strategies. As a matter of fact, increasing this parameter drives the proposed model from a cost-oriented network toward a more mortality-oriented network. Accordingly, although the objective function of the proposed model is apparently a business objective, adjusting the unit unmet demand cost can satisfy mortalityrelated concerns. In this regard, a possibility could be separating the two parts of the objective function (i.e., distribution cost and amount of unmet demand) and optimizing them through a biobjective model. The gain of such a model could be generating a set of non-dominated solutions and providing health policymakers with a portfolio of different vaccination strategies.

Conclusions
The most effective intervention strategy to control the COVID-19 pandemic, which requires careful preparation due to its availability, is vaccination, the principal hope of today's communities. In this regard, we proposed a novel exhaustive scheme that contains three phases of demand regulating, location-inventoryrouting planning, and solution approach, wherein a broad range of concerns of the COVID-19 vaccine procurement and distribution have been considered. In the first stage, we proposed a modified SIR epidemiological model for handling the dynamicity of the vaccine demands, where the optimal control problem was involved in determining the number of vaccines that need to be purchased to control the pandemic situation. In the second stage, we have presented and formulated a multi-period, multi-vaccine, multi-depot location-inventory-routing problem in the form of a mixed-integer programming model, where various realistic concerns of the COVID-19 vaccine distribution, ranging from multidose vaccine regimen to prioritizing age groups, have been reflected. Finally, we offered an accelerated Benders decomposition algorithm to solve sizable test problems in the third stage. Although the primary aim of the current research was to cope with the concerns of the COVID-19 vaccine procurement and distribution, the soundness outcomes revealed that it could indeed be applied to cope with various vaccine distribution problems. In fact, due to the complexity of the considerations that had to be re-flected in the distribution planning of the COVID-19 vaccines, the proposed optimization model factored in almost all essential concerns in the vaccine distribution. As a result, it could be fruitful to distribute existing or forthcoming vaccines. The computational results demonstrated that the proposed framework could provide a practical plan for the preparation and distribution of vaccines. In this regard, we provided operational planning to prioritize vaccination target groups, and we were also able to schedule a time interval between injections of different vaccines if possible to control the pandemic better.
For future research, we can mention a number of interesting cases that can be considered in both levels of the presented framework. Regarding the first level of the framework (i.e., the adjusted SIR model), considering the efficacy of various vaccines would be interesting to be considered in the model. Various facets of this problem make it worthwhile to study. If immunizations are only moderately effective, the booster dose should be administered sooner. Furthermore, if the vaccine's efficiency is excellent, the option of delaying the second dosage of the vaccine rises, allowing a more significant percentage of the target population to get vaccinated. The second involves anti-vaccination individuals. These people increase the possibility of infection and raise mistrust among others who have not yet been vaccinated. Regarding the second phase of the framework (i.e., vaccine distribution network), the proposed optimization model can be integrated with the adjusted SIR model in an interactive way to develop a myopic model. In this model, the SIR model and the optimization model interact in each period, and the vaccination rate u in Section 1 is considered as a decision variable, and the vaccination decisions in each period directly affect the estimated demand of the further periods.

Appendix A. Proof of Theorem 1
Proof. The optimal control intervention exists when the following conditions are satisfied: • The solution space of system dynamics (1) with control variable u in = φ.
• The mentioned set is closed, and convex, and the state system is represented with a linear function of the control variable where coefficients depend on time and also on state variables.