COVID-19 optimal vaccination policies: A modeling study on efficacy, natural and vaccine-induced immunity responses.

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.

Introduction

In late December 2019, the public health officials of Wuhan City, China, reported the emergence of a pneumonia illness of an unknown etiology (COVID-19). The virus (caused by SARS-CoV2) rapidly spread through many countries around the world, causing severe problems on their health systems. The first-line defense against the virus included several non-pharmaceutical interventions (NPIs) such as quarantine, isolation, and social distancing, being the main ones. Additionally, some drugs such as baricitinib, combined with remdesivir have been used in the U.S. to treat suspected or laboratory-confirmed COVID-19 hospital patients [1]. Despite such measures, the pandemic has not been controlled in most places. According to the WHO Dashboard [2], more than 111.5 million people have been infected, and 2.4 million individuals died around the world at the end of February 2021.

After COVID-19 became a Pandemic in March 2020, the World Health Organization (WHO) organized a group of specialists dedicated to immunization against COVID-19: the SAGE Working Group on COVID-19 vaccines [3]. In July 2020, WHO SAGE published in [3] the document entitled “Prioritized infectious disease and economic modeling questions”.

Here we attempt to answer some of the questions regarding the design of vaccination policies: when and how to administrate vaccine doses optimally.

Fortunately, the unprecedented efforts of many scientists have succeeded in developing vaccines to protect against COVID-19. When writing these sentences (February 2021), Mexico signed contracts that promise vaccine supply doses from Pfizer-BioNtech, AstraZeneca, CanSinoBio firms, and Gamaleya center. Besides, other vaccines probably get approval this year by the health authorities.

Although there are approved COVID-19 vaccine developments, its effective and fair administration implies enormous challenges. Health authorities will administrate more than one development. Each vaccine development requires different logistics, management, and administration training. For example, the Pfizer-BioNtech vaccine requires a high-tech cold chain for its transportation. In the same way, while most vaccines require two doses, others like Johnson & Johnson requires only one dose [4].

The significant barrier at the moment is to produce enough doses to vaccinate more than half the worldwide population in record time. Vaccines production demands complex processes, so we expect more situations like Pfizer’s and AstraZeneca’s announced delays. Consequently, health authorities would review and recalibrate its vaccination policies according to dynamic information.

In mathematical epidemiology, the modeling of vaccination policies reached vast and impressive advances. Some approaches range from deterministic to stochastic, discrete or continuous, and based on ordinary or partial derivatives [5], [6]. The new tendency also points towards statistical data analysis, optimization, and combinations of all mentioned techniques [7]. However, the previous studies provide limited insights into the particularities of SARS-CoV-2.

Despite new information emerges as the current pandemic evolves, the immunological responses of COVID-19 remains poorly understood. Essential and conclusive information about natural and vaccine-induced immunity remains under development [8], [9]. The variants of the SARS-CoV-2 would impact the vaccine efficacy—as in United Kingdom’s reported case with AstraZeneca [10]. Further, Johnson & Johnson reports in [4] that its vaccine efficacy differs across geographical regions.

Clearly, COVID-19 vaccination policies should endure complex and high uncertain issues. Thus, optimize the impact of the scarce vaccine supply is mandatory.

Our main objective is the formulation of optimal schedules for vaccine administration that ensure:

To cover a target population fraction in a fixed time horizon.

To minimize COVID-19 burden quantified in Disability-Adjusted Life Years (DALYs).

To preserve hospitalization occupancy below a required bound.

According to Toner et al. [11], the last influenza vaccination program may share some similarities with COVID-19 in the USA, but the latter demands new requirements. Further, health authorities would need to revise and adapt current vaccination policies as the pandemic evolves. Consequently, vaccination policies design and calibration should consider other NPIs in parallel.

Iboi et al. [12] model the effect of the combined NPIs strategies and vaccination and conclude that both are essential to mitigate the actual state of COVID-19 in the USA. Recently, [13] discusses this conducted strategy across groups differentiated according to mask-wearing. Makhoul et al. [14] analyze the impact of SARS-CoV-2 vaccination in age-stratified groups and vaccination strategies, according to reduction of susceptibility, infectiousness, duration of infection, and mitigation of severe cases.

Previous COVID-19 Kermack-McKendrick type models, as reported in [15], [16], aim to forecast the number of cases, deaths, or hospital occupancy. Other COVID-19 models explore the impact of NPIs as lockdown and exit strategies [17], [18], [19], [20], [21], or the effect of optimal serological tests as shielding [22].

To the best of our knowledge, the most popular way to model disease control strategies with optimal control is the so-called Lenhart’s approach (see, for example, [23]). Mainly, the literature report models that differ only in the compartmental model or functional cost. Although optimization problems in Engineering with mixed constraints—as boundary conditions or path restrictions—are routine, in Mathematical Epidemiology are uncommon [24].

The preprints [25], [26], [27] address the prioritization of the COVID-19 vaccine across five or more risk groups and optimize vaccine allocation. Their approaches face different prioritization policies according to vaccine efficacy and its availability. Buckner et al. [27] apply dynamic optimization with a time resolution of one month. Previous models based on dynamical optimization mainly focus on NPIs [28], [29] and just a few address optimal vaccination schedules [30].

According to real pharmaceutical profiles, our contribution employs dynamic optimization to compute optimal vaccination policies with a time-dependent vaccination rate. We model vaccination strategies as an optimal control problem with mixed restrictions. Our setup allows us to minimize the burden of COVID-19 in years of life lost and satisfying constraints of hospitalization occupancy, vaccination coverage, and time horizon of doses administration. Such strategies are aligned to the policies of the WHO strategic advisory group of experts (SAGE) on COVID-19 vaccination [3]. We run numerical experiments to explore hypothetical scenarios conforming to four main topics: optimal schedules, vaccine efficacy, natural and vaccine-induced immunity.

The present work explores hypothetical scenarios conforming to:

optimal vaccine administration schedules modulated by time-dependent vaccination rate

different vaccine profiles as efficacy and immunization periods

natural and induced-vaccine immunity responses

The present work helps to set the basis to investigate possible future scenarios with vaccination policies. For example, the implementation of optimal strategies when considering more than one vaccination dose. Particularly to explore the consequences in the vaccination policies when varying the inter-dose interval. Another implication is related to find optimal strategies for a second immunization program once the induced immunization period, of the first vaccination campaign, has expired.

Section 2 formulates our vaccination model. Section 3 discusses a vaccine reproduction number. In Section 4, we establish our optimal control problem. Section 5 displays numerical experiments regarding optimal vaccination policies. Section 6 presents our discussions. We end with a conclusions section.

Formulation of mathematical model

We use an extension of the classical Kermack–McKendrick model. Our formulation considers vaccination and vital dynamics. To this end, for fixed time , we split the total population , according to the following compartments: susceptible (), exposed (), symptomatic infectious (), asymptomatic infectious (), recovered (), dead () and vaccinated (). Our formulation requires the following hypotheses:

The vaccine is administered to all individuals exempting those with symptoms. Therefore, only individuals in the , , and classes are candidates for vaccination.

The vaccine only has effects on the susceptible individuals. Thus, susceptible individuals become vaccinated at rate.

The vaccine only protects against COVID-19.

Individuals get vaccinated only once during the epidemic.

Once an inoculated individual gets vaccine-induced immunity, returns to the class after a period of time (waning immunity period).

The vaccine is imperfect. A fraction of individuals in may become infected, with a lower probability than those in the class.

After a natural immunity period, the recovered population returns to the susceptible class.

Since we will explore disease dynamics that lasts from six months to one year, the model includes vital dynamics. We consider a constant population . Thus, we assume that birth and natural death rates are the same and represented by . All births lie into the class and all but class experience natural death. Class does not intervene in the transmission dynamics and counts reported deaths.

The infection dynamics are as follows: Susceptible individuals become infected, but not infectious, when in contact with infectious individuals and . Exposed individuals () remain in their class until they become infectious and move to either or . Individuals in class either die by disease complications or recover, whereas individuals in class move to the class after a some time. Finally, as the vaccine is considered imperfect, individuals in move to the class by interacting with infectious individuals and at a lower rate than individuals. Fig. 1 shows the compartmental model diagram which summarizes hypotheses mentioned above.

Fig. 1

Compartmental diagram of COVID-19 transmission dynamics which including vaccination dynamics. Here, there are seven different classes: Susceptible , exposed , symptomatic infected , asymptomatic infected , recovered , death and vaccinated individuals.

The model is given by the following ordinary differential equations system where the infection force is defined by Here, . For system in Eq. (1) all the variables are taken normalized by the constant total population . Therefore . Let The dynamics in Eq. (1) is positively-invariant on (see Appendix B). Additionally, the equations count the cumulative administered vaccines doses until time as the product . Also, the cumulative incidence of reported cases is given by .

Vaccine is administered to individuals in classes , , and . The amount of given vaccines is quantified by Eq. (3). However, as we assume the vaccine has preventive nature only, there is no change from classes , , and to due to vaccination.

The parameters of model Eq. (1) are described in Table 1.

Table 1

Parameters definition of system in Eq. (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 vaccine-induced 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

Calibration of baseline parameters

Multiple COVID-19 studies have shown important differences in transmission contact rates values across countries. For this reason, we establish a set of baseline parameter values for our geographic region of interest: the area constituted by Mexico-City and Mexico-State.

Mexico’s COVID-19 database provides detailed information on reported cases, hospitalized, ambulatory and deaths. Following the ideas of  [15], we consider COVID-19 confirmed deaths data to calibrate both 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 distribution with mean given by where represents the daily deaths incidence at the th day, and is the solution of the sixth equation of the system in Eq. (1) without vaccination dynamics at the th day. Fig. 2 shows fitting curves with their respective confidence bands. Estimation process considers data from February 19, 2020, to October 31, 2020. Like other studies [17], [18], it is considered perturbations on both transmission contact rates due to implementing or breaking mitigation measures. For more information about the parameter estimation process, see Appendix A. Table 2 summarizes our parameter calibration.

Fig. 2

Fitting death curve of the COVID-19 outbreak in Mexico-City and Mexico-State. Panel A shows new reported deaths per day. Panel B represents cumulative deaths per day. Reported deaths data are shown in blue bars from February 19, 2020, to October 31, 2020.

Table 2

Estimated range for some parameters of system in Eq. (1) without vaccination dynamics.

Parameter	Estimated range	Calibrated
β_S	[0.058 262,0.544 492]	0.363282
β_A	[0.101 754,0.441 215]	0.251521
p	[0.111 348,0.249 985]	0.1213

Vaccine reproduction number

In this section, we present a formulation for the vaccine reproduction number. Using the definitions of Van den Driessche and Watmough [31], the basic reproductive number for our dynamics without vaccination () and the vaccine reproduction number () are calculated for Eq. (1). Since can be rewritten in terms of , 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 represents the contribution of the symptomatic and asymptomatic infectious, respectively, to the spread of the disease.

On the other hand, the vaccine reproduction number for Eq. (1) is given by (see Appendix B) where, Let 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 as Note that . Observe that if , then . Otherwise, if and the following conditions hold then value is lower than one. Thus, for a waning rate given, there exist an adequate efficacy and vaccination rate to control the disease. However, if the conditions in (8) is not satisfied, it will not be possible to reduce the value of below 1.

To illustrate the aforementioned, Fig. 3 shows the regions where it is possible to reduce the value of . In this case, we set all the system parameters as given in Table 3 and fixed , leaving and free.

Fig. 3

Feasibility region for vaccine reproduction number. The vaccine reproduction number is plotted as a function of vaccine efficacy () and vaccination rate (). Gray shaded region, corresponds to . White region, denotes when . Red region is biologically unfeasible.

Table 3

Fixed parameters values of system in Eq. (1). The parameters corresponding to vaccination are established in each scenario.

Parameter	Value	95 %CI	Reference/Source
α_S	0.092 506 9	[0.043 233,0.181 159]	[17]
α_A	0.167 504	[0.086 355,0.198 807]	[17]
σ_E	0.196 078	[0.182 815,0.231 481]	[33]
β_S	0.363 282	a	Calibrated
β_A	0.251 521	a	Calibrated
p	0.1213	a	Calibrated
σ_R	0.002 739 73		Assumed
μ	0.000 039 138 9		Assumed
θ	0.11	b	From data

Confidence intervals for parameters , and in Table 2 and Appendix A.

Median estimated from Mexico data.

Baseline vaccination rate

Vaccination policies to reach a given coverage of a certain percentage of the population in a given period is of great importance. In this sense, we refer to this vaccination constant rate as the base vaccination rate, denoted by .

Let be the normalized unvaccinated population at time . If no individual has been vaccinated at , then . Assuming that we vaccinate individuals at a constant rate , proportional to the actual population, we have that satisfies the equation that is . It implies that the number of vaccinated individuals at time is given by . Therefore, to vaccinate a fraction of a population in the given time horizon , it follows that satisfies the equation thus

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

Fig. 4 shows the contour curves for as a function of the vaccine efficacy and vaccination rate . The green line corresponds to vaccination rate equal to 0.000 611. With this vaccination rate, it is not possible to reduce the value of below one. The intersection of the red lines correspond a vaccine efficacy equals to 0.8 and the corresponding vaccination rate such that . Note that this vaccination rate values below 0.007 implies greater than 1. Thus, vaccination rate has to be greater than 0.007 in order drives lower to 1. Further, a vaccine efficacy of 50 % or more is required so that, with an adequate vaccination rate, the value can be reduced below one.

Fig. 4

Contour plot of , as a function of vaccine efficacy () and vaccination rate () for the case where the average vaccine-induced immunity period is six month. Dark green line represents the value of , corresponding to a coverage and a time horizon days. Red lines show a scenario in which it is possible to reduce the value below one, considering a vaccine efficacy of 0.8 and a vaccination rate of 0.007.

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

Optimal vaccination policies

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

Constant Vaccination Policy

Consider the model in Eqs. (1), (3). A constant vaccination policy (CP) is a policy where the vaccination rate remains constant for all time . Thus the number of administered vaccine doses at time with this CP results

Our main idea is taking as Eq. (9) and modulating it additively by a time function . We impose that , , , then term amplifies or attenuates the constant vaccination rate . If , then control signal attenuates the vaccination rate . Meanwhile, if , then control signal amplifies this vaccination rate.

We modify components equations corresponding to , , in Eqs. (1), (3) by Then our controlled dynamics reads

Formally we define a controlled vaccination policy as follows.

Controlled Vaccination Policy

Conforming the model in Eq. (13) we say that is a controlled vaccination policy (CVP). Then, denote the number of doses at time according to the modulated vaccination rate .

We aim to obtain time-control functions that hold natural modeling constraints—as a fixed bound for hospitalized prevalence and coverage at the final time—and optimize a conveniently cost functional. To this end, we have to assure our optimal controlled model solution, so we consider the functional space

Let and control signal . Following the guidelines of WHO-SAGE modeling questions [3], we quantify the burden of COVID-19, according to the DALY indicator. Adapting DALY’s definition reported in [34], we optimize the number of years of life lost with a controlled vaccination policy. Our formulation calculates a minimum of the penalization functional Here, and are parameters related to the definition of the Years of Life Lost (YLL) due to premature mortality and the Years Lost due to Disability (YLD). We estimate as the average remaining life expectancy at the age of death, and according to the union of Mexico-City and Mexico-State data, we set . Parameter is the product of a disability weight (DW) and the average duration of cases until remission or death in years, that is, . Here we postulate the disability weight as the arithmetic average of disability weight regarding comorbidities reported in [35]. Our simulations employ . Thus, functional penalizes the pandemic burden—in Years of Life Lost—due to mortality or disability.

To describe vaccination coverage, we ask the terminal condition That is, given the time horizon , we set the vaccination coverage to 20 % or 50 % of the total population, 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 rate, and is the load capacity of the health system.

Admissible Control Vaccination Policy

Let be a pair satisfying the ODE (13). Consider as in (14). If

,

holds, then the CVP is admissible.

In other words, an admissible vaccination policy (AVP) is a CVP that satisfies the coverage and hospitalization constraints imposed on model (13). Furhter, if an AVP optimizes functional cost (15), then this AVP is an optimal vaccination policy (OVP). Formally we have the following definition.

Optimal Vaccination Policy

Let a pair that satisfies the ODE (13) such that (AC-1)–(AC-4) of Definition 3 holds. Let cost functional as in (15). If then is an optimal vaccination policy.

Optimal vaccination amplifies or attenuates the estimated baseline in an interval to optimize functional —minimizing symptomatic incidence and death reported cases in DALYs and satisfying hospitalization occupancy and coverage constraints.

We aim to minimize the cost functional (15)—over an appropriated space—subject to the dynamics in Eq. (13), coverage related to the boundary condition (16), and path constraints (17). We call this kind of policies as optimal vaccination policies (OVP). That is, we seek vaccination policies that solve the following problem.

Optimal Control Problem (OCP): Find the optimal vaccination rate such that, Fig. 5 illustrates the main ideas of the above discussion. Table 4 enclose parameter information of the functional cost and constraints.

Fig. 5

Compartmental diagram of COVID-19 transmission dynamics that includes optimal vaccination dynamics (19), penalization and a path constraint.

Table 4

Parameters regarding the constraints conditions and cost functional of the OCP (19).

Symbol	Description	Value	Ref
a_D	Penalization weight due to premature mortality (YLL) and estimated from Mexico-City an Mexico-State data	7.5years	[34], [37]
a_S	Penalization weight due to disability (YLD)	0.008 418 473years	[35]
x_coverage	Covering constraint at time horizon T		[3]
κ	Hospitalization rate	0.05	Estimated
B	Health service capacity in number of beds	9500	Estimated

Existence of solution to our (OCP) in Eq. (19) drops in the theory developed by Francis Clark [see e.g. 36, Thm. 23.11]. Since we aim to simulate hypothetical scenarios, we omit here a rigorous proof. Instead, we refer interested readers to [23] and the reference therein.

Numerical experiments

Methodology

We apply the so-called transcript method to solve our (OCP). This method transforms the underlying problem of optimizing functional governed by a differential equation into a finite-dimensional optimization problem with restrictions. To fix ideas, let , respectively denote state and control, and consider the optimal control problem Then, transcription methods transform this infinite-dimensional optimization problem into a finite dimension problem (NLP) via discretization of dynamics, state, and control. For example, if we employ the Euler method with a discretization of constant steps with size , then we can solve where , in the grid Let . Thus Eq. (20) defines a nonlinear programming problem on the discretized state and control variables of the form The numerical analysis and design of transcript methods is a well-established and active research numerical field. There is a vast literature about robust methods, and recently, implementations have been developed in vogue languages like Julia  [38], Python [39], and others.

Our simulations rely on the Bocop package  [40] to solve our OCP. Bocop is part of the development of the INRIA-Saclay initiative for open source optimal control toolbox and supported by the team Commands. BOCOP solves the NLP problem in Eq. (21) by the well-known software Ipopt and using sparse exact derivatives computed by ADOL-C.

We provide in [41] a GitHub repository with all regarding Bocop sources.

Hypothetical scenarios

Following the guidelines reported by the WHO Strategic Advisory Group of Experts (SAGE) on Immunization Working Group on COVID-19 Vaccines modeling questions presented in [3]. We simulate scenarios to illustrate vaccination policies’ response with a preventive vaccine. We aim to contrast the impact on the burden of COVID-19, according to

Optimal versus constant vaccination policies

Vaccine efficacy

Induced vaccine immunity

Natural immunity

We consider vaccine efficacies and doses compatible with the firms CansinoBio and Johnson & Johnson [42]. Although other firms such as Pfizer-BioNTech, Moderna, Astra-Zeneca, Gamaleya Research Institute, Sinovac Biotech, among others, require two doses [10], [43], [44], as a first approximation, our model simulates their total efficacy in a single dose. Further, since reinfection and vaccine-induced immunity parameters remain unavailable, we see it pertinent to explore the effect of plausible settings.

Optimal vaccination policy implies that number of doses per unit time described by mitigates the outbreak in optimal form, where optimal is defined in terms of function (see Eq. (15)), that is minimizing years of life lost in DALYs. Counterfactual scenarios implies . Constant vaccination policy means .

We assume positive prevalence of all epidemiological classes according to the following hypothesis:

The implemented initial conditions for our numerical experiments are hypothetical and do not reflect the actual data reported by Mexico-City and Mexico-State health authorities. The initial conditions are taken such that an outbreak is on its growth phase.

We suppose that around of 30 % of the population is under Lockdown and is enclosed along with recovered class . That is, encircle mainly children, senior, home office, people with low mobility and COVID-19 recovered individuals.

Our numerical results are of qualitative nature and not necessarily sustain forecasting or follows the actual profile of the underlying COVID-19 pandemic data.

Table 5 encloses a brief description and parameter values regarding each scenario. The reader can also access the web Chart Studio Graph of each figure regarding data and plotly [45] visual representation.

Table 5

Setup parameters for counterfactual and response scenarios. See Table 2 for the rest of parameters.

Simulation scene	Description	Set-up  (x_coverage,T,ε_V,ωV−1,σR−1)
(SCN-1)	Likening between optimal and constant vaccination policies.	(20  % , 180  days , 70  % , 730  days , lifelong )
(SCN-2)	Vaccine efficacy blow	(50  % , 365  days , 50%,70%,90%, 730  days , 180  days )
(SCN-3)	Vaccine-induced immunity period	(50  % , 365  days , 90  % , 180days,365days,730days, 365  days )
(SCN-4)	Natural immunity period	(50  % , 365  days , 90  % , 730  days , 90days,180days,365days)

To perform the simulations corresponding to the scenarios presented in Table 5, we fix the parameter values as in Table 6.

Table 6

Fixed parameters values of system in Eq. (19).

Parameters values
	(SCN-1)–(SCN4)
β_S, β_A, α_S, α_A, σ_E, μ, θ, p	Table 2
a_D, a_S, κ, B	Table 4
	(SCN-1)	(SCN-2)–(SCN-4)
ψ_V	0.001 239 69	0.001 899 03
u_min	−0.5ψ_V	−0.5ψ_V
u_max	5ψ_V	2.5ψ_V

Optimal versus constant vaccination policies: (SCN-1)

To fix ideas, we display in Fig. 6, Fig. 7 the counterfactual scenario regarding no intervention, constant vaccination policy (CP), and optimal vaccination policy (OVP) with a vaccine profile of efficacy , vaccine-induced immunity and a campaign for 20 % of coverage at 180 days. Fig. 6 suggests that the OVP improves CP response and counterfactual scenario given by the disease burden quantified in DALYs. Fig. 7 confirms this improvement by comparing disease dynamics with and without vaccination. We observe CP and OVP reduce the symptomatic prevalence and the cumulative deaths. Note that despite OVP requires less doses amount than CP at the time horizon (see Fig. 6B), the OVP performs better than CP (see Fig. 7).

Fig. 6

Effect of the vaccination policy on the burden COVID-19 for a 20% coverage at time horizon of half year. (A) Vaccination policies’ response regarding constant () and optimal () vaccination rates in the burden of COVID-19 quantified in DALYs. (B) Evolution of the vaccination covering according to each policy. (C) Vaccination schedule for each vaccination policy. Blue line corresponds to policies with constant vaccination rate 0.001 239 69. Green line corresponds to optimal vaccination policy. For counterfactual reference (panel A), black line represents the burden of COVID-19 without vaccination. See https://plotly.com/ MAAZ/366/ for plotly visualization and data.

Fig. 7

Effect of the vaccination policy on outbreak evolution. Optimal policy versus no vaccination (first column), constant policy versus no vaccination (second column) and optimal versus constant policy (third column). Upper row shows the symptomatic prevalence per 100 000 inhabitants. Lower row illustrates cumulative deaths (per 100 000 inhabitants). The shaded area represents the improvement of one policy over its corresponding. Data and web visualization in https://plotly.com/ MAAZ/474/.

Vaccine efficacy (SCN-2)

In February 2021, multiple vaccines have been rolled out to prevent COVID-19 disease. Table 7 summarizes some of these vaccine efficacies. To encompass these scenarios we focus on the 50 %, 70 % and 90 % vaccine efficacies.

Table 7

Vaccine efficacy of some of the approved developments for emergency use.

Developer	Vaccine name	Vaccine efficacy %, (95% CI)	Reference
Pfizer-BioNTech	BNT162b2	95(90.3 –97.6 )	[46]
Gamaleya Institute	Sputnik V	91.6(85.6 –95.2 )	[44]
Oxford University-AztraZeneca	AZD1222	74.6(41.6 -88.9 )	[10]
Johnson & Johnsona	Ad26.COV2.S	57 %, 66  % or 72  %	[4]
Sinovac Biotecha	CoronaVac	50.4 %	[42]

No available information about the confidence intervals.

Fig. 8, Fig. 9 display the optimal vaccination policy’s response according to the vaccine efficacies above mentioned. Fig. 8A displays COVID-19 burden in DALYs for the different vaccine efficacies and without vaccination dynamics. Clearly, the higher the vaccine efficacy, the lower the DALYs become.

Fig. 8

The response of COVID-19 burden on vaccine efficacy. (A) COVID-19 burden response quantified in DALYs per 100 000 inhabitants to vaccines with efficacy of 50 % (blue), 70 % (red) and 90 %(green). (B) Coverage evolution to reach 50 % of the total population vaccinated. (C) Optimal vaccination doses schedule according to the different efficacies. See https://plotly.com/ MAAZ/358/ for visualization and data.

Fig. 9

Effect of vaccine-efficacy on hospital occupancy and on the number of saved lives compared to no vaccination dynamics. Vaccine-efficacy of 50 % (blue), 70 % (red) and 90 % (green). Dotted red line represents 90 % of hospital beds per 100 000 inhabitants. See https://plotly.com/ MAAZ/470/ for data and visualization.

Figs. 8B–C also illustrate the effect of vaccine-efficacies on coverage and optimal vaccination policy, respectively. We observe that vaccine efficacy influences the design of an optimal policy. For this scenario the higher vaccine efficacy, the earlier is its massive implementation. According to 50 % coverage at a time horizon of 1 year, Fig. 9 displays an improvement from 1.5 to 3 times regarding the counterfactual scenario for the cumulative deaths at time horizon. Also, there is a notorious improvement in the hospital occupancy. For the no vaccination dynamics, the hospital occupancy exceeds the 90% of beds (dashed red line in Fig. 9). In contrast, for the considered controlled dynamics, the higher hospital occupancy lowers to about 63% of beds (see Fig. 9A).

Vaccine-induced immunity (SCN-3)

Vaccine response and its induced immunity are strongly related on the mitigation prevalence. However, vaccine-induced immunity period remains poorly understood [8]. Here, we contrast three vaccines with different induced-immunity. Let denote by , , vaccines with an induced-immunity period of a half, one, and two years, respectively, and common efficacy of 90 %. Consider a vaccine campaign of time horizon of one year and 50 % coverage. Taking the same initial conditions, and parameters values, as in Table 2, we explore the vaccines’ consequences versus the counterfactual scenario with . Vaccine reproduction number, at the initial time of the optimal policy, for each vaccine results in , , for vaccine-induced immunity periods of a half, one, and two years, respectively. We display in Fig. 10 the response of the vaccines , and . Since optimal vaccination policies are similar, Fig. 10C suggests that the vaccine-induced immunity rate is not determinant in the vaccination schedule design. Fig. 11 shows a wide reduction of prevalence and cumulative deaths concerning the uncontrolled outbreak. This reduction is consistent with the vaccine reproductive number values corresponding to each vaccine-induced immunity. Thus, vaccine reproductive number values highlight the importance of having an adequate vaccine profile and adequate vaccine schedule to control the epidemic outbreak.

Fig. 10

Vaccine-induced immunity effect on the COVID-19 burden. (A) Effect on the burden of COVID-19 quantified in DALYs per 100 000 inhabitants due to vaccine-induced immunity of 180 days (green), 365 days (red) and 730 days (blue). (B) Coverage evolution to reach 50 % of the total population vaccinated. (C) Optimal vaccination doses schedule according to the different vaccine-induced immunity periods. Visualization and data in https://plotly.com/ MAAZ/407/.

Fig. 11

Effect of vaccine-induced immunity on mitigation and saved lives of COVID-19 outbreak. Upper row shows the improvement of optimal vaccination policy over no vaccination dynamics on the symptomatic prevalence per 100 000 inhabitants. Lower row shows a comparison of the saved lives (per 100 000 inhabitants) between optimal policies and no vaccination dynamics. Optimal policy with vaccine-induced immunity of a half year versus no vaccination dynamics (first column), optimal policy with vaccine-induced immunity of a year versus no vaccination dynamics (second column), and optimal policy with vaccine-induced immunity of two years versus no vaccination dynamics (third column). See https://plotly.com/ MAAZ/465/ for data and visualization.

Natural immunity hypothesis (SCN-4)

“Reinfections raise questions about long-term immunity to COVID-19 and the prospects for a vaccine”, reported Heidi Ledford in [47]. Following this line, we display in Fig. 12, Fig. 13 the vaccine’s response with 90 % efficacy and contrasting with natural immunity periods of 90 days, 180 days, and 365 days. Here, the adjective “natural” denotes the immunity that an individual develops after recovering from a previous bout of COVID-19 without vaccination. Fig. 12A shows that, if natural immunity lasts one year, then burden of COVID-19 falls until around 120 DALYs. We confirm this fall in the reduction of the symptomatic prevalence cases and cumulative deaths, as displayed in Fig. 13. When natural immunity is 365 days, the benefit in reducing symptomatic prevalence concerning a natural immunity of 90 days is at least 100 times. The number of deaths with a natural immunity of 90 days reaches 845 cases per 100 000 inhabitants, in contrast, to 206 cases when natural immunity is 365 days. Thus, this scenario suggests that natural immunity plays an essential role in controlling the outbreak, consistent with the conclusions reported in [8].

Fig. 12

Effect of natural immunity on the burden of COVID-19. (A) Effect on the burden of COVID-19 quantified in DALYs per 100,000 inhabitants due to natural immunity of 90 days (red), 180 days (blue) and 365 days (green). (B) Coverage evolution to reach 50 % of the total population vaccinated. (C) Optimal vaccination doses schedule according to the different natural immunities. https://plotly.com/ MAAZ/402/.

Fig. 13

Effect of natural immunity on COVID-19 outbreak. (A) Effect of natural immunity over the symptomatic prevalence per 100 000 inhabitants. (B) Improvement of natural immunity over the cumulative deaths per 100 000 inhabitants. Plotly visualization and data in https://plotly.com/ MAAZ/451/.

Discussion

In February 2021, at least four vaccination developments support the primary pharmaceutical measure to recover life’s style before the pandemic. Although the leading pharmaceutics firms report high vaccine efficacy, its implementation implies new challenges as its distribution, stocks, politics, logistic, among others, emerge. A fair distribution and application strategy is imperative to manage the available resources. Despite that vaccine-immunization response remains under study, vaccination campaigns recently started. Thus evaluating the impact of different vaccine profiles and natural immune responses would be crucial to calibrate vaccination policies. Further, since the vaccine doses are scarce at the beginning of the COVID-19 vaccination campaigns optimize their administration subject to its availability, minimizing deaths and preserving health systems is mandatory.

We aim to compare different vaccine profiles in computing optimal vaccination policies for COVID-19 and evaluate the impact of hypothetical reinfection and immunization responses. Our simulations suggest that an optimal strategy design is influenced by vaccine profile and natural immunity period. Likewise, we observe that natural and vaccine-induced immunity would play an essential role in reducing COVID-19 disease mortality and prevalence levels. Poland et al. [48] stress that understanding immune responses to SARS-CoV-2 are crucial to developing vaccines. Our simulation faces different vaccine-induced and natural immunization profiles, allowing the study of these variables’ effects and developing optimal dose administration schedules. Further, our optimal solution satisfies modeling constraints to obtain the desired vaccination coverage and avoid health system overload. We optimized vaccination strategies by minimizing 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. Thus, the average remaining expectancy of lives could be explained by this population. It is well know that stratification becomes important when vaccination strategies prioritize risk groups to reduce burden. Although our model lacks an age-stratified structure, these results would apply to situations where the population has homogeneous characteristics. This kind of population’s designs arises, for example, in places like factories, offices, or even schools. The question “who vaccinate first¿‘ has been faced as an optimal allocation problem in [25], [26], [27]. Our optimal solutions have a complementary objective—while we answer when and how to administrate vaccine doses, the mentioned references answer who vaccinate first.

In [30] the authors answered when to vaccinate optimally. They used a signal that switches between not vaccination and immunization of all susceptible individuals. Since the vaccine is scarce, this control signal would be unrealistic in some situations. Our results improve this strategy by modulating a base vaccination rate in a more realistic bound. Further, we consider more detailed dynamics: our model includes exposed and dead classes and differences between symptomatic infectious kinds.

We relate vaccination rate with the available vaccine stock and health system response—limited by logistics, management, politics, and other issues related to each vaccine development. So, our contribution can help decision-makers to design vaccination schedules of homogeneous populations.

When writing this article, the vaccine-induced immunity period and SARS-CoV-2 transmission capacity from a vaccinated individual remain unexplained [48]. The conclusive answer inevitably will drive the development of COVID-19 vaccines and their policies. Experts estimate that around 7 800 billion immunization doses are required to reach herd immunity. This number, with the approved developments, is far away to accomplish in one year. However, we expect that more vaccines finish their third trial phase and also get approval. Thus, vaccination policies also have to contemplate two o more vaccine profiles, perhaps with different efficacy o number of recommended doses.

These issues also have to balance the health services and other resources optimally.

Conclusions

The absence of proven highly effective treatments against COVID-19 has led the world’s population to pin their hopes on vaccines. However, due to the limited stock of vaccines, the question arises as to which is the best vaccination strategy. Thus, different researchers have sought to answer this question from different perspectives. In the present work, we address this question through an optimal control approach. Our results show the importance of how and when to vaccinate on reducing disease levels.

The scarce information regarding immunity periods and the wide efficacy range of the approved vaccines highlights the importance of exploring scenarios with different vaccination profiles. This study observes the importance of knowing the vaccine profile and the natural immunity period to design an optimal vaccination strategy. For example, SCN-2 and SCN-4 scenarios show that an optimal vaccination strategy does not necessarily imply massive vaccination at the beginning. For the former, the vaccine efficacy level (Fig. 8) is a responsible factor to determine the beginning time of the massive vaccination strategy whereas for the latter is the natural immunity period (Fig. 12).

Nowadays, several countries have problems with vaccine stocks. In addition, the uncertainty regarding induced immunity and natural immunity has forced governments to implement vaccination programs aimed at specific populations such as medical personnel, elderly population, teachers, etc., with similar characteristics to each other. Our work can be useful in decision-making for the design of vaccination programs targeted to this populations types, since it considers stock constraints and explores different vaccination profiles and natural immunity to design optimal vaccination strategies.

At the date of writing this manuscript, several questions about vaccines and COVID-19 disease have not been addressed thoroughly. How often shall we need to get vaccinated, once in a lifetime or yearly? Do all the vaccines’ acquired immunity last equally? If the latter is false, what are the consequences of implementing campaigns with different vaccines? Would the vaccines be produced on a much larger scale than now? If the latter is true, when, how, and to whom should vaccines be given to obtaining greater reductions on incidence levels? What effect will the new COVID-19 variants have on the vaccines developed? Would the vaccines’ effectiveness significantly decrease? If the latter is true, then would be necessary more doses for each vaccine? All these questions will eventually become daily concerns in society as time passes by. Thus, optimal vaccination strategies take more relevance to address some of these questions.

Motivated by the above, as future work, we analyze extending our study when considering a stratified population problem and more than one vaccine with different efficacies and number of doses. We aim to answer some of the previously raised questions over non-homogeneous populations.

CRediT authorship contribution statement

Manuel Adrian Acuña-Zegarra: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Resources, Visualization, Project administration, Writing - original draft, Writing - review & editing. Saúl Díaz-Infante: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Visualization, Supervision, Writing - original draft, Writing - review & editing. David Baca-Carrasco: Conceptualization, Methodology, Formal analysis, Writing - original draft, Writing - review & editing. Daniel Olmos-Liceaga: Conceptualization, Methodology, Formal analysis, Writing - original draft, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table A.8

Fixed parameters values of system in Eq. (A.1).

Parameter	Value	References
σE−1	5.1days	[33]
αS−1	5.97days	[17]
αA−1	10.81days	[17]
σR−1	365days	Assumed
μ^−1	70years	Assumed

Table A.9

Prior distributions for each parameter and initial conditions.

Parameter	Prior distribution	Initial condition	Prior distribution
βˆ_S	N(0.5,0.1)	E(0)	U(1,15)
βˆ_A	N(0.5,0.1)	I_S(0)	U(1,5)
p	U(0.05,0.25)	I_A(0)	U(1,15)
ϵ_1	U(0.05,1)
ϵ_2	U(0.05,1)

Table A.10

Estimated range for some parameters of system in Eq. (A.1).

Parameter	Estimated range
βˆ_S	[0.243439,0.814272]
βˆ_A	[0.416429,0.774116]
p	[0.111348,0.249985]
ϵ_1	[0.503394,0.793054]
ϵ_2	[0.3371,0.412239]

Table A.11

Estimated range for some parameters of system in Eq. (A.1).

Parameter	Estimated range
ϵ_1βˆ_S	[0.147767,0.544492]
ϵ_1βˆ_A	[0.294708,0.441215]
ϵ_2ϵ_1βˆ_S	[0.058262,0.193831]
ϵ_2ϵ_1βˆ_A	[0.101754,0.171101]