Stability and estimation problems related to a stage-structured epidemic model

In this work, we consider a class of stage-structured Susceptible-Infectious (SI) epidemic models which includes, as special cases, a number of models already studied in the literature. This class allows for n diﬀerent stages of infectious individuals, with all of them being able to infect susceptible individuals, and also allowing for diﬀerent death rates for each stage—this helps to model disease induced mortality at all stages. Models in this class can be considered as a simpliﬁed modelling approach to chronic diseases with progressive severity, as is the case with AIDS for instance. In contradistinction to most studies in the literature, we consider not only the questions of local and global stability, but also the observability problem. For models in this class, we are able to construct two diﬀerent state-estimators: the ﬁrst one being the classical high-gain observer, and the second one being the extended Kalman ﬁlter. Numerical simulations indicate that both estimators converge exponentially fast, but the former can have large overshooting, which is not present in the latter. The Kalman observer turns out to be more robust to noise in measurable data.


Introduction
The dynamics of chronic diseases that are directly transmitted appears to be very simple at a first glance: in the presence of infectious individuals, susceptible ones get infected according to some assumed force of 1 infection -in some cases there might be an exposed or latent stage. Thus, SI and SEI models (or their modified versions with more general incidence terms) seem to be very simple and appropriate descriptions of such dynamics. However, the development of such diseases is seldom as simple as described above. In many cases, an infected individual goes through different phases of disease severity, with possibly diverse forces of infection and disease-induced mortality rates -e.g. AIDS.
In this work, we analyse a class of stage-structured SI models with n infectious stages: Let x = (S, I 1 , . . . , I n ) T ,and consider the following dynamics: . . γ n−2 I n−2 − (µ n−1 + γ n−1 )I n−1 γ n−1 I n−1 − µ n I n Most studies in the literature of mathematical epidemiology modelling focus on the understanding of both local and global stability features. The rationale behind this choice is that while local stability of the disease free equilibrium (DFE) characterises the invasibility dynamics of the disease, global stability describes the long term behaviour of the disease dynamics. In many models, these features depend only on the basic reproduction number of the model -the so-called R 0 -with the DFE being both locally and globally asymptotically stable when R 0 ≤ 1, and an endemic equilibrium (usually unique) satisfying these properties when R 0 > 1. This led Shuai and van den Driessche to introduce the concept of sharp threshold property (STP), cf. [1]. On the other hand, except perhaps for the simplest models, even when the STP is satisfied, the time-scale between departure of the population state from the DFE and convergence to the EE is not necessarily short and the transient dynamics can be of significant interest. Therefore, while combined local-global analysis is mathematically correct, it can lead to an oversimplified view of the true dynamics. This change of focus brings in additional issues: under the dichotomy of models with the STP, one expects to find the population state either near the DFE or near the EE. When considering the transient dynamics, this is no longer true. Moreover, when performing a long-term analysis the only information about the initial condition is that it has infectious individuals -for the transient dynamics, the knowledge of initial condition is important. Nevertheless, it is usually not possible to measure the whole state of the system at any time -in particular, the initial condition is usually unknown. Indeed, often one has access only to a partial information of the state. These new issues lead quite naturally to the so-called observation problem: given some measurable information about a population, we want to estimate its complete state. More precisely, let x(t) and y(t) be the population state and the measurable information about this at time t -we will assume that y(t) = h(x(t)), where h is a smooth function. Let X(s) be a Lipschitz vector field that generates a global solution and consider the following problem: Give an initial condition, x(t 0 ) say, the quantities x(t) and y(t) are uniquely determined. On the other hand, assume that only y(t) is known for t ∈ [t 1 , t 2 ]. We are interested in identifying what kind of properties the vector field X, and the function h should satisfy in order that is possible to obtain an asymptotic estimatê x(t) of the state x(t). In addition, we also want to investigate howx(t) can be computed using the vector field X and the measurable output y(t) -recall that we do not know the initial condition x(t 0 ), and hence we cannot integrate the differential equationẋ = X(x). We assume that the number of infected individuals in the last stage can be measured. Therefore, we have y = h(x) = C.x = (0, . . . , 0, 1).x = I n .
Our goal is the following: knowing I n (t) for all t ≥ 0, how to derive estimatesŜ(t),Î i (t) satisfying lim t→∞Ŝ (t)− S(t) = 0 and lim t→∞Îi (t) − I i (t) = 0 and the convergence should be as fast as possible.

Outline
The remaining of the paper is organised as follows: Section 2 is devoted to a quick review of the asymptotic behaviour of the model: we give the expression of the basic reproduction number R 0 and we study the stability properties of the equilibria. In Section 3, we will show that the system with the output y = I n is observable and we give some observers (state-estimators) that allow to give dynamical estimates of the unmeasured state variables S(t), I 1 (t), . . . , I n−1 (t). More precisely, we construct two different observer: the classical high-gain and the extended Kalman filter. Next, in section 4, we present numerical experiments that indicate that convergence of the estimatesx(t) delivered by these observers towards the state x(t) is exponentially fast. However, the high-gain observer can display significant overshooting in the initial stages, while the Kalman filter observer does not present this behaviour. The latter also is robust when measurements are noisy, whereas the former might fail to converge in this case.
We finish with a discussion on the results in Section 5.
2 Revisiting a class of stage-structured epidemic models We consider the following epidemic model with multiple infectious stages: . . .
where: Λ is the recruitment, and β i is the per capita contact in the compartment I i . The function ϕ is assumed to be continuous, positive and increasing that models the exposure of susceptible individuals to contacts with infectious ones. In particular, it allows for the specification of a more general incidence rate than the standard bilinear one -for instance, one can use ϕ(S) = S p or ϕ(S) = S 1 + aS (with a > 0) to take into account saturation effects. µ 0 is the natural death rate of the susceptible individuals and µ i is the death rate of the infected individuals in stage i, in general µ i = µ 0 + d i , d i being the additional disease induced mortality rate. γ i is the transition rate from stage i to stage i + 1. The stability of this model with ϕ(S) = S has already been studied in ( [2,3]).

The equilibria and the basic reproduction number
It is easy to show that the set Ω = (S, I 1 , . . . , I n ) ∈ IR n+1 is a positively invariant set for the system.
The basic reproduction number represents the number of secondary cases produced by one infective host in entirely susceptible population. Following [4], we define by F i (S, I) the rate of appearance of new infections in compartment i, and by V i (S, I) the rate of transfer individuals in and out the compartment i by all others means. The matrix F and V are given by: Note that our V is the opposite of the same used in [4]. The Jacobian matrices at the disease free equilibrium are: The basic reproduction number is then given by R 0 = ρ(−F V −1 ), and it is equal to: Proposition 2.1. If R 0 > 1, then system (2) has a unique endemic equilibrium.
An endemic equilibrium (S * , I * 1 , . . . , I * n ) satisfy : Then for i = 2 to n − 1 and Replacing the expressions of I * i in the second equation of (2) gives:  Thus, , and since ϕ is a continuous and increasing function, it follows that there exists a unique S * ∈ (0, Λ/µ 0 ) satisfying equation (5), more precisely, Furthermore, by replacing S * , ϕ(S * ) and I * i for i = 2 to n by their expressions, we get the expression of I * 1 . After this, the expressions of I * i , i ∈ {2, . . . , n} are also uniquely determined from (3) and (4).

Global stability of the disease free equilibrium
Theorem 2.1. If R 0 ≤ 1, the disease free equilibrium is globally asymptotically stable.
Proof. Let us consider: On the other hand: and so on, until the last component of b, This shows that V is positive in Ω. Its derivative is: The terms containing I i , for i ∈ {1, . . . , n − 1}, and β n I n ϕ(S) cancel (thanks to the choice of b), and we obtainV

5
This expression is equal to: Since ϕ is an increasing positive function and R 0 ≤ 1, we haveV ≤ 0. By LaSalle Invariance Principle [5], the disease free equilibrium is globally asymptotically stable.

Global stability of the endemic equilibrium
Theorem 2.2. If R 0 > 1 then the unique endemic equilibrium EE is globally asymptotically stable.
Proof. We consider the following candidate Lyapunov function We have The derivative of V along the solutions is: On one hand, we have On the other hand, The monotonicity of ϕ together with the inequality between the arithmetic and geometric means show thatV is negative definite with respect to EE. It follows that the endemic equilibrium EE is globally asymptotically stable on the positive orthant IR n+1 + minus the stable manifold of the DFE. 7 3 Observation problem

Preliminary results
The state of System (2) will be denoted x(t) = (S(t), I 1 (t), . . . , I n (t)) T . We suppose that only the number of infected in the last stage is available for measurement, i.e, the measurable output of the system is y(t) = I n (t). Therefore we are interested in the observability properties of the following system where the expression of f is given in (2) and h(x) = x.(0, . . . , 1).
To study the observability of (6), we introduce (as usual) the following map using the successive time derivatives of the output function along the trajectories of the system: The map ψ is a diffeomorphism from from • Ω (the interior of Ω) onto ψ( • Ω), and hence system (6) is observable.
Proof. For k ∈ {1, . . . , n − 1}, the k th derivative of y is given by the following expression: Denote by c k In−p the coefficient corresponding to the variable I n−p . The expression for y (k) can be rewritten as: Moreover, the expression of y (n) can be easily deduced from the expression of y (n−1) : The Jacobian matrix of ψ at x = (S, I 1 , I 2 , . . . , I n ) is: In−2 An inductive computation shows that every c i In−i , for i = 1 to n − 1, is uniquely determined, constant and not equal to zero. The last coefficient c n S , corresponding to the variable S in x, is given by: The absolute value of the determinant of the Jacobian matrix is det for all i ∈ {1, . . . , n − 1}, it follows that det ∂ψ ∂x = 0 if and only if (I 1 , . . . , I n ) = (0, . . . , 0). Therefore, the rank of ∂ψ ∂x is then n + 1 on • Ω the interior of Ω. Now, consider x and x such that ψ(x) = ψ(x ).
Then the equality of the first component of the vectors gives I n = I n , and therefore the equality of the second component of ψ(x) and ψ(x ) gives that I n−1 = I n−1 . By repeating this reasoning until the n th − 1 component, we get I k = I k , for k ∈ {1, . . . , n}. Using these equalities in the last equation and the fact that ϕ is strictly increasing and continuous, we obtain that S = S . Then x = x , which prove that ψ is injective.
The injectivity of ψ and the full rank of ∂ψ ∂x show that ψ is a diffeomorphism from • Ω onto ψ( • Ω).
We give also the expression of Σ θ : In the original coordinates, the expression of the observer is: This can be written: , , For System (9), it is possible to construct another type of observer called extended Kalman filter (see, for instance, [8]): where r, θ are positive real numbers, Q θ = θ ∆ −1 θ Q∆ −1 θ , Q is a given symmetric positive definite n × n matrix, and ∆ θ = diag 1, 1 θ , . . . 1 θ n−1 . A (ẑ) is the Jacobian matrix of F (z) evaluated at z =ẑ, i.e, Once again, for θ ≥ 1 and large enough, the system (12) is an exponential observer for the system (9). The difference with the Luenberger extended observer (10) is that the gain 1 r Σ(t)C T used in the correction term is not constant but it is dynamically computed as a solution of a Riccati matrix differential equation. In the original coordinates, the observer is (with Σ(t) given by (12)) , , This observer works very efficiently and, unlike observer (11), there is no initial overshoot of the estimate.
Remark 3.1. It must be emphasized that use of the globally Lipschitz prolongationsΨ andψ is not only a mathematical delicacy, but it is also essential for applications. Absence of such global Lipschitz prolongations can prevent convergence of the observer, as it is illustrated in [9]. However, it is not in general an easy task to construct explicitly such globally Lipschitz prolongations.

Simulations
We simulated the various observers in the example case of a two-stage infectious class -the model given by System (8). The parameters used in these simulations were as follows: β 1 = 0.01; β 2 = 0.15; γ = 0.02; Λ = 40; µ 0 = µ 1 = 0.01; µ 2 = 0.025. The evolution of S(t) and the corresponding estimates using the Kalman observer -Equation (13) -can be seen in Figure 1. These estimates are robust as it can be seen from the simulations in Figure 2, where the measurable output y(t) = I 2 (t) was corrupted with noise. In Figure 4 we plotted the evolution of I 1 (t) together with the Kalman observer estimates. Figure 4 displays the corresponding noisy simulations.  The evolution of S(t) and its estimateŜ(t) delivered by the observer (13) when the measurable output is corrupted by noise.
13 Figure 3: The evolution of I 1 (t) and its estimateÎ 1 (t) delivered by the observer (13). Figure 4: The evolution of I 1 (t) and its estimateÎ 1 (t) delivered by the observer (13) when the measurable output is corrupted by noise.
The estimates of S(t) and I 1 (t) provided by the High Gain observer are displayed in Figures 5 and 6. In this case there is an initial overshoot of the estimate. When the measurable output, y(t), was corrupted by noise, we were not always able to verify the convergence of this observer, and hence we do not display any noisy results for the high gain observer.

S(t)
S  (t)  Figure 5: The evolution of S(t) and its estimateŜ(t) delivered by the high-gain observer (11). One can notice a very important initial overshoot (for t < 1) of the estimate. I 1 (t) Figure 6: The evolution of I 1 (t) and its estimateÎ 1 (t) delivered by the high-gain observer (11). Once again there is an initial overshoot of the estimate.

Discussion
In this work we considered and analysed a class of stage-structured SI models, which included a number of previously studied models in the literature -cf. [6,8].
We begin by identifying an expression for the basic reproductive number (R 0 ) for these models, and we also show that they satisfy the Sharp Threshold Property [1]. We then proceed to present some preliminary results on the observability of these models, and after that we derive two state-estimators. These estimators are then put to work in Section 4 where various numerical experiments are presented. Both estimators show exponential convergence, but the extended Luenberger estimator usually displays some initial overshooting in contradistinction to the extended Kalman filter. Indeed, the behaviour of the latter is very accurate, even when the measurable output y(t) = I 2 (t) is corrupted by noise. On the other hand, the high-gain Luenberger observer is very sensitive to noisy measurements, and it might fail to converge in this case.