17 min

Where to Vaccinate First

Published on November 13, 2020

I show that as a rough heuristic, a limited amount of intervention should first be used in regions where the reproduction number is approximately \(R = 1 + \Delta\), where \(\Delta\) is the uncertainty in the estimate of reproduction number. We also reason as to why the optimal allocation involves triage. It brings some countries down to a common, lower reproduction number while at the same time, other regions receive no help at all.    

A Model for Regions

In this note, I made some attempt to improve my own intuition for what will no doubt be a tremendously controversial topic in months to come. I am not an epidemiologist, nor claim any expertise in ethics, but here I present some mathematical commentary in order to provoke some debate, and model improvements.  

I adopt a very simple model in which a number of different populations are considered non-interacting. We assume a one-off intervention where the infection rate is changed once and for all, but thereafter the epidemic runs it course independently in each population.

This stylized model, certainly open to critique, might give us some handle on the "where" but not the "who", for it is silent on the topic of which cohorts (elderly, schoolchildren, etc.) should be vaccinated first. It also ignores contagion from one region to another. 

Because the model does not distinguish between different types of people, we will assume our aim is to minimize the total number of people who are eventually infected. Other more complex models may apply different utilitarian measures, such as weighting for years of life lost. 

Vaccination and Infection Rate 

While this post is prompted by the possibility of vaccination, the outline applies to any intervention that might change the course of an epidemic. For example, a marketing spend could be allocated across different populations. The campaigns might encourage mask use, hand washing and social distancing. 

Vaccination is similar to a lowering of the infection rate - the proportional rate at which the infected population infects the as-yet-susceptible. It lowers the probability that someone you bump into has the virus. This is true whether the vaccine is 10% effective or 90%.

In the toy simulation below, blue and green people infect others. However the yellow dots might as well be non-interacting neutrinos. These are the inert recovered people - never changing to blue (infected) or causing anyone else to change to blue. To administer a vaccine in this simulation is to take a marker and re-color some purple people yellow. 


An enterprising reader might want to use the code (see pandemic library) to consider strategies for inoculation taking local effects into account. Perhaps we should vaccinate those dots in dense areas, or dots that move faster and commute more. However in this note I will content myself with a simple compartmental model that treats sub-populations as homogeneous (and fully mixing). 

With that simplification we can exploit analytic solutions. Inoculation can be viewed as moving people from the cohort of susceptible people to the cohort of recovered (assumed immune) people. But if the intervention is small, we can model a partially inoculated population as one with nobody removed, just with different parameters. Think of giving everyone a one percent effective vaccine, rather than one percent of people a perfect one. It isn't identical, but it's close.

Now let us imagine that in place of one population, there are actually hundreds of them - representing far flung regions on Earth that do not interact. We are able to airdrop in vaccines and inoculate some or all of the populations - but we have a limited supply. So we have to make the terrible decision about which populations to vaccinate the most and which to vaccinate the least. 

In this stylized setup, it may seem as if there are all sorts of topological considerations being left out which are analogous to fire-fighting strategy - it may be best to administer vaccine at the airport. That's mostly a flaw although you can imagine that some of these cells represent regions near the virus "front" and some far from it - so we are not losing absolutely everything with this approximate mental model. 

Total Infected

There is a formula that approximates the total number of people who are going to eventually become infected. The critical ratio is: \begin{equation} r= \frac{Infection\ rate}{Recovery\ rate} \end{equation} and if \(s=S_0\) denotes the initial fraction of the population that is susceptible, the fraction of the population that will catch the virus eventually is given by: \begin{equation} R_{\infty}(s;r) = 1 + \frac{ W \left( - s r e^{-r} \right)} {r} \end{equation} See, for example, Frank Wang's article Application of the Lambert W Function to the SIR Epidemic Model (pdf). I use little \(r,s\) in place of the customary \(S_0, R_0\) to declutter the formulas. 

The function \(W\) is the amazing Lambert W Function, useful for computing everything from the amount of fuel consumed by a plane, to the solution of Wien's displacement law or, if you want \begin{equation} x^{x^{x^{x^{\dots}}}} \end{equation} as high as you wish. I've found it useful over the years to spot equations that can be solved by it, and there are some examples here for those who are interested.

The function is named in honor of the equally amazing John Heinrich Lambert, a self-taught polymath who speculated on Black Holes and non-Euclidean geometry well before they were in vogue, derived beautiful cosmological results, and was able to prove that \(\pi\) is irrational. Lambert's name also lives on in cartography (Lambert projection).

We can view the Lambert W Function as the inverse of the function 

\begin{equation} w \mapsto w e^{w} \end{equation} although that is only true for \(w \ge -1\). As noted by Wang, this means that \(W\left(-re^{-r}\right) = -r\) when the virus is under control at the outset. We assume the initial infected population is a tiny fraction of the total population, which is not necessary, just neater looking. 

Impact of a Change in Infection Rate

If \(r \ge 1\), no simplification to the Lambert W Function applies, unfortunately, but one is free to wield the function as freely as \(\sin\) or \(cosh\) or any other legitimate special function (e.g. scipy.special.lambertw for the Pythonistas). Here's a plot of \(1+W(-r e^{-r})/r\) making the fine line between success and failure apparent. 


In a thought experiment, we are able to move backwards and forward along the x-axis by modifying the infection rate (or the recovery rate, though that seems harder) though some kind of intervention. The cost of administering, we shall assume, is linear in the population size and constant in \(r\) value. Thus reducing \(r\) from \(2.1\) to \(2\) in a population of \(1,000,000\) costs as much as reducing \(r\) from \(2.1\) to \(1.1\) in a population of \(100,000\). This linearity is not entirely realistic, but we must start somewhere. 

It follows that if we have a whole collection of regions, each at different points along the curve, then one reasonable approach would be to focus effort on where the gradient of this line is the largest - right after the bifurcation point at \(r=0\). This suggests that if we have just one single vial of vaccine, and it can be administered to one lucky person somewhere on the globe, we should not give it to someone in the most ravaged location, or the safest, but in fact to someone living in a region that is just beyond the \(r=1\) cusp. For in these locations, there is a chance of radically changing the path of the virus. 


Above I show the derivative of total cases as a function of \(r\). You'll notice that with the assumption \(s=S_0 \rightarrow 1\) there is no benefit whatsoever for intervening if the virus will die out on its own accord. That's counter-intuitive but remember that the initial population, and thus the number of all future infected, is considered negligible in any scenario in which the virus is in decline. 

A somewhat related critique: any given region might itself better be thought of as mixtures of reproduction numbers - or with convexity adjustments as I suggest here . Suffice to say, there are numerous ways to improve this analysis.

Intervention Under Uncertainty

The bigger problem is that life is messy and we don't really know the reproduction number. Suppose we intervene because we believe the reproduction number is \(r=1.1\) and it turns out that it was actually \(r=0.9\). A utilitarian logic suggests our efforts could have been better placed elsewhere. 

Instead, we might rank locations for intervention not by the sharply peaked plot above but rather, by some smoothing of the same that takes into account the error measuring \(r\), as shown in blue below. To compute that line, we assume that the reproduction number is known only within a standard error of \(0.5\). 

The orange line is proportional to the derivative of the blue line, and represents the marginal benefit of intervention. It suggests that regions where \(r \approx 4/3\) should take priority over those with \(r \approx 1\). Now if you worry about uncertainty of uncertainty, perhaps aim for \(r \approx 1+\Delta\) where \(\Delta\) is some philosophically robust estimate of your error. 


This is the main message. This simplified utilitarian analysis suggests that we should start with a person in a region where you are pretty sure the reproduction error exceeds the critical runaway value \(r=1\) - but not by too much. 


So much for determining a fortunate first patient. However, as we continue to administer vaccine, the infectiousness and \(r\) will continue to drop in that region. As can be seen from the orange curve, the marginal benefits are decreasing. It is reasonable to ask what the program should look like if there are \(n\) regions and a total ability to move \(r\) values is limited by some fixed budget. For example suppose we introduce a constraint \begin{equation} \sum_{i=1}^n \delta r_i = c \end{equation} where \(c\) represents our total ability to change reproduction numbers, \( \delta r_i \) is the improvement in reproduction number defined as \begin{equation} \delta r_i = r_i - r_i^* \end{equation} due to intervention. Here \(r_i\) denotes the reproduction number for region \(i\) before intervention and \(r_i^*\) the reproduction number afterwards. For simplicity we assume all populations are the same size and our task is to choose the \( \delta r_i \gt 0 \) to minimize total cases. For this we might try a greedy algorithm, which is to say that we could administer a small amount where it is needed the most, and repeat until we run out of vaccine.

But interestingly, this won't lead to an optimal solution. A region whose reproduction rate has been lowered from \(2.5\) to \(1.5\) now has a lot to gain on the margin - much more than at the outset. A greedy algorithm might not see that, and thus there is a notion of investment in a region. We might want to administer even if the marginal benefit right now isn't myopically better than a competing use.

It turns out that naive hill climbing won't quite cut it for this reason. However we can reason our way to what the optimal solution will look like. The key is to recognize that it will involve a triage - something that won't be surprising to the medical profession.

To see where the triage occurs, one can visualize interventions that move points near the hump in the orange plot. We know we want to move points near the hump because that is where the effect is high. However these interventions move them beyond the hump towards the benign regime \(r<1\) and at that point there is no investment effect. Future gains will be smaller. 

Also, by symmetry, all regions receiving assistance end up in the same place. 

To be mildly formal, and without loss of generality, we can assume the initially measured \(r\) values were in increasing order. Also loss of generality, we can assume that all regions remain in order (i.e. they can merge, but not overtake). Then my statement is that \begin{equation} r^*_1 = r^*_2 = \dots = r^*_k \end{equation} where \(r^*_k\) denotes the reproducibility number after optimal intervention. On the other hand, \begin{eqnarray*} r^*_{k+1} & = & r_{k+1} \\ \vdots & = & \vdots \\ r^*_{n} & = & r_{n+1} \end{eqnarray*} where \(r_k\) denotes the pre-treatment reproducibility number. Regions \(k+1\) through \(n\) have been triaged.

To continue the argument, assume we have found the optimal allocation. At this point advancing a single, less well-off region all the way to the top of the derivative hill must be less worthwhile than moving all \(k\) treated regions a smaller distance towards \(r=1\), for otherwise we could improve the solution. This dictates where we draw the triage line - and it will be \(k\) multiplied by the marginal benefit of further intervention in already treated regions. All regions with reproducibility numbers in excess of this threshold must, logic demands, be left to fend for themselves. 


An example with six regions is shown in the diagram above. Initially, they have equally spaced \(r\) values ranging from \(r=2\) to \(r=5\). The first three are given the vaccine, which brings them to an identical place with the virus well and truly under control. The last three regions go without any vaccine at all. 

I decided to exercise some care in debugging my thinking on this one, by the way, by using Simplicial Homological Global Search implemented by Stefan Endres. That optimizer did well on other problems, as noted in my survey, but more importantly, it is guaranteed to find all global minima (minor technicality, total cases is clearly a Lipschitz function of intervention). See the notebook for an example you can modify, if you don't buy my argument. 

Orange States

Mindful of the fact that intervention can comprise of things besides inoculation, let us roll the tape back to March 2020. Nature magazine provides a plot of estimated \(r\) values, at the time. 


If we assume, bravely, that the error in estimation is on the order of \(0.25\) or \(0.5\) this suggests that the best bang for the buck was to be found not in New York, Florida, New Mexico or other places on a steep trajectory - nor in Minnesota, Wisconsin, Arkansas or Tennessee where the virus was seemingly in retreat.

Instead, middling performers such as California, Oregon and much of the South were, in theory, prime candidates for intervention. It seems unlikely (from the perspective of health leadership) and perhaps unconscionable that the "best" strategy at that time would have involved triaging the red states (those with \(r \ge 2\) in the plot - I'm not referring to politics). The optimal strategy might have diverted all public education dollars to try to save half the country, where \(r \approx 1.5\), through behavioral changes. The mathematics above suggests exactly that, in broad brush terms, though of course this fails to take into account travel. 

Blue Countries

Moving to the global picture this map, published in July 2020 (article) suggests that the greatest marginal benefit at that time existed in (dark) blue countries. However, given the methodological difficulties, it is hard to read too much into any one estimate.  


One can also debate the morality of this allocation, naturally, for all manner of reasons including moral hazard. 

I'm sure it goes without saying that the solution implied by this utilitarian model, or any, is largely academic. We have no reason to believe actual distribution of vaccine will in any way resemble it, and nor would we want it to without further critique of the assumptions. For example one obvious flaw with the analysis is that the impact of a high case load in some regions might be much more severe than elsewhere - due to differing standards of medical treatment. The hospital overloading effect might also dictate a different choice of region dependent objective function. 

However, I hope that this can provide a straw man and that the intuition from this simple model - which I don't think is entirely obvious up front - is more valuable than the flaws in its construction. There is nothing new about the shape taken by the solution. Plenty of other situations, where there is non-linear response to intervention, will have similar characteristics. 

If you'd like to play more with this morbid calculus, or debug it, see the notebook