Adaptive and optimized COVID-19 vaccination strategies across geographical regions and age groups

We evaluate the efficiency of various heuristic strategies for allocating vaccines against COVID-19 and compare them to strategies found using optimal control theory. Our approach is based on a mathematical model which tracks the spread of disease among different age groups and across different geographical regions, and we introduce a method to combine age-specific contact data to geographical movement data. As a case study, we model the epidemic in the population of mainland Finland utilizing mobility data from a major telecom operator. Our approach allows to determine which geographical regions and age groups should be targeted first in order to minimize the number of deaths. In the scenarios that we test, we find that distributing vaccines demographically and in an age-descending order is not optimal for minimizing deaths and the burden of disease. Instead, more lives could be saved by using strategies which emphasize high-incidence regions and distribute vaccines in parallel to multiple age groups. The level of emphasis that high-incidence regions should be given depends on the overall transmission rate in the population. This observation highlights the importance of updating the vaccination strategy when the effective reproduction number changes due to the general contact patterns changing and new virus variants entering.

where v(t) is the overall national number of available vaccine doses on day t, w 1 , w 2 , and w 3 are tunable weight parameters of the strategy ( i w i = 1), I D k (t) is the number of new infections, and H D k (t) is the total hospital occupation in region k over the last D days. In our work we set D = 14 to capture the changes over two weeks starting from the initial date. The total number of new infections in region k in the last D days is computed by and similarly the number of hospitalized individuals is computed by H w kg (t) + H c kg (t) + H r kg (t) dt.
Note that, since we want H D k (t) to reflect the total hospital occupation at time t for region k, an individual may be counted more than once. One at day t, another at t − 1 and so on. As long as the individuals remain in the hospital they will be counted.
Different vaccination strategies can be obtained by changing the weights w i , e.g. setting w 1 = 1 and w 2 = w 3 = 0 corresponds to the baseline strategy Pop where vaccines are equally distributed according to the population density.

Numerical discretization of optimal control problems
In this section we describe the numerical approach to solving a general optimal control problem via Pontryagin's Maximum Principle [1,2].
T ∈ R n be a state vector and u(t) = [u 1 (t), u 2 (t), . . . , u r (t)] T ∈ R r be a control vector. Consider the following optimal control problem: find u(t) to minimize subject to the state equationsẋ In order to state the Maximum Principle we define a Hamiltonian as where q(t) are the time-dependent Lagrange multipliers [1]. The goal now is to find an optimal trajectory x(t), an optimal control u(t) and an optimal set of Lagrange multipliers q(t) so that to minimize the objective function in (2.1). From the Hamiltonian (2.3) and Pontryagin's Maximum Principle, we obtain the following theorem [1].

(2.4)
For a given initial value x 0 ∈ R n , the numerical approach now consists of finding functions x : [0, t f ] → R n , u : [0, t f ] → R n and q : [0, t f ] → R n satisfying the optimality system (2.4). The numerical algorithm consists of the following steps: Step 1: Subdivide the interval [0, t f ] into N equal sub-intervals and assume a piecewiseconstant control u (0) (t) = u (0) (t k ), t ∈ [t k , t k+1 ], k = 0, 1, . . . , N − 1 Step 2: Integrate the state equations forward in time for the assumed control u (i) and store the trajectory x (i) Step 3: Compute q (i) by solving the second equation in (2.4)

backwards in time
Step 4: Compute a new control u i+1 by solving a finite-dimensional nonlinear optimization problem using a sequential least squares programming algorithm Step 5: Compute x (i+1) and q (i+1) for the new control variable as in Steps 2 and 3 Step 6: Compute the values J (i) (u (i) , x (i) ) and Step 7: If stop the iterative procedure. Here ϵ is a small positive constant used as a tolerance.
and return to Step 4.

Different type of vaccines
In this section we investigate the sensitivity of the optimization algorithm to different levels of the vaccine efficacy e, reduction in susceptibility ω and protection against developing severe illness π.
We set e = 0.5, 0.7, 0.9, ω = 0, 0.2, 0.6 and π = 0, 0.2, 0.6 and test the robustness of the optimized strategies with respect to different combinations of these parameters. Let be a Cartesian product representing the set of all combinations for the different values of the parameters e, ω, π. Let S k be the optimized strategy obtained for P k ∈ P , k = 1, 2, . . . , 27. We investigate the difference between the optimized vaccination strategies, i.e., where D is the total number of deaths. Further, we choose the value of the mobility parameter τ = 0.5.Taking the max norm of the Σ = (σ kl ) ∈ R 27×27 matrix, we get Hence the difference in the total number of deaths is less than 1 individual and the performance of the optimized strategies is the same for combinations for the different values of the parameters. Further, Figs M and N in S2 file verify that the optimization algorithm is robust to different values of the parameters e, ω, π since the optimized strategies are similar. The algorithm is expected to be robust against changes in the model parameters concerning the efficacy of the vaccine or protection of the vaccine against severe illness since the gradient of the optimization algorithm does not explicitly depend on these parameters.

Data and parameters
4.1 Demographic data for Finland

Initial conditions
We obtain the initial conditions from data trying to mimic the pandemic situation in Finland as of 18 April 2021. More specifically, we calculate the initial conditions for the compartments in Table D  0-9   10-19  20-29  30-39  40-49  50-59  60-69  70-79  80+  Total   HYKS  221613  238313  272674  316173  285988  289128  256006  212089  106198  2198182  TYKS  82812  93001  103572  106093  101979  111874  113383  99917  56373  869004  TAYS  88071  100864  105275  112809  106951  115157  117896  100045  55613  902681  KYS  71910  84213  92466  91390  85302  103387  119723  95591  53252  797234  OYS  80308  91471  84511  88448  82348  91225  100322  75669  42261  736563   Total  544714  607862  658498  714913  662568  710771  707330  583311 313697 5503664  as follows where v kg is the cumulative number of people who have received the first dose of any vaccine until 18 April 2021, i r kg stands for the estimated number of real infectious, r r kg represents the number of real recovered people as of 18 April, h k is the reported individuals in hospital ward, c k is the reported individuals in critical care units, and H g and C g are the proportions of people at ward Table C: Finnish regional morning (between 6:00-11:59) mobility, averaged over March-May 2019. Rows represent origins and columns represent destinations [5]. and critical care, respectively. The estimation of real infectious individuals at any day t is derived directly from data as follows: where the number of undetected infectious people i u kg (t) come from upscaling the number of detected individuals i d kg (t) by a factor that depends on the index of age group g, i.e., The number of detected infectious people is calculated by summing the reported cases over the last is the number of cases in region k reported by THL (Finnish Institute of Health and Welfare) at day t, and I w g (t) is the proportion of infected people in age group g. We do not have daily counts as THL does not provide these on infected people per age group. We have chosen 18 April 2021 as the start day since it is a Sunday, and the weekly proportion I w g (t) is the same for all the sums (T I + T E = 7 days, 1 week), which gives The numerical value for I w g can be found in Table E and for the result of the summation ω i d k (ω) see Table F. The estimation of real recovered people at day t is similar, in which ω = 0 marks the beginning of the coronavirus epidemic in Finland. The estimated values at 18 April 2021 of the real infected people i r kg and real recovered people r r kg can be found in Tables  I and J, Table J: Estimated number of real recovered people in Finland by region with 9 age groups and 10y age resolution as of 18 April 2021. The entry on row g and column k indicates the number of individuals who are recovered in age group g and region k.