Threshold dynamics of an uncertain SIRS epidemic model with a bilinear incidence

1 Introduction

As a major public problem worldwide, the epidemics of infectious diseases in history have brought great disaster to the survival of human beings, national economy and person’s livelihood. Over the past century, the population have suffered from influenza, AIDS, SARS tuberculosis and other infectious diseases that have killed hundreds of millions of people. The sudden outbreak of COVID-19 at the end of 2019 and its worldwide spreading not only affect the health of individuals, but also do great harm to economic and social development. At present, mathematical models have been widely employed in the study of infectious diseases. The most classic disease dynamics model is the Kermack-Mckendrick’s SIR compartmental model [1] (1927) which divided the total population into three compartments: susceptible, infected and recovered, and marked the number of people at time t as S _t , I _t and R _t , respectively. Kermack W. and McKendrick A. [2] (1932) proposed SIS compartmental model and puted forward the “threshold theory" to distinguish between disease prevalent or not. On this basis, Allen L. [3] (1994) proposed some discrete-time SI, SIR, and SIS epidemic models. After that, many scholars have derived various ordinary differential models, such as SEIR model, SEAIR model, SIRI model, and SIRS model. The research work of these typical deterministic epidemic models can be found in the literature [4–8].

In addition to deterministic models, based on the consideration that the spread of disease is subject to stochastic perturbations, stochastic models are established. Gray A. et al. [9] (2011) studied the stationary distribution and ergodicity of a random SIS epidemic model by using Has’ minskii’s ergodicity theory. Ji C. et al. [10] (2012) discussed the stationary distribution and ergodicity of a random multi-group SIR model. Recently, various random SIRS models have been studied. El Fatini M. et al. [11] (2019) considered a stochastic SIRS epidemic model with general incidence and studied its stationary distribution and threshold dynamics. Zhao X. et al. [12] (2020) formulated a stochastic SIRS model and given the sensitivity analysis of the basic reproduction number under different semi-Markov chains. The stochastic epidemic model has been studied deeply in all aspects and has made great achievements.

Nevertheless, at the beginning of epidemic’s outbreak, there are not enough samples to estimate the probabilistic properties of the random model, then some dominate expert’s knowledge is required to assess the belief degree that each event will happen. At this point it is not appropriate to consider the propagation process subject to random perturbations based on probability theory. For dealing with uncertainties, Zadeh L. A. [13] (1965) proposed Fuzzy Sets, where a fuzzy set is defined by its membership function. and then Granular computing have been studied explicitly or implicitly in many fields. The approximation space in rough set theory is important for dealing with uncertainties. Xu, W. et al. [14] (2022) mainly investigated a dynamic approximation update mechanism for multigranulation data from a local viewpoint. Granular computing (GrC) and two-way concept learning (TCL) are influential studies of knowledge processing and cognitive learning. Xu, W. et al. [15] (2022) designed a fuzzy-based progressive learning mechanism within this framework to provide a novel and convenient method for researching two-way learning and concept-cognitive learning under dynamic environments. However, this article deals with the issue of belief degree from the perspective of uncertainty theory, which was established by Liu [16] (2007). Liu introduced an uncertain measure based on normality, duality, subadditivity and product axioms. Liu [17] (2008) introduced a type of differential equation driven by Liu process and then employed the uncertainty theory into various practical field [18, 19]. Chen X. and Liu B. [20] (2010) proposed the existence and uniqueness theorem of solution for the uncertain differential equation. Yao K. and Chen X. [21] (2013) proposed the α-path methods to obtain the numerical solutions of uncertain differential equations. Li Z. et al. [22] (2017) established an uncertain SIS model with the bilinear incidence rate via uncertain differential equations and obtained the analytical solution and the corresponding α-paths of this model. On the basis, Li Z. et al. [23] (2018) further studied uncertain SIS epidemic model with standard incidence and analyzed the extinction and permanence of the disease for a threshold quantity. Furthermore, Li Z. et al. [24] (2018) discussed the comparison with the deterministic, stochastic, and uncertain SIS models. Jia L. et al. [25] (2020) proposed an uncertain SEIAR model and estimated the parameters via moment estimation method to predict the possible numbers of active cases for COVID-19 in China. Liu Z. [26] (2021) studied the cumulative cases of COVID-19 based on the uncertain growth model. Chen X. et al. [27] (2021) proposed an uncertain SIR model to describe the development trend of COVID-19. Given that viruses can mutate over time and human resistance to drugs may increase, the risk of a decrease or even gradual disappearance of immunity is significant. Consequently, studying uncertain SIRS epidemic models is crucial in accurately describing epidemics with temporary immunity and informing effective control strategies.

Our main findings are summarized as follows. Firstly, uncertain processes are utilized to describe the non-deterministic factors in the epidemic transmission process, and an uncertain SIRS epidemic model in non-constant population size with constant inputs for only susceptible persons is established. Secondly, the concept of basic reproduction number is introduced, and the uncertain threshold function is defined to investigate the asymptotic stability of disease extinction and persistence in an uncertain environment.

The rest of this paper is structured as follows. In Section 2, the SIRS model under uncertainty theory is established. Section 3 and Section 4 show the expression of uncertain threshold function R 0 * by use of the Yao-Chen formula to establish α-path of the proposed uncertain SIRS model and then the asymptotic behavior of equilibria states are discussed. In Section 5, some numerical experiments are given. Finally, some conclusions are drawn in Section 6.

2 Model formulation

In this section, we give an uncertain SIRS epidemic model with constant input and exponential death, that is, the total population is non-constant, only susceptible persons have constant input. For this model, we divide the population being studied into three categories labeled S, I, and R,respectively. Meanwhile the numbers of individuals in the three groups at the time t are marked by S _t , I _t and R _t , respectively, where S _t represents the number of individuals susceptible to the disease, that is, who are not (yet) infected at time t, and I _t represents the number of infected who are presumed to be infectious and able to spread the disease by contact with susceptibles, and R _t represents the number of individuals who have been infected and then recovered from the possibility of re-infection or transmission of infection. The general SIRS model studied by Muroya Y. et al. [7] (2014) is given by ordinary differential equation as follows:

{ d S = ( A - ν S - β SI + δ R ) d t d I = ( β SI - ( r + ν + μ ) I ) d t d R = ( rI - ( δ + ν ) R ) d t(1) where A is the recruitment rate, β is the disease transmission coefficient and r represents the rate of recovery from infection, ν indicates the natural mortality rate, μ represents the disease-related mortality coefficient. Moreover, δ denotes the coefficient of loss of immunity.

It can be seen from the Figure 1 that the schematic diagram of SIRS epidemic spreading.

OPEN IN VIEWER

Fig. 1

The schematic diagram of SIRS epidemic spreading.

It is based on the following assumptions:

Assume that there is entry and exit in the total population so that the total population is non-constant. Let N be the total population, then N = S + I + R. Adding the three equations of the model (1), we obtain that N meets the equation d N d t = A - ν N - μ I ;

In real life, actually, the spread of a pandemic is affected by environmental noise. Based on the deterministic SIRS model, this paper considers that the contact transmission coefficient β will be affected by uncertain factors in the environment during the infectious disease spreading. Due to the influence of human factors such as travel restrictions and isolation policy in the process of disease transmission, the prediction of the stochastic model constructed with random disturbances is not accurate enough. Therefore, Liu process is introduced to depict uncertain factors in the environment and transform the rate β into β _t with uncertain perturbations that capture variability in transmission. The model (1) is considered to replace the deterministic transmission rate with a uncertain one. This is done by letting β t d t ⟶ β d t + σ d C t , where σ represents the intensity of Liu process, C _t is said to be a Liu process, which is a normal uncertain variable with expected value 0 and variance t², namely, C t ~ N ( 0 , t ) . Moreover C _t has an uncertainty distribution Φ t ( x ) = ( 1 + exp ( - π x 3 t ) ) - 1 .

Replace the β in the original deterministic SIRS model with β _t and convert it into a system of uncertain differential equations, which is called an uncertain SIRS epidemic model. The resulting uncertain model is introduced as follows:

{ d S t = ( A - ν S t - β S t I t + δ R t ) d t - σ S t I t d C t d I t = ( β S t I t - ( r + ν + μ ) I t ) d t + σ S t I t d C t d R t = ( r I t - ( δ + ν ) R t ) d t .(2)

The basic reproduction number R 0 [28], defined as the expected number of secondary infections caused by a single infected individual in a fully susceptible population during the average duration of their illness, plays an important role in studying the existence and asymptotic stability of equilibrium points in epidemic models. Specifically, R 0 provides a measure of the average number of individuals that a patient can infect when the entire population is initially susceptible to the disease. For the general SIRS model (1), the following conclusions hold:

The basic reproduction number R 0 = ( β ) ( A / ν ) r + ν + μ .

Considering the uncertain SIRS epidemic model (2), we plan to introduce Yao-Chen formula [21], which relates uncertain differential equations and ordinary differential equations, to further discuss an uncertain threshold function R 0 * which determines there is a pandemic.

The α-path of the uncertain model (2) is defined based on each uncertain differential equation. Yao K. and Chen X. [21] suggested that a solution of an uncertain differential equation is a contour process. Thus, the solutions S _t , I _t and R _t of uncertain differential system (2) are contour processes. The α-path (0 < α <1) of the uncertain model (2) with initial values S₀, I₀ and R₀ correspond to deterministic functions S t α , I t α and R t α with respect to t respectively. The corresponding ordinary differential equations are deduced,

{ d S t 1 - α = ( A - ν S t 1 - α - β S t 1 - α I t α + δ R t α ) d t + | - σ S t 1 - α I t α | Φ - 1 ( 1 - α ) d t d I t α = ( β S t 1 - α I t α - ( r + ν + μ ) I t α ) d t + | σ S t 1 - α I t α | Φ - 1 ( α ) d t d R t α = ( r I t α - ( δ + ν ) R t α ) d t(3) where Φ^-1 (α) is the inverse uncertainty distribution of standard normal uncertain variable, i.e., Φ - 1 ( α ) = 3 π ln α 1 - α , 0 < α < 1 .

Consider system (2), we have

d N t = ( A - ν N t - μ I t ) d t ,(4) It can be seen from equation (4) that when the disease does not exist (I _t  = 0), the total population N _t eventually tends to a constant A/ν. So the solution of system (2) eventually fall into or stay in the area D = { | S t ≥ 0 , I t ≥ 0 , R t ≥ 0 S t + I t + R t ≤ A ν | ∀ t > 0 } .

Theorem 1. Let R 0 * = R 0 + ( σ Φ - 1 ( ( α ) ) ( A / ν ) r + ν + μ and σ > 0 with α ∈ (0, 1). Then for any non-negative initial values (S₀, I₀, R₀), we have

Suppose that R 0 * ≤ 1 , the system (3) has an unique disease-free equilibrium state E₀ (A/ν, 0, 0).

Proof 1. At first, we need to discuss the existence of steady state points of system (3). Set x = S t 1 - α , y = I t α , z = R t α , we can find two equilibrium points under the solutions of following equations: A - ν x - β xy + δ z + σ xy Φ - 1 ( 1 - α ) = 0(5) β xy - ( r + ν + μ ) y + σ xy Φ - 1 ( α ) = 0(6) ry - ( δ + ν ) z = 0 .(7)

It is easy to get one of solution of system (5-7), that is the disease-free equilibrium state E₀ (A/ν, 0, 0). After that, we try to solve the positive epidemic solution E₁ (x^*, y^*, z^*). From the equation (7), we get the expression of z in terms of y:

z * = r δ + ν y * .(8) For y ≠ 0, using the equation (6) dividing by y, we have

x * = r + ν + μ β + σ Φ - 1 ( α ) .(9) Summing the three equations (5-7), we have

A - ν ( x + y + z ) + μ y + σ xy ( Φ - 1 ( 1 - α ) + Φ - 1 ( α ) ) = 0 ,(10) where Φ^-1 (1 - α) + Φ^-1 (α) =0. Substituting (8) and (9) in the above equation (10), we can easily get the expression of y:

y * = ( 1 - 1 R 0 * ) ( 1 + r δ + ν + μ ν ) - 1 .(11)

If R 0 * ≤ 1 , then y^* ≤ 0 . Since y ∈ D, we deduce y^* = 0 for the equation (11). There is a contradiction with the premise y ≠ 0. Therefore the system (3) has an unique disease-free equilibrium state without any positive endemic equilibria.

If R 0 * > 1 , we have y^* > 0, which means the model (3) has positive endemic equilibrium state E₁ (S^*, I^*, R^*). Thus, in addition to E₀, the model (3) has a series of positive equilibrium state E₁ for 0 < α < 1.

In this section, we discusses the dynamics of the system (3) by analyzing the asymptotic behaviors of two steady states.

Theorem 2. Suppose that R 0 * ≤ 1 , then the unique disease-free equilibrium state E₀ (A/ν, 0, 0) of the system (3) is globally asymptotically stable.

Proof 2. If R 0 * ≤ 1 , then γ + ν + μ ≥ (β + σΦ^-1 (α)) (A/ν). Define the non-negative Lyapunov function [29, 30] V t = I t α and calculate the derivative with respect to t, then we get V t ʹ = d I t α d t = β S t α I t α − ( r + ν + μ ) I t α + σ S t α I t α Φ − 1 ( α ) = β + σ Φ − 1 ( α ) S t α − ( r + ν + μ ) I t α ≤ β + σ Φ − 1 ( α ) A d − ( r + ν + μ ) I t α = ( r + ν + μ ) R 0 * − 1 I t α n ≤ 0. Obviously, equation holds if and only if I t α = 0 . Substitute I t α = 0 into the system (3), we have lim t ⟶ ∞ S t α = A / ν and lim t ⟶ ∞ R t α = 0 . By LaSalle’s invariance principle, the disease-free equilibrium of system (3) is globally asymptotically stable. Hence, the proof is complete.

In the following part, we study the local asymptotic stability near the endemic equilibrium state E₁.

Theorem 3. Suppose that R 0 * > 1 and β + σΦ^-1 (α) >0, then the equilibrium state E₁ is locally asymptotically stable.

Proof 3. To study the local behavior near the steady state E₁, we linearize model (3) in the neighborhood of E₁. The Jacobian matrix associated with system (5-7) evaluated at (x, y, z) is given by

J = ( - ν - β y + σ Φ - 1 ( 1 - α ) y - β x + σ Φ - 1 ( 1 - α ) x δ β y + σ Φ - 1 ( α ) y β x - ( r + ν + μ ) + σ Φ - 1 ( α ) x 0 0 r - ( δ + ν ) ) .

Substitute E₁ (x^*, y^*, z^*) to the above matrix, then the jacobian matrix at E₁ is J ( E 1 ) = ( J ik ) 1 ≤ i , k ≤ 3 ( - ν - ( β + σ Φ - 1 ( α ) ) y * - ( β + σ Φ - 1 ( α ) ) x * δ ( β + σ Φ - 1 ( α ) ) y * 0 0 0 r - ( δ + ν ) ) .

Then, the characteristic polynomial associated to J (E₁) is given by λ 3 + a 1 λ 2 + a 2 λ + a 3 = 0 , where λ represents the eigenvalues of the jacobian matrix, and a 1 = - tr ( J ( E 1 ) ) = - ( J 11 + J 22 + J 33 ) , a 2 = ( J 22 J 33 - J 23 J 32 ) + ( J 11 J 33 - J 13 J 31 ) + ( J 11 J 22 - J 12 J 21 ) , a 3 = - det ( J ( E 1 ) ) . That is a 1 = ν + ( β + σ Φ - 1 ( α ) ) y * + ( δ + ν ) , a 2 = ( ν + ( β + σ Φ - 1 ( α ) ) y * ) ( δ + ν ) + ( β + σ Φ - 1 ( α ) ) 2 x * y * , a 3 = ( β + σ Φ - 1 ( α ) ) y * · ( ( δ + ν ) ( β + σ Φ - 1 ( α ) ) x * - r δ ) . For the above three equations, it is obvious that a₁, a₂ are positive under the condition β + σΦ^-1 (α) >0. Using equation (9) to simplify the third of above equations, we get a 3 = ( β + σ Φ - 1 ( α ) ) y * ( ( δ + ν ) ( r + ν + μ ) - r δ ) > 0 . Furthermore, we have a 1 a 2 - a 3 = ( ν + ( β + σ Φ - 1 ( α ) ) y * ) · [ ( ν + ( β + σ Φ - 1 ( α ) ) y * ) ( δ + ν ) + ( β + σ Φ - 1 ( α ) ) ( r + ν + μ ) y * ] + ( δ + ν ) · [ ( ν + ( β + σ Φ - 1 ( α ) ) y * ) ( δ + ν ) ] + ( β + σ Φ - 1 ( α ) ) y * r δ > 0 .

According to Routh-Hurwitz criteria [29, 30], we may conclude that since all the eigenvalues of J (E₁) have negative real part, the endemic steady state E₁ exists and has local asymptotic stability for R 0 * > 1 .

5 Numerical algorithm and some examples

In this section, we perform a computing algorithm to obtain numerical solutions of the uncertain SIRS model (2) on the basis of Euler method and give some examples to illustrate theoretical results proposed in this paper.

Algorithm: To solve the model (2), it is suggested to obtain spectra of α-paths of S _t , I _t and R _t when analytic solutions are extremely difficult to obtain.

Step 1. Fix α on (0, 1) and take time t ∈ (0, T) with final time T = 200. Set the initial values S₀, I₀ and R₀ and various parameters A, β, ν, μ, r, δ and σ.

Step 2. Solve the system (3) and obtain the α-paths S t α , I t α and R t α , respectively. It is solved by recursion formula as follows: { S t + 1 1 - α = S t 1 - α + ( A - ν S t 1 - α - β S t 1 - α I t α + δ R t α ) + σ S t 1 - α I t α Φ - 1 ( 1 - α ) h I t + 1 α = I t α + ( β S t 1 - α I t α - ( r + ν + μ ) I t α ) + σ S t 1 - α I t α Φ - 1 ( α ) h R t + 1 α = R t α + ( r I t α - ( δ + ν ) R t α ) where h is the step length.

Step 3. The α-paths S t α , I t α and R t α are obtained, which give the approximate uncertainty distribution of S _t , I _t and R _t . Then M { S t ≤ S t α } = α , M { I t ≤ I t α } = α and M { R t ≤ R t α } = α .

Example 1. This example studies the dynamic changes of the uncertain SIRS model (2) with the same parameter values as shown in Table 1, but with different initial values. It can be seen from the Figure 2 that the results of the numerical simulation are consistent with the theoretical analysis, that is endemic equilibria asymptotically stable for any initial values (S₀, I₀, R₀).

OPEN IN VIEWER

Fig. 2

Trajectories of S t α , I t α and R t α in Example 1.

In the subfigure (a), the solutions of system (3) converge to the disease-free equilibrium state E₀ (3, 0, 0) with different initial values, that is, the infection will die out. Here R 0 * ≈ 0.6777 and R 0 * ≈ 0.9938 correspond to α = 0.2, α = 0.9, respectively. In the subfigure (b), we can see that I t α does not tend to 0, that is the disease will persist in mean. Here R 0 * ≈ 2.4204 and R 0 * ≈ 3.5494 correspond to α = 0.2, α = 0.9, respectively.

Example 2. The model parameter values are shown in Table 2. For any given α, the α-paths of S _t , I _t , and R _t can be algorithmically obtained using the uncertain SIRS epidemic model(2). In Figure 3, we show the trajectories of S t α , I t α and R t α with initial value (8, 1, 0) and σ = 0.004 with α = 0.1, 0.2, ⋯ , 0.9.

OPEN IN VIEWER

Fig. 3

Trajectories of S t α , I t α and R t α in Example 2.

In the left subfigure, R 0 * < 1 for different α ∈ (0, 1). It is clear that the trajectories of α-path S t α , I t α and R t α converge to E₀. In the right one, R 0 * > 1 for each α, we see that the solutions converge to E₁. By numerical simulation, we can see when R 0 * > 1 , E₁ may also be globally asymptotically stable.

Example 3. In the general SIRS model (1) and uncertain SIRS model (2), take A = 0.1, ν = 0.02, μ = 0.005, β = 0.02, r = 0.05, δ = 0.005 and σ = 0.04 with initial value (9, 1, 0). For the ordinary differential equation (1), R 0 ≈ 1.333 > 1 . Thus, the endemic equilibrium state (3.75, 0.385, 0.769) in the model (1) is asymptotically stable. By Theorem 2, the values of R 0 * for α = 0.1, 0.2, ⋯ , 0.9 is provided in Table 3. Then the curves of S t α , I t α and R t α are illustrated in Figure 4.

OPEN IN VIEWER

Fig. 4

Trajectories of S t α , I t α and R t α in example 3.

In the above figure, R 0 * have derived different values for α = 0.1, 0.2, ⋯ , 0.9 respectively. The trajectories of α-path S t α , I t α and R t α converge to E₀ (5, 0, 0) for R 0 * < 1 . Otherwise, we can see that lim t ⟶ ∞ I t α = 0.333 ( 1 - 1 R 0 * ) , It is obvious that the uncertainty of extinction and permanence of the epidemic mainly up to the value α.

Uncertain differential equation play a significant role in handling dynamic systems in epidemics. This paper concerned an uncertain SIRS epidemic model with bilinear incidence and demography, in which the uncertain disturbance was described by the Liu process. By considering the relationship between different compartments, we first established the uncertain SIRS model. Additionally, under the corresponding α-path, this paper defined an uncertain threshold function R 0 * which determines the extinction and persistence of an endemic. The result of stability analysis shows that when R 0 * ≤ 1 the disease dies out almost surely, while when R 0 * > 1 , the endemic equilibrium state is locally asymptotically stable. At last, the theoretical results proposed in this paper are verified by numerical simulation. There are some related problems which are yet not studied, such as globally asymptotically stable of endemic equilibria E₁ and uncertain parameter estimation for the uncertain SIRS infectious disease model,. These consideration will be taken as further investigation.