key: cord-1007333-80q9l7sk
authors: Acuña-Zegarra, Manuel Adrian; Díaz-Infante, Saúl; Baca-Carrasco, David; Liceaga, Daniel Olmos
title: COVID-19 optimal vaccination policies: A modeling study on efficacy, natural and vaccine-induced immunity responses
date: 2021-05-04
journal: Math Biosci
DOI: 10.1016/j.mbs.2021.108614
sha: 71331f43f3db2a53b36099cff1213dafa4cfe958
doc_id: 1007333
cord_uid: 80q9l7sk

About a year into the pandemic, COVID-19 accumulates more than two million deaths worldwide. Despite non-pharmaceutical interventions as social distance, mask-wearing, and restrictive lockdown, the daily confirmed cases remain growing. Vaccine developments from Pfizer, Moderna, and Gamaleya Institute reach more than 90% efficacy and sustain the vaccination campaigns in multiple countries. However, natural and vaccine-induced immunity responses remain poorly understood. There are great expectations, but the new SARS-CoV-2 variants demand to inquire if the vaccines will be highly protective or induce permanent immunity. Further, in the first quarter of 2021, vaccine supply is scarce. Consequently, some countries that are applying the Pfizer vaccine will delay its second required dose. Likewise, logistic supply, economic and political implications impose a set of grand challenges to develop vaccination policies. Therefore, health decision-makers require tools to evaluate hypothetical scenarios and evaluate admissible responses. Following some of the WHO-SAGE recommendations, we formulate an optimal control problem with mixed constraints to describe vaccination schedules. Our solution identifies vaccination policies that minimize the burden of COVID-19 quantified by the number of disability-adjusted years of life lost. These optimal policies ensure the vaccination coverage of a prescribed population fraction in a given time horizon and preserve hospitalization occupancy below a risk level. We explore via simulation plausible scenarios regarding efficacy, coverage, vaccine-induced, and natural immunity. Our simulations suggest that response regarding vaccine-induced immunity and reinfection periods would play a dominant role in mitigating COVID-19.

In late December 2019, the public health officials of Wuhan City, China, reported the emergence of a 30 pneumonia illness of an unknown etiology . The virus (caused by SARS-CoV2) rapidly spread 31 through many countries around the world, causing severe problems on their health systems. The first-32 line defense against the virus included several non-pharmaceutical interventions (NPIs) such as quarantine, 33 isolation, and social distancing, being the main ones. Additionally, some drugs such as baricitinib, combined 34 with remdesivir have been used in the U.S. to treat suspected or laboratory-confirmed COVID-19 hospital 35 Since we will explore disease dynamics that lasts from six months to one year, the model includes vital 133 dynamics. We consider a constant population (N (t) = N ). Thus, we assume that birth and natural death 134 rates are the same and represented by µ. All births lie into the S class and all but class D experience natural 135 death. Class D does not intervene in the transmission dynamics and counts reported deaths. 136 The infection dynamics are as follows: Susceptible individuals (S) become infected, but not infectious, 137 when in contact with infectious individuals I S and I A . Exposed individuals (E) remain in their class until 138 they become infectious and move to either I S or I A . Individuals in class I S either die by disease complications 139 or recover, whereas individuals in class I A move to the R class after a some time. Finally, as the vaccine is 140 considered imperfect, individuals in V move to the E class by interacting with infectious individuals I S and 141 I A at a lower rate than S individuals. Figure 1 shows the compartmental model diagram which summarizes 142 hypotheses mentioned above. The model is given by the following ordinary differential equations system 144

where the infection force is defined by

Here, N (t) = S(t) + E(t) + I S (t) + I A (t) + R(t) + V (t). 

count the cumulative administered vaccines doses until time t as the product N ×X(t). Also, the cumulative 149 incidence of reported cases is given by N × Y IS (t). 150 Remark 1. Vaccine is administered to individuals in classes S, E, I A and R. The amount of given vaccines 151 is quantified by Equation (3). However, as we assume the vaccine has preventive nature only, there is no 152 change from classes E, I A , and R to V due to vaccination. 153 The parameters of model Equation (1) are described in Table 1 .

Parameter Description µ Natural death rate β S (β A ) Symptomatic (Asymptomatic) transmission contact rate ψ V Vaccination rate ω V Waning rate of vaccine. 1/ω V is the average time to lose vaccineinduced immunity ε V Vaccine efficacy σ E Latency rate. 1/σ E is the average latency period p Exposed individuals' fraction who become symptomatic infectious α S Transition rate from symptomatic to recover or death. 1/α S is the average output time of symptomatic individuals class θ Proportion of symptomatic individuals who die due to the disease α A Recovery rate of asymptomatic individuals. 1/α A is the average time which asymptomatic individuals leave being infectious σ R Rate of loss of natural immunity. 1/σ R is the natural immunity period transmission contact rates and exposed individuals proportion who become symptomatic infectious.

To address this problem, we employ a MCMC method. As observation model, we use a negative binomial 162 distribution with mean given by in terms of R 0 , this representation will allow us to formulate disease control strategies.

Considering dynamics (1) without vaccination, the basic reproductive number results:

Note that each term of R 0 represents the contribution of the symptomatic and asymptomatic infectious, 176 respectively, to the spread of the disease.

On the other hand, the vaccine reproduction number for Equation (1) is given by (see Appendix B)

where,

This factor encloses parameters corresponding to the vaccine profile (efficacy, waning), and vaccination rate.

According to this factor and following the ideas of Alexander et al.

[32], the expression (6) can be rewritten

then R V value is lower than one. Thus, for a waning rate given, there exist an adequate efficacy and 186 vaccination rate to control the disease. However, if the conditions in (8) is not satisfied, it will not be 187 possible to reduce the value of R V below 1.

To illustrate the aforementioned, Figure 3 shows the regions where it is possible to reduce the value of 189 R V . In this case, we set all the system parameters as given in Table 3 Vaccination policies to reach a given coverage of a certain percentage of the population in a given period 193 is of great importance. In this sense, we refer to this vaccination constant rate as the base vaccination rate, 194 denoted by ψ V base .

Let W (t) be the normalized unvaccinated population at time t. If no individual has been vaccinated at t = 0, then W (0) = 1. Assuming that we vaccinate individuals at a constant rate ψ V base , proportional to the actual population, we have that W (t) satisfies the equatioṅ

that is W (t) = e −ψV base t . It implies that the number of vaccinated individuals ( W (t)) at time t is given by W (t) = 1 − e −ψV baset . Therefore, to vaccinate a fraction z of a population in the given time horizon T , it follows that ψ V base satisfies the equation Figure 3 : Feasibility region for vaccine reproduction number. The vaccine reproduction number R V is plotted as a function of vaccine efficacy (ε V ) and vaccination rate (ψ V ). Gray shaded region, corresponds to R V > 1. White region, denotes when R V < 1. Red region is biologically unfeasible.

Observe that in the calculation of ψ V base , it is considered all the population to be vaccinated. For our 197 study, vaccination is not applied to symptomatic infectious individuals. Therefore, Equation (9) represents 198 an approximation of our vaccination coverage (x coverage ) at constant rate.

199 Figure 4 shows the contour curves for R V as a function of the vaccine efficacy (ε V ) and vaccination rate order drives R V lower to 1. Further, a vaccine efficacy of 50 % or more is required so that, with an adequate 205 vaccination rate, the R V value can be reduced below one.

In the next section, the optimal control theory will be applied to propose optimal vaccination policies 207 that minimize the COVID-19 burden. 

In the remains of this manuscript, we use the following definitions. 

Our main idea is taking ψ V as Equation (9) and modulating it additively by a time function u V (t). We

amplifies or attenuates the constant vaccination rate ψ V . If m 1 ∈ (0, 1], then control signal u V (t) attenuates 216 the vaccination rate ψ V . Meanwhile, if m 2 > 0, then control signal u V (t) amplifies this vaccination rate. 217 We modify components equations corresponding to S, V , X in Equations (1) and (3) by 

Formally we define a controlled vaccination policy as follows.

Definition 2 (Controlled vaccination policy). Conforming the model in Equation (13) we say that

is a controlled vaccination policy (CVP). Then,

denote the number of doses at time t according to the modulated vaccination rate (ψ V + u V (t)). 221 We aim to obtain time-control functions u V (·) that hold natural modeling constraints-as a fixed bound 

Here, a S and a D are parameters related to the definition of the Years of Life Lost (YLL) due to prema- To describe vaccination coverage, we ask the terminal conditon 233 ϕ(x(T )) = X(T ),

That is, given the time horizon T , we set the vaccination coverage to 20 % or 50 % of the total population, 234 and the rest of final states free. Likewise, we impose the path constraint

to ensure that critical symptomatic cases will not overload healthcare services. Here κ denotes hospitalization 236 rate, and B is the load capacity of the health system.

Definition 3 (Admissible control vaccination policy). Let (x(·), u V (·)) be a pair satisfying the ODE 238 (13). Consider U[0, T ] as in (14). If

In other words, an admissible vaccination policy (AVP) is a CVP that satisfies the coverage and hospitaliza- 

is an optimal vaccination policy.

Remark 2. Optimal vaccination amplifies or attenuates the estimated baseline

to optimize functional J(·)-minimizing symptomatic incidence and death reported cases in DALYs and sat-252 isfying hospitalization occupancy and coverage constraints. 253 We aim to minimize the cost functional (15)-over an appropriated space-subject to the dynamics in 254 Equation (13), coverage related to the boundary condition (16) , and path constraints (17). We call this 255 kind of policies as optimal vaccination policies (OVP). That is, we seek vaccination policies that solve the 256 following problem.

Optimal Control Problem (OCP): Find the optimal vaccination rate (ψ V + u * V ) such that, Figure 5 illustrates the main ideas of the above discussion. Functional cosṫ 

where such that

The numerical analysis and design of transcript methods is a well-established and active research numerical derivatives computed by ADOL-C. 281 We provide in [41] a GitHub repository with all regarding Bocop sources. To perform the simulations corresponding to the scenarios presented in Table 5 , we fix the parameter 313 values as in Table 6 . Table 3 for the rest of parameters.

Parameters values To fix ideas, we display in Figures 6 and 7 DALYs. Figure 7 confirms this improvement by comparing disease dynamics with and without vaccination. 320 We observe CP and OVP reduce the symptomatic prevalence and the cumulative deaths. Note that despite 321 OVP requires less doses amount than CP at the time horizon (see Figure 6B ), the OVP performs better 322 than CP (see Figure 7 ). 

In February 2021, multiple vaccines have been rolled out to prevent COVID-19 disease. Table 3 , we explore the vaccines' consequences versus the counterfactual scenario with R 0 = 1.794 93.

Vaccine reproduction number, at the initial time of the optimal policy, for each vaccine results in R and two years, respectively. We display in Figure 10 the response of the vaccines vax 1 , vax 2 and vax 3 .

Since optimal vaccination policies are similar, Figure 10C suggests that the vaccine-induced immunity rate 349 is not determinant in the vaccination schedule design. Figure 11 shows a wide reduction of prevalence Figure 12A shows that, if natural immunity lasts one year, then burden of 360 COVID-19 falls until around 120 DALYs. We confirm this fall in the reduction of the symptomatic prevalence 361 cases and cumulative deaths, as displayed in Figure 13 . When natural immunity is 365 days, the benefit 362 in reducing symptomatic prevalence concerning a natural immunity of 90 days is at least 100 times. The natural immunization profiles, allowing the study of these variables' effects and developing optimal dose 383 administration schedules. Further, our optimal solution satisfies modeling constraints to obtain the desired 384 vaccination coverage and avoid health system overload. We optimized vaccination strategies by minimizing 385 the COVID-19 burden quantified in DALYs.

Note that the elderly population has been the most affected in terms of deaths due to the disease. rate in a more realistic bound. Further, we consider more detailed dynamics: our model includes exposed 398 and dead classes and differences between symptomatic infectious kinds. 399 We relate vaccination rate with the available vaccine stock and health system response-limited by [49], and the code is available in [41] .

In the estimation process, we use model in Equation (1) 

where, the index • runs over S, and A. The relaxation times T 1 , T 2 and T 3 corresponds to February 19, (Equation (4) ). That is, we count the number of new deaths 618 26 J o u r n a l P r e -p r o o f per day. In addition to the above, we assign prior probability distributions to each parameter and the initial conditions for the above mentioned classes. Finally, it is important to mention that our results were implemented considering a reduction in the 627 baseline transmission contact rates. However, implementation/breaking mitigation measures (among oth-628 ers) do not allow knowing with precision what the value of these parameters was at the beginning of the 629 vaccination campaign. For this reason, for each transmission contact rate, we join the ranges given in Ta-630 ble A.11. The resulting ranges and fitting curves with their respective bands are shown in Table 2 Now, continuing with the analysis of our model, it is easy to prove that the disease-free equilibrium is given by the point at X 0 ∈ Ω of the form

On the other hand, following [31], the next generation matrix for this model, evaluated at the disease where S * = (µ+ωV ) µ+ωV +ψV and V * = λ µ+ωV +ψV . Then, the spectral radius of K is

R V is the vaccine reproduction number. 

Efficacy of ChAdOx1 nCoV-19 (AZD1222) Vaccine

Interim framework for COVID-19 vaccine 497 allocation and distribution in the United States

Will an imperfect vaccine curtail the COVID-19 pandemic in the

Mathematical assessment of the roles of vaccination and non-503 pharmaceutical interventions on COVID-19 dynamics: a multigroup modeling approach, medRxiv

Abu-Raddad, Epidemio-506 logical Impact of SARS-CoV-2 Vaccination: Mathematical Modeling Analyses

IHME COVID-19 health service utilization forecasting team, Forecasting COVID-19 impact on hospital bed-days, ICU-509 days, ventilator-days and deaths by US state in the next 4 months

Forecasting hospital demand in metropolitan areas during the current 511

COVID-19 pandemic and estimates of lockdown-induced 2nd waves

Modeling behavioral change and COVID-19 con-514 tainment in Mexico: A trade-off between lockdown and compliance

Lifting mobility restrictions and the effect of 517 superspreading events on the short-term dynamics of COVID-19

Lockdown, relaxation, and acme period in COVID-19: A study of disease dynamics in 521

Staggered release policies for covid-19 control: Costs and benefits of relaxing restrictions by age and 523 risk

Mathematical assessment 525 of the impact of non-pharmaceutical interventions on curtailing the 2019 novel Coronavirus

Modeling shield immunity to reduce COVID-19 epidemic 529 spread

Optimal control applied to biological models

Optimal control of normalized SIMR models with vaccination and 533 treatment

Vaccine optimization for COVID-19, who to vaccinate first?

Model-informed COVID-19 537 vaccine prioritization strategies by age and serostatus

Optimal Dynamic Prioritization of Scarce COVID-19 Vaccines

Optimal Control of the COVID-19 Pandemic with non-pharmaceutical Interventions

Modeling the impact of non-pharmaceutical interventions on the dynamics of novel coronavirus 543 with optimal control analysis with a case study

Determination of an optimal control strategy for 546 vaccine administration in COVID-19 pandemic treatment

Reproduction numbers and sub-threshold endemic equilibria for compartmental 549 models of disease transmission

A vaccination model for trans-551 mission dynamics of influenza

An investigation of 553 transmission control measures during the first 50 days of the COVID-19 epidemic in China

World of Health Organization, WHO methods and data sources for global burden of disease estimates

Functional analysis, calculus of variations and optimal control

Datos abiertos

Jump: A modeling language for mathematical optimization

Bocop: an open source toolbox for optimal control

Source code for the manuscript Optimal Vac-569 cine Policies for COVID-19

Coronavirus Vaccine Tracker

Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine

Safety and efficacy of an rAd26 and rAd5 vector-579 based heterologous prime-boost COVID-19 vaccine: an interim analysis of a randomised

The Lancet

Collaborative data science

BNT162b2 mRNA Covid-19 Vaccine in a Nationwide Mass Vaccination Setting

Coronavirus reinfections: three questions scientists are asking

SARS-CoV-2 immunity: review and applications to phase 3 vaccine 588 candidates

Contemporary statistical inference for infectious 590 disease models using stan

CoV-2 mortality during the early stages of an epidemic: A modeling study in Hubei, China, and six regions in Europe

The authors whose names are listed immediately below certify that they have NO affi liations with or involvement in any organization or entity with any fi nancial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-fi nancial interest (such as personal or professional relationships, affi liations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.

The authors whose names are listed immediately below report the following details of affi liation or involvement in an organization or entity with a fi nancial or non-fi nancial interest in the subject matter or materials discussed in this manuscript. Please specify the nature of the confl ict on a separate sheet of paper if the space below is inadequate.Author names: