Article On a Discrete SEIR Epidemic Model with Two-Doses Delayed Feedback Vaccination Control on the Susceptible Manuel De la Sen 1, * , Santiago Alonso-Quesada 1 , Asier Ibeas 2 and Raul Nistal 1 1 Faculty of Science and Technology, Institute of Research and Development of Processes IIDP, University of the Basque Country, Barrio Sarriena, 48940 Leioa, Spain; santiago.alonso@ehu.eus (S.A.-Q.); raul.nistal@gmail.com (R.N.) 2 Department of Telecommunications and Systems Engineering, Universitat Autònoma de Barcelona, UAB, 08193 Barcelona, Spain; Asier.Ibeas@uab.cat * Correspondence: manuel.delasen@ehu.eus Abstract: A new discrete susceptible-exposed-infectious-recovered (SEIR) epidemic model is pre- sented subject to a feedback vaccination effort involving two doses. Both vaccination doses, which are subject to a non-necessarily identical effectiveness, are administrated by respecting a certain mutual delay interval, and their immunity effect is registered after a certain delay since the second dose. The delays and the efficacies of the doses are parameters, which can be fixed in the model for each concrete experimentation. The disease-free equilibrium point is characterized as well as its stability properties, while it is seen that no endemic equilibrium point exists. The exposed subpopulation is supposed to be infective eventually, under a distinct transmission rate of that of the infectious   subpopulation. Some simulation examples are presented by using disease parameterizations of the COVID-19 pandemic under vaccination efforts requiring two doses. Citation: De la Sen, M.; Alonso-Quesada, S.; Ibeas, A.; Nistal, R. On a Discrete SEIR Keywords: discrete epidemic model; delayed feedback vaccination control; vaccination doses; Epidemic Model with Two-Doses COVID-19 pandemic Delayed Feedback Vaccination Control on the Susceptible. Vaccines 2021, 9, 398. https://doi.org/ 10.3390/vaccines9040398 1. Introduction Typical formulations used to describe epidemic models are based on either differential Academic Editors: Khalid Hattaf, or difference equations. In that way, the basic reproduction number and its physical and Bapi Pahar and Giuseppe La Torre biological insight are discussed in [1] which is related to pertussis and measles descriptions. In addition, feedback vaccination laws have been developed using techniques such as Received: 9 March 2021 sliding-mode control or linear or impulsive feedback vaccination [2,3]. The transient Accepted: 15 April 2021 evolution of epidemic diseases is also an important issue for properly describing the day-to- Published: 18 April 2021 day time-transmission levels and the appropriate eventual interventions to perform since the stability properties are more related to the stationary states, typically the disease-free Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in and the endemic ones. For instance, in [4], an analytic methodology is given to predict and published maps and institutional affil- monitor the dates of maximum hospital occupancy of beds. The differential and difference iations. models have also been corrected with other powerful analysis techniques. In that way, the bifurcation analysis and the stability of a fractional order susceptible-infectious-recovered (SIR) epidemic model with delay has been discussed in [5]. On the other hand, discretized and discrete-time epidemic models have been proposed in the background literature. The different approaches can basically consist of the discretization of continuous-time- Copyright: © 2021 by the authors. based models by numerical methods or in the development of discrete models based in Licensee MDPI, Basel, Switzerland. difference equations. See, for instance, [6–12] and some of the references therein. It turns This article is an open access article distributed under the terms and out that since the relevant time of the dynamics evolution in epidemics is relatively long, conditions of the Creative Commons typically for instance, on the orders of days or weeks, discrete-oriented epidemic models Attribution (CC BY) license (https:// can be found to be appropriate for describing and monitoring the infection evolution creativecommons.org/licenses/by/ through time. An extended SEIR model which incorporates its usual subpopulations, 4.0/). the asymptomatic, and the dead-infective subpopulations has been proposed in [6]. That Vaccines 2021, 9, 398. https://doi.org/10.3390/vaccines9040398 https://www.mdpi.com/journal/vaccines Vaccines 2021, 9, 398 2 of 26 model has been proposed as appropriate for the Ebola disease. On the other hand, a multistaged SIR was discussed in [7]. Such a model considers several coupled layers of infectious subpopulation for a coupled disease transmission from each layer to the adjacent ones. Moreover, discrete susceptible-infectious-recovered-susceptible (SIRS) models have also been studied in heterogeneous networks [8]. The vaccination effort can be considered either as an external forcing term or as a generator of a new subpopulation, the so-called vaccinated one. Such a subpopulation becomes dynamically coupled to the remaining ones in the model rather than as a specific forcing control [9]. Other types of epidemic models, such as, for instance, discretized susceptible-infectious-recovered (SIR)-type ones, or susceptible-infectious-susceptible (SIS)-type ones, have been proposed in [10–13] and some of the references therein. Recently, a lot of research is being dedicated to studying and monitoring the new COVID-19 pandemic using registered data, like infection and detection tests, hospital bed occupancies, and mortality-related records. Normally, such data are updated in discrete time, typically, day to day or week to week. Therefore, discrete epidemic models have been found appropriate for processing such data. See, for instance, [14–36] and some related references therein. Studies on particular data for different countries or regions can be found in the literature related to COVID-19, sometimes related to public interventions, such as quarantines, isolation measures, lockdowns, use of masks, social distance rules, etc. See, for instance, those concerned with Saudi Arabia [16,17], Madrid capital town, metropolitan area and surrounding administrative area [18,19], India [23,24], Italy [25], United States [26], Canada and several of its provinces [28], Switzerland, [29], Brazil [30], etc. In addition, the analysis of data has been sometimes accompanied with mathematical analysis techniques on the pandemic evolution related to public interventions or mathematically founded analysis of the obtained data. In that way, the impact of lockdowns is investigated in [17], while the effects of total or partial quarantines are investigated in [18] for a SEIAR model, which incorporates the asymptomatic subpopulation to the typical SEIR model by consid- ering the isolated population as removed either from the infectious individuals or from the susceptible ones. In [19], a more general model of potential usefulness in the description and monitoring of COVID-19 has been proposed, discussed, and tested with recorded real data. Such a model includes three different infectious subpopulations, namely, the slightly infections, the hospitalized ones and the ones staying at the intensive care units are considered. On the other hand, the implementation of control rules oriented either toward reducing the number of exposed individuals or toward increasing the number of treated individuals is proposed and discussed in [20] while impulsive optimal control tech- niques are developed in [21]. In particular, the analysis technique proposed in [20] relies on the fact that the epidemics is endemic. In [23], the inadequacy of the implementation of open-loop (i.e., without using feedback) controls is emphasized contrarily to the use of closed-loop controls like, for instance, sliding mode-based control or other feedback laws. Moreover, it has been proved in [27] that the suppression strategies are appropriate, provided that they are sufficiently strong while taken through prompt decisions, whereas the mitigation strategies can fail because of eventual unfavorable combinations of delays, unstable dynamics, and uncertainties in the model. This paper proposes and investigates a new discrete SEIR model subject to a linear two-stage delayed feedback vaccination effort having in mind that some of the recently ap- proved vaccines for COVID-19 require two doses for increasing their average effectiveness. Such doses are administered to the susceptible subpopulation with a delay period, and their potential benefits on immunity appear several days after the administration of the second doses. In fact, the model considers the injection of two doses with different delays and eventual average different effectiveness. The proposed model considers also that the exposed have a transmission rate exposed-susceptible which may be eventually distinct from the infectious-susceptible transmission rate. The underlying idea is that in some infectious diseases, such as, for instance, in the COVID-19 pandemic, there are contagions from both exposed to susceptible and from infectious to susceptible. The delays of the Vaccines 2021, 9, 398 3 of 26 vaccination and expected achievable immunity as well as the transmission rates are model parameters which can be updated for different experiments, [37–39]. The paper is organized as follows: Section 2 is devoted to present the new mentioned proposed discrete susceptible-exposed-infectious-recovered (SEIR) epidemic model with delayed double-dose linear feedback vaccination. Section 3 discusses the non-negativity and boundedness of the solution under any given finite non-negative initial conditions, as well as the existence and components of the disease-free equilibrium point and its stability properties. It is also proved that the endemic equilibrium point does not exist for the proposed model in Section 4. Section 5 presents and discusses some examples of the proposed model related to the evolution of the COVID-19 pandemic. Finally, a set of illustrative concluding remarks ends the paper. Some auxiliary technical results are proved in Appendix A which are supported by general necessary mathematical results given or proved in [40–44]. 2. The Discrete SEIR Epidemic Model Subject to Two Vaccination Doses Note that the rationale of the sampling period interpretation is unity, typically one day or one week for a correct practical use of the model. The model parameters should be expressed in values of dimensionality being the inverse of the sampling period units. The SEIR epidemic model equations may be rewritten equivalently as follows: Sk+1 = ( ak − β k Ik − βek Ek )Sk − V2,k−d2 − V1,k−d1 −d2 (1) Ek+1 = (1 − µ) Ek + ( β k Ik + βek Ek )Sk (2) Ik+1 = (1 − γ) Ik + µEk (3) Rk+1 = Rk + γIk + V2,k−d2 + V1,k−d1 −d2 (4) V2,k−d2 = Kk−d2 ρ2 Sk−d2 (5) V1,k−d1 −d2 = Kk−d1 −d2 ρ1 Sk−d1 −d2 (6) for any integer k ∈ Z0+ = Z+ ∪ {0} and any given finite initial conditions S0 ≥ 0, E0 ≥ 0, I0 ≥ 0 and R0 ≥ 0 and Sk = Ek = Ik = Rk = 0 for k ∈ Z− , where S, E, I and R are the susceptible, exposed, infectious and recovered subpopulations, respectively. The forcing terms of Equations (5) and (6) are the two doses of vaccination on the susceptible which are generated via linear feedback and which are subject to integer delays d1 + d2 (first dose) and d2 (second dose), respectively, with respect to each current sampling instant, where min(d1 , d2 ) > 0. In the above model: – ak is the average recruitment rate proportional to the susceptible at the kth sampling instant related; for instance, to the rates of births and a is a constant reference value for the above sequence; for instance, its average over the whole time period under study and typically it might be unity. – β k and βek are, respectively, the average transmission rates of the infectious and exposed subpopulations at the k − th sampling instant. – γ is the average recovery rate. – µ is the average incubation rate. – Kk is the vaccination rate (a feedback control gain) which can be eventually depending on the sampling instants. It is assumed in the sequel that {Kk }0∞ ⊂ [0, 1]. – ρ1 , ρ1 + ρ2 are parameters in (0, 1] which quantify the average effectiveness (or efficiency) of the respective doses. In particular, ρ2 gives the extra effectiveness obtained from the injection of the second dose. In this context, ρ1 + ρ2 ∈ [0, 1] and ρ1 + ρ2 = 1 refer to the ideal situation, unattainable in practice, of 100 percent effectiveness of the combined injection of the two doses. Note that ρ1 = ρ2 = ρ1 + ρ2 = 0 refers to the worst case where the vaccination is fully superfluous. Typically, we can consider two different transmission rates for the exposed and in- fectious since the exposed are not usually identified to be allocated under quarantine or Vaccines 2021, 9, 398 4 of 26 isolation, and furthermore, it has been argued that the infective periods are, in general, of distinct time length for both stages in the case of COVID-19 pandemic. It turns out that these transmission rates can be time varying, in general, since the transmission rates can depend on the intervention measures and on the social customs in the geographic area under study. Note that, even if the vaccination is performed under the same gain K for both doses, it can be considered that the effectiveness if only the first injected dose is smaller than if both of them are injected. Therefore, it can be typically argued that ρ2 > 0 leads to ρ1 + ρ2 > ρ1 . It can also be pointed out that the proposed discrete-time model corresponds to the backward Euler discretization on the continuous –time SEIR model given by: . S(t) = ν(t)S(t) − ( β0 I (t) + β0 e E(t))S(t) − V1 0 (t − d1 − d2 ) − V2 0 (t − d2 ) . E(t) = −µ0 E(t) + ( β0 I (t) + β0 e E(t))S(t) . I (t) = µ0 E(t) − γ0 I (t) . R(t) = γ0 I (t) + V1 0 (t − d1 − d2 ) + V2 0 (t − d2 ) whose discretization is given by Sk+1 = Sk + hνk Sk − (hβ0 Ik + hβ0 e Ek )Sk − hV1,k−d1 −d2 0 − hV2,k−d2 0 = (1 + hν)Sk − (hβ0 Ik + hβ0 e Ek )Sk − hV1,k−d1 −d2 0 − hV2,k−d2 0 Ek+1 = Ek − hµ0 Ek + (hβ0 Ik + hβ0 e Ek )Sk = (1 − hµ0 ) Ek + (hβ0 Ik + hβ0 e Ek )Sk Ik+1 = Ik + hµ0 Ek − hγ0 Ik = (1 − hγ0 ) Ik + hµ0 Ek Rk+1 = Rk + hγ0 I (t) + hV1,k−d1 −d2 0 + hV2,k−d2 0 so that the discretized and the proposed discrete-time model are equivalent if the following correspondence is performed ak = 1 + hνk β = hβ0 , βe = hβ0 e V1,k−d1 −d2 = hV1,k−d1 −d2 0 , V2,k−d2 = hV2,k−d2 0 µ = hµ0 , γ = hγ0 However, the model is originally set up in discrete-time instead as the discretization of a continuous-time one. One of the daunting challenges to counteract the epidemic spreading is the design of the vaccination functions. The control theory provides a frame- work to systematically design them in order to fulfill with a prescribed control objective. Thus, an output feedback approach is adopted in this work and discussed in the sequel. However, control theory has developed many analytical tools that allows one to face this problem by using other approaches such as state feedback and feedback linearization [2], output-controllability [19], optimal control [21], impulsive control [22], sliding mode [23] and lockdown and quarantine [18,28], to cite just a few. The model (1)–(6) may be compacted after replacing the vaccination controls (5)–(6) in the susceptible and recovered subpopulations dynamics, which result in:  Sk+1 = aSk − K ρ2 Sk−d2 + ρ1 Sk−d1 −d2 + uk (7) Ek+1 = (1 − µ) Ek + ( β k Ik + βek Ek )Sk (8) Ik+1 = (1 − γ) Ik + µEk (9)  Rk+1 = Rk + γIk + K ρ2 Sk−d2 + ρ1 Sk−d1 −d2 − gk (10) uk = [ ak − a − β k ( Ik + λek Ek )]Sk + gk (11)   gk = K − K k − d2 ρ 2 Sk − d2 + K − K k − d1 − d2 ρ 1 Sk − d1 − d2 (12) for any integer k ∈ Z0+ , where λek = βek /β k and K and a are reference values for the vaccination rate gain and recruitment, for instance, the estimated average values of their Vaccines 2021, 9, 398 5 of 26 asymptotic limits if they converge. In particular, a scalar discrete function of the form of the susceptible subpopulation (7) is discussed in detail in Appendix A from the points of view of stability and convergence to a limit. Those auxiliary results based on a similar discrete Equation (7) subject to (11) are then used to discuss the stability of the proposed epidemic model. 3. Non-Negativity, Stability and Disease-Free Equilibrium Point The following result which relies on boundedness, convergence and non-negativity of the solution is direct by summing up (7) to (10) (or (1) to (4)) by considering (7) to (11) (or (1) to (6)): Theorem 1. The following properties hold: (i) Assume that N0 = S0 + E0 + I0 + R0 = 1. Then, the total population Nk = Sk + Ek + Ik + Rk = 1; ∀k ∈ Z0+ if ak ≡ 1. k −1 (ii) Nk = 1 + ∑ j=0 a j − 1 S j ; ∀k ∈ Z0+  (iii) {Sk }∞ k=0 is non-increasing, and then bounded and convergent, if V2,k−d2 + V1,k−d1 −d2 = Kk−d2 ρ2 Sk−d2 + Kk−d1 −d2 ρ1 Sk−d1 −d2 ≥ ( ak − β k ( Ik + λek Ek ) − 1)Sk ; ∀k ∈ Z0+ (13) leading to k −1 k −1 1 + Sk−1 ∑ j=0 a j − 1 ≤ Nk ≤ 1 + S0 ∑ j=0 a j − 1 ; ∀k ∈ Z0+   (14) and, if the average value of the sequence { ak }∞ k=0 is unity, then there exists the limit lim Nk = k→∞ N0 = 1. {Sk }∞ k =0 is strictly decreasing if the inequality in (13) is strict. (iv) {Sk }∞ k =0 is non-negative if and only if V2,k−d2 + V1,k−d1 −d2 = Kk−d2 ρ2 Sk−d2 + Kk−d1 −d2 ρ1 Sk−d1 −d2 ≤ ( ak − β k ( Ik + λek Ek ))Sk ; ∀k ∈ Z0+ (15) ak provided that the infectious transmission rate is sufficiently small according to β k ≤ Ik +λek Ek ; ak ∀k ∈ Z0+ which is guaranteed if β k ≤ ; ∀k ∈ Z0+ . (1+λek ) Nk (v) Assume that µ ∈ [0 , 1), γ ∈ [0 , 1) and that (15) holds. Then, any sequence trajectory solution of (1)–(4) subject to (5)–(6) is non-negative. Proof. If ak = 1; ∀k ∈ Z0+ then N0 = S0 + E0 + I0 + R0 = 1 implies directly that Nk = Sk + Ek + Ik + Rk = N0 = 1; ∀k ∈ Z0+ by summing up (1) to (4). The same calculations prove Property (ii) in the general case. To prove, Property (iii), note from (3) that, under the given constraint, {Sk }∞ k =0 is non-increasing so that (14) holds directly from Property (ii) since: S k +1 V2,k−d2 + V1,k−d1 −d2 = ak − β k ( Ik + λek Ek ) − ≤ 1; ∀k ∈ Z0+ (16) Sk Sk Note that {Sk }∞ k =0 is also bounded and convergent since it is non-increasing and its initial condition is finite by hypothesis. Property (iv) follows directly by replacing (15) with (1) subject to S0 ≥ 0 and to β k ≤ I +aλke E (guaranteed if β k ≤ ak ); ∀k ∈ Z0+ k k k (1+λek ) Nk since Vk ≥ 0; ∀k ∈ Z0+ . If, addition, µ ∈ [0 , 1), γ ∈ [0 , 1) then from (2)–(4) and { Sk } ∞ ∞ ∞ ∞ k =0 ⊂ R0+ = R+ ∪ {0} yields { Ek }k =0 ⊂ R0+ , { Ik }k =0 ⊂ R0+ and { Rk }k =0 ⊂ R0+ , since the initial conditions of all the subpopulations are non-negative. Property (v) is proved.  Vaccines 2021, 9, 398 6 of 26 Remark 1. The non-negativity property of the susceptible of Theorem 1 (iv) is a necessary constraint for a well-posedness of the model. Note that the condition of sufficiently smaller transmission rate ak of the Infectious β k ≤ ; ∀k ∈ Z0+ , which guarantees according to Theorem 1 (iv) (1+λek ) Nk that the sequence of susceptible is non-negative if S0 ≥ 0 and the vaccination satisfies (15), is reasonable since the transmission rate decreases as the total population increases in the popular class of true-mass action epidemic models. The property of the susceptible sequence being non-increasing [Theorem 1 (iii)] describes appropriately the disease growing since it starts until its first maximum peak. Remark 2. It turns out the need of the assumption {Kk }0∞ ⊂ [0, 1] for the vaccination gain sequence since the daily vaccination (or, in general, the one for the used sampling period in the model parameterization) is proportional to the susceptible subpopulation. Therefore, note that the condition (13) for {Sk }∞ k=0 to be non-increasing requires the necessary condition: ρ2 Sk−d2 + ρ1 Sk−d1 −d2 ≥ ( ak − β k ( Ik + λek Ek ) − 1)Sk ; ∀k ∈ Z0+ (17) to be fulfilled by the vaccination gain unity, which is the maximum possible one, provided that the existing stock vaccines covers such a need. Note that the above condition always holds if { Sk } ∞ k =0 is non-negative and ak ≤ e ∞  β k Ik + λk Ek + 1; ∀k ∈ Z0+ and it also holds if {Sk }k=0 is e non-negative, ak > β k Ik + λk Ek + 1; ∀k ∈ Z0+ . This second condition is guaranteed by under unity vaccination gain:  (ρ1 + ρ2 )min Sk−d2 , Sk−d1 −d2 ≥ 1; ∀k ∈ Z0+ (18) ak − β k Ik + λek Ek − 1 Sk     min Sk−d2 , Sk−d −d2 1 In practice, Sk is slightly greater than unity if the susceptible sequence is decreasing since d1 and d1 + d2 are typically delays of one and three weeks, in the case of SARS- CoV-2 which do not affected significantly the variation of the susceptible population levels. In addition, ak − β k Ik + λek Ek − 1 is typically close to zero for small infection levels compared to the susceptible population values. Therefore, for reasonable vaccine efficacies, (18) holds if e  ak > β k Ik + λk Ek + 1; ∀k ∈ Z0+ . Note that the condition (15) for the susceptible sequence {Sk }∞ k =0 to be non-negative requires the necessary condition to be fulfilled by some maximum possible vaccination gain rate K ≤ 1 (since this feature would ensure that it holds for any lower effort as well), that is, ak − β k Ik + λek Ek Sk ≥ 0; ∀k ∈ Z0+ which holds if S0 ≥ 0 if β k ≤ I +aλke E ; ∀k ∈ Z0+ k k k ak which is guaranteed if β k ≤ ; ∀k ∈ Z0+ . (1+λek ) Nk The existence of a disease- free equilibrium point is discussed in the subsequent result: ak ak Theorem 2. Assume that ρ1 + ρ2 ≤ 1, β k ≤ Ik +λek Ek (guaranteed if β k ≤ ); ∀k ∈ Z0+ , (1+λek ) Nk { ak }∞ ∞ k =0 → a (≤ 1) , { Kk }k =0 (⊂ [0, 1]) → K and assume also that the constraints (13) and (15)  T hold. Then, there is a disease-free equilibrium point xd f = Sd f , 0, 0, Rd f = Nd f − Sd f , being in general dependent on the initial conditions, with Sd f ≥ 0 and Rd f = Nd f − Sd f , where { Nk }∞ a k =0 → Nd f , such that Sd f = 0 if K 6 = ρ1 +ρ2 . Proof. Note that {Kk }∞ k =0 (⊂ [0, 1]) → K implies that K ∈ [0, 1]. Note also that Ed f = 0 and Id f = 0 are trivial solutions of (2) and (3) if {Sk }∞ k =0 is bounded. It follows furthermore that {Sk }∞ k =0 is non-increasing and bounded from Theorem 1 (iii), { ak }∞ k =0 → a (≤ 1) and ∞ (13) holds. As a result, if furthermore, (15) holds then {Sk }k=0 (⊂ R0+ ) → Sd f (≥ 0) with Sd f = 0 if K 6= ρ +a ρ2 subject to the assumed constraint ρ1 + ρ2 ≥ 1 since K ∈ [0, 1]. Note 1 that Nk+1 = ak Sk + Ek + Ik + Rk = Nk + ( ak − 1)Sk (19) Vaccines 2021, 9, 398 7 of 26 implies that Nk+1 − Nk = ( ak − 1)Sk and { ak }∞ ∞ k =0 → 1 implies that { Nk }k =0 → Nd f and ∞ ∞ ∞ { Rk }k=0 → Rd f = Nd f − Sd f since { Ek }k=0 → Ed f = 0 and { Ik }k=0 → Id f = 0 .  The subsequent result gives conditions for both the infection and the susceptibility to asymptotically vanish implying convergence of the solution to a disease-free equilibrium point without susceptible individuals. Some supportive auxiliary results used in its proof are given in Appendix A. Theorem 3. Assume that the constraint (15) holds, so that any solution sequence is non-negative for non-negative initial conditions. Assume also that the first assumptions below, and, furthermore, at least one of the assumptions 2 to 5 below also holds: ak ak (1) µ ∈ [0, 1), γ ∈ [0, 1), ρ1 + ρ2 ≤ 1, Kk ≤ Ik +λek Ek (guaranteed if β k ≤ (1+λe ) Nk ); k ∀k ∈ Z0+ , { ak }∞ ∞ e ∞ e ∞  k =0 → a (≤ 1) , { β k }k =0 → β , λk k =0 → λ , { Kk }k =0 → K . (2) Either the constraint (13) holds with strict inequality for all k ∈ Z0+ and K 6= ρ +a ρ2 . 1 (3) K ∈ (0, min(K̄, 1)), where: a−1     1+a K̄ = sup y ∈ , : y < 1/ f (y) (20) ρ1 + ρ2 ρ1 + ρ2 v u (ρ2 (cos(θd2 ) − 1) + ρ1 (cos(θ (d1 + d2 )) − 1))2 + (ρ2 sin(θd2 ) + ρ1 sin(θ (d1 + d2 )))2 u f (y) = sup t (21) θ ∈(0,2π ) 1 + ( a − (ρ1 + ρ2 )y)2 − 2( a − (ρ1 + ρ2 )y)cosθ and     |Kk − K | = o k x̂k k−1 ; | ak − a| = o k x̂k k−1 (22) T where x̂k = xk , xk−1 , . . . , xk−d2 , . . . , xk−d1 −d2 . " ! ρ2 e−iθd2 + ρ1 e−iθ (d1 +d2 ) (4) a ∈ (0, 1), K ∈ 0, 1/ sup (23) θ ∈[0,2π ) eiθ − a and (22) holds.    − d1  c1 − a c1 (5) a ∈ (0, 1), K ∈ 0,  (24) ρ1 + ρ2 c2 Then, there exists a globally asymptotically stable disease-free equilibrium point xd f = T and { Nk }∞  0, 0, 0, Rd f = Nd f k =0 → Nd f . Proof. The above assumptions 1 and 2 guarantee that {Sk }∞ ∞ k =0 → 0 , and that { Sk }k =0 is bounded, from Theorem 1 and Theorem 2 for any given non-negative finite initial conditions. Note also that, from (2) and (3), { Ek }∞ ∞ k =0 → 0 and { Ik }k =0 → 0 since µ ∈ [0, 1) and γ ∈ [0, 1). In addition, any of the assumptions 3 or 4, together with the properties { Ek }∞ k =0 → 0 and { Ik }∞ k =0 → 0 guarantees that { S } ∞ k k =0 → 0 , and { S } ∞ k k =0 is bounded, for any given non-negative finite initial conditions from Theorem A1 [(i)–(ii)] and Corollary A1 of Appendix A. The same result also holds under the assumption 5 according to Theorem A3 (i) of Appendix A.  T Thus, the disease-free equilibrium point xd f = 0, 0, 0, Rd f = Nd f exists under the joint assumptions 1 and 2 and under the assumption 1 jointly with one of the assumptions 3 to 5. Now, note that Since { β k Sk }∞ k =0 → βSd f = 0 , Equations (2) and (3) can be compacted as follows: Vaccines 2021, 9, 398 8 of 26 1 − µ + λek β k Sk      Ek+1 β k Sk Ek = Ik+1 µ 1−γ Ik  "  e    # 1 − µ + λe βSd f β λ k Sk − λ e Sd f    βSd f Ek β Sk − Sd f Ek = + (25) µ 1−γ Ik 0  0 Ik βλek Sk βSk      1−µ 0 Ek Ek = + ; k ∈ Z0+ µ 1−γ Ik 0 0 Ik   1−µ 0 Note that the spectral radius r of the matrix is max(1 − µ, 1 − γ) < 1 µ 1−γ since µ ∈ [0, 1) and γ ∈ [0, 1) and the above matrix is a convergent matrix. Since there is some matrix norm of arbitrary close value to the spectral radius there  exist some positive  1−µ 0 real constants ε ≤ δ0 − 1 + min(µ, γ) and r ≤ δ0 < 1 such that ≤ δ0 . µ 1−γ In addition, since { β k Sk }∞ e ∞ e  k =0 → βSd f = 0 and λk k =0 → λ and a decreasing nonnega- o∞ tive real sequence {δ1k }∞ n k =0 → 0 having a strictly decreasing subsequence δ1Nn →0 k k =0 such that for any given δ1Nn , there exists a non-negative integer Nnk such that k βλek Sk βSk    ∞ ≤ δ1Nk ; j = Nnk , Nnk +1 , . . . , Nnk+1 −1 such that Nnk k=0 is strictly 0 0 increasing. Thus, one gets from (25) that: n k +1 −1 " # " # ENn Nn − Nnk ENn −1− j ∑ n k +1 ≤ δ0 k+1 k + δ1Nk δ0 k+1 ; k ∈ Z0+ (26) INn INn j=nk k +1 k ( " # )∞ ENn so that k has a strictly decreasing subsequence which converges to zero INn k k =0 for any given finite non-negative initial conditions. In view of (2)–(3), that property also implies that { Ek }∞ ∞ k =0 → 0 and { Ik }k =0 → 0 for any given non-negative finite initial conditions and then (25) is globally asymptotically stable. Since it has been already proved T that {Sk }∞  k =0 → 0 then xd f = 0, 0, 0, Rd f = Nd f is globally asymptotically stable and { Nk }∞ k =0 → Nd f .  The next result gets condition for which the stationary model has a stable infective substate, irrespective of the susceptible being zero at the equilibrium point or not. It is found that there exists a critical transmission rate under which the stationary substate of the infective model is globally asymptotically stable. Under the extra conditions from Theorem 3, the critical transmission rate becomes infinity since the susceptible subpopulation is zero at the disease-free equilibrium point. As a result, the disease-free equilibrium point is globally asymptotically stable for any value of the transmission rate, and there is no endemic attractor. Theorem 4. Assume µ ∈ [0, 1) and γ ∈ [0, 1). Then, the stationary infective substate (25) µγ has stable characteristic roots if β ∈ [0, β c ), where β c = (µ+λe γ) S so that it is globally df asymptotically stable. In addition, if the conditions of Theorem 2 hold with K 6= ρ +a ρ2 then β c = +∞. 1 The stable solution sequence is guaranteed to be, furthermore, non-negative for any non- negative finite initial conditions if all the conditions of Theorem 3 hold.  e ∞ Proof. If { β k Sk }∞ e k =0 → βSd f = s and λk k =0 → λ , then the characteristic Equation of (25) is: λ2 + bλ + c = 0 (27) Vaccines 2021, 9, 398 9 of 26 where b = b ( s ) = µ + γ − (2 + λ e s ) (28) c = c(s) = (1 − µ + λe s)(1 − γ) − µs (29) Note that b(0) = µ + γ − 2 < 0 and c(0) = (1 − µ)(1 − γ) > 0. The coefficients of the polynomial defining the characteristic equation for s = 0, ordered in decreasing order, have two changes of sign. According to Descartes’ rule of signs, such a polynomial has either two or zero positive real roots. Note that by simple inspection or, equivalently, according to the unforced part of (25), the characteristic Equation (27) has in fact two positive real characteristic roots for s = 0 within the unit complex circle centered at zero which are λ1 = 1 − µ and λ2 = 1 − γ. Now, by inspection, b(s) = b(0) − λe s < 0 for all s ≥ 0. It is now discussed a valid range of non-negative values of s which guarantees that c(s) = c(0) + (λe (1 − γ) − µ)s remains to be positive. Write by convenience s = µ+λe γ . σµγ Then, c(s) > 0 if and only if ĉ(σ ) > 0, where: ĉ(σ ) = (µ + λe γ)c(s) = (µ + λe γ)c(0) − σµγ(µ − λe (1 − γ)) > 0 for any real σσ if λe ≥ µ/(1 − γ) with the eventual negative real values having no interest since s ≥ 0 and, if λe < µ/(1 − γ) then ĉ(σ ) > 0 if µ + λe γ (1 − µ)(1 − γ) µ + λe γ 1−µ−γ   σ< = + 1 = 1 + σ0 µ − λ e (1 − γ ) µγ µ + λe γ − λe µγ for some real σ0 > 0. In both cases, that is, irrespective of λe , c(s) > 0 if s = µ+λe γ and σµγ σ ∈ [0, 1]. Since there are two real positive characteristic roots for s = 0 from Descartes rule of signs (since b(0) < 0 and c(0) > 0) and since the characteristic roots are continuous σµγ functions of the argument s and b(s) < 0 and c(s) > 0 for s = µ+λe γ with σ ∈ [0, 1], µγ one concludes that the characteristic roots are real (and positive) if s = βSd f < µ+λe γ . It remains now to prove that those real positive roots are within the open unit circle. The two (real) within the complex open circle or radius unity centered at the origin of the complex plane if and only if both roots of (27), subject to (28)–(29), are within the open unit circle centered at zero in the complex plane, namely if and only if: q 2− µ − γ + λ e s ± (2−µ−γ+λe s)2 +4(µs−(1−µ+λe s)(1−γ)) −1 < λ1,2 = 2 <1 (30) equivalently, if and only if, q µ + γ − λe s − 4< − λe s2 + (γ − µ)2 + 4µs + 2λe s(γ − µ) 2 q (31) ≤ λe2 s2 + (γ − µ)2 + 4µs + 2λe s(γ − µ) < µ + γ − λe s equivalently, if and only if, q λe s2 + (γ − µ)2 + 4µs + 2λe s(γ − µ) < min(µ + γ − λe s , 4 + λe s − µ − γ) = µ + γ − λe s 2 (32) since the last equality follows directly since µ, γ ∈ (0, 1). Taking the squares in both sides of (32) yields that it is equivalent to 2 2 λe s2 + (γ − µ)2 + 4µs + 2λe s(γ − µ) < λe s2 + (µ + γ)2 − 2λe s(µ + γ) (33) and also equivalent to (γ − µ)2 + 4µs + 2λe s(γ − µ) < (µ + γ)2 − 2λe s(µ + γ) (34) and, furthermore, equivalent to (γ − µ)2 + 4µs + 4λe sγ < (µ + γ)2 (35) Vaccines 2021, 9, 398 10 of 26 which holds if and only if µγ s = βSd f < (36) µ + λe γ Since s = βSd f ≥ 0, the constraint (36) implies that the characteristic Equation (27) of µγ the stationary infective substate (25) has stable roots if β ∈ [0, β c ), where β c = (µ+λe γ) S . df Now, assume that all the conditions of Theorem 2 hold. Since Sd f = 0 if K 6= ρ +a ρ2 1 then β c = +∞ if K 6= ρ +a ρ2 and the stationary infective substate (25) has stable roots if 1 β ∈ [0, +∞ ). In addition, any solution is non-negative for all time if (15) and assumption 1 plus any of the four remaining assumptions 2 to 5 of Theorem 3 holds (according to Theorem 3). Thus, the joint non-negativity and global asymptotic stability of the disease-free stationary solution, under any finite non-negative initial conditions, are guaranteed for β ∈ [0, +∞ ) since Sd f = 0 (Theorem 3).  Note that the Descartes rule of signs used in the proof of Theorem 4 is supported by the fact that for s = 0 both characteristic roots are real, positive and stable, so that the possibility of having no positive root is excluded by such a value of the argument s, and the fact that the characteristic roots are continuous functions of s, so that they continue to be within the positive region of the unity circle centered at zero as the argument increases. Note also that the proposed model is claimed for its usefulness for short-term predic- tions in the evolution phase when the disease is blowing up. It is neither considered that vaccination is available for use nor that there is immunity lost allowing to increase again the susceptible numbers after a certain delay. Therefore, the evolution of the susceptible subpopulation is given by a decreasing sequence. It is now proved that the proposed model does not have an endemic equilibrium point. 4. Nonexistence of Endemic Equilibrium Point It is now discussed under which reasonable conditions an endemic equilibrium point xend = (Send , Eend (6= 0), Iend (6= 0) , Rend ) T either exists or it does not exist. As a final result, it is concluded that the endemic equilibrium point cannot exist. Consider the following possible cases: Case 1 (Send > 0). Combining (2) and (3) for stationary limiting model parameters for an assumed to exist endemic equilibrium point leads to provided that Eend 6= 0 and Send 6= 0: µ Iend = Eend (37) γ γµ Send = >0 (38) β(µ + λe γ) which, replaced in (1) and after cancelling common factors with Send > 0, results in the resulting constraint: Send = [ a − β( Iend + λe Eend ) − K (ρ1 + ρ2 )] Send (39) leads to a − 1 − K ( ρ1 + ρ2 ) Eend = γ (40) β(µ + γλe ) if a 6= β( Iend + λe Eend ) + K (ρ1 + ρ2 ) (41) a −1 so that Eend > 0 (and also Iend > 0 from (37)) if and only if K < ρ1 + ρ2 ≤ 1 and 1 < a −1 a ≤ 1 + ρ1 + ρ2 ≤ 2. If K = and 1 ≤ a ≤ 1 + ρ1 + ρ2 then Eend = Iend = 0 so that ρ1 + ρ2 the equilibrium point is, in fact, a disease-free one. By considering also the stationary Vaccines 2021, 9, 398 11 of 26 recovered, one gets after cancelling the stationary recovered Rk+1 = Rk = Rend in both sides of (4) that the subsequent constraint should hold: γIend + K (ρ1 + ρ2 )Send = µEend + K (ρ1 + ρ2 )Send = 0 (42) which contradicts that Send > 0 and Eend > 0. As a result, no endemic equilibrium point exists with Send > 0. Since (42) can only hold with Eend = Send = 0 then, it only could eventually be true for the disease-free equilibrium point. Thus, Case 1 is unfeasible for the existence of an endemic equilibrium point. Case 2 (Send ≥ 0). Note that (40) holds with, which includes also Send > 0 of Case 1 if   µ a = β( Iend + λe Eend ) + K (ρ1 + ρ2 ) = β + λe Eend + K (ρ1 + ρ2 ) (43) γ provided that a − K ( ρ1 + ρ2 ) Eend = γ (44) β(µ + γλe ) provided that a K< ≤ 1; a ≤ ρ1 + ρ2 (45) ρ1 + ρ2 with Send = 0, and µ a − K ( ρ1 + ρ2 ) Iend = Eend = µ (46) γ β(µ + γλe ) However, again, after cancelling Rk+1 = Rk = Rend in both sides of (4), one gets that the subsequent constraint should hold: γIend + K (ρ1 + ρ2 )Send = µEend + K (ρ1 + ρ2 )Send = 0 (47) which implies that Send = Eend = Iend = 0 (which only holds if K = ρ +a ρ2 ), and it is in fact 1 a disease- free equilibrium point. In this case, one also has that Nend = Rend . Thus, Case 2 is unfeasible for the existence of an endemic equilibrium point. Remark 3 (basic reproduction number). In biological terms, it is possible to re-interpret the condition of asymptotic stability around the disease-free equilibrium point in terms of the basic reproduction number, which indicates the number of secondary infectious individuals generated from one primary infectious one, defined by β β ( µ + λ e γ ) Sd f R= = (48) βc µγ It turns out that R < 1 is exactly identical to β < β c , that is, if the transmission rate is smaller than its critical value β c then, equivalently, the basic reproduction number is smaller than unity. That asymptotic stability condition, in the two mentioned equivalent forms, are related to the asymptotic extinction of the disease. It can be commented that such a condition is not related to the transient disease evolution but to its steady-state which is the disease-free equilibrium in this case. We can also to define a sample-dependent effective reproduction number as follows: β µ + λek γ Sk  Rk = ; k ∈ Z0+ (49) µγ whose meaning and evolution should not be confused with that of the basic reproduction number. This number can be greater than one at the beginning of the disease transmission (even if R < 1) but Rk → R as k → ∞ . Note also from (1) that Rk decreases faster (since Sk decreases faster) if vaccination is programmed than in the vaccination-free case. The effective reproduction number is usually periodically checked to elucidate the particular intervention measures to be taken depending Vaccines 2021, 9, 398 12 of 26 on the disease transmission evolution. It can be also pointed out that R = 1 is related, in common situations associated with epidemic models, to the confluence of the disease-free equilibrium point with the endemic one while, for R < 1, the disease-free equilibrium point is unstable and the endemic equilibrium point is an attractor for the solutions. According to the former discussion in this subsection, this model, which is proposed for short-term predictions, does not evaluate the influence possible endemic steady states of the disease so that the central discussion is related to the case when R < 1. 5. Simulation Results This section contains some simulation examples devoted to study the application of the two-doses vaccination to counteract the spread of the COVID-19 pandemic. To this end the case of Italy borrowed from [45], where the discrete-model dynamics is confronted with actual data from Italy, will be considered. Thus, we have a discrete-time COVID-19 model serving as benchmark for the vaccination control. The model is parameterized by: β 1 β = 0.2, βe = , µ = , γ = 0.04 1.3 6 so that λe = 0.7692. The sampling time is one day so that the units of the parameters are in days−1 . The initial conditions are S0 = 0.9999, E0 = 0.0001, I0 = 0, R0 = 0 implying that the total population is normalized to unity, without loss of generality. Notice that almost all the population is susceptible and a small fraction of the population is exposed at the beginning. The parameters of the vaccination are d1 + d2 = 21 days and d2 = 7 days; therefore, separation between the two doses is of two weeks. The values of the doses effectiveness are given by ρ1 = 0.66 and ρ2 = 0.3 in such a way that the total effectiveness of the two doses ρ = ρ1 + ρ2 is 96%, in accordance with the average effectiveness of available vaccines. The natural recruitment rate is a = 1, since the natural growth of the population Vaccines 2021, 9, x FOR PEER REVIEW may be rejected when it comes to the epidemic spreading description due to the 14 ofsmall 29 number of children affected. The Figure 1 displays the dynamics of the model without vaccination. 1 0.9 0.8 S E 0.7 I R 0.6 Populations 0.5 0.4 0.3 0.2 0.1 0 0 50 100 150 time (days) Figure Figure 1. Dynamics 1. Dynamics of the of the model model for for COVID-19 COVID-19 pandemic. pandemic. It isItremarkable is remarkable in Figure in Figure 1 how 1 how thethe number number of infectious of infectious reaches reaches a peak a peak at around at around thethe 50%50% of the of the totaltotal population. population. NowNow vaccination vaccination is applied is applied with awith gain aofgain of K cor- K = 0.01, = 0.01, corresponding to the 1% of the susceptible vaccinated every day. The vaccination responding to the 1% of the susceptible vaccinated every day. The vaccination is applied is applied starting at different moments of the epidemic spreading in order to observe the effect of starting at different moments of the epidemic spreading in order to observe the effect of vaccination in the epidemic dynamics. Thus, Figure 2 displays the effect of vaccination vaccination in the epidemic dynamics. Thus, Figure 2 displays the effect of vaccination when it is applied since the beginning of the spreading. when it is applied since the beginning of the spreading. 1 0.9 S E It is remarkable in Figure 1 how the number of infectious reaches a peak at around the 50% of the total population. Now vaccination is applied with a gain of K = 0.01, cor- responding to the 1% of the susceptible vaccinated every day. The vaccination is applied starting at different moments of the epidemic spreading in order to observe the effect of Vaccines 2021, 9, 398 vaccination in the epidemic dynamics. Thus, Figure 2 displays the effect of vaccination 13 of 26 when it is applied since the beginning of the spreading. 1 0.9 S 0.8 E I R 0.7 0.6 Populations 0.5 0.4 0.3 0.2 0.1 0 0 50 100 150 time (days) Figure Figure 2. Dynamics 2. Dynamics of the of the model model whenwhen two-doses two-doses vaccination vaccination is applied is applied fromfrom the the beginning. beginning. It isItobserved is observed fromfrom Figures Figures 1 and 1 and 2 how 2 how thethe application application of vaccination of vaccination reduces reduces thethe peak of the infectious while increasing the pace at which susceptible peak of the infectious while increasing the pace at which susceptible become immune, as become immune, as expected. expected. In addition, In addition, the Figure the Figure 3 illustrates 3 illustrates theof the effect effect of starting starting the vaccination the vaccination at dif- at ferent moments of the spreading ranging from the beginning, 29 days and 59 days 59 different moments of the spreading ranging from the beginning, 29 days and days after the first day of simulation. It can be seen in Figure 3 how the effect of vaccination in in after the first day of simulation. It can be seen in Figure 3 how the effect of vaccination spreading is higher the sooner vaccination is applied when cases are detected. Thus, vacci- spreading is higher the sooner vaccination is applied when cases are detected. Thus, nation is useful for preventing new outbreaks and reducing their severity. If vaccination is vaccination is useful for preventing new outbreaks and reducing their severity. If vac- applied when there is a relatively large number of cases among the population, the peak cination is applied when there is a relatively large number of cases among the popula- is not reduced ostensibly. Figure 4 displays the values of vaccination actions V1 and V2 tion, the peak is not reduced ostensibly. Figure 4 displays the values of vaccination ac- corresponding to this situation. It is deduced from Figure 4 that the feedback vaccination tions V1 and V2 corresponding to this situation. It is deduced from Figure 4 that the provides a lower number of vaccinated individuals when it starts at advanced stages of the feedback spreadvaccination provides due to a smaller numbera lower number of of susceptible vaccinatedTherefore, individuals. individuals thiswhen kind it of starts action is especially recommendable when the spread is at its initial steps. Figure 4 also shows that Vaccines 2021, 9, x FOR PEER REVIEW vaccination V1 has larger values than V2 due to the differences in effectiveness 15 of 29 defined by the parameters ρ1 and ρ2 . Finally, Figure 5 displays the dynamics of the model for different values of the vacci- at nation advanced gain. As itofcould stages have been the spread due to expected, a smaller the larger of number thesusceptible vaccination action is, the smaller individuals. Therefore, this kind of peak is attained. action is especially Naturally, recommendable the reduction in the peakwhen the spread is achieved is at at the its ini- of a higher expense tialvaccination steps. Figure 4 also effort as shows Figurethat vaccination 6 shows. 1 hasbe It canValso larger valuesinthan observed V2 due Figures to the 4 and 6 the effect of differences delay inin theeffectiveness defined control action, by the where parameters a value of zero and  2 . 1 vaccination is provided during the first days of simulation due to the delay. 0.5 Start time = 0 days 0.45 Start time = 29 days Start time = 59 days 0.4 0.35 0.3 Infectious, I 0.25 0.2 0.15 0.1 0.05 0 0 50 100 150 time (days) Figure Figure 3. Effect 3. Effect of vaccination of vaccination starting starting time ontime on the infectious. the infectious. -3 10 -3 10 7 3 Start time = 0 days Start time = 0 days Start time = 29 days Start time = 29 days Start time = 59 days 6 Start time = 59 days 2.5 0.05 0 0 50 100 150 time (days) Vaccines 2021, 9, 398 14 of 26 Figure 3. Effect of vaccination starting time on the infectious. -3 10 -3 10 7 3 Start time = 0 days Start time = 0 days Start time = 29 days Start time = 29 days Start time = 59 days 6 Start time = 59 days 2.5 5 2 Vaccination, V 2 Vaccination, V 1 4 1.5 3 1 2 0.5 1 0 Vaccines 2021, 09,0 x FOR PEER REVIEW 50 100 150 0 50 100 16 150 of 29 time (days) Vaccines 2021, 9, x FOR PEER REVIEW time (days) 16 of 29 Figure 4. Vaccination actions corresponding to different starting times. Figure 4. Vaccination actions corresponding to different starting times. 0.35 0.35 K = 0.005 Finally, Figure 5 displays the dynamics Kof the 0.010 model for different values of the vac- K= =0.005 0.3 K K= =0.010 0.015 cination 0.3 gain. As it could have been expected, K = 0.015the larger the vaccination action is, the smaller 0.25 peak is attained. Naturally, the reduction in the peak is achieved at the expense of0.25 a higher vaccination effort as Figure 6 shows. It can also be observed in Figures 4 and 6 the0.2 effect of delay in the control action, where a value of zero vaccination is provided Infectious, I Infectious, I 0.2 during the first days of simulation due to the delay. 0.15 0.15 0.1 0.1 0.05 0.05 0 0 50 100 150 0 0 50 time (days) 100 150 time (days) Figure 5. Dynamics of the infectious for different values of the vaccination gain. Figure Figure 5. Dynamics 5. Dynamics of the of the infectious infectious forfor different different values values of the of the vaccination vaccination gain. gain. -3 10 4.5 0.01 -3 10 K = 0.005 K = 0.005 4.5 0.01 4 KK == 0.010 0.005 0.009 KK = 0.010 = 0.005 KK == 0.015 0.010 0.009 KK = 0.015 = 0.010 4 K = 0.015 0.008 K = 0.015 3.5 0.008 3.5 0.007 3 0.007 Vaccination, V 1 Vaccination, V 2 3 0.006 Vaccination, V 1 Vaccination, V 2 2.5 0.006 0.005 2.5 0.005 2 0.004 2 0.004 1.5 0.003 1.5 0.003 1 0.002 1 0.002 0.001 0.5 0.001 0.5 0 0 0 50 100 150 0 0 50 100 150 0 50 time (days) 100 150 0 0 50 time (days) 100 150 time (days) time (days) Figure 6. Vaccination actions for different values of the control gain K. Figure 6. Vaccination actions for different values of the control gain K. Figure 7 shows Figure 6. Vaccination the actions foreffect of changing different thecontrol values of the interval of K. gain time between the application of the Figure two doses. It is assumed 7 shows the effectthat the effectivity of changing does not the interval ofchange with thethe time between period in between application of Figure them. The 7 shows the vaccination effect gain K of=changing 0.01 is the used interval while theof time initial between conditions the are the two doses. It is assumed that the effectivity does not change with the period in be- application those of indicated the two them. tween doses.TheIt isvaccination assumed that gainthe K effectivity doeswhile = 0.01 is used not change with the initial the period conditions areinthose be- tween them. The vaccination gain K = 0.01 is used while the initial conditions indicated at the beginning of this section. Moreover, d2 = 7 days and d1 + d2 range from 10 are those indicated at It to 21 days. theis beginning observed in of Figure this section. 7 that Moreover, if the secondd2 =dose 7 days and d1 +effectivity maintains d2 range from 10 regard- toless 21 the days. It is observed in Figure 7 that if the second dose maintains effectivity dose sparing, it is better to administer it as soon as possible. However, in the real regard- less the dose sparing, it is better to administer it as soon as possible. However, in the real Vaccines 2021, 9, 398 15 of 26 at the beginning of this section. Moreover, d2 = 7 days and d1 + d2 range from 10 to 21 days. Vaccines 2021, 9, x FOR PEER REVIEW It is observed in Figure 7 that if the second dose maintains effectivity regardless the dose 17 of 29 sparing, it is better to administer it as soon as possible. However, in the real world, the Vaccines 2021, 9, x FOR PEER REVIEW effectivity of the second dose may depend on when this is applied. There are no17 of 29on data how the effectivity depends on the dose sparing or these are scarce, [46]. Therefore a more accurate 0.25 simulation on the effect of dose sparing is not carried out in this work. d 1 +d 2 = 10 days 0.25 d 1 +d 2 = 14 days dd1 +d +d2 == 10 17 days days 0.2 1 2 dd1 +d +d == 14 21 days days 1 22 d +d = 17 days 0.2 1 2 d +d = 21 days 1 2 0.15 Infectious, I 0.15 Infectious, I 0.1 0.1 0.05 0.05 0 0 50 100 150 200 250 time (days) 0 0 50 100 150 200 250 Figure 7. Dynamics of the infectious for different time periods between the two doses. time (days) Figure Figure 7. Dynamics 7. Dynamics of the of the infectious infectious forfor different different time time periods periods between between thethe twotwo doses. doses. Moreover, Figure 8 compares the evolution of the infectious when two doses with a spacingMoreover, of 14 daysFigure 8 compares is applied (K = 0.1)theinevolution contrast ofof applying the infectious when a single twowith dose doses with a a dou- Moreover, Figure 8 compares the evolution of the infectious when two doses with a spacing of 14 days ble-vaccination is applied gain (K = 2x0.01(K = 0.1) The = 0.02). in contrast Figure of applying 8 shows thata with singlerespect dose with a double- to the evo- spacing of 14 days is applied (K = 0.1) in contrast of applying a single dose with a dou- vaccination lution of the infectious × 0.01 gain (K = 2the = 0.02). Theof administration Figure 8 shows a single dosethat to awith respect broader to the evolution population alle- ble-vaccination gain (K = 2x0.01 = 0.02). The Figure 8 shows that with respect to the evo- of the viates infectious the the administration peak of infectious of a single in comparison dose to aapplying to promptly broader population the second alleviates dose in or-the lution of the infectious the administration of a single dose to a broader population alle- derpeak of infectious to increase in comparison the overall to promptly effectiveness of the applying the second vaccine. This behavior dose is in orderthe behind to increase deci- viates the peak of infectious in comparison to promptly applying the second dose in or- theofoverall sion the UK effectiveness governmentofofthe vaccine. This prioritizing behavior the first dose is to behind as many thepeople decision of the UK as possible dergovernment to increase theprioritizing overall effectiveness of to theasvaccine. This behavior is behind the deci-the while delayingof the second dose, the [46]. first dose many people as possible while delaying sion of the UK government second dose, [46]. of prioritizing the first dose to as many people as possible while delaying the second dose, [46]. 0.25 double-dose single-dose 0.25 double-dose 0.2 single-dose 0.2 0.15 Infectious, I 0.15 Infectious, I 0.1 0.1 0.05 0.05 0 0 50 100 150 200 250 time (days) 0 0 50 100 150 200 250 Figure Figure 8. Comparison 8. Comparison in the in the evolution evolution time of the (days) of the number number of infectious of infectious between between thethe administration administration of aof a double dose (K = 0.01) or a single dose to a broader population (K double dose (K = 0.01) or a single dose to a broader population (K = 0.02). = 0.02). Figure 8. Comparison in the evolution of the number of infectious between the administration of a double dose (K = 0.01) The basic or a single properties of thedose to a broader model population discussed (K = 30.02). in Sections and 4 can also be corroborated The basic properties of the model discussed in Sections 3 and 4 can also be corrobo- in the rated in previous examples. the previous Thus,Thus, examples. the number of exposed the number and infectious of exposed always converges and infectious always The basic properties of converges the model to discussed in Sections 3 and 4 can also be corrobo- converges to zero so that the model converges to the disease-free equilibriumThis to zero so that the model the disease-free equilibrium point. happens point. This rated in the previous examples. Thus, the number ofdiscussed exposed in and infectious always happens because there is no endemic point for this model as discussed in Section Figure because there is no endemic point for this model as Section 4. The 4. The 9 converges to zero so that the model converges to the disease-free equilibrium point. This Figure 9 displays a longer simulation for the model with K = 0.001 and different initial happens because there is no endemic point for this model as discussed in Section 4. The conditions. It is shown in Figure 9 that the model always converges to the globally-stable Figure 9 displays a longer simulation for the model with K = 0.001 and different initial disease-free equilibrium point given by S  Edf  I df  0, Rdf  1 and all the subpopula- conditions. It is shown in Figure 9 that the dfmodel always converges to the globally-stable Vaccines 2021, 9, x FOR PEER REVIEW 18 of 29 Vaccines 2021, 9, 398 16 of 26   rems 1 and 2 since  ,   0,1 and   0.2  1.66  max  ak  is small enough; there-  I  e E  displays a longer simulation for the model with K = 0.001 k k k  and different initial conditions. fore, (15) holds, and the model trajectory solution is non-negative It is shown in Figure 9 that the model always converges to the globally-stable according todisease-free Theorem 1 equilibrium and converges to the disease-free equilibrium point according to point given by Sd f = Ed f = Id f = 0, Rd f = 1 and all the subpopulations Theorem 2 since a  0.001 non-negative Kremain  1.04  at. all In time. These addition, properties Theorem 3 alsoare indeed ensures guaranteed that by Theorems the disease-free equi- 1    a and 2 since µ, γ ∈1 [0, 21) and β = 0.2 < 1.66 = max I +λe E is small enough; therefore, k k k k librium point is globally stable since Condition 1 from Theorem 3 is satisfied, Condition (15) holds, and the model trajectory solution is non-negative according to Theorem 1 and a 2 converges also holdstoasthe  0.001  1.04equilibrium K disease-free  point and according Equation (13)toholds Theorem K= 2 since 10 as Figure 0.001 6= depicts. 1.04 = ρ +a ρ2 . In addition, Theorem 1 3also 2 ensures that the disease-free equilibrium point 1 Therefore, is globally the disease-free stable equilibrium since Condition 1 frompoint is concluded Theorem to be globally 3 is satisfied, Condition asymptotically 2 also holds as stable. a K = 0.001 6= 1.04 = ρ +ρ2 and Equation (13) holds as Figure 10 depicts. Therefore, the 1 disease-free equilibrium point is concluded to be globally asymptotically stable. 1 1 S 0.9 E 0.9 S I E R I 0.8 0.8 R 0.7 0.7 0.6 0.6 Populations Populations 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 50 100 150 200 250 0 50 100 150 200 250 time (days) time (days) 1 1 S S 0.9 E 0.9 E I I 0.8 R R 0.8 0.7 0.7 0.6 0.6 Populations Populations 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 50 100 150 200 250 0 50 100 150 200 250 time (days) time (days) Figure 9. Dynamics of the model for different initial conditions and K = 0.001. Figure 9. Dynamics of the model for different initial conditions and K = 0.001. Vaccines 2021, Vaccines 9, x9,FOR 2021, 398 PEER REVIEW 19 of 17 29 of 26 0.005 0 V +V 1 2 -0.005 e (a - (I + E )-1)S k k k k k k -0.01 -0.015 -0.02 -0.025 -0.03 -0.035 0 50 100 150 200 250 time (days) Figure Figure 10.10. Satisfaction Satisfaction of Condition of Condition 3 from 3 from Theorem Theorem 3. 3. It can It can be be pointed pointed out out that that differentalternative different alternativetheoretic theoreticmethods methodsfor forstability stabilityanaly- analysis of dynamic systems in the presence of delays can be found, for instance, sis of dynamic systems in the presence of delays can be found, for instance, in [47–49]in [47–49] and some of the references therein. and some of the references therein. 6. Concluding Remarks 6. Concluding Remarks A new discrete SEIR model has been presented in this paper subject to two delayed A new discrete SEIR model has been presented in this paper subject to two delayed doses of feedback vaccination controls on the susceptible. It is also assumed that there exists doses of feedback vaccination controls on the susceptible. It is also assumed that there a transmission rate from the exposed to the susceptible and that the transmission rates exists a transmission rate from the exposed to the susceptible and that the transmission exposed-susceptible and infectious-susceptible may eventually be distinct. The second rates ideaexposed-susceptible is based on the knowledgeand infectious-susceptible that the later period may in eventually the asymptomaticbe distinct. The incubation second idea is based on the knowledge that the later period in the asymptomatic phase of the infection in the COVID-19 pandemic is also infective as it is the first phase incuba- tion of phase of the infection the infectious period. in thefirst The COVID-19 pandemic idea of using two is also infective vaccination doses as isit based is the onfirstthe phase administration mechanisms of some of the recently approved vaccines for COVID-19. on of the infectious period. The first idea of using two vaccination doses is based Such thetwo administration doses are applied mechanisms of some with a certain of theand interval, recently approved the immunity vaccines effect for COVID- is delayed in respect 19.toSuch eachtwo doses current are applied picked- withon up sample a the certain interval, model, while and the immunity it is also assumed that effect bothis doses de- layed in respect to each current picked- up sample on the model, while can have, in general, different effectiveness. It is assumed, in particular, that the second it is also as- sumed that both doses can have, in general, different effectiveness. It is assumed, dose has an incremental benefit on the injection of the first single dose. The effectiveness in par- ticular, degreesthatofthe second both dosesdose hasas as well anthe incremental delays canbenefit on theininjection be adjusted of the firsttests the experimental single since dose. they Theareeffectiveness degrees of model parameters. Theboth doses as wellproperties non-negativity as the delays cansolution of the be adjusted in the under finite experimental non-negative tests sinceconditions, initial they are model parameters. and its boundedness The the non-negativity disease-freeproperties equilibrium of the point solution as wellunder finite non-negative as its stability properties areinitial alsoconditions, investigated.itsThe boundedness non-negativity and the of thedisease- solution free canequilibrium become lost point as well as itstransmission for large-enough stability properties areitalso rates, but investigated. is kept for low and The non- moderate negativity of the solution can become lost for large-enough transmission values of the transmission rates. The proposed model has been numerically tested through rates, but it is kept for low examples numerical and moderate values ofparameterizations. for COVID-19 the transmission rates. The proposed model has been numerically tested through numerical examples for COVID-19 parameterizations. Author Contributions: Conceptualization, M.D.l.S., R.N. and A.I.; methodology, M.D.l.S. and A.I.; Author Contributions: software, Conceptualization, A.I.; validation, S.A.-Q. and A.I.; M.D.L.S., R.N. and formal analysis, A.I.; investigation, M.D.l.S.; methodology,M.D.l.S., M.D.L.S.A.I. andand A.I.; software, S.A.-Q.; A.I.; validation, resources, M.D.l.S. andS.A-Q. A.I.; and dataA.I.; formal curation, analysis, A.I.; M.D.L.S.; investigation, writing—original M.D.L.S., draft preparation, M.D.l.S. A.I.and andA.I.; S.A-Q.; resources, M.D.L.S. writing—review and A.I and editing, .; data curation, M.D.l.S., A.I.; writing—original A.I. and S.A.-Q.; draftsupervision, visualization, A.I.; prepara- tion, M.D.L.S. S.A.-Q., and R.N.; A.I.; writing—review project and editing, administration, M.D.l.S.; M.D.L.S., funding A.I. andM.D.l.S. acquisition, S.A-Q.;All visualization, authors haveA.I.; read supervision, and agreedS.A-Q., R.N.; project to the published administration, version M.D.L.S..; funding acquisition, M.D.L.S. All au- of the manuscript. thors have read and agreed to the published version of the manuscript. Funding: This research was funded by MCIU/AEI/FEDER, UE, grant number RTI2018-094336-B-I00; Funding: SpanishThis research Institute was funded of Health by grant Carlos III, MCIU/AEI/FEDER, UE, grantand number COV 20/01213 number BasqueRTI2018-094336-B- Government, grant I00;number SpanishIT1207-19. Institute of TheHealth APC wasCarlos III, grant funded number COV 20/01213 by MCIU/AEI/FEDER, UE. and Basque Government, grant number IT1207-19 .The APC was funded by MCIU/AEI/FEDER, UE. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Vaccines 2021, 9, 398 18 of 26 Data Availability Statement: All data used are in the references and properly cited within the manuscript. Acknowledgments: The authors are grateful to the Spanish Government for Grant RTI2018-094336- B-I00 (MCIU/AEI/FEDER, UE) and to the Basque Government for Grant IT1207-19. They also thank the Spanish Institute of Health Carlos III for its support through Grant COV 20/01213. Conflicts of Interest: The authors declare that they do not have any competing interests. Appendix A Appendix A.1 Some Auxiliary Technical Results Theorem A1. Consider the difference equation:  xk+1 = axk − K ρ2 xk−d2 + ρ1 xk−d1 −d2 + uk ; k ∈ Z0+ (A1) with ρ1 , ρ2 ∈ R+ , a ∈ R, K ∈ R, nonnegative(integer) delays d1 > 0 and d1 + d2 > d1 , and initial conditions x−d1 −d2 = x−d2 = 0 and x0 ∈ R and eventually subject to a bounded forcing sequence {uk }∞ k =0 ⊂ R. Then, the following properties hold: (i) The unforced difference Equation (A1), i.e., if {uk }∞k =0 ≡ 0, is globally asymptotically ∞ ∞ stable,   so that { x } k k =0 is bounded  and { x }  k k=0 → 0 for any given finite x0 ∈ R, if K ∈ a −1 + a , K̄ max −K̄ , min ρ11 + ρ1 + ρ2 , ρ2 and, in particular, K ∈ [0, K̄ ) under the restriction K ∈ R0+ = R+ ∪ {0}, where: a−1     1+a K̄ = sup y ∈ , : y < 1/ f (y) with f : R0+ → R0+ de f ined by : ρ1 + ρ2 ρ1 + ρ2 v u (ρ2 (cos(θd2 ) − 1) + ρ1 (cos(θ (d1 + d2 )) − 1))2 + (ρ2 sin(θd2 ) + ρ1 sin(θ (d1 + d2 )))2 u f (y) = sup t θ ∈(0,2π ) 1 + ( a − (ρ1 + ρ2 )y)2 − 2( a − (ρ1 + ρ2 )y)cosθ (ii) The unforced difference Equation (A1) is globally asymptotically stable for any given finite x0 ∈ R, so that { xk }∞ ∞ k =0 is bounded and { xk }k =0 → 0 for any finite x0 ∈ R, if | a | < 1 and ρ2 e−iθd2 +ρ1 e−iθ (d1 +d2 ) K ∈ (−K̄, K̄ ), were K̄ = 1/ sup eiθ − a . θ ∈[0,2π ) (iii) Assume that either the conditions of Property (i) or those of Property (ii) hold and that {uk }∞ ∞ ∞ k =0 → 0 . Then, { xk }k =0 is bounded and { xk }k =0 → 0 for any given finite initial condition x0 . (iv) Assume that either the conditions of Property (i) or those of Property (ii) hold and that {uk }∞ ∞ k =0 is bounded. Then, { xk }k =0 is bounded for any given finite initial conditions. Proof. By using the one-step-ahead and the one-step-backward operators q and q−1 defined, respectively, by qxk = xk+1 and q−1 xk+1 = xk , one can compactly rewrite (A1) as follows by defining K1 = Kρ1 and K2 = Kρ2 :   q − a + K2 q−d2 + K1 q−(d1 +d2 ) xk = uk (A2) The characteristic equation of the unforced system (A1) is derived from the left-hand- side of (A2) by taking into account the formal analogy between the z-transform operator and the one-step-ahead one and also multiplying the result by the factor zd1 +d2 to avoid Vaccines 2021, 9, 398 19 of 26 the concourse of mixed contributions of the forward and backward operators. As a result, the characteristic equation of (A2) becomes:   p0 (z)= zd1 +d2 z − a + K2 z−d2 + K1 z−(d1 +d2 )      (A3) = zd1 +d2 z − a + K1 + K2 + K2 z−d2 − 1 + K1 z−(d1 +d2 ) − 1 = 0 The polynomial p0 (z) is stable if its zeros are within |z| < 1 and the circle is centered at the origin of unity radius. In view of (A.3), since zd1 +d2 has a zero of multiplicity d1 + d2 in the center z = 0 of the open circle |z| < 1, the stability condition holds if     p(z) = z − a + K1 + K2 + K2 z−d2 − 1 + K1 z−(d1 +d2 ) − 1 (A4) has all its zeros in the open complex circle |z| < 1. Note that ρi = Ki /K (i = 1, 2) so that one gets from (A4) that: p(z) = p0 (z) + pe(z) (A5)      p0 (z) = z − a + (ρ1 + ρ2 )K ; pe(z) = K ρ2 z−d2 − 1 + ρ1 z−(d1 +d2 ) − 1 (A6) is stable if −1 < a − (ρ1 + ρ2 )K < 1 and if, in addition, | pe(z)| < | p0 (z)| for any z ∈ C on the circumference |z| = 1, i.e., on the boundary of the circle |z| ≤ 1: This follows from the Rouché0 s theorem of zeros [40–43], which establishes that if all the zeros of p0 (z) lie in |z| < 1 (i.e., they are stable) and | pe(z)| < | p0 (z)| on |z| = 1, which is trivially a Jordan curve, then all the zeros of q(z) = p0 (z) + pe(z) are also in |z| < 1 since q(z) has the same number of zeros in |z| < 1 as p0 (z). The condition | pe(z)| < | p0 (z)| on |z| = 1 is identical to the following implicit constraint:     ρ2 e−iθd2 − 1 + ρ1 e−iθ (d1 +d2 ) − 1 |K |−1 > sup (A7) θ ∈[0,2π ) eiθ − a + (ρ1 + ρ2 )K after writing any complex number z with |z| = 1 as z = eiθ for some θ ∈ [0, 2π ). By inspecting (A7), we note that, for θ = 0, pe(0) = 0 so that θ = 0 may be removed from the test (A7) so that by using cos(−θ ) = cosθ and sin(−θ ) = −sinθ; ∀θ ∈ [0, 2π ), we get: v u (ρ2 (cos(θd2 ) − 1) + ρ1 (cos(θ (d1 + d2 )) − 1))2 + (ρ2 sin(θd2 ) + ρ1 sin(θ (d1 + d2 )))2 u 1 > |K | sup t (A8) θ ∈(0,2π ) (cosθ − a + (ρ1 + ρ2 )K )2 + sin2 θ Then, Property (i) follows from (A8) combined with −1 < a − (ρ1 + ρ2 )K < 1, which jointly ensure that the unforced difference Equation (A1) has the property that { xk }∞ k =0 → 0 for any finite x0 by defining f : R0+ → R0+ by v u (ρ2 (cos(θd2 ) − 1) + ρ1 (cos(θ (d1 + d2 )) − 1))2 + (ρ2 sin(θd2 ) + ρ1 sin(θ (d1 + d2 )))2 u f (y) = sup t (A9) θ ∈(0,2π ) 1 + ( a − (ρ1 + ρ2 )y)2 − 2( a − (ρ1 + ρ2 )y)cosθ and by noting that sup{y ∈ R0+ : y < (1/ f (y))} 6= ∅. The vaccination gain constraint is simplified  to K ∈ [0, K̄) in the  case that K ∈ R0+ since the general constraint  particular a −1 1+ a − a , that is, K̄ > a−1 . K ∈ max ρ1 + ρ2 , −K̄ , min implies that −K̄ < ρ1+ ρ1 + ρ2 , K̄ ρ1 + ρ2 n o 1 ρ2 a +1 a +1 Since, furthermore, K̄ < ρ +ρ2 , it follows that min ρ +ρ2 , K̄ = K̄ and if K ∈ R0+ then the 1 1 vaccination gain constraint is simplified to K ∈ [0, K̄ ) . Moreover, { xk }∞ k =0 is bounded since any convergent sequence is bounded. Since { xk }∞ k =0 is bounded and { xk }∞ k =0 → 0 for any finite x0 then the unforced difference Equation (A1) is globally asymptotically stable. Vaccines 2021, 9, 398 20 of 26 To prove Property (ii), rewrite (A5) as p(z) = p01 (z) + pe1 (z) where   p01 (z) = z − a; pe1 (z) = K ρ2 z−d2 + ρ1 z−(d1 +d2 ) (A10) By using the same arguments that those in the proof of Property (i), it is found that the zeros of p(z) are in |z| < 1 if −1 < a < 1 and ρ2 e−iθd2 + ρ1 e−iθ (d1 +d2 ) | K | −1 > sup (A11) θ ∈(0,2π ) eiθ − a after the replacement of (A7) by (A11). As a result, the unforced difference Equation (A1) is globally asymptotically stable. To prove, Property (iii), we first obtain from (A1) the following extended dynamic system of dimension d1 + d2 + 1: x̂k+1 =  x̂k + ûk (A12) subject to initial conditions x−i = 0 for i = 1, 2, . . . , d1 + d2 , where T x̂k = xk , xk−1 , . . . , xk−d2 , . . . , xk−d1 −d2 ; ûk = (uk , 0, . . . , 0, . . . , 0) T (A13)   a 0...0 −Kρ2 0 . . . 0 −Kρ1   1 0 0 0 0    =  .. .. .. .. ..  (A14)   . . . . .   0 0 0 1 0 First note that all the eigenvalues of  are stable under the conditions of either Property (i) or Property (ii) since if { xk }∞ ∞ ∞ ∞ k =0 → 0 if { uk }k =0 ≡ 0 then { x̂k }k =0 → 0 if { ûk }k =0 ≡ 0 so that the unforced discrete dynamic system (A12) is globally asymptotically stable. Now, from recursive calculation using (A12), one gets: k + N −1 x̂k+ N =  N x̂k + ∑ Âk+ N −i−1 ûi ; k, N ∈ Z0+ = Z+ ∪ {0} (A15) i =k If {uk }∞k =0 → 0 then, for any given ε ∈ R+ , there exists N = N ( ε ) such that k ûi k = |ui | ≤ ε; ∀i (∈ Z0+ ) ≥ N. Then, for some α̂ ∈ (0, 1) such that the complex circle centered at z = 0 of radius α̂ contains all the eigenvalues of Â, one gets: ∞ εk k x̂k+ Nk k ≤ Mα̂ Nk k x̂k k + ∑i=k Mα̂k+ Nk −i−1 kûi k ≤ Mα̂ Nk k x̂k k + M (A16) 1 − α̂ for some (norm-dependent) real constant M ≥ 1, for any given k ∈ Z0+ and any given strictly decreasing sequence {ε k }∞ k =0 ⊂ R0+ , with ε 0 = ε, and some corresponding strictly increasing sequence of integers { Nk }∞ ∞ k =0 such that Nk = Nk ( ε k ). Since { ε k }k =0 → 0 and α̂ Nk → 0 as k → ∞ then there exists the limit lim x̂k+ Nk = 0, so that lim xk+ Nk = 0 as k→∞ k→∞ well, if x0 is finite, and then { xk }∞ k =0 → 0 and, since this sequence is convergent, it is also bounded. Property (iii) has been proved. Property (iv) follows by noting that the first inequality of (A16) partially holds if {uk }∞ k =0 is bounded in the weakened form: ∞ Mε C k x̂k k ≤ Mα̂k k x̂0 k + ∑i=0 Mα̂k−i−1 kûi k ≤ Mα̂k k x̂0 k + ≤ Mk x̂0 k + < +∞; ∀k ∈ Z0+ (A17) 1 − α̂ 1 − α̂ Mε lim supk x̂k k ≤ < +∞ (A18) k→∞ 1 − α̂ Vaccines 2021, 9, 398 21 of 26 for any given k ∈ Z0+ and some real constant +∞ > C ≥ supkûi k = sup|ûi |.  i ≥0 i ≥0 Corollary A1. Consider the difference equation:  xk+1 = ak xk − Kk ρ2 xk−d2 + ρ1 xk−d1 −d2 ; k ∈ Z0+ (A19) subject to initial conditions x0 ∈ R and x−i = 0 for i = 1, 2, . . . , d1 + d2 and assume that the limiting difference equation:  xk+1 = axk − K ρ2 xk−d2 + ρ1 xk−d1 −d2 (A20) resulting from (A19) as { ak }∞ ∞ k =0 → a and { Kk }k =0 → K satisfies the conditions of the unforced Equation of system (A1) of Theorem A1 (i) and Theorem A1 (ii). Assume also that     |Kk − K | = o k x̂k k−1 ; | ak − a| = o k x̂k k−1 (A21) Then, { xk }∞ ∞ k =0 is bounded and { xk }k=0 → 0 for any given finite initial conditions so that the difference Equation (A19) is globally asymptotically stable. Proof. Equation (A19) is of the form of Equation (A1) with  u k = ( a k − a ) x k − ( K k − K ) ρ 2 x k − d2 + ρ 1 x k − d1 − d2 (A22) and {uk }∞ k=0 → 0 under (A21). Thus, the proof is direct from Theorem A1 (iii).  Corollary A2. Consider the difference Equation (A19) with x0 ≥ 0 and { ak }∞ k =0 ⊂ [0, 1], define αk = 1 − ak and assume that {αk }∞ ∞ ∞ µk k =0 → 0 , ∑ α k =0 k = + and K k = K − ρ2 x k − d2 + ρ1 x k − d − d2 ; 1 ∀k ∈ Z0+ , with {µk }∞ ρ2 xk−d2 + ρ1 xk−d1 −d2 K , +∞ being subject to the subsequent    k =0 ⊂ constraints: ∞ µ − ρ2 xk−d2 + ρ1 xk−d1 −d2 K < +∞; µk = 0 if xk−d2 xk−d1 −d2 = 0, and µk =  ∑ k =0 k   . Then, { xk }∞ ∞  −1  o ρ 2 x k − d2 + ρ 1 x k − d1 − d2 k =0 is bounded and { xk }k=0 → 0 for any given finite nonnegative initial conditions so that the difference Equation (A19) is globally asymptotically stable. Proof. The difference Equation (A19) is written equivalently as:  x k +1 = (1 − α k ) x k + µ k − K ρ 2 x k − d2 + ρ 1 x k − d1 − d2 (A23) with αk = 1 − ak ; ∀k ∈ Z0+ , { ak }∞ ∞ k =0 → 1 since { αk }k =0 → 0 and αk ≥ 0, µk ≥ 0; ∀k ∈ Z0+ , ∞ ∞  −1  ∑k=0 αk = +∞ and ∑k=0 µk < +∞ since µk = o ρ2 xk−d2 + ρ1 xk−d1 −d2 . Thus, it follows from Venter theorem [44] that { xk }∞ k =0 → 0 , and the sequence boundedness follows as a consequence of its convergence.  Theorem A2. Consider the difference equation: x k + 1 = a k x k + bk (A24) with x0 ≥ 0, where the following constraints hold: ak = 1 + e bk ; b > 0; ∀k ∈ Z0+ a k , bk = b + e (A25)  ∞ n o∞ a j j=0 ⊂ [0, e e a] ⊂ [0, +∞), e bj ⊂ [0, +∞) (A26) j =0 ak ∞ ; ∀k ∈ Z0+ , ∑k=0 ebk < +∞ e a k +1 < e (A27) ak + bk /xk Vaccines 2021, 9, 398 22 of 26 Thus, the following properties hold: (i) If x0 6= 0 then { xk }∞k =0 is unbounded. If b = 0 and the remaining above conditions hold, then { xk }∞ k =0 is bounded.  ∞ n o∞ |e a | (ii) If x0 6= 0, e a j j=0 ⊂ [−1, 0], |eak+1 | > a +bk /x , b > 0, e bj ⊂ [0, +∞) and k k k j =0 ∑∞ k = 0 bk e < +∞ then a necessary condition for { xk }∞ k =0 → 0 is ∑kj=0 e a j = −∞. The result still holds if { xk }∞ k =0 → x (> 0) . Proof. One gets recursively from (A24) and (A25) that   k x k +1 − x 0 − ( k + 1 ) b = ∑ j =0 eaj xj + e bj (A28) From the first inequality of (A27), (A24) and D0 Alembert criterion for convergence of series on non-negative terms, one gets: ak xk a x k < k k x k +1 = e ak xk ⇒ ∑ j=0 ea j x j < +∞; ∀k ∈ Z0+ e e a k +1 x k +1 < x (A29) a k x k + bk k +1 e x k +1 ak xk → 0 as k → ∞ . Now, (A29) together with the second condition of (A27) and so that e (A28) leads to: 0 ≤ xk+1 − x0 − (k + 1)b < +∞; ∀k ∈ Z0+ (A30) the non-negative lower-bound in (A30) arising from the non-negativity of the right-hand- side of the equality (A28). The relations (A30) lead to 0 ≤ lim sup( xk+1 − x0 − (k + 1)b) < +∞ (A31) k→∞ and also (A24) subject to (A26) and x0 ≥ 0 implies that { xk }∞ k =0 ⊂ clR0+ = R0+ ∪ {+ ∞ }. This concludes that +∞ ≤ lim sup xk = +∞ from (A31) and b > 0, which implies that there k→∞ exists the limit lim xk = +∞ (otherwise, lim sup xk would be finite) and then the sequence k→∞ k→∞ { xk }∞ k =0 is unbounded. If b = 0 and the remaining constraints (A25)–(A27) hold then (A31) is modified to 0 ≤ lim sup( xk+1 − x0 ) < +∞ implying that x0 ≤ lim sup( xk+1 ) < +∞ k→∞ k→∞ which concludes that { xk }∞k =0 is bounded. Property (i) has been proved. To prove Property (ii), note that under the alternative given conditions, if { xk }∞ k =0 → 0 , then one gets from (A28), by taking into account that ∑∞ b k =0 k e = Cb < + ∞ and D0 Alembert’s criterion for divergence of series, that ∞ ∞ ∞   −∞ = − lim ( x0 + kb) = ∑k=0 e ak xk + ebk ≤ Cb + ∑k=0 ea k x k ⇒ ∑ k =0 |e ak | xk = +∞ k→∞ (A32) ∞ and also (∑∞ k =0 | a e k |) sup x k = + ∞. Since { x k } k =0 is still non-negative and convergent, then k ≥0 it is bounded, and 0 < sup xk < +∞ if ∑∞ ak | = +∞ and the result is easily seen to still k =0 |e k ≥0 hold if { xk }∞ k =0 → x (> 0) . The proof of Property (ii) is complete.  Appendix A.2 A More Detailed Expansion of the Squared Numerator of (A9) Note that Equation (A9) is of the form: f (y) = sup g(θ, d1 , d2 , ρ1 , ρ2 , y) (A33) θ ∈(0,2π ) Vaccines 2021, 9, 398 23 of 26 for each given quadruple (d1 , d2 , ρ1 , ρ2 ), where: s n(θ, d1 , d2 , ρ1 , ρ2 ) g(θ, d1 , d2 , ρ1 , ρ2 , y) = (A34) d(θ, ρ1 , ρ2 , y) whose squared right-hand-side numerator n(θ, d1 , d2 , ρ1 , ρ2 ) is expanded as follows: n(θ, d1 , d2 , ρ1 , ρ2 ) 2 2 = (ρ2 (cos(θd 2 ) − 1) + ρ1 cos( θ ( d1 + d2 )) − 1) + ( ρ2 sin( θd2 ) + ρ1 sin( θ ( d1 + d2 ))) = 2 ρ22 + ρ21 − 2ρ22 cos(θd2 ) − 2ρ21 cos(θ (d1 + d2 )) + 2ρ1 ρ2 (1 − cos(θd2 ) − cos(θ (d1 + d2 )))  (A35) +2ρ1 ρ2 [sin(θd2 )(sin(θd2 )cos(θd1 ) + cos(θd2 )sin(θd1 )) + cos(θd2 )(cos(θd2 )cos(θd1 ) − sin(θd2 )sin(θd1 ))] = 2 ρ22 (1 − cos(θd2 )) + ρ21 (1 − cos(θ (d1 + d2 ))) + 2ρ1 ρ2 (1 + cos(θd1 ) − cos(θd2 ) − cos(θ (d1 + d2 ))) which leads to the following nested particular cases: (a) If d1 = d2 then, one gets by using 1 − cos(2θd2 ) = 2sin2 θd2 that:   n(θ, d1 , d1 , ρ1 , ρ2 ) = 2 ρ22 (1 − cos(θd2 )) + 2ρ1 (ρ1 + ρ2 )sin2 (θd1 ) (b) If d1 = d2 and ρ1 = ρ2 then   n(θ, d1 , d1 , ρ1 , ρ1 ) = 2 ρ21 (1 − cos(θd2 )) + 4ρ21 sin2 (θd2 ) (c) If d1 = d2 = 0 or ρ1 = ρ2 = 0 then n(θ, 0, 0, ρ1 , ρ1 ) = n(θ, d1 , d1 , 0, 0) = 0 (d) If θ = 0 then n(0, d1 , d1 , ρ1 , ρ1 ) = 0 which is not evaluated in the supremum over (0, 2π ) in (A.33). Note that the constraint y < 1/ f (y) to calculate K in Theorem A1 always holds since 1/ f (y) = +∞ for θ = 0, 2π since n(0, d1 , d1 , ρ1 , ρ1 ) = 0 so that it has not to be accounted for in the supremum evaluation as Theorem A1 formally establishes. Appendix A.3 A Fast and Simple Delay-Dependent Stability Test Based on Rouché0 s Theorem of Zeros Within Open Circles Contained in the Unit Circle Theorem A3. The following properties hold: (i) Let a ∈ (0, 1) and consider a circle |z| ≤ σ of center z = 0 and of radius σ ∈ ( a , a + ε] for ( σ − a ) σ d1 + d2 some real ε > 0. Assume that 0 < K < . Then, the unforced difference Equation ρ1 + ρ2 σ d1 (A1) has all its zeros in a < |z| < σ. As a result, if a < 1 and 0 < ε ≤ 1 − a then the unforced difference Equation (A1) has all its zeros in a < |z| < 1 so that it is stable. (ii) Let ci = ci (d1 , d2 ) ∈ R+ for i = 1, 2, be chosen such that ad1 < c1 ≤≤ σd1 ≤ c2 (implying that lnc lnσ − d2 < d1 ≤ 1 lnc2 lnσ . Then, the unforced difference Equation (A1) has its zeros in |z| < σ if   − d1 c1 − a c1 0≤K< . ρ1 + ρ2 c2 Proof. One has from the first identity of (A3) that: p0 (z) = p00 (z) + pe0 (z) (A36) Vaccines 2021, 9, 398 24 of 26 where:   p00 (z) = zd1 +d2 (z − a); pe0 (z) = ρ2 zd1 + ρ1 K (A37) Consider a circle |z| ≤ σ of center z = 0 and of radius σ ∈ ( a , a + ε] for some real ε > 0. Any point of the circumference |z| = σ, i.e., iθ √ the boundary of the circle |z| ≤ σ, is given by z = σe for some θ ∈ [0 , 2π ) with i = −1 being the complex unity. It turns out that all the zeros of p00 (z), which are z = a and z = 0 (with multiplicity d1 + d2 ), are in |z| < σ. From Rouché0 s theorem, all the zeros of p0 (z) are in a < |z| < σ if q q 2 σ d1 + d2 (σcosθ − a)2 + σ2 sin2 θ > K ρ2 σd1 cosd1 θ + ρ1 + ρ22 σ2d1 sin2 d1 θ; (A38) θ ∈ [0, 2π ) or,     σ2(d1 +d2 ) σ2 + a2 − 2aσcosθ > K2 ρ21 + ρ22 σ2d1 + 2ρ1 ρ2 σd1 cos(d1 θ ) ; θ ∈ [0, 2π ) (A39) The following cases can arise: Case a) If d1 is even then for θ ∈ [0, 2π ):     σ 2( d1 + d2 ) σ 2 + a 2 > K2 max ρ21 + ρ22 σ2d1 + 2 ρ1 ρ2 σd1 + aσ2(d1 +d2 ) cos(θ )  θ∈[0,2π )   = K2 ρ21 + ρ22 σ2d1 + 2 ρ1 ρ2 σd1 + aσ2(d1 +d2 ) (A40)  2 = K2 ρ1 + ρ2 σd1 + 2aσ2(d1 +d2 ) or, equivalently,   K ρ1 + ρ 2 σ d1 < σ d1 + d2 ( σ − a ) Case b) If d1 is odd then for θ ∈ [0, 2π ):     σ 2( d1 + d2 ) σ 2 + a 2 > K2 max ρ21 + ρ22 σ2d1 + 2 ρ1 ρ2 σd1 − aσ2(d1 +d2 ) cos(d1 θ )  θ∈[0,2π )  (A41) = K ρ21 + ρ22 σ2d1 + 2 ρ1 ρ2 σd1 − aσ2(d1 +d2 ) 2 which is still guaranteed under the above condition which is equivalent to 0 < σ d1 + d2 ( σ − a ) K< . The proof of Property (i) is complete. To prove Property (ii), note ρ1 + ρ2 σ d1 d that if c1 ≤ σd1 +d2 ≤ σd1 ≤ c2 then | a| < 1/c11 < σ. Those constraints guarantee  −d  c2 1 − a c1 σ d1 + d2 ( σ − a ) that 0 < K < ρ1 + ρ2 c2 implies that 0 < K < which proves Property ρ1 + ρ2 σ d1 (ii) as a direct consequence of Property (i).  References 1. Hethcote, H.W. The Mathematics of Infectious Diseases. SIAM Rev. 2000, 42, 599–653. [CrossRef] 2. Jiao, H.; Shen, Q. Dynamics Analysis and Vaccination-Based Sliding Mode Control of a More Generalized SEIR Epidemic Model. IEEE Access 2020, 8, 174507–174515. [CrossRef] 3. De la Sen, M.; Ibeas, A.; Alonso-Quesada, S. On vaccination controls for the SEIR epidemic model. Commun. Nonlinear Sci. Numer. Simul. 2012, 17, 2637–2658. [CrossRef] 4. De La Sen, M.; Nistal, R.; Ibeas, A.; Garrido, A.J. On the Use of Entropy Issues to Evaluate and Control the Transients in Some Epidemic Models. Entropy 2020, 22, 534. [CrossRef] [PubMed] 5. Liu, F.; Huang, S.; Zheng, S.; Wang, O.H. Stability Analysis and Bifurcation Control for a Fractional Order SIR Epidemic Model with Delay. In Proceedings of the 39th Chinese Control Conference (CCC), Shenyang, China, 27–30 July 2020; IEEE: Piscataway, NJ, USA, 2020; pp. 724–729. 6. Nistal, R.; De La Sen, M.; Alonso-Quesada, S.; Ibeas, A. On a New Discrete SEIADR Model with Mixed Controls: Study of Its Properties. Mathematics 2018, 7, 18. [CrossRef] Vaccines 2021, 9, 398 25 of 26 7. Ibeas, A.; De La Sen, M.; Alonso-Quesada, S.; Nistal, R. Parameter Estimation of Multi-Staged SI(n)RS Epidemic Models. In Proceedings of the UKACC 12th International Conference on Control (CONTROL), Sheffield, UK, 5–7 September 2018; IEEE: Piscataway, NJ, USA, 2018; pp. 456–461. 8. Wang, X.; Lu, J.; Wang, Z.; Li, Y. Dynamics of discrete epidemic models on heterogeneous networks. Phys. A Stat. Mech. Appl. 2020, 539, 122991. [CrossRef] 9. Parsamanesh, M.; Mehrshad, S. Stability of the equilibria in a discrete-time sivs epidemic model with standard incidence. Filomat 2019, 33, 2393–2408. [CrossRef] 10. Darti, I.; Suryanto, A.; Hartono, M. Global stability of a discrete SIR epidemic model with saturated incidence rate and death induced by the disease. Commun. Math. Biol. Neurosci. 2020, 2020, 33. [CrossRef] 11. Suryanto, A.; Darti, I. On the nonstandard numerical discretization of SIR epidemic model with a saturated incidence rate and vaccination. AIMS Math. 2021, 6, 141–155. [CrossRef] 12. Shalan, R.N.; Shireen, R.; Lafta, A.H. Discrete an SIS model with immigrants and treatment. J. Interdiscip. Math. 2020, 1–6. [CrossRef] 13. Córdova-Lepe, F.; Gutiérrez, R.; Vilches-Ponce, K. Analysis of two discrete forms of the classic continuous SIR epidemiological model. J. Differ. Equ. Appl. 2019, 26, 1–24. [CrossRef] 14. Radha, S.; Mahalakshmi, K.; Kumar, S.; Saravanakumar, A.R. E-learning during lockdown of COVID-19 pandemic: A global perspective. Int. J. Control Autom. 2020, 13, 1088–1089. 15. Cao, H.; Wu, H.; Wang, X. Bifurcation analysis of a discrete SIR epidemic model with constant recovery. Adv. Differ. Equ. 2020, 2020, 1–20. [CrossRef] 16. Youssef, H.M.; Alghamdi, N.A.; Ezzat, M.A.; El-Bary, A.A.; Shawky, A.M. A new dynamical modeling SEIR with global analysis applied to the real data of spreading COVID-19 in Saudi Arabia. Math. Biosci. Eng. 2020, 17, 7018. [CrossRef] [PubMed] 17. Alrashed, S.; Min-Allah, N.; Saxena, A.; Ali, I.; Mehmood, R. Impact of lockdowns on the spread of COVID-19 in Saudi Arabia. Inform. Med. Unlocked 2020, 20, 100420. [CrossRef] 18. De La Sen, M.; Ibeas, A.; Agarwal, R.P. On Confinement and Quarantine Concerns on an SEIAR Epidemic Model with Simulated Parameterizations for the COVID-19 Pandemic. Symmetry 2020, 12, 1646. [CrossRef] 19. De la Sen, M.; Ibeas, A. On a controlled SE (Is)(Ih)(Iicu) AR epidemic model with output controllability issues to satisfy hospital constraints on hospitalized patients. Algorithms 2020, 13, 322. [CrossRef] 20. Zhai, S.; Gao, H.; Luo, G.; Tao, J. Control of a multigroup COVID-19 model with immunity: Treatment and test elimination. Nonlinear Dyn. 2020, 1–15. [CrossRef] 21. Abbasi, Z.; Zamani, I.; Mehra, A.H.A.; Shafieirad, M.; Ibeas, A. Optimal Control Design of Impulsive SQEIAR Epidemic Models with Application to COVID-19. Chaos Solitons Fractals 2020, 139, 110054. [CrossRef] [PubMed] 22. Etxeberria-Etxaniz, M.; Alonso-Quesada, S.; De La Sen, M. On an SEIR Epidemic Model with Vaccination of Newborns and Periodic Impulsive Vaccination with Eventual On-Line Adapted Vaccination Strategies to the Varying Levels of the Susceptible Subpopulation. Appl. Sci. 2020, 10, 8296. [CrossRef] 23. Rohith, G.; Devika, K.B. Dynamics and control of COVID-19 pandemic with nonlinear incidence rates. Nonlinear Dyn. 2020, 101, 2013–2026. [CrossRef] [PubMed] 24. Singh, R.K.; Drews, M.; De La Sen, M.; Kumar, M.; Singh, S.S.; Pandey, A.K.; Srivastava, P.K.; Dobriyal, M.; Rani, M.; Kumari, P.; et al. Short-Term Statistical Forecasts of COVID-19 Infections in India. IEEE Access 2020, 8, 186932–186938. [CrossRef] 25. Tuite, A.R.; Ng, V.; Rees, E.; Fisman, D. Estimation of COVID-19 outbreak size in Italy. Lancet Infect. Dis. 2020, 20, 537. [CrossRef] 26. Dutta, A. Stabilizing COVID-19 Infections in US by Feedback Control based Test and Quarantine. In Proceedings of the 2020 IEEE Global Humanitarian Technology Conference (GHTC), Seattle, WA, USA, 29 October–1 November 2020; IEEE: Piscataway, NJ, USA, 2020; pp. 1–6. 27. Casella, F. Can the COVID-19 pandemic be controlled on the basis of daily test reports? IEEE Control Syst. Lett. 2021, 5, 1079–1084. [CrossRef] 28. Paul, S.; Lorin, E. Lockdown: A non-pharmaceutical policy to prevent the spread of COVID-19. Math. Modeling Comput. 2020. [CrossRef] 29. Balabdaoui, F.; Dirk, M. Age-stratified discrete compartment model of the COVId-19 epidemic with application to Switzerland. Sci. Rep. 2020, 10, 2306. [CrossRef] 30. Costa, J.A., Jr.; Martinez, A.C.; Geromel, J.C. On an alternative susceptible-infected removed epidemic model in discrete-time. Soc. Bras. Autom. 2020, 2, 2020. 31. Cooper, I.; Mondal, A.; Antonopoulos, C.G. A SIR model assumption for the spread of COVID-19 in different communities. Chaos Solitons Fractals 2020, 139, 110057. [CrossRef] [PubMed] 32. De la Sen, M.; Ibeas, A. On a Sir epidemic model for the COVID-19 pandemic and the logistic equation. Discret. Dyn. Nat. Soc. 2020, 2020, 1382870. 33. Ng, K.Y.; Gui, M.M. COVID-19: Development of a robust mathematical model package with consideration for ageing populatio- nand time delay for control action and resusceptibility. Phys. D Nonlinear Phenom. 2010, 411, 132599. [CrossRef] 34. Ng, M.K.Y. SEIRS-Based COVID-19 Simulation Packa ge. Available online: https://www.markusng.com/COVIDSiM/ (accessed on 2 January 2021). Vaccines 2021, 9, 398 26 of 26 35. Anderez, D.O.; Kanjo, E.; Pogrebna, G.; Kaiwartya, O.; Johnson, S.D.; Hunt, J.A. A COVID-19-Based Modified Epidemiologi- calModel and Technological Approaches to Help Vulnerable Individuals Emerge from the Lockdown in the UK. Sensors 2020, 20, 4967. [CrossRef] [PubMed] 36. Dobbs, P.S.; Watts, D.J. A generalized model of social and biological contagion. J. Theor. Biol. 2005, 232, 587–604. 37. Wang, Y.; Wei, Z.; Cao, J. Epidemic dynamics of influenza-like diseases spreading in complex networks. Nonlinear Dyn. 2020, 101, 1–20. [CrossRef] 38. Wang, Y.; Cao, J. Final size of network epidemic models: Properties and connections. Sci. China Inf. Sci. 2020, 64, 1–3. 39. Chinwenyi, H.C.; Ibrahim, H.D.; Adekunle, J.O. A study of two disease models: With and without incubation period. Int. J. Math.Comput. Sci. 2019, 13, 74–79. 40. Conway, A. Functions of One Complex Variable I; Springer: Berlin/Heidelberg, Germany, 1978; ISBN 0-471-99672-6. 41. Beardon, A. Complex Analysis: The Argument Principle in Analysis and Topology; John Wiley & Sons: New York, NY, USA, 1979. 42. De la Sen, M.; Barcena, R. Stability results for point time-delay systems obtained via Rouche’s theorem. Appl. Anal. 2004, 83, 157–183. [CrossRef] 43. Khanh Tu, N.; Erickson, K.T. Stability of discrete-time matrix polynomials. IEEE T. Automat. Contr. 1997, 42, 538–542. 44. Mendel, J.M. Discrete Techniques of Parameter Estimation: The Equation Error Formulation; Marcel Dekker Inc.: New York, NY, USA, 1973. 45. De la Sen, M.S.; Alonso-Quesada, A.I. On a Discrete SEIR Epidemic Model with Exposed Infectivity, feedback vaccination and partial delayed re-Susceptibility. Mathematics 2021, 9, 520. [CrossRef] 46. Niamsup, P.; Ratchagit, K.; Phat, V.N. Novel criteria for finite-time stbilization ans guaranteed cost control of delay systems delayed neural networks. Neurocomputing 2015, 160, 281–286. [CrossRef] 47. Phat, V.N.; Khongtham, Y.; Ratchagit, K. LMI approach to exponential stability of linear systems with interval time-varying delays. Linear Algebra Appl. 2012, 436, 243–251. [CrossRef] 48. Delasen, M. On some structures of stabilizing control laws for linear and time-invariant systems with bounded point delays and unmeasurable states. Int. J. Control. 1994, 59, 529–541. 49. Iacobucci, G.; Mahase, E. Covid-19 vaccination: “What’s the evidence for extending the dosing interval”? BMJ 2021, 372. [CrossRef]