Equalizing access to pandemic influenza vaccines through optimal allocation to public health distribution points

Vaccines are arguably the most important means of pandemic influenza mitigation. However, as during the 2009 H1N1 pandemic, mass immunization with an effective vaccine may not begin until a pandemic is well underway. In the U.S., state-level public health agencies are responsible for quickly and fairly allocating vaccines as they become available to populations prioritized to receive vaccines. Allocation decisions can be ethically and logistically complex, given several vaccine types in limited and uncertain supply and given competing priority groups with distinct risk profiles and vaccine acceptabilities. We introduce a model for optimizing statewide allocation of multiple vaccine types to multiple priority groups, maximizing equal access. We assume a large fraction of available vaccines are distributed to healthcare providers based on their requests, and then optimize county-level allocation of the remaining doses to achieve equity. We have applied the model to the state of Texas, and incorporated it in a Web-based decision-support tool for the Texas Department of State Health Services (DSHS). Based on vaccine quantities delivered to registered healthcare providers in response to their requests during the 2009 H1N1 pandemic, we find that a relatively small cache of discretionary doses (DSHS reserved 6.8% in 2009) suffices to achieve equity across all counties in Texas.

pair (i, j). The optimization model is as follows: s.t. (i,j):(i,j,k)∈G If all location-priority group pairs are equally important (i.e., all weights w ij are set to one) then the objective function in (1a) is the sum of the squares of vaccine shortfalls, weighted by the location-priority group population sizes across all such pairs. In the optimization model, we minimize this sum of population-weighted shortages. Below, we provide sufficient conditions under which model (1)'s optimal allocations maximize proportional fairness.
Constraint (1b) ensures that the allocated doses of vaccine type k do not exceed the available doses; constraint (1c) indicates that coverage cannot exceed one for any eligible locationpriority group pair; and, constraint (1d) precludes negative allocations.
We add weights (w ij ≥ 1) for each location-priority group pair to acknowledge their relative importance, where the value of the least important location-priority group pair under consideration is one. When the weights differ, instead of seeking equal coverage for priority groups in each location, we seek equal coverage when weighted by the proportionality constants, w ij . So, if one location-priority group pair has twice the weight of another, we seek twice the coverage rate for the higher priority pair. As Theorem 1 shows, this can indeed be achieved for two location-priority group pairs under the following conditions: (i) both of their final coverage rates, after discretionary doses are allocated, are less than one, and (ii) these two pairs have at least one common type of available discretionary vaccines allocated in an optimal solution.
With the proper inputs, model (1) provides an optimal allocation, denoted Q * ijk , and, in turn, the optimal coverage rate of each pair, denoted g ij , i.e., Here, we show that the objective function in (1a), when combined with the constraints (1b)-(1d), seeks to provide an allocation of available discretionary doses in a manner of proportional fairness.
Theorem 1. By using model (1) to allocate available discretionary vaccine doses, the optimal coverage rates of two location-priority group pairs, (i, j) and (i , j ), denoted g ij and g i j , are proportional to their weights, w ij and w i j , if (i) both g ij and g i j as defined in equation (2) are less than 1; and, (ii) (i, j) and (i , j ) have a common type of available discretionary vaccine allocated in an optimal solution.
Let F (Q) denote the objective function in (1a). Then, we have As a result, the Karush-Kuhn-Tucker optimality conditions for the convex quadratic program (1) (see, e.g., [1]) are then: Primal feasibility: constraints (1b), (1c), and (1d) Dual feasibility: Complementary slackness: Consider two (i, j) and (i , j ) pairs that satisfy hypotheses (i) and (ii). By hypothesis (i) and complementary slackness condition (4b), we have ν * ij = ν * i j = 0. By hypothesis (ii), we have Q * ijk > 0 and Q * i j k > 0 for some k and hence by (4c) for that k: Using the definition of g ij from equation (2), we have

Secondary optimization model
We index geographic health service regions by h ∈ H. We form three more index sets: and P , where I h is the subset of counties in region h; O is the subset of priority group-vaccine type pairs (O ⊆ J × K) for which priority group j is suitable for vaccine type k; and, P is the subset of region, priority group, vaccine type triples (P ⊆ H × J × K), which indicates priority group-vaccine type pairs in region h with (j, k) in set O.
We define decision variables V jk , Y hjk , and Y jk : V jk is a binary variable, which takes value 1 if we allocate discretionary doses of type k to priority group j and takes value 0 otherwise; Y hjk is the regional coverage of priority group j in region h from discretionary doses of type k; and, Y jk is the average regional coverage of priority group j from discretionary doses of type k. Instead of requiring the optimal coverage levels from the primary model, we allow the final coverage of a location-priority group pair to be at most ε away from optimality to provide more opportunity to improve secondary objectives. In our analysis presented in the main text, we use ε = 0.001.
The secondary optimization model is as follows: The objective function in (5a) consists of two terms. The first term represents the number of vaccine type-priority group pairs that receive doses, and the second term measures the variation of vaccine types allocated across health service regions. In order to understand the second term, fix a (j, k) pair for the moment. As defined mathematically in constraints (5g) and (5h), Y hjk is the regional coverage of priority group j in region h that comes from discretionary doses of vaccine type k, and Y jk is its average across all health service regions.
The average absolute deviation, represents the dispersion vaccine type k contributes to the coverage of priority group j among regions. By summing this average absolute deviation across j ∈ J, we obtain the dispersion vaccine type k contributes to the coverage among regions. Lastly, we sum across k ∈ K to form a measure of variation of vaccine types allocated among health service regions. We could assign different weights to the two terms in (5a) to adjust the relative importance of policy simplicity and geographic equity. However, we obtain desirable results, as we describe in the main text, with equal weights.
In addition to constraints (5b)-(5d), which replicate constraints from the primary model, we have five other sets of constraints to satisfy. Constraint (5e) ensures that the final coverage of a location-priority group pair is no more than ε away from the optimal proportionally fair coverage of the primary model. Constraints (5f) and (5i) calculate the number of vaccine type-priority group pairs that receive doses using the "big M " method, where M is a constant that is large enough so that constraint (5f) is vacuous if V ik = 1. Constraints (5g) and (5h) define Y hjk and Y jk for measuring the variation of vaccine types allocated among health service regions.
Model (5) takes as input the optimal proportionally fair coverage rates (g ij ) from the primary model. Then, model (5) provides an optimal allocation, denoted Q * * ijk , and the associated coverage rate of each location-priority group pair, which considers proportional fairness and the secondary objectives.

Post-processing for integer-valued allocations
By using the primary and secondary models, we allocate available discretionary doses to eligible location-priority group pairs in a proportionally fair manner with two secondary objectives. However, these two models ignore integrality of vaccine doses. Hence, we construct a post-processing step to find a near-optimal solution that allocates integer-valued discretionary doses. First, for each vaccine type k, we round each Q * * ijk from model (5) down to the nearest integer, denoted Q * * ijk . Next, we calculate the difference between the sum of Q * * ijk across all (i, j) pairs in subset L and that of Q * * ijk . This difference represents extra doses from the fractional part of each Q * * ijk . We sort eligible (i, j) pairs in descending order of their fractions and add one dose to each Q * * ijk by this order until all the extra doses are allocated.

Computation
The primary and secondary optimization models are solved using the commercial software modeling language GAMS [2], which calls the CPLEX [3] optimization solver. We can run the software to solve the models on a server at http://flu.tacc.utexas.edu/ or on a laptop. On a 1.7 GHz Intel Core i7 MacBook Air using GAMS version 24.7.3 and CPLEX version 12.6.3, the three-step solution process takes about 4 seconds for a single value of the discretionary reserve (e.g., 6.8%). This computation time includes the software creating a