# A Class of Deterministic and Stochastic Fractional Epidemic Models with Vaccination

Authors: Tingting Xue; Xiaolin Fan; Jian Zhu
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

In this paper, a class of fractional deterministic and stochastic susceptible-infected-removed- susceptible (SIRS) epidemic models with vaccination is proposed. For the fractional deterministic SIRS epidemic model, the existence of solution and the stability of equilibrium points are analyzed by using dynamic method. Then, the appropriate controls are established to effectively control the disease and eliminate it. On this basis, the fractional stochastic SIRS epidemic model with vaccination is further considered, and a numerical approximation method is proposed. The correctness of the conclusion is verified by numerical simulation.


---

## Body

## 1. Introduction

The integration of the global economy has brought about the accelerated development of international trade, frequent contact and communication between people, the aggravation of environmental pollution, and increased resistance of pathogens and vectors. As a result, many infectious diseases (such as venereal diseases, tuberculosis, schistosomiasis, and dengue fever) that had previously been extinct and controlled have reemerged and are spreading. Some newly emerging infectious diseases also come fierce, endangering human health. Therefore, it is of great significance to study epidemic dynamics model. The study of mathematical models of the spread of infectious diseases began with En’ko, and the work as a cornerstone was that of Kermark and Mekendrick in 1927. They divided the total population into three categories: susceptible (S), infected (I), and recovered (R). The SIR infectious disease model was established by using kinetic method. The spreading law and prevalence trend were also studied. In the past 20 years, the research on the dynamics of infectious diseases has developed rapidly in the world, and a large number of mathematical models have been used to analyze various infectious diseases. Most of these mathematical models are applicable to the study of general laws of various infectious diseases [1–6]. There are also some models for many specific diseases, such as measles, malaria, tuberculosis, influenza, smallpox, gonorrhea, and AIDS. A lot of achievements have been made in qualitative and quantitative analysis of infectious diseases through mathematical models, mainly focusing on the judgment and prediction of the development trend of the disease. In recent years, people are not only limited to the research and prediction of infectious diseases but also began to pay attention to the control of infectious diseases.Recently, the application of fractional differential and integral operators in mathematical models has become increasingly popular. In addition, the fractional differential equation has become the focus of many scholars [7, 8]. It frequently appears in a variety of mathematical applications, such as fluid mechanics, viscoelastic mechanics, economics, biology, physics, and engineering. With the development of fractional differential equation, the fractional model is closer to reality. Through the study of fractional model, the biological system can be studied in more detail and in depth. In particular, fractional differential equations themselves have memory and genetic properties, and memory is an important feature of immune response. Therefore, it is more reasonable to introduce fractional differential equation into infectious disease model. For example, the authors [9] proposed several fractional epidemic models. At the same time, a parameter estimation method was proposed to fit real data with fractional model. In [10], the authors studied the existence and uniqueness of solutions of a fractional SIRS model for HRSV diseases using fixed point theory. The fractional SIRS model was simulated and compared with the real experimental data in Florida. The results showed that, compared with classical models, fractional modeling and optimal control methods have lower cost and more effective effect in reducing the number of HRSV infected individuals. Article [11] proposed an operation matrix based on Bernstein wavelet to solve the fractional SIR model with unknown parameters. The fractional-order problems were transformed into algebraic equations by combining operation matrices with configuration method. This paper also discussed the Adams-Bashforth-Moulton prediction correction scheme to solve this problem. In [12], a mathematical model of hepatitis B virus with nonsingular nucleus fractional derivative of Caputo-Fabrizio was established. An iterative method for solving the model was proposed in this paper. The fixed point theorem was used to verify the uniqueness of the model solutions. For convenience, we give the definition of fractional derivative. Let f∈H1a,b, b>a, and 0<α<1. A fractional derivative with Caputo-Fabrizio version is defined as
(1)CFDtαft=∫0tf′xexp−αt−x/1−α1−αdx.A fractional derivative with Atangana-Baleanu version is defined by(2)ABCDtαft=Bα1−α∫0tf′xEα−αt−xα1−αdx,ABRDtαft=Bα1−αddt∫0tfxEα−αt−xα1−αdx.On the other hand, there are all kinds of unpredictability and randomness in the actual spread of an epidemic. In the real ecosystem, environmental disturbance is everywhere. The process of epidemic transmission is usually associated with changes in environmental factors, such as temperature, humidity, and other climatic characteristics. Environmental factors are understood as any external factors that influence epidemic model parameters. Therefore, it is of great significance to study the influence of random factors on infectious disease model. Stochastic differential equation model is a method to study this problem. In recent years, many scholars have introduced parameter perturbation into epidemic models and studied its dynamics. For example, Gray et al. [13] extended the classical infectious disease model from a deterministic framework to a stochastic model. The unique global positive solution of stochastic model was studied. Some conditions were then established for disease extinction and persistence. Perturbation was discussed by means of random noise. In the case of persistence, the author proved the existence of stationary distribution and obtained its mean and variance. The article [14] studied the dynamics of a stochastic SIRS infectious disease model with saturated incidence. When the noise was low, the authors obtained a threshold for the stochastic system that determines the extinction and persistence of infectious diseases. In addition, the authors found that loud noise suppressed the spread of the epidemic. In [15], by introducing random fluctuations, the authors extended the classical SIRS infectious disease model with intervention strategies from a deterministic framework to a stochastic differential equation model. Using Markov semigroup theory, the authors presented conditions for disease extinction and stochastical persistence. The authors found that random fluctuations are important in keeping disease outbreaks in check. In [16], the SEIR model of COVID-19 was established. The authors used particle swarm optimization to estimate the system parameters based on the measured data from Hubei Province. They found that the parameters of the SEIR model were different in different scenarios. Considering the seasonality and randomness of the parameters, the authors proposed control strategies for COVID-19 based on the structure and parameters of the model. In [17–19], a stochastic SIS infectious disease model with vaccination was studied. Applying Markov semigroup theory, the authors presented a random regeneration number, which was regarded as a threshold parameter for identifying random extinction and persistence of diseases. The authors found that loud ambient noise suppressed disease outbreaks.Inspired by the above literature, we will study a class of fractional deterministic and stochastic infectious disease models with vaccination. According to literature review, it is rare to introduce fractional derivatives and random perturbations into infectious disease models with vaccination. Therefore, the work of this paper fills this gap. Specifically, this paper extends the deterministic model of [20] to the SIRS infectious disease deterministic model with the Caputo type fractional derivative and the stochastic fractional SIRS infectious disease model and studies their dynamic behavior and control measures, respectively. Different from other articles, this paper also studies the Mittag-Leffler stability at the equilibrium point, and the stability at the equilibrium point of fractional stochastic infectious disease model. This part of work is relatively novel. The paper is arranged as follows: In Section 2, the fractional SIRS epidemic mathematical model is described. Some preliminary results, such as existence, uniqueness, non-negativity, boundedness, local stability, global stability, and Mittag-Leffler stability of equilibrium points, are presented. In Section 3, the control is exerted on the infected person, and the Lyapunov method is used to design the controller, so that the disease can be eliminated. In Section 4, the deterministic fractional SIRS model is extended to the stochastic fractional SIRS model. The stability result of the stochastic fractional SIRS model at the equilibrium point is given. The numerical approximation method for fractional stochastic SIRS epidemic model is proposed. The correctness of the conclusion is verified by numerical simulation in each section. Finally, the conclusions are given in Section 5.


## 2. Dynamical Analysis for Fractional SIRS Epidemic Model with Vaccination

The paper is interested in studying a model of infectious disease with vaccination. This kind of model comes from the article of Zhen and Ma [20]. For ease of reading, we first introduce the SIRS epidemic model with vaccination as follows:
(3)dSdt=bN−λSIN+θI+eR−μ+pS,dIdt=λSIN−μ+ε+c+θI,dRdt=cI+pS−μ+eR,with initial conditions
(4)S0=S0≥0,I0=I0≥0,R0=R0≥0,where S, I, and R, respectively, represent the number of susceptible, infected, and removed individuals at time t, S+I+R=N represents the number of total population at time t, b is the natural birth rate, λ is the transmission rate of the disease, θ is the transfer rate from I to S, c represents the treatment rate, e represents the loss of immunity rate, μ represents the natural death rate, ε represents the disease-related death rate, p represents vaccination rate. New positive variables are defined in this paper
(5)Λ=bN,a1=λN,a2=μ+p,a3=μ+ε+c+θ,a4=μ+e.Then, the model (3) is transformed into the following form:
(6)dSdt=Λ−a1SI+θI+eR−a2S,dIdt=a1SI−a3I,dRdt=cI+pS−a4R.Next, we express derivatives on left-hand-side of model (6) by Caputo fractional derivatives. Therefore, we write the fractional SIRS model as
(7)CDtαS=Λ−a1SI+θI+eR−a2S,CDtαI=a1SI−a3I,CDtαR=cI+pS−a4R,where 0<α≤1.


### 2.1. Preliminary Knowledge

For convenience, the relevant definitions and lemmas of fractional calculus in this paper are given below.Definition 1 (see [21]).
The Riemann-Liouville fractional integral of orderαα>0 of function f:0,∞⟶ℝ is defined as:
(8)Itαft=∫0tt−xα−1fxΓαdx.Definition 2 (see [21]).
The Caputo fractional derivative of functionf:0,∞⟶ℝ is defined as:
(9)CDtαft=∫0tt−x−αΓ1−αdfxdxdx,0<α≤1.Lemma 3 (see [22]).
Suppose that the vector functionft,X:ℝ+×ℝ3⟶ℝ3 satisfies the following conditions:
(1)
The functionft,X is Lebesgue measurable with respect to t, t∈ℝ+(2)
The functionft,X is continuous on ℝ3 with respect to X(3)
∂ft,X/∂X is continuous on ℝ3 with respect to X(4)
ft,X≤ω1+ω2X, ∀t∈ℝ+, and X∈ℝ3, where ω1 and ω2 are two positive constants
Then, the initial value problems(10)CDtαX=ft,X,0<α≤1,Xt0=X0,has a unique solution.Lemma 4 (see [22]).
Supposeft∈Ca,b and CDtαft∈Ca,b, 0<α≤1. If CDtαft≥0, ∀t∈a,b, then for ∀t∈a,b, ft is a nondecreasing function; if CDtαft≤0, ∀t∈a,b, then for ∀t∈a,b, ft is a nonincreasing function.Lemma 5 (see [23]).
Letxt∈ℝ be a continuous and derivable function. Then, for any t≥0, one has
(11)12DCtαx2t≤xtCDtαxt,∀α∈0,1.Lemma 6 (see [24]).
Letxt∈ℝ+ be a continuous and derivable function. Then, for any t≥0, one has
(12)CDtαxt−x∗−x∗lnxtx∗≤1−x∗xtCDtαxt,∀α∈0,1,x∗∈ℝ+.Lemma 7 (see [25]).
Supposef∈C0,∞,ℝ+,F∈C0,∞,ℝ, there exist χ,χ0 and T such that for ∀t≥T, the following inequality
(13)lnft≤χt−χ0∫0tχfsds+Fta.s.,and limt⟶∞Ft/t=0a.s. are valid; then, the following inequality holds:
(14)limt⟶∞1t∫0tfsds≤χχ0a.s..


### 2.2. Existence and Uniqueness, Nonnegativity, and Boundedness

Theorem 8.
Model (7) has a unique solution Xt=S,I,RT∈ℝ+3,t≥0 under the initial conditions (4). And all the solutions of model (7) are uniformly bounded and nonnegative.Proof.
Let us first prove that the model (7) has a unique solution. According to Lemma 3, the vector function of the model (7) is as follows
(15)ft,X=Λ−a1SI+θI+eR−a2Sa1SI−a3IcI+pS−a4R.Obviously, the vector functionf satisfies conditions (1)-(3) of Lemma 3. In order to prove the existence and uniqueness of the model solution, we only need to prove that the vector function f satisfies the condition (4) of Lemma 3. Let x1t=St,x2t=It,x3t=Rt,x10=S0=S0,x20=I0=I0,x30=R0=R0,Xt=x1tx2tx3t,A0=Λ00,A1=−a2θe0−a30pc−a4,A2=0−a10000000,andA3=000a100000. Then, the model (7) can be reduced to
(16)CDtαX=A0+A1Xt+x1tA2Xt+x2tA3Xt.Setft,Xt=A0+A1Xt+x1tA2Xt+x2tA3Xt. Then,
(17)ft,Xt=A0+A1Xt+x1tA2Xt+x2tA3Xt≤A0+A1X+hMA2X+A3X=ω1+ω2X,where hM=maxx1t,x2t,ω1=A0,ω2=A1+hMA2+A3,X=suptSt+suptIt+suptRt. Therefore, according to Lemma 3, the model (7) has a unique solution.Secondly, it is proved that all the solutions of fractional model (7) are uniformly bounded and nonnegative. By model (7) and Nt=St+It+Rt, we have
(18)CDtαNRt=Λ+Sp−a2+Iθ−a3+c+Re−a4=Λ+S−μ+I−μ−ε+R−μ=Λ−μNRt−εI.Thus,CDtαNRt+μNRt=Λ−εI≤Λ. According to Lemma 9 in [26], we obtain
(19)0≤NRt≤N0Eα−μtα+ΛtαEα,α+1−μtα,where Eα is the Mittag-Leffler function. By Lemma 5 and Corollary 6 in [26], we get 0≤NRt≤Λ/μ,t⟶∞. Thus, all the solutions of fractional model (7) that start in ℝ+3 are uniformly bounded in the region Π=S,I,R∈ℝ+3:NRt≤Λ/μ+ε,ε>0. From the first equation of model (7), one has
(20)CDtαS=Λ−a1SI+θI+eR−a2S≥−a1I+a2S≥−a1Λμ+a2S.By the fractional comparison theorem in [26], we have
(21)St≥S0Eα,1−a1Λμ+a2tα.BecauseS0≥0 and Eα,1>0, so St≥0. From the second equation of model (7), we have CDtαI=a1SI−a3I≥−a3I. Thus, It≥I0Eα,1−a3tα. Because I0≥0 and Eα,1>0, so It≥0. From the third equation of model (7), we have
(22)CDtαR=cI+pS−a4R≥−a4R.Thus,Rt≥R0Eα,1−a4tα. Because R0≥0 and Eα,1>0, so Rt≥0. Hence, all the solutions of fractional model (7) are nonnegative. That is, all solution of model (7) is in region D, where
(23)D=S,I,R∈ℝ+3:0<S+I+R≤Λμ.


### 2.3. Disease-Free Equilibrium Point, Basic Reproduction Number, and Stability

In order to study the equilibrium points of model (7), let CDtαS=0,CDtαI=0,CDtαR=0, then model (7) has a disease-free equilibrium point
(24)E0=S0,I0,R0=Λa4a2a4−ep,0,Λpa2a4−ep,where a2a4−ep=μμ+p+e>0. To calculate the basic reproduction number of model (7), we reorder the state variables of model (7) by setting X=S,I,R and rewrite in the matrix form
(25)CDtαX=FX−VX,where
(26)FX=Λ+θI+eRa1SIcI+pS,VX=a1SI+a2Sa3Ia4R.The Jacobian of Equation (26) around E0=S0,I0,R0. Then, the basic reproduction number R0 is the spectral radius of FV−1, which takes the form R0=a1S0/a3.Theorem 9.
IfR0<1, then the disease-free equilibrium point E0 of model (7) is locally asymptotically stable, while if R0>1, then E0 is unstable.Proof.
The Jacobian matrix of model (7) around E0 is
(27)JE0=−a2−a1S0+θe0a1S0−a30pc−a4.After calculation,JE0 has three eigenvalues: λ1=−μ, λ2=−μ+p+e, and λ3=a3R0−1. Obviously, λ1<0,λ2<0, and the sign of λ3 depends on R0. If and only if the R0<1, all of eigenvalues of JE0 are negative. Hence, model (7) is locally asymptotically stable around E0. If R0>1, then λ3 is positive and λ1,λ2<0, so E0 is an unstable.Theorem 10.
IfR0<1, then the disease-free equilibrium point E0 of model (7) is globally asymptotically stable in D.Proof.
Let us consider the following positive definite Lyapunov function(28)VS,I,R=12S−S0+I+R−R02+εa1I.By Lemma5, we obtain
(29)CDtαV≤S−S0+I+R−R0CDtαS+CDtαI+CDtαR+εa1CDtαI=S−S0+I+R−R0−μS−S0+I+R−R0−εI+εa1a1S−S0+S0I−a3I=−μS−S0+I+R−R02−εI2−εIR−R0+εa1Ia1S0−a3.SinceR0<1, so a1S0−a3<0. Hence, CDtαV<0. CDtαV=0 if and only if S=S0,I=0,R=R0. Set M=S,I,R∈DDCtαVS,I,R=0. When t⟶+∞, M⟶E0. So, E0 is the unique greatest positive invariant of M. According to generalized Lyapunov-Lasalle’s invariance principle [22], E0 is globally asymptotically stable in D if R0<1.Next, to prove the Mittag-Leffler stability of the solution, we first give the following definition. For more details, please refer to literature [27].Definition 11.
The trivial solution of the following initial value problem(30)CDtαX=ft,X,0<α≤1,Xt0=X0,is said to be Mittag-Leffler stable if the function ft,X:ℝ+×ℝ3⟶ℝ3 is continuous locally Lipschitz function, ft,0=0, and X≤mXt0Eα−λt−t0αυ, where t0 is the initial time, λ≥0, υ>0, m0=0,mX≥0, and mX is locally Lipschitz on X∈ℝ3 with the Lipschitz constant m0.Theorem 12.
LetX=S,I,R. If there exists a continuously differentiable function Vt,X:0,+∞×ℝ+3⟶ℝ, such that it is a locally Lipschitz with respect to X and a class K function ψ satisfying:
(31)k1Xk0≤Vt,X≤k2ψX,CDtαVt,X≤−k3ψX,where k0,k1,k2,k3>0, then the disease-free equilibrium point E0 of model (7) is Mittag-Leffler stable. If (31) holds globally on ℝ+3, then E0 is globally Mittag-Leffler stable.Proof.
By (31), one has CDtαVt,X≤−k3/k2Vt,X. According to Lemma 9 in [26], we have Vt,X≤V0,X0Eα−k3/k2tα. Combining the above formula with (31), we obtain
(32)k1Xk0≤V0,X0Eα−k3k2tα.
So, if we simplify this, we get(33)X≤V0,X0k1Eα−k3k2tα1/k0.WhenmX0=V0,X0/k1,λ=k3/k2,andυ=1/k0, by Definition 11, it follows that the disease-free equilibrium point E0 of model (7) is Mittag-Leffler stable.


### 2.4. Endemic Equilibrium Point and Stability

The endemic equilibrium is equilibrium of model (7) in which the infected component of model is nonzero. We obtain the endemic equilibrium point E∗=S∗,I∗,R∗, where
(34)S∗=a3a1,I∗=a3a2a4−ep1−R0ce+a4θ−a3a1,R∗=ca1I∗+pa3a1a4.Theorem 13.
IfR0>1, then there exists the endemic equilibrium point E∗; otherwise, it does not exist.Proof.
After simple calculation,a2a4−ep=μμ+p+e>0 and ce+a4θ−a3=−μ+eμ+ε−μc<0.
IfR0>1, then 1−R0<0. Thus, I∗>0, thereby R∗>0. Consequently, the endemic equilibrium point exists, if R0>1.The Jacobian matrix of model (7) around E∗ is
(35)JE∗=−a1I∗−a2−a3+θea1I∗00pc−a4.Through simple calculation, the characteristic equation ofJE∗ can be obtained as follows:
(36)Pλ=λ3+A1λ2+A2λ+A3=0,where A1=a1I∗+a2+a4>0,A2=a1I∗2μ+ε+c+e+μμ+p+e>0,(37)A3=a3μμ+p+eR0−1>0R0>1.Theorem 14.
If all the eigenvalues of Jacobian matrixJE∗ satisfy argλ>απ/2, then the endemic equilibrium point E∗ of model (7) is locally asymptotically stable.LetDP=18A1A2A3+A1A22−4A3A13−4A23−27A32. According to [28], we have the following conclusions.Theorem 15.
IfR0>1, and assuming that any of the following three conditions holds, then the endemic equilibrium point E∗ of model (7) is locally asymptotically stable:
(38)DP>0,A1A2>A3,DP<0,α<23,DP<0,A1A2=A3,α∈0,1.Theorem 16.
IfR0>1, then the disease-free equilibrium point E∗ of model (7) is globally asymptotically stable in D.Proof.
Let us consider the following positive definite Lyapunov function(39)VS,I,R=12S−S∗+I−I∗2+τ1I−I∗−I∗lnII∗+12τ2R−R∗2,where τ1=a2+a3−θ/a1=2μ+ε+c+p/a1>0,τ2>0. Then, the function is positive. By Lemmas 5 and 6, we obtain CDtαV≤S−S∗+I−I∗CDtαS+CDtαI+τ11−I∗/ICDtαI+τ2R−R∗CDtαR=θ−a3I−I∗+eR−R∗−a2S−S∗×S−S∗+I−I∗+τ1a1I−I∗S−S∗+τ2R−R∗cI−I∗+pS−S∗−a4R−R∗=−a2S−S∗2−a3−θI−I∗2−τ2a4R−R∗2+S−S∗I−I∗θ−a3−a2+τ1a1+e+τ2pS−S∗R−R∗+e+τ2cI−I∗R−R∗. Because τ1=a2+a3−θ/a1, so θ−a3−a2+τ1a1=0. Thus,
(40)CDtαV≤−a2S−S∗2−a3−θI−I∗2−τ2a4R−R∗2+e+τ2pS−S∗R−R∗+e+τ2cI−I∗R−R∗≤−a2S−S∗2−a3−θI−I∗2−τ2a4R−R∗2+3a24S−S∗2+e+τ2p23a2R−R∗2+3a3−θ4I−I∗2+e+τ2c23a3−θR−R∗2=−a24S−S∗2−a3−θ4I−I∗2−τ2a4−e+τ2p23a2−e+τ2c23a3−θR−R∗2,where a3−θ=μ+ε+c>0. Let τ2>0, such that τ2a4−e+τ2p2/3a2−e+τ2c2/3a3−θ≥0. Thus, CDtαV≤0. CDtαV=0 if and only if S=S∗, I=I∗, and R=R∗. Set M=S,I,R∈DDCtαVS,I,R=0. When t⟶+∞, M⟶E∗. So, E∗ is the unique greatest positive invariant of M. Hence, the disease-free equilibrium point E∗ of model (7) is globally asymptotically stable in D if R0>1.


### 2.5. Numerical Simulation

The following examples are used to verify the previous theoretical results.Example 1.
In model (7), let Λ=0.2, a1=0.4, a2=0.26, a3=0.5, a4=0.34, θ=0.05, e=0.14, c=0.15, p=0.06, S0=1, I0=0.4, and R0=0. In this case, after the calculation, we have E0=0.85,0,0.15 and R0=0.68<1. Then, according to Theorem 9, it can be concluded that E0 is locally asymptotically stable.Example 2.
The values of the parameters excepta1=0.8 are the same as those in Example 1. In this case, after calculation, we have E∗=0.625,0.1364,0.17, R0=1.36>1, A1=0.70912, A2≈0.1662, A3=0.0144, A1A2≈0.1178557>A3, and DP≈0.000839>0. By Theorem 15, it can be concluded that E∗ is locally asymptotically stable.To support our results, we present computer simulations shown in Figures1–6.Figure 1
Stability analysis ofSt with a1=0.4.Figure 2
Stability analysis ofIt with a1=0.4.Figure 3
Stability analysis ofRt with a1=0.4.Figure 4
Stability analysis ofSt with a1=0.8.Figure 5
Stability analysis ofIt with a1=0.8.Figure 6
Stability analysis ofRt with a1=0.8.Figures1–3 correspond to Example 1. By observing Figures 1–3, the system is locally asymptotically stable at the disease-free equilibrium point E0=0.85,0,0.15. Figures 4–6 correspond to Example 2. By observing Figures 4–6, the system is locally asymptotically stable at the endemic equilibrium point E∗=0.625,0.1364,0.17.


## 2.1. Preliminary Knowledge

For convenience, the relevant definitions and lemmas of fractional calculus in this paper are given below.Definition 1 (see [21]).
The Riemann-Liouville fractional integral of orderαα>0 of function f:0,∞⟶ℝ is defined as:
(8)Itαft=∫0tt−xα−1fxΓαdx.Definition 2 (see [21]).
The Caputo fractional derivative of functionf:0,∞⟶ℝ is defined as:
(9)CDtαft=∫0tt−x−αΓ1−αdfxdxdx,0<α≤1.Lemma 3 (see [22]).
Suppose that the vector functionft,X:ℝ+×ℝ3⟶ℝ3 satisfies the following conditions:
(1)
The functionft,X is Lebesgue measurable with respect to t, t∈ℝ+(2)
The functionft,X is continuous on ℝ3 with respect to X(3)
∂ft,X/∂X is continuous on ℝ3 with respect to X(4)
ft,X≤ω1+ω2X, ∀t∈ℝ+, and X∈ℝ3, where ω1 and ω2 are two positive constants
Then, the initial value problems(10)CDtαX=ft,X,0<α≤1,Xt0=X0,has a unique solution.Lemma 4 (see [22]).
Supposeft∈Ca,b and CDtαft∈Ca,b, 0<α≤1. If CDtαft≥0, ∀t∈a,b, then for ∀t∈a,b, ft is a nondecreasing function; if CDtαft≤0, ∀t∈a,b, then for ∀t∈a,b, ft is a nonincreasing function.Lemma 5 (see [23]).
Letxt∈ℝ be a continuous and derivable function. Then, for any t≥0, one has
(11)12DCtαx2t≤xtCDtαxt,∀α∈0,1.Lemma 6 (see [24]).
Letxt∈ℝ+ be a continuous and derivable function. Then, for any t≥0, one has
(12)CDtαxt−x∗−x∗lnxtx∗≤1−x∗xtCDtαxt,∀α∈0,1,x∗∈ℝ+.Lemma 7 (see [25]).
Supposef∈C0,∞,ℝ+,F∈C0,∞,ℝ, there exist χ,χ0 and T such that for ∀t≥T, the following inequality
(13)lnft≤χt−χ0∫0tχfsds+Fta.s.,and limt⟶∞Ft/t=0a.s. are valid; then, the following inequality holds:
(14)limt⟶∞1t∫0tfsds≤χχ0a.s..


## 2.2. Existence and Uniqueness, Nonnegativity, and Boundedness

Theorem 8.
Model (7) has a unique solution Xt=S,I,RT∈ℝ+3,t≥0 under the initial conditions (4). And all the solutions of model (7) are uniformly bounded and nonnegative.Proof.
Let us first prove that the model (7) has a unique solution. According to Lemma 3, the vector function of the model (7) is as follows
(15)ft,X=Λ−a1SI+θI+eR−a2Sa1SI−a3IcI+pS−a4R.Obviously, the vector functionf satisfies conditions (1)-(3) of Lemma 3. In order to prove the existence and uniqueness of the model solution, we only need to prove that the vector function f satisfies the condition (4) of Lemma 3. Let x1t=St,x2t=It,x3t=Rt,x10=S0=S0,x20=I0=I0,x30=R0=R0,Xt=x1tx2tx3t,A0=Λ00,A1=−a2θe0−a30pc−a4,A2=0−a10000000,andA3=000a100000. Then, the model (7) can be reduced to
(16)CDtαX=A0+A1Xt+x1tA2Xt+x2tA3Xt.Setft,Xt=A0+A1Xt+x1tA2Xt+x2tA3Xt. Then,
(17)ft,Xt=A0+A1Xt+x1tA2Xt+x2tA3Xt≤A0+A1X+hMA2X+A3X=ω1+ω2X,where hM=maxx1t,x2t,ω1=A0,ω2=A1+hMA2+A3,X=suptSt+suptIt+suptRt. Therefore, according to Lemma 3, the model (7) has a unique solution.Secondly, it is proved that all the solutions of fractional model (7) are uniformly bounded and nonnegative. By model (7) and Nt=St+It+Rt, we have
(18)CDtαNRt=Λ+Sp−a2+Iθ−a3+c+Re−a4=Λ+S−μ+I−μ−ε+R−μ=Λ−μNRt−εI.Thus,CDtαNRt+μNRt=Λ−εI≤Λ. According to Lemma 9 in [26], we obtain
(19)0≤NRt≤N0Eα−μtα+ΛtαEα,α+1−μtα,where Eα is the Mittag-Leffler function. By Lemma 5 and Corollary 6 in [26], we get 0≤NRt≤Λ/μ,t⟶∞. Thus, all the solutions of fractional model (7) that start in ℝ+3 are uniformly bounded in the region Π=S,I,R∈ℝ+3:NRt≤Λ/μ+ε,ε>0. From the first equation of model (7), one has
(20)CDtαS=Λ−a1SI+θI+eR−a2S≥−a1I+a2S≥−a1Λμ+a2S.By the fractional comparison theorem in [26], we have
(21)St≥S0Eα,1−a1Λμ+a2tα.BecauseS0≥0 and Eα,1>0, so St≥0. From the second equation of model (7), we have CDtαI=a1SI−a3I≥−a3I. Thus, It≥I0Eα,1−a3tα. Because I0≥0 and Eα,1>0, so It≥0. From the third equation of model (7), we have
(22)CDtαR=cI+pS−a4R≥−a4R.Thus,Rt≥R0Eα,1−a4tα. Because R0≥0 and Eα,1>0, so Rt≥0. Hence, all the solutions of fractional model (7) are nonnegative. That is, all solution of model (7) is in region D, where
(23)D=S,I,R∈ℝ+3:0<S+I+R≤Λμ.


## 2.3. Disease-Free Equilibrium Point, Basic Reproduction Number, and Stability

In order to study the equilibrium points of model (7), let CDtαS=0,CDtαI=0,CDtαR=0, then model (7) has a disease-free equilibrium point
(24)E0=S0,I0,R0=Λa4a2a4−ep,0,Λpa2a4−ep,where a2a4−ep=μμ+p+e>0. To calculate the basic reproduction number of model (7), we reorder the state variables of model (7) by setting X=S,I,R and rewrite in the matrix form
(25)CDtαX=FX−VX,where
(26)FX=Λ+θI+eRa1SIcI+pS,VX=a1SI+a2Sa3Ia4R.The Jacobian of Equation (26) around E0=S0,I0,R0. Then, the basic reproduction number R0 is the spectral radius of FV−1, which takes the form R0=a1S0/a3.Theorem 9.
IfR0<1, then the disease-free equilibrium point E0 of model (7) is locally asymptotically stable, while if R0>1, then E0 is unstable.Proof.
The Jacobian matrix of model (7) around E0 is
(27)JE0=−a2−a1S0+θe0a1S0−a30pc−a4.After calculation,JE0 has three eigenvalues: λ1=−μ, λ2=−μ+p+e, and λ3=a3R0−1. Obviously, λ1<0,λ2<0, and the sign of λ3 depends on R0. If and only if the R0<1, all of eigenvalues of JE0 are negative. Hence, model (7) is locally asymptotically stable around E0. If R0>1, then λ3 is positive and λ1,λ2<0, so E0 is an unstable.Theorem 10.
IfR0<1, then the disease-free equilibrium point E0 of model (7) is globally asymptotically stable in D.Proof.
Let us consider the following positive definite Lyapunov function(28)VS,I,R=12S−S0+I+R−R02+εa1I.By Lemma5, we obtain
(29)CDtαV≤S−S0+I+R−R0CDtαS+CDtαI+CDtαR+εa1CDtαI=S−S0+I+R−R0−μS−S0+I+R−R0−εI+εa1a1S−S0+S0I−a3I=−μS−S0+I+R−R02−εI2−εIR−R0+εa1Ia1S0−a3.SinceR0<1, so a1S0−a3<0. Hence, CDtαV<0. CDtαV=0 if and only if S=S0,I=0,R=R0. Set M=S,I,R∈DDCtαVS,I,R=0. When t⟶+∞, M⟶E0. So, E0 is the unique greatest positive invariant of M. According to generalized Lyapunov-Lasalle’s invariance principle [22], E0 is globally asymptotically stable in D if R0<1.Next, to prove the Mittag-Leffler stability of the solution, we first give the following definition. For more details, please refer to literature [27].Definition 11.
The trivial solution of the following initial value problem(30)CDtαX=ft,X,0<α≤1,Xt0=X0,is said to be Mittag-Leffler stable if the function ft,X:ℝ+×ℝ3⟶ℝ3 is continuous locally Lipschitz function, ft,0=0, and X≤mXt0Eα−λt−t0αυ, where t0 is the initial time, λ≥0, υ>0, m0=0,mX≥0, and mX is locally Lipschitz on X∈ℝ3 with the Lipschitz constant m0.Theorem 12.
LetX=S,I,R. If there exists a continuously differentiable function Vt,X:0,+∞×ℝ+3⟶ℝ, such that it is a locally Lipschitz with respect to X and a class K function ψ satisfying:
(31)k1Xk0≤Vt,X≤k2ψX,CDtαVt,X≤−k3ψX,where k0,k1,k2,k3>0, then the disease-free equilibrium point E0 of model (7) is Mittag-Leffler stable. If (31) holds globally on ℝ+3, then E0 is globally Mittag-Leffler stable.Proof.
By (31), one has CDtαVt,X≤−k3/k2Vt,X. According to Lemma 9 in [26], we have Vt,X≤V0,X0Eα−k3/k2tα. Combining the above formula with (31), we obtain
(32)k1Xk0≤V0,X0Eα−k3k2tα.
So, if we simplify this, we get(33)X≤V0,X0k1Eα−k3k2tα1/k0.WhenmX0=V0,X0/k1,λ=k3/k2,andυ=1/k0, by Definition 11, it follows that the disease-free equilibrium point E0 of model (7) is Mittag-Leffler stable.


## 2.4. Endemic Equilibrium Point and Stability

The endemic equilibrium is equilibrium of model (7) in which the infected component of model is nonzero. We obtain the endemic equilibrium point E∗=S∗,I∗,R∗, where
(34)S∗=a3a1,I∗=a3a2a4−ep1−R0ce+a4θ−a3a1,R∗=ca1I∗+pa3a1a4.Theorem 13.
IfR0>1, then there exists the endemic equilibrium point E∗; otherwise, it does not exist.Proof.
After simple calculation,a2a4−ep=μμ+p+e>0 and ce+a4θ−a3=−μ+eμ+ε−μc<0.
IfR0>1, then 1−R0<0. Thus, I∗>0, thereby R∗>0. Consequently, the endemic equilibrium point exists, if R0>1.The Jacobian matrix of model (7) around E∗ is
(35)JE∗=−a1I∗−a2−a3+θea1I∗00pc−a4.Through simple calculation, the characteristic equation ofJE∗ can be obtained as follows:
(36)Pλ=λ3+A1λ2+A2λ+A3=0,where A1=a1I∗+a2+a4>0,A2=a1I∗2μ+ε+c+e+μμ+p+e>0,(37)A3=a3μμ+p+eR0−1>0R0>1.Theorem 14.
If all the eigenvalues of Jacobian matrixJE∗ satisfy argλ>απ/2, then the endemic equilibrium point E∗ of model (7) is locally asymptotically stable.LetDP=18A1A2A3+A1A22−4A3A13−4A23−27A32. According to [28], we have the following conclusions.Theorem 15.
IfR0>1, and assuming that any of the following three conditions holds, then the endemic equilibrium point E∗ of model (7) is locally asymptotically stable:
(38)DP>0,A1A2>A3,DP<0,α<23,DP<0,A1A2=A3,α∈0,1.Theorem 16.
IfR0>1, then the disease-free equilibrium point E∗ of model (7) is globally asymptotically stable in D.Proof.
Let us consider the following positive definite Lyapunov function(39)VS,I,R=12S−S∗+I−I∗2+τ1I−I∗−I∗lnII∗+12τ2R−R∗2,where τ1=a2+a3−θ/a1=2μ+ε+c+p/a1>0,τ2>0. Then, the function is positive. By Lemmas 5 and 6, we obtain CDtαV≤S−S∗+I−I∗CDtαS+CDtαI+τ11−I∗/ICDtαI+τ2R−R∗CDtαR=θ−a3I−I∗+eR−R∗−a2S−S∗×S−S∗+I−I∗+τ1a1I−I∗S−S∗+τ2R−R∗cI−I∗+pS−S∗−a4R−R∗=−a2S−S∗2−a3−θI−I∗2−τ2a4R−R∗2+S−S∗I−I∗θ−a3−a2+τ1a1+e+τ2pS−S∗R−R∗+e+τ2cI−I∗R−R∗. Because τ1=a2+a3−θ/a1, so θ−a3−a2+τ1a1=0. Thus,
(40)CDtαV≤−a2S−S∗2−a3−θI−I∗2−τ2a4R−R∗2+e+τ2pS−S∗R−R∗+e+τ2cI−I∗R−R∗≤−a2S−S∗2−a3−θI−I∗2−τ2a4R−R∗2+3a24S−S∗2+e+τ2p23a2R−R∗2+3a3−θ4I−I∗2+e+τ2c23a3−θR−R∗2=−a24S−S∗2−a3−θ4I−I∗2−τ2a4−e+τ2p23a2−e+τ2c23a3−θR−R∗2,where a3−θ=μ+ε+c>0. Let τ2>0, such that τ2a4−e+τ2p2/3a2−e+τ2c2/3a3−θ≥0. Thus, CDtαV≤0. CDtαV=0 if and only if S=S∗, I=I∗, and R=R∗. Set M=S,I,R∈DDCtαVS,I,R=0. When t⟶+∞, M⟶E∗. So, E∗ is the unique greatest positive invariant of M. Hence, the disease-free equilibrium point E∗ of model (7) is globally asymptotically stable in D if R0>1.


## 2.5. Numerical Simulation

The following examples are used to verify the previous theoretical results.Example 1.
In model (7), let Λ=0.2, a1=0.4, a2=0.26, a3=0.5, a4=0.34, θ=0.05, e=0.14, c=0.15, p=0.06, S0=1, I0=0.4, and R0=0. In this case, after the calculation, we have E0=0.85,0,0.15 and R0=0.68<1. Then, according to Theorem 9, it can be concluded that E0 is locally asymptotically stable.Example 2.
The values of the parameters excepta1=0.8 are the same as those in Example 1. In this case, after calculation, we have E∗=0.625,0.1364,0.17, R0=1.36>1, A1=0.70912, A2≈0.1662, A3=0.0144, A1A2≈0.1178557>A3, and DP≈0.000839>0. By Theorem 15, it can be concluded that E∗ is locally asymptotically stable.To support our results, we present computer simulations shown in Figures1–6.Figure 1
Stability analysis ofSt with a1=0.4.Figure 2
Stability analysis ofIt with a1=0.4.Figure 3
Stability analysis ofRt with a1=0.4.Figure 4
Stability analysis ofSt with a1=0.8.Figure 5
Stability analysis ofIt with a1=0.8.Figure 6
Stability analysis ofRt with a1=0.8.Figures1–3 correspond to Example 1. By observing Figures 1–3, the system is locally asymptotically stable at the disease-free equilibrium point E0=0.85,0,0.15. Figures 4–6 correspond to Example 2. By observing Figures 4–6, the system is locally asymptotically stable at the endemic equilibrium point E∗=0.625,0.1364,0.17.


## 3. Apply Control over Model (7)

According to the previous discussion, whenR0<1, is globally asymptotically stable. That is, infectious diseases will be eliminated from the area. When R0>1, E0 is unstable, and E∗ is globally asymptotically stable and becomes endemic. In order to avoid the formation of E∗, we take measures to effectively control the disease, so that when R0>1, E0 can also be globally asymptotically stable, to eliminate the disease. To this end, we apply controller US,I so that model (7) becomes the following system:
(41)CDtαS=Λ−a1SI+θI+eR−a2S,CDtαI=a1SI−a3I+US,I,CDtαR=cI+pS−a4R.In medicine, two methods are commonly used to exert control over an infected person: one is effective medication, and the other is to isolate those who have the disease from those who are susceptible. Letk1 and k2 (k1>0,k2>0) represent isolation rate and cure rate, respectively, and consider the controller of the following form:
(42)US,I=−k1a1SI−k2I.Theorem 17.
IfUS,I=−k1a1SI−k2I and 0<k1<1,k2>0, then when R0>1, E0 is globally asymptotically stable in D.Proof.
Consider the Lyapunov functionV=I, then
(43)CDtαV=CDtαI=a1SI−a3I+US,I=a1SI−a3I−k1a1SI−k2I=1−k1a1SI−a3+k2I≤0.CDtαV=0 if and only if I=0. And E0 is the largest invariant set on I=0, so E0 is globally asymptotically stable.To support our results, we present numerical simulations shown in Figures7–12. In model (41), let Λ=0.2,a1=0.8,a2=0.26,a3=0.5,a4=0.34, θ=0.05,e=0.14, c=0.15,p=0.06,ε=0.1,μ=0.2,S0=1,I0=0.4,andR0=0. In this case, after the calculation, we have E0=0.85,0,0.15 and R0=1.36>1. Let k1 and k2 take different values that satisfy the conditions in Theorem 17. Figures 7–12 show the numerical simulation of different values of fractional orders. By observing the images, it can be seen that the higher the values of k1 and k2 are, the faster the rate of disease I tends to 0; thus, the disease is eliminated. That is, the controller selected in this part is effective.Figure 7
Infected population in case ofa1=0.8, a=1.Figure 8
Infected population in case ofa1=0.8, a=0.9.Figure 9
Infected population in case ofa1=0.8, a=0.85.Figure 10
Infected population in case ofa1=0.8, a=0.8.Figure 11
Infected population in case ofa1=0.8, a=0.7.Figure 12
Infected population in case ofa1=0.8, a=0.6.


## 4. Dynamical Analysis for Stochastic Fractional SIRS Epidemic Model with Vaccination

We extend SIRS model (7) to stochastic SIRS model as follows:
(44)CDtαS=Λ−a1SI+θI+eR−a2S+σ1Sdσ1dt,CDtαI=a1SI−a3I+σ2Idσ2dt,CDtαR=cI+pS−a4R+σ3Rdσ3dt,where σi0=0,i=1,2,3.σi>0 represents the intensity of white noise, which is an independent standard Brownian motion. The disease-free equilibrium obtained when It≡0 is the same point as deterministic model (7). The stability of stochastic model (44) can be studied by using the Lyapunov stability method suitable for stochastic model (44). Let us consider the following positive definite Lyapunov function VS,I,R=1/2S2t+1/2I2t+1/2R2t. By Lemma 5, we obtain
(45)CDtαV=12CDtαS2t+12CDtαI2t+12CDtαR2t≤StCDtαSt+ItCDtαIt+RtCDtαRt=StbS+I+R−a1SI+θI+eR−a2S+σ1Sdσ1dt+Ita1SI−a3I+σ2Idσ2dt+RtcI+pS−a4R+σ3Rdσ3dt=S2tb−a1I−a2+σ1dσ1dt+I2ta1S−a3+σ2dσ2dt+R2t−a4+σ3dσ3dt+StItb+θ+StRtb+e+p+cItRt≤S2tb−a1I−a2+σ1dσ1dt+I2ta1S−a3+σ2dσ2dt+R2t−a4+σ3dσ3dt+b+θ2S2t+I2t+b+e+p2S2t+R2t+c2I2t+R2t=S2t2b+θ+e+p2−a1I−a2+σ1dσ1dt+I2ta1S−a3+b+θ+c2+σ2dσ2dt+R2t−a4+b+e+p+c2+σ3dσ3dt.Let2b+θ+e+p/2−a1I−a2+σ1dσ1/dt<0,a1S−a3+b+θ+c/2+σ2dσ2/dt<0, −a4+b+e+p+c/2+σ3dσ3/dt<0. Then, CDtαV≤0. Hence, the stochastic SIRS model (44) is stable to E0.Fractional stochastic SIRS model (44) can be written as follows:
(46)CDtαXt=FX+GXdσdt,where X=S,I,RT,GX=σ1S,σ2I,σ3R,dσ/dt=dσ1/dt,dσ2/dt,dσ3/dtT and the mapping FX=F1X,F2X,F3XT with F1X=Λ−a1SI+θI+eR−a2S,F2X=a1SI−a3I,F3X=cI+pS−a4R. By Definition 1, we have
(47)Xt=X0+1Γα∫0tt−sα−1FXds+1Γα∫0tt−sα−1GXdσs.From the perspective of the numerical implementation, the fractional derivative definition of the Grunwald-Letnikov is the most direct. So, it is more appropriate to solve the fractional stochastic epidemic model. The Grunwald-Letnikov fractional derivativeGLDtα is defined as follows:
(48)GLDtαft=limh⟶0h−α∑j=0n−1jαjft−jh,where t=nh. For details, see [29–32]. Formula (48) can be simplified to
(49)GLDtαftn≈h−α∑j=0nwjαftn−j,where tn=nh,wjα satisfies recursive relation as follows:
(50)w0α=1,wjα=1−1+αjwj−1α,j=1,2,3,⋯.Ifft is a continuous integrable function, then the relation between Grunwald-Letnikov and Caputo fractional derivative is as follows:
(51)CDtαftn=GLDtαftn−f0tn−αΓ1−α≈h−α∑j=0nwjαftn−j−f0.For details, see [33–35]. Model (44) can be transformed into the following system through the approximate formula of Grunwald-Letnikov fractional derivative
(52)Sn=S0+hαΛ−a1Sn−1In−1+θIn−1+eRn−1−a2Sn−1+σ1Sn−1hξ1n−∑j=1nwjαSn−j−S0,In=I0+hαa1Sn−1In−1−a3In−1+σ2In−1hξ2n−∑j=1nwjαIn−j−I0,Rn=R0+hαcIn−1+pSn−1−a4Rn−1+σ3Rn−1hξ3n−∑j=1nwjαRn−j−R0,where σi are real constants, ξin represent the three-dimensional Gaussian white noise processes, i=1,2,3, and
(53)ξjt=0j=0,1,2,3,ξit1ξjt2=δijδt1−tj,where δij represents Kronecker delta and δt1−tj represents the Dirac function.In this part, we present numerical simulations of stochastic fractional SIRS epidemic model (44). Let Λ=0.2, a1=0.4, a2=0.26, a3=0.5, a4=0.34, θ=0.05, e=0.14, c=0.15, p=0.06, S0=1, I0=0.4, and R0=0. In this case, after the calculation, we have E0=0.85,0,0.15 and R0=0.68<1. Numerical simulations are performed for the different fractional order values. Figures 13–15 describe the following random density values of σ1=0.4, σ2=0.75, σ3=0.35. Figures 16–18 describe the stochastic behavior of St for different fractional derivative values 1, 0.8, and 0.6. Figures 19–21 describe the stochastic behavior of It for different fractional derivative values 1, 0.8, and 0.6. Figures 22–24 describe the stochastic behavior of Rt for different fractional derivative values 1, 0.8, and 0.6. Observing the images, we found that different random perturbations and different fractional derivative values affect the speed of the system to the equilibrium point.Figure 13
Stochastic behavior ofSt with σ1=0.4.Figure 14
Stochastic behavior ofIt with σ2=0.75.Figure 15
Stochastic behavior ofRt with σ3=0.35.Figure 16
Stochastic behavior ofSt with α=1.Figure 17
Stochastic behavior ofSt with σ2=0.8.Figure 18
Stochastic behavior ofSt with σ2=0.6.Figure 19
Stochastic behavior ofIt with α=1.Figure 20
Stochastic behavior ofIt with α=0.8.Figure 21
Stochastic behavior ofIt with α=0.6.Figure 22
Stochastic behavior ofRt with α=1.Figure 23
Stochastic behavior ofRt with α=0.8.Figure 24
Stochastic behavior ofRt with α=0.6.


## 5. Conclusion

This paper investigates a class of deterministic and stochastic fractional SIRS epidemic models with vaccination. The models studied in this paper are more general. Specifically, this paper extends the deterministic model of [20] to the SIRS infectious disease deterministic model with the Caputo type fractional derivative and the stochastic fractional SIRS infectious disease model and studies their dynamic behavior and control measures, respectively. Different from other articles, this paper also studies the Mittag-Leffler stability at the equilibrium point and the stability at the equilibrium point of fractional stochastic infectious disease model. This part of work is relatively novel. For the fractional deterministic SIRS model, the existence and uniqueness of solution, nonnegativity, and boundedness of solution, equilibrium points, and stability analysis (local stability, global stability, and Mittag-Leffler stability of equilibrium points) are given by using existing techniques. Then, the article takes measures to effectively control the disease and control the infected. By using Lyapunov method, a controller is designed to make the model’s disease-free equilibrium globally asymptotically stable when R0>1, so that the disease can be eliminated. On the other hand, by introducing noise into the disease transmission term, a fractional stochastic SIRS epidemic model with vaccination is further considered. The stability result of the stochastic fractional SIRS model at the equilibrium point is given. A numerical approximation method for fractional stochastic SIRS epidemic model is proposed. The correctness of the conclusion is verified by numerical simulation in each section. Our study shows that the fractional stochastic epidemic models based on virus dynamics are more realistic. This theory can provide a solid foundation for the study of similar diseases and has a wide range of applications in the biomedical field. For example, a stochastic delayed infectious disease model can be considered to study the effect of incubation periods on disease dynamics. In addition, our proposed theory can also be used to study other infectious diseases, such as HIV, COVID-19, and tuberculosis. We leave these problems to future work.



---

*Source: 1797258-2022-08-16.xml*