Open AccessFeature PaperArticle

Threshold Analysis and Stationary Distribution of a Stochastic Model with Relapse and Temporary Immunity

Peng Liu

Xinzhu Meng

and

Haokun Qi

College of Mathematics and Systems Science, Shandong University of Science and Technology, Qingdao 266590, China

Author to whom correspondence should be addressed.

Symmetry 2020, 12(3), 331; https://doi.org/10.3390/sym12030331

Submission received: 8 January 2020 / Revised: 17 February 2020 / Accepted: 19 February 2020 / Published: 26 February 2020

(This article belongs to the Special Issue Advances in Stochastic Differential Equations)

Download

Browse Figures

Versions Notes

Abstract

In this paper, a stochastic model with relapse and temporary immunity is formulated. The main purpose of this model is to investigate the stochastic properties. For two incidence rate terms, we apply the ideas of a symmetric method to obtain the results. First, by constructing suitable stochastic Lyapunov functions, we establish sufficient conditions for the extinction and persistence of this system. Then, we investigate the existence of a stationary distribution for this model by employing the theory of an integral Markov semigroup. Finally, the numerical examples are presented to illustrate the analytical findings.

Keywords:

stochastic epidemic model; extinction; persistence; stationary distribution

1. Introduction

Infectious diseases have always threatened the health of human beings. In recent decades, mathematicians have tried to study the spread of infectious diseases. Many disease models have been introduced to understand and analyze epidemics [1,2,3,4,5,6]. The simplest and best-known one is the SIR (susceptible, infections, removed) model, where the total population is divided into three compartments: susceptible (S), infections (I), and removed (R). If the epidemic models are conferred on temporary immunity, then we can establish SIS (susceptible, infections, susceptible) or SIRS (susceptible, infections, removed, susceptible) models. These models have been researched widely, from the deterministic to stochastic perspective [7,8,9,10,11,12,13,14,15,16]. However, the recovery of diseases may relapse, when latent infections reactivate and revert resulting in actively infected people. These types of diseases can be modeled by the SIRI (susceptible, infections, removed, infections) model [17,18,19,20], which reads as follows

\{\begin{matrix} \frac{d S}{d t} = μ - μ S - φ (S, I), \\ \frac{d I}{d t} = φ (S, I) - (μ + η) I + γ_{1} R, \\ \frac{d R}{d t} = η I - (μ + γ_{1}) R, \end{matrix}

(1)

where

μ

denotes the recruitment rate of humans, and we suppose that it is equal to the natural death rate of humans.

η

is the recovery rate,

γ_{1}

denotes the rate that recovered individuals are reverted to the infective state.

φ (S, I)

is the transmission function, which is important in order to analyze the stochastic properties of the model. In such a model, relapse is an important feature of these types of diseases that occurs in some animals and humans, for instance, tuberculosis and herpes.

Tudor [17] and Ding [20] established SIRI models with a bilinear incidence rate. However, in real life, the types of incidence rates of the same disease may be different as the environment changes. Nevertheless, on the basis of what we have studied, few works have been done on epidemic models with different incidence rates. Thus, we propose a model with a novel type of transmission function,

φ (S, I) = \sum_{i = 1}^{n} p_{i} φ_{i} (S, I)

, where

p_{i} (i = 1, 2, \dots, n)

denotes the probability of

φ_{i} (S, I) (i = 1, 2, \dots, n)

occurring, and

\sum_{i = 1}^{n} p_{i} = 1

Taking into account a bilinear incidence rate and the Beddington–DeAngelis incidence rate, in this regard, we assume

φ (S, I) = p β_{1} S I + (1 - p) \frac{β_{2} S I}{1 + m S + n I}

, where p represents the probability of a bilinear incidence rate occurring, the corresponding

1 - p

is the probability of the Beddington–DeAngelis incidence rate occurring, and evidently

0 < p < 1

. m and n are nonnegative constants.

β_{i} (i = 1, 2)

are positive constants, which represent the disease transmission coefficients. Hence, the corresponding model can be formulated as follows

\{\begin{matrix} d S = [μ - μ S - p β_{1} S I - (1 - p) \frac{β_{2} S I}{1 + m S + n I} + γ_{2} R] d t, \\ d I = [p β_{1} S I + (1 - p) \frac{β_{2} S I}{1 + m S + n I} - (μ + η) I + γ_{1} R] d t, \\ d R = [η I - (μ + γ_{1} + γ_{2}) R] d t, \end{matrix}

(2)

where

γ_{2}

denotes the per capita immune loss rate of removed.

Stochastic differential equations (SDEs) are ordinary differential equations (ODEs) that include random processes in their vector fields. In this paper, we establish an SDE system by introducing terms representing stochastic perturbations into the ODE system (2), which is achieved by letting

β_{i} \to β_{i} + σ_{i} B_{i} (i = 1, 2)

. Hence the stochastic epidemic model of system (2) is given by

\{\begin{cases} d S = [μ - μ S - p β_{1} S I - (1 - p) \frac{β_{2} S I}{1 + m S + n I} + γ_{2} R] d t - p σ_{1} S I d B_{1} - (1 - p) \frac{σ_{2} S I}{1 + m S + n I} d B_{2}, \\ d I = [p β_{1} S I + (1 - p) \frac{β_{2} S I}{1 + m S + n I} - (μ + η) I + γ_{1} R] d t + p σ_{1} S I d B_{1} + (1 - p) \frac{σ_{2} S I}{1 + m S + n I} d B_{2}, \\ d R = [η I - (μ + γ_{1} + γ_{2}) R] d t . \end{cases}

(3)

With the following initial conditions:

S (0) > 0

I (0) > 0

R (0) > 0

and

S (0) + I (0) + R (0) = 1

, we have

d (S + I + R) = [μ - μ (S + I + R)] d t .

We can also easily show that

S (t) + I (t) + R (t) = 1 for all t ⩾ 0 .

Therefore, in this paper we can also discuss the following system for

(I (t), R (t))

\{\begin{matrix} \begin{array}{l} d I = & [p β_{1} (1 - I - R) I + (1 - p) \frac{β_{_{2}} (1 - I - R) I}{1 + m (1 - I - R) + n I} - (μ + η) I + γ_{1} R] d t, \\ + p σ_{1} (1 - I - R) I d B_{1} + (1 - p) \frac{σ_{_{2}} (1 - I - R) I}{1 + m (1 - I - R) + n I} d B_{2}, \\ d R = & [η I - (μ + γ_{1} + γ_{2}) R] d t . \end{array} \end{matrix}

(4)

Precisely, the region

Γ_{1} = \{(S, I, R) \in R^{3}, S + I + R < 1\}

is a positively invariant set of system (3) and the region

Γ_{2} = \{(I, R) \in R^{2}, I + R < 1\}

is a positively invariant set of system (4).

The symmetric method for differential equations is a popular way to study ordinary differential equations (ODEs) and partial differential equations (PDEs). In this regard, we apply the ideas of this method shown in [21,22] to the SDE system.

Unless otherwise specified, throughout this paper, let

(Ω, F, {F}_{t \geq 0}, P)

be a complete probability space with a filtration

{\{F_{t}\}}_{t \geq 0}

satisfying the usual conditions (i.e., it is increasing and right continuous, while

F_{0}

contains all

P

-null sets). Let

R_{+}^{3} = \{x_{i} > 0, i = 1, 2, 3\}, R_{+}^{2} = \{x_{i} > 0, i = 1, 2\}

. Let

B (t)

be an

F_{t}

-adapted Brownian motion.

The rest of this paper is organized as follows. The existence and uniqueness of the positive solution of the system (3) are provided in Section 2. In Section 3 and Section 4, we explore sufficient conditions for extinction and persistence results of the epidemic. In Section 5, we use integral Markov semigroup theory to study the existence of a unique stationary distribution for system (4). Finally, we present some numerical simulations to demonstrate our analytical results.

2. Existence and Uniqueness of the Global Positive Solution

In the following, we discuss the existence and uniqueness of the global positive solution of system (3) for any initial condition on

R^{+}

. For further research, this is an important premise.

Theorem 1.

There is a unique and positive solution

(S (t), I (t), R (t))

of system (3) for any initial condition on

R_{+}

. Then the solution remains in

Γ_{1}

, almost surely.

Proof.

As the coefficients of system (3) are locally Lipschitz continuous, there exists a unique solution on

[0, τ_{e}]

, where

τ_{e}

is the explosion time.

Let

n_{0} \in N

, such that the initial values

S (0), I (0), R (0)

are all in the interval

(1 / n_{0}, n_{0})

. For any integer

n \geq n_{0}

, we consider the stop-time

τ_{n}

τ_{n} = inf \{t \in [0, τ_{e}) : min {S (t), I (t), R (t)} \leq \frac{1}{n} or max {S (t), I (t), R (t)} \geq n\} .

Clearly,

τ_{n}

is increasing as

n \to \infty

and less than

τ_{e}

for any

n \geq n_{0}

. Set

τ_{\infty} = {lim}_{n \to \infty} τ_{n}

. We have

τ_{\infty} \leq τ_{e}

. If

τ_{\infty} = \infty

a.s. holds, then

τ_{e} = \infty

a.s. Otherwise, there exist

ϵ \in (0, 1)

and

T \in (0, \infty)

, such that

P (τ_{\infty} \leq T) \geq ε

. Thus, there exists an integer

n_{1} \geq n_{0}

, such that

P (τ_{n} \leq T) \geq ϵ

for all

n \geq n_{1}

Let

N (t)

represent the population size at time t. Thus,

N (t) = S (t) + I (t) + R (t)

. For any

n \geq n_{1}

and

t \in [0, τ_{n})

, we have by derivation

d N (t) = [μ - μ N (t)] d t .

Then, by the comparison theorem, we obtain

N (t) \leq 1 + [N (0) - 1] e^{- μ t} \leq M = max {1, N (0)}, \forall t \in [0, τ_{k}) .

Define a

C^{2}

-function

V_{1}

V_{1} (S, I, R) = (S - 1 - log S) + (I - 1 - log I) + (R - 1 - log R) .

Clearly,

V_{1} (t)

is nonnegative. According to the Itô formula, one has

\begin{array}{r} L V_{1} = & (1 - \frac{1}{S}) [μ - μ S - p β_{1} S I - (1 - p) \frac{β_{2} S I}{1 + m S + n I} + γ_{2} R] \\ + (1 - \frac{1}{I}) [p β_{1} S I + (1 - p) \frac{β_{2} S I}{1 + m S + n I} - (μ + η) I + γ_{1} R] \\ + (1 - \frac{1}{R}) [η I - (μ + γ_{1} + γ_{2}) R] \\ + \frac{1}{2} p^{2} σ_{1}^{2} I^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} {(\frac{I}{1 + m S + n I})}^{2} \\ + \frac{1}{2} p^{2} σ_{1}^{2} S^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} {(\frac{S}{1 + m S + n I})}^{2} \\ = & 4 μ + η + γ_{1} + γ_{2} + p β_{1} I + (1 - p) \frac{β_{2} I}{1 + m S + n I} \\ + \frac{1}{2} p^{2} σ_{1}^{2} I^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} {(\frac{I}{1 + m S + n I})}^{2} \\ + \frac{1}{2} p^{2} σ_{1}^{2} S^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} {(\frac{S}{1 + m S + n I})}^{2} \\ - (μ S + μ I + μ R + \frac{μ}{S} + γ_{2} \frac{R}{S} + p β_{1} S + (1 - p) \frac{β_{2} S}{1 + m S + n I} \\ + γ_{1} \frac{R}{I} + η \frac{I}{R}) \\ \leq & 4 μ + η + γ_{1} + γ_{2} + p β_{1} + (1 - p) β_{2} + \frac{1}{2} p^{2} σ_{1}^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} = : Λ . \end{array}

Here

Λ

is a generic positive constant. We compute

d V_{1} = Λ d t + (\frac{1}{S} - \frac{1}{I}) p σ_{1} S I d B_{1} (t) + (\frac{1}{S} - \frac{1}{I}) (1 - p) \frac{σ_{2} S I}{1 + m S + n I} d B_{2} (t) .

Integrating over

[0, τ_{n} \land T]

and taking the expectation, then

E V_{1} (S (τ_{n} \land T), I (τ_{n} \land T), R (τ_{n} \land T)) \leq V_{1} (S (0), I (0), R (0)) + Λ T,

which leads to

\begin{array}{l} V_{1} (S (0), I (0), R (0)) + Λ T & \geq E [1_{\{τ_{n} \land T\}} V_{1} (S (τ_{n} \land T), I (τ_{n} \land T), R (τ_{n} \land T))] \\ \geq ε \{(n - 1 - log n) \land (\frac{1}{n} - 1 - log \frac{1}{n})\} . \end{array}

Let

n \to \infty

, and we obtain

\infty > V_{1} (S (0), I (0), R (0)) + Λ T = \infty,

which is a contraction, then,

τ_{\infty} = \infty

a.s. This completes the proof. □

3. Extinction of the Disease

In this section, we consider the extinction of the disease of the system (4) under some sufficient assumptions.

Theorem 2.

Assume that

R_{0} < 1,

(5)

β_{1} > p σ_{1}^{2}, β_{2} > (1 - p) σ_{2}^{2},

(6)

where

R_{0} = \frac{β_{1} p + β_{2} (1 - p)}{μ + η + \frac{1}{2} σ_{1}^{2} p^{2} + \frac{1}{2} σ_{2}^{2} {(1 - p)}^{2} - \frac{γ_{1} η}{μ + γ_{1} + γ_{2}}}

, then for any initial value

(I (0), R (0))

Γ_{2}

, the solution

(I (t), R (t))

of system (4) converges to

(0, 0)

a.s.

Proof.

Let

R_{0} < 1

and

θ > 1

, such that

- k_{θ} ≜ p β_{1} + (1 - p) β_{2} - (μ + η + \frac{1}{2} σ_{1}^{2} p^{2} + \frac{1}{2} σ_{2}^{2} {(1 - p)}^{2}) + \frac{θ γ_{1} η}{μ + γ_{1} + γ_{2}} < 0 .

Define the stop-time

τ_{1} = inf \{t \geq 0, I (t) \leq \frac{1}{θ}\} .

We suppose that

E (τ_{1}) < \infty

. Actually, we assume that

I (0) \in (\frac{1}{θ}, 1)

(if not

τ_{1} = 0

). Hence, for all

T > 0

and

t \leq T \land τ_{1}

, one has

I (t) \geq \frac{1}{θ} .

(7)

Now, employing Itô’s formula to

log I

, we have

\begin{array}{l} d log I = & [p β_{1} (1 - I - R) + (1 - p) \frac{β_{2} (1 - I - R)}{1 + m (1 - I - R) + n I} - (μ + η) + γ_{1} \frac{R}{I} \\ - \frac{1}{2} p^{2} σ_{1}^{2} {(1 - I - R)}^{2} - \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} {(\frac{1 - I - R}{1 + m (1 - I - R) + n I})}^{2}] d t \\ + p σ_{1} (1 - I - R) d B_{1} (t) + (1 - p) \frac{σ_{2} (1 - I - R)}{1 + m (1 - I - R) + n I} d B_{2} (t) \\ \leq & [- (μ + η) + β_{1} p - \frac{1}{2} p^{2} σ_{1}^{2} + β_{2} (1 - p) - \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} + θ γ_{1} R] d t \\ + p σ_{1} (1 - I - R) d B_{1} (t) + (1 - p) \frac{σ_{2} (1 - I - R)}{1 + m (1 - I - R) + n I} d B_{2} (t), \end{array}

(8)

where (6) and the monotony of the functions

x_{1} \mapsto - (μ + γ_{1} + γ_{2}) + p β_{1} x_{1} - \frac{1}{2} p^{2} σ_{1}^{2} x_{1}^{2}, x_{2} \mapsto - (μ + γ_{1} + γ_{2}) + (1 - p) \frac{β_{2} x_{2}}{1 + m x_{2} + n I} - \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} {(\frac{x_{2}}{1 + m x_{2} + n I})}^{2}

are used. Precisely, according to assumption (6), the functions

x_{1}

and

x_{2}

are both increasing on the interval

[0, 1]

. Combining (8) and the R-equation, one has

\begin{array}{l} d [log I + \frac{θ γ_{1}}{μ + γ_{1} + γ_{2}} R] \leq & [- (μ + η) + β_{1} p - \frac{1}{2} p^{2} σ_{1}^{2} + β_{2} (1 - p) - \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2} + \frac{θ γ_{1} η}{μ + γ_{1} + γ_{2}} I] d t \\ + p σ_{1} (1 - I - R) d B_{1} (t) + (1 - p) \frac{σ_{2} (1 - I - R)}{1 + m (1 - I - R) + n I} d B_{2} (t) \\ \leq & - k_{θ} d t + p σ_{1} (1 - I - R) d B_{1} (t) + (1 - p) \frac{σ_{2} (1 - I - R)}{1 + m (1 - I - R) + n I} d B_{2} (t) . \end{array}

Further,

E [log I (T \land τ_{1}) + \frac{θ γ_{1}}{μ + γ_{1} + γ_{2}} R (T \land τ_{1})] \leq log I (0) + \frac{θ γ_{1}}{μ + γ_{1} + γ_{2}} R (0) - k_{θ} E (T \land τ_{1}),

which leads to

E (T \land τ_{1}) \leq \frac{θ γ_{1}}{(μ + γ_{1} + γ_{2}) k_{θ}} - \frac{E [log I (T \land τ_{1})]}{k_{θ}} .

Then, letting

T \to \infty

and using Fatou’s lemma, we have

E (τ_{1}) \leq \frac{θ γ_{1}}{(μ + γ_{1} + γ_{2}) k_{θ}} + \frac{log θ}{k_{θ}} < \infty .

Therefore, our claim is true. In the next step, we shall show that

I (t)

goes by

\frac{1}{θ^{2}}

in a finite mean time. Therefore, we set

\begin{matrix} τ_{1}^{'} = inf \{t \geq τ_{1}, I (t) > \frac{1}{θ}\}, \\ τ_{2} = inf \{t < τ_{1}^{'}, I (t) < \frac{1}{θ^{2}}\} . \end{matrix}

Hence, for all

T > 0

and any

t \in [τ_{1}, T \land τ_{2}]

we obtain

\frac{1}{θ^{2}} \leq I (t) \leq \frac{1}{θ} .

(9)

By using (9), we can also have

\begin{array}{l} d [log I + \frac{θ^{2} γ_{1}}{μ + γ_{1} + γ_{2}} R] = & [β_{1} p + β_{2} (1 - p) - (μ + η + \frac{1}{2} p^{2} σ_{1}^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2}) + \frac{θ^{2} γ_{1} η}{μ + γ_{1} + γ_{2}} I] d t \\ + p σ_{1} (1 - I - R) d B_{1} (t) + (1 - p) \frac{σ_{2} (1 - I - R)}{1 + m (1 - I - R) + n I} d B_{2} (t) \\ \leq & - k_{θ} d t + p σ_{1} (1 - I - R) d B_{1} (t) + (1 - p) \frac{σ_{2} (1 - I - R)}{1 + m (1 - I - R) + n I} d B_{2} (t) . \end{array}

Integrating over

[τ_{1}, T \land τ_{2}]

, one can find

E [log I (T \land τ_{2})] \leq log I (τ_{1}) + \frac{θ^{2} γ_{1}}{μ + γ_{1} + γ_{2}} R (τ_{1}) - k_{θ} E (T \land τ_{2} - τ_{1}) .

Thus we have

E (τ_{2}) \leq E (τ_{1}) + \frac{θ^{2} γ_{1}}{(μ + γ_{1} + γ_{2}) k_{θ}} + \frac{2 log θ}{k_{θ}} < \infty .

By induction, we have the following definitions

\begin{matrix} τ_{n}^{'} = inf \{t \geq τ_{n - 1}, I (t) > \frac{1}{θ^{n - 1}}\}, \\ τ_{n} = inf \{t < τ_{n}^{'}, I (t) < \frac{1}{θ^{n}}\} . \end{matrix}

We obtain

E (τ_{n}) \leq E (τ_{n - 1}) + \frac{θ^{n} γ_{1}}{(μ + γ_{1} + γ_{2}) k_{θ}} + \frac{n log θ}{k_{θ}} < \infty .

(10)

Clearly, (

τ_{n}

) is an increasing sequence. Therefore,

τ_{n} \to τ_{\infty}

a.s. We denote

Ω^{'} = ⋂_{n = 1}^{\infty} (τ_{n} < \infty)

, by (10) and we obtain

P (Ω^{'}) = P (⋂_{n = 1}^{\infty} (τ_{n} < \infty)) = lim_{n \to \infty} P (τ_{n} < \infty) = 1 .

(11)

We write

Ω^{'}

in the following form

Ω^{'} = ((τ_{\infty} < \infty) \cap Ω^{'}) \cup ((τ_{\infty} = \infty) \cap Ω^{'}) = (τ_{\infty} < \infty) \cup ((τ_{\infty} = \infty) \cap Ω^{'}) .

We show that

P (τ_{\infty} < \infty) = 0

. We assume that

P (τ_{\infty} < \infty) > 0

. Let

ω \in (τ_{\infty} < \infty)

. We obtain for all

n \in N^{*}

I (τ_{n} (ω)) = \frac{1}{a^{n}} .

By extending n to ∞, we obtain

I (τ_{\infty} (ω)) = 0

. Hence,

τ_{0} (ω) \leq τ_{\infty} (ω) < \infty

, where

τ_{0} = inf \{t > 0, I (t) = 0\}

. Then

P (τ_{0} < \infty) \geq P (τ_{\infty} < \infty) > 0,

which is a contradiction with

P (τ_{0} = \infty) = 1

. Then,

P (τ_{\infty} < \infty) = 0

, furthermore

P ((τ_{\infty} = \infty) \cap Ω^{'}) = 1 .

Finally, we set

ω \in (τ_{\infty} = \infty) \cap Ω^{'}

ϵ > 0

t > 0

and

n_{0} = [- log ϵ {(log θ)}^{- 1}]

. For

t > τ_{n_{0}}

, there is n such that

n \geq n_{0}

and

τ_{n} \leq t \leq τ_{n + 1} .

Then

\frac{1}{θ^{n + 1}} \leq I (t) \leq \frac{1}{θ^{n}} \leq \frac{1}{θ^{n_{0} + 1}} \leq ϵ,

which suggests

lim_{t \to \infty} I (t) = 0 .

(12)

This indicates that

I (t)

converges to 0 with probability 1. Using Fatou’s lemma and the last equation of system (4), we obtain

\underset{t \to \infty}{lim sup} R (t) \leq η \int_{0}^{\infty} \underset{t \to \infty}{lim sup} e^{- (μ + γ_{1} + γ_{2}) u} I (t - u) d u \leq \frac{η}{μ + γ_{1} + γ_{2}} \underset{t \to \infty}{lim sup} I (t),

which implies with (12) that

{lim}_{t \to \infty} R (t) = 0 .

a.s. This proof is therefore complete. □

4. Persistence in the Mean

Theorem 3.

Let

(S (0), I (0), R (0)) \in Γ_{1})

be an initial value. Assume that

R_{s}^{0} > 1

holds, where

R_{s}^{0} = \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ + η - \frac{γ_{1} η}{μ + γ_{1} + γ_{2}} + \frac{1}{2} p^{2} σ_{1}^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2}}

. Then, the epidemic of (3) is permanent in the mean. Precisely, we obtain

\underset{t \to \infty}{lim inf} \frac{1}{t} \int_{0}^{t} S (u) d u \geq \frac{μ}{μ + p β_{1} + (1 - p) β_{2} C} a . s .,

\underset{t \to \infty}{lim inf} \frac{1}{t} \int_{0}^{t} I (u) d u \geq \frac{μ}{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}} (1 - \frac{1}{R_{s}^{0}}) a . s .,

\underset{t \to \infty}{lim inf} \frac{1}{t} \int_{0}^{t} R (u) d u \geq \frac{μ η}{(μ + γ_{1} + γ_{2}) (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n})} (1 - \frac{1}{R_{s}^{0}}) a . s .

Proof.

Integrating the first equation of (3) between 0 and t, we obtain

\begin{array}{l} μ t - (S (t) - S (0)) = & \int_{0}^{t} (μ S + p β_{1} S I + \frac{(1 - p) β_{2} S I}{1 + m S + n I} - γ_{2} R) d s \\ + \int_{0}^{t} p σ_{1} S I d B_{1} (s) + \int_{0}^{t} \frac{(1 - p) σ_{2} S I}{1 + m S + n I} d B_{2} (s) \\ \leq & \int_{0}^{t} (μ + p β_{1} I + \frac{(1 - p) β_{2} I}{1 + m S + n I}) S d s \\ + \int_{0}^{t} p σ_{1} S I d B_{1} (s) + \int_{0}^{t} \frac{(1 - p) σ_{2} S I}{1 + m S + n I} d B_{2} (s) . \end{array}

It is obvious that

p β_{1} S I < p β_{1}

\frac{S I}{1 + m S + n I} < C

, where C is a positive constant. Thus, we have

\begin{array}{l} μ t - (S (t) - S (0)) \leq & (μ + p β_{1} + (1 - p) β_{2} C) \int_{0}^{t} S d s \\ + \int_{0}^{t} p σ_{1} S I d B_{1} (s) + \int_{0}^{t} \frac{(1 - p) σ_{2} S I}{1 + m S + n I} d B_{2} (s) . \end{array}

Multiplying both sides by

\frac{1}{t}

, we obtain

\begin{array}{l} \frac{1}{t} \int_{0}^{t} S d s \geq & \frac{1}{μ + p β_{1} + (1 - p) β_{2} C} (μ - \frac{S (t) - S (0)}{t}) \\ - \frac{1}{μ + p β_{1} + (1 - p) β_{2} C} (\frac{1}{t} \int_{0}^{t} p σ_{1} S I d B_{1} (s) + \frac{1}{t} \int_{0}^{t} \frac{(1 - p) σ_{2} S I}{1 + m S + n I} d B_{2} (s)) . \end{array}

Note that

\int_{0}^{t} p σ_{1} S I d B_{1} (s)

\int_{0}^{t} \frac{(1 - p) σ_{2} S I}{1 + m S + n I} d B_{2} (s)

are continuous martingales with finite quadratic variation. By employing the strong law of large numbers for local martingales, we have

lim_{t \to \infty} \frac{S (t) - S (0)}{t} + \frac{1}{t} \int_{0}^{t} p σ_{1} S I d B_{1} (s) + \frac{1}{t} \int_{0}^{t} (1 - p) \frac{σ_{2} S I}{1 + m S + n I} d B_{2} (s) = 0 a . s .

Hence,

lim_{t \to \infty} \frac{1}{t} \int_{0}^{t} S d s \geq \frac{μ}{μ + p β_{1} + (1 - p) β_{2} C} a . s .

Then, we define a function F on

Γ_{1}

F (S, I, R) = log I + \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} S + B log R,

(13)

where B is the unique solution to the equation

- (μ + γ_{1} + γ_{2}) B + 2 \sqrt{γ_{1} η} \sqrt{B} - \frac{γ_{1} η}{μ + γ_{1} + γ_{2}} = 0,

(14)

given by

B = \frac{γ_{1} η}{{(μ + γ_{1} + γ_{2})}^{2}}

. Using

{(1 - I - R)}^{2} \leq 1

and (8) we infer that

\begin{array}{l} d log I \geq & [p β_{1} S + (1 - p) \frac{β_{2} S}{1 + m S + n I} - (μ + η + \frac{1}{2} p^{2} σ_{1}^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2}) + γ_{1} \frac{R}{I}] d t \\ + p σ_{1} S d B_{1} (t) + \frac{(1 - p) σ_{2} S}{1 + m S + n I} d B_{2} (t) . \end{array}

(15)

Additionally, we also have

\begin{array}{r} d (\frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} S) = & [p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n} - p β_{1} S - \frac{(1 - p) β_{2} S}{1 + m + n} - \frac{p^{2} β_{1}^{2} S I}{μ} - \frac{p (1 - p) β_{1} β_{2} S I \frac{S I}{1 + m + n}}{μ} \\ - \frac{p (1 - p) β_{1} β_{2} S I \frac{S I}{1 + m S + n I}}{μ} - \frac{\frac{β_{2}^{2} {(1 - p)}^{2} S I}{(1 + m + n) (1 + m S + n I)}}{μ} + \frac{γ_{2} (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}) R}{μ}] d t \\ - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} p σ_{1} S I d B_{1} (t) - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} \frac{(1 - p) σ_{2} S I}{1 + m S + n I} d B_{2} (t) \\ \geq & [p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n} - p β_{1} S - \frac{(1 - p) β_{2} S}{1 + m + n} - \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} I] d t \\ - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} p σ_{1} S I d B_{1} (t) - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} \frac{(1 - p) σ_{2} S I}{1 + m S + n I} d B_{2} (t) . \end{array}

(16)

Moreover,

d B log R = [η B \frac{I}{R} - (μ + γ_{1} + γ_{2}) B] d t .

(17)

Combining (15)–(17), we have

\begin{array}{l} d F (S, I, R) \geq & [p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n} - (μ + η + \frac{1}{2} p^{2} σ_{1}^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2}) - (μ + γ_{1} + γ_{2}) B \\ + {(\sqrt{\frac{γ_{1} R}{I}} - \sqrt{\frac{B η I}{R}})}^{2} + 2 \sqrt{γ_{1} η B} - \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} I] d t \\ + (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I) p σ_{1} S d B_{1} (t) + (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I) \frac{(1 - p) σ_{2} S}{1 + m S + n I} d B_{2} (t) \\ \geq & [p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n} - (μ + η + \frac{1}{2} p^{2} σ_{1}^{2} + \frac{1}{2} {(1 - p)}^{2} σ_{2}^{2}) - (μ + γ_{1} + γ_{2}) B \\ + 2 \sqrt{γ_{1} η B} - \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} I] d t \\ + (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I) p σ_{1} S d B_{1} (t) + (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I) \frac{(1 - p) σ_{2} S}{1 + m S + n I} d B_{2} (t) . \end{array}

By using (14), we have

\begin{array}{l} d F (S, I, R) \geq & [(p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}) (1 - \frac{1}{R_{s}^{0}}) - \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} I] d t \\ + (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I) p σ_{1} S d B_{1} (t) + (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I) \frac{(1 - p) σ_{2} S}{1 + m S + n I} d B_{2} (t) . \end{array}

(18)

Integrating the last inequality from 0 to t and multiplying by

\frac{1}{t}

, it yields that

\begin{array}{l} \frac{F (S (t), I (t), R (t)) - F (S (0), I (0), R (0))}{t} \geq & (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}) (1 - \frac{1}{R_{s}^{0}}) \\ - \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ t} \int_{0}^{t} I (u) d u + \frac{M_{1} (t)}{t} + \frac{M_{2} (t)}{t}, \end{array}

where

M_{1} (t) = p σ_{1} \int_{0}^{t} S (u) (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I (u)) d B_{1} (u),

and

M_{2} (t) = (1 - p) σ_{2} \int_{0}^{t} \frac{σ_{2} S (u)}{1 + m S (u) + n I (u)} (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I (u)) d B_{2} (u) .

Therefore

\begin{array}{l} \frac{1}{t} \int_{0}^{t} I (u) d u \geq & \frac{μ (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n})}{{(p β_{1} + (1 - p) β_{2})}^{2}} (1 - \frac{1}{R_{s}^{0}}) - \frac{μ}{{(p β_{1} + (1 - p) β_{2})}^{2}} \frac{F (S (t), I (t), R (t))}{t} \\ + \frac{μ}{{(p β_{1} + (1 - p) β_{2})}^{2}} \frac{F (S (0), I (0), R (0))}{t} \\ + \frac{μ}{{(p β_{1} + (1 - p) β_{2})}^{2}} \frac{M_{1} t}{t} + \frac{μ}{{(p β_{1} + (1 - p) β_{2})}^{2}} \frac{M_{2} t}{t} . \end{array}

Clearly,

M_{1} (t)

M_{2} (t)

are continuous martingales and

M_{1} (0) = 0

M_{2} (0) = 0

with their quadratic variation

[M_{1}] (t) = p^{2} {σ_{1}}^{2} \int_{0}^{t} {[S^{2} (u) (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I (u))]}^{2} d u < p^{2} {σ_{1}}^{2} (1 + \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ})^{2} t,

\begin{array}{l} [M_{2}] (t) = & {(1 - p)}^{2} {σ_{2}}^{2} \int_{0}^{t} (1 - \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} I (u))^{2} {(\frac{S (u)}{1 + m S (u) + n I (u)})}^{2} d u \\ < & {(1 - p)}^{2} {σ_{2}}^{2} (1 + \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ})^{2} t . \end{array}

Hence, by the strong law of large numbers for local martingales, we have

lim_{t \to \infty} \frac{1}{t} \frac{M_{1} (t)}{t} = 0 a . s .,

lim_{t \to \infty} \frac{1}{t} \frac{M_{2} (t)}{t} = 0 a . s .

Consequently, we have

lim_{t \to \infty} inf \frac{1}{t} \int_{0}^{t} I (u) d u \geq \frac{μ (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n})}{{(p β_{1} + (1 - p) β_{2})}^{2}} (1 - \frac{1}{R_{s}^{0}}) - \frac{μ}{{(p β_{1} + (1 - p) β_{2})}^{2}} lim_{t \to \infty} sup \frac{F (S (t), I (t), R (t))}{t} a . s .

(19)

Clearly,

lim_{t \to \infty} sup \frac{F (S (t), I (t), R (t))}{t} \leq lim_{t \to \infty} sup \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ t} = 0 .

Combining (19), we find the desired result.

lim_{t \to \infty} inf \frac{1}{t} \int_{0}^{t} I (u) d u \geq \frac{μ [p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}]}{{(p β_{1} + (1 - p) β_{2})}^{2}} (1 - \frac{1}{R_{s}^{0}}) a . s .

(20)

Integrating the last equation of (3) and multiplying by

\frac{1}{t}

, we obtain

\frac{μ + γ_{1} + γ_{2}}{t} \int_{0}^{t} R (s) d s = \frac{η}{t} \int_{0}^{t} I (s) d s + \frac{R (0) - R (t)}{t} .

lim_{t \to \infty} \frac{R (0) - R (t)}{t} = 0

, using (20), we deduce that

lim_{t \to \infty} inf \frac{1}{t} \int_{0}^{t} R (s) d s \geq \frac{η μ (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n})}{(μ + γ_{1} + γ_{2}) {(p β_{1} + (1 - p) β_{2})}^{2}} (1 - \frac{1}{R_{s}^{0}}) .

Hence, we complete the proof. □

5. Existence of a Stationary Distribution

In Section 4, we obtain that the persistence in the mean of the system starts from any initial condition of the invariant region

Γ_{1}

on certain conditions. In this section, to better understand the asymptotic behavior of the diseases shown in the SDE system (4), in the case of

R_{s}^{0} > 1

, we present that there exists a stationary distribution of the process

(I (t), R (t))

. The diffusion matrix of system (4) is shown as

\sum (i, r) = (\begin{matrix} {(p σ_{1} (1 - i - r) i + (1 - p) \frac{σ_{2} (1 - i - r) i}{1 + m (1 - i - r) + n i})}^{2} \\ 0 \end{matrix} \begin{matrix} 0 \\ 0 \end{matrix}) .

As the diffusion matrix is degenerate and dissatisfies the uniform ellipticity condition, we can not directly prove the asymptotic stability using Khasminskii’s theorem [23]. To handle this situation, we apply the integral Markov semigroup theory, presented in [14,24,25].

Now, we consider the space

(Γ_{2}, B (Γ_{2}), m)

, where

B (Γ_{2})

is the

σ

-algebra of Borel subsets of

Γ_{2}

and m is the Lebesgue measure on

(Γ_{2}, B (Γ_{2}))

. Throughout this section, we denote the transition probability function for the diffusion process

(I (t), R (t))

, i.e.,

P (t, i_{0}, r_{0}, B) = P ((I (t), R (t)) \in B | I (0) = i_{0}, R (0) = r_{0}) .

In Theorem 4, we show the absolute continuity of the transition function of the degenerate diffusion processes, as described by the Stratonovich equation, using the Hörmander condition, given in [26]. Let

K (t, i, r, i_{0}, r_{0})

be the density of

P (t, i_{0}, r_{0}, .)

. For any

t \geq 0

, the operator

P (t)

can be defined as follows

P (t) g (i, r) = \int_{Γ_{2}} K (t, i, r, i_{0}, r_{0}) g (i, r) m (d i, d r), g \in D .

Hence, according to system (4), we define the integral Markov semigroup

{P (t)}_{t \geq 0}

Theorem 4.

The transition probability function

P (t, i_{0}, r_{0}, .)

of the solution

(I (t), R (t))

of system (4), has a density

K (t, i, r, i_{0}, r_{0}) \in C^{\infty} (R, Γ_{2}, Γ_{2})

Proof.

Denote

a (x)

and

b (x)

two vectors fields defined on

R^{n}

, then the Lie bracket

[a, b] (x)

is presented as

{[a, b]}_{j} (x) = \sum_{k = 1}^{n} [a_{k} \frac{\partial b_{j}}{\partial x_{k}} (x) - b_{k} \frac{\partial a_{j}}{\partial x_{k}} (x)],

j = 1, \dots, n .

We write the Itô system (4) as the Stratonovitch SDE

(\begin{matrix} d I \\ d R \end{matrix}) = (\begin{matrix} a_{1} (I, R) \\ a_{2} (I, R) \end{matrix}) d t + (\begin{matrix} p σ_{1} (1 - I - R) I + (1 - p) \frac{σ_{2} (1 - I - R) I}{1 + m (1 - I - R) + n I} \\ 0 \end{matrix}) \circ d B (t),

(21)

where

\begin{matrix} a_{1} (I, R) = p β_{1} (1 - I - R) I + (1 - p) \frac{β_{2} (1 - I - R) I}{1 + m (1 - I - R) + n I} - (μ + η) I + γ_{1} R, \\ a_{2} (I, R) = η I - (μ + γ_{1} + γ_{2}) R . \end{matrix}

Let

a (I, R) = {(a_{1} (I, R), a_{2} (I, R))}^{'}

and

b (I, R) = {(p σ_{1} (1 - I - R) I + (1 - p) \frac{σ_{2} (1 - I - R) I}{1 + m (1 - I - R) + n I}, 0)}^{'}

. By direct calculation, we obtain

det ([a, b] (I, R), b (I, R)) = η {(p σ_{1} (1 - I - R) I + (1 - p) \frac{σ_{2} (1 - I - R) I}{1 + m (1 - I - R) + n I})}^{2} > 0

for any

[I, R] \in Γ_{2}

. Accordingly, the vectors

[a, b] (I, R)

and

b (I, R)

span the space

Γ 2_{}

. By employing the Hörmander theorem, the probability function

P (t, i_{0}, r_{0}, .)

has a density

K (t, i, r, i_{0}, r_{0}) \in C^{\infty} (R, Γ_{2}, Γ_{2})

. □

Now we use the approach given in [24,25] to verify that

K

is positive. We fix a function

ϕ \in C ([0, T], R)

. Consider the system of differential equations below

\{\begin{matrix} {\dot{I}}_{ϕ} (t) = a_{1} (I_{ϕ} (t), R_{ϕ} (t)) + (p σ_{1} (1 - I_{ϕ} (t) - R_{ϕ} (t)) I_{ϕ} (t) + (1 - p) \frac{σ_{2} (1 - I_{ϕ} (t) - R_{ϕ} (t)) I_{ϕ} (t)}{1 + m (1 - I_{φ} (t) - R_{φ} (t)) + n I_{φ} (t)} φ (t), \\ {\dot{R}}_{ϕ} (t) = a_{2} (I_{ϕ} (t), R_{ϕ} (t)), \end{matrix}

(22)

with the initial condition

I_{ϕ} (0) = i_{0}

R_{ϕ} (0) = r_{0}

and let

D_{i_{0}, r_{0}, ϕ}

be the Frechét derivative of the function

h \to (\begin{matrix} I_{ϕ + h} (T) \\ R_{ϕ + h} (T) \end{matrix})

from

C ([0, T], R)

R^{2}

. The derivative

D_{i_{0}, r_{0}, ϕ}

could be found as follows. Let

E (t) = a^{'} (I_{ϕ} (t), R_{ϕ} (t)) + ϕ (t) b^{'} (I_{ϕ} (t), R_{ϕ} (t))

, where

a^{'}

and

b^{'}

are the Jacobian of

a (I, R)

and

b (I, R)

, respectively. Denote

Q (t, t_{0})

, for

0 \leq t_{0} \leq t \leq T

, as the function satisfying

Q (t_{0}, t_{0}) = I

\frac{\partial Q}{\partial t} (t, t_{0}) = E (t) Q (t, t_{0})

. Then

D_{i_{0}, r_{0}, ϕ} h = \int_{0}^{T} Q (T, s) b (I_{ϕ} (s), R_{ϕ} (s)) h (s) d s .

If, for some

ϕ \in C ([0, T], R)

, the derivative

D_{i_{0}, r_{0}, ϕ}

has rank 2, then

K (t, i, r, i_{0}, r_{0}) > 0

for

i = I_{ϕ} (T)

and

r = R_{ϕ} (T)

Theorem 5.

For each

(i_{0}, r_{0}) \in Γ_{2}

and

ϕ \in C ([0, T], R)

, the derivative

D_{i_{0}, r_{0}, ϕ}

has rank 2.

Proof.

Let

ϵ \in (0, T)

and

h = 1_{[T - ϵ, T]}

. By a first-order Taylor series approximation, we obtain

Q (T, s) = 1 - E (T) (T - s) + o (T - s)

. Then

D_{i_{0}, r_{0}, ϕ} h = ϵ e - \frac{1}{2} ϵ^{2} E (T) e + o (ϵ^{2})

where

e = (\begin{matrix} 1 \\ 0 \end{matrix})

. By computation we have

E (T) e = (\begin{matrix} p (β_{1} + ϕ σ_{1}) φ_{ϕ} (t) + (1 - p) (β_{2} + ϕ σ_{2}) (\frac{φ_{ϕ} (t)}{1 + m ρ_{ϕ} (t)) + n I_{ϕ} (t)} - \frac{(n - m) ρ_{ϕ} (t) I_{ϕ} (t)}{{(1 + m ρ_{ϕ} (t)) + n I_{ϕ} (t))}^{2}} - (μ + η) \\ η \end{matrix}) .

where

φ_{ϕ} (t) = 1 - 2 I_{ϕ} (t) - R_{ϕ} (t)

and

ρ_{ϕ} (t) = 1 - I_{ϕ} (t) - R_{ϕ} (t)

. Therefore, e and

E (t) e

are linearly independent and the derivative

D_{i_{0}, r_{0}, ϕ}

has rank 2. □

Now we verify the positivity of

K

. Hence, we investigate the controllability of the region

Γ_{2}

Theorem 6.

For each

(I_{0}, R_{0}) \in Γ_{2}

and for almost every

(I_{1}, R_{1}) \in Γ_{2}

, there exists

T > 0

, such that

K (T, I_{1}, R_{1}, I_{0}, R_{0}) > 0

Proof.

Fixed

(i_{0}, r_{0})

(i_{1}, r_{1}) \in Γ_{2}

. Without a loss of generality, we assume that

r_{1} < r_{0}

. Now we show that there is a control function

ϕ

and

T > 0

such that

I_{ϕ} (0) = i_{0}

R_{ϕ} (0) = r_{0}

I_{ϕ} (T) = i_{1}

R_{ϕ} (T) = r_{1}

. We first consider the following ODEs

\{\begin{matrix} {\dot{X}}_{1} = - (μ + γ_{1} + γ_{2}) X_{1}, \\ {\dot{Y}}_{1} = η - (μ + γ_{1} + γ_{2} + η) Y_{1}, \\ X_{1} (0) = Y_{1} (0) = r_{0}, \end{matrix}

(23)

with positive solutions

X_{1}

and

Y_{1}

. We define the function

R_{1} (t)

, as follows

R_{1} (t) = X_{1} (t) + \frac{i_{0}}{1 - r_{0}} (Y_{1} (t) - X_{1} (t)) .

Thus, we have

R_{1} (0) = r_{0}, {\dot{R}}_{1} (0) = η i_{0} - (μ + γ_{1} + γ_{2}) .

(24)

Furthermore

\begin{array}{l} {\dot{R}}_{1} (t) & = - (μ + γ_{1} + γ_{2}) R_{1} (t) + \frac{λ i_{0}}{1 - r_{0}} (1 - Y_{1} (t)), \\ = η - (μ + γ_{1} + γ_{2} + η) R_{1} (t) + \frac{η (1 - i_{0} - r_{0})}{1 - r_{0}} (X_{1} (t) - 1) . \end{array}

X_{1} (t), Y_{1} (t) \in (0, 1)

for all

t \geq 0

. Then

- (μ + γ_{1} + γ_{2}) R_{1} (t) < {\dot{R}}_{1} (t) < η - (μ + γ_{1} + γ_{2} + η) R_{1} (t) .

(25)

Similarly, we consider

X_{2}

and

Y_{2}

the solutions of the ODEs

\{\begin{matrix} {\dot{X}}_{2} = - (μ + γ_{1} + γ_{2}) X_{2}, \\ {\dot{Y}}_{2} = η - (μ + γ_{1} + γ_{2} + η) Y_{2}, \\ X_{2} (0) = Y_{2} (0) = r_{1} . \end{matrix}

(26)

Then, the function

Z (t) = X_{2} (t) + \frac{i_{1}}{1 - r_{1}} (Y_{2} (t) - X_{2} (t)),

verifies the following

Z (0) = r_{1}, \dot{Z} (0) = η i_{1} - (μ + γ_{1} + γ_{2}) r_{1} .

(27)

Now, we choose a differentiable function,

R_{3} (t) = Z (t - T)

, where T is a sufficiently large number, such that

X_{1} (t) \leq R_{3} (t) \leq Y_{1} (t), \forall t \in [T - ε, ε],

(28)

where

ε

is a small enough positive number. Hence, from (27) and (28) we get that for all

t \in [T - ε, T]

R_{3} (t)

verifies

- (μ + γ_{1} + γ_{2}) R_{3} (t) < {\dot{R}}_{3} (t) < η - (μ + γ_{1} + γ_{2} + η) R_{3} (t),

(29)

R_{2} (T) = r_{1}, {\dot{R}}_{2} (T) = η i_{1} - (μ + γ_{1} + γ_{2}) r_{1} .

(30)

Finally, we define a function

R_{2} (t)

[ε, T - ε]

such that

X_{1} (t) \leq R_{2} (t) \leq Y_{1} (t)

, then we have

- (μ + γ_{1} + γ_{2}) R_{2} (t) < {\dot{R}}_{2} (t) < η - (μ + γ_{1} + γ_{2} + η) R_{2} (t) .

(31)

In addition, we choose

R_{2} (t)

appropriate to define the

C^{1}

-function

R_{ϕ}

as follows

R_{ϕ} (t) = \{\begin{matrix} R_{1} (t), 0 \leq t \leq ε, \\ R_{2} (t), ε \leq t \leq T - ε, \\ R_{3} (t), T - ε \leq t \leq T . \end{matrix}

Hence the function

R_{ϕ}

, constructed above, possesses the properties below

\begin{array}{l} R_{ϕ} (0) = r_{0}, R_{ϕ} (t) = r_{1}, \\ {\dot{R}}_{ϕ} (0) = η i_{0} - (μ + γ_{1} + γ_{2}) r_{0}, {\dot{R}}_{ϕ} (T) = η i_{1} - (μ + γ_{1} + γ_{2}) r_{1}, \\ - (μ + γ_{1} + γ_{2}) R_{ϕ} (t) < {\dot{R}}_{ϕ} (t) < η - (μ + γ_{1} + γ_{2} + η) R_{ϕ} (t) . \end{array}

(32)

Now, we define the function

I_{ϕ} (t)

on the interval

[0, T]

I_{ϕ} (t) = \frac{1}{η} ({\dot{R}}_{ϕ} (t) + (μ + γ_{1} + γ_{2}) R_{ϕ} (t)) .

Thus, we can deduce from (32) that

0 < I_{ϕ} (t) + R_{ϕ} (t) < 1

for all

t \in [0, T]

. After that, we choose the control function

ϕ (t)

as follows

ϕ (t) = \frac{{\dot{I}}_{ϕ} (t) - a_{1} (I_{ϕ} (t), R_{ϕ} (t))}{[(p σ_{1} (1 - I_{ϕ} (t) - R_{ϕ} (t)) + (1 - p) \frac{σ_{2} (1 - I_{ϕ} (t) - R_{ϕ} (t))}{1 + m (1 - I_{ϕ} (t) - R_{ϕ} (t)) + n I_{ϕ} (t)}] I_{ϕ} (t)} .

Therefore, with the above determinate continuous function

ϕ

(I_{ϕ} (t), R_{ϕ} (t))

is the solution to system (22). Furthermore, from (32) and the chosen of

I_{ϕ} (t)

, we find that this solution satisfies

(I_{ϕ} (0), R_{ϕ} (0)) = (i_{0}, r_{0}), (I_{ϕ} (T), R_{ϕ} (T)) = (i_{1}, r_{1}) .

This completes the proof. □

Theorem 7.

R_{s}^{0} > 1

, then the semigroup

{P (t)}_{t \geq 0}

is asymptotically stable. This indicates that there exists a unique density

K^{*} (i, r)

, such that

\int \int_{Γ_{2}} | K (t, i, r, i_{0}, r_{0}) - K^{*} (i, r) | d i d r = 0 f o r a l l (i_{0}, r_{0}) \in Γ_{2} .

Proof.

First, we construct a positive

C^{2}

-function

V_{2}

and a closed set

Δ \subset Γ_{2}

such that

sup_{(I, R) \in Γ_{2} \ Δ} L V_{2} (I, R) < 0,

where

L

is the differential operator related to the system (4) and

Δ = (α_{1}, α_{2}) \times (α_{3}, α_{4}) \subset Γ_{2}

in which

α

α_{1}

α_{2}

α_{3}

, and

α_{4}

are positive constants to be chosen later. Then the semigroup

{P (t)}_{t \geq 0}

is not sweeping from the set

Δ

. It is clear that the function

f_{u, v} (x) = u log x - v x, u, v > 0

has its maximum at

\frac{u}{v}

. Using the definition of

F (S, I, R)

in (13), we can easily deduce that

\begin{array}{l} F (I, R) + a log R & = log I + \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} (1 - I - R) + B log R + a log R \\ = f_{1, \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ}} (I) + f_{B + a, \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ}} (R) + \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} \\ \leq f_{1, \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ}} (\frac{μ}{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}) + f_{B + a, \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ}} (\frac{μ (B + a)}{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}) + \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} . \end{array}

Hence,

\begin{array}{l} V_{2} (I, R) = & f_{1, \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ}} (\frac{μ}{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}) + f_{B + a, \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ}} (\frac{μ (B + a)}{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}) + \frac{p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}}{μ} + 1 \\ - F (I, R) - a log R > 0 . \end{array}

Combining (18) and the differential operator

L

, we have

L V_{2} (I, R) \leq - (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}) (1 - \frac{1}{R_{s}^{0}}) + \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} I + a (μ + γ_{1} + γ_{2}) - \frac{a η I}{R} .

(33)

Let

(I, R) \in Γ_{2} \ Δ

. We consider the following four cases:

If $I < a_{1}$ . In this case, we choose $a_{1}$ and a to be sufficiently small, such that

$L V_{2} (I, R) \leq - (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}) (1 - \frac{1}{R_{s}^{0}}) + \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} a_{1} + a (μ + γ_{1} + γ_{2}) < 0 .$
If $a_{1} \leq I < 1$ and $R < a_{3}$ , we choose $a_{3}$ to be sufficiently small, such that

$L V_{2} (I, R) \leq \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} a_{1} + a (μ + γ_{1} + γ_{2}) - \frac{a η a_{1}}{a_{3}} < 0 .$
If $a_{2} \leq I < 1 - a_{3}$ and $a_{3} \leq R < 1 - a_{2}$ , we choose $a_{2} < 1$ to be sufficiently close to 1, such that

$L V_{2} (I, R) \leq \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} a_{1} + a (μ + γ_{1} + γ_{2}) - \frac{a η a_{2}}{1 - a_{2}} < 0 .$
If $a_{1} \leq I < 1 - a_{4}$ and $a_{4} \leq R < 1 - a_{1}$ , we choose $a_{4} < 1$ to be sufficiently close to 1 and a to be sufficiently small, such that

$L V_{2} (I, R) \leq - (p β_{1} + \frac{(1 - p) β_{2}}{1 + m + n}) (1 - \frac{1}{R_{s}^{0}}) + \frac{{(p β_{1} + (1 - p) β_{2})}^{2}}{μ} a_{1} (1 - a_{4}) + a (μ + γ_{1} + γ_{2}) < 0 .$

Then, with the above appropriate choice on the five parameters, we conclude that

L V_{2} (I, R) < 0, \forall (I, R) \in Γ_{2} \ Δ,

as intended. We complete the proof of this theorem. □

6. Numerical Simulations

We simulate the system (3) for various parameter sets to verify the main theoretical results raised in this paper. For simplicity, we choose the following parameters

η = 0.1

p = 0.3

Example 1.

Figure 1 presents a single path of the solution

(S (t), I (t), R (t))

to the stochastic system (3) with initial condition

(0.7, 0.2, 0.1)

for the parameters

β_{1} = β_{2} = 0.3

σ_{1} = σ_{2} = 0.5

γ_{1} = 0.04

γ_{2} = 0.01

m = 0.25

n = 0.3

, where

R_{0} = 0.8415 < 1

. Then the disease of system (4) will die out.

Example 2.

Figure 2 presents the solution

(S (t), I (t), R (t))

to system (3) with an initial conditions

(0.7, 0.2, 0.1)

for the parameters

μ = 0.3

β_{1} = β_{2} = 0.6

σ_{1} = σ_{2} = 0.2

γ_{1} = 0.2

γ_{2} = 0.1

m = 0.15

n = 0.2

, where

R_{s}^{0} = 1.2983 > 1

. Then from Theorem 4 the disease of system (3) is permanent in mean, and admits a stationary distribution.

Example 3.

Figure 3 presents the results of SDE (3) with initial conditions

(0.7, 0.2, 0.1)

for

γ_{1} = 0, 0.5, 1

, and the other parameters are the same as in Example 2. Then we have

1 < R_{s}^{0} (0) = 1.1932 < R_{s}^{0} (0.5) = 1.3794 < R_{s}^{0} (1) = 1.4437

. Therefore, a relapse can increase the incidence of the disease.

7. Conclusions

In this paper, we focus on the principle properties of a stochastic epidemic model with relapse and temporary immunity. The highlight of this paper is that, considering that the incidence rate of diseases may change as the environment changes, we put forward a new version of the transmission function

φ (S, I) = \sum_{i = 1}^{n} p_{i} φ_{i} (S, I)

\sum_{i = 1}^{n} p_{i} = 1

. We assume

φ (S, I) = p β_{1} S I + (1 - p) \frac{β_{2} S I}{1 + m S + n I}

(0 < p < 1)

, and this transmission function indicates that, in different environments, the disease has two different incidence rates, which are the bilinear incidence rate and the Beddington–DeAngelis incidence rate. Moreover, if the incidence rate of this disease is only the bilinear incidence rate, then

p = 1

. In contrast, if

p = 0

, this implies that the disease has only one incidence, which is the Beddington–DeAngelis incidence rate. Therefore, our assumptions are both mathematically and biologically reasonable. The extinction result shows that the disease will die out almost surely when

R_{0} < 1

. Then we obtain that the disease persists in the mean when

R_{s}^{0} > 1

, and there exists a stationary distribution for the stochastic model in the meantime (see Figure 2). We also mention that a relapse can increase the risk of diseases, as

R_{0}

increases with the relapse rate

γ_{1}

(see Figure 3).

There are also many related problems to be solved with further investigation. For example:

-: The long-time behavior of the system (3) when $R_{s}^{0} = 1$ .
-: If there exists a threshold $R^{*}$ , when $R_{s}^{0} < R^{*} < R_{0}$ , for the behavior of the system (3).

In order to make this model more reasonable, we can also investigate delay stochastic differential equation models, control stochastic differential equation models, and impulsive stochastic differential equation models of system (2) for further work. The approaches are shown in [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41].

Author Contributions

All authors contributed equally and significantly in writing this article. All authors have read and agreed to the published version of the manuscript.

Funding

This project was funded by the Research Fund for the Taishan Scholar Project of Shandong Province of China, Shandong Provincial Natural Science Foundation of China (No. ZR2019MA003), and SDUST Innovation Fund for Graduate Students (No. SDKDYC190119).

Conflicts of Interest

The authors declare no conflict of interest.

References

Anderson, R.; May, R.; Medley, G. A preliminary study of the transmission dynamics of the human immunodeficiency virus (HIV), the causative agent of AIDS. IMA J. Math. Appl. Med. Biol. 1986, 3, 229–263. [Google Scholar] [CrossRef] [PubMed]
Ma, Z.E.; Zhou, Y.C.; Wu, J.H. Modeling and Dynamics of Infectious Diseases; Higher Education Press: Beijing, China, 2009. [Google Scholar]
Herbert, H.W. The mathematics of infectious diseases. SIAM Rev. 2000, 42, 599–653. [Google Scholar]
Brauer, F.; Chavez, C.C. Mathematical Models in Population Biology and Epidemiology; Springer: New York, NY, USA, 2001. [Google Scholar]
Liu, K.Y.; Zhang, T.Q.; Chen, L.S. State-dependent pulse vaccination and therapeutic strategy in an SI epidemic model with nonlinear incidence rate. Comput. Math. Methods Med. 2019, 2019, 3859815. [Google Scholar] [CrossRef] [PubMed]
Zhao, W.C.; Liu, J.L.; Chi, M.N.; Bian, F.F. Dynamics analysis of stochastic epidemic models with standard incidence. Adv. Differ. Equ. 2019, 2019, 22. [Google Scholar] [CrossRef] [Green Version]
Qi, H.K.; Leng, X.N.; Meng, X.Z.; Zhang, T.H. Periodic Solution and Ergodic Stationary Distribution of SEIS Dynamical Systems with Active and Latent Patients. Qual. Theory Dyn. Syst. 2019, 18, 347–369. [Google Scholar] [CrossRef]
Wang, X.H.; Wang, Z.; Shen, H. Dynamical analysis of a discrete-time SIS epidemic model on complex networks. Appl. Math. Lett. 2019, 94, 292–299. [Google Scholar] [CrossRef]
Korobeinikov, A. Lyapunov functions and global stability for SIR and SIRS epidemiological models with nonlinear transmission. Bull. Math. Biol. 2006, 30, 615–636. [Google Scholar] [CrossRef] [Green Version]
Liu, Q.; Jiang, D.Q.; Hayat, T.; Alsaedi, A. Threshold behavior in a stochastic delayed SIS epidemic model with vaccination and double diseases. J. Franklin Inst. 2019, 356, 7466–7485. [Google Scholar] [CrossRef]
Meng, X.Z.; Zhao, S.N.; Feng, T.; Zhang, T.H. Dynamics of a novel nonlinear atochastic SIS epidemic model with double epidemic hypothesis. J. Math. Anal. Appl. 2016, 433, 227–242. [Google Scholar] [CrossRef]
Chang, Z.B.; Meng, X.Z.; Zhang, T.H. A new way of investigating the asymptotic behaviour of a stochastic SIS system with multiplicative noise. Appl. Math. Lett. 2019, 87, 80–86. [Google Scholar] [CrossRef]
Gao, N.; Song, Y.; Wang, X.Z.; Liu, J.X. Dynamics of a stochastic SIS epidemic model with nonlinear incidence rates. Adv. Differ. Equ. 2019, 2019, 41. [Google Scholar] [CrossRef] [Green Version]
Fatini, M.E.; Khalifi, M.E.; Gerlach, R.; Laaribi, A.; Taki, R. Stationary distribution and threshold dynamics of a stochastic SIRS model with a general incidence. Phys. A 2019, 534, 120696. [Google Scholar] [CrossRef]
Song, Y.; Miao, A.Q.; Zhang, T.Q. Extinction and persistence of a stochastic SIRS epidemic model with saturated incidence rate and transfer from infectious to susceptible. Adv. Differ. Equ. 2018, 2018, 293. [Google Scholar] [CrossRef] [Green Version]
Liu, Q.; Jiang, D.Q.; Hayat, T.; Alsaedi, A. Threshold dynamics of a stochastic SIS epidemic model with nonlinear incidence rate. Phys. A 2019, 526, 120946. [Google Scholar] [CrossRef]
Tudor, D. A deterministic model for herpes infections in human and animal polulations. SIAM Rev. 1990, 32, 130–139. [Google Scholar] [CrossRef]
Zhang, W.W.; Meng, X.Z.; Dong, Y.L. Periodic Solution and Ergodic Stationary Distribution of Stochastic SIRI Epidemic Systems with Nonlinear Perturbations. J. Syst. Sci. Complex. 2019, 32, 1104–1124. [Google Scholar] [CrossRef]
Blower, S. Modeling the genital herpes epidemic. Herpes 2004, 11 (Suppl. 3), 138A. [Google Scholar]
Ding, S.S.; Wang, F.J. Sili epidemiological model with nonlinear incidence rates. J. Biomath. 1994, 9, 1–59. [Google Scholar]
Dorodnitsyn, V. Applications of Lie Groups to Difference Equations; Chapman and Hall/CRC: Boca Raton, FL, USA, 2010. [Google Scholar]
Stephani, H. Differential Equations: Their Solution Using Symmetries; Cambridge University Press: Cambridge, UK, 1989. [Google Scholar]
Khasminskii, R. Stochastic Stability of Differential Equations; Springer: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
Rudnicki, R. Long-time behaviour of a stochastic prey-predator model. Stoch. Proc. Appl. 2003, 108, 93–107. [Google Scholar] [CrossRef] [Green Version]
Rudnicki, R.; Pichor, K. Influence of stochastic perturbation on prey-predator systems. Math. Biosci. 2007, 206, 108–119. [Google Scholar] [CrossRef]
Bell, D.R. The Malliavin Calculus; Dover Publications: New York, NY, USA, 2006. [Google Scholar]
Qi, H.K.; Meng, X.Z.; Chang, Z.B. Markov semigroup approach to the analysis of a nonlinear stochastic plant disease model. Electron. J. Differ. Equ. 2019, 2019, 1–19. [Google Scholar]
Yang, Y.; Xia, J.W.; Zhao, J.L.; Li, X.D.; Wang, Z. Multiobjective nonfragile fuzzy control for nonlinear stochastic financial systems with mixed time delays. Nonlinear Anal. Model. Control 2019, 24, 696–717. [Google Scholar] [CrossRef]
Hou, J.; Zhao, Y. Some remarks on a pair of seemingly unrelated regression models. Open Math. 2019, 17, 979–989. [Google Scholar] [CrossRef]
Wang, F.; Liu, Z.; Zhang, Y.; Chen, C.L.P. Adaptive finite-time control of stochastic nonlinear systems with actuator failures. Fuzzy Sets Syst. 2019, 374, 170–183. [Google Scholar] [CrossRef]
Wang, Y.F.; Cheng, H.D.; Li, Q.J. Dynamic analysis of wild and sterile mosquito release model with Poincare map. Math. Biosci. Eng. 2019, 16, 7688–7706. [Google Scholar] [CrossRef] [PubMed]
Liu, G.D.; Chang, Z.B.; Meng, X.Z. Asymptotic analysis of impulsive dispersal predator-prey systems with Markov switching on finite-state space. J. Funct. Spaces 2019, 2019, 8057153. [Google Scholar] [CrossRef] [Green Version]
Zhang, J.; Xia, J.W.; Sun, W.; Zhuang, G.M.; Wang, Z. Finite-time tracking control for stochastic nonlinear systems with full state constraints. Appl. Math. Comput. 2018, 338, 207–220. [Google Scholar] [CrossRef]
Zhao, Y.; Zhang, T.L.; Fu, Y.; Ma, L.M. Finite-Time Stochastic H_∞ Control for Singular Markovian Jump Systems With (x,v)-Dependent Noise and Generally Uncertain Transition Rates. IEEE Access 2019, 7, 64812–64826. [Google Scholar] [CrossRef]
Wang, F.; Chen, B.; Sun, Y.M.; Lin, C. Finite time control of switched stochastic nonlinear systems. Fuzzy Sets Syst. 2019, 365, 140–152. [Google Scholar] [CrossRef]
Shi, Z.Z.; Li, Y.N.; Cheng, H.D. Dynamic analysis of a pest management smith model with impulsive state feedback control and continuous delay. Mathematics 2019, 7, 591. [Google Scholar] [CrossRef] [Green Version]
Li, F.; Zhang, S.Q.; Meng, X.Z. Dynamics analysis and numerical simulations of a delayed stochastic epidemic model subject to a general response function. Comput. Appl. Math. 2019, 38, 95. [Google Scholar] [CrossRef]
Liu, Q.; Jiang, D.Q. The threshold of a stochastic delayed SIR epidemic model with vaccination. Phys. A 2016, 461, 140–147. [Google Scholar] [CrossRef]
Gao, S.J.; Chen, L.S.; Nieto, J.J.; Torres, A. Analysis of a delayed epidemic model with pulse vaccination and saturation incidence. Vaccine 2006, 24, 6037–6045. [Google Scholar] [CrossRef]
Li, Y.J.; Meng, X.Z. Dynamics of an impulsive stochastic nonautonomous chemostat model with two different growth rates in a polluted environment. Discret. Dyn. Nat. Soc. 2019, 2019, 15. [Google Scholar] [CrossRef] [PubMed]
Zhu, C.; Yin, G. Asymptotic properties of hybrid diffusion systems. SIAM J. Control Optim. 2007, 46, 1155–1179. [Google Scholar] [CrossRef]

Figure 1. Paths of the solution and the kernel density functions to the system (3) using the data of Example 2.

Figure 2. Paths of the solution to the system (3) using the data of Example (1).

Figure 3. Paths of the solution to the system (3) using the data of Example 3.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Liu, P.; Meng, X.; Qi, H. Threshold Analysis and Stationary Distribution of a Stochastic Model with Relapse and Temporary Immunity. Symmetry 2020, 12, 331. https://doi.org/10.3390/sym12030331

AMA Style

Liu P, Meng X, Qi H. Threshold Analysis and Stationary Distribution of a Stochastic Model with Relapse and Temporary Immunity. Symmetry. 2020; 12(3):331. https://doi.org/10.3390/sym12030331

Chicago/Turabian Style

Liu, Peng, Xinzhu Meng, and Haokun Qi. 2020. "Threshold Analysis and Stationary Distribution of a Stochastic Model with Relapse and Temporary Immunity" Symmetry 12, no. 3: 331. https://doi.org/10.3390/sym12030331

APA Style

Liu, P., Meng, X., & Qi, H. (2020). Threshold Analysis and Stationary Distribution of a Stochastic Model with Relapse and Temporary Immunity. Symmetry, 12(3), 331. https://doi.org/10.3390/sym12030331

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Threshold Analysis and Stationary Distribution of a Stochastic Model with Relapse and Temporary Immunity

Abstract

1. Introduction

2. Existence and Uniqueness of the Global Positive Solution

3. Extinction of the Disease

4. Persistence in the Mean

5. Existence of a Stationary Distribution

6. Numerical Simulations

7. Conclusions

Author Contributions

Funding

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI