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

1.   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 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, \ldots, I_n)^T$, and consider the following dynamics:

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\leq1$, 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\in [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 estimate $\hat x(t)$ of the state $x(t)$. In addition, we also want to investigate how $\hat x(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 $\dot x = 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, \dots, 0, 1). x = I_n$.

Our goal is the following: knowing $I_n(t)$ for all $t\geq 0$, how to derive estimates $\hat S(t), \hat I_i(t)$ satisfying $\lim_{t\rightarrow\infty} \hat S(t) - S(t) = 0$ and $\lim_{t\rightarrow\infty} \hat I_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 $\mathcal{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), \ldots, 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 estimates $\hat x(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:

$\Lambda$ is the recruitment, and $\beta_i$ is the per capita contact in the compartment $I_i$. The function $\varphi$ 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 $\varphi(S) = S^p$ or $\varphi(S) = \dfrac{S}{1+a S}$ (with $a > 0$) to take into account saturation effects. $\mu_0$ is the natural death rate of the susceptible individuals and $\mu_i$ is the death rate of the infected individuals in stage $i$, in general $\mu_i = \mu_0 +d_i$, $d_i$ being the additional disease induced mortality rate. $\gamma_i$ is the transition rate from stage $i$ to stage $i+1$. The stability of this model with $\varphi (S) = S$ has already been studied in (^[2,3]).

2.1.   The equilibria and the basic reproduction number

It is easy to show that the set

is a positively invariant set for the system.

The disease free equilibrium is $DFE = (\dfrac{\Lambda}{\mu_0}, 0, \ldots, 0)$.

The basic reproduction number represents the number of secondary cases produced by one infective host in entirely susceptible population. Following ^[4], we define by $\mathcal{F}_i(S, I)$ the rate of appearance of new infections in compartment $i$, and by $\mathcal{V}_i(S, I)$ the rate of transfer individuals in and out the compartment $i$ by all others means. The matrix $\mathcal{F}$ and $\mathcal{V}$ are given by:

and

Note that our $\mathcal{V}$ is the opposite of the same used in ^[4]. The Jacobian matrices at the disease free equilibrium are:

and

The basic reproduction number is then given by $\mathcal{R}_0 = \rho(-F V^{-1})$, and it is equal to:

Proposition 2.1. If $\mathcal{R}_0 > 1$, then system (2.1) has a unique endemic equilibrium.

An endemic equilibrium $(S^*, I_1^*, \ldots, I_n^*)$ satisfy :

Then for $i = 2$ to $n-1$

and

Replacing the expressions of $I_i^*$ in the second equation of (2.1) gives:

Thus,

Therefore, if $\mathcal{R}_0 > 1$ then $\dfrac{\varphi(\Lambda/\mu_0)}{\mathcal{R}_0} < \varphi(\Lambda/\mu_0)$, and since $\varphi$ is a continuous and increasing function, it follows that there exists a unique $S^* \in (0, \; \Lambda/\mu_0)$ satisfying equation (2.4), more precisely, $S^* = \varphi^{-1}\left(\dfrac{\varphi(\Lambda/\mu_0)}{\mathcal{R}_0}\right)$.

Furthermore, by replacing $S^*$, $\varphi(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\in\{2, \ldots, n\}$ are also uniquely determined from (2.2) and (2.3).

2.2.   Global stability of the disease free equilibrium

Theorem 2.1. If $\mathcal{R}_0\leq 1$, the disease free equilibrium is globally asymptotically stable.

Proof. Let us consider:

where

and

On the other hand:

Since $\mathcal{R}_0\leq 1$, we have $\mu_{1}+\gamma_1-\beta_1\varphi(\Lambda/\mu_0) > 0$, and also

and so on, until the last component of $b$,

This shows that $V$ is positive in $\Omega$. Its derivative is:

The terms containing $I_i$, for $i\in\{1, \ldots, n-1\}$, and $\beta_n I_n \varphi(S)$ cancel (thanks to the choice of $b$), and we obtain

This expression is equal to:

Since $\varphi$ is an increasing positive function and $\mathcal{R}_0\leq 1$, we have $\dot{V}\leq 0$. By LaSalle Invariance Principle ^[5], the disease free equilibrium is globally asymptotically stable.

2.3.   Global stability of the endemic equilibrium

Theorem 2.2. If $\mathcal{R}_0 > 1$ then the unique endemic equilibrium EE is globally asymptotically stable.

Proof. We consider the following candidate Lyapunov function

where

We have

The derivative of $V$ along the solutions is:

On one hand, we have

On the other hand,

Therefore,

with

Hence,

The monotonicity of $\varphi$ together with the inequality between the arithmetic and geometric means show that $\dot V$ is negative definite with respect to EE. It follows that the endemic equilibrium EE is globally asymptotically stable on the positive orthant ${\bf R}^{n+1}_{+}$ minus the stable manifold of the DFE.

3.   Observation problem

3.1.   Preliminary results

The state of System (2.1) will be denoted $x(t) = (S(t), I_1(t), \ldots, 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.1) and $h(x) = x.(0, \ldots, 1)$.

To study the observability of (3.1), we introduce (as usual) the following map using the successive time derivatives of the output function along the trajectories of the system:

Lemma 3.1. The map $\psi$ is a diffeomorphism from from $\overset{\circ}{\Omega}$ (the interior of $\Omega$) onto $\psi(\overset{\circ}{\Omega})$, and hence system (3.1) is observable.

Proof. For $k\in\{1, \ldots, n-1\}$, the $k^{th}$ derivative of $y$ is given by the following expression:

Denote by $c_{I_{n-p}}^k$ 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 $\psi$ at $x = (S, I_1, I_2, \ldots, I_n)$ is:

An inductive computation shows that every $c_{I_{n-i}}^{i}$, for $i = 1$ to $n-1$, is uniquely determined, constant and not equal to zero. The last coefficient $c_S^n$, corresponding to the variable $S$ in $x$, is given by:

The absolute value of the determinant of the Jacobian matrix is $\left|\text{det}\left(\dfrac{\partial \psi}{\partial x}\right)\right| = c_S^n \prod_{i = 1}^{n-1}c_{I_{n-i}}^{i}$. Since $c_{I_{n-i}}^{i} \neq 0$ for all $i\in \{1, \ldots, n-1\}$, it follows that $\text{det}\left(\dfrac{\partial \psi}{\partial x}\right) = 0$ if and only if $(I_1, \ldots, I_n) = (0, \ldots, 0)$. Therefore, the rank of $\dfrac{\partial \psi}{\partial x}$ is then $n+1$ on $\overset{\circ}{\Omega}$ the interior of $\Omega$. Now, consider $x$ and $x'$ such that $\psi(x) = \psi(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 $\psi(x)$ and $\psi(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\in\{1, \ldots, n\}$. Using these equalities in the last equation and the fact that $\varphi$ is strictly increasing and continuous, we obtain that $S = S'$. Then $x = x'$, which prove that $\psi$ is injective. The injectivity of $\psi$ and the full rank of $\dfrac{\partial \psi}{\partial x}$ show that $\psi$ is a diffeomorphism from $\overset{\circ}{\Omega}$ onto $\psi(\overset{\circ}{\Omega})$.

3.2.   Some state estimators

We are now able to construct an observer for (3.1).

We begin by performing variable change $z = \psi(x)$, and we determine $\psi^{-1}$. Through the identification $z_i = \psi_i(x)$, we obtain that:

repeating this until the $(n-1)$-th component, we show that:

and the $n+1$ th component gives that:

recall that $\varphi$ is continuous and strictly increasing.

With the new coordinate $z$, system (3.1) is given by the following:

We denote the matrix in (3.2) by $A$ and $C = (1, \ldots, 0)$. It is easily shown that $\Psi(z) = y^{(n+1)}(\psi^{-1}(z))$ is $k$-Lipschitzian on $\Omega$. The observer of (3.2) is given by (^[7,6]):

where $\Sigma_\theta$ is the solution of the equation:

The resolution of the equation leads to:

The system

is an observer for (3.1).

Application: We apply the result for $n = 2$ and $\varphi(x) = x$. The corresponding system is:

The output is $y = I_2$ and

$\psi(S, I_1, I_2) = \left(\begin{array}{c} y\\ \dot{y}\\ \ddot y \end{array} \right)$ = $\left(\begin{array}{c} I_2\\ \gamma I_1-\mu_{2} I_2\\ \gamma(\beta_1 I_1+\beta_2 I_2) S-\gamma(\mu_{1}+\mu_{2}+\gamma) I_1+\mu_{2}^2 I_2 \end{array} \right)$.

Its Jacobian can be written as:

The determinant of the Jacobian is equal to $\det\left(\dfrac{\partial \psi}{\partial x}\right) = \gamma^2(\beta_1 I_1+\beta_2 I_2)$, which is positive in the positively invariant open set $\Omega_0 = \left\{ S > 0, I_1 > 0, I_2 > 0, S+I_1+I_2 < \dfrac{\Lambda}{\mu_0}\right\}$.

We give also the expression of $\Sigma_\theta$:

With the variable change, we have: $z = \psi(S, I_1, I_2)$, which leads to:

and when we replace $I_1$ and $I_2$ by their expression, we get:

The expression of $\psi^{-1}$ is:

According to the expression of $\psi^{-1}$ we deduce $\Psi(z)$:

With the variable $z$, system (3.3) is given by:

Let $\tilde\Psi$ be a globally Lipschitz map on the whole ${\bf R}^3$ whose restriction to $\Omega_0$ is $\Psi$, and $\tilde\psi$ an extension of the diffeomorphism $\psi$ to the whole ${\bf R}^3$.

The observer of (3.4) is simply:

In the original coordinates, the expression of the observer is:

This can be written:

For System (3.4), it is possible to construct another type of observer called extended Kalman filter (see, for instance, ^[8]):

where $r$, $\theta$ are positive real numbers, $Q_{\theta} = \theta\, \Delta_{\theta}^{-1} Q \Delta_{\theta}^{-1}$, $Q$ is a given symmetric positive definite $n\times n$ matrix, and $\Delta_{\theta} = \text{diag}\left(1, \dfrac{1}{\theta}, \ldots\dfrac{1}{\theta^{\, n-1}}\right)$. $A^{\star}(\hat{z})$ is the Jacobian matrix of $F(z)$ evaluated at $z = \hat{z}$, i.e, $A^{\star}(\hat{z}) = \left.\dfrac{\partial}{\partial z}(F(z)) \right|_{z = \hat z}$.

Once again, for $\theta \geq 1$ and large enough, the system (3.7) is an exponential observer for the system (3.4). The difference with the Luenberger extended observer (3.5) is that the gain $\dfrac{1}{r}\Sigma(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 $\Sigma(t)$ given by (3.7))

This observer works very efficiently and, unlike observer (3.6), there is no initial overshoot of the estimate.

Remark 3.1. It must be emphasized that use of the globally Lipschitz prolongations $\tilde\Psi$ and $\tilde\psi$ 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.

4.   Simulations

We simulated the various observers in the example case of a two-stage infectious class – the model given by System (3.3). The parameters used in these simulations were as follows: $\beta_1 = 0.01; \beta_2 = 0.15; \gamma = 0.02; \Lambda = 40; \mu_0 = \mu_1 = 0.01; \mu_2 = 0.025$.

The evolution of $S(t)$ and the corresponding estimates using the Kalman observer — Equation (3.8) — 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 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.

5.   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.

Acknowledgments

We are grateful to two anonymous referees for valuable comments and suggestions that led to an improvement of this paper.

This work was partially funded by Inria and CAPES-Brazil in the framework of the program STIC AmSud (project MOSTICAW). MOS was partially supported by CNPq under grant # 309079/2015-2 and by CAPES — Finance Code 001.

Conflict of interest

All authors declare no conflicts of interest in this paper.