Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 2, Number 4, December 2021, pp.273-289 https://doi.org/10.5206/mase/14233 AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION IN MEXICO: ANALYSIS AND FORECAST ÁNGEL G. C. PÉREZ AND DAVID A. OLUYORI Abstract. In this study, we propose and analyse an extended SEIARD model with vaccination. We compute the control reproduction number Rc of our model and study the stability of equilibria. We show that the set of disease-free equilibria is locally asymptotically stable when Rc < 1 and unstable when Rc > 1, and we provide a sufficient condition for its global stability. Furthermore, we perform numerical simulations using the reported data of COVID-19 infections and vaccination in Mexico to study the impact of different vaccination, transmission and efficacy rates on the dynamics of the disease. 1. Introduction The Coronavirus Disease 2019 (COVID-19) pandemic, caused by the Severe Acute Respiratory Syn- drome Coronavirus 2 (SARS-CoV-2) has caused a worldwide crisis due to its effects on society and global economy. Due to the absence of specific anti-COVID-19 medical treatments, most countries had been relying on non-pharmaceutical interventions, such as wearing of face masks, social/physical dis- tancing, partial/total lockdown, travel restrictions, and closure of schools and work centres, in order to curtail the spread of the disease before December 2020. However, these measures have been insufficient to mitigate the pandemic globally as medical facilities were overstretched and death toll heightened. Vaccination has been an effective strategy in combatting the spread of infectious diseases, e.g., pertussis, measles, and influenza. Historically, the eradication of smallpox has been considered as the most remarkable success of vaccination ever recorded [40]. So far, the development and testing of vaccines against SARS-CoV-2 has occurred at an unprecedented speed and, in the last months, several vaccines have been approved for use in many countries, and their deployment is already underway. In a pandemic situation such as this, current preventive vaccines consisting of inactivated viruses do not protect all vaccine recipients equally as the protection conferred by the vaccine is dependent on the immune status of the recipient [3]. Over the past few decades, a large number of simple compartmental models with vaccination have been proposed in the literature to assess the effectiveness of vaccines in combatting the infectious diseases [19, 28, 4, 1, 16, 21, 9, 8, 37, 30]. With the recent development of anti-COVID vaccines, several models have been proposed to provide insight into the effect that inoculation of a certain portion of the population will have on the dynamics of the COVID-19 pandemic. The authors in [10] studied an SEIHRDV model and an SMEIHRDV model (the latter including a semi-susceptible class) and fitted the parameters using data from several countries to evaluate the effect of social distancing and vaccination on controlling COVID-19. Another model was studied in [25] to compare the outcomes of single-dose and two-dose anti-COVID vaccination regimes. The global stability of a two-strain COVID-19 model Received by the editors 13 September 2021; accepted 8 November 2021; published online 10 November 2021. 2010 Mathematics Subject Classification. 92D30, 34D20, 37M05. Key words and phrases. Coronavirus, vaccine efficacy, asymptomatic infection, stability. 273 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/14233 274 Á. G. C. PÉREZ AND D. A. OLUYORI with vaccination against one strain was studied in [32]. The stability analysis of an SEIR model with prophylactic and therapeutic vaccines was performed in [38]. In [20], the effect of immunity, vaccination, and reinfection with changing parameters was analysed using an SEVIS model, while in [13] an SIQRD model was used to simulate several scenarios of vaccine delivery in Indonesia. Since it is known that individuals infected with SARS-CoV-2 can transmit the virus to other peo- ple without presenting symptoms of the disease, some authors have proposed models that distinguish between symptomatic and asymptomatic infections. For instance, an age- and region-structured model was proposed in [22] to simulate the rollout of a two-dose vaccination programme in the UK using the Pfizer–BioNTech and Oxford–AstraZeneca vaccines. An agent-based transmission model was used in [34] to project the impact of a two-dose vaccination campaign with the Pfizer–BioNTech and Moderna vaccines in Ontario. On the other hand, a sensitivity analysis and uncertainty quantification of an SEIsI- aQR model with vaccination strategy was conducted in [23]. A three-patch metapopulation epidemic model structured by risk was used to investigate several vaccination scenarios in Mexico City in [27]. The threshold dynamics of a COVID-19 model combining vaccination and treatment was established in [12], while [11] considered an SE(Is)(Ih)AR model with vaccination and antiviral controls. Some other COVID-19 models with vaccination can be found in [14, 31, 5, 2, 26]. Of the works mentioned above, only [11, 12, 32, 38] studied the global stability of the equilibrium states; moreover, none of these studies considered that people who become infected after being vac- cinated may have reduced infectivity and a lower chance of showing severe symptoms. Hence, in the present work, we aim to study a model that takes into account all of these factors and perform a local and global stability analysis of the disease-free equilibria, as well as present an application of this model to the COVID-19 epidemic in Mexico. In mid-November 2020, Mexico passed the mark of 1,000,000 confirmed cases and 100,000 deaths due to COVID-19. On 11 December 2020, the Mexican government’s medical safety commission approved the emergency use of the Pfizer–BioNTech coronavirus vaccine, with the first 250,000 doses intended for health workers. The inoculation of frontline health personnel started that year on 24 December. Four other COVID-19 vaccines were approved between January and February 2021, and vaccination of people over 60-years old with the AstraZeneca vaccine began in February. The first vaccines to be applied in Mexico followed a two-dose regime, until the single-dose CanSino vaccine began to be applied to a portion of the population in April 2021. The Johnson & Johnson vaccine, also single-dose, was deployed in June 2021 to inoculate the population over age 18 in the municipalities of the northern border. The government expected to have vaccinated all adult population with at least one dose by October 2021, after which, the application of vaccines would continue for people aged 12 to 17 with comorbidities. In this paper, we propose a differential equation model to simulate the application of a two-dose vaccine against COVID-19, considering the possibility of vaccine leakiness and asymptomatic infections. The motivation of this study is derived from the work of the authors in [7, 6], who considered an SEIARD mathematical model to investigate the outbreak of COVID-19 in Mexico. Therefore, in the present work, we incorporate the vaccination component to the model in [6] to derive an extended SEIARD model to examine the effectiveness of the COVID-19 jabs which are currently being deployed to many countries to help combat the raging pandemic situation. The rest of this paper is organized as follows. In Section 2, we present the equations and assumptions of the extended SEIARD model with vaccination. In Section 3, we perform a theoretical analysis of the model, compute its control reproduction number and study the stability of the disease-free equilibria. In Section 4, we carry out numerical simulations using reported data on COVID-19 infections and vaccination in Mexico. Lastly, we provide some discussions and concluding remarks in Section 5. AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 275 Figure 1. Flow diagram of the model with vaccination. 2. Model formulation To derive the mathematical model, we divide the unvaccinated population into susceptible (S), exposed (E), symptomatic infectious (I), asymptomatic infectious (A), and recovered (R). The number of individuals in each subpopulation at time t is denoted by S(t), E(t), etc. Since the time scale we consider in our analysis is considerably shorter than the mean lifespan of individuals in the population, our model does not consider the natural death and birth rates. Susceptible individuals become exposed by contact with symptomatic infectious individuals at a rate β1 and by contact with asymptomatic infectious individuals at a rate β2. The exposed individuals become infectious at a rate w: a proportion p1 of them will show symptoms, while the rest remains asymptomatic. We assume that the symptomatic class has a disease-induced death rate, denoted by δ1. Both symptomatic and asymptomatic infectious people recover at a rate γ. We also assume that the susceptible population S is vaccinated at a rate v ≥ 0 (the number of first doses administered per day at time t is given by vS(t)). Individuals who have received only the first dose of the vaccine are included in the class V1, and they move to the class V2 upon receiving the second dose, which occurs at a rate θ. Due to vaccine leakiness, vaccination of an individual does not completely remove the risk of infection. Hence, we also assume that the vaccinated population can become exposed (EV ), symptomatic infectious (IV ) and asymptomatic infectious (AV ). Individuals in the class V1 (respectively, V2) move to the class EV due to contact with symptomatic infectious people at a rate (1 − η1)β1 (respectively, (1 − η2)β1) and by contact with the asymptomatic infectious at a rate (1 − η1)β2 (respectively, (1 − η2)β2), where η1 is the efficacy of the vaccine after one dose (η2 is the efficacy of the vaccine after two doses). The infectivity of individuals in the IV and AV classes is reduced by a factor 1 −q with respect to that of individuals in the I and A classes. The population in the class EV becomes infectious at a rate w; we assume that the proportion of people from this class who become symptomatic infectious is p2, which may be different from that of unvaccinated people due to the effect of the vaccine in reducing the severity of the infection. Likewise, the disease-induced death rate δ2 is lower for the vaccinated population. Individuals in the IV and AV classes also move to the R class upon recovery at a rate γ. We will denote by N(t) the total population at time t, which is given by N(t) = S(t) + E(t) + I(t) + A(t) + V1(t) + V2(t) + EV (t) + IV (t) + AV (t) + R(t). 276 Á. G. C. PÉREZ AND D. A. OLUYORI Hence, our model is described by the following system of differential equations: Ṡ = − S N β1(I + qIV ) − S N β2(A + qAV ) −vS, Ė = S N β1(I + qIV ) + S N β2(A + qAV ) −wE, İ = p1wE − (δ1 + γ)I, Ȧ = (1 −p1)wE −γA, V̇1 = vS − (1 −η1) V1 N β1(I + qIV ) − (1 −η1) V1 N β2(A + qAV ) −θV1, V̇2 = θV1 − (1 −η2) V2 N β1(I + qIV ) − (1 −η2) V2 N β2(A + qAV ), ĖV = (1 −η1) V1 N β1(I + qIV ) + (1 −η1) V1 N β2(A + qAV ) + (1 −η2) V2 N β1(I + qIV ) + (1 −η2) V2 N β2(A + qAV ) −wEV , İV = p2wEV − (δ2 + γ)IV ȦV = (1 −p2)wEV −γAV , Ṙ = γ(I + A + IV + AV ). (2.1) We define an additional variable D(t) that denotes the number of people deceased due to COVID-19, which is governed by the equation Ḋ = δ1I + δ2IV . (2.2) We assume that β1, β2, θ, w and γ are positive, v,δ1,δ2 ≥ 0 and q,η1,η2,p1,p2 ∈ [0, 1]. The interpretation of parameters is summarized in Table 1. The flow diagram of the model can be seen in Figure 1. 3. Theoretical analysis In this section, we will derive some theoretical results for model (2.1). Our aim is to investigate the behaviour of the epidemic in the next few years after the vaccination campaign is over. Thus, we will mainly focus on the case when the vaccination rate v is zero, but the vaccinated subpopulations (V1, V2, EV , IV , AV ) may have positive initial values. First, we will determine the disease-free equilibria of the model, i.e., the equilibria with E = I = A = EV = IV = AV = 0. In the case when v is positive, the disease-free equilibria (DFE) of system (2.1) are given by the set S = { (S,E,I,A,V1,V2,EV ,IV ,AV ,R) = (0, 0, 0, 0, 0,V ∗ 2 , 0, 0, 0,R ∗) ∈ R10 : V ∗2 > 0, R ∗ > 0 } . The equilibria in S correspond to the case when the whole population has been fully vaccinated or recovered from COVID-19, and there are no susceptible individuals. When v = 0, model (2.1) has a continuum of disease-free equilibria, given by the set S0 = { (S,E,I,A,V1,V2,EV ,IV ,AV ,R) = (S ∗, 0, 0, 0, 0,V ∗2 , 0, 0, 0,R ∗) ∈ R10 : S∗ ≥ 0, V ∗2 > 0, R ∗ > 0 } . Each equilibrium P0 ∈ S0 represents a scenario when the anti-COVID vaccination programme has ended, and a certain number of individuals V ∗2 has been vaccinated to achieve herd immunity in the population. We will compute the control reproduction number Rc of the model based on this expression for the DFE. AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 277 Using the notation in [33], we determine the matrix of new infections F and the transition matrix V , considering only the infected compartments (E, I, A, EV , IV and AV ). First, we define F =   S N β1(I + qIV ) + S N β2(A + qAV ) 0 0 (1 −η1)V1N [ β1(I + qIV ) + β2(A + qAV ) ] + (1 −η2)V2N [ β1(I + qIV ) + β2(A + qAV ) ] 0 0   and V =   wE −p1wE + (δ1 + γ)I −(1 −p1)wE + γA wEV −p2wEV + (δ2 + γ)IV −(1 −p2)wEV + γAV   . Then, the derivative of F at a disease-free equilibrium P0 ∈S0 is F =   0 S ∗ N∗ β1 S∗ N∗ β2 0 S∗ N∗ qβ1 S∗ N∗ qβ2 0 0 0 0 0 0 0 0 0 0 0 0 0 (1 −η2) V ∗2 N∗ β1 (1 −η2) V ∗2 N∗ β2 0 (1 −η2) V ∗2 N∗ qβ1 (1 −η2) V ∗2 N∗ qβ2 0 0 0 0 0 0 0 0 0 0 0 0   where N∗ = S∗ +V ∗2 +R ∗ denotes the total population at the equilibrium. The derivative of V evaluated at P0 is V =   w 0 0 0 0 0 −p1w δ1 + γ 0 0 0 0 −(1 −p1)w 0 γ 0 0 0 0 0 0 w 0 0 0 0 0 −p2w δ2 + γ 0 0 0 0 −(1 −p2)w 0 γ   . It follows that FV −1 =   A11 S∗β1 N∗(δ1+γ) S∗β2 N∗γ A14 S∗qβ1 N∗(δ2+γ) S∗qβ2 N∗γ 0 0 0 0 0 0 0 0 0 0 0 0 A41 (1−η2)V ∗2 β1 N∗(δ1+γ) (1−η2)V ∗2 β2 N∗γ A44 (1−η2)V ∗2 qβ1 N∗(δ2+γ) (1−η2)V ∗2 qβ2 N∗γ 0 0 0 0 0 0 0 0 0 0 0 0   , 278 Á. G. C. PÉREZ AND D. A. OLUYORI where A11 = S∗ N∗ [ β1p1 δ1 + γ + β2(1 −p1) γ ] , A14 = S∗ N∗ q [ β1p2 δ2 + γ + β2(1 −p2) γ ] , A41 = V ∗2 N∗ (1 −η2) [ β1p1 δ1 + γ + β2(1 −p1) γ ] , A44 = V ∗2 N∗ q(1 −η2) [ β1p2 δ2 + γ + β2(1 −p2) γ ] . The control reproduction number Rc of model (2.1) is given by Rc = ρ ( FV −1 ) , where ρ denotes the spectral radius. Hence, Rc = S∗ N∗ [ β1p1 δ1 + γ + β2(1 −p1) γ ] + V ∗2 N∗ q(1 −η2) [ β1p2 δ2 + γ + β2(1 −p2) γ ] . (3.1) The quantity Rc measures the average number of new COVID-19 cases generated by a typical infectious individual introduced into a population where a fraction V ∗2 /N ∗ has been fully vaccinated using a vaccine with efficacy η2. According to [33, Theorem 2], we can obtain the following result about the control reproduction number. Theorem 3.1. A disease-free equilibrium P0 = (S ∗, 0, 0, 0, 0,V ∗2 , 0, 0, 0,R ∗) of system (2.1) with v = 0 is locally asymptotically stable if Rc < 1, and it is unstable if Rc > 1. The epidemiological interpretation of Theorem 3.1 is that, if Rc < 1, then a small perturbation from a disease-free equilibrium will not generate an epidemic outbreak. On the other hand, if Rc > 1, the epidemic curve will initially show an exponential growth, then reach a peak and start to decrease until becoming extinct. The following theorem gives a sufficient condition for the global stability of the disease-free equilibria. Theorem 3.2. Suppose that β1p1 δ1 + γ + β2(1 −p1) γ < 1 and q(1 −η2) [ β1p2 δ2 + γ + β2(1 −p2) γ ] < 1. (3.2) Then, the disease-free equilibrium P0 = (S ∗, 0, 0, 0, 0,V ∗2 , 0, 0, 0,R ∗) of system (2.1) with v = 0 is globally asymptotically stable. Proof. Consider the following Lyapunov function: L = g1E + g2I + g3A + g4EV + g5IV + g6AV + g7V1, where g1 = γ(δ1 + γ), g2 = γβ1, g3 = (δ1 + γ)β2, g4 = γ(δ1 + γ) 1 −η2 , g5 = qγβ1(δ1 + γ) δ2 + γ , g6 = qβ2(δ1 + γ), g7 = γ(δ1 + γ) 1 −η2 . AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 279 The time derivative of L evaluated at the solutions of system (2.1) with v = 0 is given by L̇ = g1 [ S N β1(I + qIV ) + S N β2(A + qAV ) −wE ] + g2 [ p1wE − (δ1 + γ)I ] + g3 [ (1 −p1)wE −γA ] + g4 [ (1 −η1) V1 N β1(I + qIV ) + (1 −η1) V1 N β2(A + qAV ) + (1 −η2) V2 N β1(I + qIV ) + (1 −η2) V2 N β2(A + qAV ) −wEV ] + g5 [ p2wEV − (δ2 + γ)IV ] + g6 [ (1 −p2)wEV −γAV ] + g7 [ −(1 −η1) V1 N β1(I + qIV ) − (1 −η1) V1 N β2(A + qAV ) −θV1 ] . After cancelling terms and simplifying, we obtain L̇ = γβ1(δ1 + γ)(I + qIV ) ( S N + V2 N − 1 ) + γβ2(δ1 + γ)(A + qAV ) ( S N + V2 N − 1 ) + wγ(δ1 + γ) [ β1p1 δ1 + γ + β2(1 −p1) γ − 1 ] E + wγ(δ1 + γ) 1 −η2 [ qβ1p2(1 −η2) δ2 + γ + qβ2(1 −p2)(1 −η2) γ − 1 ] EV − γθ(δ1 + γ) 1 −η2 V1. Since S(t) + V2(t) ≤ N(t) for all t, we have SN + V2 N ≤ 1. Combining this with the hypothesis (3.2), we can see that L̇ ≤ 0, and L̇ = 0 if and only if E(t) = 0 and EV (t) = 0. Substituting E(t) = 0 and EV (t) = 0 in system (2.1) with v = 0 shows that (S,E,I,A,V1,V2,EV ,IV ,AV ,R) → (S∗, 0, 0, 0, 0,V ∗2 , 0, 0, 0,R ∗) as t → ∞. Hence, the largest positively invariant set where L̇ = 0 is the continuum of disease-free equilibria. Therefore, by LaSalle’s invariance principle, we conclude that the DFE is globally asymptotically stable. � 3.1. Impact of vaccination coverage on the control reproduction number. Next, we will study how the control reproduction number Rc is affected by some of the model parameters. By equation (3.1), we know that Rc does not only depend on the parameters of system (2.1), but also on the final proportions of unvaccinated susceptible people (S∗/N∗) and fully vaccinated people (V ∗2 /N ∗) at the time when vaccines are no longer being deployed to the population. We recall that a disease-free equilibrium takes the form P0 = (S ∗, 0, 0, 0, 0,V ∗2 , 0, 0, 0,R ∗), where the total population is N∗ = S∗ + V ∗2 + R ∗. If we define x = V ∗2 N∗ as the proportion of fully vaccinated people and y = R∗ N∗ as the proportion of people recovered from COVID-19, we can rewrite the expression for the control reproduction number as Rc(x,y) = [ β1p1 δ1 + γ + β2(1 −p1) γ ] (1 −x−y) + q(1 −η2) [ β1p2 δ2 + γ + β2(1 −p2) γ ] x. Figure 2 depicts the value of Rc as function of the proportions x and y, using several values for the transmission rates and efficacy after the second vaccine dose. Other parameter values were taken as in Table 1. We can see that an increase in either x or y contributes to reducing the reproduction number, and therefore, is helpful towards achieving herd immunity. 280 Á. G. C. PÉREZ AND D. A. OLUYORI Figure 2. Value of the control reproduction number as function of the proportion of fully vaccinated individuals (horizontal axis) and recovered individuals (vertical axis). Herd immunity occurs when a large portion of the population has become immune to the disease due to vaccination or natural recovery, which makes spread of the disease difficult. Thus, the minimal level of vaccination coverage that is required to achieve herd immunity (that is, making Rc < 1) will also depend on the percentage of the population that has been infected and then successfully recovered. Comparing the different panels of Figure 2, we can see that increasing the vaccine efficacy η2 reduces the vaccination coverage needed to make Rc < 1 for a fixed proportion of recovered people. However, this reduction is small compared to the effect gained by decreasing the transmission rate. For example, when η2 = 0.65 and the recovered population is close to zero, it is necessary to vaccinate 46% of population to obtain Rc = 1 in the case of 120% transmission rate, 33% in the case of baseline transmission rate, and only 11% of population in the case of 80% transmission rate (bottom row of Figure 2). 4. Numerical simulations In this section, we perform some numerical simulations for model (2.1) to provide estimates for the evolution of the COVID-19 outbreak in Mexico. 4.1. Data fitting and estimation of parameters. We used cumulative data provided by the Johns Hopkins University repository [18] to fit the parameters of model (2.1) in the absence of vaccination. We considered the data for reported COVID-19 infections, deaths and recovered cases during the period from 12 November 2020 to 24 December 2020, which is before the vaccination programme in Mexico began. For this part, we considered system (2.1) with v = 0 and the vaccinated subpopulations V1, V2, EV , AV and IV equal to zero. For the other variables, we used the initial conditions S(0) = 1.1938 × 108, E(0) = 1.58582 × 105, I(0) = 1.58582 × 105, A(0) = 1.1629 × 106, R(0) = 6.134975 × 106, D(0) = 97056. AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 281 Table 1. Baseline values for the parameters used in the simulations. Parameter Description Value Source β1 Transmission rate by contact with I and IV classes 0.2 day −1 Fitted to data β2 Transmission rate by contact with A and AV classes 0.0330 day −1 Fitted to data q Relative infectivity of individuals in IV and AV classes 0.52 [17] v Vaccination rate Variable Assumed θ Application rate of second doses 1/70 day−1 Assumed η1 Vaccine efficacy rate after one dose 0.463 [35] η2 Vaccine efficacy rate after two doses 0.557 [35] w Transfer rate from exposed to infectious 0.25 day−1 [29] p1 Portion of people in E class that develop symptoms 0.12 [6] p2 Portion of people in EV class that develop symptoms 0.089 Estimated δ1 Disease-induced death rate of I class 3.2135 × 10−3 day−1 Fitted to data δ2 Disease-induced death rate of IV class 0 Assumed γ Recovery rate of infectious individuals 3.6987 × 10−2 day−1 Fitted to data The values for I(0), R(0) and D(0) were chosen based on the reported data for 12 November 2020. The value of A(0) was chosen so the symptomatic infections represent 0.12 times the total infections, i.e., I(0) = 0.12 [A(0) + I(0)]. The value of E(0) was assumed equal to I(0), and S(0) was estimated as S(0) = N −E(0) − I(0) −A(0) −R(0) −D(0), where N is the population of Mexico. We regarded as fixed parameters w = 0.25, which corresponds to a latent period of 4 days [29], and a proportion p1 = 0.12 of symptomatic infections [6]. The set of differential equations was solved using Matlab 2016b with the ode45 solver. The values for β1, β2, δ1 and γ were estimated by minimizing the Sum of Squared Errors (SSE), which is calculated as follows. For a given vector of parameters x, we compute numerically the I(t) and D(t) components of the solutions for our model, as well as the estimated number of people recovered from symptomatic infections RI(t) and the cumulative number of symptomatic infected cases C(t), defined by ṘI(t) = γI(t), C(t) = I(t) + RI(t) + D(t). Then the Sum of Squared Errors is given by SSE(x) = n∑ i=1 [ k1 (C(ti) −C exp i ) 2 + k2 (D(ti) −D exp i ) 2 + k3 (RI(ti) −R exp i ) 2 ] , where C exp i , D exp i and R exp i denote the experimental data for cumulative infections, deaths and re- coveries, respectively, reported for day ti (i = 1, . . . ,n), while k1, k2 and k3 are coefficients used to compensate the order of magnitude for the data. In our simulations, we used k1 = 20, k2 = 10 and k3 = 1. The global minimum of the SSE function was obtained by applying three consecutive searches: a gradient-based method, a gradient-free algorithm and, again, a gradient-based method. The best fit values for β1, β2, δ1 and γ are shown in Table 1. Figure 3 depicts a comparison between the model solutions and the observed cumulative COVID-19 data before the vaccination period. 4.2. Simulations for the model with vaccination. We will now simulate the solutions to model (2.1) to assess the impact of the vaccination campaign that started in Mexico in December 2020 to combat the COVID-19 pandemic. As of August 2021, seven COVID-19 vaccines have received Emergency Use Authorization for their deployment in Mexico: • BNT162b2 (Pfizer–BioNTech), 282 Á. G. C. PÉREZ AND D. A. OLUYORI Nov 12 Nov 19 Nov 26 Dec 03 Dec 10 Dec 17 Dec 24 2020 0.8 1 1.2 1.4 P e o p le 10 6 Cumulative number of symptomatic cases Model solutions Reported data Nov 12 Nov 19 Nov 26 Dec 03 Dec 10 Dec 17 Dec 24 2020 0.8 1 1.2 1.4 P e o p le 10 5 Death toll Model solutions Reported data Nov 12 Nov 19 Nov 26 Dec 03 Dec 10 Dec 17 Dec 24 2020 6 8 10 12 P e o p le 10 5 Recovered from symptomatic infection Model solutions Reported data Figure 3. Reported cumulative number of symptomatic cases, COVID-19 deaths and recovered cases in Mexico for the pre-vaccination period, and simulations using model (2.1) with the parameters in Table 1 and v = 0. • AZD1222 (Oxford–AstraZeneca), • Sputnik V (Gamaleya Institute), • CoronaVac (Sinovac), • BBV152 (Covaxin), • Ad5-nCoV (CanSino), and • Ad26.COV2-S (Johnson & Johnson). The first five of these vaccines require two doses, while CanSino and Johnson & Johnson are single-dose vaccines [15]. Efficacy estimates for each vaccine based on data from clinical trials are subject to change with the emergence of new analyses. An interim analysis for the Oxford–AstraZeneca vaccine [35] estimated an efficacy against infection (symptomatic or asymptomatic) of 46.3% (31.8%–57.8%), considering people who had a nucleic acid amplification test (NAAT)-positive swab more than 21 days after a single dose, and 55.7% (41.1%–66.7%) for people who tested positive more than 14 days after a second dose of the vaccine. However, a more recent study [36] estimated an efficacy of 63.9% (46.0%–76.9%) after one dose and 59.9% (35.8%–75.0%) after two standard doses given 12 or more weeks apart. Due to longer dose intervals being associated with greater efficacy against symptomatic infection, the WHO has recommended to administer the Oxford–AstraZeneca vaccine with an interval of 8 to 12 weeks between first and second doses [39]. Based on the above, we will assume in our simulations an average length of 1/θ = 70 days for the inter-dose period, and we will use η1 = 0.463 and η2 = 0.557 as baseline values for the efficacy parameters. Furthermore, we assume a reduction of 48% in the infectivity of individuals becoming infected after being vaccinated (i.e., q = 0.52), following the estimations in [17]. AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 283 Table 2. Estimated values for the vaccination rate v. Date Value (day−1) 24 Dec 2020 – 11 Jan 2021 4.0 × 10−5 12 Jan 2021 – 15 Jan 2021 7.9 × 10−4 16 Jan 2021 – 14 Feb 2021 6.0 × 10−5 15 Feb 2021 – 7 Mar 2021 7.3 × 10−4 8 Mar 2021 – 14 Mar 2021 0.0021 15 Mar 2021 – 31 Mar 2021 0.0017 1 Apr 2021 – 15 Apr 2021 0.0035 16 Apr 2021 – 15 May 2021 0.0017 16 May 2021 – 27 Jul 2021 0.0060 Jan 2021 Mar 2021 May 2021 Jul 2021 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 p e o p le 10 7 0% 5% 10% 15% 20% 25% 30% 35% % o f M e x ic a n p o p u la ti o n Vaccinated with at least one dose Model solutions Reported data Jan 2021 Mar 2021 May 2021 Jul 2021 0 0.5 1 1.5 2 2.5 p e o p le 10 7 0% 2% 4% 6% 8% 10% 12% 14% 16% 18% % o f M e x ic a n p o p u la ti o n Fully vaccinated Figure 4. COVID-19 vaccination coverage in Mexico from 24 December 2020 to 27 July 2021. X represents real data, continuous lines represent model simulations using the vaccination rate in Table 2. For computing the proportion of infectious vaccinated individuals that show symptoms of the disease (p2), we follow [35], who reported 37 cases of symptomatic COVID-19 disease out of a total of 68 NAAT-positive swabs in the group of people vaccinated with AZD1222, and 112 symptomatic cases out of 153 NAAT-positive cases in the control group. This yields a reduction from 0.732 to 0.544 in the symptomatic proportion after vaccination. Since we have chosen p1 = 0.12, we will take p2 = 0.089 so that p1 : p2 = 0.732 : 0.544. Furthermore, we assume that the death rate δ2 of infectious vaccinated people is zero since it is widely accepted that current anti-COVID vaccines provide full protection against severe infections. We used the daily data on COVID-19 vaccinations in Mexico obtained from [24] to estimate the value of the vaccination rate v over nine different date ranges, as shown in Table 2. We plot in Figure 4 a comparison of the reported number of vaccinated people and the simulations obtained with model (2.1) for the period 24 December 2020 – 27 July 2021. In these graphs, we considered the total population of Mexico as 127 090 000 people. In order to obtain long-term projections for the vaccination coverage in Mexico, we simulated two different scenarios. First, we assumed that the vaccination rate is kept constant at its baseline value on 27 July 2021 (0.6% of susceptible population per day, which equals roughly 294 000 first doses per day). Second, we assumed that the vaccination rate increases to twice its baseline value starting on September 284 Á. G. C. PÉREZ AND D. A. OLUYORI Jan 2021 Jul 2021 Jan 2022 Jul 2022 Jan 2023 0 10 20 30 40 50 60 70 80 90 m il li o n s o f p e o p le 0% 10% 20% 30% 40% 50% 60% 70% % o f M e x ic a n p o p u la ti o n Vaccinated with at least one dose Reported data Baseline vaccination rate 200% vaccination rate Jan 2021 Jul 2021 Jan 2022 Jul 2022 Jan 2023 0 10 20 30 40 50 60 70 80 90 m il li o n s o f p e o p le 0% 10% 20% 30% 40% 50% 60% 70% % o f M e x ic a n p o p u la ti o n Fully vaccinated Figure 5. Long-term projections of COVID-19 vaccination coverage in Mexico. X represents real data, continuous lines represent simulations using the baseline vacci- nation rate, and dashed lines represent simulations using 200% vaccination rate from September 2021 onwards. 2021. Figure 5 shows that, if vaccines continue to be delivered at their baseline rate, the vaccination coverage with at least one dose will have reached 70% of Mexican population by January 2023. On the other hand, if the vaccination rate is doubled, the same coverage with at least one dose could be achieved by May 2022, and 70% of the Mexican population could be fully vaccinated by September 2022. 4.2.1. Assessing the effect of vaccination and different transmission rates. We will next compute the solutions of model (2.1) to simulate the evolution of the pandemic in Mexico as the vaccination campaign takes place. We consider the initial date for simulations as 24 December 2020. Based on the results obtained in Subsection 4.1, we use the initial conditions S(0) = 1.1622 × 108, E(0) = 3.4415 × 105, I(0) = 2.1247 × 105, A(0) = 1.6521 × 106, R(0) = 8.5421 × 106, D(0) = 1.2128 × 105, and V1(0) = V2(0) = EV (0) = IV (0) = AV (0) = 0. In these subsection, we will consider different values for the transmission rates β1 and β2 to account for the possibility that the number of infectious contacts between people may increase or decrease due to resumption of economical activities, compliance with social/physical distancing, wearing of face masks, etc. Hence, we consider three cases: when β1 and β2 are kept with the values in Table 1, when both of them decrease to an 80% of these values, and when they increase to a 120%. The values for other parameters are fixed as in Tables 1 and 2. Figure 6 depicts the time evolution of the number of infectious COVID-19 cases with symptoms (I(t)+ IV (t)) and the death toll (D(t)) for each of the above cases. In each graph, we have plotted the solutions assuming the baseline vaccination rate and the 200% vaccination rate, as well as a counterfactual case with no vaccination. Figure 6(a) shows that, in the case of low transmission rate, the number of active cases would start to decrease in the early months of 2021, and the epidemic would be almost extinguished by January 2022. In the cases with higher transmission rate (Figures 6(b) and (c)), the epidemic curve would reach its peak around May 2021, and the number of active symptomatic cases would be less than 1000 by May 2022. Figures 6(d)–(f) show that the cumulative number of deaths would be around 250 000 for low transmission, 390 000 for baseline transmission, and 580 000 for high transmission rates. AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 285 Jan 2021 Jan 2022 Jan 2023 0 2 4 6 8 10 I( t) + I V (t ) 10 5 80% transmission rate No vaccination Baseline vaccination rate 200% vaccination rate Jan 2021 Jan 2022 Jan 2023 2 4 6 8 D (t ) 10 5 Jan 2021 Jan 2022 Jan 2023 0 2 4 6 8 10 I( t) + I V (t ) 10 5 Baseline transmission Jan 2021 Jan 2022 Jan 2023 2 4 6 8 D (t ) 10 5 Jan 2021 Jan 2022 Jan 2023 0 2 4 6 8 10 I( t) + I V (t ) 10 5 120% transmission rate Jan 2021 Jan 2022 Jan 2023 2 4 6 8 D (t ) 10 5 (f) (c)(b)(a) (e)(d) Figure 6. Simulations of model (2.1) using different values for the transmission rates. (a) and (d): 80% transmission rate. (b) and (e): baseline transmission rate. (c) and (f): 120% transmission rate. Top row: number of active symptomatic infectious cases. Bottom row: Cumulative number of deaths. We can also see that an increase in the vaccination rate to double its baseline value does not result in a considerable change in the number of infections or deaths, although the vaccination scenarios result in around 250 000 less deaths compared with the case with no vaccination. On the other hand, comparing Figures 6(d) and (e) shows that more than 130 000 deaths can be avoided by reducing the transmission rate to 80%, while a 20% increase in the transmission rate would result in almost 200 000 additional deaths (Figure 6(f)). This suggests that decreasing the number of infectious contacts by complying with preventive measures is more effective than simply accelerating the deployment of vaccines. 4.2.2. Assessing the effect of different vaccine efficacy rates. Given that there is still uncertainty re- garding the efficacy of anti-COVID vaccines against infection, including asymptomatic cases, we will also simulate the solutions of model (2.1) using different values for the parameters η1 and η2. Figure 7 shows the number of active infectious cases with symptoms (I(t) + IV (t)) and without symptoms (A(t) + AV (t)), as well as the death toll (D(t)), using different efficacy rates: in addition to the baseline case (η1 = 0.463, η2 = 0.557), we include a case with lower efficacy (η1 = 0.4, η2 = 0.45) and a case with higher efficacy (η1 = 0.6, η2 = 0.65). Here, we have plotted all solutions using the baseline vaccination rate. The simulations show that lower efficacy results in an additional 7 429 symptomatic cases and 79 624 asymptomatic cases at the peak of the infection curve, compared with the case with higher efficacy. However, this does not significantly affect the time when the peak occurs. Moreover, lower efficacy also results in 4 953 additional deaths. 5. Conclusion In this work, we studied a model for COVID-19 with vaccination. Our work was based on the SEIARD model proposed in [6], which included an exposed (latent) compartment and different transmission rates for the symptomatic and asymptomatic infectious individuals; we extended this model by incorporating 286 Á. G. C. PÉREZ AND D. A. OLUYORI Jan 2021 Jan 2022 Jan 2023 0 1 2 3 4 5 I( t) + I V (t ) 10 5 Symptomatic infectious 1 =0.4, 2 =0.45 1 =0.463, 2 =0.557 1 =0.6, 2 =0.65 Jan 2021 Jan 2022 Jan 2023 0 0.5 1 1.5 2 2.5 3 3.5 A (t )+ A V (t ) 10 6 Asymptomatic infectious Jan 2021 Jan 2022 Jan 2023 1 1.5 2 2.5 3 3.5 4 D (t ) 10 5 Death toll Figure 7. Simulations of model (2.1) using different values for the vaccine efficacy rates. Left panel: number of active symptomatic infectious cases. Central panel: number of active asymptomatic infectious cases. Right panel: cumulative death toll. vaccinated compartments and considering a two-dose vaccination regime. Although several COVID-19 models with vaccination have been proposed in the literature, most studies have focused only on carrying out numerical simulations, while our work shed some light on the theoretical properties of model (2.1). When compared to other models analysed in the literature, we can see that the novelty of our model lies in the inclusion of a parameter q representing the reduction in infectivity due to vaccination, as well as the different probabilities (p1 and p2) of developing COVID-19 symptoms depending on whether the infected person has been vaccinated. We showed that our model has multiple disease-free equilibria and computed the control reproduction number Rc using the next-generation matrix method. We established that the set of disease-free equilibria is locally asymptotically stable when Rc < 1 and unstable when Rc > 1. Furthermore, we determined a condition that guarantees the global asymptotic stability of the DFE. We performed a numerical simulation on our model using repository data on the outbreak of COVID- 19 in Mexico and the daily data on COVID-19 vaccinations to estimate the value of the vaccination rate over nine different date ranges. We used the efficacy estimates based on data from clinical trials of the Oxford–AstraZeneca vaccine, which is the one that is being more widely distributed in Mexico at the time of this writing. We remark that, in this article, we considered vaccine efficacy in the sense of protection against COVID-19 infection (symptomatic or asymptomatic), while other works consider efficacy as protection against symptomatic infection only. We simulated two different scenarios to obtain projections for the vaccination coverage in the next few years. First, we assumed that the vaccination rate is kept constant by vaccinating the same proportion of susceptible individuals per day, and secondly, we assumed that the vaccination rate increases to twice its baseline value from September 2021 onwards. Our study showed that if vaccines continue to be delivered at their baseline rate, by January 2023 the first dose will be applied to 90 million people, which represents roughly the total Mexican population over age 18. On the other hand, if the vaccination rate is doubled, the total adult population could be partially vaccinated by May 2022 and fully vaccinated by September 2022. In the case of low transmission rate, the number of active cases would start to decrease in the early months of 2021, and the epidemic would be almost eradicated in early 2022, while in the cases with medium to high transmission rate the epidemic curve would reach its peak around May 2021 and would be close to zero by mid-2022. Our simulations show that keeping a low transmission rate (by wearing face masks, complying with social/physical distancing, etc.) is the most effective method to reduce the death toll. For example, reducing the transmission rate to 80% its baseline value results in 130 000 less deaths, while doubling the AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 287 vaccination rate does not yield a significant reduction in death toll. Also, decreasing the transmission rate is more effective to reduce the control reproduction number and achieve herd immunity than deploying vaccines with higher efficacy rates. Our model has certain limitations that could affect the results presented in this study. Firstly, our model assumes that the time between first and second doses and the efficacy rate are the same for all vaccines applied in the population; in reality, several vaccines (including single-dose vaccines) may be applied simultaneously in the same country, which would yield different parameter values for each of them. Secondly, we assumed that the protection against COVID-19 provided by vaccines does not wane over time since research about the immunity waning rate is still ongoing. Thirdly, we assume that the contact rate remains constant throughout the simulations, when it could actually change at different times due to lifting of restrictions and lockdown in the population. Lastly, the emergence of different variants of SARS-CoV-2 could be taken into account to consider the increased transmissibility and possible reduction of vaccine efficacy against these variants. These shortcomings will be addressed in future studies. References [1] M. E. Alexander, C. Bowman, S. M. Moghadas, R. Summers, A. B. Gumel and B. M. Sahai, A vaccination model for transmission dynamics of influenza, SIAM Journal on Applied Dynamical Systems 3 (2004), no. 4, 503–524. [2] M. M. Alvarez, S. Bravo-González and G. Trujillo-de Santiago, Modeling vaccination strategies in an Excel spread- sheet: increasing the rate of vaccination is more effective than increasing the vaccination coverage for containing COVID-19, PLoS One 16 (2021), no. 7, e0254430. [3] R. M. Anderson and R. M. May, Infectious diseases of humans: Dynamics and control, Oxford University Press, 1991. [4] J. Arino, C. C. McCluskey and P. van den Driessche, Global results for an epidemic model with vaccination that exhibits backward bifurcation, SIAM Journal on Applied Mathematics 64 (2003), no. 1, 260–276. [5] E. Aruffo, P. Yuan, Y. Tan, E. Gatov, E. Gournis, S. Collier, N. Ogden, J. Bélair and H. Zhu, Community structured model for vaccine strategies to control COVID19 spread: a mathematical study, medRxiv (2021). [6] U. Avila-Ponce de León, Á. G. C. Pérez and E. Avila-Vales, A data driven analysis and forecast of an SEIARD epidemic model for COVID-19 in Mexico, Big Data and Information Analytics 5 (2020), no. 1, 14–28. [7] U. Avila-Ponce de León, Á. G. C. Pérez and E. Avila-Vales, An SEIARD epidemic model for COVID-19 in Mexico: mathematical analysis and state-level forecast, Chaos, Solitons & Fractals 140 (2020), 110165. [8] B. Buonomo and R. Della Marca, Oscillations and hysteresis in an epidemic model with information-dependent imperfect vaccination, Mathematics and Computers in Simulation 162 (2019), 97–114. [9] L.-M. Cai, Z. Li and X. Song, Global analysis of an epidemic model with vaccination, Journal of Applied Mathematics and Computing 57 (2018), 605–628. [10] M. Dashtbali and M. Mirzaie, A compartmental model that predicts the effect of social distancing and vaccination on controlling COVID-19, Scientific Reports 11 (2021), 8191. [11] M. De la Sen and A. Ibeas, On an SE(Is)(Ih)AR epidemic model with combined vaccination and antiviral controls for COVID-19 pandemic, Advances in Difference Equations 2021 (2021), 92. [12] M. L. Diagne, H. Rwezaura, S. Y. Tchoumi and J. M. Tchuenche, A mathematical model of COVID-19 with vacci- nation and treatment, Computational and Mathematical Methods in Medicine 2021 (2021), 1250129. [13] A. Fuady, N. Nuraini, K. K. Sukandar and B. W. Lestari, Targeted vaccine allocation could increase the COVID-19 vaccine benefits amidst its lack of availability: a mathematical modeling study in Indonesia, Vaccines 9 (2021), no. 5, 462. [14] J. R. Gog, E. M. Hill, L. Danon and R. N. Thompson, Vaccine escape in a heterogeneous population: insights for SARS-CoV-2 from a simple model, Royal Society Open Science 8 (2021), no. 7, 210530. [15] Government of Mexico, Vacuna Covid, http://vacunacovid.gob.mx/wordpress/informacion-de-la-vacuna/, ac- cessed 11 August 2021. [16] A. B. Gumel, C. C. McCluskey and J. Watmough, An SVEIR model for assessing potential impact of an imperfect anti-SARS vaccine, Mathematical Biosciences & Engineering 3 (2006), no. 3, 485. [17] R. J. Harris, J. A. Hall, A. Zaidi, N. J. Andrews, J. K. Dunbar and G. Dabrera, Effect of vaccination on household transmission of SARS-CoV-2 in England, New England Journal of Medicine (2021), NEJMc2107717. http://vacunacovid.gob.mx/wordpress/informacion-de-la-vacuna/ 288 Á. G. C. PÉREZ AND D. A. OLUYORI [18] Johns Hopkins CSSE, 2019 Novel Coronavirus COVID-19 (2019-nCoV) Data Repository, https://github.com/ CSSEGISandData/COVID-19/, accessed 26 February 2021. [19] C. M. Kribs-Zaleta and M. Martcheva, Vaccination strategies and backward bifurcation in an age-since-infection structured model, Mathematical Biosciences 177 (2002), 317–332. [20] Y. Li, L. Ge, Y. Zhou, X. Cao and J. Zheng, Toward the impact of non-pharmaceutical interventions and vaccination on the COVID-19 pandemic with time-sependent SEIR model, Frontiers in Artificial Intelligence 4 (2021), 648579. [21] M. Martcheva, On the mechanism of strain replacement in epidemic models with vaccination, Current Developments in Mathematical Biology (K. Mahdavi, R. Culshaw and J. Boucher, eds.), World Scientific, 2007, pp. 149–172. [22] S. Moore, E. M. Hill, M. J. Tildesley, L. Dyson and M. J. Keeling, Vaccination and non-pharmaceutical interventions for COVID-19: a mathematical modelling study, The Lancet Infectious Diseases 21 (2021), no. 6, 793–802. [23] A. Olivares and E. Staffetti, Uncertainty quantification of a mathematical model of COVID-19 transmission dynamics with mass vaccination strategy, Chaos, Solitons & Fractals 146 (2021), 110895. [24] M. Roser, H. Ritchie, E. Ortiz-Ospina and J. Hasell, Coronavirus pandemic (COVID-19), Our World in Data (2020). [25] C. M. Saad-Roy, S. E. Morris, C. J. E. Metcalf, M. J. Mina, R. E. Baker, J. Farrar, E. C. Holmes, O. G. Pybus, A. L. Graham, S. A. Levin, B. T. Grenfell and C. E. Wagner, Epidemiological and evolutionary considerations of SARS-CoV-2 vaccine dosing regimes, Science 372 (2021), no. 6540, 363–370. [26] G. A. Salcedo-Varea, F. Peñuñuri, D. González-Sánchez and S. Dı́az-Infante, Optimal piecewise constant vaccination and lockdown policies for COVID-19, medRxiv (2021). [27] F. Saldaña and J. X. Velasco-Hernández, The trade-off between mobility and vaccination for COVID-19 control: a metapopulation modelling approach, Royal Society Open Science 8 (2021), no. 6, 202240. [28] A. Scherer and A. McLean, Mathematical models of vaccination, British Medical Bulletin 62 (2002), no. 1, 187–199. [29] H. Sjödin, A. Wilder-Smith, S. Osman, Z. Farooq and J. Rocklöv, Only strict quarantine measures can curb the coronavirus disease (COVID-19) in Italy, 2020, Eurosurveillance 25 (2020), no. 13, 2000280. [30] F. Sulayman, F. A. Abdullah and M. H. Mohd, An SVEIRE model of tuberculosis to assess the effect of an imperfect vaccine and other exogenous factors, Mathematics 9 (2021), no. 4, 327. [31] B. Tang, X. Zhang, Q. Li, N. Bragazzi, D. Golemi-Kotra and J. Wu, The minimal COVID-19 vaccination coverage and efficacy to compensate for potential increase of transmission contacts and increased transmission probability of the emerging strains, Research Square (2021). [32] S. Y. Tchoumi, H. Rwezaura and J. M. Tchuenche, Dynamic of a two-strain COVID-19 model with vaccination, Research Square (2021). [33] P. van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compart- mental models of disease transmission, Mathematical Biosciences 180 (2002), no. 1-2, 29–48. [34] T. N. Vilches, K. Zhang, R. Van Exan, J. M. Langley and S. M. Moghadas, Projecting the impact of a two-dose COVID-19 vaccination campaign in Ontario, Canada, Vaccine 39 (2021), no. 17, 2360–2365. [35] M. Voysey, S. A. Costa Clemens, S. A. Madhi, L. Y. Weckx, P. M. Folegatti, P. K. Aley, B. Angus, V. L. Baillie, S. L. Barnabas, Q. E. Bhorat, et al., Safety and efficacy of the ChAdOx1 nCoV-19 vaccine (AZD1222) against SARS-CoV-2: an interim analysis of four randomised controlled trials in Brazil, South Africa and the UK, The Lancet 397 (2021), no. 10269, 99–111. [36] M. Voysey, S. A. Costa Clemens, S. A. Madhi, L. Y. Weckx, P. M. Folegatti, P. K. Aley, B. Angus, V. L. Baillie, S. L. Barnabas, Q. E. Bhorat, et al., Single-dose administration and the influence of the timing of the booster dose on immunogenicity and efficacy of ChAdOx1 nCoV-19 (AZD1222) vaccine: a pooled analysis of four randomised trials, The Lancet 397 (2021), no. 10277, 881–891. [37] P. Widyaningsih, P. Candrawati, Sutanto and D. R. S. Saputro, Maternal antibody Susceptible Vaccinated Infected Recovered (MSVIR) model for tetanus disease and its applications in Indonesia, Journal of Physics: Conference Series 1306 (2019), 012002. [38] P. Wintachai and K. Prathom, Stability analysis of SEIR model related to efficiency of vaccines for COVID-19 situation, Heliyon 7 (2021), no. 4, e06812. [39] World Health Organization, The Oxford/AstraZeneca COVID-19 vaccine: what you need to know, https://www.who.int/en/news-room/feature-stories/detail/the-oxford-astrazeneca-covid-19-vaccine-what- you-need-to-know, accessed 25 March 2021. [40] World Health Organization, WHO advisory committee on variola virus research: report of the thirteenth meeting, Tech. report, World Health Organization, 2011. https://github.com/CSSEGISandData/COVID-19/ https://github.com/CSSEGISandData/COVID-19/ https://www.who.int/en/news-room/feature-stories/detail/the-oxford-astrazeneca-covid-19-vaccine-what- you-need-to-know AN EXTENDED SEIARD MODEL FOR COVID-19 VACCINATION 289 Corresponding author, Facultad de Matemáticas, Universidad Autónoma de Yucatán, Mérida, Yucatán, Mexico E-mail address: agcp26@hotmail.com Department of Mathematics, School of Physical Science, Ahmadu Bello University, Zaria, Kaduna State, Nigeria E-mail address: oluyoridavid@gmail.com 1. Introduction 2. Model formulation 3. Theoretical analysis 3.1. Impact of vaccination coverage on the control reproduction number 4. Numerical simulations 4.1. Data fitting and estimation of parameters 4.2. Simulations for the model with vaccination 5. Conclusion References