mathematics Article On a New Discrete SEIADR Model with Mixed Controls: Study of Its Properties Raul Nistal 1, *, Manuel de la Sen 1 , Santiago Alonso-Quesada 1 and Asier Ibeas 2 1 Department of Electricity and Electronics, University of Basque Country, UPV/EHU, 48940 Leioa, Spain; manuel.delasen@ehu.eus (M.d.l.S.); santiago.alonso@ehu.eus (S.A.-Q.) 2 Department of Telecommunications and Systems Engineering, Universitat Autonoma de Barcelona, 08193 Barcelona, Spain; asier.ibeas@uab.cat * Correspondence: raul.nistal@gmail.com  Received: 15 October 2018; Accepted: 20 December 2018; Published: 25 December 2018  Abstract: A new discrete SEIADR epidemic model is built based on previous continuous models. The model considers two extra subpopulation, namely, asymptomatic and lying corpses on the usual SEIR models. It can be of potential interest for diseases where infected corpses are infectious like, for instance, Ebola. The model includes two types of vaccinations, a constant one and another proportional to the susceptible subpopulation, as well as a treatment control applied to the infected subpopulation. We study the positivity of the controlled model and the stability of the equilibrium points. Simulations are made in order to provide allocation and examples to the different possible conditions. The equilibrium point with no infection and its stability is related, via the reproduction number values, to the reachability of the endemic equilibrium point. Keywords: vaccination controls; nonlinear dynamics and control systems; epidemic models; discrete models 1. Introduction We present an SEIADR epidemic model (Susceptible-Exposed-symptomatic Infectious- Asymptomatic Infectious- Dead infected-Recovered) based on typical descriptions of the spread of an infectious disease made in the background literature [1–6]. We will adapt the continuous time model from previous papers, such as [7], into a discrete model. We incorporate two subpopulations, namely, asymptomatic and dead infective corpses to the SEIR module. In this way, the model can be useful for diseases where the lying corpses are infective like, for example, Ebola virus disease, with eventual asymptomatic infectious. Our motivations for the discretization are an improvement in the manipulation of the dynamic equations and a major control in the adjustments of the sampling period, which will facilitate the adaptation of the model into a real situation, which this work is ultimately aimed at. The model has four types of infected subpopulations so that each of the possible stages of the disease is characterized by a different subpopulation. This makes the model interesting for its application for studies in the Ebola disease and similar diseases in which the re-emergences of the presence of infectious individuals cannot be explained otherwise [8]. The infected stage of the alive population is split into three different subpopulations, namely, the exposed, the symptomatic infectious and the asymptomatic infectious subpopulation, where an individual may be contagious even if he/she/it does not show the symptoms of the disease [9,10]. Furthermore, we have designed our model in order to take into account the deceased infectious individuals as a new infectious subpopulation, so we can prove their importance in dynamics of the population as the contagiousness is still present on some of the corpses, that must be disposed with great care [9,11–13]. Thus, the final stages of the illness under study can be split into live and immune individuals, i.e., recovered (R) subpopulation, and Mathematics 2019, 7, 18; doi:10.3390/math7010018 www.mdpi.com/journal/mathematics Mathematics 2019, 7, 18 2 of 19 infectious and dead individuals, i.e., dead (D) subpopulation. Given this, it has been widely studied the non negativity of the solutions [14–16] and the use of perturbations and nonlinear incidence rates (as seen in [17–19] and [19–22], respectively). Reducing the disease within the population should be main objective of the control measures proposed i.e., vaccination strategies and the medical treatment of the infectious subpopulation through antiviral treatment or care of the symptoms, which we will also include in our equations. This type of control strategies have been proposed in many models in the literature [3,9,14,15,23–27]. The stability analysis and control optimization of a a subpopulation of infectious individuals in treatment with vaccination is addressed in [28], and a time-delay model with a strategy of impulsive vaccination and a saturated incidence rate is addressed in [29]. For both the continuous and the discrete-time dynamic system with closed-loop stabilizing controllers, the properties of convergence of the state-trajectory to the equilibrium points and their global stability are very important tools for analyzing and designing those models. See, for instance [30–33], and some of the references therein. Therefore, we pay special attention to those properties in the proposed and studied controlled discrete epidemic model subject to a feedback vaccination control applied on the susceptible and to a feedback control on the infected subpopulation and to the relevance of the reproduction number, linked to the controller gains in the rates of convergence. Another relevant pillar of the study performed in this paper is the relevance of the discretization concerned to the positivity of the solution sequences and the attainability or un-attainability of the endemic equilibrium point depending on the value of the reproduction number. On the other hand, and concerned with the advantages statement of discrete epidemic models, we might point out that the performance of their continuous-time counterparts can be essentially kept if the discretization process is performed with good adjustments of the sampling period and the auxiliary discretization parameters in such a discrete-model. In particular, the computing time and the computer memory storage needs for control implementation can be reduced. A discrete epidemic model has been proposed, analyzed and discussed in [34] involving two interacting populations, the vector and the avian populations. The stability analysis and the stability properties of the equilibrium points of a pair of SIRS models subject to two viruses has been performed in [35]. Also, an important computational work has been implemented to illustrate the validity of the obtained results. On the other hand, the design of the stabilizing controllers have been recently invoked and developed with a well-worked mathematical rigor in [36,37], respectively. Finally, a study of the nonlinear phenomena of bifurcation and chaos has been detailed in [38] for a discrete SI epidemic model with fractional order while [39] gives a relevance of the fractional framework to the statement of fractional-type discrete epidemic models especially through the list of commented listed “ad hoc” references. In this paper, a discrete SEIADR model is proposed as a more generic variant of the SEIR one. A multiple set of possible control parameters is proposed, and all the different possible predictions about the population dynamics are studied through simulations in order to confirm the performed study. The main novelties of this paper are: • The characterization of the relations between the stability of the disease free equilibrium (DFE) point and the reachability of the endemic (END) one for the discrete SEIADR model under positive conditions. • The study of the stability and positivity properties and the equilibrium points and their properties. • The study through numerical examples of the influence of the controller gain in the equilibrium points and in the rates of convergence even under a similar reproduction number. The paper is organized as follows. Section 2 describes the SEIADR model with vaccination campaigns and antiviral treatments. A study of the disease free and the endemic equilibrium points is made. The discretization of the continuous model and the conditions for the positivity of the subpopulations are described in detail in Section 3. Section 4 gives and discusses the stability of the equilibrium points of this discrete model through the use of the next generation matrix applied to the disease free equilibrium point. Some results of simulations of the dynamics with a reproduction Mathematics 2019, 7, 18 3 of 19 number above and under one are presented at Section 5. Conclusions based on these results are presented in Section 6. 2. The Continuous SEIADR Model The dynamics of the SEIADR model is defined by the following equations:   βI(t) + β A A(t) + β D D(t) S(t) = b1 − V + ηR(t) − S(t) + b2 + KV , (1) N (t)   βI(t) + β A A(t) + β D D(t) E( t ) = S( t ) − (b2 + γ)E(t), (2) N (t) I(t) = pγE − (b2 + τ0 + Kξ + α)I(t), (3) A(t) = (1 − p)γE − (b2 + τ0 )A(t), (4) D(t) = αI(t) + (I(t) + A(t))b2 − µD(t), (5) R(t) = V + Kξ I(t) + KV S(t) + τ0 (A(t) + I(t)) − (b2 + η )R(t), (6) with N (t) = S(t) + E(t) + I(t) + A(t) + D(t) + R(t) being the total population, where S, E, I, A, D and R are the Susceptible, Exposed, Symptomatic Infectious, Asymptomatic Infectious, Dead Infectious corpses and Recovered subpopulations respectively. See [7] for a discussion of the continuous-time model. The parameters of the model are: • b1 is the recruitment/birth rate, • b2 is the natural average death rate, • β, β A , β D are the disease transmission coefficients from the susceptible to the symptomatic and asymptomatic infectious, and to the infective corpses subpopulations, respectively, • 1/η is the average duration of the immunity period which reflects a transition state from the recovered to the susceptible, • γ is the transition rate from the exposed to the symptomatic and asymptomatic infectious, • α is the extra average mortality being associated with the disease which affects to the symptomatic infectious subpopulation, • τ0 is the natural recovery rate for the whole infectious subpopulation (i.e., A + I ), • p is the exposed subpopulation fraction which becomes symptomatic infectious, • 1−p is the exposed subpopulation fraction which becomes asymptomatic infectious, • 1/µ is the average time of infectiousness after death, • V, KV and Kξ are the constant vaccination gain, the proportional vaccination control gain and the antiviral treatment control gain respectively. The constant vaccination is bounded such that V ∈ [0, b1 ], so a fraction V/b1 ∈ [0, 1] of the new individuals of the system (newborn, immigrants) is vaccinated. Remark 1. We define a set of auxiliary parameters to be then used as follows: BE = b2 + γ, BI = b2 + τ0 + Kξ + γ (1− p ) α, BA = b2 + τ0 , BR = b2 + η, f = BI , f A = BA , f D = f A bµ2 + f b2µ+α , G1 = BR + BE − b2 + µ f D BRb−b2 , γp 2 F = f β + f A β A + f D β D and G2 = ( BR + KV )(b2 f D − α f ). After using these auxiliary parameters in the Equations (1)–(6) the resulting equilibrium points of the solutions are: • A disease-free equilibrium point given by (SDFE , 0, 0, 0, 0, RDFE ) where BR NDFE − V K N +V SDFE = , RDFE = V DFE (7) BR + KV BR + KV Mathematics 2019, 7, 18 4 of 19 with the total population at this point defined as NDFE = bb1 . 2 • An endemic one given by (SEND , EEND , IEND , AEND , DEND , REND ) where: G2 SDFE + b2 G1 NDFE (S − SEND )( BR + KV ) SEND = BE , EEND = DFE , G2 BE + b2 FG1 G1 IEND = f EEND , AEND = f A EEND , DEND = f D EEND , (SDFE − SEND )( BE − (b2 + µ f D (1 + KV /b2 ))) REND = RDFE + , (8) G1 with the total population at this point given by: SEND F NEND = SEND + EEND + IEND + AEND + DEND + REND = . BE 3. Discretization of the SEIADR In some epidemic models, the vaccination and treatment control actions can be exerted with feedback information of the subpopulations, commonly the susceptible and infectious as the one presented in this paper. Since amounts of data on the subpopulation levels and parametrical computations to exert the controls through time have to be processed and monitored, the use of computer-controlled actions may be necessary. Generally speaking, the discretization of dynamic systems offers several potential computational advantages towards the controller synthesis, namely: 1. Normally feedback control actions are exerted by discrete-time controllers, especially, if the volume of data to be processed is relevant since the computational load has to be supported by a computer. Therefore, it can be preferred to start with a discrete-time model of the process, which then generates discrete sequences of measurable data, than a continuous-time one since then discrete control sequences are directly generated by processing the available sequence of discrete measurable data. Note also that the discretization of a continuous-time model towards the use of a computer for taking actions is always an approximation of the continuous time-model. However the sampling period of a discrete-time model is a design parameter, which does not imply an approximation when running the model. 2. It could be argued that in fact the use of a continuous-time model can be used for control generations through a computer but, in this case, the discretization period has to be very small in order to consider approximately valid the continuous-time control generation from continuous-time data. That is, there is no freedom to select the discretization sampling period. Note that if a discrete controller is accommodated to a discrete-time model then there is an important freedom in the choice of the sampling period, which takes the role of an extra control parameter which can be eventually time-varying, if it is compatible with the stability and bandwidth. 3. There is an important saving in data memory storage needs when implementing control actions, since only a discretized sequence of measurements needs to be stored and the control actions can be exerted along a set of time instants while the computer can exert alternative monitoring or computation actions. Thus, the computing time for control implementation is reduced. So, the use of a discrete-time feedback controller being accommodated to a discrete-time model does offer relevant advantages over related to a fast discretization of a closed-loop tandem of continuous-time and controller configuration in order to make an “ad hoc” computational control implementation. A discrete SEIADR model is formulated based on previous examples, as Mathematics 2019, 7, 18 5 of 19 those of [34,39,40], i.e., on the approximation of the subsequent derivative of the vector x(t) = (S(t), E(t), I(t), A(t), D(t), R(t)) x( t + T ) − x( t − T ) ẋ(t) → + a(x(t + T ) − 2x (t) + x(t − T )) (9) 2T where xk is defined as xk = x(kT ) so that the derivatives at the moments xk is implemented in the equation: x (1 + 2aT ) − 4aTxk − xk (1 − 2aT ) ẋk = k 2T The parameters “a” and T are real positive constants, being T the sampling period that we choose in order to adjust our model as desired, and “a” a modulating parameter introduced in order to guarantee the system’s positivity further described later on Proposition 1. Thus, the discrete equations result to be   ! 2T βIk + β A Ak + β D Dk Sk +1 = 1+2aT 2a − b2 + KV + Nk Sk + b1 − V + ηRk   (10) −2aT S + 11+ 2aT k −1 ,     2T Sk −2aT E Ek +1 = 1+2aT ( βIk + β A Ak + β D Dk ) N + ( 2a − BE )Ek + 11+ 2aT k −1 , (11) k   2T 1−2aT Ik +1 = 1+2aT (γpEk + (2a − BI )Ik ) + Ik −1 , 1+2aT (12)   1−2aT Ak+1 = 1+2T 2aT ( γ ( 1 − p ) E k + ( 2a − B A ) A k ) + 1+2aT Ak −1 , (13)   1−2aT Dk+1 = 1+2T2aT ( b2 Ak + ( b2 + α )Ik + (2a − µ )Dk ) + 1+2aT Dk −1 , (14)   1−2aT  Rk +1 = 1+2T 2aT V + K V S k + τ 0 ( A k + I k ) + K I ξ k + ( 2a − BR ) R k + 1+2aT Rk −1 , (15) Positivity of the Solution In this section, we will prove that the subpopulations affected by the disease always remain non-negative, provided that the initial condition is non-negative. The parameter “a” will be used to tune the discrete model so that its positivity is guaranteed. In the following, the notation for a vector x ≥ 0 means that all its components are non negative. We refer to a model with nonnegative solutions as a positive system. Proposition 1. Assume that b1 ≥ V and that for any k ≥ 1, xk ≥ 0, xk−1 ≥ 0 1 1 > a > max [ BE , BI , BR , µ, max [ β, β A , β D ] + b2 + KV ] (16) 2T 2 1 for which a necessary condition is T ≤ [ B ,B ,B ,µ,max[ β,β . Then xk+1 ≥ 0. As a consequence, if E I R A ,β D ]+ b2 + KV ] x0 ≥ 0 and x−1 ≥ 0 then xk ≥ 0; ∀k ≥ 0. Mathematics 2019, 7, 18 6 of 19 Proof. Assume that xk = (Sk , Ek , Ik , Ak , Dk , Rk ) ≥ 0. Then, from the discrete dynamic Equations (10)–(15) the positivity for the values of xk+1 = (Sk+1 , Ek+1 , Ik+1 , Ak+1 , Dk+1 , Rk+1 ) is studied as:   2T ( βIk + β A Ak + β D Dk ) Sk +1 ≥ 2a − (b2 + KV ) − Sk (17) 1 + 2aT Nk 2T = (2a − (b2 + KV ) − ( βni + β A n a + β D nd )) Sk (18) 1 + 2aT 2T ≥ (2a − (b2 + KV + max [ β, β A , β D ])) Sk ≥ 0, (19) 1 + 2aT I A D Since N + N + N = ni + n a + nd ≤ 1. Also 2T Ek +1 ≥ (2a − BE ) Ek ≥ 0, (20) 1 + 2aT 2T Ik +1 ≥ (2a − BI ) Ik ≥ 0, (21) 1 + 2aT 2T Ak +1 ≥ (2a − BA ) Ak ≥ 0, (22) 1 + 2aT 2T Dk +1 ≥ (2a − µ) Dk ≥ 0, (23) 1 + 2aT 2T Rk +1 ≥ (2a − BR ) Rk ≥ 0. (24) 1 + 2aT So xk+1 = (Sk+1 , Ek+1 , Ik+1 , Ak+1 , Dk+1 , Rk+1 ) is non negative provided that xk ≥ 0. As a result via compute induction if x0 ≥ 0 and x−1 ≥ 0 then xk ≥ 0, ∀k ≥ 0. Remark 2. Note that the modulating discretization parameter “a" can grow at most linearly with the half of the inverse of the sampling period. As a result of proposition, if x0 ≥ 0 and x−1 = 0 then xk ≥ 0 ∀k ≥ 0. 4. Equilibrium Points: Positivity and Stability From the discrete dynamic Equations (10)–(15), we calculate the equilibrium points such that xk+1 = xk = xk−1 leading to: ( βI∗ + β A A∗ + β D D∗ )     2T 1 + 2aT ∗ S∗ = − b2 + KV + S ∗ + b1 − V + ηR ∗ + S , (25) 1 + 2aT N∗ 2T S∗     2T 1 + 2aT E∗ = ( βI∗ + β A A∗ + β D D∗ ) ∗ + − BE E∗ , (26) 1 + 2aT N 2T     2T 1 + 2aT I∗ = γpE∗ + − BI I∗ , (27) 1 + 2aT 2T     2T 1 + 2aT A∗ = γ (1 − p )E∗ + − BA A∗ , (28) 1 + 2aT 2T     2T 1 + 2aT D∗ = b2 A∗ + (b2 + α) I∗ + − µ D∗ , (29) 1 + 2aT 2T     2T 1 + 2aT R∗ = V + KV S∗ + τ0 (A∗ + I∗ ) + Kξ I∗ + − BR R∗ . (30) 1 + 2aT 2T The resulting equilibrium points are equal to the continuous-time model counterpart. The disease free equilibrium point is defined as (SDFE , 0, 0, 0, 0, RDFE ), with SDFE and RDFE defined as in Equation (7), while the endemic equilibrium point (SEND , EEND , IEND , AEND , DEND , REND ) is defined as in Equation (8). Observe that, from Proposition 1, for a non-negative initial state, the system remains positive, so if SDFE − SEND < 0 the values at the endemic equilibrium of the infectious subpopulations, defined in Equation (8), are negative. Therefore, the endemic equilibrium is not reachable. Mathematics 2019, 7, 18 7 of 19 4.1. Local Asymptotic Stability of the DFE Point In this section, we will study the disease free equilibrium point (DFE) xDFE and obtain the Reproduction number, which would reveal us if it is locally stable or unstable. The stability of the DFE point is studied by constructing the next generation matrix. First we define the vector yk as yk = [ Ek , Ik , Ak , Dk , Ek −1 , Ik −1 , Ak −1 , Dk −1 , Sk , Rk , Sk −1 , Rk −1 ] and then linearize around the DFE point (xk = xk−1 = xDFE ), so we get the constant Jacobian matrix: ! ∂yk+1 F−Σ 0 J = = . ∂yk xDFE ∆ C We focus on the the 8 × 8 submatrix F − Σ, related to the infectious subpopulations E, I, A, and D, while F represents the appearance of new set of infections and Σ would describe the transitions between those infectious subpopulations. Thus, we define F as   0 β βA βD 0 0 0 0   0 0 0 0 0 0 0 0   0 0 0 0 0 0 0 0       2SDFE T  0 0 0 0 0 0 0 0  F=   NDFE (1 + 2aT )   0 0 0 0 0 0 0 0   0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0 and Σ as 2T ( BE −2a)   2aT −1 2aT +1 0 0 0 2aT +1 0 0 0  −2T pγ 2T ( BI −2a) 2aT −1    2aT +1 2aT +1 0 0 0 2aT +1 0 0   −2T (1− p)γ 2T ( BA −2a) 2aT −1   2aT +1 0 2aT +1 0 0 0 2aT +1 0    −2T (b2 +α) −2Tb2 2T (µ−2a) 2aT −1  Σ=  0 2aT +1 2aT +1 2aT +1 0 0 0 2aT +1 .    −1 0 0 0 0 0 0 0   0 −1 0 0 0 0 0 0     −1    0 0 0 0 0 0 0  0 0 0 −1 0 0 0 0 The next generation matrix (NGM) will be defined as: F ( I + Σ ) −1 = 1−2aT 1−2aT 1−2aT 1−2aT   a11 a12 a13 a14 1+2aT a11 1+2aT a12 1+2aT a13 1+2aT a14   0 0 0 0 0 0 0 0   0 0 0 0 0 0 0 0         0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0   0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0 with SDFE F SDFE ( βµ+ β D (b2 +α)) SDFE ( β A µ+ β D b2 ) SDFE β D a11 = NDFE BE , a12 = NDFE BI , a13 = NDFE BA µ , a14 = NDFE µ , Mathematics 2019, 7, 18 8 of 19 The spectral radius of the NGM, defined as the maximum of the absolute value of all possible eigenvalues of the NGM, ρ( F ( I + Σ)−1 ) will correspond to the reproduction number SDFE F S N R0 = = DFE END . (31) NDFE BE SEND NDFE We have concluded from the above discussion the following result: Theorem 1. The DFE point is locally asymptotically stable if R0 < 1, while if R0 > 1, then the DFE point is unstable. 4.2. Conditions of Positivity of the Equilibrium Points Here we will establish the conditions for the positivity, and thus, the potential reachability for any non-negative initial conditions, of the equilibrium points. Note that any non-negative solutions under Proposition 1 can reach an equilibrium point only if that has non-negative components. 4.2.1. DFE Point The following result is direct: Proposition 2. xDFE ≥ 0 if V ≤ b1 (1 + µ/b2 ). Proof. Note that IDFE = ADFE = EDFE = DDFE = 0. Also note that SDFE ≥ 0 if V ≤ b1 (1 + η/b2 ) from (7) and NDFE = b1 /b2 > 0. Then, again from (7), RDFE ≥ 0. As a result, xDFE ≥ 0. 4.2.2. END Point The positivity of the END point relies on the positivity of ∆S = (SDFE − SEND ) as it is seen in (8) to guarantee that EEND ≥ 0. Let us define: NDFE − NEND ∆N G2 x≡− = = . SDFE − SEND ∆S b2 G1 It is direct to see that the reproduction number, defined in (31), can be equivalently written as ∆S ∆S R0 (∆S, x ) = (1 + )(1 + x ), (32) SEND NDFE so if x ≥ 0, this implies that the reachability of the END point is given a reproduction number R0 > 1, as ∆S > 0 ⇔ R0 > 1, and ∆S < 0 ⇔ R0 < 1. The case x = 0 is not relevant for the analysis. b22 BI +b2 (ατ −b2 Kξ ) p Proposition 3. The non negativity of x is guaranteed if µ < BA αp . G2 Proof. The parameter x can be written as x = b2 G1 , so that the positivity of G1 is easily proven as G1 = b2 + η + γ + µ f D η/b2 > 0 b22 BI − p( BA αµ+b2 (b2 Kξ −τα)) Thus sign( G2 ) = sign(b2 f D − α f ) = sign(γ BA BI µ ), so b22 BI − p( BA αµ + b2 (b2 Kξ − τα)) γ > 0, i f BA BI µ b22 BI − p( BA αµ + b2 (b2 Kξ − τα)) > 0, b22 BI + b2 (ατ0 − b2 Kξ ) p G2 > 0 i f µ < . (33) BA αp Mathematics 2019, 7, 18 9 of 19 The subsequent result establishes that the existence of a reachable END point, which is defined by the positivity of x, will be also characterized by the reproduction number. Proposition 4. If b2 < µ, then the END point is reachable for any non-negative initial conditions if and only if R0 ≥ 1 ⇔ ∆S ≥ 0. If R0 = 1 then the END point is coincident with the DFE one. Proof. From (32), one gets: SDFE ( NDFE + x∆S) R0 (∆S, x ) = . (34) NDFE (SDFE − ∆S) We study first the solution for the critic value of the reproduction number R0 = 1. It follows from (32) that R0 = 1 if either ∆S = 0 or x = −SNDFE DFE . The derivative at the critic points will be d(R0 ) S (N + xSDFE ) R00 (∆S, x ) = = DFE DFE (35) d(∆S) NDFE (SDFE − ∆S)2 which gives us 1. R00 (∆S, − NDFE /SDFE ) = 0, ∀∆S ∈ R, and 2. for x 6= − N DFE 0 SDFE ⇒ R0 ( ∆S = 0, x ) = SDFE + 1 x NDFE so that if ( x > − NDFE /SDFE ⇔ R00 (∆S = 0, x ) > 0), then ∀∆S > 0 ⇒ R0 > 1, since R0 (∆S, x ) is non-decreasing with respect to ∆S for all real x > −N DFE SDFE . We know that x depends on the parameter p (i.e., the fraction of exposed which become symptomatic) as G2 ( p) x = x ( p) = (36) b2 G1 ( p) γ ( BR + KV ) − pBA (b2 (α + b2 ) − αµ) + b22 BI ( p − 1)  = , b2 BI γηµ( p − 1) − µBA (γη p (α + b2 ) + b2 BI (b2 + γ + η )) so that the derivative with respect to p is dx ( p) b2 (αγη +b2 (b2 +γ+η )(α+k ξ )) x0 ( p) = dp = (µ − b2 ) (αγη +αBµA( B(b2)+γ+η )) + µ ( BP ) (37) P with (b22 BA (α+ BA +kξ )+b2 ( BA2 (γ+η )+ BA (η (α+γ)+αγ+(γ+η )kξ )−γη ( p−1)(α+kξ ))+αγη pBA )2 BP = > 0. b2 γBA (b2 +η +KV )(α+ BA +k ξ ) Since µ > b2 , then x 0 ( p) > 0, ∀ p ∈ [0, 1], so the minimum value of x ( p) would be x (0) and then: −γb2 (η + b2 + KV ) x ( p ) ≥ x (0) = . (38) (b2 + τ )(b2 + η + γ)µ + γηµ Thus, if NDFE x ≥ x (0) > − , (39) SDFE (so that R0 (∆S, x ) is non-decreasing with respect to ∆S) ⇒ ∆S > 0 ⇔ R0 > 1. Mathematics 2019, 7, 18 10 of 19 The inequality from (39) becomes through Equation (7): NDFE b + η + KV (b2 + η + KV )γb2 /µ ≥ 2 > − x (0) = , (40) SDFE b2 + η γη + (b2 + τ )(b2 + η + γ) b2 /µ which can be reduced to 1 > , which is true if b2 < µ. 1 + b τ+τ + BγA 2 On the other hand, the value of the susceptible subpopulation in the END point is defined from (8) as N + xSDFE ( NDFE + xSDFE ) SEND = DFE = , (41) F/BE + x (R0 + xSDFE /NDFE ) NDFE /SDFE which will be positive if x > − NDFE /SDFE and R0 ≥ 1. Note that for R0 = 1 the endemic equilibrium point is coincident with the DFE one since ∆S = 0, so it is reachable. Then R00 (∆S) > 0 and the existence of the END point is characterized by R0 ≥ 1. 4.3. Global Stability We know that the reproduction number of (31) can be rewritten as      η p(α+b2 ) β D + βµp ( p−1)(µβ S +b2 β D ) b2 γ b1 b + 1 − V α+b +k +τ − b2 +τ0 2 2 ξ 0 R0 = (42) b1 µ (b2 + γ) (b2 + η + KV ) and by using the relative disease transmission coefficients β Ar = β A /β and β Dr = β D /β we can rewrite the reproduction number as      η p(α+b2 ) β Dr +µp ( p−1)(µβ Ar +b2 β Dr ) b2 γ b1 b + 1 − V α+b +k +τ − b2 +τ0 2 2 ξ 0 R0 = β . (43) b1 µ (b2 + γ) (b2 + η + KV ) Furthermore, the Jacobian matrix J on the equilibrium points is equal to ! A B J= (44) I6 0 1−2aT 2T with B = 1+2aT I6 and A = 1+2aT Q, being   q11 q12 q13 q14 q15 q16   q21 q22 q23 q24 q25 q26    0 pγ 2a − (b2 + τ0 + k ξ + α) 0 0 0  Q= ,    0 γ − pγ 0 2a − (b2 + τ0 ) 0 0     0 0 α + b2 b2 2a − µ 0  KV 0 k ξ + τ0 τ0 0 2a − η − b2 with BE E∗ ( N ∗ −S∗ ) BE E∗ q11 = 2a − b2 − KV − N ∗ S∗ q12 = N∗  BE E∗ − βS∗  q13 = N∗ q21 = BE E∗ S1∗ − N1∗ BE E∗ βS∗ − BE E∗ q22 = 2a − γ − b2 − N∗ q23 = N∗ BE E∗ −S∗ β A BE E∗ −S∗ β D (45) q14 = N∗ q15 = N∗ ∗ β A S∗ − BE E∗ q16 = η + BNE E∗ q24 = N∗ β D S − BE E∗ ∗ ∗ q25 = N∗ q26 = − BNE E∗ Mathematics 2019, 7, 18 11 of 19 and { N ∗ , S∗ , E∗ } → { NDFE , SDFE , 0} for the Jacobian matrix associated to the DFE point JDFE , and { N ∗ , S∗ , E∗ } → { NEND , SEND , EEND } for the Jacobian matrix associated to the END point JEND . The matrix Q corresponding to JDEF , namely, QDFE of the JDFE will then be equal to: β A SDFE D SDFE   2a − b2 − KV 0 − βS DFE NDFE − NDFE − βN DFE η  βSDFE β A SDFE β D SDFE    0 2a − γ − b2 NDFE NDFE NDFE 0   QDFE  = 0 pγ 2a − (b2 + τ0 + k ξ + α) 0 0   0 (1 − p ) γ 2a − (b2 + τ0 ) 0 0     2a − µ    0 0 α + b2 b2 0  KV 0 k ξ + τ0 τ0 0 2a − (η + b2 ) βSDFE β A SDFE β D SDFE   b2 + KV 0 NDFE NDFE NDFE −η β A SDFE D SDFE − βSDFE − βN   0 γ + b2 NDFE − NDFE 0    DFE  QDFE  = 2aI6 −  0 − pγ b2 + τ0 + k ξ + α 0 0   (46) 0 ( p − 1) γ b2 + τ0 0 0     −α − b2 − b2    0 0 µ 0  − KV 0 −k ξ − τ0 − τ0 0 η + b2 Assume that (16) hold. Then, the following result guarantees the global stability under positivity conditions of Proposition 1. Theorem 2. The following properties hold: (i) The total population Nk ; ∀k ≥ 0 is positive and bounded for any given non-negative initial conditions. (ii) The discrete SEIADR epidemic model is globally Lyapunov’s stable for any given finite non-negative initial conditions irrespective of the value of the reproduction number. (iii) if R0 ≤ 1 then the DFE point is the unique reachable equilibrium which is globally asymptotically stable. Proof. Since Nk = SK + EK + IK + AK + DK + RK ; ∀k ≥ 0, subject to any finite initial conditions satisfying min(S0 , E0 , I0 , A0 , D0 , R0 ) ≥ 0, it turns out that Nk ≥ 0; ∀k ≥ 0 since all the subpopulations are non-negative for any sample from Proposition 1. Now, by summing up (10)–(15), one gets the following evolution equation for the total population: 1 − 2aT 2T Nk+1 = N + [2aNk + b1 − µDK − b2 (EK + RK + SK )]; ∀k ≥ 1 (47) 1 + 2aT k−1 1 + 2aT One gets from (47) that: 2Tν 2Tb1 N̂k+1 + (DK + EK + RK + SK )e1  A N N̂k + e ; ∀k ≥ 1 (48) 1 + 2aT 1 + 2aT 1 where ν = min(b2 , µ) and ! 4aT/(1 + 2aT ) (1 − 2aT )/(1 + 2aT ) AN = , (49) 1 0 with N̂k = [ Nk Nk−1 ] T ; e1 = (1, 0) T and the notation “” in M = ( Mij )  N = ( Nij ) with M and N being of the same order means that Mij ≤ Nij ; ∀i, j, that is, the inequalities holds for all the corresponding matrix entries of M and N. It turns out by direct inspection that A N is positive (in the sense that it has no negative entries) and convergent (i.e., stable in the discrete sense) since its eigenvalues are within the unity ratio circle in the complex plane. Now, assume that the sequence { N̂k }∞k =0 is unbounded to prove its boundedness by contradiction arguments. Then, Mathematics 2019, 7, 18 12 of 19 it has a subsequence { N̂kn }∞ ∞ ∞ n=0 ⊆ { N̂k }k =0 which is strictly increasing with { k n }n=0 being a strictly ∞ increasing sequence of non-negative integers, that is, { N̂kn }n=0 → ∞ as n → ∞ and one gets via recursive calculations from (48) and Proposition 1 that: k −1 0  N̂kn+1 + Σi=n+k1n (A N )(kn+1 −(i+1)) (Di + Ei + Ri + Si )e1 , (50) k −1 2Tb1  (A N )(kn+1 −kn ) N̂kn + Σi=n+k1n (A N )(kn+1 −(i+1)) e ; ∀k ≥ 1, (51) 1 + 2aT 1 since A N is a convergent matrix so that spectral radius r (A N ) = in f kA N k = |λmax (A N )| < 1 , A N  0 and { N̂kn }∞ n=0 is a non-negative strictly increasing sequence. Then one has from (51): k −1 k N̂kn+m+1 + Σi=n+kmn +1 (A N )(kn+m+1 −(i+1)) (Di + Ei + Ri + Si )e1 k k −1 ≤ k(A N )(kn+m+1 −kn ) N̂kn k + kΣi=n+kmn +1 (A N )(kn+m+1 −(i+1)) 12Tb +2aT e1 k; ∀ k ≥ 1 1 2Tb1 (1−ρk n+m+1 −k n ) ≤ ρkn+m−1 −kn k N̂kn k + (1+2aT )ρ ∀m ≥ 0, ∀k ≥ 1 (52) for some matrix norm k•k such that kA N k ≤ ρ = r (A N ) + e < 1 for some given real constant e fulfilling 0 ≤ e < 1 − r (A N ). From the definition of the spectral radius and the fact that A N is convergent such a norm upper-bounded by such a ρ parameter always exists. By talking limits in (52) as n, m → ∞, so that (k n+m+1 − k n ) → ∞, one gets the contradiction +∞ ≤ (1+2Tb1 2aT )ρ . This happens since N̂kn+m+1 → +∞ as (n + m) → ∞ by the unboundedness assumption on Nk as k → ∞ and ρ < 1, DK + EK + RK + SK ≥ 0 and Nk ≥ 0; ∀k ≥ 0 from the non-negativity of the subpopulations with the conditions of Proposition 1. Then, Nk ; ∀k ≥ 0 is bounded for any given non-negative initial conditions and Property (i) follows. Property (ii) is a direct consequence of Property (i) and the non-negativity of all the subpopulations for any sample and given finite non-negative initial conditions. Property (iii) follows since if R0 < 1 the END point is not reachable and if R0 = 1 it is coincident with the DFE one. Therefore for R0 ≤ 1 the unique reachable equilibrium point from non-negative initial conditions is the DFE one which is globally asymptotically stable since it is locally asymptotically stable from Theorem 1. 5. Numerical Simulations We start the simulation with the subpopulations in an initial conditions next to the DFE point, concretely: ∗ ∗ (S(0), E(0), I(0), A(0), D(0), R(0)) = (SDFE , 0.05, 0.05, 0.05, 0.05, RDFE ) With the susceptible and recovered DFE values as defined in Equation (7). Then, our parameters are set in order to obtain a reproduction number less than unity. The value of the parameters relative to the host population are based in what is expected to be the average rates of growth in a human society, and the ones related to the disease are based on previous models describing Ebola disease [7]. In this way, the mortality and natality rates will be b1 = b2 = 1/(365 × 70) days−1 , while the recovery rates from the infection are τ0 = 1/12 days−1 , µ = 1/3 days−1 , the parameters of the disease are defined by the infectivity rates β = 0.15 days−1 , β A = 0.02 days−1 and β D = 0.15 days−1 , the probability of transition p = 0.75, the mortality of the infected subpopulation α = 0.1 days−1 , the transition from infected to recovered γ = 1/300 days−1 and the recovered to susceptible rate is defined by η = 0.15 days−1 . On the other hand, the parameters for the controller methods are defined by the vaccination constants V = b1 /100 days−1 , KV = b2 /2 days−1 , and the antiviral constant Kξ = 0.01 days−1 . Given these conditions, we obtain a reproduction number smaller than unity 1, namely, R0 = 0.81 and the initial conditions based on the DFE point RDFE = 1.3 × 10−4 and SDFE = 0.99987. It can be seen in Figures 1 and 2 the time evolution of the subpopulations given these parameters: The infectious subpopulations decrease exponentially as the susceptible and recovered Mathematics 2019, 7, 18 13 of 19 subpopulations continuously approaches to the values corresponding the DFE point. Note that the simulations are performed to display the trajectory evolution tendencies through time rather than describe real amount of individuals. Figure 1. Time evolution of the different subpopulations for a R0 = 0.81. Another simulation is made with the subpopulations starting at the same initial condition that the used in the previous example, but in this time with the reproduction number exceeding unity. The parameters, as in the previous simulation, are based on the average rates of growth in a human society afflicted by an Ebola-type disease. However, the value of the main infectivity rate β is changed from β = 0.15 days−1 to β = 0.30 days−1 so the reproduction number goes from R0 = 0.81 to R0 = 1.39. Observe that the change of the infectibity does not change the values of the parameters SDFE and RDFE , so the initial conditions remain the same. We can see at Figures 3 and 4 the time evolution of the subpopulations reaching the END point. Observe in Figures 4 and 5 that even though given the parameters used in this simulation system takes quite long to stabilize, the subpopulations tend to the endemic values (red dotted). It is worth mentioning that the purpose of running the simulations at such a long time is not with the intention of showing the desired final state of our population, as it would be ridiculous. We want however, to show the general tendency of the dynamic of the subpopulations when the reproduction number exceeds one ( the infectious subpopulations will be relevant for an indefinite time) and when it doesn’t (the infectious subpopulations will decrease exponentially over time). More numerical experimentation has been performed by modifying the control gains while maintaining the reproduction number. The purpose is to study the different velocities at which the subpopulations reach the equilibrium point for a given reproduction number. In this case, the stability properties of the DFE point hold intact depending on the value of R0 being smaller or greater than unity, but the values of the subpopulations corresponding to the equilibrium points of the model and the rate of convergence can be modified. In Figure 6 we can see the dynamics characterized by the parameters of the previous simulations for the infectivity rates β = 0.92 days−1 , β A = 0.2 and β D = 0.5, and two different sets of vaccination control gains, one corresponding to the constant vaccination strategy and the other one corresponding to the vaccination strategy based on the feedback of the susceptible subpopulation. Both present the same reproduction number and the same values for the subpopulations in the DFE point. However, it is seen that the disease decreases exponentially in both examples, as it is expected when the disease is being eradicated, and that there is a noticeable difference in the presence of the sick individuals. The infectious subpopulations for the set of parameters 2 are larger during the first days, while in the long run they decrease faster than the ones from set 1. Mathematics 2019, 7, 18 14 of 19 Figure 2. A detailed display of the reaction of the different subpopulations during the first days after the introduction of the disease given a R0 = 0.81. Figure 3. Time evolution of the different subpopulations given a R0 = 1.39. Mathematics 2019, 7, 18 15 of 19 Figure 4. A better display of the different subpopulations for a R0 = 1.39. Figure 5. A close look at the beginning of the dynamic simulation of the different subpopulations for a R0 = 1.39. Mathematics 2019, 7, 18 16 of 19 Figure 6. Time evolution of the different subpopulations given a R0 = 0.81 and two different sets of control vaccination gains. Set 1 shares common parameters with set 2 except for V = 0 and KV = 0.89, while set 2 presents V = 0.13 and KV = 0. In Figure 7 we show a similar disposition than the previous graphic. With the same parameters as before and V = 0, we run three different simulations with the same parameters and reproduction number R0 = 0.81, varying the proportional vaccination control gain KV and the antiviral treatment control gain Kξ . Figure 7. Time evolution of the different subpopulations given a R0 = 0.81 and three different sets of control vaccination gains. Set 1 shares common parameters with set 2 and 3, except for KV = 0 and Kξ = 1.29, while set 2 presents KV = 0.45 and Kξ = 0.14 and set 3 KV = 0.90 and Kξ = 0. 6. Conclusions A functional discrete model based on classic models [1,2] describing an Ebola-type disease spreading through a population has been constructed and the control gains applied to synthesize the vaccination and treatment signals are also put into the model successfully. The discrete model has been properly created and the boundaries of the parameters of the simulation are set so that the non-negativity of the system remains as long as the initial conditions are non-negative, as it should be in a model that depicts real populations. Thus, the conditions that guarantee the boundaries of the total population stability of the END and the DFE points have been tested and studied. The simulations presented in this paper also show that the values assigned to the control gain can substantially affect Mathematics 2019, 7, 18 17 of 19 the evolution of the infection. This happens even though the final equilibrium state is the same, which can be used to pick the best strategy when facing humanitarian crisis situations, or when trying to eradicate an endemic disease happening in a region. It can be seen that the increase of the control gains V, Kv and Kξ leads to a computable decrease of the reproduction number (see Equations (31) and (42). It can be seen that the DFE point can be modified with the vaccination control gains, in particular, the susceptible decrease with the constant vaccination V, while the recovered increase with such a control. It can be seen that the disease free susceptible subpopulation decreases as KV increases, while the recovered increases since the total population remains constant. Also, we have observed the influence of the vaccination control gains in the dynamics even when the reproduction number stays invariant under certain changes of the control gains. This suggests an interesting point when building the control gains of the epidemic models in future works, as it shows that the reproduction number is not the only relevant parameter to be taken into account. Author Contributions: Conceptualization, R.N. and M.d.l.S.; methodology, R.N., M.d.l.S. and S.A.-Q.; software, R.N.; validation, S.A.-Q. and M.d.l.S.; formal analysis, M.d.l.S., S.A.-Q., A.I. and R.N.; investigation, R.N. and M.d.l.S.; resources, M.d.l.S. and R.N.; data curation, M.d.l.S. and R.N.; writing—original draft preparation, R.N.; writing—review and editing, S.A.-Q., A.I. and M.d.l.S.; visualization, R.N.; supervision, M.d.l.S., A.I. and S.A.-Q.; project administration, M.d.l.S.; funding acquisition, M.d.l.S. Funding: This research was funded by the Spanish Government and by the European Fund of Regional Development FEDER through Grant DPI2015-64766-R and by UPV/EHU by Grant PGC 17/33. Conflicts of Interest: The authors declare that there is no conflict of interests regarding this manuscript. References 1. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser. A 1927, 115, 700–721. [CrossRef] 2. Hethcote, H.W. The Mathematics of Infectious Diseases. SIAM Rev. 2000, 42, 599–653. [CrossRef] 3. Keeling, M.J.; Rohani, P. Modeling Infectious Diseases in Humans and Animals; Princeton University Press: Princeton, NJ, USA; Oxford, UK, 2008. 4. Hethcote, H.W. Three basic epidemiological models. Appl. Math. Ecol. 1989, 18, 119–144. 5. Daley, D.; Gani, J. Epidemic Modeling: An Introduction; Cambridge Studies in Mathematical Biology: 15; Cambridge University Press: New York, NY, USA, 2005. 6. Khan, H.; Mohapatra, R.; Vajravelu, K.; Liao, S.J. The explicit series solution of SIR and SIS epidemic models. Appl. Math. Comput. 2009, 215, 653–669. [CrossRef] 7. Nistal, R.; De la Sen, M.; Alonso-Quesada, S.; Ibeas, A. A supervised multi-control for monitoring the antiviral treatment strategy for an SEIADR epidemic model. In Proceedings of the 2018 5th International Conference on Control, Decision and Information Technologies (CoDIT), Thessaloniki, Greece, 10–13 April 2018. 8. Glynn, J.; Bower, H.; Johnson, S.; Houlihan, C.; Montesano, C.; Scott, J.T.; Semple, M.G.; Bangura, M.S.; Kamara, A.J.; Kamara, O.; et al. Asymptotic infection and unrecognised Ebola virus disease in Ebola-affected households in Sierra Leone: A cross-sectional study using new non-invasive assay for antibodies to Ebola virus. Lancet Infect. Dis. 2017, 6, 645–653. [CrossRef] 9. de Pinho, M.; Maurer, H.; Kornienko, I. Optimal control of a SEIR model with mixed constraints and L1 cost. In Proceedings of the 11th Portuguese Conference on Automatic Control Lecture Notes in Electrical Engineering, Porto, Portugal, 21–23 July 2014. 10. Leroy, E.; Baize, S.; Volchkov, V.; Fisher-Hoch, S.P.; Georges-Courbot, M.C.; Lansoud-Soukate, J.; Capron, M.; Debre, P.; Georges, A.J.; McCormick, J.B. Human asymptomatic Ebola infection and strong inflammatory response. Lancet 2000, 355, 2210–2215. [CrossRef] 11. Bellan, S.; Pulliam, J.; Dushoff, J.; Meyers, L. Ebola control: Effect of asymptomatic infection and acquired immunity. Lancet 2014, 384, 1499–1500. [CrossRef] 12. Santermans, E.; Robesyn, E.; Sudre, T.G.B.; Faes, C.; Quinten, C.; Van Bortel, W.; Haber, T.; Kovac, T.; Reeth, F.V.; Testa, M.; et al. Spatiotemporal evolution of Ebola disease at sub-national level during the 2014 West Africa epidemic: Model scrutiny and data meagerness. PLoS ONE 2016, 11, e0147172. [CrossRef] Mathematics 2019, 7, 18 18 of 19 13. Al-Darabsah, I.; Yuan, Y. A time-delayed epidemic model for Ebola disease transmission. Appl. Math. Comput. 2016, 290, 307–325. [CrossRef] 14. De la Sen, M.; Agarwal, R.P.; Ibeas, A.; Alonso-Quesada, S. On the existence of equilibrium points, boundedness, oscillating behaviour and positivity of a SVEIRS epidemic model under constant and impulsive vaccination. Adv. Differ. Equ. 2011, 2011, 748608. [CrossRef] 15. De la Sen, M.; Alonso-Quesada, S. Vaccination strategies based on feedback control techniques for a SEIR-epidemic model. Appl. Math. Comput. 2011, 218, 3888–3904. [CrossRef] 16. Wei, Z.; Le, M. Existence and Convergence of the Positive Solutions of a Discrete Epidemic Model. Discret. Dyn. Nat. Soc. 2015, 2015, 434537. [CrossRef] 17. Wang, L.; Liu, Z.; Zhang, X. Global dynamics of an SVEIR epidemic model with distributed delay and nonlinear incidence. Appl. Math. Comput. 2016, 2016, 47–65. [CrossRef] 18. Wang, X. An SIRS Epidemic Model with Vital Dynamics and a Ratio-Dependent Saturation Incidence Rate. Discret. Dyn. Nat. Soc. 2015, 2015, 720682. [CrossRef] 19. Fengying, W.; Chen, F. Stochastic permanence of an SIQS epidemic model with saturated incidence and independent random perturbations. Phys. A Stat. Mech. Its Appl. 2016, 453, 99–107. 20. De la Sen, M.; Alonso-Quesada, S.; Ibeas, A. On the stability of an SEIR epidemic model with distributed time-delay and a general class of feedback vaccination rules. Appl. Math. Comput. 2015, 270, 953–976. [CrossRef] 21. Shaikhet, L. Stability of equilibrium states for a stochastically perturbed exponential type system of differential equations. J. Comput. Appl. Math. 2015, 290, 92–103. [CrossRef] 22. Shaikhet, L.; Korobeinikov, A. Stability of a stochastic model for HIV-1 dynamics within a host. Appl. Anal. 2016, 95, 1228–1238. [CrossRef] 23. De la Sen, M.; Agarwal, R.P.; Ibeas, A.; Alonso-Quesada, S. On a generalized time-varying SEIR epidemic model with mixed point and distributed time-varying delays and combined regular and impulsive vaccination. Adv. Diff. Equ. 2010, 2010, 281612. [CrossRef] 24. Tripathi, J.; Abbas, S. Global dynamics of autonomous and nonautonomous SI epidemic models with nonlinear incidence rate and feedback controls. Nonlinear Dyn. 2016, 86, 337–351. [CrossRef] 25. Buonomo, B.; Lacitignola, D.; de Leon, C.V. Qualitative analysis and optimal control of an epidemic model with vaccination and treatment. Math. Comput. Simul. 2014, 100, 88–102. [CrossRef] 26. Ling, L.; Jiang, G.; Long, T. The dynamics of an SIS epidemic model with fixed-time birth pulses and state feedback pulse treatments. Appl. Math. Model. 2015, 39, 5579–5591. [CrossRef] 27. He, Y.; Gao, S.; Xie, D. An SIR epidemic model with time-varying pulse control schemes and saturated infectious force. Appl. Math. Model. 2013, 37, 8131–8140. [CrossRef] 28. Sharma, S.; Samanta, G. Stability analysis and optimal control of an epidemic model with vaccination. Int. J. Biomath. 2015, 8, 1550030. [CrossRef] 29. Samanta, G. A delayed hand-foot-mouth disease model with pulse vaccination strategy. Comput. Appl. Math. 2015, 34, 1131–1152. [CrossRef] 30. De la Sen, M.; Ibeas, A. On the global asymptotic stability of switched linear time-varying systems with constant point delays. Discret. Dyn. Nat. Soc. 2008, 2008, 231710. [CrossRef] 31. Bilbao-Guillerna, A.; De la Sen, M.; Ibeas, A.; Alonso-Quesada, S. Robustly stable multiestimation scheme for adaptive control and identification. Discret. Dyn. Nat. Soc. 2005, 2015, 31–67. [CrossRef] 32. Bilbao-Guillerna, A.; De la Sen, M.; Alonso-Quesada, S. Multimodel discrete control with online updating of the fractional order hold gains. Cybern. Syst. 2007, 38, 249–274. [CrossRef] 33. Herrera, J.; Ibeas, A.; Alcantara, S.; De la Sen, M. Multimodel-based techniques for the identification and adaptive control of delayed multi-input multi-output systems. IET Control Theory Appl. 2011, 5, 188–202. [CrossRef] 34. Jang, S. On a discrete west Nile epidemic model. Comput. Appl. Math. 2007, 6, 397–414. [CrossRef] 35. Zhao, J.; Wang, L.; Han, Z. Stability analysis of two new SIRS models with two viruses. Int. J. Comput. Math. 2018, 95, 2026–2035. [CrossRef] 36. Degue, K.; Ny, J.L. An interval observer for discrete-time SEIR epidemic model. In Proceedings of the 2018 American Control Conference (ACC), Milwaukee, WI, USA, 27–29 June 2018; pp. 5934–5939. Mathematics 2019, 7, 18 19 of 19 37. Chang, C.; Jing, Y.; Zhu, B. Modeling and control for a descriptor epidemic system with nonlinear incidence rate. In Proceedings of the 2018 Chinese Control Conference (CCDC), Shenyang, China, 27–29 June 2018; pp. 2183–2188. 38. Abdelaziz, M.; Ismail, A.I.; Abdullah, F.; Mohd, M.H. Bifurcations and chaos in a discrete SI epidemic model with fractional order. Adv. Diff. Equ. 2018, 44, 1–19. [CrossRef] 39. Chiranjeevi, T.; Biswas, B. Discrete-time fractional optimal control. Mathematics 2017, 5, 25. [CrossRef] 40. De la Sen, M. Preserving positive realness through discretization. Positivity 2002, 6, 31–45. [CrossRef] c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).