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

: 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 sufﬁcient 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 ﬁndings.


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 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) = ∑ n i=1 p i ϕ i (S, I), where p i (i = 1, 2, . . . , n) denotes the probability of ϕ i (S, I)(i = 1, 2, . . . , n) occurring, and ∑ n i=1 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 SI + (1 − p) β 2 SI 1 + mS + nI , 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 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 Hence the stochastic epidemic model of system (2) is given by With the following initial conditions: We can also easily show that Therefore, in this paper we can also discuss the following system for (I(t), R(t)) Precisely, the region Γ 1 = (S, I, R) ∈ R 3 , S + I + R < 1 is a positively invariant set of system (3) and the region Γ 2 = (I, R) ∈ 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≥0 , P ) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous, while F 0 contains all P -null sets). Let R 3 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 Sections 3 and 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.

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(t) represent the population size at time t. Thus, N(t) = S(t) + I(t) + R(t). For any n ≥ n 1 and t ∈ [0, τ n ), we have by derivation Then, by the comparison theorem, we obtain Clearly, V 1 (t) is nonnegative. According to the Itô formula, one has Here Λ is a generic positive constant. We compute Integrating over [0, τ n ∧ T] and taking the expectation, then EV 1 (S (τ n ∧ T) , I (τ n ∧ T) , R (τ n ∧ T)) ≤ V 1 (S(0), I(0), R(0)) + ΛT, which leads to Let n → ∞, and we obtain which is a contraction, then, τ ∞ = ∞ a.s. This completes the proof.

Extinction of the Disease
In this section, we consider the extinction of the disease of the system (4) under some sufficient assumptions.
Define the stop-time We suppose that E(τ 1 ) < ∞. Actually, we assume that Hence, for all T > 0 and t ≤ T ∧ τ 1 , one has Now, employing Itô's formula to log I, we have where (6) and the monotony of the functions Further, which leads to Then, letting T → ∞ and using Fatou's lemma, we have Therefore, our claim is true. In the next step, we shall show that I(t) goes by 1 θ 2 in a finite mean time. Therefore, we set By using (9), we can also have By induction, we have the following definitions We obtain Clearly, (τ n ) is an increasing sequence. Therefore, τ n → τ ∞ a.s. We denote Ω = ∞ n=1 (τ n < ∞), by (10) and we obtain We write Ω in the following form We show that P(τ ∞ < ∞) = 0. We assume that P(τ ∞ < ∞) > 0. Let ω ∈ (τ ∞ < ∞). We obtain for all n ∈ N * I(τ n (ω)) = 1 a n .
By extending n to ∞, we obtain which is a contradiction with P(τ 0 = ∞) = 1. Then, P(τ ∞ < ∞) = 0, furthermore This indicates that I(t) converges to 0 with probability 1. Using Fatou's lemma and the last equation of system (4), we obtain which implies with (12) that lim t→∞ R(t) = 0. a.s. This proof is therefore complete.

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

Proof.
Integrating the first equation of (3) between 0 and t, we obtain It is obvious that pβ 1 SI < pβ 1 , Multiplying both sides by 1 t , we obtain 1+mS+nI dB 2 (s) are continuous martingales with finite quadratic variation. By employing the strong law of large numbers for local martingales, we have Then, we define a function F on Γ 1 where B is the unique solution to the equation given by B = γ 1 η (µ+γ 1 +γ 2 ) 2 . Using (1 − I − R) 2 ≤ 1 and (8) we infer that Additionally, we also have Moreover, Combining (15)- (17), we have By using (14), we have Integrating the last inequality from 0 to t and multiplying by 1 t , it yields that F(S(t), I(t), R(t)) − F(S(0), I(0), R(0)) t where Clearly, M 1 (t),M 2 (t) are continuous martingales and M 1 (0) = 0, M 2 (0) = 0 with their quadratic variation Hence, by the strong law of large numbers for local martingales, we have Clearly, Combining (19), we find the desired result.
Integrating the last equation of (3) and multiplying by 1 t , we obtain Hence, we complete the proof.

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 0 s > 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 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 [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)) ∈ 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 ≥ 0, the operator P(t) can be defined as follows Hence, according to system (4), we define the integral Markov semigroup {P(t)} t≥0 .
Proof. Denote a(x) and b(x) two vectors fields defined on R n , then the Lie bracket [a, b](x) is presented as . . , n. We write the Itô system (4) as the Stratonovitch SDE where  ) 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 ) ∈ C ∞ (R, Γ 2 , Γ 2 ). Now we use the approach given in [24,25] to verify that K is positive. We fix a function φ ∈ C([0, T], R). Consider the system of differential equations below 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 The derivative D i 0 ,r 0 ,φ could be found as follows. Let , where a and b are the Jacobian of a(I, R) and b(I, R), respectively. Denote Q(t, t 0 ), for 0 ≤ t 0 ≤ t ≤ T, as the function satisfying Q(t 0 , t 0 ) = I, ∂Q ∂t (t, t 0 ) = E(t)Q(t, t 0 ). Then If, for some φ ∈ 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).
Proof. Let ∈ (0, T) and h = 1 [T− ,T] . By a first-order Taylor series approximation, we obtain By computation we have 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 ) ∈ Γ 2 and for almost every (I 1 , R 1 ) ∈ Γ 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 ) ∈ Γ 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 with positive solutions X 1 and Y 1 . We define the function R 1 (t), as follows . FurthermoreṘ Similarly, we consider X 2 and Y 2 the solutions of the ODEs Then, the function verifies the following Now, we choose a differentiable function, R 3 (t) = Z(t − T), where T is a sufficiently large number, such that where ε is a small enough positive number. Hence, from (27) and (28) we get that for all t ∈ [T − ε, T], Finally, we define a function In addition, we choose R 2 (t) appropriate to define the C 1 -function R φ as follows Hence the function R φ , constructed above, possesses the properties below Now, we define the function I φ (t) on the interval [0, T] by Thus, we can deduce from (32) that 0 < I φ (t) + R φ (t) < 1 for all t ∈ [0, T]. After that, we choose the control function φ(t) as follows . 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. If R 0 s > 1, then the semigroup {P(t)} t≥0 is asymptotically stable. This indicates that there exists a unique density K * (i, r), such that Proof. First, we construct a positive C 2 -function V 2 and a closed set ∆ ⊂ Γ 2 such that where L is the differential operator related to the system (4) and ∆ = (α 1 , α 2 ) × (α 3 , α 4 ) ⊂ Γ 2 in which α, α 1 , α 2 , α 3 , and α 4 are positive constants to be chosen later. Then the semigroup {P(t)} t≥0 is not sweeping from the set ∆. It is clear that the function f u,v (x) = u log x − vx, u, v > 0 has its maximum at u v . Using the definition of F(S, I, R) in (13), we can easily deduce that Hence, Combining (18) and the differential operator L, we have Let (I, R) ∈ Γ 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 • If a 1 ≤ I < 1 and R < a 3 , we choose a 3 to be sufficiently small, such that • If a 2 ≤ I < 1 − a 3 and a 3 ≤ R < 1 − a 2 , we choose a 2 < 1 to be sufficiently close to 1, such that • If a 1 ≤ I < 1 − a 4 and a 4 ≤ R < 1 − a 1 , we choose a 4 < 1 to be sufficiently close to 1 and a to be sufficiently small, such that Then, with the above appropriate choice on the five parameters, we conclude that as intended. We complete the proof of this theorem.

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) = ∑ n i=1 p i ϕ i (S, I), ∑ n i=1 p i = 1. We assume ϕ(S, I) = pβ 1 SI + (1 − p) β 2 SI 1 + mS + nI (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 0 s > 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 0 s = 1. -If there exists a threshold R * , when R 0 s < R * < R 0 , for the behavior of the system (3).
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.

Conflicts of Interest:
The authors declare no conflict of interest.