key: cord-0792649-uep62a4g
authors: Acuna-Zegarra, M. A.; Diaz-Infante, S.; Baca-Carrasco, D.; Olmos Liceaga, D.
title: COVID-19 optimal vaccination policies: a modeling study on efficacy, natural and vaccine-induced immunity responses
date: 2020-11-20
journal: nan
DOI: 10.1101/2020.11.19.20235176
sha: f3f82d3e23717807bb0fed26b2709663e8694522
doc_id: 792649
cord_uid: uep62a4g

At the date, Europe and part of North America face the second wave of COVID-19, causing more than 1 300 000 deaths worldwide. Humanity lacks successful treatments, and a sustainable solution is an effective vaccine. Pfizer and the Russian Gamaleya Institute report that its vaccines reach more than 90 % efficacy in a recent press release. If third stage trial results favorable, pharmaceutical firms estimate big scale production of its vaccine candidates around the first 2021 quarter and the World Health organization fix as objective, vaccinate 20 % of the whole population at the final of 2021. However, since COVID-19 is new to our knowledge, vaccine efficacy and induced-immunity responses remain poorly understood. There are great expectations, but few think the first vaccines will be fully protective. Instead, they may reduce the severity of illness, reducing hospitalization and death cases. Further, logistic supply, economic and political implications impose a set of grand challenges to develop vaccination policies. For this reason, health decision-makers require tools to evaluate hypothetical scenarios and evaluate admissible responses. Our contribution answers questions in this direction. According to the WHO Strategic Advisory Group of Experts on Immunization Working Group on COVID-19 Vaccines, we formulate an optimal controlled model to describe vaccination policies that minimize the burden of COVID-19 quantified by the number of disability-adjusted years of life lost. Additionally, we analyze the reproductive vaccination number according to vaccination profiles depending on coverage, efficacy, horizon time, and vaccination rate. We explore scenarios regarding efficacy, coverage, vaccine-induced immunity, and natural immunity via numerical simulation. Our results suggest that response regarding vaccine-induced immunity and natural immunity would play a dominant role in the vaccination policy design. Likewise, the vaccine efficacy would influence the time of intensifying the number of doses in the vaccination policy.

In late December 2019, a new virus's appearance is reported in Wuhan City, Hubei Province, China. Called SARS-CoV2, it is the virus that causes the 2019 coronavirus disease and that, very quickly since its appearance, has spread throughout much of the world, causing severe problems to health systems of all the countries in which it is present [1] .

In this section, we formulate our baseline mathematical model, which additionally to the transmission dynamics, includes vaccination. In order to build our model, we follow the classical Kermack-McKendrick approach. Figure 1 shows the compartmental diagram of our mathematical model. Figure 1 : Compartmental diagram of COVID-19 transmission dynamics which including vaccination dynamics. Here, there are seven different classes: Susceptible (S), exposed (E), symptomatic infected (I S ), asymptomatic infected (I A ), recovered (R), death (D) and vaccinated (V ) individuals. It is important to mention that I S represents the proportion of symptomatic individuals who will later report to some health medical center.

Information about reinfection dynamics on COVID-19 disease is unclear to date. However, to explore some scenarios related to this dynamic, we assume that reinfection is possible after a period of time. On the other hand, we model the vaccination process considering some assumptions: i) Vaccination is applied to all the alive individuals except those in the symptomatic class. In this situation, vaccines are applied indiscriminately over individuals on the S, E, I A , and R classes; ii) the vaccine has preventive nature, that is, only reflected in the susceptible individuals (S); iii) people will only get one vaccine during the campaign, and iv) vaccines do not necessarily have a hundred percent of effectivity, which implies that some vaccinated people can get the disease. We denote the effectivity rate by . Based on these assumptions our model becomes

whereN (t) = S(t) + E(t) + I S (t) + I A (t) + R(t) + V (t) and N =N + D. Additionally, we include the 3 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint equations X (t) = λ V (S + E + I A + R)

where X(t) and Y I S (t) represent the cumulative doses at time t, and the cumulative incidence at time t, respectively. Parameters description of system in Equation (1) is provided in Table 1 .

Parameter Description µ Natural death rate β S Infection rate between susceptible and symptomatic infected β A Infection rate between susceptible and asymptomatic infected λ V Vaccination rate δ −1 

It is now necessary to define a set of baseline parameter values to explore some scenarios of interest. In the present work, we consider Mexico City plus Mexico state as our study region and use COVID-19 data to estimate some parameter values. We follow a Bayesian approach to address this problem. We use a negative binomial distribution as an observation model, in which the mean parameter is given by

where I SA (k) represents the incidence per day of infected symptomatic individuals at the k − th day. The parameter estimation splits into two-stage: before and after mitigation measures were implemented, and we only focus on the early phase of the COVID-19 outbreak. Figure 2 shows fitting curves with their respective confidence bands for both stages. Here, we observe that our estimations follow the growth profile of the epidemic curve. For more information about the parameter estimation process, see Appendix A. available, then we assume that our scenarios start on the growth stage of a second outbreak (see Figure 3 ) with the objective of starting under a plausible scenario. Thus, for numerical results, initial conditions 4 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.19.20235176 doi: medRxiv preprint 

In this section, we use the estimated parameters from the previous section to analyze model in Equation (1) and the effect of the vaccine. Some plausible scenarios are presented, depending on the effectiveness of the vaccine, as well as the rate of vaccination.

In the first instance, it has been shown that the interest region of the state variables of the system is positively invariant and the proof can be found in Appendix B.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. Now, an important concept in the analysis of the spread of diseases is the basic reproductive number, defined as the number of secondary infections produced by a typical infected individual, throughout his infectious period, when in contact with a totally susceptible population. Following ideas of [25] , the basic reproduction number for system in Equation (1) is (see Appendix B)

Note that each sum of R V represents the contribution of the symptomatic and asymptomatic infected, respectively, to the spread of the disease. Following the ideas of Alexander et. al. [7] , expression for R V can be rewritten as

where

is the basic reproduction number of the system without vaccine. Note that 1 − λ V (µ+δ V +λ V ) < 1. Therefore, this factor, which enclose the parameters corresponding to the application of the vaccine, will allow us to modulate the value of R 0 . In the first instance, if R 0 < 1, then R V < 1. But, if R 0 > 1, we ask if the application of the vaccine can lower R V value below 1. In this sense, it is easy to prove that, if

for > 1 − (1/R 0 ), it is possible to reduce the value of R V below one. That is, there is a region in the parameter space in which it is possible to reduce the value of R V below one, considering adequate efficacy, vaccination rate and duration of the effect of the vaccine. However, if the Inequality (6) is not satisfied, it will not be possible to reduce the value of R V below 1.

To illustrate the aforementioned, Figure 4 shows the regions where it is possible to reduce the value of R V . In this case, we set all the system parameters as given in Table 3 and with δ V = 1/180, leaving and λ V free.

Parameter Value β S 0.363 282 β A 0.251 521 α S 0.092 506 9 α A 0.167 504 δ E 0.196 078 δ R 0.002 739 73 µ 0.000 039 138 9 θ 0.11 p 0.1213 Table 3 : Fixed parameters values of system in Equation (1) . The parameters corresponding to vaccination are established in each scenario studied. See Table 2 and Appendix A for more details. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.19.20235176 doi: medRxiv preprint Figure 4 : Vaccine efficacy versus vaccination rate feasibility. In the gray shaded region R V > 1 and in the white region R V < 1. Note that, for our scenario, with a 50 percent vaccine efficacy and an adequate vaccination rate, it is possible to reduce the R V value below one. The orange region is unfeasible.

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 λ V base .

Let W (t) be the normalized unvaccinated population at time t, we consider that at t = 0 no person has been vaccinated, which implies W (0) = 1. Then, by assuming that we vaccinate individuals at a constant rate λ V base proportional to the actual population, we have that W (t) satisfies the equatioṅ

or W (t) = e −λ V base t . It implies that the number of vaccinated individuals at time t is given by V (t) = 1 − e −λ V base t . Then, if we look to reach a coverage x coverage at a horizon time T , it follows that λ V base satisfies the equation

Observe that in the calculation of λ V base , it is considered all the population to be vaccinated. Vaccination is not applied to infective symptomatic individuals. Therefore, Equation (7) represents an approximation of the vaccination scheme at constant rate. Figure 5 , shows the contour curves for R 0 considering it as a function of the efficacy of the vaccine ( ) and of the vaccination rate (δ V ), considering an immunity period induced by the vaccine of half year . The blue line, correspond to the values of λ V base and we can see that with this vaccination rate, no matter how effective the vaccine is, it is not possible to reduce the value of R V below one. Purple lines show a scenario in which it is possible to reduce the R V value below one, considering a vaccine efficacy of 0.8 and a vaccination rate of 0.7. The figure shows plausible combinations of and λ V values in order to reduce the value of R V below one.

Note that a vaccine efficacy of 50 % or more is required so that, with an adequate 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 dynamics that minimize the number of cases of symptomatic infection and deaths due to the disease. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; 

According to dynamics in Equations (1) and (2), we modulate the vaccination rate by a time-dependent control signal u V (t) to achieve an imposed vaccine coverage. That is, according to the components S, V , X of Equations (1) and (2), we modulate the vaccination rate λ V by an additive control u V (t). Thus, we modify components equations related to S, V , X as

To assure the solution of our controlled model we consider the functional space

such that u V (·) bounded and piecewise continuous} . 

8 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ;

Here, a S and a D are parameters related to the definition of the unit of the Years of Life Lost (YLL) due to premature mortality and the Years Lost due to Disability (YLD). We estimate a D as the mean life expectancy at the age of death, and according to Mexico City+Mexico state data, we handle a D = 7.5 years. Parameter a S is the product of a disability weight (DW) and the average duration of cases until remission or death in years, that is, a S = DW ×α −1 S . Here we postulate the disability weight as the arithmetic average of disability weight regarding comorbidities reported in [26] . Thus, our simulations employ a S = 0.008 418 473 years. Thus, functional J penalizes the pandemic burden-in Years of Life Lost-due to mortality or disability. We display in Table 4 parameters regarding the optimal controlled model.

Since we aim to simulate vaccination policies contra factual scenarios following the SAGE modeling guidelines reported in [23] , we impose the vaccination counter state's horizon time condition X(T )

x(T ) = (·, ·, ·, ·, ·, X(T )) ∈ Ω,

Thus, given the time horizon T , we set the last fraction of vaccinated population corresponds to 20% or 50%, and the rest of final states are free. We also impose the path constraint

to ensure that critical symptomatic cases will not overload healthcare services. Here κ denotes hospitalization rate, and B is the load capacity of a health system. We illustrate the main ideas of the above discussion in Figure 6 . Given a fixed time horizon and vaccine with efficacy , we estimate the constant vaccination rate λ V as the solution of Equation (7) . That is, λ V estimates the constant rate to cover-with a vaccine dose per individual-population fraction x coverage in time horizon T . Thus, according to this vaccination rate, we postulate a policy u V that modulates vaccination rate according to λ V as a baseline.

Then, optimal vaccination amplifies or attenuates this estimated baseline λ V in an interval [λ min V , λ max V ] to optimize functional J(·)-minimizing symptomatic, death reported cases and optimizing resources. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ;

We aim to minimize the cost functional (9)-over an appropriated space-subject to the dynamics in Equations (1) and (8), boundary conditions related to (10) , and path constrain (11) . That is, we seek for vaccination policies u V (·), that solves the following optimal control problem (OCP) Existence of solution to our (OCP) in Equation (12) drops in the theory developed by Francis Clark [see e.g. 27, Thm. 23.11] . Since our aim is the simulation of hypothetical scenarios, we omit here a rigorus proof, instead we refer interested readers to [28, 29] is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.19.20235176 doi: medRxiv preprint

We apply the so-called transcript method to solve our (OCP). This schemes' main idea consists of transforming the underlying problem of optimizing functional governed by a differential equation into a finite-dimensional optimization problem with restrictions. To fix ideas, let x, u denote state and control and consider the optimal control problem min J(x(·), u(·)) = g 0 (T, x(T ))

Functional cosṫ

Boundary conditions.

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 N constant steps with size h, then we can solve

where 

The numerical analysis and design of transcript methods is a well established and active research numerical field. There is a baste literature about robust methods and resently it apperars implementations in vogue languages like Julia [31, 32] , Python [33] , Matlab [34] , and others. We refer the reader to [35, 36] for a more systematic discussion. Our simulations rely on the Bocop package [37, 38] 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 Equation (14) by the well known software Ipopt and using sparse exact derivatives computed by ADOL-C.

We provide in [39] a GitHub repository with all regarding R and Bocop sources for the sake of reproductivity. This repository also encloses data sources and python code to reproduce all reported figures.

We follow the guidelines reported by the WHO Strategic Advisory Group of Experts (SAGE) on Immunization Working Group on COVID-19 Vaccines modeling questions presented in [23] . According to this SAGE's document, we simulate scenarios to illustrate vaccination policies' response with a preventive vaccine. We aim to contrast the impact of the burden of COVID-19 mitigation regarding 11 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.19.20235176 doi: medRxiv preprint (SCN-1) Optimal versus constant vaccination policies (SCN-2) Vaccine efficacy (SCN-3) Induced vaccine immunity (SCN-4) Natural immunity

We consider vaccine profiles-efficacy and induced vaccine immunity in concordance with the expected but still unconfirmed data-from the firms Cansino-Biologics, Astra Zeneca, and Pfizer. Further, since reinfection and induced vaccine immunity parameters remain unavailable, we see pertinent explore the effect of plausible settings. Remark 1. We mean an optimal vaccination policy as the number of doses per unit time described by

Counterfactual scenarios implies u V (t) = 0. Without vaccination scenarios we means λ V = u V (t) = 0. Table 3 for the rest of parameters.

To perform the simulations corresponding to the scenarios presented in Table 5 , we fix the parameter values as in Table 6. 12 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. 

To fix ideas, we display in Figures 7 and 8 the counterfactual scenario regarding no intervention, constant vaccination policy (CP), and optimal vaccination policy (OP) with a vaccine profile of efficacy = 70 %, and induced immunity δ −1 V = 180 days and over a campaign for 20 % of coverage at 180 days. Figure 7 suggests that the OP improves CP vaccination policy response according to the burden disease due to mortality, morbidity, and coverage time. Figure 8 confirms this improvement by comparing symptomatic reported cases, saving lives, and the disease dynamics without vaccination. Although both campaigns use the same number of vaccine doses and the same vaccine profile, we observe that OP implies fewer deaths and symptomatic cases. Figure 7 , shows a scenario where R 0 > 1. Although R V remains below but close to R 0 , vaccination reproductive number R V , could explain vaccination response with constant policies according to the mitigation factor

That is, disease mitigation is strongly related to vaccine efficiency and vaccination rate λ V . Further, given a dynamic with not vaccine intervention and R 0 > 1, R V suggests a minimal vaccination rate to drive this dynamic to the disease-free state but subject to vaccines with particular efficacy.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; 

In the author's words of press releases 14 .

CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ;  "Pfizer says early analysis shows its Covid-19 vaccine is more than 90 % effective [41] ." "Russia says its Sputnik V COVID-19 vaccine is 92 % effective [42] ." "Moderna's coronavirus vaccine is 94.5 % effective, according to company data [43] ." This is a game changer fact-FDA would accept a vaccine of 50 %. Following this idea, Figures 9 and 10 display the response of the optimal vaccination policy according to three vaccines with different efficacy. Figure 9 -A displays the burden COVID-19 in DALYs for vaccines with efficacy of 50 %, 70 %, 90 %. According with a time horizon of 1 year and coverage of 50 %, Figure 10 displays an improvement of at least three times in the prevalence of symptomatic cases and saved lives concerning the uncontrolled outbreak. Figure 10 reflects this effect in the prevalence of symptomatic cases (A) and in the number of saved lives shaded by the translucent green color (B). is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; 

Vaccine response is also strongly related to its induced immunity -parameter that remains poorly understood [44] . Here we contrast two vaccines with different induced immunity. Let denote by vax 1 , vax 2 vaccines with an induced-immunity capacity of one and two years, respectively, and common efficacy of 70 %. Consider a vaccine camping of time horizon of one year and 50 % coverage. Taking the same dynamics parameters, that is initial conditions, and base line parameters as in Table 3 we explore a contra factual scenario with an uncontrolled outbreak of R 0 = 1.794 93 and two controlled dynamics according to vaccines vax 1 , vax 2 . Thus, according to these immunity parameters and factor defined in Equation (15) = 0.867 56 for vaccine immunities periods of one and two years. We display in Figure 11 the response of the vaccines vax 1 and vax 2 . Since in this set time horizon is of one year, the optimal policies follows similar schedules and implies similar gains in the number of years of life lost. Despite this similarities, Figure 12 displays in panel A a dramatically gain respect to the uncontrolled outbreak-since R

[vax1] V is near to one, prevalence fall-down more than five times and because R

[vax2] V is lest that one, prevalence of symptomatic cases tends to zero with damped oscillations. Figure 12 also endorses this gain, note that saved lives is represented by shaded region with translucent and overlapped red blue colors.

. CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ;

"Second infections raise questions about long-term immunity to COVID-19 and the prospects for a vaccine", reported Heidi Ledford in [45] . Following this line, we display in Figures 13 and 14 the vaccine's response with 90 % efficacy and contrasting with natural immunity periods of 90 days, 180 days, 365 days. Here, the adjective "natural" denotes the immunity that an individual develops after recovering from a previous bout of COVID-19. When natural immunity lasts one year, the burden of COVID-19 fall-down until around 120 DALYs. We confirm this behavior in the prevalence of symptomatic cases and cumulative deaths, as displayed in Figure 14 . When natural immunity is 365 days the gain in mitigation concerning a natural immunity of 90 days is at least and 100 times, while the number of deaths with a natural immunity of 90 days reach 845 cases per 100 000 inhabitants, in contrast, of 206 when natural immunity is 365 days. Thus, this simulation suggests that natural immunity plays a vital role in the controlled outbreak's behavior, which is consistent with the conclusions reported in [44] . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.19.20235176 doi: medRxiv preprint 

At the date of writing this article, humankind lacks strategies to eradicate COVID-19. Although NPIs implemented in most countries prevent citizens from being infected, these strategies leave them susceptible-people can not develop immunity to face futures waves. Thus, vaccination becomes the primary pharmaceutical measure to recover life's style before the pandemic. However, this vaccine has to be effective and well implemented in global vaccination programs. Thus new challenges as its distribution, stocks, politics, vaccination efforts, among others, emerge. A fair distribution and application strategy is imperative to manage the available resources, especially in developing countries.

We established an optimal control problem to design vaccination strategies where vaccination modulates dynamics susceptibility through an imperfect vaccine. We aimed to provide vaccination policies that minimize the lost life years due to disability or premature death by COVID-19, determined by cumulative deaths and cumulative incidence. Policies' acts in the minimization of infected people's prevalence and the number of deaths.

Our simulations suggest a better response with optimal vaccination policies than policies with a constant vaccination rate. For example, the optimal policy schedule in scenario SCN-1 increases the number of doses in the scheme's initial stage. This vaccination scheme improves the mitigation of the symptomatic prevalence, the incidence of deaths, and in consequence, the years of life lost quantified in DALYs.

Emerging press releases reported that Pfizer's, Russian Sputnik V, and Moderna's coronavirus vaccines reach efficacy over 90 % [41] [42] [43] . However, this information remains under development. Thus vaccine efficacy scenarios of 50 %, 70 %, and 90 % (SCN-2) illustrate the effect on optimal vaccination policies' schedule by pointing when to intensify the number of doses. According to the time horizon of one year and coverage of 50%, our numerical experiments suggest that 90 % vaccine efficacy reduces around three times the number of deaths regarding the dynamics without vaccination. Likewise, these vaccines reach a gain of eighteen times in the years of life lost compared to the without vaccination scenario.

Our numerical experiments also illustrate vaccine-induced immunity's relation between the reproductive vaccination number R V and vaccination policies (SCN-3) . Considering an outbreak with a reproductive 19 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.19.20235176 doi: medRxiv preprint number R 0 of 1.794 93, vaccine-induced immunity of 365 or 730 days implies a reduction of R 0 -dropping its value respectively to 1.139 130 and 0.867 56. Likewise, optimal policies linked to vaccine-induced immunities enhance symptomatic prevalence mitigation and the number of saved lives. Moreover, according to the initial number of deaths, the scenario without vaccination accumulates 503 deads compared to 211 and 206 deads of the underlying dynamics with vaccine-induced immunities.

Barbosa et al. recently report in [22] a simpler model about modeling of COVID-19 vaccination with a similar approach. Although they establish a less detailed model -they do not distinguish between symptomatic and asymptomatic infected individuals-its control problem return policies according to multiobjective policies. Their optimal policy is pragmatic but, in our opinion, not necessarily practical for large populations. Further, our model extends the result of [22] by a more detailed vaccine profile. Thus we can evaluate how vaccine efficacy, vaccine-induced, and natural immunity parameters impact the mitigation of an optimal vaccination schedule.

Perkins and España report in [20] a vaccination model with optimal control, but they approach to optimize the NPIs. The methodology presented here is similar, but aim very different. However, we want to stress the relevance of also including NPIs effects.

Since any vaccine's efficacy will be subject to uncertainty and immunization regarding COVID-19 remains under development, policymakers need better modeling tools to design fair vaccination programs. We faced this problem by simulation.

According to DALYs definition, segregation as age, comorbidities, and other risk groups is imperative to design more realistic vaccination policies. Moreover, it is well-known that various vaccines platforms and strategies are developing in parallel, and the most recent advance is with vaccines that require two doses. From [20] , we can deduce that NPIs, together with vaccination, would constitute a better description of COVID-19 control. We will direct our attention to extend this work according to the segregation and optimization of NPIs-vaccination controls.

All code is available at https://github.com/SaulDiazInfante/ NovelCovid19-ControlModelling/ tree/master/ UNISON-ITSON-VACCINATON-PRJ

The first and second authors named are leads. The second named author is the corresponding author. All other authors are listed in alphabetical order.

Manuel Adrian Acuña-Zegarra: is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint For the first stage, the following system is considered:

Here, we consider COVID-19 data from the first day of symptoms onset reported (February 19) until March 23, 2020. We also assume that θ = 0 because the first reported death was on March 18, and there were three reported deaths until March 23. The initial values of recovered and dead people are set to zero. Symptomatic class initial value was fixed in one individual, while E(0) and I A (0) were estimated. Thus, S(0) = N − (E(0) + I A (0) + 1), where N = 26446435 [48] . For the STAN implementation, we employ a negative-binomial model as the likelihood function with the mean parameter given by incidence solution per day. In addition to the above, we assign prior probability distributions to each parameter and the exposed and asymptomatic classes' initial conditions. Thus, we propose thatβ A andβ S follow a normal distribution with parameters µ = 1 and σ 2 = 0.13. Then, p follows a uniform distribution in (0, 0.25), and E(0) and I A (0) also follow a uniform distribution in (2, 20) and (2, 10) , respectively. When employing our STAN implementation, we run 5 chains with 100,500 iterations each, discard the first 500, and use 10,000 samples to generate estimates of parametersβ A ,β S and p. For the second stage, we took a complete month starting the day when mitigation measures were implemented, that is, from March 23 to April 23, 2020. Now, we consider parameter ξ to model the implementation 23 . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ;  of non-pharmaceutical measures. Thus, system in Equation (A.1) becomes:

whereN (t) = S(t) + E(t) + I S (t) + I A (t) + R(t) and N =N + D. At this stage, we consider that θ = 0.11. Here, our objective is to estimate the value of parameter ξ. To do this, we use the median posterior of all the estimated parameters from the first stage (see Table A .8).

Other parameter values are given in Table A Similar to the first stage, we consider a negative-binomial model as the likelihood function with the mean parameter given by incidence solution per day, while that we postulate a uniform distribution in (0.25, 0.75) as a prior probability distribution for the parameter ξ. For the second stage, we run 5 chains with 100,500 iterations each, discard the first 500, and use 10,000 samples to generate estimates of parameters ξ. Finally, it is important to mention that our results were implemented considering that the effective transmission contact rates (β • ) were equal to ξβ • . This last means that our scenarios consider the first reduction in the effective transmission contact rates by NIPs. Using values in Tables A.8 and A.9, we build confidence intervals for ξβ • . These results are shown in Table A is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted November 20, 2020. ;

World of Health Organization, Novel Coronavirus (2019-nCoV)-SITUATION REPORT -1

Isolation, quarantine, social distancing and community containment: pivotal role for old-style public health measures in the novel coronavirus (2019-ncov) outbreak

The effect of control measures on COVID-19 transmission in Italy: Comparison with Guangdong province in China

A review on Promising vaccine development progress for COVID-19 disease

COVID-19 Vaccine: A comprehensive status report

Unwavering Regulatory Safeguards for COVID-19 Vaccines

A vaccination model for transmission dynamics of influenza

Modelling the evolution trajectory of COVID-19 in Wuhan, China: experience and suggestions

Modeling and forecasting the COVID-19 pandemic in India

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

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

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

A mathematical model of COVID-19 using fractional derivative: outbreak in India with dynamics of transmission and control

The COVID-19 pandemic: model-based evaluation of nonpharmaceutical interventions and prognoses

Optimal control of vaccine distribution in a rabies metapopulation model

Vaccination models and optimal control strategies to dengue

Optimal control with multiple human papillomavirus vaccines

Optimal control of vaccination dynamics during an influenza epidemic

Controlling the Spread of COVID-19: Optimal Control Analysis

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

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

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

World Health Organization

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

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

The Burden of Disease due to COVID-19 in Korea Using Disability-Adjusted Life Years

Functional analysis, calculus of variations and optimal control

A survey of the maximum principles for optimal control problems with state constraints

Optimal control applied to biological models

Datos abiertos

Jump: A modeling language for mathematical optimization

Computing in operations research using julia

Mathworks, fmincon interior point algorithm

Practical Methods for Optimal Control Using Nonlinear Programming, Third Edition

Goddard problem in presence of a dynamic pressure limit

Team Commands, Bocop: an open source toolbox for optimal control

Bocop -A collection of examples

Source code for the manuscript Optimal Vaccine Policies for COVID-19

Collaborative data science

Pfizer says early analysis shows its Covid-19 vaccine is more than 90% effective

Russia says its sputnik v covid-19 vaccine is 92% effective

Moderna's coronavirus vaccine is 94.5% effective, according to company data

Immunological considerations for COVID-19 vaccine strategies

Coronavirus reinfections: three questions scientists are asking

Contemporary statistical inference for infectious disease models using stan

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

Parameter estimation Mathematical models for COVID-19 have shown that the parameters' values are not necessarily the same in each country. We use COVID-19 data from Mexico City plus Mexico state to follow the epidemic curve's initial growth in this work. Consequently, we estimate some parameter values of system in Equation

i) before and ii) after mitigation measures were implemented. For both stages, we use model in Equation (1) with no vaccination dynamics (λ V = 0 and V (0) = 0), and STAN R-package. This package is used for statistical inference by the Bayesian approach. For the code implementation of our system, we follow ideas of [46], and it is made freely available at [39]. For this section

The authors acknowledge support from grant DGAPA-PAPIIT IV100220 and the Laboratorio Nacional de Visualización Científica UNAM. MAAZ acknowledges support from PRODEP Programme (No. 511-6/2019-8291). DBC acknowledges support from PRODEP Programme (No. 511-6/2019-8022). We thank

The authors have no competing interests.Proof. Let Ω = {(S, E, I S , I A , R, D, V ) ∈ R 7 + : S + E + I S + I A + R + D + V = N }. First, note that for this model we have a closed population, which allows the solutions to be bounded superiorly by the total population.On the other hand, to show the positivity of the solutions with initial conditions 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 formOn the other hand, following ideas of [25] , the next generation matrix for this model, evaluated in the disease equilibrium point, is given byThen, the spectral radius of K is