On a Discrete SEIR Epidemic Model with Two-Doses Delayed Feedback Vaccination Control on the Susceptible

Abstract

A new discrete susceptible-exposed-infectious-recovered (SEIR) epidemic model is presented subject to a feedback vaccination effort involving two doses. Both vaccination doses, which are subject to a non-necessarily identical effectiveness, are administrated by respecting a certain mutual delay interval, and their immunity effect is registered after a certain delay since the second dose. The delays and the efficacies of the doses are parameters, which can be fixed in the model for each concrete experimentation. The disease-free equilibrium point is characterized as well as its stability properties, while it is seen that no endemic equilibrium point exists. The exposed subpopulation is supposed to be infective eventually, under a distinct transmission rate of that of the infectious subpopulation. Some simulation examples are presented by using disease parameterizations of the COVID-19 pandemic under vaccination efforts requiring two doses.

1. Introduction

Typical formulations used to describe epidemic models are based on either differential or difference equations. In that way, the basic reproduction number and its physical and biological insight are discussed in [1] which is related to pertussis and measles descriptions. In addition, feedback vaccination laws have been developed using techniques such as sliding-mode control or linear or impulsive feedback vaccination [2,3]. The transient evolution of epidemic diseases is also an important issue for properly describing the day-to-day time-transmission levels and the appropriate eventual interventions to perform since the stability properties are more related to the stationary states, typically the disease-free and the endemic ones. For instance, in [4], an analytic methodology is given to predict and monitor the dates of maximum hospital occupancy of beds. The differential and difference models have also been corrected with other powerful analysis techniques. In that way, the bifurcation analysis and the stability of a fractional order susceptible-infectious-recovered (SIR) epidemic model with delay has been discussed in [5]. On the other hand, discretized and discrete-time epidemic models have been proposed in the background literature. The different approaches can basically consist of the discretization of continuous-time-based models by numerical methods or in the development of discrete models based in difference equations. See, for instance, [6,7,8,9,10,11,12] and some of the references therein. It turns out that since the relevant time of the dynamics evolution in epidemics is relatively long, typically for instance, on the orders of days or weeks, discrete-oriented epidemic models can be found to be appropriate for describing and monitoring the infection evolution through time. An extended SEIR model which incorporates its usual subpopulations, the asymptomatic, and the dead-infective subpopulations has been proposed in [6]. That model has been proposed as appropriate for the Ebola disease. On the other hand, a multistaged SIR was discussed in [7]. Such a model considers several coupled layers of infectious subpopulation for a coupled disease transmission from each layer to the adjacent ones. Moreover, discrete susceptible-infectious-recovered-susceptible (SIRS) models have also been studied in heterogeneous networks [8]. The vaccination effort can be considered either as an external forcing term or as a generator of a new subpopulation, the so-called vaccinated one. Such a subpopulation becomes dynamically coupled to the remaining ones in the model rather than as a specific forcing control [9]. Other types of epidemic models, such as, for instance, discretized susceptible-infectious-recovered (SIR)-type ones, or susceptible-infectious-susceptible (SIS)-type ones, have been proposed in [10,11,12,13] and some of the references therein.

Recently, a lot of research is being dedicated to studying and monitoring the new COVID-19 pandemic using registered data, like infection and detection tests, hospital bed occupancies, and mortality-related records. Normally, such data are updated in discrete time, typically, day to day or week to week. Therefore, discrete epidemic models have been found appropriate for processing such data. See, for instance, [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36] and some related references therein. Studies on particular data for different countries or regions can be found in the literature related to COVID-19, sometimes related to public interventions, such as quarantines, isolation measures, lockdowns, use of masks, social distance rules, etc. See, for instance, those concerned with Saudi Arabia [16,17], Madrid capital town, metropolitan area and surrounding administrative area [18,19], India [23,24], Italy [25], United States [26], Canada and several of its provinces [28], Switzerland, [29], Brazil [30], etc. In addition, the analysis of data has been sometimes accompanied with mathematical analysis techniques on the pandemic evolution related to public interventions or mathematically founded analysis of the obtained data. In that way, the impact of lockdowns is investigated in [17], while the effects of total or partial quarantines are investigated in [18] for a SEIAR model, which incorporates the asymptomatic subpopulation to the typical SEIR model by considering the isolated population as removed either from the infectious individuals or from the susceptible ones. In [19], a more general model of potential usefulness in the description and monitoring of COVID-19 has been proposed, discussed, and tested with recorded real data. Such a model includes three different infectious subpopulations, namely, the slightly infections, the hospitalized ones and the ones staying at the intensive care units are considered. On the other hand, the implementation of control rules oriented either toward reducing the number of exposed individuals or toward increasing the number of treated individuals is proposed and discussed in [20] while impulsive optimal control techniques are developed in [21]. In particular, the analysis technique proposed in [20] relies on the fact that the epidemics is endemic. In [23], the inadequacy of the implementation of open-loop (i.e., without using feedback) controls is emphasized contrarily to the use of closed-loop controls like, for instance, sliding mode-based control or other feedback laws. Moreover, it has been proved in [27] that the suppression strategies are appropriate, provided that they are sufficiently strong while taken through prompt decisions, whereas the mitigation strategies can fail because of eventual unfavorable combinations of delays, unstable dynamics, and uncertainties in the model.

This paper proposes and investigates a new discrete SEIR model subject to a linear two-stage delayed feedback vaccination effort having in mind that some of the recently approved vaccines for COVID-19 require two doses for increasing their average effectiveness. Such doses are administered to the susceptible subpopulation with a delay period, and their potential benefits on immunity appear several days after the administration of the second doses. In fact, the model considers the injection of two doses with different delays and eventual average different effectiveness. The proposed model considers also that the exposed have a transmission rate exposed-susceptible which may be eventually distinct from the infectious-susceptible transmission rate. The underlying idea is that in some infectious diseases, such as, for instance, in the COVID-19 pandemic, there are contagions from both exposed to susceptible and from infectious to susceptible. The delays of the vaccination and expected achievable immunity as well as the transmission rates are model parameters which can be updated for different experiments, [37,38,39].

The paper is organized as follows: Section 2 is devoted to present the new mentioned proposed discrete susceptible-exposed-infectious-recovered (SEIR) epidemic model with delayed double-dose linear feedback vaccination. Section 3 discusses the non-negativity and boundedness of the solution under any given finite non-negative initial conditions, as well as the existence and components of the disease-free equilibrium point and its stability properties. It is also proved that the endemic equilibrium point does not exist for the proposed model in Section 4. Section 5 presents and discusses some examples of the proposed model related to the evolution of the COVID-19 pandemic. Finally, a set of illustrative concluding remarks ends the paper. Some auxiliary technical results are proved in Appendix A which are supported by general necessary mathematical results given or proved in [40,41,42,43,44].

2. The Discrete SEIR Epidemic Model Subject to Two Vaccination Doses

Note that the rationale of the sampling period interpretation is unity, typically one day or one week for a correct practical use of the model. The model parameters should be expressed in values of dimensionality being the inverse of the sampling period units. The SEIR epidemic model equations may be rewritten equivalently as follows:

$$S_{k+1}=\left(a_k-{\beta }_k{I}_k-{\beta }_k^e{E}_k\right)S_k-V_{2,k-d_2}-V_{1,k-d_1-d_2}$$

(1)

$$E_{k+1}=\left(1-\mu \right)E_k+\left({\beta }_k{I}_k+{\beta }_k^e{E}_k\right)S_k$$

(2)

$$I_{k+1}=\left(1-\gamma \right)I_k+\mu E_k$$

(3)

$$R_{k+1}=R_k+\gamma I_k+V_{2,k-d_2}+V_{1,k-d_1-d_2}$$

(4)

$$V_{2,k-d_2}=K_{k-d_2}{\rho }_2{S}_{k-d_2}$$

(5)

$$V_{1,k-d_1-d_2}=K_{k-d_1-d_2}{\rho }_1{S}_{k-d_1-d_2}$$

(6)

for any integer $k\in {\mathit{Z}}_{0+}={\mathit{Z}}_{+}\cup \left\{0\right\}$ and any given finite initial conditions $S_0\ge 0$, $E_0\ge 0$, $I_0\ge 0$ and $R_0\ge 0$ and $S_k=E_k=I_k=R_k=0$ for $k\in {\mathit{Z}}_{-}$, where $S$, $E$, $I$ and $R$ are the susceptible, exposed, infectious and recovered subpopulations, respectively. The forcing terms of Equations (5) and (6) are the two doses of vaccination on the susceptible which are generated via linear feedback and which are subject to integer delays $d_1+d_2$ (first dose) and $d_2$ (second dose), respectively, with respect to each current sampling instant, where $\mathit{min}\left(d_1,d_2\right)>0$. In the above model:

–

$a_k$ is the average recruitment rate proportional to the susceptible at the $kth$ sampling instant related; for instance, to the rates of births and a is a constant reference value for the above sequence; for instance, its average over the whole time period under study and typically it might be unity.

–

${\beta }_k$ and ${\beta }_k^e$ are, respectively, the average transmission rates of the infectious and exposed subpopulations at the $k-th$ sampling instant.

–

$\gamma$ is the average recovery rate.

–

$\mu$ is the average incubation rate.

–

$K_k$ is the vaccination rate (a feedback control gain) which can be eventually depending on the sampling instants. It is assumed in the sequel that ${\left\{K_k\right\}}_0^{\infty }\subset \left[0 , 1\right]$.

–

${\rho }_1$, ${\rho }_1+{\rho }_2$ are parameters in $\left(0 , 1\right]$ which quantify the average effectiveness (or efficiency) of the respective doses. In particular, ${\rho }_2$ gives the extra effectiveness obtained from the injection of the second dose. In this context, ${\rho }_1+{\rho }_2\in \left[0 , 1\right]$ and ${\rho }_1+{\rho }_2=1$ refer to the ideal situation, unattainable in practice, of 100 percent effectiveness of the combined injection of the two doses. Note that ${\rho }_1={\rho }_2={\rho }_1+{\rho }_2=0$ refers to the worst case where the vaccination is fully superfluous.

Typically, we can consider two different transmission rates for the exposed and infectious since the exposed are not usually identified to be allocated under quarantine or isolation, and furthermore, it has been argued that the infective periods are, in general, of distinct time length for both stages in the case of COVID-19 pandemic. It turns out that these transmission rates can be time varying, in general, since the transmission rates can depend on the intervention measures and on the social customs in the geographic area under study. Note that, even if the vaccination is performed under the same gain $K$ for both doses, it can be considered that the effectiveness if only the first injected dose is smaller than if both of them are injected. Therefore, it can be typically argued that ${\rho }_2>0$ leads to ${\rho }_1+{\rho }_2>{\rho }_1.$ It can also be pointed out that the proposed discrete-time model corresponds to the backward Euler discretization on the continuous –time SEIR model given by:

$$\begin{array}{l}\dot{S}\left(t\right)=\nu \left(t\right)S\left(t\right)-\left({\beta }^{\prime }I\left(t\right)+{\beta }^{\prime }{}^{e}E\left(t\right)\right)S\left(t\right)-{V_1}^{\prime }\left(t-d_1-d_2\right)-{V_2}^{\prime }\left(t-d_2\right)\\ \dot{E}\left(t\right)=-{\mu }^{\prime }E\left(t\right)+\left({\beta }^{\prime }I\left(t\right)+{\beta }^{\prime }{}^{e}E\left(t\right)\right)S\left(t\right)\\ \dot{I}\left(t\right)={\mu }^{\prime }E\left(t\right)-{\gamma }^{\prime }I\left(t\right)\\ \dot{R}\left(t\right)={\gamma }^{\prime }I\left(t\right)+{V_1}^{\prime }\left(t-d_1-d_2\right)+{V_2}^{\prime }\left(t-d_2\right)\end{array}$$

whose discretization is given by

$$\begin{array}{l}{S}_{k+1}=S_k+h{\nu }_k{S}_k-\left(h{\beta }^{\prime }{I}_k+h{\beta }^{\prime }{}^e{E}_k\right)S_k-h{V_{1,k-d_1-d_2}}^{\prime }-h{V_{2,k-d_2}}^{\prime }\\ =\left(1+h\nu \right)S_k-\left(h{\beta }^{\prime }{I}_k+h{\beta }^{\prime }{}^e{E}_k\right)S_k-h{V_{1,k-d_1-d_2}}^{\prime }-h{V_{2,k-d_2}}^{\prime }\\ E_{k+1}=E_k-h{\mu }^{\prime }{E}_k+\left(h{\beta }^{\prime }{I}_k+h{\beta }^{\prime }{}^e{E}_k\right)S_k=\left(1-h{\mu }^{\prime }\right)E_k+\left(h{\beta }^{\prime }{I}_k+h{\beta }^{\prime }{}^e{E}_k\right)S_k\\ I_{k+1}=I_k+h{\mu }^{\prime }{E}_k-h{\gamma }^{\prime }{I}_k=\left(1-h{\gamma }^{\prime }\right)I_k+h{\mu }^{\prime }{E}_k\\ R_{k+1}=R_k+h{\gamma }^{\prime }I\left(t\right)+h{V_{1,k-d_1-d_2}}^{\prime }+h{V_{2,k-d_2}}^{\prime }\end{array}$$

so that the discretized and the proposed discrete-time model are equivalent if the following correspondence is performed

$$\begin{array}{l}{a}_k=1+h{\nu }_k\\ \beta =h{\beta }^{\prime },{\beta }^e=h{\beta }^{\prime }{}^e\\ V_{1,k-d_1-d_2}=h{V_{1,k-d_1-d_2}}^{\prime },V_{2,k-d_2}=h{V_{2,k-d_2}}^{\prime }\\ \mu =h{\mu }^{\prime },\gamma =h{\gamma }^{\prime }\end{array}$$

However, the model is originally set up in discrete-time instead as the discretization of a continuous-time one. One of the daunting challenges to counteract the epidemic spreading is the design of the vaccination functions. The control theory provides a framework to systematically design them in order to fulfill with a prescribed control objective. Thus, an output feedback approach is adopted in this work and discussed in the sequel. However, control theory has developed many analytical tools that allows one to face this problem by using other approaches such as state feedback and feedback linearization [2], output-controllability [19], optimal control [21], impulsive control [22], sliding mode [23] and lockdown and quarantine [18,28], to cite just a few.

The model (1)–(6) may be compacted after replacing the vaccination controls (5)–(6) in the susceptible and recovered subpopulations dynamics, which result in:

$$S_{k+1}=a{S}_k-K\left({\rho }_2{S}_{k-d_2}+{\rho }_1{S}_{k-d_1-d_2}\right)+u_k$$

(7)

$$E_{k+1}=\left(1-\mu \right)E_k+\left({\beta }_k{I}_k+{\beta }_k^e{E}_k\right)S_k$$

(8)

$$I_{k+1}=\left(1-\gamma \right)I_k+\mu E_k$$

(9)

$$R_{k+1}=R_k+\gamma I_k+K\left({\rho }_2{S}_{k-d_2}+{\rho }_1{S}_{k-d_1-d_2}\right)-g_k$$

(10)

$$u_k=\left[a_k-a-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)\right]S_k+g_k$$

(11)

$$g_k=\left(K-K_{k-d_2}\right){\rho }_2{S}_{k-d_2}+\left(K-K_{k-d_1-d_2}{\rho }_1\right)S_{k-d_1-d_2}$$

(12)

for any integer $k\in {\mathit{Z}}_{0+}$, where ${\lambda }_k^e={\beta }_k^e/{\beta }_k$ and $K$ and $a$ are reference values for the vaccination rate gain and recruitment, for instance, the estimated average values of their asymptotic limits if they converge. In particular, a scalar discrete function of the form of the susceptible subpopulation (7) is discussed in detail in Appendix A from the points of view of stability and convergence to a limit. Those auxiliary results based on a similar discrete Equation (7) subject to (11) are then used to discuss the stability of the proposed epidemic model.

3. Non-Negativity, Stability and Disease-Free Equilibrium Point

The following result which relies on boundedness, convergence and non-negativity of the solution is direct by summing up (7) to (10) (or (1) to (4)) by considering (7) to (11) (or (1) to (6)):

Theorem 1.

The following properties hold:

(i)

Assume that$N_0=S_0+E_0+I_0+R_0=1$. Then, the total population$N_k=S_k+E_k+I_k+R_k=1$; $\forall k\in {\mathit{Z}}_{0+}$if$a_k\equiv 1$.

$$(\mathrm{ii}) N_k=1+{{\displaystyle \sum }}_{j=0}^{k-1}\left(a_j-1\right)S_j; \forall k\in {\mathit{Z}}_{0+}$$

(iii)

${\left\{S_k\right\}}_{k=0}^{\infty }$is non-increasing, and then bounded and convergent, if

$$V_{2,k-d_2}+V_{1,k-d_1-d_2}=K_{k-d_2}{\rho }_2{S}_{k-d_2}+K_{k-d_1-d_2}{\rho }_1{S}_{k-d_1-d_2}\ge \left(a_k-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)-1\right)S_k; \forall k\in {\mathit{Z}}_{0+}$$

(13)

leading to

$$1+S_{k-1}{{\displaystyle \sum }}_{j=0}^{k-1}\left(a_j-1\right)\le N_k\le 1+S_0{{\displaystyle \sum }}_{j=0}^{k-1}\left(a_j-1\right);\forall k\in {\mathit{Z}}_{0+}$$

(14)

and, if the average value of the sequence${\left\{a_k\right\}}_{k=0}^{\infty }$is unity, then there exists the limit$\underset{k\to \infty }{lim}{N}_k=N_0=1$. ${\left\{S_k\right\}}_{k=0}^{\infty }$is strictly decreasing if the inequality in (13) is strict.

(iv)

${\left\{S_k\right\}}_{k=0}^{\infty }$is non-negative if and only if

$$V_{2,k-d_2}+V_{1,k-d_1-d_2}=K_{k-d_2}{\rho }_2{S}_{k-d_2}+K_{k-d_1-d_2}{\rho }_1{S}_{k-d_1-d_2}\le \left(a_k-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)\right)S_k ; \forall k\in {\mathit{Z}}_{0+}$$

(15)

provided that the infectious transmission rate is sufficiently small according to${\beta }_k\le \frac{a_k}{I_k+{\lambda }_k^e{E}_k{}^{}}$; $\forall k\in {\mathit{Z}}_{0+}$which is guaranteed if${\beta }_k\le \frac{a_k}{{\left(1+{\lambda }_k^e\right)}^{}{N}_k}$; $\forall k\in {\mathit{Z}}_{0+}$.

(v)

Assume that$\mu \in \left[0 , 1\right)$, $\gamma \in \left[0 , 1\right)$and that (15) holds. Then, any sequence trajectory solution of (1)–(4) subject to (5)–(6) is non-negative.

Proof. 

If $a_k=1$; $\forall k\in {\mathit{Z}}_{0+}$ then $N_0=S_0+E_0+I_0+R_0=1$ implies directly that $N_k=S_k+E_k+I_k+R_k=N_0=1$; $\forall k\in {\mathit{Z}}_{0+}$ by summing up (1) to (4). The same calculations prove Property (ii) in the general case. To prove, Property (iii), note from (3) that, under the given constraint, ${\left\{S_k\right\}}_{k=0}^{\infty }$ is non-increasing so that (14) holds directly from Property (ii) since:

$$\frac{S_{k+1}}{S_k}=a_k-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)-\frac{V_{2,k-d_2}+V_{1,k-d_1-d_2}}{S_k}\le 1; \forall k\in {\mathit{Z}}_{0+}$$

(16)

Note that ${\left\{S_k\right\}}_{k=0}^{\infty }$ is also bounded and convergent since it is non-increasing and its initial condition is finite by hypothesis. Property (iv) follows directly by replacing (15) with (1) subject to $S_0\ge 0$ and to ${\beta }_k\le \frac{a_k}{I_k+{\lambda }_k^e{E}_k{}^{}}$ (guaranteed if ${\beta }_k\le \frac{a_k}{{\left(1+{\lambda }_k^e\right)}^{}{N}_k}$); $\forall k\in {\mathit{Z}}_{0+}$ since $V_k\ge 0$; $\forall k\in {\mathit{Z}}_{0+}$. If, addition, $\mu \in \left[0 , 1\right)$, $\gamma \in \left[0 , 1\right)$ then from (2)–(4) and ${\left\{S_k\right\}}_{k=0}^{\infty }\subset {\mathit{R}}_{0+}={\mathit{R}}_{+}\cup \left\{0\right\}$ yields ${\left\{E_k\right\}}_{k=0}^{\infty }\subset {\mathit{R}}_{0+}$, ${\left\{I_k\right\}}_{k=0}^{\infty }\subset {\mathit{R}}_{0+}$ and ${\left\{R_k\right\}}_{k=0}^{\infty }\subset {\mathit{R}}_{0+}$, since the initial conditions of all the subpopulations are non-negative. Property (v) is proved. □

Remark 1.

The non-negativity property of the susceptible of Theorem 1 (iv) is a necessary constraint for a well-posedness of the model. Note that the condition of sufficiently smaller transmission rate of the Infectious${\beta }_k\le \frac{a_k}{{\left(1+{\lambda }_k^e\right)}^{}{N}_k}$;$\forall k\in {\mathit{Z}}_{0+}$, which guarantees according to Theorem 1 (iv) that the sequence of susceptible is non-negative if$S_0\ge 0$and the vaccination satisfies (15), is reasonable since the transmission rate decreases as the total population increases in the popular class of true-mass action epidemic models. The property of the susceptible sequence being non-increasing [Theorem 1 (iii)] describes appropriately the disease growing since it starts until its first maximum peak.

Remark 2.

It turns out the need of the assumption${\left\{K_k\right\}}_0^{\infty }\subset \left[0 , 1\right]$for the vaccination gain sequence since the daily vaccination (or, in general, the one for the used sampling period in the model parameterization) is proportional to the susceptible subpopulation. Therefore, note that the condition (13) for${\left\{S_k\right\}}_{k=0}^{\infty }$to be non-increasing requires the necessary condition:

$${\rho }_2{S}_{k-d_2}+{\rho }_1{S}_{k-d_1-d_2}\ge \left(a_k-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)-1\right)S_k; \forall k\in {\mathit{Z}}_{0+}$$

(17)

to be fulfilled by the vaccination gain unity, which is the maximum possible one, provided that the existing stock vaccines covers such a need. Note that the above condition always holds if${\left\{S_k\right\}}_{k=0}^{\infty }$is non-negative and$a_k\le {\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)+1$;$\forall k\in {\mathit{Z}}_{0+}$and it also holds if${\left\{S_k\right\}}_{k=0}^{\infty }$is non-negative,$a_k>{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)+1$;$\forall k\in {\mathit{Z}}_{0+}$. This second condition is guaranteed by under unity vaccination gain:

$$\frac{\left({\rho }_1+{\rho }_2\right)\mathit{min}\left(S_{k-d_2} , S_{k-d_1-d_2}\right)}{\left(a_k-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)-1\right)S_k}\ge 1;\forall k\in {\mathit{Z}}_{0+}$$

(18)

In practice,$\frac{\mathit{min}\left(S_{k-d_2} , S_{k-d_1-d_2}\right)}{S_k}$is slightly greater than unity if the susceptible sequence is decreasing since$d_1$and$d_1+d_2$are typically delays of one and three weeks, in the case of SARS-CoV-2 which do not affected significantly the variation of the susceptible population levels. In addition,$a_k-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)-1$is typically close to zero for small infection levels compared to the susceptible population values. Therefore, for reasonable vaccine efficacies, (18) holds if$a_k>{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)+1$;$\forall k\in {\mathit{Z}}_{0+}$.

Note that the condition (15) for the susceptible sequence${\left\{S_k\right\}}_{k=0}^{\infty }$to be non-negative requires the necessary condition to be fulfilled by some maximum possible vaccination gain rate$K\le 1$(since this feature would ensure that it holds for any lower effort as well), that is,

$\left(a_k-{\beta }_k\left(I_k+{\lambda }_k^e{E}_k\right)\right)S_k\ge 0$;$\forall k\in {\mathit{Z}}_{0+}$which holds if$S_0\ge 0$if${\beta }_k\le \frac{a_k}{I_k+{\lambda }_k^e{E}_k{}^{}}$;$\forall k\in {\mathit{Z}}_{0+}$which is guaranteed if${\beta }_k\le \frac{a_k}{{\left(1+{\lambda }_k^e\right)}^{}{N}_k}$;$\forall k\in {\mathit{Z}}_{0+}$.

The existence of a disease- free equilibrium point is discussed in the subsequent result:

Theorem 2.

Assume that${\rho }_1+{\rho }_2\le 1$,${\beta }_k\le \frac{a_k}{I_k+{\lambda }_k^e{E}_k{}^{}}$(guaranteed if${\beta }_k\le \frac{a_k}{{\left(1+{\lambda }_k^e\right)}^{}{N}_k}$);$\forall k\in {\mathit{Z}}_{0+}$,${\left\{a_k\right\}}_{k=0}^{\infty }\to a\left(\le 1\right)$,${\left\{K_k\right\}}_{k=0}^{\infty }\left(\subset \left[0 , 1\right]\right)\to K$and assume also that the constraints (13) and (15) hold. Then, there is a disease-free equilibrium point$x_{df}={\left(S_{df}, 0,0,R_{df}=N_{df}-S_{df}\right)}^T$, being in general dependent on the initial conditions, with$S_{df}\ge 0$and$R_{df}=N_{df}-S_{df}$, where${\left\{N_k\right\}}_{k=0}^{\infty }\to N_{df}$, such that$S_{df}=0$if$K\ne \frac{a}{{\rho }_1+{\rho }_2}$.

Proof. 

Note that ${\left\{K_k\right\}}_{k=0}^{\infty }\left(\subset \left[0, 1\right]\right)\to K$ implies that $K\in \left[0, 1\right]$. Note also that $E_{df}=0$ and $I_{df}=0$ are trivial solutions of (2) and (3) if ${\left\{S_k\right\}}_{k=0}^{\infty }$ is bounded. It follows furthermore that ${\left\{S_k\right\}}_{k=0}^{\infty }$ is non-increasing and bounded from Theorem 1 (iii), ${\left\{a_k\right\}}_{k=0}^{\infty }\to a\left(\le 1\right)$ and (13) holds. As a result, if furthermore, (15) holds then ${\left\{S_k\right\}}_{k=0}^{\infty }\left(\subset {\mathit{R}}_{0+}\right)\to S_{df}\left(\ge 0\right)$ with $S_{df}=0$ if $K\ne \frac{a}{{\rho }_1+{\rho }_2}$ subject to the assumed constraint ${\rho }_1+{\rho }_2\ge 1$ since $K\in \left[0, 1\right]$. Note that

$$N_{k+1}=a_k{S}_k+E_k+I_k+R_k=N_k+\left(a_k-1\right)S_k$$

(19)

implies that $N_{k+1}-N_k=\left(a_k-1\right)S_k$ and ${\left\{a_k\right\}}_{k=0}^{\infty }\to 1$ implies that ${\left\{N_k\right\}}_{k=0}^{\infty }\to N_{df}$ and ${\left\{R_k\right\}}_{k=0}^{\infty }\to R_{df}=N_{df}-S_{df}$ since ${\left\{E_k\right\}}_{k=0}^{\infty }\to E_{df}=0$ and ${\left\{I_k\right\}}_{k=0}^{\infty }\to I_{df}=0$. □

The subsequent result gives conditions for both the infection and the susceptibility to asymptotically vanish implying convergence of the solution to a disease-free equilibrium point without susceptible individuals. Some supportive auxiliary results used in its proof are given in Appendix A.

Theorem 3.

Assume that the constraint (15) holds, so that any solution sequence is non-negative for non-negative initial conditions. Assume also that the first assumptions below, and, furthermore, at least one of the assumptions 2 to 5 below also holds:

(1)

$\mu \in \left[0 , 1\right)$,$\gamma \in \left[0 , 1\right)$,${\rho }_1+{\rho }_2\le 1$,$K_k\le \frac{a_k}{I_k+{\lambda }_k^e{E}_k}$(guaranteed if${\beta }_k\le \frac{a_k}{{\left(1+{\lambda }_k^e\right)}^{}{N}_k}$);$\forall k\in {\mathit{Z}}_{0+}$,${\left\{a_k\right\}}_{k=0}^{\infty }\to a\left(\le 1\right)$,${\left\{{\beta }_k\right\}}_{k=0}^{\infty }\to \beta$,${\left\{{\lambda }_k^e\right\}}_{k=0}^{\infty }\to {\lambda }^e$,${\left\{K_k\right\}}_{k=0}^{\infty }\to K$.

(2)

Either the constraint (13) holds with strict inequality for all$k\in {\mathit{Z}}_{0+}$and$K\ne \frac{a}{{\rho }_1+{\rho }_2}$.

(3)

$K\in (0 , \mathit{min}( \stackrel{—}{K} , 1) )$, where:

$$\stackrel{—}{K}=\mathit{sup}\left\{y\in \left(\frac{a-1}{{\rho }_1+{\rho }_2} , \frac{1+a}{{\rho }_1+{\rho }_2}\right) :y<1/f\left(y\right)\right\}$$

(20)

$$f\left(y\right)=\underset{\theta \in \left(0,2\pi \right)}{sup}\sqrt{\frac{{\left({\rho }_2\left(\mathit{cos}\left(\theta d_2\right)-1\right)+{\rho }_1\left(\mathit{cos}\left(\theta \left(d_1+d_2\right)\right)-1\right)\right)}^2+{\left({\rho }_2\mathit{sin}\left(\theta d_2\right)+{\rho }_1\mathit{sin}\left(\theta \left(d_1+d_2\right)\right)\right)}^2}{1+{\left(a-\left({\rho }_1+{\rho }_2\right)y\right)}^2-2\left(a-\left({\rho }_1+{\rho }_2\right)y\right)\mathit{cos}\theta }}$$

(21)

and

$$\left|K_k-K\right|=o\left(\Vert {\widehat{x}}_k\Vert {}^{-1}\right); \left|a_k-a\right|=o\left(\Vert {\widehat{x}}_k\Vert {}^{-1}\right)$$

(22)

where${\widehat{x}}_k={\left(x_k, x_{k-1},\dots ,x_{k-d_2},\dots ,x_{k-d_1-d_2}\right)}^T$.

$$(4) a\in \left(0 , 1\right), K\in \left[0 , 1/\underset{\theta \in \left[0,2\pi \right)}{sup}\left|\frac{{\rho }_2{e}^{-i\theta d_2}+{\rho }_1{e}^{-i\theta \left(d_1+d_2\right)}}{e^{i\theta }-a}\right|\right)$$

(23)

and (22) holds.

$$(5) a\in \left(0 , 1\right),K\in \left[0 , \frac{\left(c_1^{-d_1}-a\right)c_1}{{\rho }_1+{\rho }_2{c}_2}\right)$$

(24)

Then, there exists a globally asymptotically stable disease-free equilibrium point$x_{df}={\left(0, 0,0,R_{df}=N_{df}\right)}^T$and${\left\{N_k\right\}}_{k=0}^{\infty }\to N_{df}$.

Proof. 

The above assumptions 1 and 2 guarantee that ${\left\{S_k\right\}}_{k=0}^{\infty }\to 0$, and that ${\left\{S_k\right\}}_{k=0}^{\infty }$ is bounded, from Theorem 1 and Theorem 2 for any given non-negative finite initial conditions. Note also that, from (2) and (3), ${\left\{E_k\right\}}_{k=0}^{\infty }\to 0$ and ${\left\{I_k\right\}}_{k=0}^{\infty }\to 0$ since $\mu \in \left[0, 1\right)$ and $\gamma \in \left[0, 1\right)$.

In addition, any of the assumptions 3 or 4, together with the properties ${\left\{E_k\right\}}_{k=0}^{\infty }\to 0$ and ${\left\{I_k\right\}}_{k=0}^{\infty }\to 0$ guarantees that ${\left\{S_k\right\}}_{k=0}^{\infty }\to 0$, and ${\left\{S_k\right\}}_{k=0}^{\infty }$ is bounded, for any given non-negative finite initial conditions from Theorem A1 [(i)–(ii)] and Corollary A1 of Appendix A. The same result also holds under the assumption 5 according to Theorem A3 (i) of Appendix A.

Thus, the disease-free equilibrium point $x_{df}={\left(0, 0,0,R_{df}=N_{df}\right)}^T$ exists under the joint assumptions 1 and 2 and under the assumption 1 jointly with one of the assumptions 3 to 5. Now, note that

Since ${\left\{{\beta }_k{S}_k\right\}}_{k=0}^{\infty }\to \beta S_{df}=0$, Equations (2) and (3) can be compacted as follows:

$$\begin{array}{ll}\left[\begin{array}{c}{E}_{k+1}\\ I_{k+1}\end{array}\right]& =\left[\begin{array}{cc}1-\mu +{\lambda }_k^e{\beta }_k{S}_k& {\beta }_k{S}_k\\ \mu & 1-\gamma \end{array}\right]\left[\begin{array}{c}{E}_k\\ I_k\end{array}\right]\\ & =\left[\begin{array}{cc}1-\mu +{\lambda }^e\beta S_{df}& \beta S_{df}\\ \mu & 1-\gamma \end{array}\right]\left[\begin{array}{c}{E}_k\\ I_k\end{array}\right]+\left[\begin{array}{cc}\beta \left({\lambda }_k^e{S}_k-{\lambda }^e{S}_{df}\right)& \beta \left(S_k-S_{df}\right)\\ 0& 0\end{array}\right]\left[\begin{array}{c}{E}_k\\ I_k\end{array}\right]\\ & =\left[\begin{array}{cc}1-\mu & 0\\ \mu & 1-\gamma \end{array}\right]\left[\begin{array}{c}{E}_k\\ I_k\end{array}\right]+\left[\begin{array}{cc}\beta {\lambda }_k^e{S}_k& \beta S_k\\ 0& 0\end{array}\right]\left[\begin{array}{c}{E}_k\\ I_k\end{array}\right]; k\in {\mathit{Z}}_{0+}\end{array}$$

(25)

Note that the spectral radius $r$ of the matrix $\left[\begin{array}{cc}1-\mu & 0\\ \mu & 1-\gamma \end{array}\right]$ is $\mathit{max}\left(1-\mu , 1-\gamma \right)<1$ since $\mu \in \left[0, 1\right)$ and $\gamma \in \left[0, 1\right)$ and the above matrix is a convergent matrix. Since there is some matrix norm of arbitrary close value to the spectral radius there exist some positive real constants $\epsilon \le {\delta }_0-1+\mathit{min}\left(\mu ,\gamma \right)$ and $r\le {\delta }_0<1$ such that $\Vert \left[\begin{array}{cc}1-\mu & 0\\ \mu & 1-\gamma \end{array}\right]\Vert \le {\delta }_0$. In addition, since ${\left\{{\beta }_k{S}_k\right\}}_{k=0}^{\infty }\to \beta S_{df}=0$ and ${\left\{{\lambda }_k^e\right\}}_{k=0}^{\infty }\to {\lambda }^e$ and a decreasing nonnegative real sequence ${\left\{{\delta }_{1k}\right\}}_{k=0}^{\infty }\to 0$ having a strictly decreasing subsequence ${\left\{{\delta }_{1N_{n_k}}\right\}}_{k=0}^{\infty }\to 0$ such that for any given ${\delta }_{1N_{n_k}}$, there exists a non-negative integer $N_{n_k}$ such that $\Vert \left[\begin{array}{cc}\beta {\lambda }_k^e{S}_k& \beta S_k\\ 0& 0\end{array}\right]\Vert \le {\delta }_{1N_k}$; $j=N_{n_k},N_{n_k+1},\dots ,N_{n_{k+1}-1}$ such that ${\left\{N_{n_k}\right\}}_{k=0}^{\infty }$ is strictly increasing. Thus, one gets from (25) that:

$$\Vert \left[\begin{array}{c}{E}_{N_{n_{k+1}}}\\ I_{N_{n_{k+1}}}\end{array}\right]\Vert \le {\delta }_0^{N_{n_{k+1}}-N_{n_k}}\Vert \left[\begin{array}{c}{E}_{N_{n_k}}\\ I_{N_{n_k}}\end{array}\right]\Vert +{\delta }_{1N_k}{{\displaystyle \sum }}_{j=n_k}^{n_{k+1}-1}{\delta }_0^{n_{k+1}-1-j}; k\in {\mathit{Z}}_{0+}$$

(26)

so that ${\left\{\Vert \left[\begin{array}{c}{E}_{N_{n_k}}\\ I_{N_{n_k}}\end{array}\right]\Vert \right\}}_{k=0}^{\infty }$ has a strictly decreasing subsequence which converges to zero for any given finite non-negative initial conditions. In view of (2)–(3), that property also implies that ${\left\{E_k\right\}}_{k=0}^{\infty }\to 0$ and ${\left\{I_k\right\}}_{k=0}^{\infty }\to 0$ for any given non-negative finite initial conditions and then (25) is globally asymptotically stable. Since it has been already proved that ${\left\{S_k\right\}}_{k=0}^{\infty }\to 0$ then $x_{df}={\left(0, 0,0,R_{df}=N_{df}\right)}^T$ is globally asymptotically stable and ${\left\{N_k\right\}}_{k=0}^{\infty }\to N_{df}$. □

The next result gets condition for which the stationary model has a stable infective substate, irrespective of the susceptible being zero at the equilibrium point or not. It is found that there exists a critical transmission rate under which the stationary substate of the infective model is globally asymptotically stable. Under the extra conditions from Theorem 3, the critical transmission rate becomes infinity since the susceptible subpopulation is zero at the disease-free equilibrium point. As a result, the disease-free equilibrium point is globally asymptotically stable for any value of the transmission rate, and there is no endemic attractor.

Theorem 4.

Assume$\mu \in \left[0 , 1\right)$and$\gamma \in \left[0 , 1\right)$. Then, the stationary infective substate (25) has stable characteristic roots if$\beta \in \left[0, {\beta }_c \right)$, where${\beta }_c=\frac{\mu \gamma }{{\left(\mu +{\lambda }^e\gamma \right)}^{}{S}_{df}}$so that it is globally asymptotically stable.

In addition, if the conditions of Theorem 2 hold with$K\ne \frac{a}{{\rho }_1+{\rho }_2}$then${\beta }_c=+\infty$.

The stable solution sequence is guaranteed to be, furthermore, non-negative for any non-negative finite initial conditions if all the conditions of Theorem 3 hold.

Proof. 

If ${\left\{{\beta }_k{S}_k\right\}}_{k=0}^{\infty }\to \beta S_{df}=s$ and ${\left\{{\lambda }_k^e\right\}}_{k=0}^{\infty }\to {\lambda }^e$, then the characteristic Equation of (25) is:

$${\lambda }^2+b\lambda +c=0$$

(27)

where

$$b=b\left(s\right)=\mu +\gamma -\left(2+{\lambda }^{e}s\right)$$

(28)

$$c=c\left(s\right)=\left(1-\mu +{\lambda }^{e}s\right)\left(1-\gamma \right)-\mu s$$

(29)

Note that $b\left(0\right)=\mu +\gamma -2<0$ and $c\left(0\right)=\left(1-\mu \right)\left(1-\gamma \right)>0$. The coefficients of the polynomial defining the characteristic equation for $s=0$, ordered in decreasing order, have two changes of sign. According to Descartes’ rule of signs, such a polynomial has either two or zero positive real roots. Note that by simple inspection or, equivalently, according to the unforced part of (25), the characteristic Equation (27) has in fact two positive real characteristic roots for $s=0$ within the unit complex circle centered at zero which are ${\lambda }_1=1-\mu$ and ${\lambda }_2=1-\gamma$. Now, by inspection, $b\left(s\right)=b\left(0\right)-{\lambda }^{e}s<0$ for all $s\ge 0$. It is now discussed a valid range of non-negative values of $s$ which guarantees that $c\left(s\right)=c\left(0\right)+\left({\lambda }^e\left(1-\gamma \right)-\mu \right)s$ remains to be positive. Write by convenience $s=\frac{\sigma \mu \gamma }{\mu +{\lambda }^e\gamma }$. Then, $c\left(s\right)>0$ if and only if $\widehat{c}\left(\sigma \right)>0$, where: $\widehat{c}\left(\sigma \right)=\left(\mu +{\lambda }^e\gamma \right)c\left(s\right)=\left(\mu +{\lambda }^e\gamma \right)c\left(0\right)-\sigma \mu \gamma \left(\mu -{\lambda }^e\left(1-\gamma \right)\right)>0$ for any real $\sigma \sigma$ if ${\lambda }^e\ge \mu /\left(1-\gamma \right)$ with the eventual negative real values having no interest since $s\ge 0$ and, if ${\lambda }^e<\mu /\left(1-\gamma \right)$ then $\widehat{c}\left(\sigma \right)>0$ if

$$\sigma <\frac{\mu +{\lambda }^e\gamma }{\mu -{\lambda }^e\left(1-\gamma \right)}\frac{\left(1-\mu \right)\left(1-\gamma \right)}{\mu \gamma }=\frac{\mu +{\lambda }^e\gamma }{\mu +{\lambda }^e\gamma -{\lambda }^e}\left(\frac{1-\mu -\gamma }{\mu \gamma }+1\right)=1+{\sigma }_0$$

for some real ${\sigma }_0>0$. In both cases, that is, irrespective of ${\lambda }^e$,$c\left(s\right)>0$ if $s=\frac{\sigma \mu \gamma }{\mu +{\lambda }^e\gamma }$ and $\sigma \in \left[0 , 1\right]$. Since there are two real positive characteristic roots for $s=0$ from Descartes rule of signs (since $b\left(0\right)<0$ and $c\left(0\right)>0$) and since the characteristic roots are continuous functions of the argument $s$ and $b\left(s\right)<0$ and $c\left(s\right)>0$ for $s=\frac{\sigma \mu \gamma }{\mu +{\lambda }^e\gamma }$ with $\sigma \in \left[0 , 1\right]$, one concludes that the characteristic roots are real (and positive) if $s=\beta S_{df}<\frac{\mu \gamma }{\mu +{\lambda }^e\gamma }$. It remains now to prove that those real positive roots are within the open unit circle. The two (real) within the complex open circle or radius unity centered at the origin of the complex plane if and only if both roots of (27), subject to (28)–(29), are within the open unit circle centered at zero in the complex plane, namely if and only if:

$$-1<{\lambda }_{1,2}=\frac{2-\mu -\gamma +{\lambda }^{e}s\pm \sqrt{{\left(2-\mu -\gamma +{\lambda }^{e}s\right)}^2+4\left(\mu s-\left(1-\mu +{\lambda }^{e}s\right)\left(1-\gamma \right)\right)}}{2} <1$$

(30)

equivalently, if and only if,

$$\begin{array}{cc}\mu +\gamma -{\lambda }^{e}s-4& <-\sqrt{{\lambda }^{e^2}{s}^2+{\left(\gamma -\mu \right)}^2+4\mu s+2{\lambda }^{e}s\left(\gamma -\mu \right)}\\ & \le \sqrt{{\lambda }^{e^2}{s}^2+{\left(\gamma -\mu \right)}^2+4\mu s+2{\lambda }^{e}s\left(\gamma -\mu \right)}<\mu +\gamma -{\lambda }^{e}s\end{array}$$

(31)

equivalently, if and only if,

$$\sqrt{{\lambda }^{e^2}{s}^2+{\left(\gamma -\mu \right)}^2+4\mu s+2{\lambda }^{e}s\left(\gamma -\mu \right)}<\mathit{min}\left(\mu +\gamma -{\lambda }^{e}s , 4+{\lambda }^{e}s-\mu -\gamma \right)=\mu +\gamma -{\lambda }^{e}s$$

(32)

since the last equality follows directly since $\mu ,\gamma \in \left(0 , 1\right)$. Taking the squares in both sides of (32) yields that it is equivalent to

$${\lambda }^{e^2}{s}^2+{\left(\gamma -\mu \right)}^2+4\mu s+2{\lambda }^{e}s\left(\gamma -\mu \right)<{\lambda }^{e^2}{s}^2+{\left(\mu +\gamma \right)}^2-2{\lambda }^{e}s\left(\mu +\gamma \right)$$

(33)

and also equivalent to

$${\left(\gamma -\mu \right)}^2+4\mu s+2{\lambda }^{e}s\left(\gamma -\mu \right)<{\left(\mu +\gamma \right)}^2-2{\lambda }^{e}s\left(\mu +\gamma \right)$$

(34)

and, furthermore, equivalent to

$${\left(\gamma -\mu \right)}^2+4\mu s+4{\lambda }^{e}s\gamma <{\left(\mu +\gamma \right)}^2$$

(35)

which holds if and only if

$$s=\beta S_{df}<\frac{\mu \gamma }{\mu +{\lambda }^e\gamma }$$

(36)

Since $s=\beta S_{df}\ge 0$, the constraint (36) implies that the characteristic Equation (27) of the stationary infective substate (25) has stable roots if $\beta \in \left[0, {\beta }_c \right)$, where ${\beta }_c=\frac{\mu \gamma }{{\left(\mu +{\lambda }^e\gamma \right)}^{}{S}_{df}}$. Now, assume that all the conditions of Theorem 2 hold. Since $S_{df}=0$ if $K\ne \frac{a}{{\rho }_1+{\rho }_2}$ then ${\beta }_c=+\infty$ if $K\ne \frac{a}{{\rho }_1+{\rho }_2}$ and the stationary infective substate (25) has stable roots if $\beta \in \left[0, +\infty \right)$.

In addition, any solution is non-negative for all time if (15) and assumption 1 plus any of the four remaining assumptions 2 to 5 of Theorem 3 holds (according to Theorem 3). Thus, the joint non-negativity and global asymptotic stability of the disease-free stationary solution, under any finite non-negative initial conditions, are guaranteed for $\beta \in \left[0, +\infty \right)$ since $S_{df}=0$ (Theorem 3). □

Note that the Descartes rule of signs used in the proof of Theorem 4 is supported by the fact that for $s=0$ both characteristic roots are real, positive and stable, so that the possibility of having no positive root is excluded by such a value of the argument $s$, and the fact that the characteristic roots are continuous functions of $s$, so that they continue to be within the positive region of the unity circle centered at zero as the argument increases.

Note also that the proposed model is claimed for its usefulness for short-term predictions in the evolution phase when the disease is blowing up. It is neither considered that vaccination is available for use nor that there is immunity lost allowing to increase again the susceptible numbers after a certain delay. Therefore, the evolution of the susceptible subpopulation is given by a decreasing sequence. It is now proved that the proposed model does not have an endemic equilibrium point.

4. Nonexistence of Endemic Equilibrium Point

It is now discussed under which reasonable conditions an endemic equilibrium point $x_{end}={\left(S_{end}, E_{end}\left(\ne 0\right), I_{end}\left(\ne 0\right) , R_{end}\right)}^T$ either exists or it does not exist. As a final result, it is concluded that the endemic equilibrium point cannot exist.

Consider the following possible cases:

Case 1$(S_{end}>0)$. Combining (2) and (3) for stationary limiting model parameters for an assumed to exist endemic equilibrium point leads to provided that $E_{end}\ne 0$ and $S_{end}\ne 0$:

$$I_{end}=\frac{\mu }{\gamma }{E}_{end}$$

(37)

$$S_{end}=\frac{\gamma \mu }{\beta \left(\mu +{\lambda }^e\gamma \right)}>0$$

(38)

which, replaced in (1) and after cancelling common factors with $S_{end}>0$, results in the resulting constraint:

$$S_{end}={\left[a-\beta \left(I_{end}+{\lambda }^e{E}_{end}\right)-K\left({\rho }_1+{\rho }_2\right)\right]}^{}{S}_{end}$$

(39)

leads to

$$E_{end}=\frac{a-1-K\left({\rho }_1+{\rho }_2\right)}{\beta \left(\mu +\gamma {\lambda }^e\right)}\gamma$$

(40)

if

$$a\ne \beta \left(I_{end}+{\lambda }^e{E}_{end}\right)+K\left({\rho }_1+{\rho }_2\right)$$

(41)

so that $E_{end}>0$ (and also $I_{end}>0$ from (37)) if and only if $K<\frac{a-1}{{\rho }_1+{\rho }_2}\le 1$ and $1<a\le 1+{\rho }_1+{\rho }_2\le 2$. If $K=\frac{a-1}{{\rho }_1+{\rho }_2}$ and $1\le a\le 1+{\rho }_1+{\rho }_2$ then $E_{end}=I_{end}=0$ so that the equilibrium point is, in fact, a disease-free one. By considering also the stationary recovered, one gets after cancelling the stationary recovered $R_{k+1}=R_k=R_{end}$ in both sides of (4) that the subsequent constraint should hold:

$$\gamma I_{end}+K\left({\rho }_1+{\rho }_2\right)S_{end}=\mu E_{end}+K\left({\rho }_1+{\rho }_2\right)S_{end}=0$$

(42)

which contradicts that $S_{end}>0$ and $E_{end}>0$. As a result, no endemic equilibrium point exists with $S_{end}>0$. Since (42) can only hold with $E_{end}=S_{end}=0$ then, it only could eventually be true for the disease-free equilibrium point. Thus, Case 1 is unfeasible for the existence of an endemic equilibrium point.

Case 2$\left(S_{end}\ge 0\right)$. Note that (40) holds with, which includes also $S_{end}>0$ of Case 1 if

$$a=\beta \left(I_{end}+{\lambda }^e{E}_{end}\right)+K\left({\rho }_1+{\rho }_2\right)=\beta \left(\frac{\mu }{\gamma }+{\lambda }^e\right)E_{end}+K\left({\rho }_1+{\rho }_2\right)$$

(43)

provided that

$$E_{end}=\frac{a-K\left({\rho }_1+{\rho }_2\right)}{\beta \left(\mu +\gamma {\lambda }^e\right)}\gamma$$

(44)

provided that

$$K<\frac{a}{{\rho }_1+{\rho }_2}\le 1; a\le {\rho }_1+{\rho }_2$$

(45)

with $S_{end}=0$, and

$$I_{end}=\frac{\mu }{\gamma }{E}_{end}=\frac{a-K\left({\rho }_1+{\rho }_2\right)}{\beta \left(\mu +\gamma {\lambda }^e\right)}\mu$$

(46)

However, again, after cancelling $R_{k+1}=R_k=R_{end}$ in both sides of (4), one gets that the subsequent constraint should hold:

$$\gamma I_{end}+K\left({\rho }_1+{\rho }_2\right)S_{end}=\mu E_{end}+K\left({\rho }_1+{\rho }_2\right)S_{end}=0$$

(47)

which implies that $S_{end}=E_{end}=I_{end}=0$ (which only holds if $K=\frac{a}{{\rho }_1+{\rho }_2}$), and it is in fact a disease- free equilibrium point. In this case, one also has that $N_{end}=R_{end}$. Thus, Case 2 is unfeasible for the existence of an endemic equilibrium point.

Remark 3

(basic reproduction number). In biological terms, it is possible to re-interpret the condition of asymptotic stability around the disease-free equilibrium point in terms of the basic reproduction number, which indicates the number of secondary infectious individuals generated from one primary infectious one, defined by

$$R=\frac{\beta }{{\beta }_c}=\frac{\beta \left(\mu +{\lambda }^e\gamma \right)S_{df}}{\mu \gamma }$$

(48)

It turns out that$R<1$is exactly identical to$\beta <{\beta }_c$, that is, if the transmission rate is smaller than its critical value${\beta }_c$then, equivalently, the basic reproduction number is smaller than unity. That asymptotic stability condition, in the two mentioned equivalent forms, are related to the asymptotic extinction of the disease. It can be commented that such a condition is not related to the transient disease evolution but to its steady-state which is the disease-free equilibrium in this case. We can also to define a sample-dependent effective reproduction number as follows:

$$R_k=\frac{\beta \left(\mu +{\lambda }_k^e\gamma \right)S_k}{\mu \gamma }; k\in {\mathit{Z}}_{0+}$$

(49)

whose meaning and evolution should not be confused with that of the basic reproduction number. This number can be greater than one at the beginning of the disease transmission (even if$R<1$) but$R_k\to R$as$k\to \infty$. Note also from (1) that$R_k$decreases faster (since$S_k$decreases faster) if vaccination is programmed than in the vaccination-free case. The effective reproduction number is usually periodically checked to elucidate the particular intervention measures to be taken depending on the disease transmission evolution. It can be also pointed out that$R=1$is related, in common situations associated with epidemic models, to the confluence of the disease-free equilibrium point with the endemic one while, for$R<1$, the disease-free equilibrium point is unstable and the endemic equilibrium point is an attractor for the solutions. According to the former discussion in this subsection, this model, which is proposed for short-term predictions, does not evaluate the influence possible endemic steady states of the disease so that the central discussion is related to the case when$R<1$.

5. Simulation Results

This section contains some simulation examples devoted to study the application of the two-doses vaccination to counteract the spread of the COVID-19 pandemic. To this end the case of Italy borrowed from [45], where the discrete-model dynamics is confronted with actual data from Italy, will be considered. Thus, we have a discrete-time COVID-19 model serving as benchmark for the vaccination control. The model is parameterized by:

$$\beta =0.2,{\beta }^e=\frac{\beta }{1.3},\mu =\frac{1}{6},\gamma =0.04$$

so that ${\lambda }^e=0.7692$. The sampling time is one day so that the units of the parameters are in days⁻¹. The initial conditions are ${\mathrm{S}}_0=0.9999,{\mathrm{E}}_0=0.0001,{\mathrm{I}}_0=0,{\mathrm{R}}_0=0$ implying that the total population is normalized to unity, without loss of generality. Notice that almost all the population is susceptible and a small fraction of the population is exposed at the beginning. The parameters of the vaccination are ${\mathrm{d}}_1+{\mathrm{d}}_2=21$ days and ${\mathrm{d}}_2=7$ days; therefore, separation between the two doses is of two weeks. The values of the doses effectiveness are given by ${\mathsf{\rho }}_1=0.66$ and ${\rho }_2=0.3$ in such a way that the total effectiveness of the two doses $\rho ={\rho }_1+{\rho }_2$ is 96%, in accordance with the average effectiveness of available vaccines. The natural recruitment rate is a = 1, since the natural growth of the population may be rejected when it comes to the epidemic spreading description due to the small number of children affected. The Figure 1 displays the dynamics of the model without vaccination.

It is remarkable in Figure 1 how the number of infectious reaches a peak at around the 50% of the total population. Now vaccination is applied with a gain of K = 0.01, corresponding to the 1% of the susceptible vaccinated every day. The vaccination is applied starting at different moments of the epidemic spreading in order to observe the effect of vaccination in the epidemic dynamics. Thus, Figure 2 displays the effect of vaccination when it is applied since the beginning of the spreading.

It is observed from Figure 1 and Figure 2 how the application of vaccination reduces the peak of the infectious while increasing the pace at which susceptible become immune, as expected. In addition, the Figure 3 illustrates the effect of starting the vaccination at different moments of the spreading ranging from the beginning, 29 days and 59 days after the first day of simulation. It can be seen in Figure 3 how the effect of vaccination in spreading is higher the sooner vaccination is applied when cases are detected. Thus, vaccination is useful for preventing new outbreaks and reducing their severity. If vaccination is applied when there is a relatively large number of cases among the population, the peak is not reduced ostensibly. Figure 4 displays the values of vaccination actions $V_1$ and $V_2$ corresponding to this situation. It is deduced from Figure 4 that the feedback vaccination provides a lower number of vaccinated individuals when it starts at advanced stages of the spread due to a smaller number of susceptible individuals. Therefore, this kind of action is especially recommendable when the spread is at its initial steps. Figure 4 also shows that vaccination $V_1$ has larger values than $V_2$ due to the differences in effectiveness defined by the parameters ${\rho }_1$ and ${\rho }_2$.

Finally, Figure 5 displays the dynamics of the model for different values of the vaccination gain. As it could have been expected, the larger the vaccination action is, the smaller peak is attained. Naturally, the reduction in the peak is achieved at the expense of a higher vaccination effort as Figure 6 shows. It can also be observed in Figure 4 and Figure 6 the effect of delay in the control action, where a value of zero vaccination is provided during the first days of simulation due to the delay.

Figure 7 shows the effect of changing the interval of time between the application of the two doses. It is assumed that the effectivity does not change with the period in between them. The vaccination gain K = 0.01 is used while the initial conditions are those indicated at the beginning of this section. Moreover, d₂ = 7 days and d₁ + d₂ range from 10 to 21 days. It is observed in Figure 7 that if the second dose maintains effectivity regardless the dose sparing, it is better to administer it as soon as possible. However, in the real world, the effectivity of the second dose may depend on when this is applied. There are no data on how the effectivity depends on the dose sparing or these are scarce, [46]. Therefore a more accurate simulation on the effect of dose sparing is not carried out in this work.

Moreover, Figure 8 compares the evolution of the infectious when two doses with a spacing of 14 days is applied (K = 0.1) in contrast of applying a single dose with a double-vaccination gain (K = 2 × 0.01 = 0.02). The Figure 8 shows that with respect to the evolution of the infectious the administration of a single dose to a broader population alleviates the peak of infectious in comparison to promptly applying the second dose in order to increase the overall effectiveness of the vaccine. This behavior is behind the decision of the UK government of prioritizing the first dose to as many people as possible while delaying the second dose, [46].

The basic properties of the model discussed in Section 3 and Section 4 can also be corroborated in the previous examples. Thus, the number of exposed and infectious always converges to zero so that the model converges to the disease-free equilibrium point. This happens because there is no endemic point for this model as discussed in Section 4. The Figure 9 displays a longer simulation for the model with K = 0.001 and different initial conditions. It is shown in Figure 9 that the model always converges to the globally-stable disease-free equilibrium point given by $S_{df}=E_{df}=I_{df}=0,R_{df}=1$ and all the subpopulations remain non-negative at all time. These properties are indeed guaranteed by Theorems 1 and 2 since $\mu ,\gamma \in \left[0,1\right)$ and $\beta =0.2<1.66=\mathit{max}\left(\frac{a_k}{I_k+{\lambda }_k^e{E}_k}\right)$ is small enough; therefore, (15) holds, and the model trajectory solution is non-negative according to Theorem 1 and converges to the disease-free equilibrium point according to Theorem 2 since $K=0.001\ne 1.04=\frac{a}{{\rho }_1+{\rho }_2}$. In addition, Theorem 3 also ensures that the disease-free equilibrium point is globally stable since Condition 1 from Theorem 3 is satisfied, Condition 2 also holds as $K=0.001\ne 1.04=\frac{a}{{\rho }_1+{\rho }_2}$ and Equation (13) holds as Figure 10 depicts. Therefore, the disease-free equilibrium point is concluded to be globally asymptotically stable.

It can be pointed out that different alternative theoretic methods for stability analysis of dynamic systems in the presence of delays can be found, for instance, in [47,48,49] and some of the references therein.

6. Concluding Remarks

A new discrete SEIR model has been presented in this paper subject to two delayed doses of feedback vaccination controls on the susceptible. It is also assumed that there exists a transmission rate from the exposed to the susceptible and that the transmission rates exposed-susceptible and infectious-susceptible may eventually be distinct. The second idea is based on the knowledge that the later period in the asymptomatic incubation phase of the infection in the COVID-19 pandemic is also infective as it is the first phase of the infectious period. The first idea of using two vaccination doses is based on the administration mechanisms of some of the recently approved vaccines for COVID-19. Such two doses are applied with a certain interval, and the immunity effect is delayed in respect to each current picked- up sample on the model, while it is also assumed that both doses can have, in general, different effectiveness. It is assumed, in particular, that the second dose has an incremental benefit on the injection of the first single dose. The effectiveness degrees of both doses as well as the delays can be adjusted in the experimental tests since they are model parameters. The non-negativity properties of the solution under finite non-negative initial conditions, its boundedness and the disease-free equilibrium point as well as its stability properties are also investigated. The non-negativity of the solution can become lost for large-enough transmission rates, but it is kept for low and moderate values of the transmission rates. The proposed model has been numerically tested through numerical examples for COVID-19 parameterizations.