Lyapunov stability of an SIRS epidemic model with varying population

TL;DR

This paper analyzes the stability of an SIRS epidemic model with density-dependent mortality, demonstrating global stability of disease-free and endemic states using Lyapunov functions.

Contribution

It introduces a Lyapunov function approach to establish global stability for both disease-free and endemic equilibria in a generalized SIRS model.

Findings

Global stability of disease-free equilibrium proven.

Global stability of endemic state demonstrated.

Lyapunov function constructed for broad conditions.

Abstract

In this paper we consider an SIRS epidemic model under a general assumption of density-dependent mortality. We prove the global stability of the disease-free equilibrium and propose a Lyapunov function that allows to demonstrate the global stability of the (unique) endemic state under broad conditions.

Equations

$$\left\{\begin{array}[c]{l}X^{\prime}=b(N-pY)-\mu X-\beta X\dfrac{Y}{N}+\alpha Z\\ Y^{\prime}=bpY+\beta X\dfrac{Y}{N}-\left(\mu+\nu+\delta\right)Y\\ Z^{\prime}=\nu Y-\left(\alpha+\mu\right)Z\end{array}\right.$$

$$\beta \left(b-\mu \right)\left(\alpha +\mu +\nu \right)=\delta \left(\alpha +\mu \right)\left(\left(pb+\beta \right)-\left(\mu +\nu +\delta \right)\right)$$

$$N^{\prime }=\left(b-\mu \left(N\right)\right)N-\delta Y.$$

$$\mathcal{D}=\left\{\left(X,Y,Z,N\right):X\ge 0,Y\ge 0,Z\ge 0,X+Y+Z=N\le N^{\ast }\right\}.$$

$$\left\{\begin{array}[c]{l}S^{\prime}=b(1-S-pI)-\left(\beta-\delta\right)SI+\alpha R\\ I^{\prime}=\left(\beta S-\left(\left(1-p\right)b+\nu+\delta\right)+\delta I\right)I\\ R^{\prime}=\nu I-\left(b+\alpha-\delta I\right)R\end{array}\right.$$

$$N^{\prime }=N\left(b-\mu \left(N\right)\right)-\delta I.$$

$${\mathcal{D}}^{\text{fra}}=\left\{\left(S,I,R\right):S\ge 0,I\ge 0,R\ge 0,S+I+R=1\right\}$$

$$R_0=\frac{\beta }{\gamma },\gamma =\left(1-p\right)b+\nu +\delta$$

$$\gamma \left(1-R_0\right)S+b\left(1-S-pI\right)+\left(\alpha +\beta S\right)R=0$$

$$\left\{\begin{array}[c]{l}I^{\prime}=\left(\beta\left(1-I-R\right)-\gamma+\delta I\right)I\\ R^{\prime}=\nu I-\delta\left(\rho-I\right)R\end{array}\right.$$

$${\overline{\mathcal{D}}}_0^{\text{fra}}=\left\{\left(I,R\right):I\ge 0,R\ge 0,I+R\le 1\right\}\backslash \left\{\left(0,0\right)\right\}$$

$$P\left(I\right)=-\delta \left(\beta -\delta \right)I^2+\left(\delta \rho \left(\beta -\delta \right)+\left(\beta -\gamma \right)\delta +\beta \nu \right)I-\delta \rho \left(\beta -\gamma \right)=0.$$

$$\underset{I\to {\rho }^{-}}{\lim }P\left(I\right)=\beta \nu \rho >0$$

$$I^{\prime }=\left(\gamma \left(R_0-1\right)-\beta R-\left(\beta -\delta \right)I\right)I.$$

$$I^{\prime }=\left(-\beta R-\left(\delta -\beta \right)(1-I)-\left(\left(1-p\right)b+\nu \right)\right)I\le 0.$$

$$\left\{\begin{array}[c]{l}\beta\left(1-I_{e}-R_{e}\right)-\gamma+\delta I_{e}=0\\ \nu I_{e}-\delta\left(\rho-I_{e}\right)R_{e}=0\end{array}\right.$$

$$\left\{\begin{array}[c]{l}I^{\prime}=-\left(\left(\beta-\delta\right)\left(I-I_{e}\right)-\beta\left(R-R_{e}\right)\right)I\\ R^{\prime}=\left(\nu+\delta R_{e}\right)\left(I-I_{e}\right)-\delta\left(\rho-I\right)\left(R-R_{e}\right)\end{array}\right.$$

$$L_1\left(t\right)=I-I_e-I_e\ln I;L_2\left(t\right)=\frac{\beta }{2\left(\nu +\delta R_e\right)}{\left(R-R_e\right)}^2.$$

$$\frac{d{L}_1\left(t\right)}{dt}$$

$$\frac{d{L}_2\left(t\right)}{dt}$$

$$L_{\text{EE}}^{\prime }\left(t\right)=-\left(\beta -\delta \right){\left(I-I_e\right)}^2-\delta \left(\rho -I\right){\left(R-R_e\right)}^2.$$

$$0<I_u=\frac{\beta -\gamma }{\beta -\delta }=\frac{\beta -\left(\left(1-p\right)b+\nu +\delta \right)}{\beta -\delta }=1-\frac{\left(1-p\right)b+\nu }{\beta -\delta }<1.$$

$$\beta -\gamma -\rho \left(\beta -\delta \right)\le 0.$$

$$I^{\prime }\le \beta -\gamma -\left(\beta -\delta \right)I<\beta -\gamma -\rho \left(\beta -\delta \right)$$

$$\beta -\gamma -\rho \left(\beta -\delta \right)>0.$$

$$\beta >\gamma +\frac{\rho \left(\gamma -\delta \right)}{1-\rho }.$$

$$N^{\prime }=\left(b-\mu \left(N\right)-\delta I\right)N.$$

$$N^{\prime }\to (b-\delta I_e-\mu (N))N$$

$$b\le \delta I_e+\mu (0).$$

Full text

Key words: Lyapunov function, SIRS models, variable population, global asymptotic stability, endemic equilibrium.

1 Introduction

Since the seminal work initiated a century ago by Kermack and McKendrik (see [9], [5]), Mathematical Epidemiology has undertaken an extraordinary development to the point that mathematical models nowadays represent a key support of public policies aimed to control infectious diseases. All this is based on refinements of a few basic deterministic models with a simple structure, where the population is divided into states or “compartments” representing the status with respect to the infection such as e.g., susceptible, infectious, and immune, as in classical SIR and SIRS models (see [5]), which individuals can visit according to simple transition rules. These simple structures can be generalized allowing to (i) model any type of infectious diseases, ranging from vaccine preventable to vertically or sexually transmitted and to vector-borne [1], to (ii) include any type of intermediate determinants of epidemiological outputs such as e.g., the role of individual’s age in transmission, geographic structures, population dynamics, and control variables [1], such as vaccination or treatments, up to the psychological dimension of human behavior [18], and to (iii) include different kinds of nonlinearities in social contacts and transmission processes [5].

The basic models of infection spread are deterministic, continuous time, and expressed by a system of nonlinear ordinary differential equation (ODEs). The main dynamic features of these models are the existence of a disease-free equilibrium and, provided an appropriate parameter representing the reproduction of infection is above a threshold, of an endemic equilibrium where the infection persists. A key issue regards the (local and/or global) stability of these equilibria.

The involved mathematics follow two alternative approaches. The first one is based on Poincaré-Bendixon theory for planar systems and its recent multidimensional extension (see [16] and [3] for an application). The second uses the classical Lyapunov direct method, still widely applied (see [10], [11], [12], [14], [13], [15], [17], [20] and [22]).

In [20], the following SIRS model was proposed

[TABLE]

where: $b>0$ and $\mu>0$ are the birth and death rates, $\beta>0$ is the transmission rate, $\nu>0$ the rate of recovery from infection, $\delta>0$ the disease-specific mortality rate, $p\in\left(0,1\right)$ the fraction of vertically infected newborn and $\alpha>0$ the rate of return to susceptibility by loss of immunity. The stability analysis of the endemic equilibrium was performed by Lyapunov direct method under the very special assumption of constant population size $N$. However, it is easy to verify that $N=X+Y+Z$ is constant under the highly special condition on model parameters

[TABLE]

which makes this case totally uninteresting as missing any non trivial dynamics.

The aim of this paper is therefore to re-analyze the stability of steady states of model (1) in non trivial conditions by appropriately reformulating the dynamics of the population under a general assumption of density-dependent mortality. We demonstrate the global stability of the disease-free equilibrium and propose a Lyapunov function that allows to demonstrate the stability of the (unique) endemic state under broad conditions.

2 Equilibria and their stability

Following [8], [7], we assume that $\mu=\mu\left(N\right)$ is strictly increasing, with $b>\mu\left(0\right)$ and $b<\mu\left(+\infty\right)$. The resulting population dynamics is

[TABLE]

Note that, since the equation $y^{\prime}=\left(b-\mu\left(y\right)\right)y$ admits a unique globally asymptotically stable (GAS) equilibrium $y^{\ast}$, our variant to model (1) can be studied on the positively invariant set

[TABLE]

Since $N$ is not constant, it is convenient to pass to fractions $S=X/N$, $I=X/N$ and $R=Z/N$, obtaining the same epidemiological structure of system (1):

[TABLE]

completed by the following equation for the population:

[TABLE]

It is simple to verify that the region

[TABLE]

is positive invariant and attractive. Moreover, model (3) always admits the disease-free equilibrium $E_{0}=\left(1,0,0\right)\in\mathcal{D}^{\text{fra}}$. By setting

[TABLE]

it results that $E_{0}$ is the unique equilibrium when $R_{0}\leq 1$, where $R_{0}$ represents the appropriate reproduction number for system (3). Indeed, following [4], by conditions $S^{\prime}=0$, $I^{\prime}=0$ and $S=1-I-R$, for $S>0$ and $I>0$, one easily obtain the equality

[TABLE]

that shows that no solutions can exist when $R_{0}\leq 1$.

Next, we prove that, if $R_{0}>1$, a unique endemic equilibrium $E_{1}=\left(S_{e},I_{e},R_{e}\right)\in\mathcal{D}^{\text{fra}}$ for (3) exists. Since $S+I+R=1$, we can consider the reduced system

[TABLE]

where $\rho=\left(b+\alpha\right)/\delta$, on the invariant set

[TABLE]

and show that (5) has a unique positive solution $\widehat{E}_{1}=\left(I_{e},R_{e}\right)\in\overline{\mathcal{D}}_{0}^{\text{fra}}$ when $R_{0}>1$ (equivalent to $\beta>\gamma$). Note preliminarily that condition $R_{0}>1$ implies $\beta>\delta$ and it can occur $R^{\prime}=0$ only if $I<\rho$. By solving $I^{\prime}=R^{\prime}=0$, we obtain the quadratic equation in $I$

[TABLE]

Since ($P\left(0\right)<0$) and

[TABLE]

there is only one solution of $P\left(I\right)=0$ smaller than $\rho$. Therefore, if $\rho\leq 1$ the proof immediately follows, while if $\rho>1$ note that $P\left(1\right)=\delta\left(\rho-1\right)\left(\gamma-\delta\right)+\beta\nu>0$, and again the claim follows.

We first obtain the global stability of the disease-free equilibrium by adopting the Lyapunov function $I(t)$.

Theorem 1

The disease-free equilibrium $E_{0}\in\mathcal{D}^{\text{fra}}$ is GAS in $\mathcal{D}^{\text{fra}}$ if and only if $R_{0}\leq 1$, and unstable for $R_{0}>1$.

Proof. By linearization it is easy to see that $E_{0}$ is locally asymptotically stable when $R_{0}<1$, and unstable when $R_{0}>1$. We assume in the following $R_{0}\leq 1$ and we show that $L_{\text{DFE}}\left(t\right)=I\left(t\right)$ is a Lyapunov function. In fact, as $S=1-I-R$, we can write

[TABLE]

Then, by $R_{0}\leq 1$ and $\beta>\delta$, we immediately obtain $I^{\prime}\leq 0$. If $\delta\geq\beta$, for $R_{0}<1$, as $\gamma\left(R_{0}-1\right)=\beta-\gamma$, it follows

[TABLE]

Since the DFE is the only positively invariant subset of $\left\{\left(S,I,R\right):I^{\prime}=0\right\}$, by LaSalle Invariance Principle, we conclude that $E_{0}$ is GAS for $R_{0}\leq 1$.

We show now that, under suitable assumptions, there exists a Lyapunov function for the endemic equilibrium $\widehat{E}_{1}$ of system (5).

Theorem 2

If $R_{0}>1$ and $\beta\leq\gamma+\dfrac{\rho\left(\gamma-\delta\right)}{1-\rho}$, the unique endemic equilibrium $\widehat{E}_{1}$ of system (5) is GAS on $\overline{\mathcal{D}}_{0}^{\text{fra}}$.

Proof. Since, by definition, at $\widehat{E}_{1}$ it holds:

[TABLE]

and $IR-I_{e}R_{e}=I\left(R-R_{e}\right)+R_{e}\left(I-I_{e}\right)$, we can rewrite system (5) as:

[TABLE]

Let us now consider the positive functions on $\overline{\mathcal{D}}_{0}^{\text{fra}}\backslash\left\{\left(I_{e},R_{e}\right)\right\}$

[TABLE]

Along the solutions of (5) we have

[TABLE]

Therefore, function $L_{\text{EE}}=L_{1}+L_{2}$ is positive for $\left(I,R\right)\neq\left(I_{e},R_{e}\right)$ and

[TABLE]

Note that if $\rho\geq 1$, then $L_{\text{EE}}$ is a Lyapunov function and the endemic equilibrium is GAS. In passing, we also note that (see the second equation in (6) ) it holds $I_{e}<\rho$, which means that $L_{\text{EE}}$ acts as a local Lyapunov function for system (5).

Consider now case $\rho<1$. Observe that the set $\left\{\left(I,R\right)\in\overline{\mathcal{D}}_{0}^{\text{fra}}:I^{\prime}=0\right\}$ is a straight line that intersects the line $I=0$ at a point $I_{u}$ fulfilling

[TABLE]

If $I_{u}\leq\rho$, the set $\Omega=\left\{\left(I,R\right)\in\overline{\mathcal{D}}_{0}^{\text{fra}}:I<\rho\right\}$ is positively invariant. In fact, it is easy to verify that $I_{u}\leq\rho$ if and only if

[TABLE]

This implies that on $\left\{\left(I,R\right)\in\overline{\mathcal{D}}_{0}^{\text{fra}}:I=\rho\right\}$ we have $I^{\prime}\leq 0$. Furthermore, $\Omega$ is attractive as

[TABLE]

for each $I\in\left[\rho,I_{u}\right]$. Since $L_{\text{EE}}$ is a Lyapunov function on $\Omega$, the endemic equilibrium $\widehat{E}_{1}$ is GAS when $I_{u}\leq\rho$.

The result just obtained straightforwardly extends to the endemic equilibrium $E_{1}$ of system (3).

Remark 3

It is possible to verify that the case $I_{u}>\rho$ in the previous proof is far from trivial. In fact, $I_{u}>\rho$ if and only if

[TABLE]

Indeed, with simple manipulations, we obtain the equivalent condition

[TABLE]

showing that this particular case deals with situations where the interplay between demographic and epidemiological parameters favour a relatively high infection transmission.

Remark 4

The present characterization of the stability of equilibria allows clear conclusions about the effects of endemicity on the dynamics of the population

[TABLE]

Indeed, if the endemic state is GAS, then

[TABLE]

which implies that the disease will bring the extinction of the population under the condition:

[TABLE]

Conversely, if the previous condition is not met, then the disease and the population will reach the equilibrium state $N_{e}=\mu^{-1}(b-\delta I_{e})$ where the persistent presence of the disease will regulate the population size [4], [8], [19] [7].

**Acknowledgments. **The financial support of Università del Piemonte Orientale is acknowledged by the authors.