Optimal Vaccination in a Stochastic Epidemic Model of Two Non-Interacting Populations

Developing robust, quantitative methods to optimize resource allocations in response to epidemics has the potential to save lives and minimize health care costs. In this paper, we develop and apply a computationally efficient algorithm that enables us to calculate the complete probability distribution for the final epidemic size in a stochastic Susceptible-Infected-Recovered (SIR) model. Based on these results, we determine the optimal allocations of a limited quantity of vaccine between two non-interacting populations. We compare the stochastic solution to results obtained for the traditional, deterministic SIR model. For intermediate quantities of vaccine, the deterministic model is a poor estimate of the optimal strategy for the more realistic, stochastic case.


Introduction
As rapid, long-range transportation becomes increasingly accessible, transmission of infectious diseases is a growing global concern. Advances in biomedical therapies and production have enabled the development of large quantities of pre-pandemic vaccine [1]. The United Kingdom, Japan, and the United States have plans to stockpile 3.3 million, 10 million, and 40 million doses, respectively, of pre-pandemic H5N1 vaccine [2]. However, in the face of a spreading pandemic, even seemingly extensive resources would be insufficient to provide global coverage, mandating the development of effective protocols for the allocation of limited vaccine [3] [4].
A starting point for many studies of disease transmission in populations is the Susceptible-Infected-Recovered (SIR) model introduced by Kermack and McKendrick [5]. In this model, at any given time each individual is in one of three states. The dynamic evolution of the population is described by two irreversible transition probabilities: one describes the rate at which a susceptible individual becomes infected, and the other describes the rate at which an infected individual recovers (or dies). The net effect of both transition rates can be described by a single number, r 0 , which characterizes how effectively the infective agent moves through the population. On average, r 0 describes the number of susceptible individuals infected by a single infected individual in a population of susceptible individuals. Reducing the number of susceptible individuals in a population, via vaccination for example, decreases the effective reproductive number r eff . When r eff < 1, the number of infected individuals tends to decline, whereas if r eff > 1, the number tends to grow.
The concept of a reproductive number for an infectious disease can be generalized to more complex models of epidemiology. Previous work on developing optimal vaccination strategies typically focus on minimizing r eff , either by proactive dispersal of vaccine before the infection reaches a population [6], or reactive dispersal [7] [8], after infection has been detected in a group. In both cases, the overall size of the epidemic, measured by the total number of individuals who have been infected throughout the course of the epidemic, is lowered. Vaccination accomplishes this by removing an initial number of susceptible individuals and, thereby, also suppressing the rate of infection.
Numerous computational studies of large-scale veterinary infections, such as foot-andmouth disease [9] [10], Johne's Disease [11], as well as human infections like measles [12] and SARS [13], have been performed. In models that aim to capture field observations, detailed, case-specific information such as demographics of the population, timing and logistics for vaccine deployment, delays associated with the immune response, and overall vaccine efficacy are often essential to the investigation. In all cases, there are tradeoffs between complexity and realism, and between computational viability and the generality of the results.
In this paper, we abstract the geographic, demographic, and disease specific information, and instead focus on the fundamental problem of stochastic SIR dynamics with prophylactic vaccination in two non-interacting populations (e.g. two well-separated cities). Previously, Keeling and Shattock [2] considered the deterministic SIR model in this scenario, and obtained striking results. As the total amount of available vaccine is increased, the allocation of vaccine that minimizes the total number of infected individuals can undergo discontinuous transitions. With a small amount of vaccine, the optimal strategy involves ensuring that the smaller population was well protected foremost. However, with enough vaccine, the optimal strategy switches abruptly to protecting the large population, leaving the smaller population entirely unprotected. These results were well explained in terms of a phenomenon referred to as "herd immunity," whereby immunization of a fraction of a population protects even those who are not vaccinated, by reducing the effective reproductive number r eff to a value below unity. Vaccination removes susceptible individuals from the population. If there are less susceptible individuals in the population, on average, an infected individual will infect fewer individuals. Herd immunity occurs when r eff < 1, i.e., on average, at the start of the epidemic each infected individual transmits his disease to less than one person. Keeling and Shattock explained the sharp transitions in the optimal strategy as arising from a strategy that aims to induce herd immunity in the largest population possible.
While the deterministic SIR model is characterized by two coupled ordinary differential equations, the stochastic SIR model involves a high dimensional state space with probabilistic transitions between partitions of the overall population, characterized by the number of individuals in each state. Stochasticity leads to noteworthy differences in the epidemic size. When r eff > 1 the probability distribution for the total epidemic size is bimodal [14], comprised of a roughly Gaussian peak centered at the deterministic epidemic size, as well as a second peak for small, "terminal infections," describing the likelihood the disease will fail to propagate significantly during the initial phase of infection. The peaks of the distribution are well separated when the population size is large, so that if the number of infected individuals exceeds a critical size, the epidemic progresses to a large size, characterized on average by the deterministic results. However, the non-negligible probability that the disease will fail to propagate in a given population results in significant differences for optimal allocation of vaccine over a wide range of parameters.
The rest of this paper is organized as follows. In Methods we review the deterministic SIR model, and its stochastic generalization. We approach the stochastic problem using a master equation for the time evolution of the complete probability distribution for the number of individuals in each state. Building on the computationally efficient algorithm recently developed by Jenkinson and Goustias [15], we introduce a modification which leads to even greater computational efficiency. In Results we compute the probability distribution for the final epidemic size for a range of parameters to identify regimes for which the stochastic and deterministic models differ most significantly. We compute the optimal allocation of vaccine between two non-interacting populations, and compare our results with the deterministic case. Stochastic effects are most pronounced in situations involving an intermediate amount of resource availability. We conclude with a discussion of our results and future directions.

Methods
We briefly review the deterministic SIR model [5], a system of coupled differential equations for modeling the growth of an epidemic in a well-mixed population within which all agents interact equally with all other agents. The model describes a population of N individuals divided into three classes, susceptible S, infected I, and recovered R: We may omit the equation for the recovered class because we can always deduce the number of recovered individuals from the fact that the total number of individuals in the population, N, is fixed, so that R(t) = N-S(t)-I(t). Equations 1-3 can be thought of as a mean field theory where the continuous variables S(t) and I(t) are the average values (over many iterations) of two discrete integer-valued variables S and I. At any time then, the system can be characterized as being in a state (S, I), which can undergo one of two transitions: ðS; IÞ ! ðS À 1; I þ 1Þat rate bSI; ðS; IÞ ! ðS; I À 1Þ at rate gI: The parameters β and γ can be defined in terms of physical observables, the average number of contacts each person makes per day c, the probability of infection through contact p, and the characteristic duration of the infection T: Optimal Vaccination in a Stochastic Epidemic Model For each set of parameters β and γ, the reproductive number, r 0 is defined to be: where S 0 is the initial number of susceptible individuals, S(t = 0). Here r 0 can be interpreted as the average number of new infections a single infected individual will produce in a completely susceptible population. Thus if r 0 < 1, in the deterministic model dI/dt < 0, the number of infected individuals will decline from the initial seed value, I(t) I 0 = I(t = 0), and no epidemic will occur. In our numerical simulations, the value of r 0 is tuned by varying β. Because β is inversely proportional to N, transmission is frequency dependent, and the rate at which each individual makes contacts with others, c, is independent of the population size N.
In this paper we investigate the effects of prophylactic vaccination. Vaccinating V individuals proactively corresponds to removing V susceptible individuals before the epidemic begins, thus lowering the effective reproductive number. This assumes that the vaccine is completely effective. Let r 0 denote the reproductive number prior to vaccination, and r eff denote the effective reproductive number which is achieved after vaccinating V individuals: If a sufficient number of individuals V are vaccinated, r eff may be reduced to a value below unity, so that d I/dt < 0, and, as a result, the epidemic will not grow. Thus the entire population will be safeguarded without vaccinating the entire population. This phenomenon is known as herd immunity.
In the stochastic SIR model, the infection and recovery reactions are modeled as continuous-time Markovian processes. Let ϕ S,I (t) be the probability at time t of a population with S susceptible individuals and I infected individuals with N = S 0 +I 0 where S 0 and I 0 are the initial values of S and I. The evolution of ϕ S,I (t) in time is then governed by: where the first two terms on the right hand side correspond to transitions into the state (S, I) by a susceptible individual becoming infected or an infected individual recovering, respectively, and the third term corresponds to the probability of leaving the state (S, I) through infection or recovery. While the deterministic model tracks the time evolution of two ensemble averaged variables S(t) and I(t), the stochastic model has up to (S 0 + I 0 + 1)(S 0 + 1) * N 2 possible states (I can take values 0 to S 0 + I 0 and S can take values 0 to S 0 ). All the probabilities ϕ S,I (t) can be assembled into a vector and the entire system of equations can then be written in matrix form. We computationally integrate this system of equations using a modified version of Jenkinson and Goutsias's method [15] of Implicit-Euler integration. The matrix A consists of the coefficients of Equation 8, and describes the transition probabilities: The above equation is discretized by introducing a time step, which controls the accuracy of the method: The ordering of the components of the vector in such a way so that A is lower triangular reduces the number of computations needed to solve Equation 10 from OðK 3 Þ to OðK 2 Þ, where K is the length of the vector and scales with the system size. [15] We made modifications to the way the algorithm counts states, enabling considerably faster computational speed, especially as the population size increases. Where Jenkinson and Goutsias take the approach of counting the so-called "degree of advancement," a scenario in which each state corresponds to a specific sequence of reactions, we instead take the "population process" approach by enumerating all states of the system without tracking which reactions might have led the system to the state in question. In both methods one begins with (S 0 + I 0 + 1)(S 0 + 1) states. In our method, we remove those states that have zero probability of occurring, but are included in the original computational algorithm. For example, many states where S + I > N are retained in the degree of advancement procedure but are explicitly excluded in our method. The result is that we track [(S 0 + 1)(I 0 + 1) + (S 0 + 1)(S 0 )/2] states, which in the limit I 0 ( S 0 is approximately *N 2 /2. As the system size N grows, the difference in the total number of states between the two methods can significantly impact the time it takes to integrate the system of equations. The system is initialized with a population size N, I 0 infected individuals, and an initial reproductive number r 0 . Thus, at time t = 0 the probability of state ϕ N − I 0 ,I 0 (0) = 1 and the probability of all other states equals zero. The collection of probabilities of all accessible states is then evolved forward in time until the distribution reaches a stationary state where the probability of having any state (S, I) where I > 0 is vanishingly small. At that point, all individuals in the initial population of size N, have either been infected, and are now recovered, or remain susceptible. For the parameters considered here, we observe that an integration time of t = 200 is sufficient in all cases. Once the simulation is complete, we define the final epidemic size E as: N À SðtÞ in the stochastic model; Fig. 1 illustrates numerical results for the epidemic size distribution P(E), describing the probability of having a total of E individuals infected over the course of the entire simulation period. We observe that for r 0 > 1 the probability distribution consists of two parts. On the left side of Fig. 1A, there is peak describing small, "terminal infections," which fail to propagate significantly in the population (i.e. the infection terminates before a large number of individuals are impacted). The peak describing terminal infections decays approximately exponentially from the peak value at P(I 0 ). In the stochastic model, there is always a nonzero probability the infection will end without becoming a large scale epidemic.
When r 0 > 1, the distribution P(E) exhibits a second peak towards the right side of Fig. 1A, describing "large-scale epidemics." This peak is approximately centered at the epidemic size predicted by the deterministic SIR model, illustrated for each value of r 0 by the corresponding vertical dashed line in Fig. 1A. The size of the large-scale epidemic scales with the size of the population N, resulting in increasing separation of the peaks for increasing population sizes. To quantify the likelihood of a terminal infection versus a large-scale epidemic by the relative weight associated with each of the peaks, numerically we define the point separating the terminal infection and the large-scale epidemic to correspond to the local minimum in probability that exists between the two peaks. The likelihood of a terminal infection, represented by the total weight in the terminal infection peak, decreases with increasing values of r 0 and I 0 . As r 0 approaches unity from above, the large-scale epidemic progressively decreases in mean size, but increases in variance. Eventually the distinction between terminal infections and largescale epidemics vanishes (the local minimum in P(E) ceases to exist). This is associated with a critical phase transition [16], and occurs at a value of r 0 that approaches unity as the population size N tends to infinity. When r 0 1, the probability distribution is described only by terminal infections.
The cumulative epidemic size distribution P(E < E max ) describes the probability of having an epidemic of size less than E max , and is shown in Fig. 1B. The extended flat portions of the curves indicate that a population of N = 500 individuals is well within the large population limit, defined by a large separation between the terminal infection and large-scale epidemic, with little probability of observing an epidemic size in between the two. Fig. 1B also illustrates how the total probability is distributed between the terminal infection and the large-scale epidemic. The smaller the value of r 0 , the greater the likelihood that the initial seed population of infected individuals will fail to spread the disease.

Results
Our aim is to highlight key differences between stochastic and deterministic approaches to developing a framework for the optimal allocation of vaccine between two non-interacting populations. We begin by observing how a single population reacts to different levels of vaccination in the stochastic and deterministic SIR models; the results provide the basis for the optimization process. Subsequently, we determine the optimal deterministic and stochastic solutions that minimize the average epidemic size between two populations, where one is twice the size of the other, and contrast their properties. We then demonstrate the robustness of these results by computing the corresponding optimal solutions over a range of alternative parameters, that include variations in the ratio of population sizes, and increasing the number of individuals who are initially infected. Finally, we consider an alternative optimization based on imposing a  Figure B illustrates the corresponding cumulative probability distribution P(E < E max ) or the probability of an epidemic of size less than E max . This also shows the relative weight in the terminal infection and in the large-scale epidemic. In each case N = 500 and I 0 = 1. maximum tolerance for the epidemic size and show that the stochastic optimal solution better fulfills this measure than the deterministic optimum.

Impact of Vaccination on the Epidemic Size of a Single Population
We first consider how the epidemic size within a single population decreases as a function of increasing vaccine allocation. Vaccine allocation V removes V susceptible individuals from the initial state (S 0 , I 0 ) ! (S 0 − V, I 0 ) after which the stochastic SIR model evolves according to Equation 8 (Equations 1-3 for the deterministic SIR model). The resulting dynamics determine the size of the epidemic according to Equation 11. Decreasing the initial number of susceptible individuals S 0 by V will not in general lead to a corresponding reduction V in the final epidemic size. An important quantity for optimizing the allocation is the incremental reduction in the expected epidemic size per incremental increase in the allocation.
In Fig. 2 we illustrate the numerical results for a population of N = 500 individuals with different initial numbers of infected individuals I 0 and different reproductive numbers r 0 . An amount of vaccine V (0 V N − I 0 ) is given to the population and we compute the average final epidemic size hEi = R P(E) × E dE as a function of V, where P(E) is computed as in Fig. 1A. We also plot the corresponding deterministic curve in each case, where P(E) here is described by a δ-function at the deterministic epidemic size δ(E). In the stochastic model, the quantity hEi depends on the statistics of both the terminal infection and the large-scale epidemic; hEi may not correspond to an epidemic size that is likely to be observed because there may be a large separation between the observed sizes of terminal infections and large-scale epidemics, with the mean size lying somewhere in between.
Herd immunity occurs in the deterministic SIR model when the initial effective growth rate of the number of infected individuals in the population becomes less than unity (r eff < 1) [17], and is achieved at a value of V determined by Equation 7, i.e. when V/S 0 = V/(N − I 0 ) = 1 − r 0 −1 . In the limit of large populations, the fraction that must be vaccinated to achieve herd immunity approaches 1 − r 0 −1 . Thus for a population N = 500, herd immunity occurs approximately when V = 250 for r 0 = 2, V = 400 for r 0 = 5, and V = 450 for r 0 = 10. Approaching this value, the incremental reduction in expected epidemic size per increase in vaccine allocation increases monotonically. Note that the peak epidemic reduction rate occurs for a slightly smaller V when N is finite, compared to the N ! 1 limit. This is due to the non-negligible (compared to N) contribution of the initial seed population of infected individuals I 0 in the definition of the herd immunity threshold In the stochastic model, the corresponding transition is subtler. Increasing the vaccine allocation has three effects on P(E). It decreases the mean size hEi and increases the variance of the large-scale epidemic, and also increases the relative likelihood of terminal infections. We associate the onset of "effective herd immunity" in the stochastic model with the value of V for which the distinction between terminal infections and large-scale epidemics ceases to exist, as measured by the existence of a local minimum in P(E). Because of the probability of terminal infections, this generally occurs at a value of V which is smaller than that of the herd immunity transition in the deterministic model. Furthermore, unlike the deterministic model, in the stochastic case approaching the onset of effective herd immunity does not coincide with a specific value of r eff and is not generally the point of maximum impact per vaccine in the allocation (as measured by reduction in the average epidemic size).
Before the population has reached the deterministic herd immunity transition, i.e. when r eff > 1, the deterministic epidemic size E (dashed lines of Fig. 2A, B, C) is generally larger than the average stochastic epidemic size hEi (solid lines). While the maximum size of the large-scale epidemic can be greater than the deterministic epidemic E, the average epidemic size hEi is smaller due to the fact that the stochastic model includes the possibility of a terminal infection.
When sufficient vaccine is available to establish herd immunity in the deterministic model, the situation is reversed, and the deterministic size E is generally smaller than the average stochastic epidemic size hEi. When r eff < 1, in the deterministic model dI/dt < 0 at t = 0, and the initial number of infected individuals decreases. On the other hand, stochastically there is always a possibility that the initial number of infected individuals will grow. Hence, beyond the herd immunity threshold, the deterministic epidemic E is smaller than the average stochastic epidemic hEi.
The size of the average stochastic epidemic hEi also approaches the deterministic outcome E as both I 0 and r 0 become large. A larger value of r 0 causes each infected individual to infect more susceptible individuals, while larger values of I 0 makes it less likely for every member of the initial group of infected individuals to recover before spreading their disease. Both of these effects decrease the probability of a terminal infection. Fig. 2D shows the standard deviation of the large-scale epidemic. While this quantity does not factor independently into any of the optimization problems considered in this paper, the variation of the standard deviation of the large-scale epidemic with r 0 , I 0 , and V illustrates several key features of the stochastic model that differ from the deterministic case where the standard deviation is a priori zero.
Firstly, the standard deviation for large epidemics is independent of the initial number of infected individuals (i.e. I 0 = 1, 2, 5), as long as I 0 ( N. If the infection grows into a large-scale epidemic, reaching a size comparable to the system size N, the impact of the original number of initial infected individuals I 0 on the SIR dynamics becomes negligible. The standard deviation does, however, increase with decreasing r 0 . This is illustrated directly for several values of r 0 in Fig. 1A. In each case, a smaller value of r 0 implies each infected individual on average infects fewer susceptible individuals. This increases the variability in the outcome arising from stochastic effects. The standard deviation in the large-scale epidemic is also a function of V, the amount of vaccine allocated to the population. When V is small, the standard deviation increases with increasing V. This is attributed to the fact that having more vaccinated individuals reduces r eff . This effect is balanced by the fact that more vaccinated individuals results in fewer available configurations (S, I) for the system to transition into. Hence as V increases further, the standard deviation eventually peaks and then drops sharply to zero. The value of zero corresponds to the disappearance of the local minimum in P(E) separating terminal infections from largescale epidemics, coinciding with our definition of effective herd immunity. The observation that the standard deviation peaks at a value of V just below the onset of effective herd immunity indicates that the largest uncertainty in the size of the large-scale epidemic is expected for allocations just below effective herd immunity.
An important quantity in determining the optimal allocation of vaccine is the "gain" G, which corresponds to the incremental reduction in the expected epidemic size hEi per incremental increase in the allocation V: In the stochastic model, the dependence of the gain on r 0 can be separated into three distinct cases: the subcritical case r 0 < 1 (not shown), the large r 0 case (r 0 ≳ 2.5 for the other parameters considered here), and the intermediate case where r 0 is greater than the critical value of unity, but below the large r 0 limit (1 < r 0 ≲ 2.5). For r 0 ≳ 2.5, the gain initially increases (at a smaller rate than the corresponding deterministic curve), peaks, and then declines to zero. In this case then, when optimizing vaccine allocation, there is a value of V prior to reaching herd immunity where the gain from vaccination peaks. This is shown for r 0 = 5 in Fig. 3A. For 1 < r 0 ≲ 2.5, the gain G instead declines continuously from a maximum value at V = 1. This is illustrated explicitly for r 0 = 2 in Fig. 3B. This behavior implies that the larger the vaccine allocation V given to the population, the smaller the benefits of even larger V. When r 0 < 1 (not shown), there is essentially zero probability of a large-scale epidemic, thus very little decrease in the epidemic size hEi per increase in allocation V, and the gain is effectively zero. In the deterministic model, for all r 0 > 1 the gain curve follows the same qualitative pattern as the large r 0 case in the stochastic model. Initially, the gain G increases, rising sharply prior to herd immunity, and then falling sharply after the vaccine exceeds the herd immunity point. This means there is a significant increase in the gain G from vaccination as the level of vaccine in the population nears the herd immunity threshold.
What is notably different between the stochastic and deterministic models is the sharpness of the peak and the rate of decline that follows. This is apparent when comparing the stochastic and deterministic curves of Fig. 3A. This illustrates that beyond a threshold level of vaccination (r eff < 1 in the deterministic model) there is almost no reduction in E by increasing the allocation V to the population. In the stochastic model there is not as definitive a threshold level of vaccination. The gain G in the stochastic model begins to decline well before effective herd immunity is reached. Thus in the stochastic model, the point of diminishing returns from vaccination will generally take place at smaller vaccine allocations V compared to the deterministic model.

Optimal Vaccination Allocation for Two Populations
Next we consider the problem of vaccine allocations for two non-interacting populations (e.g., two well-separated cities). This scenario isolates a fundamental tradeoff in resource management, whereby allocating vaccine to one population occurs at the expense of the other.
Unless otherwise specified, we identify properties specific to each population with superscripts 1 and 2. We assume in this section one population is relatively small (N 1 = 500 individuals), and the other is relatively large (N 2 = 1000 individuals). Both populations are initialized with a single infected individual I 1 0 ¼ I 2 0 ¼ 1. In this scenario, both are also exposed to an infection with the same reproductive number r 1 0 ¼ r 2 0 , so that β 1 = 2β 2 , as per Equation 4. A fixed total amount of vaccine V (0 V 1498, where the maximum value of V is given by N 1 þ N 2 À I 1 0 À I 2 0 ¼ 1498; accounting for one seed infected individual in each population) can be distributed between the two populations, so that the small population receives V 1 , and the large population receives V 2 = V − V 1 . We define the optimal allocation to be the partition of V into Optimal Vaccination in a Stochastic Epidemic Model where 0 < E N 1 + N 2 = N = 1500. In this scenario, the cost of producing and distributing vaccine is not taken into account, so it is always beneficial to use all of the available vaccine.
Our objective is to determine the optimal solution as a function of V for the stochastic SIR model. We compare our results to the corresponding optimal solution for the deterministic SIR model, which we also compute. For both models, we sample the space of all possible allocations in order to find the exact optimal solution. This scenario was considered previously for the deterministic case in the limit of large population sizes (N 1 = 100,000 and N 2 = 200,000, using our notation) by Keeling and Shattock [2], who found that for a wide range of values of the reproductive number r 0 , the optimal solution as a function of increasing V was governed by the ability to induce herd immunity in the smaller population (small V), then in the larger population (intermediate V), and finally in both (large V). The colormap of Figs. 4A and 4B illustrates the average epidemic size hEi corresponding to a particular allocation of vaccine, quantified in the figure by the amount of vaccine in the smaller population V 1 . The optimal vaccine allocation is illustrated by the black line, along which hEi is minimized. The range in V 1 is limited by the constraints V 1 V À I 1 0 ¼ V À 1, and V 1 ! V − 999, i.e. neither population receives more vaccine than the number of initial susceptible individuals in that population. This results in the limiting diagonals in the colormap of Figs. 4A and 4B.

Switching
In the deterministic model, Keeling and Shattock [2] found that the optimal solution exhibited "switching" behavior, in which the optimal vaccine allocation makes a significant, discontinuous, change when the total amount of vaccine V exceeds a threshold size. When the amount of vaccine V is below this threshold size, the majority of the vaccine is optimally allocated to the smaller population. When the amount of vaccine V is above this size, all of it is optimally allocated to the larger population. This behavior persists for a wide range of reproductive numbers r 0 > 1 in the deterministic SIR model.
The stochastic optimal solution exhibits switching behavior only for larger values of r 0 (r 0 ≳ 2.9, for the other, fixed parameters considered here). This is demonstrated for r 0 = 5 in Fig. 4A. Switching occurs first at V = 474 and then again at V = 780. The first switching point, above which all vaccine is optimally allocated to the larger population, is present in both the stochastic and deterministic models, although the switching point of the stochastic model occurs at a smaller amount of vaccine V. The second switching point is absent in the deterministic model. The presence of the second switch in the stochastic model is explained in terms of the relative heights of the peaks of the gain curves in the next subsection.
For intermediate values of r 0 (1 < r 0 ≲ 2.9) in the stochastic model, there is no switching behavior, which is in contrast to the results of the deterministic model. It is instead optimal to distribute any given total amount of vaccine V approximately in proportion to the sizes of the populations themselves. This is shown for r 0 = 2 in Fig. 4B. The continuous transition between large r 0 values where switching does take place and small r 0 values where it does not is illustrated in Fig. 4C. As r 0 is decreased, the region between the two instances of switching behavior, where all vaccine is taken out of the smaller population and V 1 = 0, becomes narrower and disappears completely between r 0 = 2 and r 0 = 3 (at r 0 % 2.9).
The conclusion for the stochastic optimal solution is that for intermediate values of r 0 (1 < r 0 ≲ 2.9), the optimal solution is to approximately distribute vaccine in proportion to population size. For large r 0 (r 0 ≳ 2.9), two switches take place. In the deterministic optimal solution, a single switch takes place for all values of r 0 > 1. Note that the stochastic gain curves exhibit peaks only when r 0 ≳ 2.5, while the deterministic gain curves always exhibit peaks. The relationship between switching and the presence of peaks in the gain curve will be discussed in the following subsection.

Understanding the Optimal Stochastic Solution
Next we examine the optimal stochastic solutions of Fig. 4A and 4B in closer detail. Figs. 5A and 5B show the optimal stochastic solution as a solid line and the optimal deterministic solution as a dashed line, with both solutions represented by the fraction of the amount of vaccine V 1 /V given to the smaller population. We seek to characterize the different strategies employed by the optimal solution in the different resource regimes, and also to quantify why the allocation transitions abruptly from one population to another. More broadly, we explain why the optimal stochastic solution, which minimizes the average epidemic size hEi, differs from the optimal deterministic solution, which minimizes the characteristic epidemic size E.
Much of our insight comes from analysis of the gain curves, which are shown in Fig. 3A and 3B for r 0 = 5 and r 0 = 2. By Equation 6, the area under each gain curve G(V) up to a particular value of V is the decrease in the epidemic size due to that amount of vaccine V. We begin with the case when r 0 = 5, which is representative of the large r 0 regime (r 0 ≳ 2.9) where switching does take place. Fig. 5A shows that with small amounts of vaccine V, all of the vaccine is optimally allocated to the smaller population. In the deterministic case, this strategy persists for larger amounts of vaccine V, unless there is enough vaccine for the small population to achieve herd immunity, at V = 400. For a range of V greater than V = 400, herd immunity is preserved in the small population, and the remaining vaccine is optimally allocated to the large population. For the stochastic case, vaccine allocation to the larger population begins for a smaller V, V = 324, above which the optimal solution is to maintain V 1 = 324 and devote the remainder of the vaccine to V 2 .
This difference in strategy can be attributed to the fact that in the stochastic model, the gain G begins to decline well before the onset of effective herd immunity. This is evident in both the N 1 and N 2 solid curves of Fig. 3A. In contrast, the gain in the deterministic model peaks very close to herd immunity at V = 400, as the dashed curves of Fig. 3A show. Compared to the deterministic model, one can attribute this earlier decline in the average epidemic size hEi as being due to the probability of a terminal infection, which significantly lowers the average hEi. More quantitatively, in Fig. 3A, V = 324 is the point at which the stochastic curve for N 1 crosses the initial value of curve N 2 .
As V increases further in Fig. 5A, there is a sharp transition, indicating that if more vaccine exists than V = 474 in the stochastic model or V = 657 in the deterministic model, all vaccine should optimally be allocated to the large population. This is the first switch noted earlier that occurs in both models. As with the earlier transition, the first switch takes place at a smaller amount of vaccine V in the stochastic model than in the deterministic model. This sacrifices herd immunity that could have been achieved in the small population, in favor of relatively larger gains in protection that can be achieved with this level of vaccine in the large population. Quantitatively it is clear from Fig. 3A that beyond a certain amount of vaccine, the stochastic gain curve G(V) begins to decline for N 2 while the stochastic curve for N 1 is still relatively large and constant. Thus around this level of vaccine, all the available vaccine should optimally be switched into the larger population. This same behavior is observed for the deterministic gain curves G(V) for correspondingly larger values of V.
Complete resource allocation to the large population continues until the large population achieves herd immunity, at which point a fraction of the vaccine is allocated to the smaller population. For the deterministic case, the optimal solution retains herd immunity for the large population, and increasingly allocates resources to the small population, until both populations achieve herd immunity. After that point, the optimal solution plateaus. For the deterministic model, the epidemic never progresses (I(t) I 0 ). Because there is no cost for vaccination, remaining resources are allocated based solely on the relative population sizes (i.e. 1/3 for the small population and 2/3 for the large population). For the deterministic model, this corresponds to a situation with excess vaccine, since both populations are fully protected once each has sufficient resources to insure herd immunity.
For the stochastic model, once there is sufficient vaccine to induce effective herd immunity in the large population, at around V = 660, vaccine is once again allocated to the small population. However, unlike the deterministic case, for the stochastic model, there is a second abrupt shift in resources around V = 780, resulting in a cusp in the optimal V 1 /V, with the optimal solution approaching the final population based plateau value V 1 /V = 1/3 from above. This is due to the fact that in the stochastic model, for large r 0 , the smaller population N 1 has a greater peak in gain. Thus if there is enough vaccine available, there is a benefit to removing some vaccine from the large population in order to take advantage of the higher gain in the smaller population. This second switch does not occur in the deterministic model because the opposite is true, the peak of deterministic curve N 2 for the larger population is always higher than the peak of deterministic curve N 1 for the smaller population in Fig. 3A.
For the stochastic model, the gain curves of Fig. 3B can also be used to explain the absence of switching behavior for r 0 = 2, which is generally observed for lower values of r 0 (1 < r 0 ≲ 2.9). A significant difference in this case is that the gain decreases continuously. Due to the absence of peaks in the gain curve, the second switch observed for the large r 0 stochastic model, does not occur for small values of r 0 . The absence of the first switch is more subtle and depends on more than just the presence of a peak which exists when r 0 ≳ 2.5 as discussed previously. For the first switch to occur, the peak of the curve G(V) must be large enough to offset the declines in the gain that first population N 1 exhibits. Hence the first switch takes place for a more restrictive set of r 0 , and only when the peak in the gain is sufficiently large (r 0 ≳ 2.9).
In summary, the switching behavior of the optimal vaccination allocation are due, firstly, to the presence of peaks in the gain curves, and secondly, due to the relative heights of these peaks. This explains why in the stochastic model, switching occurs only for large values of r 0 , while in the deterministic model, it occurs for all values of r 0 > 1. Fundamentally, this difference arises from the bimodal nature of the epidemic size distribution P(E).

Conditioning on the Large-Scale Epidemic
In the previous sections we identified and analyzed the optimal vaccine allocations, obtained by minimizing the average epidemic size hEi, when both populations are initialized with a single infected individual I 1 0 ¼ I 2 0 ¼ 1. For such small initial infection exposures, the average epidemic size hEi has significant contributions from the terminal infection and the large-scale epidemic of both populations.
Alternatively, a policymaker may be interested in a more cautious approach which ignores contributions from terminal infections in the optimization, focusing instead on optimal allocation assuming large-scale epidemics are likely to develop. This scenario is achieved in the stochastic model by initializing the system with a relatively large number of initial infected individuals I 1 0 and I 2 0 . Results are shown below for two populations with N 1 = 500 and N 2 = 1000 individuals and I 1 0 ¼ I 2 0 ¼ I 0 initially infected individuals. Both Figs. 6A and 6B illustrate that as the initial number of infected individuals I 0 increases, the stochastic optimal distribution curves (solid lines) are increasingly similar to the deterministic optimal protocol (dotted line), both qualitatively and quantitatively. In the r 0 = 5 case, the second vaccine switch at V = 780, characteristic of the stochastic optimum, disappears between I 0 = 1 and I 0 = 2. Similarly for r 0 = 2, the optimal distribution changes incrementally from an approximately proportional distribution when I 0 = 1 (blue line), to a distribution with a single switch when I 0 = 5 (green line) that has features similar to those of the optimal distribution curve for the deterministic model.
For the case of two populations, N 1 and N 2 , we define the overall gain G to be the magnitude of the incremental decrease in the average combined epidemic size hEi per incremental increase in the total amount of vaccine V. Figs. 6C and 6D illustrate the overall gain curves associated with each of the optimal solutions illustrated in Figs. 6A and 6B, respectively. Note that the overall gain for the two population optimization problem generalizes the notion of gain for an individual population that was introduced in Equation 12 and Fig. 3. For a single population, gain varies smoothly with increasing V. In contrast, for two populations, the overall gain curve incorporates portions of the gain curves for the individual populations, and may exhibit sharp kinks and discontinuous jumps, reflecting changes, such as switching, which occur at transition points where the allocation changes discontinuously from one population to the other.
Comparing the overall gain curves in Figs. 6C an 6D with the corresponding optimal protocols in Figs. 6A and 6B, we see that switching (i.e. discontinuous jumps in the optimal protocol) coincides with a discontinuous increase in gain. Kinks (slope discontinuities) in the optimal  Figures B and D, and populations of sizes N 1 = 500 individuals and N 2 = 1000 individuals. As the initial seed population of infected individuals is increased, from I 1 0 ¼ I 2 0 ¼ 1 (as in Fig. 4), to I 1 0 ¼ I 2 0 ¼ 5, both populations are increasingly likely to experience a large-scale epidemic. Specifically, in the stochastic model, the probability of having large-scale epidemics in both populations for I 0 = 1, 2, 5 is, respectively, 0.6392 (blue), 0.9210 (red), and 0.9993 (green) for r 0 = 5 and 0.2468, 0.5559, and 0.9334 for r 0 = 2. While both the optimal allocation and gain in the stochastic model becomes more similar to the deterministic limit with increasing I 0 , it does not converge to the deterministic limit for finite r 0 due to the non-negligible width of the distribution of large-scale epidemics that is observed in the stochastic case. Optimal Vaccination in a Stochastic Epidemic Model solution (i.e. points where the protocol shifts from complete allocation to one population, gradually increasing the allocation of the other population) coincide with kinks at local minima of the overall gain. Peaks in the overall gain occur at intermediate points, rather than turning points, in the protocol. As the initial number of infected individuals I 0 increases in each of the populations, the amplitude of the peaks in the overall gain curves become increasingly sharp, becoming more similar, yet not identical to, the corresponding overall gain of the deterministic optimal solution.
For finite values of r 0 such as those considered here, the optimal distribution curves also approach but never fully converge to the deterministic optimal solution. Comparing the I 0 = 5 (green lines) and deterministic curves (dotted lines), appreciable differences remain apparent, even though the probability of a large-scale epidemic exceeds 99.9% for r 0 = 5, and 93.3% for r 0 = 2. Even after conditioning on large-scale epidemics, the non-zero width of the probability distribution for large-scale epidemics that is present in the stochastic models, and absent in the deterministic case, factors nontrivially into the optimization problem. This is the primary factor contributing to the difference between the stochastic and deterministic optimal solutions.
We tested even larger values of I 0 (e.g. I 0 = 10, not shown) to verify the non-convergence between the stochastic and deterministic protocols in the large I 0 limit. However, the system states generated by increasing I 0 , leaving S 0 + I 0 = N fixed, become unrealistic if I 0 is taken to be too large. In fact, the optimal allocation curves for I 0 = 10 lie farther from the deterministic curves (dotted lines) than the I 0 = 5 curves (green lines). This is due to the fact that there is only a very small probability of a system initialized with I 0 = 1 and S 0 = 999 individuals, as in the deterministic model, to later transition into a state with I 0 = 10 and S 0 = 990 individuals without the recovery of any infected individuals at intermediate times. Thus, while increasing the number of initially infected individuals I 0 brings the system closer to the deterministic limit at first, initializing the system with too many initial infected I 0 creates an unrealistic scenario.

The Range of Outcomes
The differences in the optimal vaccination protocols between the stochastic and deterministic models can lead to substantial differences in the observed outcomes. The optimal protocols for the stochastic and deterministic models coincide for small quantities of vaccine, where in both cases it is optimal to allocate all vaccines to the smaller population. The optimal solutions also coincide in the limit of large quantities of vaccine, where it is optimal to allocate vaccines in proportion to the population size. The range of possible outcomes for different vaccination strategies indicates that optimization is most important when intermediate amounts of vaccine are available. One way of understanding the potential impact associated with the optimal stochastic and deterministic protocols is by comparing their projected outcomes when applied to the presumably more realistic stochastic SIR model. In this scenario, there is considerable difference between the best and worst possible outcomes and a significant but smaller difference between the stochastic and deterministic optimal solutions.
The dashed gold lines of Figs. 5A and 5B illustrate the difference in the resultant average epidemic size hEi between the stochastic and deterministic optimal protocols. Both protocols are the same in resource rich and resource poor regimes, and hence yield identical results. Figs. 7A and 7B illustrate hEi for both the stochastic and deterministic optimal protocols as well as the worst case allocation.
We define the "worst case" protocol as that which maximizes hEi within the range of allowed allocations illustrated in Fig. 4A and 4B. Together the stochastic optimal solution and worst case allocation define the possible range of hEi at a given value of V. Fig. 7 shows that the difference between the outcome of the worst case protocol, and either the optimal stochastic and deterministic protocols, is substantially larger than the difference between the stochastic and deterministic cases. This is particularly pronounced for smaller values of r 0 , i.e. r 0 = 2. The worst case protocol would involve continuing to place vaccine in a population even after it is near or has reached herd immunity. This is represented by the plateaus where the average epidemic size hEi is not significantly lowered by further vaccinating members of the population. Deterministically, this is evident from the fact that d I/dt < 0 as soon as the herd immunity threshold has been reached. The deterministic herd immunity threshold, serves as an approximate guide for when to stop vaccinating even in the stochastic case.
The differences between the stochastic and deterministic protocols have a complex r 0 and I 0 dependence. The effect on the difference in hEi between the stochastic and deterministic models that is caused by increasing I 0 is different for small compared to large reproductive numbers r 0 . With a small reproductive number, e.g., r 0 = 2, the difference in the average epidemic size between the stochastic and deterministic optimal protocols is largest at an intermediate value of I 0 , I 1 0 ¼ I 2 0 ¼ 2 for the case illustrated in Fig. 7B. In contrast for large reproductive number, e.g. r 0 = 5, Fig. 7A illustrates that the difference in the average epidemic size between the stochastic and deterministic optimal protocols is maximized for I 1 0 ¼ I 2 0 ¼ 1 and decreases steadily as I 0 is increased.

Comparison with Proportional Distribution
In the previous section we compared average epidemic sizes hEi obtained in the stochastic model when the stochastic, deterministic, and worst case outcome protocols are applied. Another relevant protocol for comparison is one where vaccine is always distributed in a manner that is proportional to the population sizes, i.e.
Politically, a proportional distribution of vaccine is expected to be much easier to implement with the public, compared to a protocol that allocates vaccine to one population at the expense of another, even if the expected epidemic size is reduced in the skewed distribution. Fig. 8 illustrates the increase in the average epidemic size hEi for a proportional distribution of vaccine in the stochastic model, compared to the optimal stochastic solution for different Optimal Vaccination in a Stochastic Epidemic Model values of r 0 and I 0 . Alternatively, the curves can be interpreted as the average reduction in the epidemic size when the distribution changes from a proportional to the stochastic optimal distribution. The most significant differences occur for large values of r 0 and I 0 (i.e. r 0 = 5), which is the regime in which the stochastic model is most similar to the deterministic model. In contrast, when r 0 and I 0 are relatively small (i.e. r 0 = 2), the difference between the stochastic and proportional distributions is much less significant. Additionally, for all values of the parameters, when there is sufficient vaccine available to ensure effective herd immunity in both populations, the optimal protocol approaches the proportional distribution, and the difference in the average epidemic size for the two cases approaches zero.

Population Size Variations
Thus far we have considered the case of two noninteracting populations, where one population is twice as large as the other, i.e., N 1 = 500 and N 2 = 1000. In this section, we show that the qualitative characteristics of the optimal stochastic and deterministic protocols are robust to variations in the ratio of the population sizes N 1 /N 2 . For this investigation, the value of N 2 = 1000 is held constant, and we compute optimal stochastic protocols by minimizing hEi for N 1 = 250, 500, 750, and 1000. Fig. 9 illustrates the optimal stochastic (Figs. 9A and 9B) and deterministic (9C and 9D) protocols for r 0 = 5 (9A and 9C) and r 0 = 2 (9B and 9D) for varying values of N 1 . The same qualitative behavior observed in Figs. 4 and 5 for N 1 = 500 is illustrated here for different ratios of the population sizes. For r 0 = 5 the stochastic and deterministic models exhibit switching in a manner that is qualitatively similar to our previous results (two switches in the stochastic case, and one in the deterministic case). And as before, when r 0 = 2, switching is absent in the stochastic model, but preserved, with one switch, in the deterministic case.
The only exception arises when N 1 = 1000, i.e. N 1 /N 2 = 1. In this case, there is no broken symmetry in the population sizes. In the deterministic model, for both r 0 = 5 and r 0 = 2, the optimal protocol allocates vaccine to one population (either population can be chosen), until herd Optimal Vaccination in a Stochastic Epidemic Model immunity is reached. At higher levels of vaccine, the deterministic protocol allocates the remaining vaccine to the other population until it also achieves herd immunity. Subsequently, vaccine is divided proportionally. For the stochastic model, when r 0 = 5 and N 1 = 1000 (Fig. 9A), a protocol that is similar to the deterministic case is optimal, but involves effective herd immunity, and transitions earlier and more sharply to a proportional distribution. However, when r 0 = 2 in the stochastic model (Fig. 9B), the optimal solution transitions into exactly the proportional distribution beginning at a very small value of V.

Alternate Cost Functions
So far, we have defined the optimal allocation as that which minimizes the average epidemic size hEi, a quantity that contains contributions from both terminal infections and large-scale epidemics, but is not necessarily representative of any specific epidemic size that is likely to be observed because of the gap in the size distribution P(E) (Fig. 1A). Choosing to minimize the deterministic result, which is the same as the average large-scale epidemic size, might potentially be viewed as a conservative approach that safeguards against the case in which both populations will experience large-scale epidemics.
Other criteria for optimization may be considered within this framework. For example, one possibility would be to maximize the likelihood of having no (or very few) infections. We computed the distribution of vaccine that would maximize the probability of having no further infection beyond the initial infectious seed. In this case, the optimal protocol allocates all vaccine to the smaller population until every individual is vaccinated, only allocating vaccine to the larger population when V exceeds N 1 .
A cost function that may be of particular interest to policymakers is one that illustrates how, for a given quantity of vaccine V, to allocate vaccine so as to minimize the probability of having an epidemic greater than a particular size. To address this, in the same scenario of two non-interacting populations (e.g. two well-separated cities) with N 1 = 500 individuals and N 2 = 1000 individuals and both populations initialized with a single initial infected I 1 0 ¼ I 2 0 ¼ 1, here we alternatively consider the probability that the epidemic is below some particular threshold tolerance size E max .
A policymaker may be interested in how much vaccine V would be necessary and how it must be allocated between two populations in order to keep the total epidemic below some size E max . We compute the best achievable probability of having an epidemic below a given size E max given a total amount of vaccine V. The results are shown in Fig. 10A.
The sharp color contrast of the diagonal bands in Fig. 10A are associated with step-like changes in probability, arising from the bimodal nature of the epidemic size distributions P(E). Because there is very little probability for an event in the size range between the large-scale epidemic and the terminal infection peaks, when the threshold E max passes through the large-scale epidemic size (which depends on the vaccine allocation) in the small population, the large population, or the sum of the two, nearly discrete steps in probability are observed.
The allocation that maximizes this probability is shown in Fig. 10B, and is a function of both the amount of vaccine V and also E max . Unlike our previous optimization based on expected size (where the corresponding plot depends only on V), here the solution is extremely complex, switching discontinuously and frequently depending on both V and E max , as indicated by sharp grey scale contrasts reflecting boundaries between high and low allocations to the small population. In the resource poor regime (small V, corresponding to the lower horizontal boundary of the color plot) the solution switches from full allocation to the small population, to full allocation to the large population, back to full allocation to the small population. The lower left white triangle in Fig. 10B corresponds to the situation with few resources, and minimal tolerance for the epidemic size. As in the previous stochastic and deterministic solutions aimed at minimizing the average epidemic size, here the optimal solution allocates all resources to the smaller population. In the E max dependent resource rich regime, corresponding to points above the highest diagonal, the maximum achievable probability in Fig. 10A is near unity, and the optimal allocation simplifies to depend only on V (corresponding to horizontal bands in Fig. 10B). However, in intermediate cases, where tradeoffs are most critical, the structure of the resulting solution is much too subtle to be realistically implemented for real populations given a limited amount of vaccine V.
For comparison, we evaluate the corresponding probabilities based on our previous stochastic and deterministic optimal protocols. While both solutions are suboptimal for this alternative criterion, the stochastic solution comes close to the optimal case. Fig. 10C shows this result for the stochastic optimal solution, which replicates much of the green and blue high probability regions above the intermediate reference line. It does a suboptimal job for relatively smaller epidemics in the regions where the amount of vaccine ranges from V = 400 to V = 1000. Fig. 10D illustrates the corresponding results when the optimal deterministic protocol is applied. In maximizing P(E < E max ), the deterministic protocol underperforms compared to the protocols of both Fig. 10A and 10C. Comparatively, the deterministic protocol minimizes the area of the high probability (blue) regions. It does slightly better than the stochastic optimum in roughly the same regions where the stochastic optimum fails compared to the best possible result, from about V = 400 to V = 750.
This shows that the situation does indeed become more complicated when one looks beyond optimizing the average epidemic size hEi. If the goal is to keep the epidemic below some size, given some amount of vaccine, there are indeed regions where the deterministically Optimizing the probability of having an epidemic less than a given size. Figure A shows in color, the optimal (largest) probability of having an epidemic less than some given size (x-axis), given some amount of vaccine (y-axis). Figure B shows the fraction of vaccine V 1 /V in the smaller population that corresponds to optimal probability shown in Figure  optimal solution will yield slightly better results. Most of the time however, optimizing the average stochastic epidemic size gives a result closer to the best possible one of Fig. 10A. These figures thus indicate that the average epidemic size is a potentially useful metric for gauging the effects of stochasticity and will most of the time yield a solution that is preferable to the deterministic optimum.

Discussion
This paper illustrates the viability and power of developing the exact numerical solution of the master equations, done here for the stochastic SIR model. We modified the computational algorithm developed by Jenkinson and Goutsias [15] to obtain even greater numerical efficiency. This is accomplished by removing excess states of the system which have no probability of occuring, but are naturally included in the original algorithm. Even more significantly, our work and that of Jenkinson and Goutsias [15] provide proof of concept for obtaining accurate, exact solutions for SIR-type models, rather than relying on sampling methods [18] or approximations to the master equation [19] [20]. Furthermore, Black and Ross recently demonstrate an efficient and numerically stable computational methodology that computes the final epidemic size distribution, for a broad range of Markovian SIR-type models, without integration using the jump chain [21].
Our analysis focuses on the fundamental tradeoff involving allocation of vaccine between two non-interacting communities of different size. Our procedure involved three steps. First, for each population we separately calculate the probability distribution of epidemic sizes for a given amount of vaccine. Second, we evaluate the expected epidemic size as a function of the amount of vaccine in each population. Third, we impose a constraint on the total amount of vaccine to distribute between the two populations, and determine the optimal allocation which minimizes the expected combined epidemic size of the two populations.
We obtain several results that serve to elaborate and refine principles first identified by Keeling and Shattock [2], who considered the corresponding tradeoff in the context of the deterministic SIR model. Where the deterministic SIR model predicts a definite epidemic size for any given set of parameters, the stochastic SIR model produces a distribution, the characteristics of which significantly impact protocols for optimal allocation of vaccine. Under conditions that promote spread of the epidemic (i.e., the reproductive number r 0 > 1), the distribution of epidemic sizes obtained from the stochastic SIR model is bimodal [14] in the limit of large population sizes, consisting of a peak describing terminal infections, that fail to propagate significantly in the population, and a peak describing large-scale epidemics, which have a mean size well approximated by the deterministic size. For finite population sizes, the distinction between terminal infections and large-scale epidemics vanishes at a value of r 0 that approaches unity as N ! 1.
Both the possibility of a terminal infection and the width of the distribution of the largescale epidemic sizes contribute significantly to differences in the optimal allocation of vaccine for the stochastic model compared to the deterministic case. The differences are most significant for intermediate ranges of vaccine. In contrast, for both the stochastic and deterministic cases, when vaccine is severely limited or abundant, there is little or no difference in the optimal allocation of vaccine between the two models.
Differences in optimal allocations are amplified for intermediate amounts of vaccine because of the strong switching behavior of the optimal strategy. This switching can arise in both the stochastic and deterministic models, but at different points quantitatively, and is not always observed in the stochastic case. If the deterministic protocol is applied to the more realistic stochastic description of the epidemic evolution in the two populations, the performance is suboptimal, leading to a greater average epidemic size than would occur using the stochastic protocol. The difference is most significant for smaller values of r 0 where there is the most significant probability of a terminal infection. The dependence on I 0 , the number of infected individual, is more complex and depends on r 0 , but in the limit where both r 0 and I 0 are large, the results converge to those of the deterministic SIR model. In the absence of vaccine, these quantities both increase the relative weight in the peak describing terminal infections.
Keeling and Shattock [2] attribute the switching behavior to the property of herd immunity, which occurs when the amount of vaccine is sufficient to prevent the epidemic from spreading significantly in the population. Herd immunity occurs in the deterministic SIR model when the initial effective growth rate of the number of infected individuals in the population becomes less than unity [17]. While the optimal deterministic solution approximately distributes vaccine in a manner that achieves herd immunity in the largest possible population, this is not exactly the case. More precisely, the sharp transitions in both the deterministic and stochastic models arise from optimizing the overall impact of the vaccine in reducing the joint epidemic size, which we attribute to maximizing the overall gain. In the deterministic model, the impact of vaccine on epidemic size reduction is maximized as herd immunity is approached. For the stochastic model, the maximal impact typically occurs earlier, and in some cases, there is no sharp, intermediate transition.
Policies involving strong switching may be difficult to implement publically, as one community could be reluctant to voluntarily sacrifice their entire vaccine allocation to another community in favor of a reduction in the overall epidemic size. In contrast, allocations that are proportional to population size are likely to be less controversial to implement. Our results demonstrate that in certain scenarios in the stochastic model, a proportional distribution can be justified as nearly optimal. For intermediate values of r 0 (1 < r 0 ≲ 2.9) and small initial infected populations, we find that switching behavior is absent in the optimal stochastic protocol. In these situations, the optimal solution is reasonably well approximated by a proportional distribution. In contrast, in situations where r 0 and I 0 are both large, and large-scale epidemics are likely, we find that the optimal stochastic protocol is more similar to the deterministic case, and proportional distribution results in significant increases in the overall epidemic size. However, the reduction in magnitude of the gain peaks for the stochastic model in Fig. 3, compared to the deterministic case, indicate that the overall magnitude of the benefit (as measured by reduction of the epidemic size), is less sensitive to the precise details of the allocation in the stochastic model than it is in the corresponding deterministic case.
Interestingly, our analysis reveals that compared to the deterministic protocol, the stochastic protocol that minimizes the expected epidemic size, also overall better approximates an alternative target based on specifying a maximum tolerance (or threshold) for the overall epidemic size. This result is somewhat surprising. One might have expected the deterministic model to be more accurate in this case, because it predicts a large-scale epidemic whenever r 0 > 1, and as such might have captured a threshold criterion more accurately. The fact that the stochastic protocol continues to outperform the deterministic counterpart provides additional impetus to include the more complete and accurate stochastic dynamics of epidemic evolution in further studies.
This paper isolates the tradeoff in vaccination allocation between two non-interacting populations, prior to the onset of widespread disease, in order to illustrate the significance of the full stochastic solution compared to deterministic case. Our analysis relies on some strong assumptions, particularly the assumption of non-interacting populations. The extreme switching behavior in the deterministic case results from this non-interaction. It is less clear what the optimal policy will be for the case of weakly interacting populations in the stochastic model.
One might speculate that, in the deterministic limit, the presence of even a modest amount of interaction yields dynamics of a single population.
Our conceptual framework and methods can potentially be generalized to include increasingly realistic situations, including interacting populations and real time allocation of vaccine as the epidemic evolves. In these scenarios, we anticipate detailed monitoring of stochastic effects, as well as incorporation of delays associated with transportation and the onset of immunity, will play a critical role in determining the optimal dynamic protocol, and we expect that the critical differences between the stochastic and deterministic SIR models illustrated here will have an increasingly significant impact in identifying protocols that aid in minimizing the overall epidemic size.
The hope is that the systematic study of such tradeoffs will shed light on the development of effective policies. For example, in the case of epidemic outbreak in a localized geographic region, government officials might have to decide whether to allocate scare vaccination doses exclusively to that region or to allocate the vaccine proportionately for population as a whole. In situations where vaccine doses have been prepositioned geographically, the question of "giving away" vaccines from one region to another will be of intense debate. Thus, issues of fairness will complicate decisions even more. Identifying policies that are close-to-optimal and can actually be implemented is an important topic for future research.