# On the Stochastic Dynamics of a Social Epidemics Model

Authors: Xun-Yang Wang; Peng-Zhan Zhang; Qing-Shan Yang
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

Alcohol abuse is a major social problem, which has caused a lot of damages or hidden dangers to the individual and the society. In this paper, with random factors of alcoholism considered in mortality rate of compartment populations, we formulate a stochastic alcoholism model according to compartment theory of infectious disease. Based on this model, we investigate the long-term stochastic dynamics behaviors of two equilibria of the corresponding deterministic model and point out the effect of random disturbance on the stability of the system. We find that whenR0≤1, we get the estimation between the trajectory of stochastic system and E0=(Π/μs,0,0,0) in the average in time with respect to the disturbance intensity, while when R0>1, stochastic system is ergodic and has the unique stationary distribution. Finally, we carry out numerical simulations to support the corresponding theoretical results.


---

## Body

## 1. Introduction

Along with the rapid development of human civilization, the pressure of human life is becoming more and more large, especially in the spirit and emotion. In order to release the pressure, many people will lose their rationality and take some unwise actions; for example, they will take addictive drugs or indulge in long-term alcohol abuse or heavy smoking [1]. These extreme behaviors not only seriously hurt the people’s own health but also cause serious security risks to the total society; for example, drunk driving and drunk dangerous sexual behaviors have frequently occurred [2–4]. Sociologists have noticed these undesirable negative phenomena, which are called social diseases with great destructive and infectivity due to mental or habitual addiction [5]. In the above social phenomena, alcoholism is particularly serious, showing the trend of lowering ages and femininity [6].Actually, similar to some typical infectious diseases, alcohol abuse also has a strong infectiousness, especially for the social crowd in frequent contact; therefore, it is often called social infectious diseases [5, 7]. In view of this, in the last ten years, many applied mathematicians and interdisciplinary workers are committed to formulating mathematical models to study the spreading behavior of alcohol abuse and predict the trend of alcoholism from a mathematical point of view [8–14].Sanchez et al. [8] earlier proposed a simple and basic mathematical model with relapse to describe the spread and infection of alcohol abuse. In their model, according to the progress of alcoholism, the population is divided into three compartments, that is, normal persons S(t) who do not drink or drink moderately; the alcoholics D(t), who will infect the people around to become new alcoholics; and the persons who temporarily quit drinking, denoted by R(t). The model in [8] is as follows:(1)S′=Λ-βSD-μS,D′=βSD+ρRD-μ+ϕD,R′=ϕD-ρRD-μR,based on which, the authors calculated the basic reproduction number of alcoholism R0 and proved the existence and stability of the two equilibria by virtue of R0 and, simultaneously, analyzed the influence of the parameters on the stability of the system. However, this model is somewhat rough and simple. On the basis of this model, many authors have built new mathematical models from different angles to study the behavior of alcoholism. Taking the different stages of alcohol abuse into account, [9, 10] proposed a two-stage alcoholism model, while [10] also investigated different infection rates according to the two stages. References [12] considered the warning function of information publicity and education on persons against alcohol abuse. Besides the awareness of information, [14] also considered the time delay between the contact and infection during the course of alcohol abuse. To highlight the negative effects of alcohol abuse, [13] proposed a rather complex model to investigate the effect of alcohol on HIV infection.Due to the fact that alcoholism can be prevented and treated [15], based on an earlier model of alcoholism [8], Wang et al. [16] proposed an SATR-type alcoholism model, and it is presented as follows:(2)S′=μN-1-u1tβSAN-μS,A′=1-u1tβSAN+ξT-u2t+μA,T′=u2tA-μ+ξ+δT,Q′=δT-μQ,in which the persons in treatment are considered as a population compartment separately and put forward effective measures to control alcohol abuse. Specifically, in model (2), the total population is partitioned into four compartments: S(t), which refers to the persons who never drink or drink moderately without affecting the physical health; A(t), which refers to the persons who drink heavily; T(t), which refers to the persons being treated by taking pills or other medical interventions after alcoholism; Q(t), which refers to the persons who recover from alcoholism after treatment and never drink again.Recently, Wang et al. [17] continue to consider model (2) with a distributed time delay between contact and infection during the course of alcoholism. Furthermore, in order to make the model have wider applicability in different environments, they generalize the incidence function from standard case to the abstract one.Actually, in one hand, alcoholism is a social behavior; therefore, it is only infected and transmitted in a limited population. To depict this attribute, we will adopt saturation generating functionβS(t)A(t)/(1+γA(t)) as [18] does. In the other hand, alcoholism will inevitably be affected by many uncertainties, such as emotions and environment. However, the models in the above-mentioned references did not take these random factors into account. In this paper, based on model (2), we will consider the random disturbances of the mortality in all populations due to the effect of alcohol abuse. To do that, we let (Ω,F,{Ft}t≥0,P) be a complete probability space with a filtration {Ft}t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F0 contains all P-null sets) and let Bi(t),1≤i≤4 be 4 independent standard Brownian motions. In practice, we usually estimate a parameter by an average value plus an error term. In this case, the parameters μs,μa,μt, and μq in (2) change to random variables μs~,μa~,μt~, and μq~, respectively, such that (3)μs~=μs+error1,μa~=μa+error2,μt~=μt+error3,μq~=μq+error4.By the central limit theorem, the error terms may be approximated by normal distributions with mean zero and variance σi2dt,1≤i≤4, respectively, so we represent them as follows: (4)erroridt=σidBit,1≤i≤4.With the above crucial factors considered, in this paper, based on model (2), we can get a new stochastic alcoholism dynamic model which is characterized by(5)dSt=Π-βStAt1+γAt-μsStdt+σ1StdB1t,dAt=βStAt1+γAt-ξ1At+ξ2Tt-μaAtdt+σ2AtdB2t,dTt=ξ1At-ξ2Tt-ηTt-μtTtdt+σ3TtdB3t,dQt=ηTt-μqQtdt+σ4QtdB4t.For the parameters in (5) and their explanations, please see Table 1. We will utilize stochastic analysis method to investigate the dynamic behavior of system (5).Table 1
The parameters description of model (5).
Parameter Description μs,μa,μt,μq Death rate inS(t),A(t),T(t),Q(t) accordingly β Transmission coefficient betweenS(t) and A(t) Π The recruitment constant of population γ Half saturation coefficient ξ1 The ratio of alcoholics in the treatment ξ2 The ratio of the treatment population who are failed and return to be alcoholics η The ratio of successfully treated population and never drink hereafter σ i, i=1,2,3,4 Disturbance intensity according toμs,μa,μt,μq, respectively B i ( t ), i=1,2,3,4 Brown movement followed byμs,μa,μt,μq, respectivelyThis paper is arranged as follows. In Section2, we put forward preliminaries including some tools for stochastic analysis and basic properties of deterministic model corresponding to (5) with σi=0, i=1,2,3,4. In Section 3, we discuss the existence and uniqueness of the positive solution of (5). In Section 4, we discuss the stochastic stability of alcohol-free equilibrium point E0. In Section 5, we discuss the stochastic stability of the internal alcoholism equilibrium point E∗. In Section 6, we carry out some simulations to support our theoretical results, In the last section, we give some conclusions to end this paper.


## 2. Preliminaries

Now, we give some criteria on the ergodic property. Denote(6)R+l=x∈Rl:xi>0,∀1≤i≤l.In general, let X be a regular temporally homogeneous Markov process in El⊂R+l described by the stochastic differential equation(7)dXt=bXtdt+∑r=1dσrXtdBrt,with initial value X(t0)=x0∈El and Br(t), 1≤r≤d, are standard Brownian motions defined on the above probability space. The diffusion matrix is defined as follows: (8)Ax=Aijx1≤i,j≤l,Aijx=∑r=1dσrixσrjx.Define the differential operator L associated with (7) by (9)L=∑i=1lbix∂∂xi+12∑i,j=1lAijx∂2∂xi∂xj.If L acts on a function V∈C2,1(El×R+;R), then (10)LVx=∑i=1lbix∂V∂xi+12∑i,j=1lAijx∂2V∂xi∂xj,where Vx=(∂V/∂x1,…,∂V/∂xl) and Vxx=∂2V/∂xi∂xjl×l. By Itô’s formula, we have (11)dVXt=LVXtdt+∑r=1dVxXtσrXtdBrt.Lemma 1 (see [19]).
We assume that there exists a bounded domainU⊂El with regular boundary, having the following properties:(B.1)
In the domainU and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix A(x) is bounded away from zero.(B.2)
Ifx∈El∖U, the mean time τ at which a path issuing from x reaches the set U is finite, and supx∈KExτ<∞ for every compact subset K⊂El.
Then, the Markov processX(t) has a stationary distribution υ(·) with density in El such that for any Borel set B⊂El(12)limt→∞⁡Pt,x,B=υB,PxlimT→∞⁡1T∫0Tfxtdt=∫Elfxυdx=1,for all x∈El and f(x) being a function integrable with respect to the probability measure υ.Remark 2.
(i) The existence of the stationary distribution with density is referred to Theorem 4.1 on page 119 and Lemma 9.4 on page  138 in [19] while the ergodicity and the weak convergence are referred to Theorem 5.1 on page 121 and Theorem 7.1 on page  130 in [19].
(ii) To verify Assumptions (B.1) and (B.2), it suffices to show that there exists a bounded domainU with regular boundary and a nonnegative C2-function V such that A(x) is uniformly elliptical in U and, for any x∈El∖U, LV(x)≤-C for some C>0 (see, e.g., [20, page 1163]).


### 2.1. Basic Properties of Deterministic Model

To compare some results between the deterministic model and stochastic model, firstly, we letσi=0, i=1,2,3,4 to get the corresponding deterministic model as follows:(13)S′t=Π-βStAt1+γAt-μsSt,A′t=βStAt1+γAt-ξ1At+ξ2Tt-μaAt,T′t=ξ1At-ξ2Tt-ηTt-μtTt,Q′t=ηTt-μqQt.We can derived the alcohol-free equilibriumE0 of the deterministic model corresponding to (13); that is, E0=(Π/μs,0,0,0). By using the method of next generation matrix [21], we can calculate to get the fundamental reproduction number of alcoholism which is (14)R0=βΠμsξ2+η+μtξ1η+μt+μaξ2+η+μt.The biological meaning ofR0 is as follows. Once an alcoholic is placed in the environment full of healthy persons with the initial population Π/μs, during the course of alcohol transmission time 1/(ξ1+μa), the alcoholic will infect the healthy persons. R0 denotes the number of healthy persons who are successfully infected.DenotingN(t)=S(t)+A(t)+T(t)+Q(t) and adding up the four equations in (13), we can get (15)N′t=Π-μsS-μaA-μtT-μqQ≤Π-μNt,where μ=min⁡{μs,μa,μt,μq}. Obviously, as t→∞, N(t)≤Π/μ. Hence, all the solutions of (13) will fall in (16)Ω=S,A,T,Q∈R+4∣S+A+T+Q≤Πμ,which is easily proven to be an invariant set. Based on the analysis of [16], we list the relevant results of (13) as follows.Theorem 3.
For system (13), when R0<1, only an alcohol-free equilibrium E0=(Π/μs,0,0,0) exists and it is globally asymptotically stable; when R0>1, besides E0 (it is not stable), there exists an internal equilibrium E∗=(S∗,A∗,T∗,Q∗), which is globally asymptotically stable.


## 2.1. Basic Properties of Deterministic Model

To compare some results between the deterministic model and stochastic model, firstly, we letσi=0, i=1,2,3,4 to get the corresponding deterministic model as follows:(13)S′t=Π-βStAt1+γAt-μsSt,A′t=βStAt1+γAt-ξ1At+ξ2Tt-μaAt,T′t=ξ1At-ξ2Tt-ηTt-μtTt,Q′t=ηTt-μqQt.We can derived the alcohol-free equilibriumE0 of the deterministic model corresponding to (13); that is, E0=(Π/μs,0,0,0). By using the method of next generation matrix [21], we can calculate to get the fundamental reproduction number of alcoholism which is (14)R0=βΠμsξ2+η+μtξ1η+μt+μaξ2+η+μt.The biological meaning ofR0 is as follows. Once an alcoholic is placed in the environment full of healthy persons with the initial population Π/μs, during the course of alcohol transmission time 1/(ξ1+μa), the alcoholic will infect the healthy persons. R0 denotes the number of healthy persons who are successfully infected.DenotingN(t)=S(t)+A(t)+T(t)+Q(t) and adding up the four equations in (13), we can get (15)N′t=Π-μsS-μaA-μtT-μqQ≤Π-μNt,where μ=min⁡{μs,μa,μt,μq}. Obviously, as t→∞, N(t)≤Π/μ. Hence, all the solutions of (13) will fall in (16)Ω=S,A,T,Q∈R+4∣S+A+T+Q≤Πμ,which is easily proven to be an invariant set. Based on the analysis of [16], we list the relevant results of (13) as follows.Theorem 3.
For system (13), when R0<1, only an alcohol-free equilibrium E0=(Π/μs,0,0,0) exists and it is globally asymptotically stable; when R0>1, besides E0 (it is not stable), there exists an internal equilibrium E∗=(S∗,A∗,T∗,Q∗), which is globally asymptotically stable.


## 3. The Existence and Uniqueness of the Positive Solution of (5)

For the sake of biological and practical meanings, in this section, we will investigate the long-term behaviors of system (5); that is, no matter how much the disturbance intensity is, it has a unique global positive solution. According to classical theory of the stochastic differential equations, in order to make sure the existence and uniqueness of the global solution of system (5) are under the given initial conditions (i.e., it will not blow up in a limited time), we generally require that the coefficient of stochastic system (5) satisfies the linear growth condition and the local Lipschitz condition [22]. However, the coefficient of system (5) does not satisfy the linear growth condition; therefore, the solution is likely to blow up in a finite time. In the following, we will use the analysis method of Lyapunov functional to prove global existence of positive solution.Theorem 4.
For arbitrary initial value(S(0),A(0),T(0),Q(0))∈R+4, as t≥0, system (5) shows a unique positive solution (S(t),A(t),T(t),Q(t)), which fall in R+4 almost with probability 1; that is, (S(t),A(t),T(t),Q(t))∈R+4,a.s.Proof.
We firstly prove that system (5) shows a unique local solution. Noticing that the right hand of system (5) does not satisfy the local Lipschitz condition, we introduce variable transformations S=eu, A=eυ; thus, system (5) can be converted into the one whose coefficients satisfy the local Lipschitz condition.
Next, we prove that the only solution is positive and nonexplosive. We assume that system (5) only shows the unique local solution (S(t),A(t),T(t),Q(t)), t∈[0,τe), where τe means the explosive time. To show that the solution is global, we should prove τe=∞, a.e. There exists a sufficient large m0 such that (S(0),A(0),T(0),Q(0))∈[1/m0,m0]. For arbitrary positive integer m≥m0, we define the stopping time (17)τm=inf⁡t∈0,τe:min⁡St,At,Tt,Qt≤1mormax⁡St,At,Tt,Qt≥m,where inf⁡ϕ=∞ (ϕ means empty set). Obviously, as m→∞, τm is monotonically increasing. Due to τm≤τe, sequence {τm} is convergent. Let τ∞=limm→∞τm, so τ∞≤τe a.s. If we can show τ∞=∞, so τe=∞ is obviously established.
We take apagoge to prove this conclusion. Assumeτ∞≠∞ (i.e., τ∞<∞); then there must be a large enough positive constant T and a small enough positive constant ϵ∈(0,1) such that (18)Pτ∞≤T>ϵ.
Hence, there exists a positive integerm1≥m0 such that(19)Pτm≤T≥ϵ,∀m≥m1.We define a Lyapunov functional V:R+4→R+ as follows: (20)VS,A,T,Q=S-1-ln⁡S+A-1-ln⁡A+T-1-ln⁡T+Q-1-ln⁡Q.By the simple inequality of u-1-ln⁡u≥0, ∀u>0, it is easy to know that V(S,A,T,Q) is nonnegative. In virtue of Itô formula, we calculate the stochastic differential of V(S,A,T,Q) along with (5) as follows: (21)dV=1-1SdS+12S2dS2+1-1AdA+12A2dA2+1-1TdT+12T2dT2+1-1QdQ+12Q2dQ2=Π-μsS-ΠS+βA1+γA+μs+σ122-μaA-βS1+γA+ξ1-ξ2TA+μa+σ222-μtT-ξ1AT+ξ2+η+μt+σ322-μqQ-ηTQ+μq+σ422dt+1-1Sσ1SdB1+1-1Aσ2AdB2+1-1Tσ3TdB3+1-1Qσ4QdB4≤1+1μΠ+βγ+μs+μa+μt+μq+ξ1+ξ2+η+σ12+σ22+σ32+σ422dt+σ11-1SSdB1+σ21-1AAdB2+σ31-1TTdB3+σ41-1QQdB4≔LVdt+σ11-1SSdB1+σ21-1AAdB2+σ31-1TTdB3+σ41-1QQdB4,where(22)LV=1+1μΠ+βγ+μs+μa+μt+μq+ξ1+ξ2+η+σ12+σ22+σ32+σ422≔K.
Integrating the two sides ofdV above from 0 to τm∧T, we get (23)∫0τm∧TdVSl,Al,Tl,Ql≤∫0τm∧TKdl+∫0τm∧Tσ11-1SlSldB1l+σ21-1AlAldB2l+σ31-1TlTldB3l+σ41-1QlQldB4l.
Take expectation of the two sides in the above expression to get(24)EVSτm∧T,Aτm∧T,Tτm∧T,Qτm∧T≤VS0,A0,T0,Q0+E∫0τm∧TKdl≤VS0,A0,T0,Q0+KT.
LetΩm={τm≤T}; there exists m1>0, for ∀m>m1; by (19), we know P(Ωm)≥ϵ.
∀ ω ∈ Ω m, at least one of S(τm),A(τm),T(τm) and (τm) equals m or 1/m. Specifically, when S(τm,ω)=m or 1/m, A(τm,ω)=m or 1/m, T(τm,ω)=m or 1/m, or Q(τm,ω)=m or 1/m, there will always be (25)VSτm∧T,Aτm∧T,Tτm∧T,Qτm∧T≥m-1-ln⁡m∧1m-1-ln⁡1m.
Comprehensive consideration of (19) and (24), we can get (26)VS0,A0,T0,Q0+KT≥E1ΩmωVSτm∧T,Aτm∧T,Tτm∧T,Qτm∧T≥ϵm-1-ln⁡m∧1m-1-ln⁡1m,where 1Ωm denotes the indicator function of Ωm. By letting m→∞, it is easy to get the following contradiction: (27)∞>VS0,A0,T0,Q0+KT=∞.Thus, τ∞=∞ a.e. The proof is complete.Remark 5.
This theorem implies that, besides the alcohol free equilibrium, system (5) will always show a unique positive solution, no matter how much the disturbance intensity is.


## 4. Stochastic Stability of Alcohol-Free Equilibrium PointE0

Seen from Theorem4, for deterministic system (13), as R0<1, alcohol-free equilibrium point E0=(Π/μs,0,0,0) is globally asymptotically stable, which means the alcohol population will disappear in the appropriate conditions finally, and the alcohol behavior will be effectively controlled. However, there is not any equilibrium point in the stochastic system (5); therefore, in this section, we will discuss the disturbance behavior of system (5) by investigating the stability and the ergodic property of E0.Theorem 6.
Assume(S(t),A(t),T(t),Q(t)) is the solution of system (5) with the initial value (S(0),A(0),T(0),Q(0))∈R+4. If(28)R0=βΠμs·ξ2+η+μtξ1η+μt+μaξ2+η+μt≤1,then(29)limsupT→∞⁡1T∫0Tμs-c1St-S02+μa2A2-c2A2t+η+μt2-c3T2tdt≤c1S02,where ci, i=1,2,3 are defined as follows: (30)c1=12μsμs+μa22μa+μs+η+μt22η+μt+12σ12,c2=σ222,c3=12+μa+η+μt2ξ1σ32.Proof.
Define nonnegativeC2-function V1:R+→R+ as (31)V1S=S-S022;then(32)dV1=S-S0dS+12dS2=S-S0Π-βSA1+γA-μsSdt+σ1SdB1t+12σ1S2dt.Hence,(33)LV1=S-S0Π-βSA1+γA-μsS+12σ1S2=S-S0μsS0-βSA1+γA-μsS+12σ1S2≤-μsS-S02-βS0S-S0A1+γA+12σ1S2.Let(34)V2A,T=A+ξ2η+μt+ξ2T;then(35)dV2=dA+ξ2η+μt+ξ2dT=βSA1+γA-ξ1A+ξ2T-μaAdt+σ2AdB2t+ξ2η+μt+ξ2ξ1A-ξ2T-ηT-μtTdt+σ3TdB3t.Hence,(36)LV2=βSA1+γA-ξ1A+ξ2T-μaA+ξ2η+μt+ξ2ξ1A-ξ2T-ηT-μtT=βS-S0A1+γA+βS0A1+γA-ξ1+μaA+ξ1ξ2Aη+μt+ξ2=βS-S0A1+γA+AβS01+γA-ξ1+μa-ξ1ξ2η+μt+ξ2≤βS-S0A1+γA+AβS0-ξ1η+μtη+μt+ξ2-μa≤βS-S0A1+γA.The last inequality is attributed to the following relation: (37)βS0-ξ1η+μtη+μt+ξ2-μa≤0⟺R0=βΠμs·ξ2+η+μtξ1η+μt+μaξ2+η+μt≤1.
Consider linear combinationV1+S0V2; then (38)LV1+S0V2≤-μsS-S02+12σ1S2.Let(39)V3=T22;then(40)dV3=TdT+12dT2=Tξ1A-ξ2T-ηT-μtTdt+σ3TdB3t+12σ3T2dt,LV3=Tξ1A-ξ2T-ηT-μtT+12σ3T2=ξ1AT-η+μt+ξ2T2+12σ3T2.Let(41)V4=S-S0+A+T22;then(42)dV4=S-S0+A+TdS+dA+dT+dS2+dA2+dT22=S-S0+A+TΠ-μsS-μaA-μtT-ηTdt+σ1SdB1t+σ2AdB2t+σ3TdB3t+σ1S2+σ2A2+σ3T22dt,LV4=S-S0+A+TΠ-μsS-μaA-μtT-ηT+σ1S2+σ2A2+σ3T22=-μsS-S02-μaA2-η+μtT2-μs+μaS-S0A-μs+η+μtS-S0T-μa+η+μtAT+σ1S2+σ2A2+σ3T22≤-μsS-S02-μa2A2-η+μt2T2+μs+μa22μaS-S02+μs+η+μt22η+μtS-S02-μa+η+μtAT+σ1S2+σ2A2+σ3T22=-μs+μs+μa22μa+μs+η+μt22η+μtS-S02-μa2A2-η+μt2T2-μa+η+μtAT+σ1S2+σ2A2+σ3T22.
In order to eliminate the cross termAT, we consider linear combination of V3 and V4 as ((μa+η+μt)/ξ1)V3+V4; then (43)Lμa+η+μtξ1V3+V4≤-μs+μs+μa22μa+η+μt+μs22η+μtS-S02-μa2A2-η+μt2T2+σ1S2+σ2A2+σ3T22+μa+η+μtξ1σ3T22.
At last, to eliminate the term(S-S0)2, we consider linear combination (44)1μsμs+μa22μa+η+μt+μs22η+μtV1+S0V2+μa+η+μtξ1V3+V4.
LetV=(1/μs)[(μs+μa)2/2μa+(μs+η+μt)2/2η+μt](V1+S0V2)+((μa+η+μt)/ξ1)V3+V4. Integrating all the formulas of LVi, then (45)LV≤-μsS-S02-μa2A2-η+μt2T2+12μsμs+μa22μa+μs+η+μt22η+μt+12σ1S2+12σ2A2+12+μa+η+μt2ξ1σ3T2≔-μsS-S02-μa2A2-η+μt2T2+c1S2+c2A2+c3T2,where ci(i=1,2,3) is defined as follows: (46)c1=12μsμs+μa22μa+μs+η+μt22η+μt+12σ12,c2=σ222,c3=12+μa+η+μt2ξ1σ32.Seen from the definition of Vi, i=1,2,3,4, (47)V=1μsμs+μa2μa+η+μt+μs2η+μtV1+S0V2+η+μa+μtξ1V3+V4is obviously positive definite; it is easy to know that (48)limT→∞⁡1T∫0Tμs-c1St-S02+μa2A2-c2A2t+η+μt2-c3T2tdt≤c1S02.The proof is complete.Remark 7.
Besides,Q0 becomes the disease-free equilibrium of system (1) as σ1=0. In the proof, we see(49)LV≤-μsS-S02-μa2-c2A2-η+μt2-c3T2,
Thus, the solution of system is stochastically asymptotically stable in the large [22] as the conditions μa/2>c2, (η+μt)/2>c3 are satisfied.


## 5. Stochastic Stability of the Internal Equilibrium PointE∗

In this section, we will investigate the stochastic behavior of (5) around the internal equilibrium point E∗ in deterministic model (13) from two aspects; that is, one is concerned with the stability of its distribution, the other is concerned with ergodic property of equilibrium point E∗.Theorem 8.
If the following conditions are satisfied,(50)μs>CS∗2+c1,μa>2CA∗2+c2,η+μt>2CT∗2+c3,where(51)c1=σ12+μs+μa22μa+μs+μt22η+μtσ12μs,c2=σ22,c3=σ32+μa+μtξ1σ12,c4=1μsμs+μa22μa+μs+μt22η+μtS∗1+γA∗σ222+ξ2T∗2ξ1A∗S∗1+γA∗σ322,C=c1S∗2+c2A∗2+c3T∗2+c4,then (5) is ergodic and positive recurrent.Proof.
Construct Lyapunov functionU1 as (52)U1=S-S∗22and calculate to get (53)LU1=S-S∗Π-βSA1+γA-μsS+σ122S2=S-S∗βS∗A∗1+γA∗-βSA1+γA-μsS-S∗+σ122S2=βS-S∗S∗A∗1+γA∗-SA1+γA-μsS-S∗2+σ122S2=βS-S∗S∗A∗11+γA∗-11+γA+S∗A∗-SA1+γA-μsS-S∗2+σ122S2=βγS∗A∗A-A∗S-S∗1+γA∗1+γA-βS∗S-S∗A-A∗1+γA-βAS-S∗21+γA-μsS-S∗2+σ122S2≤-βS∗S-S∗A-A∗1+γA1+γA∗+σ1S22-μsS-S∗2.Next, we let (54)U2=A-A∗-A∗ln⁡AA∗,and calculate to get (55)LU2=1-A∗AβSA1+γA-ξ1A+ξ2T-μaA+σ222=1-A∗AβSA1+γA-ξ1+μaA+1-A∗Aξ2T+σ222=A-A∗βS1+γA-ξ1+μa+ξ2T-ξ2TA∗A+σ222=A-A∗βS1+γA-βS∗1+γA∗+ξ2T∗A∗+ξ2T-ξ2TA∗A+σ222=βA-A∗S1+γA-S∗1+γA∗+ξ2T∗A∗A-A∗+ξ2T-ξ2TA∗A+σ222=βA-A∗S-S∗1+γA+βS∗A-A∗11+γA-11+γA∗+ξ2T∗A∗A-A∗+ξ2T-ξ2TA∗A+σ222≤βA-A∗S-S∗1+γA+ξ2T∗A∗A-A∗+ξ2T-ξ2TA∗A+σ222.Then, it is easy to know (56)LU1+S∗1+γA∗U2≤-μsS-S∗2+σ12S22+S∗1+γA∗ξ2T∗A∗A-A∗+S∗1+γA∗ξ2T-S∗1+γA∗ξ2TA∗A+S∗1+γA∗σ222=-μsS-S∗2+σ12S22+S∗1+γA∗ξ2T∗A∗A+S∗1+γA∗ξ2T-T∗-S∗1+γA∗ξ2TA∗A+S∗1+γA∗σ222=-μsS-S∗2+σ12S22-S∗1+γA∗ξ2T∗AA∗+TT∗+TT∗A∗A-1+σ12S22+S∗1+γA∗σ222≤-μsS-S∗2-S∗1+γA∗ξ2T∗AA∗-logA∗A-1+TT∗-log⁡T∗T-1+σ12S22+S∗1+γA∗σ222,where the last inequality is derived from x≥log⁡x+1, x>0.
Let(57)U3=S-S∗+A-A∗+T-T∗22;we calculate to get(58)LU3=S-S∗+A-A∗+T-T∗Π-μsS-μaA-η+μtT+σ122S2+σ222A2+σ322T2=S-S∗+A-A∗+T-T∗-μsS-S∗-μaA-A∗-η+μtT-T∗+σ122S2+σ222A2+σ322T2=-μsS-S∗2-μaA-A∗2-η+μtT-T∗2-μs+μaS-S∗A-A∗-μs+μtS-S∗T-T∗-μa+μtA-A∗T-T∗+σ122S2+σ222A2+σ322T2.Let(59)U4=T-T∗22;similar to the technique in computing LU1, we calculate to get (60)LU4=T-T∗ξ1A-A∗-ξ2+η+μtT-T∗+σ122T2=ξ1A-A∗T-T∗-ξ2+η+μtT-T∗2+σ122T2.Considering the linear combination of U3 and U4 as (61)U3+μa+μtξ1U4,then (62)LU3+μa+μtξ1U4=-μsS-S∗2-μa2A-A∗2-η+μt2+μa+μtξ1ξ2+η+μtT-T∗2+μs+μa22μaS-S∗2+μs+μt22η+μtS-S∗2+μa+μtξ1σ122T2+σ122S2+σ222A2+σ322T2.
Next, we define(63)U5=T-T∗-T∗log⁡TT∗;then ξ1A∗=(ξ2+η+μt)T∗ implies (64)LU5=1-T∗Tξ1A-ξ2+η+μtT+σ32T∗2=ξ1A∗AA∗-AA∗T∗T-TT∗+1+σ32T∗2≤ξ1A∗AA∗-log⁡AA∗-1-TT∗-log⁡TT∗-1+σ32T∗2.
We let(65)V=U1+S∗1+γA∗U2+ξ2T∗ξ1A∗S∗1+γA∗U5;it is easy to prove that (66)LV≤-μsS-S∗2+σ12S22+S∗1+γA∗σ222+ξ2T∗2ξ1A∗S∗1+γA∗σ322.
We consider the linear combination of(67)μs+μa22μa+μs+μt22η+μtVμs+U3+μa+μtξ1U4;then(68)Lμs+μa22μa+μs+μt22η+μtVμs+U3+μa+μtξ1U4≤-μsS-S∗2-μa2A-A∗2-η+μt2+μa+μtξ1ξ2+η+μtT-T∗2+c1S22+c2A22+c3T22+c4≤-μs-c1S-S∗2-μa2-c2A-A∗2-η+μt2-c3T-T∗2+C,where ci, i=1,2,3,4 is represented as (69)c1=σ12+μs+μa22μa+μs+μt22η+μtσ12μs,c2=σ22,c3=σ32+μa+μtξ1σ12,c4=1μsμs+μa22μa+μs+μt22η+μtS∗1+γA∗σ222+ξ2T∗2ξ1A∗S∗1+γA∗σ322,C=c4+c1S∗2+c2A∗2+c3T∗2
Choose appropriate parameters in system (2) to satisfy (70)μs-c1S∗2>C,μa2-c2A∗2>C,η+μt2-c3T∗2>C;that is,(71)μs>CS∗2+c1,μa>2CA∗2+c2,η+μt>2CT∗2+c3.
Hence, if we set(72)U=μs+μa22μa+μs+μt22η+μtVμs+U3+μa+μtξ1U4>0,then LU≤-C for some C>0. The proof is complete.Next, we will discuss the mean convergence of system (5) based on the condition in Theorem 6.Theorem 9.
If the condition in Theorem6 holds, then the following conclusions are established: (73)limT→∞⁡1T∫0TβSuAu1+γAudu+limT→∞⁡1T∫0TμsSudu=Π,a.e.,limT→∞⁡1T∫0TβSuAu1+γAudu+limT→∞⁡1T∫0Tξ2Tudu=limT→∞⁡1T∫0Tξ1Au+μaAudu,a.e,limT→∞⁡1T∫0Tξ1Audu=limT→∞⁡1T∫0Tξ2Tu+ηTu+μtTudu,a.e.Proof.
By ergodic theorem, for arbitrary functionf, if it satisfies (74)∫fx,y,zmdx,dy,dz<+∞,then(75)limT→∞⁡1T∫0TfSu,Au,Tudu⟶(76)∫fx,y,zmdx,dy,dz<+∞,a.e.Note C=min{μs,μa,η+μt}, since (77)dSt+At+Tt=Π-μsSt-μaAt-η+μtTt+σ1StdB1t+σ2AtdB2t+σ3TtdB3t;we take stochastic integral of the above formula over [0,t], 0<t<T, (78)St+At+Tt=S0+A0+T0+Πt-∫0tμsSu+μaAu+η+μtTudu+∫0tσ1SudB1u+∫0tσ2AudB2u+∫0tσ3TudB3u,ESt+At+Tt≤S0+A0+T0-C∫0tESu+Au+Tudu;by Gronwall inequality, we can derive (79)0≤ESt+At+Tt<M<+∞.Furthermore, (80)limT→∞⁡1T∫0TSt+At+Tt∧ndt=∫x+y+z∧nmdx,dy,dz.After taking expectation of this formula, by the control convergence theorem, we know (81)∫R+3x+y+z∧nmdx,dy,dz≤M<+∞.
Letn→+∞; then (82)∫R+3x+y+zmdx,dy,dz<+∞.Therefore, if we define f=βSA/(1+γA) or S or A or T, respectively, we can prove the corresponding convergence.
By the use of (5), we can get (83)St=∫0tΠ-βSuAu1+γAu-μsSudu+σ1SudB1u.Hence,(84)Stt=1t∫0tΠ-βSuAu1+γAu-μsSudu+1t∫0tσ1SudB1u.By the law of large numbers of martingale, (85)limt→∞⁡1t∫0tσ1SudB1u⟶0.Considering that limt→∞(S(t)/t) exists, then (86)Elimt→∞⁡Stt=limt→∞⁡EStt.Since ES(t)≤M,∀t, then (87)Elimt→∞⁡Stt=0.Thus, we proved (88)Stt⟶0,a.e.
Similarly, we can prove(89)Att⟶0,Ttt⟶0,a.e.From (84), we can conclude (90)limT→∞⁡1T∫0TΠ-βSuAu1+γAu-μsSudu=0.The others conclusions can be similarly proved. The proof is completed.


## 6. Numerical Simulations

In this section, to verify the results of stochastic stability, we use Milstein’s numerical method to make corresponding simulations based on stochastic model (5) with some given initial values and parameters [23]. Thus, discreting system (5), we get the following iterative scheme: (91)Si+1=Si+Π-βSiAi1+γAi-μsSiΔt+σ1Sic1,iΔt+0.5σ12Sic1,i2-1Δt,Ai+1=Ai+βSiAi1+γAi-ξ1Ai+ξ2Ti-μaAiΔt+σ2Aic2,iΔt+0.5σ22Aic2,i2-1Δt,Ti+1=Ti+ξ1Ai-ξ2Ti-ηTi-μtTiΔt+σ3Tic3,iΔt+0.5σ32Tic2,i2-1Δt,Qi+1=Qi+ηTi-μqQiΔt+σ4Qic4,iΔt+0.5σ42Qic4,i2-1Δt,where c(l,i), l=1,2,3,4; i=1,2,…,n, are independent Gaussian random variables N(0,1) and σi,i=1,2,3,4, are intensities of white noise.We choose appropriate initial values(S(0),A(0),T(0),Q(0))=(6000,70,20,5) and system parameter values as γ=0.2 as well as population mortality are μs=0.001, μa=0.005, μt=0.003, and μq=0.0015, respectively, by the latest report from WHO [7]. In addition, some parameters are reasonably estimated as η=0.3, ξ1=ξ2=0.2. In our simulations, in view of the fact that drinking is related to the people’s habit and thus is a long-term behavior, for the sake of simplicity, we will only consider the evolution curves of S(t) and A(t). Furthermore, we will select different critical parameters to discuss stability of the stochastic solution of (5) as follows.(1) We choose population recruitment rate asΠ=0.1 and alcohol infection rate as β=0.001 to make R0=0.797<1. To compare we let disturbance intensity be σi=0.05 and0.01(i=1,2,3,4.), respectively (see Figures 1 and 2).Figure 1
R 0 < 1, population evolution curve of S(t) with disturbance intensity σi=0.05,0.01(i=1,2,3,4).Figure 2
R 0 < 1, population evolution curve of A(t) with disturbance intensity σi=0.05,0.01(i=1,2,3,4.).Synthesizing and comparing the information in Figures1 and 2, we can conclude that when R0<1, regardless of the size of the disturbance intensity, the solution of stochastic system (5) will randomly tend to the alcohol-free equilibrium E0=(Π/μs,0,0,0) of deterministic system (13). Specifically speaking, if we can take some effective measures, for example, controlling the population recruitment rate Π or reducing the infection rate β, the population in A(t),T(t),Q(t) will eventually tend to zero while the population in S(t) will eventually tend to a constant Π/μs. Figures 1 and 2 agree with the conclusions in Theorems 6 and 9.(2) We choose population recruitment rate asΠ=0.2 and alcohol infection rate as β=0.002 to make R0=1.59>1. To compare, we still let disturbance intensity be σi=0.05 and 0.01(i=1,2,3,4), respectively (see Figures 3 and 4).Figure 3
R 0 > 1, population evolution curve of S(t) with disturbance intensity σi=0.05,0.01(i=1,2,3,4.).Figure 4
R 0 > 1, population evolution curve of A(t) with disturbance intensity σi=0.05,0.01(i=1,2,3,4.).Similarly, by synthesizing and comparing the information in Figures3 and 4, we can conclude that, in one hand, when R0>1, regardless of the size of the disturbance intensity, the solution of stochastic system (5) will randomly tend to the internal equilibrium point E∗=(S∗,A∗,T∗,Q∗) of deterministic system corresponding to σi=0(i=1,2,3,4) in (5). Specifically speaking, if the population recruitment rate Π or the infection rate β is rather large, the population in A(t),T(t),Q(t) will increase with time going and eventually tends to a fixed level, while the population in S(t) will decrease with time going and eventually tends to a constant. In the other hand, the disturbance intensity of the system will affect the speed of the trend; specifically, the greater the disturbance, the slower this trend, and vice versa. Figures 3 and 4 agree with the conclusions in Theorems 8 and 9.(3) To testify the ergodic property of system (5), using the same parameters as case (2), we draw Figures 5 and 6 as follows only with relatively large disturbance intensity σi=0.05, i=1,2,3,4.Figure 5
R 0 > 1, ergodic property of S(t) with disturbance intensity σi=0.05(i=1,2,3,4.).Figure 6
R 0 > 1, ergodic property of A(t) with disturbance intensity σi=0.05(i=1,2,3,4.).Figures5 and 6 reveal the fact that if system parameters meet the conditions to make R0>1, then the model is ergodic and has a uniqueness stationary distribution which illustrates the persistence. This agrees with the conclusions in Theorem 9.


## 7. Conclusions

In this paper, with the fact that the addiction of alcohol abuse complies with the law of infectious diseases being taken into account, we formulate a stochastic alcoholism model to investigate the stochastic stability of the model. We find that, ifR0=β(Π/μs)((ξ2+η+μt)/[ξ1(η+μt)+μa(ξ2+η+μt)])<1, the alcohol-free equilibrium E0=(Π/μs,0,0,0) is stochastic asymptotic stability irrelated to system disturbance intensity, which points a way to control alcohol by limiting recruitment to a drinking environment and avoiding contact with the alcoholics as far as possible. On the contrary, when system parameters meet the conditions to make R0>1, then the model is ergodic and has a uniqueness stationary distribution which illustrates the persistence. Furthermore, it is worth noting that, in the two cases (i.e., when R0<1 or R0>1), disturbance intensity will obviously affect the trend speed; specifically, the smaller the disturbance, the faster the trend, and vice versa.



---

*Source: 6297074-2017-01-02.xml*