mathematics Article On a Discrete SEIR Epidemic Model with Exposed Infectivity, Feedback Vaccination and Partial Delayed Re-Susceptibility Manuel De la Sen 1, * , Santiago Alonso-Quesada 1 and Asier Ibeas 2 1 Institute of Research and Development of Processes, IIDP, University of the Basque Country, 48080 Biscay, Spain; santiago.alonso@ehu.eus 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 pro- posed, and its properties of non-negativity and (both local and global) asymptotic stability of the solution sequence vector on the first orthant of the state-space are discussed. The calculation of the disease-free and the endemic equilibrium points is also performed. The model has the following main characteristics: (a) the exposed subpopulation is infective, as it is the infectious one, but their re- spective transmission rates may be distinct; (b) a feedback vaccination control law on the Susceptible is incorporated; and (c) the model is subject to delayed partial re-susceptibility in the sense that a partial immunity loss in the recovered individuals happens after a certain delay. In this way, a portion  of formerly recovered individuals along a range of previous samples is incorporated again to the  susceptible subpopulation. The rate of loss of partial immunity of the considered range of previous Citation: De la Sen, M.; samples may be, in general, distinct for the various samples. It is found that the endemic equilibrium Alonso-Quesada, S.; Ibeas, A. On a point is not reachable in the transmission rate range of values, which makes the disease-free one Discrete SEIR Epidemic Model with to be globally asymptotically stable. The critical transmission rate which confers to only one of the Exposed Infectivity, Feedback equilibrium points the property of being asymptotically stable (respectively below or beyond its Vaccination and Partial Delayed value) is linked to the unity basic reproduction number and makes both equilibrium points to be Re-Susceptibility. Mathematics 2021, 9, coincident. In parallel, the endemic equilibrium point is reachable and globally asymptotically stable 520. https://doi.org/10.3390/ in the range for which the disease-free equilibrium point is unstable. It is also discussed the relevance math9050520 of both the vaccination effort and the re-susceptibility level in the modification of the disease-free equilibrium point compared to its reached component values in their absence. The influences of Academic Editors: Martha L. Abell the limit control gain and equilibrium re-susceptibility level in the reached endemic state are also and James P. Braselton explicitly made viewable for their interpretation from the endemic equilibrium components. Some Received: 26 January 2021 simulation examples are tested and discussed by using disease parameterizations of COVID-19. Accepted: 21 February 2021 Published: 2 March 2021 Keywords: discrete epidemic model; feedback vaccination control; transmission rate of the exposed; transmission rate of the infectious; delayed re-susceptibility (or immunity loss); COVID-19 pandemic Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affil- iations. 1. Introduction Epidemic models have received much attention in the last decades. Many formula- tions used to describe them are based on either differential or difference equations leading to significant numbers of existing continuous-time and discrete-time epidemic models. Copyright: © 2021 by the authors. For instance, the relevance of the basic reproduction number is discussed in [1] by ob- Licensee MDPI, Basel, Switzerland. taining related threshold theorems and applications to pertussis and measles descriptions. This article is an open access article Feedback vaccination laws have been proposed by employing feedback control concepts distributed under the terms and and techniques as sliding-mode control or linear or impulsive feedback vaccination. See, conditions of the Creative Commons for instance [2,3] and references therein. It has been described in the mentioned research Attribution (CC BY) license (https:// how the various vaccination laws can affect the basic reproduction number, or the disease creativecommons.org/licenses/by/ transmission rate, and then the stability of both the disease-free equilibrium point and 4.0/). Mathematics 2021, 9, 520. https://doi.org/10.3390/math9050520 https://www.mdpi.com/journal/mathematics Mathematics 2021, 9, 520 2 of 32 the endemic one. The transient evolution of epidemic diseases is focused on in [4] to approximately predict the dates of maximum hospital occupancy of beds according to different epidemic models of continuous type, and the dates can be monitored to some extent under vaccination controls. The stability analysis and bifurcation control for a frac- tional order Susceptible-Infectious-Recovered (SIR) epidemic model with delay is focused on in [5]. The discretization of continuous-time epidemic models and the derivation of purely discrete epidemic models has also received much attention. See, for instance, [6–12] and references therein. It can be pointed out that, since the significant periods in a dis- ease dynamics are relatively long, typically for instance on the orders of days or weeks, discrete-oriented epidemic models might be appropriate for the disease evolution effi- cient description and also to efficiently perform the computational studies to predict and eventually to control its transmission. A way of controlling the disease transmission is an appropriate programming of the vaccination campaigns. The above consideration is reinforced since the public vaccination policies are typically programmed for their man- agement and administration to the population along coming periods of days or weeks, which are considered long periods for modelling discretization purposes and concerns that facilitate the decision implementations. Discretization techniques are in fact very useful in practical situations to accommodate the data acquisition from the disease evo- lution to its concrete use, interpretation and public decisions about the disease control when running through time mathematical epidemic models. In particular, an extended Susceptible-Exposed-Infectious-Recovered (SEIR) model which incorporates asymptomatic and dead infective subpopulations to the classical SEIR model has been discussed in [6] where the properties of stability of the equilibrium points and model positivity have been investigated in detail. Regular and impulsive feedback controls have been incorporated to the basic model. Such a model has been claimed to be of potential usefulness for Ebola evolution description. A multi-staged epidemic discrete model has been discussed in [7]. Although the exposed subpopulation is not included in such a model, it considers several coupled layers of the infectious subpopulation where the infection is transmitted from each layer to the adjacent ones. Discrete Susceptible-Infectious-Recovered-Susceptible (SIRS) models have also been studied in heterogeneous networks [8] with related bifurcation and stability studies. On the other hand, some epidemic models consider the vaccination actions as a generator of a new subpopulation, the vaccinated one, which is dynamically coupled to the remaining ones in the model rather than as a specific forcing control. See, for instance, [9]. In addition, discretized Susceptible-Infectious-Recovered (SIR)-type models, or Susceptible-Infectious-Susceptible (SIS)-type ones with stability formulations have been proposed in [10–13] and some of the references therein. Over the last year, a lot of research has been developed concerning the recent COVID- 19 pandemic. Since the authorities of the different countries supply very frequently disease registered data, typically day-to-day, such as infection and detection tests, hospital bed occupancies and, in particular, intensive unity care occupancies and mortality-related data, discrete epidemic models have been seen to be appropriate in the related research. See, for instance, [14–30], and some related references therein. The existing analysis of data is abundant for different countries or regions such as, for instance, Saudi Arabia [16,17], Madrid and the surrounding administrative area [18,19], India [23,24], Italy [25], United Sates [26], Canada and several of its provinces [28], Switzerland, [29], Brazil [30], etc. See also [31,32] for other recent studies on the spread of COVID-19. 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. The COVID-19 disease is currently evolving with alternate intensive growing-up periods of time followed by damping periods and vice-versa, which have regional characteristics and are typically governed by public interventions, such as isola- tions, total or partial quarantines, movement restrictions, social distance keeping, mask use duties, lock down of commerce, teaching and non-essential activities, etc. In this way, the impact of lockdowns is investigated in [17]. The effects of total or partial quarantines is Mathematics 2021, 9, 520 3 of 32 investigated in [18] for a SEIAR model, which incorporates the asymptomatic subpopula- tion to the SEIR model by considering the isolated population as removed either from the infectious individuals or from the susceptible ones through impulsive control actions from the infection chain. The approach in [18] is formulated in the continuous-time framework rather than in a discrete-time one. The whole extended SEIR model under consideration includes the asymptomatic individuals as a separate subpopulation with its own transition from the exposed individuals. In addition, in [19], a general model with three different infectious subpopulations, namely, the slightly infections, the hospitalized ones and those in the intensive care unit, are considered. Each one of the above subpopulations as well as the asymptomatic ones have their own transitions from the exposed subpopulation, and vaccination and treatment controls are also considered to satisfy a priori established hospitalization constraints for the seriously infectious individuals on availability of beds and specific needed hospital means for intensive treatment. On the other hand, the im- plementation of control strategies oriented to reduce the number of exposed individuals or to increase the number of treated individuals is proposed and discussed in [20], while impulsive optimal control techniques are developed in [21]. In particular, the proposal of [20] relies on the fact that the pandemic is now endemic. In [23], the inadequacy of open-loop controls is emphasized in contrast to closed-loop controls like, for instance, slid- ing mode-based control laws. On the other hand, it is proven in [27] that the suppression strategies might work if they are sufficiently strong and taken through early decisions, while the mitigation strategies can fail because of unfavourable combinations of delays, unstable dynamics and uncertainties. The new discrete SEIR model, which is proposed and discussed in this paper, includes linear feedback vaccination efforts on the susceptible, partial re-susceptibility (or partial loss of immunity) of previously recovered individuals and infectivity (transmission rate) in the Exposed-Susceptible contagious contacts which may be potentially distinct from its counterpart related to the Infectious-Susceptible contagious contacts. It can be pointed out that in some later periods of the disease incubation the exposed individuals (still without exhibiting external symptoms) can produce contagious contacts with susceptible individuals what can translate into an uncontrolled disease transmission along that phase of the disease evolution due to the difficulty of performing tests for those individuals (basically, due to their absence of symptoms), [33–35]. This circumstance is known to happen, for instance, in the COVID-19 pandemic evolution. So, this is an important motivating reason to consider infectivity levels not only in the Infectious-Susceptible contacts but in the Exposed-Susceptible ones as well. It is also known that some diseases like influenza, common viral cough, and very probably COVID-19 do not generate permanent immunity in the recovered subpopulation that was previously infected. Therefore, the proposed discrete model considers also a partial loss of immunity in the recovered subpopulation, which increases the current susceptible levels after a certain delay from the recovery date. The incorporation of a term to take account of potential re-susceptibility, which is subtracted from the recovered subpopulation and added to the susceptible one, is made in a versatile way in the sense that the weights for each considered interval of loss of immunity can be affected from distinct weights, and such an interval itself is parameterized by modeling delays. 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 eventual linear feedback vaccination and partial re-susceptibility (or partial loss of immu- nity) subject to a certain minimum threshold discrete delay and which weights such a lost on a set of consecutive delayed samples related to the current sample and previous to such a threshold delay. Three kinds of loss of immunity are discussed that evolve, respectively, either at exponential or lower than exponential rates, and also loss of immunity is consid- ered a constant rate. Section 3 discusses the non-negativity of the solution for bounded non-negative initial conditions and the disease-free and endemic equilibrium points that are unique and calculated explicitly. It is proven that the endemic equilibrium point does Mathematics 2021, 9, 520 4 of 32 not exist for the transmission rates for which the disease-free one is locally asymptotically stable, which coincides with the basic reproduction number having values below unity. In addition, the endemic one is only reachable (i.e., all its components are non-negative) as the disease- free one is unstable. The existence of a critical transmission rate corresponding to the basic reproduction number having the unity critical value (critical basic reproduction number) is also proven, where both equilibrium points coincide. The basic reproduction number, critical transmission rate and equilibrium points are seen to be dependent on the limit vaccination gains. In particular, if the vaccination gain increases, then the basic reproduction number decreases. Section 4 is devoted to the stability properties. Essentially, it is proven that there are no limit oscillations, and that the conditions of local asymptotic stability around the equilibrium points imply that those stability properties are also global since no limit oscillations exist. It is also proven that there is only a global asymptotic attractor (being either the disease-free equilibrium point or the endemic one) depending on the transmission rate being under or over its critical value. Section 5 presents and discusses some examples of the proposed model related to the evolution of the COVID-19 pandemic. Finally, conclusions end the paper. 2. Discrete SEIR Epidemic Model with Vaccination Control and Re-susceptibility It can be first pointed out that the use of a discrete epidemic model allows a direct acquisition of recorded data, and it is easier to implement related to continuous-time models from a computational point of view. Note that the rationale of the sampling period interpretation is that it is unity, typically, either 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. Consider the following Susceptible-Exposed- Infectious-Recovered (SEIR)-type discrete model with vaccination effort and partial loss of immunity (or re-susceptibility) given by Si+1 = (1 − βIi − βe Ei ) Si − Vi + LI R̂i  (1) Ei+1 = (1 − µ) Ei + ( βIi + βe Ei )Si (2) Ii+1 = (1 − γ) Ii + µEi (3)  Ri+1 = Ri + γIi + Vi − LI R̂i (4) for any integer i ∈ Z0+ = Z+ ∪ {0} where {Si }i∞=0 , { Ei }i∞=0 , { Ii }i∞=0 and { Ri }i∞=0 are, respectively, the susceptible, exposed, infectious and recovered sequences. {Vi }i∞=0 ⊂ R0+   ∞ and LI R̂i i=0 ⊂ R0+ are, respectively, the vaccination and partial loss of immunity max(0 , i −d ) sequences, where R̂i = R j j=max(0 , i−1 d −d ) reflects that the immunity is lost for a certain  1 2 range of d2 + 1 previous recovered values ranging until the d1 back-step sample before the current sample. Consider any finite initial conditions S0 ≥ 0, E0 ≥ 0, I0 ≥ 0 and R0 ≥ 0. In the above model, note the following: (a) β and βe are, respectively, the transmission rates of the infectious and exposed sub- populations, γ is the recovery rate and µ is the incubation rate The consideration of infectivity to the Susceptible from both the Exposed and the Infectious is consid- ered since recent studies from the COVID-19 pandemic have established that the last period of the incubation stage and the former period of the infectious one are both contagious. In this paper, we consider that both transmission rates might be distinct in general. max(0 , i −d ) The term LI R̂i , where R̂i = R j j=max(0 , i−1 d −d ) , which is subtracted from the   (b) 1 2 contribution to the recovered subpopulation while added to the susceptible subpopu- lation, takes into account the partial loss of immunity in the recovered subpopulation. Its calculation involves a range of previous samples of the discrete model that can be adjusted from d1 + d2 samples to d1 samples prior to the current sampling time instant. The sampling period can be typically one day, but this is not a restriction of the model. Mathematics 2021, 9, 520 5 of 32 The delays d1 and d2 parameterize the model, and they can be varied or monitored for different model evaluations. Some ideas about how to implement the term of loss of immunity are discussed later on in Remark 2 with several possible weights that relate the loss of immunity with the time elapsed since recovery. The basic idea is the that the degree of re-susceptibility or of loss of immunity might be gradual so that the effective loss of immunity increases as time increases since recovery. This modeling proposal is introduced in the discrete epidemics framework in this paper. In the continuous framework, re-susceptibility has been considered in several papers. See, for instance, [36–39]. In particular, in [36], the contribution of re-susceptibility to the Susceptible and Recovered dynamics is formulated through integrals that involve delays in their integration limits. For generality purposes, we can consider that the period of loss of immunity can be parameterized by two delay parameters, and the effective loss of immunity along its duration period can be weighted, if desired, with different average weights for each sample depending on how far is it from the current sample under consideration. On the other hand, note that the given non-negativity conditions are set for the general mathematical setting. For a proper spread of the infection, S0 > 0, since otherwise, Si = 0; ∀i ∈ Z0+ , even in the presence of initial infection, that is, even if min( E0 , I0 ) > 0, in the event that there is no loss of immunity in the recovered subpopulation. It is argued that this model may be of interest for the evaluation of the COVID-19 middle-term or long-term disease propagation under eventual potential feedback vaccina- tion in its blowing-up phase because of the following facts: - It is of simple structure and of a discrete nature and, therefore, very appropriate for computational experiments. Furthermore, it does not need the incorporation of a modulation discretization parameter to guarantee the non-negativity of the sequence solution as other epidemic discrete models usually need. See, for instance, [1,2] and some references therein. - It has two eventual potentially distinct coefficient transmission rates βe and β for the exposed and infectious subpopulations, which allow to consider infectivity to the susceptible subpopulation from both infective subpopulations. This might be potentially advantageous for its use for description of COVID-19 since it is now known that this disease has an infective period at the end of the incubation period and another one along the first days of the symptomatic infectious period. - It might be relevant for the forthcoming studies related to COVID-19 to evaluate the possible existence of an endemic steady-state. Note that the disease is now behaving in blowing-up phases, which follow the different intervention measures of confinement, isolations and other measures like social distancing, leisure control measures, etc., which are being taken in most of the countries. However, vaccination is starting to progress, and it is foreseen that its influence can reduce the disease force and the appearance of an endemic steady state, for which infection, susceptibility and immunity levels can depend on the vaccination effort and the transmission rate. - It is seen later on from the mathematical study of this model that the susceptible subpopulation is a decreasing sequence and that the exposed and infectious sub- populations increase for consecutive samples under certain conditions of the disease parameters and upper-bounds for the susceptible subpopulation. This behavior is also of interest for the use of the model as the disease blows-up along a transient period of time. Remark 1. Note that by summing-up Equations (1)–(4), one gets directly that the total population is constant, then satisfying Ni = Si + Ei + Ii + Ri = N0 for any integer i ∈ Z0+ . As a result, if the total population is initially unity then the model remains as a normalized model for all samples with the total population remaining unity through time. Mathematics 2021, 9, 520 6 of 32 3. Non-Negativity of the Solution Sequence and Equilibrium Points The following result gives sufficiency-type conditions for the non-negativity of the solution sequence of Models (1)–(4) under non-negative bounded initial conditions. The proof is based on a complete induction method from the non-negativity assumption for the initial conditions. Theorem 1. Assume that Vi = Ki Si with {Ki }i∞=0 ⊂ [0 , 1] and, furthermore, the following constraints hold: (1) µ ∈ (0 , 1], γ ∈ (0 , 1] and β/βe ≤ γ/µ with β > 0 and βe > 0. (2) The initial conditions are subject to 0 ≤ S0 ≤ µ/(2βe ), R0 ≥ 0, 0 ≤ E0 ≤ 1/(2βe ) and 0 ≤ I0≤ 1/(2β).  (3) LI R̂i /Si ≤ Ki ≤ 1 if Si 6= 0 and if Si = 0 then LI R̂i = 0 and Ki = 0; ∀i ∈ Z0+ . Then, the vector sequence solution is non-negative for all samples. In addition, the sequences of all the subpopulations and that of the total population are bounded. Proof. Note that from Condition 3, one concludes from (1)–(4) that Si+1 = (1 − Ki − βIi − βe Ei )Si + LI R̂i ≤ (1 − βIi − βe Ei ) Si  (5)  Ri+1 = Ri + γIi + Ki Si − LI R̂i ≥ Ri + γIi (6) The proof goes ahead by complete induction. Assume that there exists an integer k ≥ 0 such that min(Si , Ei , Ii , Ri ) ≥ 0, Si ≤ µ/2βe , Ei ≤ 1/2βe , Ii ≤ 1/2β; i = 0, 1, . . . , k. Such an integer k ≥ 0 always exists since the above constraints hold from the given hypotheses on the initial conditions (at least) for k = 0. It has to be proven that the constraints still hold for k + 1. First, note that, since Ek ≤ 1/2βe , Ik ≤ 1/2β and Sk ≥ 0, it follows that Sk+1 = (1 − Kk − βIk − βe Ek )Sk + LI R̂k ≥ 0  (7) and, since Ek ≥ 0 and Ik ≥ 0, one has that Ik+1 = (1 − γ) Ik + µEk ≥ 0 (8) where γ ∈ (0 , 1] has been used. On the other hand, since Ek ≤ 1/2βe , Ik ≤ 1/2β, 0 ≤ Sk ≤ µ/(2βe ), γ ∈ (0 , 1], µ ∈ (0 , 1] and β/βe ≤ γ/µ, one has that 0 ≤ Ek+1 = (1 − µ) Ek + ( βIk + βe Ek )Sk −µ 1− µ (9) ≤ 12β µ 1 e + Sk ≤ 2βe + 2βe = 2βe 1−γ µ 1−γ γ 1 0 ≤ Ik+1 = (1 − γ) Ik + µEk ≤ + e ≤ + = (10) 2β 2β 2β 2β 2β where β/βe ≤ γ/µ has been used. In addition, since Sk ≤ µ/(2βe ) and 1 − βIk − βe Ek ≥ 0 since 0 ≤ Ek ≤ 1/2βe ; 0 ≤ Ik ≤ 1/2β, it follows that µ 0 ≤ Sk+1 = (1 − βIk − βe Ek ) Sk − Vk + LI R̂k ≤ Sk ≤  (11) 2βe Additionally, since Rk ≥ 0, Ik ≥ 0 and γ > 0 then Rk+1 ≥ 0. It has been proven that (min(Si , Ei , Ii , Ri ) ≥ 0; i = 0, 1, . . . , k) ⇒ (min(Si , Ei , Ii , Ri ) ≥ 0; i = 0, 1, . . . , k + 1) (12) (Sk ≤ µ/(2β ) , Ek ≤ 1/2βe , Ik ≤ 1/2β ; i = 0, 1, . . . , k) e ⇒ (Sk ≤ µ/(2βe ) , Ek ≤ 1/2βe , Ik ≤ 1/2β ; i = 0, 1, . . . , k + 1) Mathematics 2021, 9, 520 7 of 32 As a result, the proof is complete since (min(S0 , E0 , I0 , R0 ) ≥ 0) ⇒ (min(Si , Ei , Ii , Ri ) ≥ 0; ∀i ∈ Z+ ) (S0 ≤ µ/(2βe ), E0 ≤ 1/2βe , I0 ≤ 1/2β) (13) ⇒ (Si ≤ µ/(2βe ), Ei ≤ 1/2βe , Ii ≤ 1/2β; ∀i ∈ Z+ which implies also that the sequences of all the subpopulations and that of the total population are bounded.  Note that the above result also holds in the particular cases of absence of vaccination and re-susceptibility. In particular, it can be concluded that, in the case of∞ absence of vaccination and partial immunity loss, that is, if {Vi }i∞=0 ≡ 0 and LI R̂i i=0 ≡ 0, one has that, under Constraints 1 and 2 of Theorem 1, and because of the inequality of (5), it follows that 0 ≤ Si ≤ µ/(2βe ), Ri ≥ 0, 0 ≤ Ei ≤ 1/(2βe ) and 0 ≤ Ii ≤ 1/(2β); ∀i ∈ Z0+ so that, in view of (5)–(6), those boundedness and non-negativity solution conditions still hold in this particular case.  Three simple conditions on the evaluation of the amount LI R̂i describing the partial immunity loss, and some of their respective implications related to the vaccination gain Ki and the solution non-negativity conditions, are now discussed in the subsequent remark. Remark 2. Some claimed mechanisms of losing the temporary immunity are now discussed: i −d (1) Exponentially fast partial immunity loss with delay: if LI R̂i = ∑ j=i−1 d −d a j R j with  1 2 ai−d1 −k = a(i )k+1 for 0 ≤ k ≤ d2 with some given delays d1 > 0 and 0 ≤ d2 ≤ ∞ and suppose that a = a(i ) ∈ [0 , 1); ∀i ∈ Z0+ , which describes an exponentially fast partial immunity loss after a certain delay d1 . If the partial immunity loss parameter a is zero, this means that the immunity is permanent. Otherwise, the immunity is lost at exponential rate after the d1 preceding samples to the current one, while for the d1 − 1 immediately previous samples the immunity is maintained. Thus, one has   a 1 − a d2 +1 i − d1 ∑ aj Rj  0≤ min R j ≤ LI R̂i = 1−a max(i −d1 −d2 ,0)≤ j≤max(i −d1 ,0) j =i − d − d 1 2     a 1 − a d2 +1 a 1 − a d2 +1 ≤ max Rj ≤ R i − d1 1−a max(i −d1 −d2 ,0)≤ j≤max(i −d1 ,0) 1−a Since from Condition 3 of Theorem 1 in (4), the sequence{ Ri }i∞=0 is non-decreasing so that R i − d1 = max R j and i − d1 − d2 ≤ j ≤ i − d1   ∞ ∞ ∞ ∞ d + 2 a 1 − a d2 +1 a a 2 a + a 2 + · · · + a d2 +1 = ∑ a j − a d2 +2 ∑ a j = a ∑ a j − a d2 +2 ∑ a j = − = j =1 j =0 j =0 j =0 1 − a 1−a 1−a a(1− ad2 +1 ) Rmax (i−d , 0)   Thus, Condition 3 of Theorem 1 is guaranteed for a normalized model if Ki ∈ (1− a ) S 1 , 1 i a ( 1 − a d2 +1 ) R i − d provided that ( 1 − a ) Si 1 ≤ 1 with Ri−d1 ≤ R M = sup Ri ≤ 1; ∀i ∈ Z0+ , which 0≤ i ≤ ∞ requires in turn as necessary condition compatible with the non-negativity conditions of The- µ a ( 1 − a d2 +1 ) orem 1 that 2βe ≥ Si ≥ 1− a Rmax (i−d1 , 0) subject in turn to the necessary condition µ (1− a ) µ (1− a ) βe ≤ , which is guaranteed if βe ≤ since R M ≤ 1. If d2 = ∞ (i.e., the 2a(1− ad2 +1 ) R M 2a(1− ad2 +1 ) µ (1− a ) immunity is lost backward in time after a delay); then, the above condition reduces to βe ≤ 2a . (2) Constant partial immunity loss with delay: If ai−d1 −k = a for d1 > 0, 0 ≤ k ≤ d2 < ∞ and a ∈ [0 , 1); ∀i ∈ Z0+ then Mathematics 2021, 9, 520 8 of 32 i − d1 ∑  0 ≤ LI R̂i = a j R j ≤ a(d2 + 1) Rmax j = i − d1 − d2   a(d2 +1) Rmax (i−d , 0) Thus, Condition 3 of Theorem 1 is guaranteed if Ki ∈ Si 1 , 1 ; ∀i ∈ Z0+ under the necessary condition compatible with the non-negativity conditions of Theorem 1 µ that 2βe ≥ Si ≥ a(d2 + 1) Rmax (i−d1 , 0) ; ∀i ∈ Z0+ subject in turn to the necessary condition βe ≤ which is guaranteed if βe ≤ µ µ 2a(d2 +1) Rmax (i−d , 0) , 2a(d2 +1) . 1 (3) Non-constant partial immunity loss at a smaller rate than exponential: It can be deduced in a similar way to the above situation with ai−d1 = ai ∈ [0, 1) and ai−d1 −k = ai − ε i−d1  k =max(d2 , i −d1 ) for k ∈ [1, d2 ] ∩ Z and ε i−d1 −k k=0 being a strictly increasing finite sequence; ∀i, k(≤ i − d1 ) ∈ Z0+ . If the parameterizing sequences are constant, i.e., and ε i = ε ∈ (0, a); ∀i, k(≤ i − d1 ) ∈ Z0+ then ∑ij− d1 =i −d −d a j R j ≤ ( a ( d2 + 1) − εd2 ) Rmax . 2 1 The disease-free and endemic equilibrium points are now calculated provided that the vacci- nation gains converge asymptotically to non-negative limits that can be, in general, distinct for a ( 1 − a d2 +1 ) the disease-free and endemic equilibrium points. Note that A ≡ Ae = 1− a in the case of exponentially fast partial immunity loss with delay, A ≡ Ac = a(d2 + 1) in the case of constant partial immunity loss with delay and A ≡ Anc = ( a − ε)d2 + a in the case of non-constant immunity loss at a smaller rate than exponential if the parameterizing sequences are constant. Proposition 1. (Disease-free equilibrium point). Assume that A > 0, that there is a partial immunity loss and that Vi = Ki Si ; ∀i ∈ Z0+ , where {Ki }i∞=0 (⊂ [0, 1]) → Kd f ∈ [0 , 1] . Then, there is a unique disease- free equilibrium point ! A Kd f T xd f = , 0 , 0, A + Kd f A + Kd f Proof. Note that the following equilibrium equations hold:   Sd f = (1 − βIi − βe Ei ) Sd f − Vd f + LI Rd f (14) Ed f = Id f = 0 (15)   LI Rd f = ARd f = Vd f (16)   with Vd f = Kd f Sd f = Kd f 1 − Rd f with Kd f = lim Ki > 0. It follows that i→∞ A Kd f Sd f = ; Rd f = (17) A + Kd f A + Kd f   1− a 1− a 1 provided that a < 1, that is a < min 1, = if A = Ae and a < if 1 − a d2 +1 1 − a d2 +1 d2 +1 A = Ac .  Remark 3. The subsequent observations have an interesting biological concern: (a) Note from Proposition 1 Rd f /Sd f = Kd f /A, that is, the relation limit Recovered to Suscepti- ble improves as the vaccination gain increases, and the partial immunity loss level decreases, as expected. (b) Note also that the disease-free equilibrium point given in Proposition 1 stands for that A > 0 and Kd f ∈ [0 , 1]. In particular, if Kd f = 0 then Sd f = 1 and Rd f = 0, which implies that the Susceptible are not asymptotically removed, while the total population becomes Mathematics 2021, 9, 520 9 of 32 asymptotically susceptible if the vaccination effort is removed and there is a permanent level of partial immunity loss. On the other hand, if A = Kd f = 0 then the susceptible and recovered subpopulations at the disease-free equilibrium point are both indeterminate in the form 0/0 according to Proposition 1. However, the inspection of (1) concludes that, in this case, {Si }i∞=0 is non-increasing so that lim (Si+1 − Si ) = 0 and {Si }i∞=0 → Sd f ≥ 0 for some existing i→∞ limit Sd f , eventually a residual nonzero limit susceptibility amount, depending on each given set of non-negative initial conditions, so that Rd f = 1 − Sd f as a result. Generally speaking, the vaccination injection reduces the values of the susceptible at the disease-free equilibrium and increases the equilibrium values of the recovered. Contrarily, the equilibrium values of the susceptible increase as the re-susceptibility force increases at the disease-free equilibrium point. Proposition 2. (Endemic equilibrium point). Assume that A > 0, that there is a partial immunity loss and that Vi = Ki Si ; ∀i ∈ Z0+ , where {Ki }i∞=0 (⊂ [0, 1]) → Kend ∈ [0 , 1] . Thus, there is a unique endemic equilibrium point xend = (Send , Eend , Iend , Rend ) T , where γµ Send = (18) βµ + βe γ γ Aµ( β − γ) + γ( Aβe − µKend ) Eend = (19) γ(µ + A) + µA βµ + βe γ µ Aµ( β − γ) + γ( Aβe − µKend ) Iend = (20) γ(µ + A) + µA βµ + βe γ µγ [µ( β − γ) + γβe ] + Kend γµ(γ + µ) Rend = (21) [γ(µ + A) + µA] ( βµ + βe γ)   µγ Kend which exists, that is, its components are positive, provided that β > β c = µ+γλe 1+ A , where βA(µ+γλe ) λe = βe /β, and equivalently, if the basic reproduction number R0 defined by R0 = µγ( A+Kend ) = A β c A(µ+γλe ) Send ( A+Kend ) satisfies the constraint R0 > R0c = µγ( A+Kend ) = 1. Proof. The proof of the endemic equilibrium components (18)–(21) follows by direct calculations by noting that the following equilibrium  equations  hold. From  (3),one gets Iend = µ γ Eend which replaced in (2) yields µEend = βµ γ + βe Eend Send = β µ γ + λe Eend Send γµ γµ which yields, since Eend > 0, Send = βµ+ βe γ = β (µ+λe γ) . Now, one gets from (4) with  µ γµ LI R̂end = ARend , Iend = γ Eend , Send = βµ+ βe γ , Rend = 1 − Send − Eend − Iend and Kend = lim Ki > 0 that i→∞ Kend γµ ARend = γIend + Vend = γIend + Kend Send = µEend + βµ+ βe γ h   i (22) = A 1 − 1 + γµ Eend − βµ+ γµ e β γ By comparing the two last equalities, one gets γµ + Aγ + µA A( βµ + βe γ) − Aγµ − Kend γµ Eend = (23) γ βµ + βe γ which leads directly to (19) and is positive if the right-hand-side is positive, that is if A[µ( β− γ) + γβ e e e  ] > µγKend , which requires, in turn, if λ = β /β, that β > β c = µγ Kend βA(µ+γλe ) µ+γλe 1 + A , equivalently, that the basic reproduction number R0 = µγ( A+K ) = end A β c A(µ+γλe ) Send ( A+Kend ) > R 0c = 1 and the critical basic reproduction number is R 0c = µγ( A+Kend ) =1 Mathematics 2021, 9, 520 10 of 32 µ as a result. On the other hand, one gets (20) from (19) since Iend = γ Eend . Finally, note from (4) with Rend = Ri = Ri+1 ; ∀i ∈ Z0+ leading to γ( Aβe −µKend ) Rend = µ A Eend + Kend γµ A( βµ+ βe γ) = Aµ γ(µ+ Aγ)+µA Aµ( β−γ)+ βµ+ βe γ + Kend γµ A( βµ+ βe γ) e − µK (24) = µγ [ Aµ( β−γ)+Aγ[γ((Aβ end )]+ Kend γµ [ γ ( µ + A )+ µA ] µ+ A)+µA] ( βµ+ βe γ) which can be rewritten as (21) after simplification in the numerator and cancellation of the coefficient A in the numerator and denominator.  Note from (18)–(21) that the endemic equilibrium susceptible subpopulation gradu- ally extinguishes as the transmission rate increases irrespective of the limit vaccination control and the re-susceptibility level. However, the limit vaccination gain decreases the equilibrium exposed and infectious levels, while the re-susceptibility force value at the equilibrium increases them. At the same time, the recovered equilibrium value increases as the equilibrium vaccination effort increases, while it decreases as the re-susceptibility force at the equilibrium increases. The above qualitative results agree, which is expected from the vaccination control and the re-susceptibility levels. Note also by summing-up (18)–(21) and both expressions in (17) that Ni = Nend = Nd f = 1; ∀i ∈ Z+ if N0 = 1 as expected from the normalization of the model (see Remark 1). Proposition 2 relies on the case that the endemic equilibrium point has positive components. The next result is concerned with the case when the endemic equilibrium point coincides with the disease-free one. It is shown that this happens for the critical values of the disease- free transmission rate, equivalently, and the basic reproduction number equalizes its critical unity value. Proposition 3. Assume that there exists some partial loss of immunity level, that is A > 0, and β = β c for a vaccination law Vi = Ki Si ; ∀i ∈ Z0+ , and any limit given vaccination gain K ∈ [0 , 1], with K = Kend = Kd f , such that {Ki }i∞=0 → K . Then, the disease-free and endemic   A K T. equilibrium points are confluent, that is, xend = xd f = A+ K , 0 , 0, A+K   Proof. Note from (19)–(20) that if β = β c = µ+γλe 1 + K µγ A then Eend = Iend = 0. The   relation µEend = β γ + λe Eend Send derived in the proof of Proposition 2 for Eend > 0, µ used to calculate the susceptible subpopulation at the equilibrium, is now indeterminate since Eend = 0 and (18) is not valid. So, for Eend = Iend = 0, the equilibrium point is A K T xd f = A+K , 0 , 0, A+K for K = Kd f (Proposition 1), which implies in turn from (19)–(20) with K = Kend = Kd f that Aµ( β − γ) + γ( Aλe β − µK ) = 0 so that xend = xd f , provided   that A > 0, if and only if β = β c = µ+γλe 1 + K µγ A for any given limit vaccination gain K ∈ [0 , 1].  Remark 4. From Propositions 2 and 3, important biological  conclusions  might be derived from µγ Kend the critical disease transmission rate β c = µ+γλe 1+ A , its associate (critical) basic re- βcA(µ+γλe ) production number R0c = µγ( A+Kend ) = 1 and the basic reproduction number R0 defined by βA(µ+γλe ) A R0 = µγ( A+Kend ) = Send ( A+Kend ) as follows: (a) If β < β c , equivalently if R0 < R0c , then the endemic equilibrium point is not reachable since it has some negative components incompatible with the non-negativity of the solution sequence. (b) If β ≥ β c , equivalently if R0 ≥ R0c , then the endemic equilibrium point is reachable since it has no negative component. (c) If β = β c , equivalently if R0 = R0c , then the disease-free equilibrium point and the endemic one are confluent. Mathematics 2021, 9, 520 11 of 32 (d) The basic reproduction number has an inverse proportionality ratio (so it decreases) with the endemic susceptible level. (e) The critical disease transmission rate increases with the vaccination gain, while the critical basic reproduction number decreases in parallel. That means the vaccination effort increases the range of potential transmission rate values for which the endemic equilibrium point is not reachable, that is, the only equilibrium point for such an increased range of potential values is the disease-free one. (f ) From (18)–(21), it is seen that the susceptible endemic equilibrium value does not depend on the vaccination effort. However, the recovered endemic equilibrium value increases as the vaccination limit gain increases while the exposed and infectious endemic equilibrium values both decrease as the vaccination effort increases. On the other hand, the recovered endemic equilibrium decreases as the partial loss of immunity becomes increased, as intuitively expected. 4. Stability Results This section relies on the main stability results on the epidemic model linked to the two existing equilibrium points. 4.1. Stability around the Disease-Free Equilibrium Point The local stability of the disease-free equilibrium point for sufficiently small state perturbations is addressed in the subsequent result: Theorem 2. i(Local asymptotic stability of the disease-free equilibrium point). Assume that λ ∈ 0 , 2−γ ,Vi = Ki Si ∈ [0, 1]; ∀i ∈ Z0+ with {Ki }i∞=0 → K with K = Kd f = Kend . Then, h e µ the disease-free equilibrium point is locally asymptotically stable in some neighborhood centred at itself if β ∈ [0 , β c ), it is critically stable if β = β c and it is instable if β > β c . Proof. The linearized model around the disease-free equilibrium point of the infective Equations (2) and (3) becomes ∆Ei+1 1 − µ + βλe Sd f ∆Ei      βSd f = (25) ∆Ii+1 µ 1−γ ∆Ii where ∆Ei and ∆Ii denote the first-order incremental values of the exposed and infectious related to the equilibrium levels. The characteristic equation of (25) is     z2 + γ + µ − 2 − βSd f λe z + µ − 1 − βSd f λe (γ − 1) − µβSd f = 0 (26) whose rootsrare  2     2 + βSd f λe − γ − µ ± 2 + βSd f λe − γ − µ + 4 µβSd f − µ − 1 − βSd f λe (γ − 1) z1,2 = (27) 2 2 The radicand in (27) is β2 Sd2 f λe + 2βSd f λe (γ − µ) + (γ − µ)2 + 4µβSd f ≥ 0 so that both roots are real. It follows then that z1,2 ∈ (−1, 1) so that (25) is an asymptotically stable discrete linear system if and only if r  2       ± 2 + βSd f λe − γ − µ + 4 µβSd f − µ − 1 − βSd f λe (γ − 1) ∈ γ + µ − βSd f λe − 4, γ + µ − βSd f λe (28)   Using the trivial identity γ + µ − βSd f λe = 2 − 2 − γ − µ + βSd f λe in (28), one obtains that largest root is less than one if  2 h  i 2 − γ − µ + βSd f λe − 4(µ − 1)(γ − 1) + 4βSd f λe (γ − 1) + 4µβSd f < 2 − 2 − γ − µ + βSd f λe 2 (29) Mathematics 2021, 9, 520 12 of 32 and, equivalently, if and only if Kd f   µγ µγ β < βc = = 1 + (30) ( µ + λ e γ ) Sd f (µ + λe γ) A r On the other hand, the smaller root is larger than minus one if and only if  2     2 + βSd f λe − γ − µ + 4 µβSd f − µ − 1 − βSd f λe (γ − 1) + γ + µ − βSd f λe < 4 (31) µγ−2(µ+γ)+4 Note that, since λe ≤ µ which holds for β ∈ [0 , β c ) where β c = Sd f [µ+λe (γ−2)] . 2− γ , β ∈ [0 , βc ) so that β ∈ [0 , β c ) guarantees the stability of both z1,2 . Then, the linearized sys- tem around the disease-free equilibrium point is locally asymptotically stable if β ∈ [0 , β c ), critically stable if β = β c and instable if β > β c . It has been proven that ∆Ei → 0 and ∆Ii → 0 as i → ∞ for some sufficiently small initial perturbations of the disease- free equilibrium point. It remains to be proven that the substate of the susceptible and recovered subpop- ulations of the linearized system for small perturbations around the disease-free equi- librium point also converges asymptotically to zero. From (1) and (4), one gets for i −d LI R̂i = ∑ j=i−1 d −d a j R j , by assuming that Ki = Kd f ∈ [0 , 1] as i → ∞ , the following  1 2 incremental auxiliary discrete system: ∆Si+1 ∆Si     1−K 0 0 · · · 0 a i − d1 · · · a i − d1 − d2   ∆Ri+1 ∆Ri 0 · · · 0 − a i − d1 · · · − a i − d1 − d2     K 1 ∆Ri ∆Ri−1              ..   0 1 0 0 0   ..   + ∆vi        . = 0 0 1 0 0   . (32)         ∆Ri−d1 +1     0 0 0 0 0     ∆Ri−d1    ..   0 0 0 1 0 0   ..   .   .  0 0 0 1 0 ∆Ri−d1 −d2 +1 ∆Ri−d1 −d2 where {vi }i∞=0 is a sequence that depends on the sequences {∆Ei }i∞=0 and {∆Ii }i∞=0 via (1)–(4) such that {vi }i∞=0 → 0 . If the state vector sequence and the matrix of dynamics of the above system are denoted by {∆xi }i∞=0 and A L , it follows that m −1 ∆xn+m = Am L ∆xn + ∑ Am L − i −1 ∆vn+i (33) i =0 Now, for any given integer and real constants m and ε > 0, there exists a non- m −1 m − i −1 negative integer N = N (ε, m) such that ∆xn+m − Am L ∆xn = ∑i =0 A L ∆vn+i ≤ ε; ∀n ≥ N and there is also a strictly decreasing positive real subsequence ε nk ∞  k =0 , with ε nk = ε nk ( Nk ), ∀k ∈ Z0+ , for a strictly increasing sequence of non-negative integers { Nk }∞ k =0 such that m −1 ∆xnk +m − Am L ∆xnk = ∑ Am L − i −1 ∆vnk +i ≤ ε nk ; ∀nk ∈ [ Nk , Nk+1 ) ∩ Z0+ , ∀k ∈ Z0+ (34) i =0 m −1 lim ∆xnk +m − Am k→∞ L ∆xnk = lim k → ∞ i =0 ∑ Am L − i −1 ∆vnk +i = 0 (35) ∞ Since the solution sequence is bounded, the sequences {∆xi }i∞=0 and Am L ∆xnk k =0  cannot diverge with opposed signs so that they are either oscillatory or they converge to some limits. Note that the eigenvalues of A L are all real of values 0, 1 − K and 1, so that the limit steady-state solution sequence verifies the following: (a) It cannot diverge; and Mathematics 2021, 9, 520 13 of 32 (b) It cannot be oscillatory for any period of m samples either. This occurs in spite of the matrix of dynamics A L is critically stable with one (respectively, two) positive unity  m eigenvalue(s) if K ∈ (0 , 1] (respectively, if K = 0). Therefore,{∆xi }i∞=0 and ∞ A L ∆xnk k=0 both converge to limit values of opposed sign. Those limit values cannot be nonzero since then the disease-free equilibrium point would be non-unique, contradicting Proposition 2. As a result, those limits are both zero, and the susceptible and recovered subpopulations are locally asymptotically stable around the disease- free equilibrium point. It has been fully proven that the linearized system is locally asymptotically stable around the disease-free equilibrium point.  Remark 5. Note that Theorem 2 excludes the presence of steady-state oscillations around the disease- free equilibrium point so that it is locally asymptotically stable. The infective sub-state has been proven to be locally asymptotically stable in the first part of the proof. The second part excludes that the linearized susceptible-recovered substate might exhibit an asymptotically oscillatory steady-state since its matrix of dynamics have no complex conjugate eigenvalues and the disease-free equilibrium point is unique. Remark 6. Note that the critical value of the transmission rate β c and the corresponding unity critical basic reproduction number R0c = 1 guaranteeing the local asymptotic stability in Theorem 2 are identical to those which give the coincidence of the endemic equilibrium point and the disease-free one in Proposition 2. As a result, it is ensured that, for ranges of corresponding values smaller than those critical values, the disease-free equilibrium point is locally asymptotically stable while the endemic one is not reachable. The above asymptotic stability is also global for β < β c since the endemic equilibrium point is not reachable, so that the only attractor is the disease-free equilibrium point, and the existence of non trivial oscillations is excluded. By trivial oscillations we understand the constant solutions corresponding to the disease equilibrium points and the zero solution. A more precise formulation of this concern follows: Vi = Ki Si ∈ Theorem 3. (Global asymptotic stability of the disease-free equilibrium point). Assume that [0, 1]; ∀i ∈ Z0+ with {Ki }i∞=0 → K with K = Kd f = Kend and that lim inf Ki Si − LI R̂i ≥ 0. Then, i→∞ the disease-free equilibrium point has no non-trivial limit periodic oscillation, and it is globally asymptotically stable for any bounded initial conditions in the first orthant of the phase space if β ∈ [0 , βc ). Proof. Note from Proposition 1 that Sd f /Rd f = A/K, which implies in (1) that β ∈ [0 , β c ), {Ki }i∞=0 → K and lim inf Ki Si − LI R̂i ≥ 0 implies that {Si }i∞=0 → Sd f since {Si }i∞=0  i→∞ is non-negative and asymptotically non-increasing, so that it cannot tend asymptoti- cally to a non-trivial periodic oscillation, irrespective of the finite initial conditions in the first orthant of the state space, and then lim inf (Si+1 − Si ) ≤ 0. Note that to con- i→∞ clude this global property to hold for the susceptible, it is not considered to be global the local property of Theorem 2 that lim Ei = lim Ii = 0 for small initial conditions i→∞ i→∞ in a neighborhood of the disease-free equilibrium point. Now, assume that there is a non-trivial limit oscillation of some period p(≥ 1) in the infectious subpopulation of val- ues lim Inp+i = lim I(n− j) p+i = lim I(n− j) p+i = Îi ; ∀i ∈ {0 , 1, · · · , p − 1}. n→∞ n,j(≤n−1)→∞ n,`(≤n−1)→∞  Thus, one gets for (3) that  np+i −1 µ lim  Inp+i −  ∑ (1 − γ)np+i−1− j Ej  = 0; ∀i ∈ {0 , 1, · · · , p − 1} (36) n →∞ ,1≤n−`<∞ 1 − (1 − γ ) ` p j=(n−`) p+i Mathematics 2021, 9, 520 14 of 32    np+i −1 µ lim  Inp+i −  ∑ (1 − γ)np+i−1− j Ej  = 0 (37) n,`(≤n−1)→∞ 1 − (1 − γ ) ` p j=(n−`) p+i irrespective of ` with 1 ≤ n − ` < ∞ as n → ∞ and also as n, `(≤ n − 1) → ∞ . As a result, lim Enp+i = lim E(n−`) p+i = Êi ; ∀i ∈ {0 , 1, · · · , p − 1} concluding that if n→∞ n,`(≤n−1)→∞ { Ii }i∞=0has a non-trivial limit oscillation of any period of p samples, then { Ei }i∞=0 has also a non-trivial limit oscillation of the same period. Until now, it has been proven that {Si }i∞=0 has no non-trivial limit oscillation, and that if { Ii }i∞=0 has a non-trivial limit oscillation of any period, then { Ei }i∞=0 has also a non-trivial oscillation of the same period. By using the same argument, it might be proven from (2) that if { Ei }i∞=0 has a non-trivial limit oscillation of any period, then { Ii }i∞=0 has also a non-trivial limit oscillation of the same period since {Si }i∞=0 has no limit non-trivial oscillation, and it has a limit. Until now, we have proven that {Si }i∞=0 has no limit non-trivial oscillation and {Si }i∞=0 → Sd f and that { Ii }i∞=0 and { Ei }i∞=0 either (a) they have both jointly non-trivial limit oscillations of some period, or (b) they do not have both such jointly non-trivial limit oscillations, while they converge both asymptotically to zero as time tends to infinity for any given non-negative and not all zero finite initial conditions. In Case a, {Si }i∞=0 → Sd f = 0 from (1), but by inspecting (2) and (3) with {Si }i∞=0 → Sd f = 0 , one concludes that, for any given strictly decreasing sequence {ε i }i∞=0 ⊂ R+ , there is a strictly increasing sequence { Ni (ε i )}i∞=0 ⊂ Z0+ such that S Ni +k ≤ ε i ; ∀k ∈ Z0+ and one gets from (2) that N + k −1 Ni +k −1− j ENi +k ≤ (1 − µ)k ENi + ∑ j=i N βIj + βe Ej Ij  i (1 − µ ) Ni +k −1− j i + k −1 ≤ (1 − µ)k ENi + ε i β(1 + λe )∑ N  j= N (1 − µ )i max Ej , Ij Ni ≤ j≥ Ni +k −1 Ni +k −1− j (38) i + k −1 ≤ (1 − µ)k ENi + ε i β(1 + λe )∑ N j= N (1 − µ ) i j = (1 − µ)k ENi + ε i β(1 + λe )∑kj=−01 (1 − µ) k = (1 − µ)k ENi + ε i β(1 + λe ) 1−(1µ−µ) < ENi   since ε i << 1 so that lim ENi +k = 0 and lim sup ENi +k − (1 − µ)k ENi ≤ 0. Furthermore, i,k →∞ i→∞ 1−(1−µ)k ENi +k ≤ (1 − µ)k ENi + 1− µ ε i β (1 + λ e ) < ENi ; k = 0, 1, . . . , Ni+1 − Ni − 1; (39) ∀i ∈ Z0+ ! µENi (1−µ)µEN provided that ε i < β (1+ λ e ) ; ∀i ∈ Z0+ and, under the choice ε i < min ε i−1 ,  N +1  i e , 1 − µ i +1 (1+ βλ ) ∞ then EN0 +k k=0 → 0 . Thus, { Ei }i=0 → 0 . In the same way, { Ii }i=0 → 0 so that { Ri }i∞=0 ∞ ∞  → Rd f = 1. In Case b, one has {Si }i∞=0 → Sd f ≥ 0 , { Ei }i∞=0 → 0 and { Ii }i∞=0 → 0 so that { Ri }i∞=0 → Rd f = 1 − Sd f . In both Case a and Case b there are no non-trivial oscillations of the solution for any given non-negative finite initial conditions, and the proof is complete.  If the local asymptotic stability condition of Theorem 2, being also global asymptotic according to Theorem 3 in the bounded first orthant, is changed under an alternative simpler proof to a more restrictive one under a smaller critical transmission rate, then the disease-free equilibrium point is ensured to be globally asymptotically stable. An alternative related sufficiency-type condition for it is discussed in the subsequent result. Mathematics 2021, 9, 520 15 of 32 Theorem 4. (Stronger conditions for global asymptotic stability of the disease-free equilibrium point). The disease-free equilibrium point is globally asymptotically stable if β < βcg = sup βcg (α), for α∈(1,∞)   (1−α)µ some real constant α ∈ (0 , 1] and where βcg (α) = min γα , λe . Proof. Take any parameterized Lyapunov sequence candidate {V (α)i }i∞=0 on the infec- tive subsystem defined by Vi (α) = Ei + αIi for the parameter α ∈ (0 , 1]. Then, since {Si }i∞=0 ⊂ [0 , 1], one gets that ( βIi + βe Ei )Si ≤ βIi + βe Ei ; ∀i ∈ Z0+ so that ∆Vi (α) = Vi+1 (α) − Vi (α) = −(αγ − βSi ) Ii − (µ(1 − α) − βλe Si ) Ei (40) ≤ −(αγ − β) Ii − (µ(1 − α) − βλe ) Ei ≤ 0 with the above inequality being strict if Ei + Ii > 0 provided that β < β cg (α) for any given α > 1. The proof is completed since it suffices that β < β cg (α) holds for some α ∈ (0 , 1] to ensure that the candidate {Vi (α)}i∞=0 is a Lyapunov sequence, which ver- ifies {Vi (α)}i∞=0 → 0 , so that { Ei }i∞=0 → 0 and { Ii }i∞=0 → 0 , irrespective of the given fi- nite non-negative initial conditions in (1)–(4). It has been proven that { Ei }i∞=0 → 0 and { Ii }i∞=0 → 0 provided that β < βcg . The convergences of the sequences {Si }i∞=0 → Sd f and { Ri }i∞=0 → Rd f follow under the same condition using similar arguments to those of the second part of the proof of Theorem 2 on the uniqueness of the disease-free equilibrium point, the lack of complex eigenvalues in the dynamics of a formally similar auxiliary system to (32) by replacing for all samples ∆Ei → Ei , ∆Ii → Ii and the convergence to a limit of its forcing term as a result of the convergence to zero of the infective subpopulations.  Remark 7. In practice, Theorem 4 can be checked for convenient design values of the test param- eter α ∈ (0 , 1] giving sufficiently large critical values β cg (α) or critical values easily testable. For instance: (1) The choice α = 1 − λe . gives β < min((1 − λe )γ , µ) provided that λe ∈ (0 , 1).   1−µγ (2) If µγ ≤ 1, then the choice α = µγ yields β < µmin γ2 , λe . βe µ−µ If the assumption for non-negativity of the solution sequence λe = β > γ 0 > γ (Theorem µ (3) 1) for some real constant µ0 ∈ (0 , µ) is incorporated then β < γmin(α , 1 − α) and, under the choice α = 1 − λe , Theorem 4 is satisfied for β < γmin(λe , 1 − λe ) if λe ∈ (0 , 1). 4.2. Stability around the Endemic Equilibrium Point The subsequent result addresses the global asymptotic stability of the endemic equilib- rium point excluding the presence of oscillations when the disease-free equilibrium point is unstable. ∈ [0, 1]; Theorem 5. (Global stability of the endemic equilibrium point). Assume that Vi = Ki Si  ∀i ∈ Z0+ with {Ki }i∞=0 → K with K = Kd f = Kend and that lim inf Ki Si − LI R̂i ≥ 0. i→∞ Then, the endemic equilibrium point has no non-trivial limit periodic oscillation, and it is globally asymptotically stable for any bounded initial conditions in the first orthant of the phase space if β ∈ ( β c , ∞ ). Proof. From (1), since β ∈ ( β c , ∞) and {Ki }i∞=0 → K and lim inf Ki Si − LI R̂i  ≥ 0 i→∞ one gets that {Si }i∞=0 → Send ≥ 0 since {Si }i∞=0 is non-negative and asymptotically non- increasing, so that it cannot tend asymptotically to a non-trivial periodic oscillation, ir- respective of the finite initial conditions in the first orthant of the state space, and then lim inf (Si+1 − Si ) ≤ 0. By summing-up (1) and (2), one gets that i→∞  (Si+1 + Ei+1 ) − (Si + Ei ) = −Vi + LI R̂i − µEi ≤ 0; ∀i ∈ Z0+ (41) Mathematics 2021, 9, 520 16 of 32 and then {Si + Ei }i∞=0 is asymptotically non-increasing and {Si + Ei }i∞=0 → Send + Eend so that it cannot have an asymptotic limit oscillation. Since, furthermore, {Si }i∞=0 does not oscillate asymptotically, then { Ei }i∞=0 does not oscillate asymptotically either. Thus, { Ii }i∞=0 does not oscillate asymptotically since { Ei }i∞=0 does not oscillate asymptotically (see proof of Theorem 3). Since the total population is constant then { Ri }i∞=0 cannot have an asymptotic oscillation. As a result the endemic equilibrium point is globally asymptotically stable for any finite initial conditions in the first orthant since the disease-free one is unstable if β ∈ ( β c , ∞).  Remark 8. The absence of oscillations and the global stability of the endemic equilibrium point of Theorem 5 requiring K ∈ [0 , 1] and lim inf Ki Si − LI R̂i ≥ 0 needs not to overpass a i→∞ minimum value of the constant A affecting as a multiplicative factor the term describing the partial loss of immunity in the susceptible and recovered in view of the endemic values of the equilibrium steady state. Note from (18) that the limit vaccination effort satisfies Vend = Kend Send ≤ 1 if β(µ+λe γ) K≤ γµ . On the other hand, the constraint requires that Kend Send ≥ ARend to guarantee  lim inf Ki Si − LI R̂i ≥ 0. Then, one gets from (18) and (21) that i→∞ Aµ( β − γ) + Aγλe β A[ β(µ + γλe ) − µγ]   A(γ + µ) K 1− ≥ = (42) γ(µ + A) + µA γ(µ + A) + µA γ(µ + A) + µA since the left-hand-side is non-negative. One gets from (42) that [ β(µ + γλe ) − µγ] A   K ≥ min ,0 (43) γµ which, together with the former non-negativity constraint of the vaccination gain yields [ β(µ + γλe ) − µγ] A β(µ + λe γ)   min ,0 ≤K≤ (44) γµ γµ µγ Note that the left hand-side constraint requires β ≥ β cv f = µ+γλe , which is the critical vaccination-free transmission rate associated with the basic reproduction number in the absence of vaccination effort. 5. Numerical Examples This section contains some numerical simulation examples aimed at illustrating the theoretical results discussed in the previous sections. The simulation results are divided into two examples. The first one is a numerical theoretical example, while the second one is related to the COVID-19 pandemic. Each example is split into different cases in order to show the effect of different parameter values in the model dynamics. 5.1. Example 1. Numerical Theoretical Example The parameters of the model in the absence of vaccination and loss of immunity are given by β 1 1 β = 0.5, βe = , µ = , γ = 3 6 2 so that λe = 13 . The sampling time is one day so that the units of the above parame- ters are in days accordingly. The initial conditions are given by S0 = 0.5 , E0 = 0.3 , I0 = 0.1 , R0 = 0.1. It can be noticed that the total population has been normalized to unity, Sk = Ek = Ik + Rk = S0 + E0 + I0 = R0 = 1, for all k ≥ 0 without loss of generality. The simulation example is split into a number of cases. The above set of parameters and initial values are used in all the cases except otherwise indicated. Mathematics 2021, 9, 520 19 of 34 Mathematics 2021, 9, 520 17 of 32 simulation example is split into a number of cases. The above set of parameters and ini- tial values are used in all the cases except otherwise indicated. Case Case 1.1.1.1. There There is no is no vaccination vaccination nornor lossloss of immunity. of immunity. It can It can be be readily readily verified verified that that the defined model satisfies Theorem 1 conditions, the defined model satisfies Theorem 1 conditions, since since β β γ γ ββ μµ∈∈( 0(0, 1, ]1, ],γγ∈∈(0(0, 1,] ,1], e ≤ ≤ , , e >> 00 (45)(45) β βe μ µ β βe and and μ μ μ 0.5 = S 0 ≤ µ e = 0.5 , 0.3 = E0 ≤ µ e = 0.5 , 0.1 = I 0 ≤ µ = 0.17 , R0 ≥ 0 (46) 0.5 = S0 ≤ 2βe = 0.5, 0.3 = E0 ≤ 2βe = 0.5, 0.1 = I0 ≤ 2 β = 0.17, R0 ≥ 0 (46) 2β 2β 2β As a result of Theorem 1, all the sequences of populations are non-negative and boundedAsfora result of Theorem all samples 1, all as Figure the sequences 1 shows. of populations It is observed in Figureare non-negative 1 that and the system bounded reaches for all samples a disease-free as Figure equilibrium point1 shows. It is observed with a nonzero numberin of Figure 1 that the susceptible andsystem im- reaches mune a disease-free equilibrium point with a nonzero number of susceptible and im- individuals. mune individuals. 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 10 20 30 40 50 60 time (days) Figure Figure 1. Sequences 1. Sequences of populations of populations for for CaseCase 1.1 1.1 conditions. conditions. Case Case 1.2.1.2. There There is vaccination is vaccination but but not not lossofofimmunity. loss immunity.TheThevaccination vaccinationgain gain is is set set to K = 0.2 and K = to Kddff = 0.2 and Kddff = 0.5 to show the effect of this control action on the dynamics ofthe 0.5 to show the effect of this control action on the dynamics of system. Since Theorem 1 conditions are still met, all the sequences of populations are non- thenegative system. and Sincebounded Theoremfor 1 conditions all time, as are it canstill bemet, seenallinthe sequences Figure of populations 2 for both vaccination aregains. non-negative and bounded Moreover, it is observedfor allintime, as it2can Figure thatbethe seen in Figure model 2 for both converges vaccination to the disease-free gains. equilibrium point given by the coordinates Sd f = Ed f = Id f = 0, Rd f = 1 in both cases, as Proposition 1 predicts. It is also observed from Figures 1 and 2 that the effect of K d f = 0.2 vaccination 1 when there is no loss of immunity is to make the whole population immune S asymptotically 0.9 while accelerating the convergenceE to the disease-free equilibrium point. I The effect of changing the vaccination gain is slight; R the susceptible vanish at a higher pace 0.8 for larger values of the gain, while the shape of the exposed and infectious does not change 0.7 significantly with vaccination. Figure 3 displays the vaccination function for both gains. As0.6expected, a larger vaccination gain results in a higher vaccination control action. Populations 0.5 Furthermore, according to Theorem 3, the obtained disease-free equilibrium points are globally 0.4 asymptotically stable. Figure 4 shows the dynamics of the system for perturbed initial conditions and Kd f = 0.5. It is seen in Figure 4 how all the system trajectories 0.3 converge to the globally asymptotically stable disease-free equilibrium point regardless of the0.2initial values. d f = 0.1, and an exponentially fast partial 0.1 Case 1.3. The vaccination gain is set to K loss0 of immunity with delay is considered. The parameters of the lost are given by d1 = 10 0 10 20 30 40 50 60 and d2 = 5 days and αtime=(days) 0.001 with A = 0.001 > 0. The dynamics of the system are depicted in Figure 5. Moreover, in addition to Conditions (1) and (2) from Theorem 1 we now need  to check that Constraint (3) holds. This condition can be equivalently written as LI R̂i ≤ Ki Si = Vi . Figure 6 shows that this constraint is also met in this example so that 0 0 10 20 30 40 50 60 time (days) Figure 1. Sequences of populations for Case 1.1 conditions. Mathematics 2021, 9, 520 Case 1.2. There is vaccination but not loss of immunity. The vaccination gain is18set of 32 to K df = 0.2 and K df = 0.5 to show the effect of this control action on the dynamics of the system. Since Theorem 1 conditions are still met, all the sequences of populations are all the population non-negative sequences and bounded aretime, for all guaranteed as it cantobe remain seen innon-negative Figure 2 for and bothbounded at all vaccination times, as is shown in Figure 5. gains. K d f = 0.2 1 S 0.9 E I 0.8 R 0.7 Mathematics 2021, 9, 520 0.6 20 of 34 Populations 0.5 0.4 K d f = 0.5 0.3 1 S 0.2 0.9 E I 0.1 0.8 R Mathematics 2021, 9, 520 20 of 34 0.7 0 0 10 20 30 40 50 60 0.6 time (days) Populations 0.5 K d f = 0.5 1 0.4 S 0.9 E 0.3 I 0.8 R 0.2 0.7 0.1 0.6 Populations 0 0 10 20 30 40 50 60 0.5 time (days) 0.4 Figure 0.3 2. Sequences of populations for Case 1.2 conditions. The vaccination gains K df = 0.2 and K df 0.2= 0.5 are employed in this simulation. 0.1 Moreover, 0 10 it is20 observed 0 30 in 40Figure502 that 60the model converges to the disease-free equi- librium point given by time coordinates S df = E df = I df = 0 , Rdf = 1 in both cases, as Prop- the(days) df =d f0.2 and Figure osition Figure 2. 12. Sequences predicts. Sequences ofispopulations ofItpopulations for for also observedCase Case from 1.2 conditions. Figures 1.2 conditions. andThe 1The vaccination 2 that gains thegains vaccination Kof effect = 0.2 and Kvaccination Kd f =there when 0.5 areisemployed no loss in immunity of this simulation. is to make the whole population immune asymptoti- K df = 0.5 are employed in this simulation. cally while accelerating the convergence to the disease-free equilibrium point. The effect of changing Furthermore, the system the vaccination gain istrajectories slight; the converge susceptible to vanish the disease-free at a higher equilibrium point pace for larger Moreover, it is observed in Figure 2the that the model K given values bythe of = Sd fgain,A+ A Kd f =the while 0.0099, shape Edoff = f = 0 and Idexposed and Rdconverges infectious d f f = A+Kd f does to=the disease-free 0.9901, not equi- in accordance change signifi- librium cantly point given by with vaccination. with Proposition the 1, which coordinates Figure 3 displays is globally S df = E df asymptotically= I df the vaccination = 0 , stable R df = function 1 according in both cases, for both as Prop- gains. As to Theorem ex- 3 since osition pected, β = 0.5 1apredicts. larger < It is also β c =vaccination 25.225 from observed gain results (30). fromin aFigures 1 and 2 that higher vaccination the effect control of vaccination action. when there is no loss of immunity is to make the whole population immune asymptoti- cally while 0.1 accelerating theK dconvergence f = 0.2 to the disease-free equilibrium point. The effect of changing the vaccination gain is slight; the susceptible vanish at a higher pace for larger 0.09 values of the gain, while the shape of the exposed and infectious does not change signifi- cantly0.08with vaccination. Figure 3 displays the vaccination function for both gains. As ex- 0.07 pected, a larger vaccination gain results in a higher vaccination control action. Vaccination, V i 0.06 K d f = 0.2 0.05 0.1 0.04 0.09 0.03 0.08 0.02 0.07 0.01 Vaccination, V i 0.06 0 0.05 0 10 20 30 40 50 60 time (days) 0.04 Figure 0.03 3. Cont. 0.02 0.01 0 0 10 20 30 40 50 60 time (days) 0.25 0.2 Mathematics 2021, 9, 520 19 of 32 Vaccination, V i 0.15 Mathematics 2021, 9, 520 21 of 34 0.1 K d f = 0.5 0.25 0.05 0.2 0 0 10 20 30 40 50 60 time (days) Vaccination, V i 0.15 Figure 3. Vaccination control action for K df = 0.2 and K df = 0.5 . 0.1 Furthermore, according to Theorem 3, the obtained disease-free equilibrium points are 0.05 globally asymptotically stable. Figure 4 shows the dynamics of the system for per- turbed initial conditions and K df = 0.5 . It is seen in Figure 4 how all the system trajecto- 0 ries converge 0 10 to the 20 globally 30 time (days) asymptotically 40 50 60 stable disease-free equilibrium point re- gardless of the initial values. Vaccinationcontrol Figure3.3.Vaccination Figure controlaction actionfor Kddff == 0.2 for K 0.2and andKd K df = 0.5 . f = 0.5. 0.8 Furthermore, according to Theorem 0.4 3, the obtained disease-free equilibrium points are globally asymptotically stable. Figure 4 shows the dynamics of the system for per- turbed initial conditions and K df = 0.5 . It is seen in Figure 4 how all the system trajecto- 0.6 0.3 ries converge to the globally asymptotically stable disease-free equilibrium point re- gardless of the initial values. 0.4 0.2 E S 0.2 0.1 0.8 0.4 0 0 0.6 0 10 20 30 0.3 0 10 20 30 time (days) time (days) 0.4 0.2 E S 0.2 1 0.2 0.1 0.8 0.15 0 0 0 10 20 30 0.6 0 10 20 30 0.1 R I time (days) time (days) 0.4 0.2 0.05 1 0.2 0.8 0.15 0 0 0 10 20 30 0 10 20 30 0.6 time (days) time (days) 0.1 R I 0.4 Sequencesof Figure4.4.Sequences ofpopulations populations for Case 1.2 conditions, f = 0.5 KdK and perturbed initial values. Figure 0.05 for Case 1.2 conditions, df = 0.5 and perturbed initial values. 0.2 Case 1.4. The vaccination gain is set to Kd f = 0.1, and a constant loss of immunity 0 is setofto df =are theKlost 0.1given 0 with Case delay1.3. The vaccination is considered. gain The parameters , andby andexponentially 1 = 4 and d2 =fast partial 2 days 0 10 20 30 0 10 20 30 and a loss of = 0.0005 immunity time (days) with A = 0.0015. The dynamics of the with delay is considered. The parameters system are of the lost are given by d1it= 10 depicted in Figure 7. As happened before, in addition to Conditions (1)time and(days) (2) from Theorem 1 we need to check and d 2 = 5 days that Constraint and α = 0.001 with A = 0.001 > 0 . The dynamics of the system are (3) holds, which is performed in the same way as done for Case 1.3. Figure 8 Figure 4. Sequences of populations depicts for Case 1.2 how Condition conditions, 3 from Theorem K df 1=is0.met. 5 and Furthermore, perturbed initialthevalues. system trajectories A converge to the disease-free equilibrium point Sd f = A+K = 0.0566, Ed f = Id f = 0 and df Case K1.3. df The vaccination gain is set to K df = 0.1 , and an exponentially fast partial Rd f = A+K = 0.9434, which is globally asymptotically stable according to Theorem 3 loss of immunity df with delay is considered. The parameters of the lost are given by d1 = 10 since β = 0.5 < β c = 3.5 from (30). and d 2 = 5 days and α = 0.001 with A = 0.001 > 0 . The dynamics of the system are depicted in Figure 5. Moreover, in addition to Conditions (1) and (2) from Theorem 1 we depicted in Figure 5. Moreover, in addition to Conditions (1) and (2) from Theorem 1 we now need to check that Constraint (3) holds. This condition can be equivalently written (( )) now need to check that Constraint (3) holds. This condition can be equivalently written as LI R̂i ≤ K i S i = Vi . Figure 6 shows that this constraint is also met in this example so as LI R̂i ≤ K i S i = Vi . Figure 6 shows that this constraint is also met in this example so Mathematics 2021, 9, 520 that all the population sequences are guaranteed to remain non-negative and bounded 20 of 32 that all the population sequences are guaranteed to remain non-negative and bounded at all times, as is shown in Figure 5. at all times, as is shown in Figure 5. 1 1 S 0.9 SE 0.9 EI 0.8 IR 0.8 R 0.7 0.7 0.6 Populations 0.6 Populations 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 5 10 15 20 25 30 35 40 45 0 0 5 10 15 time 20 (days) 25 30 35 40 45 time (days) Figure 5. Sequences of populations for Case 1.3 conditions. The vaccination gain K df = 0.1 is Figure Figure 5. Sequences 5. Sequences of populations of populations for Case for Case 1.3 conditions. 1.3 conditions. The vaccination The vaccination df = K gain Kgain 0.d1f is = 0.1 is employed employed in this in simulation this with simulation exponentially with fast exponentially partial fast loss partial of of loss immunity with immunity delay. with delay. employed in this simulation with exponentially fast partial loss of immunity with delay. 0.05 0.05 LI 0.045 V=K *S LI df 0.045 V=K *S df 0.04 0.04 0.035 0.035 0.03 0.03 0.025 0.025 0.02 0.02 0.015 Mathematics 2021, 9, 520 23 of 34 0.015 0.01 0.01 0.005 0.005 0 K df Rdf = 00 5 = 0. 10 943415 15 20 , which 25 is 30 35 globally 40 45 asymptotically stable according to Theorem 3 0 A + K df5 10 time 20 (days) 25 time (days) 30 35 40 45 since β6.=Satisfaction Figure Figure 5 < β c = of 6.0.Satisfaction 3.5Condition of from (30). Condition 3 from 3 from Theorem Theorem 1 in 1 in Case Case 1.31.3 conditions. conditions. Figure 6. Satisfaction of Condition 3 from Theorem 1 in Case 1.3 conditions. Furthermore, 1 the system trajectories converge to the disease-free equilibrium point Furthermore, the system trajectories converge to the disease-free K df equilibrium point A S given0.9 by S df = A = 0. 0099 , E = I = 0E and R df = K df = 0.9901 , in accord- given by S df = A + K df = 0.0099 , E dfdf = I dfdf = 0IR and Rdf = A + K df = 0.9901 , in accord- 0.8 A + K df A + K df ance0.7with Proposition 1, which is globally asymptotically stable according to Theorem 3 ance with Proposition 1, which is globally asymptotically stable according to Theorem 3 since0.6 β = 0.5 < β c = 25.225 from (30). since β = 0.5 < β c = 25.225 from (30). Populations Case 1.4. The vaccination gain is set to K df = 0.1 , and a constant loss of immunity with 0.5 Case 1.4. The vaccination gain is set to K df = 0.1 , and a constant loss of immunity with delay is considered. The parameters of the lost are given by d1 = 4 and d 2 = 2 days and 0.4 delay0.3is considered. The parameters of the lost are given by d1 = 4 and d 2 = 2 days and a = 0.0005 with A = 0.0015. The dynamics of the system are depicted in Figure 7. As it hap- a = 0. 0.20005 with A = 0.0015. The dynamics of the system are depicted in Figure 7. As it hap- pened before, in addition to Conditions (1) and (2) from Theorem 1 we need to check that pened before, in addition to Conditions (1) and (2) from Theorem 1 we need to check that Constraint 0.1 (3) holds, which is performed in the same way as done for Case 1.3. Figure 8 de- Constraint (3) holds, which is performed in the same way as done for Case 1.3. Figure 8 de- picts 0how 0 Condition 10 20 3 from30Theorem 40 1 is met. Furthermore, the system trajectories con- picts how Condition 3 from Theorem time (days) 1 is50met. Furthermore, 60 A the system trajectories con- verge to the disease-free equilibrium point S df = A = 0.0566 , E df = I df = 0 and verge 7.toSequences theSequences disease-free equilibrium point = AThe = 0.vaccination + Kvaccination S dfconditions. 0566gain = I K= 0=and , E dfK gain df =df0.1d f is Figure Figure 7. of populations of populations for 1.4 for Case Case 1.4 conditions. df The 0.1 is A + K df employed employed in this in this simulation simulation with with constant constant lossloss of immunity of immunity withwith delay. delay. 0.05 LI 0.045 V=K *S df 0.04 0.035 0.03 0.025 0.02 0.015 0.2 0.1 0 0 10 20 30 40 50 60 time (days) Mathematics 2021, 9, 520 Figure 7. Sequences of populations for Case 1.4 conditions. The vaccination gain K df = 0.1 is 21 of 32 employed in this simulation with constant loss of immunity with delay. 0.05 LI 0.045 V=K *S df 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0 0 5 10 15 20 25 30 35 40 45 time (days) Figure Figure 8. Satisfaction 8. Satisfaction of Condition of Condition 3 from 3 from Theorem Theorem 1 in1 Case in Case 1.41.4 conditions. conditions. CaseCase1.5.1.5. There There is noisvaccination no vaccination nor of nor loss loss of immunity, immunity, but the butinitial the initial populations populations do not meet Theorem 1 conditions. This case is aimed at showing that the Theorem Theorem do not meet Theorem 1 conditions. This case is aimed at showing that the 1 condi- 1 conditions tions are sufficient are sufficient to guarantee to guarantee the boundedness the boundedness and positivity and positivity of the populations of the populations but are but are not necessary. Thus, the following initial conditions are used in this case S0 = 0.99, not necessary. Thus, the following initial conditions µare used in this case S 0 = 0.99 , E0 = 0.01, I0 = 0 and R0 = 0 for which 0.99 = S0 > 2βe = 0.5, not satisfying Condition (2). μ Mathematics 2021, 9, 520 01 , I 0 = 0 ofand = 0.dynamics E0The 0 = 0 for theRsystem arewhich 0.99in= Figure displayed S 0 > 9.e As = 0it.5can , not satisfying inCondition be observed Figure 9 all 24 of 34 2 β populations remain nonnegative and bounded for all time. In addition, the disease-free (2).equilibrium The dynamics of the point system is also are displayed globally in Figure asymptotically 9. As stable as itTheorem can be observed 3 states. in Figure 9 all populations remain nonnegative and bounded for all time. In addition, the disease-free equilibrium 1 point is also globally asymptotically stable as Theorem 3 states. 0.9 0.8 0.7 S E I 0.6 R Populations 0.5 0.4 0.3 0.2 0.1 0 0 50 100 150 time (days) Figure Figure 9. Sequences 9. Sequences of populations of populations for for Case Case 1.5 1.5 conditions. conditions. Case Case 1.6.1.6. ThereThere is no is no vaccination, vaccination, and a constant and a constant loss of immunity loss of immunity with with delay delay is is consid- considered. ered. The parameters The parameters of the lost of arethe lostby given are dgiven = by, dd1 == 180, 180 10 d2 =and days 10 days a = and 0. a =with 0007 0.0007 1 2 with A = 0.0077. The value of β is increased to β = 2 so that the endemic equilibrium point A = 0.0077. The value of β is increased to β = 2 so that the endemic equilibrium point is is now reachable and stable since β = 2 > β c = 0.25. The system trajectories are depicted now reachable and stable since β = 2 > β c = 0.25 . The system trajectories are depicted in in Figure 10. Figure 10. It is observed in Figure 10 how after 180 days the loss of immunity makes some former immune individuals become susceptible again. In the end, an endemic equilibrium point is reached. 1 From the numerical simulation we can obtain the location of the endemic equilibrium 0.9 point and compare it with the theoretical values given in Proposition 2. The comparison is performed in Table 1, which shows agreement between them. 0.8 S 0.7 E I 0.6 R Populations 0.5 0.4 0.3 0.2 0.8 0.7 S E I 0.6 R Populations 0.5 Mathematics 2021, 9, 520 22 of 32 0.4 0.3 0.2 Table 1. Comparison of the theoretical values predicted by Proposition 2 and the numerically obtained 0.1 ones for the location of the endemic point. 0 0 50 100 Theoretical 150 Numerical time (days) Send 0.1250 0.1249 Figure 9. Sequences Eendof populations for Case 1.5 conditions. 0.0381 0.0380 Iend 0.0127 0.0127 Rend is no vaccination, and a constant Case 1.6. There 0.8242 loss of immunity with delay 0.8243 is consid- ered. The parameters of the lost are given by d1 = 180 , d 2 = 10 days and a = 0.0007 with A = 0.0077. Furthermore, of β is increased The value according to 5βthe to Theorem = 2endemic so that equilibrium the endemicpoint equilibrium point is globally is asymp- nowtotically stable, reachable andas Figure stable since β = 2 >for 11 shows β c different = 0.25 . The initial conditions. system Moreover, trajectories Figure are depicted in 10 shows Figure 10. that the convergence to the endemic point exhibits a transient oscillation, but the trajectory does not have a steady-state periodic oscillation. 1 0.9 0.8 S 0.7 E I 0.6 R Populations 0.5 0.4 0.3 0.2 0.1 0 0 500 1000 1500 time (days) Figure Figure 10. 10. Sequences Sequences of populations of populations for for Case Case 1.6 1.6 conditions. conditions. Case It is 1.7. The observed in vaccination Figure 10 how is Kd180 gainafter f = days 0.1, and the aloss constant loss of immunity of immunity makes some with delay is considered. The parameters of the loss of immunity former immune individuals become susceptible again. In the end, an endemic are given by d 1 = 180, d equilib- 2 = 10 days and a = 0.0007 with A = 0.0077. The value of β = 2 is used rium point is reached. From the numerical simulation we can obtain the location of thenow as well. The system trajectories endemic are depicted equilibrium in Figure point and compare 12. itItwith is observed in Figure the theoretical 12 how values given the in vaccination Proposi- changes tion the dynamics 2. The comparison of the system is performed with1, respect in Table to theagreement which shows vaccination-free between case. them. Thus, vaccination is able to lessen the impact of loss of immunity so that the endemic equilibrium point Table is avoided and 1. Comparison the of the system converges theoretical and remains values predicted in the disease-free by Proposition equilibrium 2 and the numerically ob- one. Therefore, tained ones forthe theimplementation of adequate location of the endemic point. vaccination campaigns throughout the year may act as a powerful tool to counteract the possibility of loss of immunity. This may be the case for the COVID-19 pandemic. Thus, Example Theoretical2 is devoted to study the Numerical application of the proposed Smodelend 0.1250 to COVID-19. From a mathematical 0.1249 point of view, the effect of vaccination is to changeE end the critical beta value so that Condition 0.0381 (30) becomes β 0.0380> β c = 3.5. = 2 Consequently,I end the disease-free equilibrium point is globally asymptotically 0.0127 0.0127stable while the endemic Rend one is not reachable. Figure 13 displays the vaccination action to be applied 0.8242 0.8243 through time. It is seen that vaccination does not vanish now asymptotically since it has to cope with the new susceptible appearingdue to the loss of immunity. Finally, Theorem 3 requires the condition inf Ki Si − LI R̂i ≥ 0 asymptotically in order to guarantee the global stability. This condition is shown in Figure 14. As it is observed in Figure 14 this condition may be violated within some time intervals, but it is achieved asymptotically. Thus, all conditions of Theorem 3 are achieved, and the disease-free equilibrium point is asymptotically stable. Furthermore, according to Theorem 5 the endemic equilibrium point is globally as- ymptotically stable, as Figure 11 shows for different initial conditions. Moreover, Figure 10 Mathematics 2021, 9, 520 shows that the convergence to the endemic point exhibits a transient oscillation, but 23 of the 32 tra- jectory does not have a steady-state periodic oscillation. 1 0.9 0.8 0.7 0.6 0.5 S 0.4 0.3 0.2 0.1 0 0 500 1000 1500 time (days) 0.7 0.6 0.5 0.4 E 0.3 0.2 0.1 0 0 500 1000 1500 time (days) 0.2 0.18 0.16 0.14 0.12 0.1 I 0.08 0.06 0.04 0.02 0 0 500 1000 1500 time (days) Figure 11. Cont. Case 1.7. The vaccination gain is K df = 0.1 , and a constant loss of immunity with de- lay is considered. The parameters of the loss of immunity are given by d1 = 180 , d 2 = 10 days and a = 0.0007 with A = 0.0077. The value of β = 2 is used now as well. The sys- Mathematics 2021, 9, 520 26 of 34 Mathematics 2021, 9, 520 tem trajectories are depicted in Figure 12. It is observed in Figure 12 how the vaccination 24 of 32 changes the dynamics of the system with respect to the vaccination-free case. Thus, vac- cination is able to lessen the impact of loss of immunity so that the endemic equilibrium point is1 avoided and the system converges and remains in the disease-free equilibrium one. Therefore, 0.9 the implementation of adequate vaccination campaigns throughout the year may act as a powerful tool to counteract the possibility of loss of immunity. This may be the case 0.8 for the COVID-19 pandemic. Thus, Example 2 is devoted to study the application of the 0.7 proposed model to COVID-19. From a mathematical point of view, the effect of vaccination is to change the critical beta value so that Condition (30) becomes β = 2 > β c = 3.5 . Conse- 0.6 quently, 0.5 the disease-free equilibrium point is globally asymptotically stable while the en- R demic0.4one is not reachable. Figure 13 displays the vaccination action to be applied through time. 0.3 It is seen that vaccination does not vanish now asymptotically since it has to cope with the new susceptible appearing due to the loss of immunity. Finally, Theorem 3 requires the ( ( )) ≥ 0 asymptotically in order to guarantee the global stabil- 0.2 condition 0.1 inf K i Si − LI Rˆi ity. This 0 condition is shown in Figure 14. As it is observed in Figure 14 this condition may be 0 500 1000 1500 violated within some time intervals, time (days) but it is achieved asymptotically. Thus, all conditions of Theorem 3 are Figure 11. achieved, Global and the asymptotical disease-free stability equilibrium of the endemic point equilibrium is asymptotically point stable. for Case 1.6 conditions. Figure 11. Global asymptotical stability of the endemic equilibrium point for Case 1.6 conditions. Case 1 1.7. The vaccination gain is K df = 0.1 , and a constant loss of immunity with de- lay is0.9considered. The parameters of the loss of immunity are given by d1 = 180 , d 2 = 10 S and a = 0.0007 with A = 0.0077. The valueEI of β = 2 is used now as well. The sys- days 0.8 0.7 R tem trajectories are depicted in Figure 12. It is observed in Figure 12 how the vaccination changes 0.6 the dynamics of the system with respect to the vaccination-free case. Thus, vac- Populations cination 0.5 is able to lessen the impact of loss of immunity so that the endemic equilibrium point0.4is avoided and the system converges and remains in the disease-free equilibrium one. Therefore, 0.3 the implementation of adequate vaccination campaigns throughout the year may act as a powerful tool to counteract the possibility of loss of immunity. This may be the case 0.2 for the COVID-19 pandemic. Thus, Example 2 is devoted to study the application of the 0.1 Mathematics 2021, 9, 520 proposed model to COVID-19. From a mathematical point of view, the effect of vaccination 27 of 34 is to change 0 0 50 the100critical 150 beta 250 200 value 300 so350that400Condition 450 500 (30) becomes β = 2 > β c = 3.5 . Conse- time (days) quently, the disease-free equilibrium point is globally asymptotically stable while the en- demic Figure Figure one 12.is not reachable. Sequences Sequences Figurefor13 ofofpopulations populations for Casedisplays Case 1.7 1.7 the vaccination action to be applied through conditions. conditions. time. It is seen that vaccination does not vanish now asymptotically since it has to cope with the new susceptible appearing K d f =due 0.1 to the loss of immunity. Finally, Theorem 3 requires the ( ( )) ≥ 0 asymptotically in order to guarantee the global stabil- 0.05 condition 0.045 inf K i Si − LI Rˆi 0.04 ity. This condition is shown in Figure 14. As it is observed in Figure 14 this condition may be violated 0.035within some time intervals, but it is achieved asymptotically. Thus, all conditions of Theorem 0.03 3 are achieved, and the disease-free equilibrium point is asymptotically stable. Vaccination, V i 0.025 10.02 0.015 0.9 S 0.80.01 E I 0.005 0.7 R 0.6 0 Populations 0 50 100 150 200 250 300 350 400 450 500 time (days) 0.5 Figure Figure 13. Vaccination 0.4 13. Vaccinationaction actionforfor Case 1.71.7 Case conditions. conditions. 0.3 0.20.05 LI 0.045 0.1 V=K *S df 00.04 0 50 100 150 200 250 300 350 400 450 500 0.035 time (days) 0.005 0 0 50 100 150 200 250 300 350 400 450 500 time (days) Mathematics 2021, 9, 520 25 of 32 Figure 13. Vaccination action for Case 1.7 conditions. 0.05 LI 0.045 V=K *S df 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0 0 50 100 150 200 250 300 350 400 450 500 time (days) -3 10 9.5 LI V=K *S df 9 8.5 8 7.5 7 6.5 300 350 400 450 time (days) Figure Figure Condition to 14.14.Condition to guarantee guarantee thethe global asymptotic global stability asymptotic of the disease-free stability equilibrium of the disease-free equilibriu point for Case 1.7. point for Case 1.7. 5.2. Example 2. Application to the COVID-19 Pandemic 5.2. Example 2. Application This example is devotedto tothe COVID-19 study Pandemic the application of the proposed model to the COVID- 19 pandemic. To this end, the case of Italy will be considered. Thus, the model is parame- terized by β 1 β = 0.2 , βe = , µ = , γ = 0.04 1.3 6 e while λ = 0.7692. The sampling time is one day so that the units of the above parameters are in days−1 . The initial conditions are given by S0 = 240, 000, E0 = 0, I0 = 240 and R0 = 0. These populations have been normalized to unity in the simulation examples, without loss of generality. The number of infectious, however, will be presented with actual data. This example has been divided into a number of cases in order to illustrate the dynamics of the model under different circumstances. Case 2.1. There is no vaccination nor loss of immunity. Figure 15 displays the dynamics of the model along with the number of real cases reported in Italy from 26 February 2020 (which corresponds to the first day of simulation) for 145 consecutive days. This example for the COVID-19 pandemic has been considered previously in [31,32]. trate the dynamics of the model under different circumstances. Case 2.1. There is no vaccination nor loss of immunity. Figure 15 displays th namics of the model along with the number of real cases reported in Italy from 26 Mathematics 2021, 9, 520 ruary 2020(which corresponds to the first day of simulation) for 145 consecutive 26 of 32 This example for the COVID-19 pandemic has been considered previously in [31,32] 1 0.9 S 0.8 E I 0.7 R Populations 0.6 0.5 0.4 0.3 0.2 0.1 0 0 50 100 150 time (days) 4 10 12 I real cases 10 8 Populations 6 4 2 0 0 50 100 150 time (days) Figure Figure Dynamics of 15.Dynamics 15. ofthe themodel modelfor for the COVID-19 pandemic the COVID-19 in Italy for pandemic inCase 2.1. Italy for Case 2.1. In the absence of loss of immunity, all the population becomes immune asymptotically. TheInfollowing the absence of lossthe cases discuss ofeffects immunity, of loss ofall the population immunity becomes and vaccination on theimmune COVID- asym cally. The following cases discuss the effects of loss of immunity and vaccination o 19 dynamics. COVID-19 Case dynamics. 2.2. There is no vaccination, and an exponentially fast loss of immunity with delay is considered. The parameters of the loss of immunity are given by d1 = 180, d2 = 10 days and a = 0.009 with A = 0.0091. The dynamics of the models are depicted in Figure 16. In Figure 16, a substantial change in the system dynamics due to the loss of immunity is observed. Therefore, the asymptotically stable disease-free equilibrium point gives way to an endemic equilibrium point, whose location is given by Proposition 2 as Send = 0.1688, Eend = 0.0353, Iend = 0.1473 and Rend = 0.6486, in normalized units. Thus, in the steady- state there is always a number of infectious individuals. Moreover, the endemic equilibrium point is globally asymptotically stable since β = 0.2 > β c = 0.034 according to Theorem 5. Case 2.3. There is no vaccination, and a constant loss of immunity with delay is considered. The parameters of the loss of immunity are given by d1 = 180, d2 = 10 days and a = 0.0007 with A = 0.0077. The dynamics of the models are depicted in Figure 17. As it happened in Case 2.2, the dynamics of the system changed substantially, and an endemic equilibrium point appears. This situation is not recommendable at all, and human Mathematics 2021, 9, 520 27 of 32 Mathematics 2021, 9, 520 29 intervention is necessary to avoid it. Consequently, the next Case 2.4 will study the effect of vaccination when loss of immunity occurs. Case 2.4. The vaccination gain is Kd f = 0.1, and a constant loss of immunity with delay is considered. The parameters of the loss of immunity are given by d1 = 180, d2 = 10 Case 2.2. There is no vaccination, and an exponentially fast loss of immunity with days and a = 0.0007 with A = 0.0077. The vaccination is only applied 180 days after the is considered. first day. ThisThe parameters situation ofconsidered has been the loss ofinimmunity are given order to discuss the effect = 180 , d 2 = 10 by ofd1vaccination a = 0.009 the andto counteract with lossAof=immunity. 0.0091. The Thedynamics dynamics ofofthe thesystem modelsareare depicted depicted in Figure in Figure 18, 16. while the vaccination is displayed in Figure 19. 1 S 0.9 E I 0.8 R 0.7 0.6 Populations 0.5 0.4 0.3 0.2 0.1 0 0 500 1000 1500 2000 2500 time (days) 4 10 14 I real cases 12 10 Populations 8 6 4 2 0 0 500 1000 1500 2000 2500 time (days) Figure 16.16.Dynamics Figure Dynamics of of the modelforfor the model thethe COVID-19 COVID-19 pandemic pandemic in Italyinfor Italy Casefor 2.2.Case 2.2. There There is n is no vaccination,and vaccination, and there there isisanan exponentially exponentiallyfast loss fastofloss immunity with delay. of immunity with delay. It is deduced from Figure 18 that the vaccination, even with as low a percentage as 10%, In Figure allows avoiding 16,thea endemic substantial change state caused byinthethe losssystem dynamics of immunity duethetototal and makes thepop- loss of imm ity ulation is observed. be eitherTherefore, the susceptible or asymptotically immune asymptotically. stable disease-freeexplanation The mathematical equilibrium for point way thisto is that an vaccination endemic generates a new point, equilibrium critical value whose for beta so that βis=given location 0.2 < βby c = 0.4723. Proposition Therefore, the disease-free equilibrium point is globally asymptotically S end = 0.1688 , E end = 0.0353 , I end = 0.1473 and Rend = 0.6486 , in normalized stable, while the u endemic equilibrium point is no longer reachable. Overall, the proposed model is useful Thus, in the steady-state there is always a number of infectious individuals. Moreove endemic equilibrium point is globally asymptotically stable since β = 0.2 > β c = 0 according to Theorem 5. Case 2.3. There is no vaccination, and a constant loss of immunity with delay is co Mathematics 2021, 9, 520 28 of 32 ematics 2021, 9, 520 to describe the COVID-19 pandemic, study the effect of loss of immunity and assess the impact of vaccination in the epidemic spread. 1 S 0.9 E I 0.8 R 0.7 0.6 Populations 0.5 0.4 0.3 0.2 0.1 0 0 500 1000 1500 2000 2500 time (days) 4 10 14 I real cases 12 10 Populations 8 6 4 2 0 0 500 1000 1500 2000 2500 time (days) Figure 17. Dynamics of the model for the COVID-19 pandemic in Italy for Case 2.3. There is no Figure 17. Dynamics of the model for the COVID-19 pandemic in Italy for Case 2.3. vaccination, and there is a constant loss of immunity with delay. vaccination, and there is a constant loss of immunity with delay. As it happened in Case 2.2, the dynamics of the system changed subs an endemic equilibrium point appears. This situation is not recommenda human intervention is necessary to avoid it. Consequently, the next Case the effect of vaccination when loss of immunity occurs. Case 2.4. The vaccination gain is K df = 0.1 , and a constant loss of immun is considered. The parameters of the loss of immunity are given by d1 = 180 , and a = 0.0007 with A = 0.0077. The vaccination is only applied 180 days day. This situation has been considered in order to discuss the effect of Mathematics Mathematics 2021, 2021, 9, 9, 520 520 31 31 of of 34 34 Mathematics 2021, 9, 520 29 of 32 1 1 0.9 0.9 S 0.8 S 0.8 E E I I 0.7 R 0.7 R 0.6 Populations 0.6 Populations 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 00 100 200 300 400 500 600 0 100 200 300 400 500 600 time (days) time (days) 4 10 4 12 10 12 I I real cases real cases 10 10 8 8 Populations Populations 6 6 4 4 2 2 0 00 100 200 300 400 500 600 0 100 200 300 400 500 600 time (days) time (days) Figure 18. Figure Figure Model 18.18. Model dynamics Model dynamics dynamics with with vaccination with vaccination vaccination and aa constant andand loss a constant constant of of immunity lossloss with of immunity immunity with delay. with The delay. delay. vac- TheThe vac- cination gain is set to K = 0.1 and it is applied after 180 cination gain is set to Kdf =d f0.1 and it is applied after 180 days. vaccination gain is set to K = 0.1 and it is applied after days. 180 days. df -3 K d f = 0.1 10 -3 K d f = 0.1 8 10 8 7 7 6 6 5 Vaccination,VV i 5 i Vaccination, 4 4 3 3 2 2 1 1 0 00 100 200 300 400 500 600 0 100 200 300 400 500 600 time (days) time (days) Vaccination forCase Case2.4, 2.4,where whereaa constant constant loss Figure Figure Figure 19.19. 19. Vaccination Vaccination for for Case 2.4, where a constant lossof loss ofimmunity of immunitywith immunity delay with with is is delay delay considered. The is considered. considered. Thevaccination vaccinationgain is is gain setset to to =df0.1=and Kd f K 0.1 itand is applied it is after 180 applied days. after 180 days. The vaccination gain is set to K = 0.1 and it is applied after 180 days. df 6. 6. Conclusions Conclusions A A new new discrete discrete SEIR SEIR model model has has been been presented presented in in this this paper paper subject subject to to linear linear feed- feed- back back vaccination controls on the susceptible and delayed partial loss of immunity of vaccination controls on the susceptible and delayed partial loss of immunity of Mathematics 2021, 9, 520 30 of 32 6. Conclusions A new discrete SEIR model has been presented in this paper subject to linear feedback vaccination controls on the susceptible and delayed partial loss of immunity of previously recovered individuals. It has also been considered that there is infectivity power in the Exposed-Susceptible contacts being, in general, of distinct transmission rate than that related to the Infectious-Susceptible contacts. The above idea is based on the knowledge that the later periods in the asymptomatic incubation phase of the infection in some infectious diseases, like for instance COVID-19, are also infective. It is considered that the partial loss of immunity in the recovered subpopulation increases the susceptibility after a certain delay depending on the period of the recovery immunity maintenance. Such an immunity loss is considered to take place in a range of samples previous to the current sample being active in the model, and both the duration of such a time interval where immunity is lost and the particular force per sample of such a loss can be adjusted with different weights in the proposed model to include the effect that the immunity loss is usually gradual through time. The proposed model considers also the use of feedback vaccination on the susceptible subpopulation. Furthermore, both the disease-free and the endemic equilibrium points are calculated explicitly as being dependent on the vaccination control gains. The stability of both equilibrium points is also discussed. It is proven that when the disease-free equilibrium point is asymptotically stable, the endemic one is not attainable, and when the endemic one is reachable then the disease-free one is unstable. In this way, it is proven that there is only a global attractor depending on the reproduction number or, equivalently, on its associate critical transmission rate. The proposed model has been tested through numerical examples based on COVID-19 parameterizations. Author Contributions: M.D.l.S. contributed to the formal treatment of the material in the paper and wrote the manuscript. S.A.-Q. revised, checked and completed the proofs and supervised the article. A.I. prepared, programmed and interpreted the simulated examples. All authors have read and agreed to the published version of the manuscript. Funding: The work has been funded by Grant RTI2018-094336-B-I00 from MCIU/AEI/FEDER, UE; by Grant IT1207-19, by the Basque Government and by Grant COV 20/01213 from Spanish Institute of Health Carlos III. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Availability of Data and Material: The data used to support the findings of this study are included in the references within the article. Acknowledgments: The authors are grateful to the Spanish Government for Grant RTI2018-094336- B-I00 (MCIU/AEI/FEDER, UE), to the Basque Government for Grant IT1207-19 and to the Spanish Institute of Health Carlos III for its support through Grant COV 20/01213. The authors thank also the reviewers for their useful suggestions. Conflicts of Interest: The authors declare that they do not have any competing interests. 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.; O Wang, H. Stability Analysis and Bifurcation Control For a Fractional Order SIR Epidemic Model with Delay. In Proceedings of the 2020 39th Chinese Control Conference (CCC), Shenyang, China, 27–29 July 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] Mathematics 2021, 9, 520 31 of 32 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 2018 UKACC 12th International Conference on Control (CONTROL), Sheffield, UK, 5–7 September 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. [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. Equations 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–7044. [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. Informatics 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] 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; pp. 300–305. 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. Res. 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. [CrossRef] 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, 17. [CrossRef] 33. Wang, Y.; Wei, Z.; Cao, J. Epidemic dynamics of influenza-like diseases spreading in complex networks. Nonlinear Dyn. 2020, 101, 1–20. [CrossRef] 34. Wang, Y.; Cao, J. Final size of network epidemic models: Properties and connections. Sci. China Inf. Sci. 2020, 64, 1–3. [CrossRef] 35. 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. Mathematics 2021, 9, 520 32 of 32 36. Ng, K.Y.; Gui, M.M. COVID-19: Development of a robust mathematical model package with consideration for ageing population and time delay for control action and resusceptibility. Phys. D Nonlinear Phenom. 2010, 114, 132599. 37. Ng, M.K.Y. SEIRS-Based COVID-19 Simulation Packge. Available online: https://www.markusng.com/COVIDSiM/ (accessed on 10 February 2021). 38. Anderez, D.O.; Kanjo, E.; Pogrebna, G.; Kaiwartya, O.; Johnson, S.D.; Hunt, J.A. A COVID-19-Based Modified Epidemiological Model and Technological Approaches to Help Vulnerable Individuals Emerge from the Lockdown in the UK. Sensors 2020, 20, 4967. [CrossRef] [PubMed] 39. Dobbs, P.S.; Watts, D.J. A generalized model of social and biological contagion. J. Theor. Biol. 2005, 232, 587–604.