Optimal Vaccine Allocation for the Early Mitigation of Pandemic Influenza

With new cases of avian influenza H5N1 (H5N1AV) arising frequently, the threat of a new influenza pandemic remains a challenge for public health. Several vaccines have been developed specifically targeting H5N1AV, but their production is limited and only a few million doses are readily available. Because there is an important time lag between the emergence of new pandemic strain and the development and distribution of a vaccine, shortage of vaccine is very likely at the beginning of a pandemic. We coupled a mathematical model with a genetic algorithm to optimally and dynamically distribute vaccine in a network of cities, connected by the airline transportation network. By minimizing the illness attack rate (i.e., the percentage of people in the population who become infected and ill), we focus on optimizing vaccine allocation in a network of 16 cities in Southeast Asia when only a few million doses are available. In our base case, we assume the vaccine is well-matched and vaccination occurs 5 to 10 days after the beginning of the epidemic. The effectiveness of all the vaccination strategies drops off as the timing is delayed or the vaccine is less well-matched. Under the best assumptions, optimal vaccination strategies substantially reduced the illness attack rate, with a maximal reduction in the attack rate of 85%. Furthermore, our results suggest that cooperative strategies where the resources are optimally distributed among the cities perform much better than the strategies where the vaccine is equally distributed among the network, yielding an illness attack rate 17% lower. We show that it is possible to significantly mitigate a more global epidemic with limited quantities of vaccine, provided that the vaccination campaign is extremely fast and it occurs within the first weeks of transmission.

Vaccination Suppose that we have M doses of vaccine available for age-group i in city k. Ideally, one would like to use all the available vaccine for the susceptible people of age-group i. In practice however, it is very difficult to track susceptible and infectious individuals, especially given that a fraction of the latter are asymptomatic, so a proportion of vaccine is given to people who are already infected or recovered, and it is hence wasted. To take this into account, we use only a fraction f ki of the available vaccine to vaccinate the susceptible individuals of age-group i on day t. This fraction f ki is set to be equal to the fraction of people in age-group i who are still susceptible on day t. Formally, we have for each k ∈ {1, . . . , K} and i ∈ {1, 2}, So for example, if 40% of the children in city i are susceptible on day t, then we will use only 40% of the available vaccine to vaccinate the susceptible individuals of age-group i on day t, and consider that the remaining vaccine is not used.

Calculation of the contact rates
The paper by Wallinga et al. [3] served as a basis for our computations for the contact rates between and within the two age-groups. In this section the notation follows that used in [3]. Wallinga et al. divided the population into six age classes, a 0-5 class, 6-12 class, 13-19 class, 20-39, 40-59 and 60 or older. They gathered information on the number of conversational partners during a week and the age distribution of these partners. They arranged this information in a matrix M , which they called "social contact matrix". Each entry of the matrix m ij represents the "mean numbers of conversational partners per week in age class i as reported by a participant in class j". We used this matrix to reduce the number of age-groups from six to two by computing weighted averages among the groups.
Let w i be the projected Dutch population in age-group i, and let w tot be the projected total Dutch population on January 1, 1987.
To obtain the contact rates for our two age-groups, we collapsed three age classes into a single one. Namely, we regrouped all the children into a single class (0-19 years old) and all the adults into a single class (≥ 20). To do this, we started by collapsing two age-groups at a time and apply the algorithm recursively until we have only two age-groups as follows. Assume that we want to collapse age-group i and age-group j into a single age-group. Let this new age-group be the group ij. Furthermore, let m ij,k be the mean number of conversational partners per week in age class ij as reported by a participant in class k. Then, m ij,k needs to be a weighted average of m ik and m jk , where the weights are given by the fraction of the population in each age-group, Analogously, m k, ij is defined as Finally, m ij,ij is defined also as a weighted average After repeating this process 4 times, we obtained a 2-by-2 matrix M of mean number of conversational partners with only two age-groups, given by Using the same definition, we obtained a matrix C , with c ij = m ij w tot /w i where the weights are now given by the sums of the populations of each collapsed age-group. So, The matrix C is given by Finally, we normalize the matrix C by dividing it by the c 11 entry, The matrix Cnor is given by The matrix Cnor is the contact matrix used in the system of ODE's.

Implementation
We adapted Pyevolve [4] so that the evaluation function is set to use our objective function f . We used a crossover rate of 0.8 and a mutation rate of 0.02. We ran the genetic algorithm independently three times. In each of these, we initialized it with 50 chromosomes and let it evolve for 80 generations. For each of these three individual solutions, we ran 500 epidemics and computed the mean of the attack rate. The optimal solution was selected among the three individual solutions as the one which had the lowest mean attack rate. We then bootstrapped the 500 epidemic results of the optimal solution 1000 times to obtain 95% confidence intervals for the attack rates and the epidemic prevention potential.
We approximate the expected values in equation (??) by computing the mean over 20 runs of our model, for minimizing the function f .

Genetic algorithm
A genetic algorithm [5,6] is a heuristic algorithm that searches for nearly-optimal solutions to a minimization problem by mimicking nature's evolution. The algorithm is initialized with a set of chromosomes and evolves in generations. A chromosome is a particular solution to the given problem. Each chromosome is evaluated according to some predetermined fitness function and is assigned a score. In our case, a chromosome is a solution to our optimization problem, and we can think of a gene as a particular control vector. We initialize the genetic algorithm by randomly generating 48 feasible solutions in two ways. First, we generate half of the initial population by randomly (uniformly) selecting numbers between 0 and 1 for each city and age-group in the network, and map the resulting chromosome to the corresponding feasible solution using the transformation T defined in the main text (eq. 9). Second, we create the other half of the chromosomes by randomly selecting a city and an age-group and allocating enough vaccine to cover as many people as possible in this age-group. Then, if vaccine is still available, we select another city and age-group and use the rest of the vaccine to allocate as much as possible in this age-group and so forth until we run out of vaccine.

Reproduction number
We performed sensitivity analysis for all the possible scenarios (six different coverages and six different vaccination dates) to evaluate the robustness of our results when the epidemic was assumed to be less transmissible (basic reproduction number R 0 = 1.2), when the epidemic was assumed to be more transmissible (R 0 = 1.8). Figures S6 and S7 present the results when we assumed the virus was less transmissible. All three strategies perform better under this scenario. This is expected, as we know from previous studies [1] that for a low R 0 as this one, the threshold in deterministic models for the number of children and adults needed to be vaccinated to bring R 0 below 1 is much lower. Also, it is important to note that when R 0 = 1.2, the epidemic is much slower, so that after 250 days (the final time considered for this study) the epidemic has not yet peaked in most of the cities. Figures S8 and S9 present the results for a more transmissible virus, with R 0 = 1.8. Again, as expected, here the three strategies presented performed worse than when R 0 = 1.5. The optimal solution always outperforms the other strategies considered, and it is able to reduce the size of the epidemic by more than 40% compared with no vaccination with as few as 4 million doses. Similar to the base case where R 0 = 1.5, the biggest difference between the optimal strategy and the children-only pro rata strategy is found when vaccination occurs early, and this difference tends to decrease as vaccination occurs later. However, this difference is more prominent for R 0 = 1.8 than for R 0 = 1.5 when vaccination occurs during the first few days of the epidemic and we have more resources available (Fig. S8D-F), suggesting that when the virus is more transmissible the optimal allocation of resources is more important, at least at the beginning of the epidemic.
For the rest of the sensitivity analysis, we picked a single vaccination coverage to perform the analysis. Because the most interesting results (i.e, those where the optimal vaccine distribution showed greatest advantage with respect to the other two strategies considered) were obtained when there were five million doses of vaccine, we concentrated the analysis on this coverage but considered all the six different vaccination dates.

Probability of travel
We next analyzed the sensitivity of our results to the probability of infectious people traveling through the network. To the best of our knowledge, there is no data concerning the probability of travel if a person becomes infected and ill. For this reason, we considered that symptomatic infectious people have a 25% reduction in their probability of traveling when compared to an asymptomatic infectious traveler for the baseline cases. Figures S10 and S11 show the results when symptomatic infectious people have a 10% reduction in their probability of traveling (panel A) or a 75% reduction in their probability of traveling (panel B). A little lower reduction in the probability of traveling yields a slight increase in the attack rates ( fig. S10, panel A) and a higher reduction in the probability of traveling yields a decrease in the attack rates. Interestingly, when the probability of travel is reduced by 75%, the optimal strategy is able to prevent 44% of the epidemics ( fig. S11, panel B). This is expected, if symptomatic infectious people are traveling less, the epidemic is spreading more slowly and the optimal strategy is able to prevent the epidemics for a longer time. However, the main conclusions are insensitive to this parameter: the optimal strategy outperforms the other two strategies with the major difference seen when vaccination occurs during the first days of the epidemic.

Probability of travel for infected children
In this section we consider the probability of an infected child traveling through the network. Due to the lack of data for this region, we conservatively assumed that infected children will travel with the same probability than infected adults for the baseline cases. Figure S12 shows the attack rates (panel A) and the EPP (panel B) when the optimization was done assuming that infected children will have a 50% less probability of travel than infected adults. The results are very similar to those presented in the main text. The conclusions were not sensitive to this parameter, but, as expected, we obtained a better EPP when the probability of travel for children was reduced by 50%. As before, if fewer infected people (here, children) are traveling the network, the epidemic occurs at a slower pace and the optimal strategy is able to prevent more epidemics.

Vaccination
Figures S13 and S14 present the results when the vaccine efficacies considered are lower than those presented in the main text. We considered two new vaccine efficacy values: the first one assumed that the vaccine efficacies would be one-third of the original values (VE S = 0.13, VE I = 0.15, and VE P = 0.25, results given in panel A) and the second one assumed that vaccine efficacies would be two-thirds of their original values (VE S = 0.27, VE I = 0.30, and VE P = 0.50, results given on panel B). With very low vaccine efficacies (one third of the original values, panel A), the three strategies perform poorly, but the optimal strategy still performs better than the other two strategies for all the dates considered. With low vaccine efficacies (two thirds of the original values, panel B), the results presented in the main text are preserved: the optimal strategy outperforms the other two strategies with a maximum difference attained if vaccination occurs early (11% difference in the attack rate), but this difference quickly declines as vaccination starts later on. We also performed sensitivity analysis with respect to the vaccination timing. Figure S15 presents the results when we assumed that vaccination would be completed in 10 days. Here, the difference in the attack rates between the optimal strategy and the children-only pro rata strategy is only 5% (panel A). The optimal strategy can contain 36% of the epidemics but only if vaccination starts on day 5.

Optimization when vaccines are only given to children
In this section, we performed the optimization when vaccines can be given only to children. The results are shown in figure S16. The optimal vaccination strategy outperforms the other two strategies. Here, the attack rates are slightly lower than when the optimization is performed over the entire population, but the EPP is very similar. These results suggest that allocating vaccine to the high-transmission groups (in this case children) could be more advantageous than to allocate the vaccine to cities where the epidemic is already under way, or to cities that are well connected.