Threshold dynamics and sensitivity analysis of a stochastic semi-Markov switched SIRS epidemic model with nonlinear incidence and vaccination

In this paper, a stochastic SIRS epidemic model with nonlinear incidence and vaccination is formulated to investigate the transmission dynamics of infectious diseases. The model not only incorporates the white noise but also the external environmental noise which is described by semi-Markov process. We first derive the explicit expression for the basic reproduction number of the model. Then the global dynamics of the system is studied in terms of the basic reproduction number and the intensity of white noise, and sufficient conditions for the extinction and persistence of the disease are both provided. Furthermore, we explore the sensitivity analysis of \begin{document}$ R_0^s $\end{document} with each semi-Markov switching under different distribution functions. The results show that the dynamics of the entire system is not related to its switching law, but has a positive correlation to its mean sojourn time in each subsystem. The basic reproduction number we obtained in the paper can be applied to all piecewise-stochastic semi-Markov processes, and the results of the sensitivity analysis can be regarded as a prior work for optimal control.


1.
Introduction. Since the seminal work of Kermack and McKendrick [1], mathematical modeling has been an important tool in understanding the transmission dynamics of infectious diseases. In recent years, with the help of mathematical modeling, much significant progress has been made in epidemiology [2,3,28,38,39,47,49]. Most research on the dynamics of infectious diseases is based on the mathematical framework of deterministic differential equations. These models provide us with a far-sighted theoretical guidance for understanding the transmission mechanism of infectious diseases [4,17,18,30].
However, population dynamics is often affected by various types of external environmental noises, such as the changes in weather, temperature and climate [41,44,45]. Since these fluctuations may cause dramatic changes in some important parameters of infectious diseases (e.g., transmission rate, death rate, recovery rate), they are generally described as switching between two or more regimes of the environment. For example, the transmission rate of mosquito-borne disease is completely different in summer and winter [19]. In the literature, the switching can be generally illustrated as a continuous-time Markov chain, which drives the changes of the main parameters of population models with state switchings of homogeneous Markov chain [6,9,10,22,31], and the systems are described by the regime-switching differential equations, i.e., the piecewise-deterministic Markov process. To discuss the impact of environmental noise on the epidemic system in the host population, it is quite important to investigate the effect of random switching of environmental regimes on the spread dynamics of the disease.
In all these work mentioned above, it is assumed that the regime-switching of external environment factors follows a homogeneous continuous-time Markov chain. Li et al. [12] and Swishchuk et al. [48] claimed that such simple continuoustime Markov chain has many limitations in practical applications, because for a continuous-time Markov chain, the sojourn time in each environmental regime is necessarily exponentially distributed. In [15,16], Mollison and Vergu proved that the assumption of exponentially distributed sojourn times in infectious states is an acceptable approximation, but it is generally not realistic. For example, as shown in [13,14], it is more appropriate to model the dry spell length distribution by Pearson type III distribution, gamma distribution or Weibull distribution, where the dry spell is defined as the consecutive days with daily rain amount below some given threshold. In [12], Li et al. assumed that in each environmental state the conditional holding time distribution is a general distribution, i.e., the process used for describing environmental stochasticity is semi-Markov one, and they proposed an SIRS model whose infection coefficient is driven by semi-Markov switching. In addition to semi-Markov switching caused by changes in environmental factors, infectious diseases are also subject to environmental fluctuations, i.e., the white noise.
The main purpose of studying the dynamics of infectious diseases is to establish reasonable mathematical models, thereby providing effective strategies to control the spread of diseases. It has been verified that the sensitivity analysis of the basic reproduction number plays an important role in the formulation of disease prevention and control strategies [17,18,20,21]. As far as we know, in the current research on diseases transmission based on semi-Markov switching, the sensitivity analysis of the basic reproduction number is rarely mentioned. In this paper, we will provide a method to study the sensitivity analysis of the basic reproduction number for the stochastic epidemic system driven by semi-Markov switching.
In this paper, we will formulate a stochastic SIRS model with nonlinear incidence and vaccination, and the disease transmission rate is driven by a semi-Markov process and under the effect of white noise. The main contributions of this paper are: 1) We incorporate the nonlinear incidence and vaccination into the model proposed by Li et al. [12], where the general non-linear effective incidence rate can contain almost all of the cases raised in the known literature. Moreover, we also assume that the disease transmission rate is also affected by white noise.
2) We derive a new basic reproduction number R s 0 for the generalized stochastic epidemic model, and show that the disease extinction or persistence can be determined by the threshold number.
3) By giving the sensitivity analysis of the basic reproduction number under different semi-Markov chains, we demonstrate that the effects of different factors on the basic reproduction number will vary with different semi-Markov regimeswitching environment. 4) We use a specific example to compare the effects of the different semi-Markov switchings on disease dynamics, and a detailed flowchart of the simulation of a semi-Markov chain is given in the Appendix. This paper is organized as follows. The basic knowledge about semi-Markov process and our epidemic model are given in Section 2. Section 3 establishes the basic reproduction number R s 0 of the stochastic model and gives threshold dynamics of the system. Section 4 gives examples under different semi-Markov switchings, and uses sensitivity analysis to discuss the essence that affects the dynamics of the entire system. Finally, we conclude our paper in Section 5. The flowchart of the simulation of semi-Markov chain is given in the Appendix.

Model and preliminaries.
2.1. The semi-Markov process. Let {r(t), t ≥ 0} be a semi-Markov process taking values in state space M = {1, ..., N }. It means that the states of external environments are depended on the weather conditions or other uncertainties. Then we define that: (p i,j ) M×M is the transition probability matrix, {F i (t), t ∈ [0, ∞)} is the conditional holding time distribution of {r(t), t ≥ 0}, {τ i , i = 0, 1...} are jump times, and σ i = τ i − τ i−1 are time intervals between adjacent two jumps. The transition probability from state i to state j is given as P (r(τ n+1 ) = j, σ n+1 ≤ t|r(τ n ) = i) = p i,j · F i (t) and the embedded chain {J n := r(τ n ), n = 0, 1...} of {r(t), t ≥ 0} is Markovian with one-step transition probability matrix (p i,j ) M×M . For each t ∈ [0, ∞) and i, j ∈ M, we define that We give the following same assumption (H) as that in [7,12], which will be valid in the whole paper: (2) For each i ∈ M, F i (·) has a continuous and bounded density f i (·), and f i (t) > 0 for all t ∈ [0, ∞); (3) For each i ∈ M, there exists a constant i > 0 such that The constraint conditions of the assumption (H) are essentially very weak, but in Definition 1 of [12], Li established a phase-type distribution (PH-distribution) of a nonnegative random variable, and in Proposition 2 they proved that this PHdistribution can be chosen to approximate the original distribution in any accuracy.
Since assumption (H) holds true for PH-distribution, Li and q i,j is the transition rate of the continuous-time Markov chain {r(t), t ≥ 0}. It follows from Section 2 in [22] that the Markov chain r(t) is generated by the transition rate matrix Q := (q i,j ) M×M , i.e., where ∆t > 0 is a small time interval, q i,j ≥ 0 if i = j and Then {θ(t), t > 0} is a process associated to the semi-Markov process {r(t), t ≥ 0}, and the process {θ(t), t > 0} represents the amount of time that the process {r(t), t ≥ 0} is at the current state after the last jump.
In [23], the process {θ(t), t ≥ 0} is also known as the backward recurrence time process or "age". The introduction of the backward recurrence time process is justified by the fact that the pair {(θ(t), r(t)), t ≥ 0} satisfies the Markov property.
2.2. Stochastics SIRS epidemic model. In this subsection, we will introduce our model and the definition of each parameter. SIRS model is more often and properly used to describe the transmission dynamics of infectious diseases, such as the Respiratory Syncytial Virus and influenza A. In [11], Li et al. presented the following SIRS model with varying population size:    dS(t) = (Λ − µS(t) + λR(t) − β(r(t))S(t)G(I(t))) dt, dI(t) = (β(r(t))S(t)G(I(t)) − (µ + α + δ)I(t)) dt, dR(t) = (δI(t) − (µ + λ)R(t)) dt. (1) In system (1), Λ is the recruitment rate of the susceptible; µ is the natural death rate of the population; δ is the recovery rate of the infected individuals, and α is the disease-induced death rate of the infected individuals. Furthermore, in system (1) the transmission rate β is driven by semi-Markov switching r(t) which is induced by the environmental factors. However, the system (1) did not include the impact of vaccination and white noise disturbance, and in the epidemiology, vaccination is an important strategy for the elimination of infectious diseases [27]. The studies of vaccination, treatment, and associated behavioral changes related to disease transmission have become the subject of intense theoretical analysis [5,18,27,29,31]. Thus, in the paper we will include the vaccination into system (1), and assume that a fraction p(p ∈ [0, 1]), of newborns susceptible individuals are vaccinated [28,30]. As mentioned in the introduction, it is necessary to include the white noise into system (1), and in this paper we wil use the expression β i + σḂ(t) to replace the expression β i , wherė B(t) represents the white noise, i.e., B(t) is the Brownian motion, and σ is the intensity of the noise. In addition, in system (1) Based on system (1) and the above discussion, the model we will investigate in the paper can be described as: where ψ(S, I) is an increasing function with respect to S and I. For infectious disease systems, there always exists ψ(0, 0) = 1(more details can be found in the paper [33,34,35,36,37,38]). The definitions of all variables and parameters are summarized in Table 1.
Remark 3. The infection coefficient β(·) in the incidence rate β(r(t))SI ψ(S,I) is driven by semi-Markov switching r(t). The positive function ψ(S, I) measures the inhibition effect from the behavioral change of the susceptible individuals S when their number increases or from the crowding effect of the infective individuals I, and it is assumed to satisfy the following assumption (H1): is locally Lipschitz continuous. Throughout this paper, let (Ω, F, P) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions, and denote ω as an element of Ω, the space R k k}. For any initial value Z(0) = (S(0), I(0), R(0)) ∈ R 3,o + , Z(t, ω, Z(0)) = (S(t, ω, Z(0)), I(t, ω, Z(0))), R(t, ω, Z(0)))) is denoted as the corresponding solution of system (2).  (2) Natation Biological meanings S(t) Number of susceptibles at time t The recruitment rate of the population p The proportion of population that is vaccinated µ The death rates of susceptibles, infectives, and recovered individuals β(·) The infection coefficient α The death rate of infected individuals from disease-related causes δ The recovery rate of the infective individuals λ The recovered individuals immunity lose rate The proof of Lemma 2.1 is similar to Lemma 1 in [10] and Theorem 2.1 in [45], hence we omit it here. For convenience and simplicity, in the following we always assume that (S 0 , I 0 , R 0 ) ∈ R 3 + .
3. Threshold dynamics. In this section, the threshold condition for extinction and persistence of the disease will be derived. In order to derive the expression for the basic reproduction number, we first give the definition of mean sojourn time of {r(t), t > 0} in state i: According to the definitions of semi-Markov process in Section 2, we further assume the following assumption (H2) holds. Assumption (H2): Under Assumption (H2), we can easily see that the semi-Markov r(t) process is ergodic by Section 2 in [42]. Set J n be an embedded Markov chain related to the semi-Markov process {r(t), t > 0}. It then follows from the ergodic theorem (Theorem 2, p. 244 of [8]) that we have for any bounded measureable function f , where v = {v k ; k ∈ M} is the unique stationary positive probability distribution of J n , i.e., the invariant measure for the transition probability matrix P = (p i,j ) (M×M) with vP = v. Now we are able to define the expression for the basic reproduction number. Denote where {β k , k ∈ M} is the infection coefficient under different regimes. In the following, we will prove the disease persistence and extinction can be determined by R s 0 .
3.1. Stochastic disease-free dynamics. In this subsection, we will investigate the disease-free dynamics under each one of the following conditions.
. The main theorem in the subsection can be stated as: Proof. Using the Itô formula, we have )dB(t), . In order to prove that the assertion (7) holds, let us consider the auxiliary processR(t) which is defined by the following differential equation Since I(t) ≥ 0, it follows from the comparison theorem that By Lemma 2.1, we have lim t→+∞ [S(t) + In order to prove (6), integrating both sides of (8) yields that Then M is a square-integrable continuous martingale. Since ψ(S(t), I(t)) > 1 and (10), the quadratic variation of M can be expressed as By the strong law of large numbers for a real-valued continuous local martingales (Theorem 3.4, p.12 of [43]), we have lim sup Recalling Lemma 2.1, (13) and taking limitation on each side of (11), we have )ds a.s.
Since R s 0 < 1, i.e., By the ergodic property of r(t), we have
Since lim t→∞ I(t) = 0 a.s., it follows from the third equation in (2) that Since Since σ > Substituting (22) into (14) yields s. Then using the similar arguments as that in Case 1, we can also get the required assertion (5), (6) and (7). This completes the proof of Theorem 3.1. when the disease dies out, the population density of susceptible is positively related to the recruitment rate of the population Λ, while the population density of recovered individuals is positively related to the vaccination proportion p. Therefore, vaccination increases the immunity of the population and inhibits the occurrence of the infection.

3.2.
Persistence of the disease. In this subsection, we will discuss the persistence of the disease, which corresponds to the stochastic strong persistence in the mean defined by Wang [40] and Zhao et al. [41]. First we give a lemma which will be used in the proof of Theorem 3.3.
For convenience, set S(t) = S, I(t) = I, ψ(S(t), I(t)) = ψ. It follows from Lemma 2.1 and Remark 3 that Considering the inequality (24) and (27), we have Then direct calculation yields that Since Φ k (x) is an unary quadratic function, it takes its maximum at x max = β k σ 2 . In order to discuss the extreme value, definẽ Since (24) is always satisfied, it is easy to see that In the following let us consider two cases: It then follows that Since Φ k (x) is an unary quadratic function and it takes its maximum at x max = β k σ 2 . Substituting the above inequality into (30), we have Case (2) Ifσ > σ, it then follows from (24) and (31) that By using (31), we can easily see that Direct computation yileds that Substituting (33) and (34) into the above expression, it then follows that

THRESHOLD DYNAMICS AND SENSITIVITY ANALYSIS OF A SIRS EPIDEMIC MODEL 13
Summarizing the above discussion, we can conclude that for σ > 0. Substituting (37) into (26) yields Then it follows from the equation for S(t) in system (2) that we have By the hypotheses (H1) we have ψ(0, 0)= 1, which implies, Define the function h(I) = 1 − 1 ψ(·,I) for ψ(S, I). Since ∂ψ(S,I) ∂I ≥ 0 and ψ(S, I) is locally Lipschitz continuous, we can conclude that there exists a constant θ > 0 such that Substituting (40) into (39), and recalling Lemma 2.1 and Remark 3, we obtain (42) From (38) and (42), it follows that Set Then we can see that M 1 and M 2 are both square-integrable continuous martingales. Since ψ > 1 and lim t→∞ [S(t) + I(t) + R(t)] < Λ µ , the quadratic variactions of M 1 and M 2 are given by . By Lemma 2.1 and the strong law of large numbers for local martingales, we have For the ergodic result of the semi-Markov process {r(t), t > 0}, it implies that Therefore, for Integrating the third equation of system (2) and dividing both sides by t, we get Combining (45) and (46), together with Lemma 2.1, we get lim inf This completes the proof of Theorem 3.3.
Theorem 3.3 indicates that when the intensity of noise σ is small enough and the threshold R s 0 > 1, then the disease will strong persistent in the mean, i.e., lim inf In this case, the long-term average density of the infected is positively related to the minimum infection coefficient of each subsystemβ k . By the definitions of ψ(S, I) in Remark 3 and Lemma 3.2, if we set p = 0 (no vaccination), the disease will break out at the initial time with any initial values since ψ(S, I) > 1 always holds. 1 1−p . Moreover, when the disease persists in the mean, the population density of recovered individuals R(t) is positively related to the p. Therefore, if the model parameters satisfy all conditions in Theorem 3.3, the vaccination cannot completely eliminate the disease. But the increase in vaccination proportion p still enable more population to gain immunity and prolongs the time before disease outbreak. Next section we will give some numerical examples to verify our theorems and illustrate the sensitivity analysis of the basic reproduction number (4).
By using the threshold of [45] and our R s 0 in (49), we will analyze the disease status of each subsystem (48) and the dynamic properties of system (47) under different semi-Markov switching, respectively. Here we choose the parameters in Table 2, which are used to investigate the dynamics of Pasteurella muris in colonies of laboratory mice [25].
By using the parameter in Table 2 and the basic reproduction number in (49) , we get R s 0 = 1.3181 > 1, hence the disease will be almost surely persistent in the time mean sense (Fig.(1) (a)). However, for each subsystem, by using the threshold from [45], we can obtain that in state 1, R 2 = 3.153 > 1, so there exists stationary distribution of solution (S 1 , I 1 , R 1 ) (Fig.(1) (b)); in state 2, since R 1 = −0.001 < 0, we get that I 2 (t) tends to zero exponentially almost surely (Fig.(1) (c)). Therefore, we come to the conclusion that the infection coefficient β driven by semi-Markov switching enhances the spread of disease, but this is only in the case of Hyperexponential distribution, next we will give a counterexample.   Table 2 with initial values (5, 1, 0), the distribution function of each state of semi-Markov chain as the Gamma distribution.
In this counterexample, since we still use the parameter values in Table 2, the dynamics of the model (48) (Fig.(2) (a,b,c)).
From those two examples, we obtain a completely contradictory result of the entire system, although each subsystem is identical. In order to analyze the direct factors affecting the existence and extinction of the disease in the entire system, we give the sensitivity analysis of (49) under different distribution function in Fig.(3) by using the partial rank correlation coefficients (PRCC) method from [46]. From Fig.(3), equations (3) and (49), it is obtained that the basic reproduction number R s 0 of model (47) is completely determined by the mean sojourn times m 1 and m 2 of the semi-Markov process {r(t), t ≥ 0}, namely, the corresponding distribution F i (t) of jump time between each state change the dynamics of the entire system. By fixing other parameters, we can establish a function: {R s 0 (m 1 , m 2 ) : ∂R s 0 ∂m2 < 0} and give the following analysis in Fig.(4). From Fig.(1), Fig.(2) and Fig.(4), our simulation results prove the following conclusion: (1) If all the subsystems share the same dynamics under the same state, the different distribution functions of jump time between subsystem will still induce changes in the dynamics of the entire system. (2) The dynamics of the entire system has nothing to do with its switching law but has a positive correlation to its mean sojourn time in each subsystem. Namely, the longer the system stays in a stable subsystem, the closer the entire system becomes to stability, and vice versa.

5.
Conclusion. In this paper, we provide a theoretical framework based on the stochastic SIRS model to study how semi-Markov switching and White noise affect the dynamics of disease transmission. As important potential factors that may affect the spread of the disease, nonlinear incidence and vaccination strategies are also included in the proposed model. By using the ergodic theory of semi-Markov switching and related techniques of stochastic analysis, we provide a threshold for judging whether the disease is extinct in the sense of probability. Comparing to the classical stochastic SIRS epidemic model driven by Markov-regime switching, our model has the following improvements: (1) The transmission rate of our model is driven by semi-Markov switching, which includes the case of Markov switching and can be used to describe more random changes in the environment. The model also considers the effects of parametric disturbances on the spread of the diseases.
(2) This model contains a nonlinear incidence function which is associated with both susceptible and infected individuals. This nonlinear incidence function is more general and can be applied to a series of models.
(3) In Example we give the sensitivity analysis of our R s 0 under two semi-Markov switchings, the most important conclusion we obtain is: if we keep all parametric values unchanged, the basic reproduction number R s 0 of model (2) is completely determined by the mean sojourn times m i of the semi-Markov process {r(t), t ≥ 0}.
Some interesting topics deserve further investigation. On the one hand, one can study the unique ergodic stationary distribution of system (2), on the other hand, one can introduce the influence of colored noise such as Lévy jump on the stochastic system. The sensitivity analysis of R s 0 is a guide and a prior for optimal control. Our simulation results prove that when studying the optimal control theory of different diseases, the most effective way to suppress the spread of disease is to prolong its time to extinction. One can further obtain the results of optimal control based on our results. We leave these for further investigations.
Appendix: The Simulation of Semi-Markov chain. The Markov chain is a special type of semi-Markov chain, and its numerical simulation is relatively simple. In the existing literature, few scholars have clearly defined the simulation process of the semi-Markov chain. So in this appendix we will give a simple and understandable semi-Markov chain simulation method by two steps as follows: (1) Simulate a "continuous" time chain The semi-Markov process is a continuous process, but in numerical simulation we can only use discrete processes to approximate this continuous process. First we give the following symbols and their definitions: Table 3. Symbols used in simulation

Symbols
Definition Value T Total simulation time 10 ti(i = 0, 1, · · · ) Every checkpoint {0, 1, · · · , 10} ts Accumulation of time of per cycle --∆t Time interval between each jump 1 si(i = 1, 2) The states of the semi-Markov switching {1,2} L(t) The state of the semi-Markov chain in t 1 or 2 r(t) The state of the discrete chain in t 1 or 2 M The transition probability matrix -- The most difference between Markov chain and semi-Markov is that semi-Markov chain is not memoryless, which means that the state of semi-Markov chain at time t i will affect its state at time t j . This memorable expression in the simulation is that the transition probability matrix of semi-Markov at each moment is constantly changing. For example, since our chain has two states, we can assume the transition probability matrix like: where P i,j (t) = 1 · F i (t) means the probability of transition from state i to j at time t, F i is the distribution function of state i and lim t→∞ F i = 1. Hence we can give the step below: Start L(t) = s0,t = 0  Figure 5. In this example we set F i (t) as the Hyper-exponential distribution. In order to make the picture clear, we adopted ∆t = T 10 . One can use the smaller ∆t to ensure accuracy.
(2) Construct the discrete semi-Markov chain After obtaining a continuous-time chain, we still have to proceed to the next step, namely, select a discrete semi-Markov chain from this continuous-time chain to simulate it to ensure higher accuracy.  Figure 6. We use a new interval ∆ẗ = ∆t/2, and assume r(t) = L(t i ) if t ∈ [t i ,t i+1 ]. Obviously, the smaller ∆t we choose, the more accurately this discrete semi-Markov chain is simulated.