Stability analysis and observer design for discrete-time SEIR epidemic models

This paper applies Micken’s discretization method to obtain a discrete-time SEIR epidemic model. The positivity of the model along with the existence and stability of equilibrium points is discussed for the discrete-time case. Afterwards, the design of a state observer for this discrete-time SEIR epidemic model is tackled. The analysis of the model along with the observer design is faced in an implicit way instead of obtaining first an explicit formulation of the system which is the novelty of the presented approach. Moreover, some sufficient conditions to ensure the asymptotic stability of the observer are provided in terms of a matrix inequality that can be cast in the form of a LMI. The feasibility of the matrix inequality is proved, while some simulation examples show the operation and usefulness of the observer.


Introduction
Mathematical models have become an important tool in analyzing the causes, dynamics, and spread of epidemics []. Thus, its study is crucial in order to obtain valuable knowledge of the underlying aspects of the disease. Moreover, the analysis of mathematical models describing epidemics spreading allows authorities to make decisions regarding the vaccination policies, quarantine application and so on. These models can be classified into two main categories: continuous-time and discrete-time models.
Continuous-time models were the first proposed to describe epidemics spreading since modeling has traditionally been focused from a differential equations point of view. Therefore, many specific features regarding these models have been studied in many works such as the presence of bifurcations [, ], the construction of approximate solutions [, ], the fractional dynamics [], the existence of equilibrium points and oscillatory behavior [], and the presence of waves [-], to cite but a few. However, model stability has been by far the most important property to be studied [-].
On the other hand, discrete-time models have gained substantial importance during the last years as the increasing number of papers concerning this type of models reveals [-]. There are many reasons that explain this tendency. Among them we can find the simplicity of its simulation, the fact that it is easier to adjust system's parameters from statistical data in discrete-time models than in continuous-time ones and the fact that discrete-time models may exhibit a richer dynamic behavior than its continuous-time counterparts []. Until now, as [] states, the study of discrete epidemic models was focused on the computation of the reproduction number, the existence and stability (both local and global) of the equilibrium points, and the extinction, permanence, and persistence of the disease. In addition, these studies have been performed basically on SI, SIR, and SIRS models.
Moreover, there are two ways to face the problem of setting up a discrete-time epidemic model. The first one consists in modeling the disease spreading directly in a discrete domain. The second one faces the problem by discretizing a previously existing continuoustime model. The latter approach will be used in this paper. In this way, several discretization methods have been used in the literature including the implicit and explicit Euler, mixed, and Runge-Kutta methods [-]. However, discretization procedures do not always preserve the structural properties exhibited by continuous-time models such as the conservation of the total population or the positivity of the solutions. Therefore, much research has been done in the direction of obtaining discretization methods able to preserve all these properties. One of the most important breakthroughs in this respect is the research made by Micken on the use of nonstandard finite difference (NSFD) methods []. These methods combine implicit and explicit Euler approximations in the incidence rate function. Therefore, this paper will apply Micken's NSFD discretization method to obtain a SEIR discrete-time epidemic model while stating its main properties such as the well-posedness of the model, the positivity of solutions and the existence and stability of equilibrium points. Since most of the previous literature is confined to SI, SIR, and SIRS models, this study will contribute to the application of NSFD approaches to a wider class of models. In addition, the properties will be proved based on an implicit representation of the equations, in comparison with previous approaches [, ], where an explicit description of the model is obtained prior to analysis.
Furthermore, control theory has recently emerged as one of the proposed approaches to deal with the design of vaccination campaigns [-]. Thus, different types of vaccination laws based on control theory have appeared in the literature during the last years, such as the state-feedback control [, ], robust sliding control [], or feedback linearization []. Nevertheless, most of these works calculate the vaccination effort based on the values of the subpopulations, namely, susceptible, exposed, infectious and immune of the SEIR model. From a practical point of view only the total population and the number of infectious may be obtained from clinical data at health centers. Therefore, the number of susceptible, exposed, and immune cases may not be available to on-line compute the control vaccination law. Under these circumstances, a state observer is necessary to estimate the value of these subpopulations.
State observers are largely used in control applications to cope with the inability to measure some state variables. Thus, we can find many applications in chemical [], aeronautic [], and electrical engineering [], to cite just a few. Also, there are few previous works where a state observer is designed for epidemic models [, ]. In these works, an observer is designed for a continuous-time SEIR epidemic model whose values are then used to implement a vaccination strategy. In the present work, however, a state observer is designed for a discrete-time SEIR model obtained by a NSFD discretization method, which is for the first time being considered in the literature. The design of such a state observer must face some particularities, which make it different from previous approaches to observer design: • NSFD schemes lead to implicit systems rather than explicit. Thus, the design of the state observer must be done implicitly. This makes a big difference from traditional observer design where the problem is usually formulated in an explicit way. • Epidemic models are singular systems so that there exists an algebraic relation between the dynamic variables of the system. • The number of infectious can be measured, implying that there is no need to estimate this state component, which leads to a reduced-order observer. In this way, the combination of all the above items makes the design of state observers for discrete-time SEIR epidemic models rather different from standard approaches in the literature. In addition, the asymptotic stability of the observer will be discussed, while some simulation examples will show its operation. Some preliminary results of this paper were first presented in [] and [].
The paper is organized as follows. Section  describes the model description and problem formulation. Section  studies the main properties of the discrete model. The design of the state observer is carried out in Section . Section  is devoted to the study of the analytic properties of the observer. Finally, Section  shows some numerical simulation examples while our Conclusions end the paper.

Model description
Consider the SEIR continuous-time epidemic model described by [] where S(t), E(t), I(t), and R(t) denote the subpopulations of susceptible, exposed, infectious and immune, respectively, defining the state vector x T (t) = [S E I R]. N(t) denotes the total population at time t (i.e. N(t) = S(t) + E(t) + I(t) + R(t)), μ is the rate of deaths from causes unrelated to the infection, ν denotes the birth rate and ω is the rate of losing immunity. The term βS(t)I(t) N(t) is referred to as the disease incidence rate and defines the total number of new infections per unit of time at time t. σ - and γ - are, respectively, the average duration of the latent and infectious periods. All the above parameters are assumed to be positive so as to represent a real situation. This model is also referred to as the SEIRS model in the literature (see, for instance, []) due to the presence of the terms ±ωR(t) in () and (). The total population dynamics can be calculated by summing up all the equations ()-(), leading tȯ It can be deduced from () that the total population is constant when ν = μ, increases when ν > μ and decreases when ν < μ. This model is discretized by using Micken's NSFD method [], to get where h >  denotes the integration step while S k , E k , I k , and R k denote the value of the susceptible, exposed, infectious and immune at sampling time t k = kh which define the discrete-time state vector The dynamics of the total population in discrete-time can be obtained by summing up again all the equations ()-() to get which corresponds to the NSFD discretization of the continuous-time equation (). This fact implies that the discretization method is consistent with the total population dynamics given by the continuous-time model, property not satisfied by all discretization methods []. Equation () can be re-written as which provides the explicit value of N k+ in terms of the system's parameters and N k . For the sake of simplicity, let us rearrange ()-() in the following implicit matrix form: Note that the matrix M(x k ) depends explicitly on x k through entries A and H, since they depend on N k and I k as () shows. Equation () allows one to isolate the value of state vector x k+ , when possible, as For the SEIR model, the closed-form explicit representation given by () is intricate. Therefore, all the properties of this model will be stated starting from the implicit form () rather than from () in comparison with some previous works [, ], where an explicit (or semi-explicit) formulation of the model is obtained prior to analysis. The following section is devoted to the study of the main properties of this scheme.

Main properties of the discrete-time model
Initially, the well-posedness of () will be proved.

Lemma  M(x k ) is invertible, its inverse is given by () and all the entries of the matrix M(x k ) - are non-negative for any set of positive parameters and any
Proof Firstly, notice that, since all the parameters (including the time step h) are assumed to be positive and I k , N k ≥ , A, B, C, D, J, F, G >  and H ≥  from (). Now, the determinant of M(x k ) can be symbolically computed as Moreover, its inverse can be symbolically computed as (). Finally, it can be directly seen from () that all its entries are non-negative since all the quantities involved are non-negative from (). Thus, the lemma is proved.
Furthermore, we have the following.
Theorem  S k , E k , I k , and R k remain non-negative for all discrete time t k = kh, k ∈ N provided that S  , E  , I  , R  ≥ ,  + (μν)h > , and all the system's parameters are positive.
Proof The proof is performed by mathematical induction. Thus, let us assume that x  is non-negative. Now suppose that x k is non-negative and it will be proved that so is x k+ .
In this way, since x k is non-negative so is N k and, in particular, I k . This fact implies that M(x k ) is invertible due to Lemma . In addition, its inverse is given by (), whose entries are non-negative. Then the (k + )th state vector, x k+ , is calculated through () as Since  + (μν)h > , all the elements in the right-hand side of () are non-negative, so is x k+ and the theorem is proved.
Remark  Notice that Theorem  states that the system is always well defined, since I k ≥  and N k ≥  for all k ≥  implies that M(x k ) is always invertible due to Lemma  and makes possible the calculation of x k+ at each time step through ().
Remark  Moreover, Theorem  can also be used to establish the range of values of the time step h that guarantees the well-posedness of the discrete-time model. Thus, from  + (μν)h >  we see that h must fulfill: Since h must meet this constraint, the method is a NSFD-based one, rather a pure NSFD method.
Remark  The non-negativity and well-posedness of the discrete-time system are crucial properties for the discrete-time model []. These properties are obtained easier by using the NSFD method than by using other discretization schemes. For instance, consider the forward Euler approximation of the derivative applied to the left-hand side of ()-() [], to obtain These equations can be cast in the matrix form: If we want to apply now the same argument as in () to prove the non-negativity of the system (i.e. the non-negativity of all the elements of the right-hand side of ()) we find two main drawbacks: (i) the first one is that the vector  possesses negative components, and, secondly, (ii) all the entries of the matrix M  are non-negative provided that the following constraints hold: Thus, the simple argument used in () is not directly applicable here, while it may lead to further conditions on the sampling period h arising from ()-(): Therefore, the NSFD method reveals as a very useful one to treat the discretization problem of the SEIR model since its properties are proved easily and under mild mathematical conditions.

Existence and stability of equilibrium points
This section studies the existence and stability of equilibrium points in the discrete-time SEIR model given by ()-().

Lemma 
and there are three potential equilibrium points given by: Proof The proof of the first part is straightforward from ()-() and ()-() while the equilibrium points are calculated in []. Now, the stability of the equilibrium points is discussed. We will focus on the two nontrivial equilibrium points obtained when μ = ν: the disease-free and the endemic equilibrium points. a For this, the nonlinear system () is linearized around these equilibrium points and the first Lyapunov theorem (linearization theorem) is used. Thus, the linearized equation reads b where x * denotes the corresponding equilibrium point and v k = x kx * is the incremental variable. In this way, the results are obtained starting from the implicit equation () rather than from the explicit one () as some previous works do [, ]. The advantage of this approach is that the calculus of the Jacobian matrix ( x * is easier implicitly than explicitly. In order to perform the calculation of the Jacobian, we must firstly write M(x k ) · x k+ in terms of only x k+ since, otherwise, we would be omitting some terms in the Taylor development. Therefore Thus, the Jacobian is given by so that around the disease-free equilibrium point and the linearized system is given by Since the incremental variables satisfy v ,k + v ,k + v ,k + v ,k =  due to the constraint S k + E k + I k + R k = N we can omit one variable in our analysis, since it can be obtained from this algebraic equation. Therefore, if we omit the first one we obtain the reduced system, The linearized reduced system () will be asymptotically stable provided that all the eigenvalues of the left-hand side matrix are strictly greater than unity in absolute value. In this way, the eigenvalues are given by It can be readily seen that λ  > , while The only eigenvalue that might threaten stability is λ  , which may be re-written as Notice that the argument of the square root is positive from (). Therefore, () implies that λ  must also be positive since CD > . The stability condition then reads λ  > . Some algebra on () shows that this condition can be met, provided that which corresponds to the same stability condition found in [] for the continuoustime system. Thus, not only does the discretization method preserve positivity and total population dynamics, but also the equilibrium points and their stability. In this way, the reproduction number of the discrete-time model is given by (). When () is met, the linearized system is globally stable around the disease-free equilibrium point, v ,k , v ,k , v ,k , v ,k →  and S k → N as k → ∞, while the nonlinear system is locally stable around the equilibrium point. Thus, we have already proved the following theorem.
On the other hand, the following theorem holds when R  > .
Theorem  If μ = ν and R  >  then one of the following claims remains true:

() The solution of the discrete-time SEIR model ()-() tends to the endemic equilibrium point. () Some variables of the discrete-time model ()-() oscillate.
Proof Firstly note that the positivity of the model along with the constraint S k + E k + I k + R k = N imply that  ≤ S k , E k , I k , R k ≤ N . Therefore, all variables are globally bounded through time. Now, since the state variables cannot diverge as k → ∞, there are only two possibilities: () All the state variables converge to a constant finite value.
() Some of the state variables do not converge to a constant value, and therefore, oscillate. If all the variables converge to a constant value, it must be an equilibrium point. The trivial point is not reachable due to the constraint S k + E k + I k + R k = N while the disease-free equilibrium point is unstable when R  > . In consequence, the finite constant limit must be the endemic equilibrium point. On the other hand, if any of the variables does not converge to a finite limit it must oscillate. Thus the theorem is proved.
Theorem  represents an alternative-type theorem to establish the behavior of the solutions when R  > . The particular case corresponding to one situation or another depends on the particular values of the parameters. However, the complexity of the endemic point expressions ()-() makes it difficult to obtain analytical conditions based entirely on the parameters to distinguish in between. This result is also in line with those presented in [] for the continuous-time model, where the conditions under which the endemic equilibrium point is globally stable are presented. The important feature regarding the case R  >  is that the disease persists in the population (in a constant or oscillatory way), and this is the reason why the equilibrium point is called endemic. In this case, the application of vaccination contributes to removing the disease. Many vaccination laws make use of the values of the subpopulations to generate the vaccination effort. Nevertheless, the number of susceptible, exposed, and immune cases may be very difficult to obtain in practice since only the infectious can be registered at health centers. Therefore, a state observer is designed in the next section for the discrete-time model.

Design of the state observer
The purpose of the state observer is to obtain an estimation of the state,x k , from the unique available data, the infectious subpopulation, I k . This issue is the main feature of the presented design process. In this way, the following feasible assumptions will be used in the sequel: Assumption  The total population N k is known at all k ≥ .

Assumption  The number of infectious is available at all k ≥ .
Remark  The total population may be obtained from () if birth and mortality rates are known or by direct measurement. This assumption is particularly feasible in the case of non-lethal diseases, for which the population remains almost constant during the spreading stage.

Remark  The number of infectious may be obtained from clinical data at health centers.
Therefore, I k does not have to be estimated, since it is directly obtained from the actual data. This fact allows us to consider the design of a reduced-order observer. On the other hand, the constraint N k = S k + E k + I k + R k for all k ≥  also permits the algebraic direct calculation of one state variable in terms of the others. In this way, ifx T k = [Ŝ kÊkÎkRk ] denotes the observed state, there is no need to estimateÎ k , whileŜ k will be calculated fromŜ Thus, we just need to observe the values ofÊ k andR k . The first step is to re-organize the state vector so as to put together the variables that are to be truly estimated. To this end define the re-organized state vector so that only the first two components of z k have to be observed. With this notation, the implicit epidemic model () can be re-written as where K k = TD - M(x k )T - . Equation () will be the starting point to design the reducedorder state observer. Now, () is expanded as , z ,k = S k , and z ,k = I k . The reason to proceed with the observer design implicitly is that matrix K k can be computed symbolically in an easy way. Hence, we have The last equation of () reads K k z ,k+ + K k z ,k+ + K k z ,k+ = z ,k () and will be used to reduce the order of the observer. Since z ,k+ = S k+ = N k+ -E k+ - which finally becomes Let us analyze the meaning of (). The right-hand side of (), δ k , is known since the infectious and the total population are known by Assumptions  and . On the other hand, the left-hand side is unknown since it just involves the actual state, which is not known. This equation will be used as error signal to drive the state observer. The observer is then designed from the first equation of (): by substituting the actual state by the observed one and by adding the observation error, leading to which is the reduced-order observer to estimate the only two components that remain to be observed. L T k = [l k l k ] ∈ R  is the so-called observer gain. Note that () is an implicit observer since the observed variableẑ ,k+ appears at both sides of equation. Therefore, the observer is designed implicitly instead of explicitly. An explicit representation of the observer, more appropriate for implementation, can be obtained by re-arranging now () to Q kẑ,k+ =ẑ ,k -L k δ k -K k N k+ + (K k -K k )z ,k+ , so that whenever Q - k exists. It will be proved in the sequel that this inverse matrix exists if the observer is designed to be asymptotically stable. In this way, () and () provide together an estimation of all the state variables. The derivation of the observer in an implicit way allows the symbolic calculation of the matrix Q k above as with the parameters defined by (). The explicit calculation of this matrix will allow us to analyze the stability properties of the observer.

Stability analysis of the observer
This section is devoted to the stability analysis of the observer. For this purpose, denote the observation error as e k = z ,k -ẑ ,k . Thus, the error equation reads The stability of the observer is discussed in the following theorem, which provides sufficient conditions to guarantee the observation of the states.
holds for all k ≥ .
Proof Consider the positive-definite Lyapunov candidate sequence V k = e T k P k e k . Its time variation is given by Since V k is positive definite and V k is negative definite, the second Lyapunov theorem guarantees the convergence to zero of the observation error and the observed variables converge to the actual ones.
Remark  It will be proved below that the stability condition () is feasible so that it is always possible to design an asymptotically stable state observer.
Remark  Equation () can be expressed as a linear matrix inequality (LMI) by using Schur complements. Therefore, the observer design problem can be efficiently solved by available numerical suites.
Remark  Since e k → , we have (E k -Ê k ) →  and (R k -R k ) →  as k → ∞. As a consequence, () implies (S k -Ŝ k ) →  as k → ∞ and all immeasurable state variables are observed.
The following corollary appears from Theorem .
Corollary  If Q k satisfies the matrix equation (), then it is invertible.
Proof Equation () can be re-written as Since S k and P k+ are symmetric positive-definite matrices, v T (S k + P k+ )v >  for any real vector v =  and all k ≥ . Thus, from () we have In this way, if the gain vector L k meets Theorem 's requirements, the observer is well defined and asymptotically stable. Furthermore, we can prove now that inequality () is feasible.
Theorem  There exists at least one gain vector L T k = [l k l k ] such that () holds.
Proof It will be proved that the equality in () holds when L k =  and S k = S is symmetric and positive definite. In this case, the matrix Q k reads whose eigenvalues are given by λ  = C + βh I k N k ≥ C =  + (μ + σ )h >  and λ  = J =  + (μ + ω)h >  for all k ≥ . Since both eigenvalues are larger than unity at all time, we see that Q - k always exists and it is a convergent matrix. The proof is based on the explicit calculation of a sequence of symmetric positive-definite matrices P k satisfying (). The following notation will be used: The rationale behind the above definitions is that R multiplies the matrices by the right while L multiplies them by the left. Moreover, the following relation between the two definitions appears: With this notation, we can now explicitly construct a solution to () as Now we can check that () is a valid solution to (): • Convergence. Since Q - k is a convergent matrix for all k ≥  the series () is convergent and well defined. • Positive definiteness. Equation () can be re-written as by using (). Thus, As a consequence, there exists at least one solution to () and the observer design problem is feasible.
Remark  A further interpretation of Theorem  is that the observer is asymptotically stable in the open-loop case since L =  implies the vanishing of the error driving term in ().

Numerical examples
This section contains some simulation examples regarding the results introduced in the previous sections. The parameters of the system are assumed to be known and given by Model's parameters are taken from an outbreak of influenza in a British boarding school in the late s, [], where the μ and ν parameters have been modified in such a way that a short-term simulation is enough to show the dynamic behavior of the system. In this way, in a period of some days the effect of natural growth can be noticed. The initial values of the populations are S  = , E  = , I  = , R  =  for a total initial population of N  = ,. Initially, Figures , , , and  depict the dynamics of the continuous and discrete-time systems for time steps of h =  hours and h =  hours. These time steps are considered since they imply that data are gathered twice or once per day, respectively, which fits in a typical sampling procedure. It can be seen that the transient response depends slightly on the value of the time step, while maintaining the shape of the response and tending to be the same value in the steadystate. Moreover, since R  >  in this example, the disease is persistent and there are always a number of infectious and exposed individuals among the population. Now, the operation of the reduced state observer will be illustrated through some numerical examples. A time step of h =  hours has been used. The observer gain has been selected as L T = l[ ] for the sake of simplicity so that there is only one free parameter.     The initial values of the populations are S  = , E  = , I  = , R  =  for a total initial population of N  = . Since the observed populations are initially unknown we may choose the values ofŜ  ,Ê  , andR  as desired while satisfyingŜ  +Ê  +R  = N  -I  =  - = . However, as a systematic way to initialize the observer we suggest to use the initial valuesŜ  = N  -I  =  andÊ  =R  = . This is always a valid selection. It can be noticed in Figures - that the actual states are observed and the estimation errors vanish asymptotically. Now, we may analyze the effect of varying the observer gain l into the dynamics of the observation. Thus, l will be changed from l = -. to l = . Figure  depicts the behavior of the observer for the susceptible when the gain is varied   show the observation error for the remaining variables. As can be deduced, the larger the observer gain is the faster the estimation of the state variables is. Hence, the most advisable tuning for the observer gain is the largest constant l satisfying the stability constraint ().

Conclusions
In this paper a discrete-time SEIR model has been obtained from a continuous-time one by using Micken's discretization method. Afterwards, the well-posedness of the model along with the existence and stability of the equilibrium points have been studied. Moreover, an implicit method has been used to analyze these properties instead of the explicit one used in previous works. It has been shown that this method preserves all the structural prop-   erties exhibited by the continuous-time model such as the expression of the reproduction number and the stability properties of the equilibrium points. In addition, the design of a state observer has been introduced for a discrete-time SEIR epidemic model. The observer was designed in an implicit way since the derivation of its stability properties as well as the design process itself are easier than when proceeding with an explicit model of the system. Moreover, some sufficient conditions to ensure the asymptotic stability of the observer have been provided in terms of a matrix inequality that can be cast in the form of a LMI. The feasibility of these conditions has been proved while some simulation examples showed the operation and usefulness of the observer. The observer is shown to be capable of estimating the actual state variables, while all the observation errors are converging to zero.