# Global Stability for a Binge Drinking Model with Two Stages

Authors: Hai-Feng Huo; Na-Na Song
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

A more realistic two-stage model for binge drinking problem is introduced, where the youths
with alcohol problems are divided into those who admit the problem and those who do not
admit it. We also consider the direct transfer from the class of susceptible individuals towards the
class of admitting drinkers. Mathematical analyses establish that the global dynamics of the model
are determined by the basic reproduction number,R0. The alcohol-free equilibrium is globally asymptotically stable, and the alcohol problems are eliminated from the population if R0<1. A unique
alcohol-present equilibrium is globally asymptotically stable if R0>1. Numerical simulations are also conducted in the analytic results.


---

## Body

## 1. Introduction

Young people’s binge drinking problem is a major concern to public health. Recently, US surveys indicate that approximately 90% of college students have consumed alcohol at least once [1], and more than 40% of college students have engaged in binge drinking [2, 3]. The binge drinking refers to youths in 17–30 age group who drink a large amount of alcohol and become so drunk; they are likely to exhibit antisocial behavior [4]. Although there have been many attempts to reduce the problem, alcohol abuse by college students has persisted and in some cases increased over the past several decades [5]. Prior studies have indicated that heavy alcohol drinkers are likely to engage in risky sexual behaviours and more likely to get sexually transmitted infection than social drinkers [6, 7]. There is a strong medical evidence that treatment of individuals with alcohol problems is a major issue [8–10].Thus, it is very important to use a mathematical method to study the binge drinking problems in youths. A simple model for alcohol treatment is presented by Sanchez et al. [11]. Since then, there have been numerous studies investigating campus drinking and the associated consequences [12–15]. Manthey et al. [12] focus on a college campus, divide the student population into three classes: nondrinkers, social drinkers and problem drinkers, and show that campus alcohol abuse may be reduced by minimizing the ability of problem drinkers to directly recruit nondrinkers. Cintron-Arias et al. [13] focus on situations where relapse rates are high and conclude that the systematic removal of individuals from high-risk environments, or the development of programs that limit access or reduce the residence times in such environments (or both approaches combined), may reduce the level of alcohol abuse. Mubayi et al. [14] show that if the relative residence times of moderate drinkers are distributed randomly between low- and high-risk environments, then the proportion of heavy drinkers is likely to be higher than expected. Mulone and Straughan [15] investigate a model for binge drinking taking into account admitting and nonadmitting drinkers. But the global stability of binge drinking model is not discussed in the literature.Motivated by the binge drinking model in [15], we develop a more realistic model with two stages. Drinking is often encouraged by peer pressure. A susceptible individual acquires alcohol problems through the direct contact with the admitting drinker or the nonadmitting drinker. The new drinker can become either the admitting drinker or the nonadmitting drinker. So we consider the direct transfer from the class of susceptible individuals towards the class of admitting drinkers; furthermore, we study the global dynamics of the model. The reason to introduce this new direct transfer is that about one-third of the American population admit to drinking problems, 17.8% of the population admit to the binge drinking problem [16]. Therefore, this fact cannot be neglected in the binge drinking model.The organization of this paper is as follows. In the next section, the binge drinking model with two stages and some basic properties are derived. In Section3, the existence and the global stability of equilibria are investigated. Some numerical simulations are given in Section 4. Some discussions are given in Section 5.


## 2. The Model

### 2.1. System Description

The total population is divided into four compartments, namely, the susceptible compartment of those who do not drink or drink only moderately, denoted byS(t), those who drink heavily at least some of the time but do not admit having a problem, denoted by A1(t), those who drink heavily and admit having a problem, denoted by A2(t), and those people in treatment, denoted by R(t). The total number of population at time t is given by
(2.1)N(t)=S(t)+A1(t)+A2(t)+R(t).
The model structure is shown in Figure 1. The transfer diagram leads to the following system of ordinary differential equations:
(2.2)S˙=μN-S(βA1+γA2)N-μS,A˙1=(1-p)S(βA1+γA2)N-(α+μ)A1,A˙2=pS(βA1+γA2)N+ηR+αA1-(π+μ)A2,R˙=πA2-(μ+η)R,
where μN is the number of individuals entering into the system in a given time interval (say each year), so μ represents the rate of entry. Since we are dealing with youths, we assume that  the death rate is negligible, and so the leaving rate is also μ. A susceptible individual can be turned to drink through direct contact with an admitting drinker or a nonadmitting drinker. β is the transmission coefficient for the individuals who do not admit to have alcohol problems, γ is the transmission coefficient for the individuals who admit to have these problems. We assume that β>γ due to the fact that they are more likely to be unaware of their condition. α is the rate at which represents those with the alcohol problems admitting to have these problems and then transferring from class A1 to A2, π is that fraction of A2 who go into treatment, η is that fraction of R who relapse into admitting drinkers (we adopt a linear relapse term rather than peer pressure since we argue that the relapse is primarily due to the person), p is the probability of a new drinker to admit having the problem. We may show that  N=S+A1+A2+R is constant, and then we introduce the fractions of S, A1, A2, and R:
(2.3)s=SN,a1=A1N,a2=A2N,r=RN,
with s+a1+a2+r=1. Then the system (2.2) becomes
(2.4)s˙=μ-s(βa1+γa2)-μs,a˙1=(1-p)s(βa1+γa2)-(α+μ)a1,a˙2=ps(βa1+γa2)+ηr+αa1-(π+μ)a2,r˙=πa2-(η+μ)r.
As system (2.4) is equivalent to system (2.2), we only need to study system (2.4).Figure 1
Transfer diagram of model (2.2).


### 2.2. Basic Properties

#### 2.2.1. Invariant Region

Adding all the equations of system (2.4) gives
(2.5)dNdt=0.Thus, the total populationN is a constant. Since system (2.4) monitors human population, it is plausible to assume that all its state variables and parameters are nonnegative for all t≥0. Further, it can be shown that the region
(2.6)Ω={(s(t),a1(t),a2(t),r(t))∈R+4:s(t)+a1(t)+a2(t)+r(t)≤1}
is positively invariant. Thus, each solution of system (2.4), with initial conditions in Ω, remains there for t≥0. Therefore, the ω-limit sets of the solutions of system (2.4), in Ω, are contained in Ω. Furthermore, in Ω, the usual existence, uniqueness, and continuation results hold for the system, so that the system (2.4) is well posed mathematically and epidemiologically [17]. So we consider the dynamics of system (2.4) on the set Ω in this paper.


#### 2.2.2. Positivity of Solutions

For system (2.4), it is necessary to prove that all the state variables are positive, so that the solutions of the system with positive initial conditions remain positive for all t≥0. We thus state the following lemma.Lemma 2.1.
Ifs(0)>0,a1(0)>0,a2(0)>0, and r(0)>0, the solutions s(t), a1(t), a2(t), and r(t) of system (2.4) are positive for all t≥0.Proof.
Under the given initial conditions, it is easy to prove that the solutions of system (2.4) are positive; if not, we assume a contradiction: that there exists a first time t1 such that
(2.7)s(t1)=0,s′(t1)<0,a1(t)>0,a2(t)>0,r(t)>0,0<t<t1,
there exists a t2,
(2.8)a1(t2)=0,a1′(t2)<0,s(t)>0,a2(t)>0,r(t)>0,0<t<t2,
there exists a t3,
(2.9)a2(t3)=0,a2′(t3)<0,s(t)>0,a1(t)>0,r(t)>0,0<t<t3,
and there exists a t4,
(2.10)r(t4)=0,r′(t4)<0,s(t)>0,a1(t)>0,a2(t)>0,0<t<t4.
In the first case, we have
(2.11)s′(t1)=μ>0,
which is a contradiction, meaning that s(t)>0, t≥0. In the second case, we have
(2.12)a1′(t2)=(1-p)γa2(t)s(t)>0,
which is a contradiction, meaning that a1(t)>0, t≥0. Similarly, it can be shown that a2(t)>0 and r(t)>0 for all t≥0.Thus, the solutionss(t), a1(t),  a2(t), and r(t) of system (2.4) remain positive for all t>0.


## 2.1. System Description

The total population is divided into four compartments, namely, the susceptible compartment of those who do not drink or drink only moderately, denoted byS(t), those who drink heavily at least some of the time but do not admit having a problem, denoted by A1(t), those who drink heavily and admit having a problem, denoted by A2(t), and those people in treatment, denoted by R(t). The total number of population at time t is given by
(2.1)N(t)=S(t)+A1(t)+A2(t)+R(t).
The model structure is shown in Figure 1. The transfer diagram leads to the following system of ordinary differential equations:
(2.2)S˙=μN-S(βA1+γA2)N-μS,A˙1=(1-p)S(βA1+γA2)N-(α+μ)A1,A˙2=pS(βA1+γA2)N+ηR+αA1-(π+μ)A2,R˙=πA2-(μ+η)R,
where μN is the number of individuals entering into the system in a given time interval (say each year), so μ represents the rate of entry. Since we are dealing with youths, we assume that  the death rate is negligible, and so the leaving rate is also μ. A susceptible individual can be turned to drink through direct contact with an admitting drinker or a nonadmitting drinker. β is the transmission coefficient for the individuals who do not admit to have alcohol problems, γ is the transmission coefficient for the individuals who admit to have these problems. We assume that β>γ due to the fact that they are more likely to be unaware of their condition. α is the rate at which represents those with the alcohol problems admitting to have these problems and then transferring from class A1 to A2, π is that fraction of A2 who go into treatment, η is that fraction of R who relapse into admitting drinkers (we adopt a linear relapse term rather than peer pressure since we argue that the relapse is primarily due to the person), p is the probability of a new drinker to admit having the problem. We may show that  N=S+A1+A2+R is constant, and then we introduce the fractions of S, A1, A2, and R:
(2.3)s=SN,a1=A1N,a2=A2N,r=RN,
with s+a1+a2+r=1. Then the system (2.2) becomes
(2.4)s˙=μ-s(βa1+γa2)-μs,a˙1=(1-p)s(βa1+γa2)-(α+μ)a1,a˙2=ps(βa1+γa2)+ηr+αa1-(π+μ)a2,r˙=πa2-(η+μ)r.
As system (2.4) is equivalent to system (2.2), we only need to study system (2.4).Figure 1
Transfer diagram of model (2.2).


## 2.2. Basic Properties

### 2.2.1. Invariant Region

Adding all the equations of system (2.4) gives
(2.5)dNdt=0.Thus, the total populationN is a constant. Since system (2.4) monitors human population, it is plausible to assume that all its state variables and parameters are nonnegative for all t≥0. Further, it can be shown that the region
(2.6)Ω={(s(t),a1(t),a2(t),r(t))∈R+4:s(t)+a1(t)+a2(t)+r(t)≤1}
is positively invariant. Thus, each solution of system (2.4), with initial conditions in Ω, remains there for t≥0. Therefore, the ω-limit sets of the solutions of system (2.4), in Ω, are contained in Ω. Furthermore, in Ω, the usual existence, uniqueness, and continuation results hold for the system, so that the system (2.4) is well posed mathematically and epidemiologically [17]. So we consider the dynamics of system (2.4) on the set Ω in this paper.


### 2.2.2. Positivity of Solutions

For system (2.4), it is necessary to prove that all the state variables are positive, so that the solutions of the system with positive initial conditions remain positive for all t≥0. We thus state the following lemma.Lemma 2.1.
Ifs(0)>0,a1(0)>0,a2(0)>0, and r(0)>0, the solutions s(t), a1(t), a2(t), and r(t) of system (2.4) are positive for all t≥0.Proof.
Under the given initial conditions, it is easy to prove that the solutions of system (2.4) are positive; if not, we assume a contradiction: that there exists a first time t1 such that
(2.7)s(t1)=0,s′(t1)<0,a1(t)>0,a2(t)>0,r(t)>0,0<t<t1,
there exists a t2,
(2.8)a1(t2)=0,a1′(t2)<0,s(t)>0,a2(t)>0,r(t)>0,0<t<t2,
there exists a t3,
(2.9)a2(t3)=0,a2′(t3)<0,s(t)>0,a1(t)>0,r(t)>0,0<t<t3,
and there exists a t4,
(2.10)r(t4)=0,r′(t4)<0,s(t)>0,a1(t)>0,a2(t)>0,0<t<t4.
In the first case, we have
(2.11)s′(t1)=μ>0,
which is a contradiction, meaning that s(t)>0, t≥0. In the second case, we have
(2.12)a1′(t2)=(1-p)γa2(t)s(t)>0,
which is a contradiction, meaning that a1(t)>0, t≥0. Similarly, it can be shown that a2(t)>0 and r(t)>0 for all t≥0.Thus, the solutionss(t), a1(t),  a2(t), and r(t) of system (2.4) remain positive for all t>0.


## 2.2.1. Invariant Region

Adding all the equations of system (2.4) gives
(2.5)dNdt=0.Thus, the total populationN is a constant. Since system (2.4) monitors human population, it is plausible to assume that all its state variables and parameters are nonnegative for all t≥0. Further, it can be shown that the region
(2.6)Ω={(s(t),a1(t),a2(t),r(t))∈R+4:s(t)+a1(t)+a2(t)+r(t)≤1}
is positively invariant. Thus, each solution of system (2.4), with initial conditions in Ω, remains there for t≥0. Therefore, the ω-limit sets of the solutions of system (2.4), in Ω, are contained in Ω. Furthermore, in Ω, the usual existence, uniqueness, and continuation results hold for the system, so that the system (2.4) is well posed mathematically and epidemiologically [17]. So we consider the dynamics of system (2.4) on the set Ω in this paper.


## 2.2.2. Positivity of Solutions

For system (2.4), it is necessary to prove that all the state variables are positive, so that the solutions of the system with positive initial conditions remain positive for all t≥0. We thus state the following lemma.Lemma 2.1.
Ifs(0)>0,a1(0)>0,a2(0)>0, and r(0)>0, the solutions s(t), a1(t), a2(t), and r(t) of system (2.4) are positive for all t≥0.Proof.
Under the given initial conditions, it is easy to prove that the solutions of system (2.4) are positive; if not, we assume a contradiction: that there exists a first time t1 such that
(2.7)s(t1)=0,s′(t1)<0,a1(t)>0,a2(t)>0,r(t)>0,0<t<t1,
there exists a t2,
(2.8)a1(t2)=0,a1′(t2)<0,s(t)>0,a2(t)>0,r(t)>0,0<t<t2,
there exists a t3,
(2.9)a2(t3)=0,a2′(t3)<0,s(t)>0,a1(t)>0,r(t)>0,0<t<t3,
and there exists a t4,
(2.10)r(t4)=0,r′(t4)<0,s(t)>0,a1(t)>0,a2(t)>0,0<t<t4.
In the first case, we have
(2.11)s′(t1)=μ>0,
which is a contradiction, meaning that s(t)>0, t≥0. In the second case, we have
(2.12)a1′(t2)=(1-p)γa2(t)s(t)>0,
which is a contradiction, meaning that a1(t)>0, t≥0. Similarly, it can be shown that a2(t)>0 and r(t)>0 for all t≥0.Thus, the solutionss(t), a1(t),  a2(t), and r(t) of system (2.4) remain positive for all t>0.


## 3. Analysis of the Model

There are one alcohol-free equilibriumE0 and one alcohol-present equilibrium E* for system (2.4).


### 3.1. Alcohol-Free Equilibrium and the Reproduction Number

The model has an alcohol-free equilibrium given by(3.1)E0=(1,0,0,0).
In the following, the basic reproduction number of system (2.4) will be obtained by the next-generation matrix method formulated in [18].Letx=(a1,a2,r,s)T, then system (2.4) can be written as
(3.2)dxdt=ℱ(x)-𝒱(x),
where
(3.3)ℱ(x)=((1-p)(βa1+γa2)sp(βa1+γa2)s00),𝒱(x)=((α+μ)a1(π+μ)a2-αa1-ηr(η+μ)r-πa2μs+(βa1+γa2)s-μ).
The Jacobian matrices of ℱ(x) and 𝒱(x) at the alcohol-free equilibrium E0 are, respectively,
(3.4)Dℱ(E0)=(F3×3000),D𝒱(E0)=(V3×30ββ0μ),
where
(3.5)F=((1-p)β(1-p)γ0pβpγ0000),V=(α+μ00-απ+μ-η0-πη+μ).
The model reproduction number, denoted by R0, is thus given by
(3.6)R0=ρ(FV-1)=(1-p)βμ(π+η+μ)+(1-p)γα(η+μ)+pγ(α+μ)(η+μ)μ(α+μ)(π+η+μ).
Following Theorem 2 of [18], we have the following result on the local stability of E0.Theorem 3.1.
The alcohol-free equilibriumE0 is locally asymptotically stable for R0<1 and unstable otherwise.


### 3.2. Global Stability ofE0

Theorem 3.2.
For system (2.4), the alcohol-free equilibrium E0 is globally asymptotically stable if R0<1.Proof.
We use the comparison theorem to prove the global stability of the alcohol-free equilibrium. The rate of change of the variables(a1,a2,r) of system (2.4) can be rewritten as
(3.7)(a˙1a˙2r˙)=(F-V)(a1a2r)-(1-s)((1-p)β(1-p)γ0pβpγ0000)(a1a2r),
where F and V are defined in (3.5). Since s≤1 for all t≥0 in Ω, then
(3.8)(a˙1a˙2r˙)≤(F-V)(a1a2r).
Since the eigenvalues of the matrix F-V all have negative real parts (this comes from the local stability results in Lemma 1 in [18]), then system (2.4) is stable whenever R0<1. So, (a1,a2,r)→(0,0,0) as t→∞. By the comparison theorem [19], it follows that (a1,a2,r)→(0,0,0) and s→1 as t→∞. The (s,a1,a2,r)→E0 as t→∞. So, E0 is globally asymptotically stable for R0<1.


### 3.3. Alcohol-Present Equilibrium

#### 3.3.1. Existence of the Alcohol-Present Equilibrium

IfR0>1, system (2.4) has a unique alcohol-present equilibrium E*(s*,a1*,a2*,r*), where
(3.9)s*=1R0,a1*=μ(1-p)α+μ(1-1R0),a2*=p(α+μ)(η+μ)+α(1-p)(η+μ)(1-p)μ(π+η+μ)a1*,r*=πη+μa2*.


#### 3.3.2. Global Stability of the Alcohol-Present Equilibrium

Theorem 3.3.
IfR0>1, the alcohol-present equilibrium E* is globally asymptotically stable.Proof.
To study the global stability of the alcohol-present equilibrium, motivated by [20–22], we use a Lyapunov function V as follows:
(3.10)V=x1(s-s*lns)+x2(a1-a1*lna1)+x3(a2-a2*lna2)+x4(r-r*lnr).
Applying the identity μ=s*(βa1*+γa2*)+μs*, the derivative of V is given by
(3.11)V˙=x1(1-s*s)s˙+x2(1-a1*a1)a˙1+x3(1-a2*a2)a˙2+x4(1-r*r)r˙=x1[μ-s(βa1+γa2)-μs-μs*s+s*(βa1+γa2)+μs*]+x2[(1-p)s(βa1+γa2)-(α+μ)a1-(1-p)βa1sa1*a1-(1-p)γa2sa1*a1+(α+μ)a1*]+x3[ps(βa1+γa2)+ηr+αa1-(π+μ)a2-pβa1sa2*a2-pγa2sa2*a2-ηra2*a2-αa1a2*a2+(π+μ)a2*]+x4[πa2-(η+μ)r-πa2r*r+(η+μ)r*]=x1μs*(2-ss*-s*s)+[x1s*(βa1*+γa2*)+x2(α+μ)a1*+x3(π+μ)a2*+x4(η+μ)r*]-[x1(s*)2(βa1*+γa2*)s+x2(1-p)βsa1*+x2(1-p)γa2sa1*a1+x3pβa1sa2*a2+x3pγsa2*+x3ηra2*a2+x3αa1a2*a2+x4πa2r*r(s*)2(βa1*+γa2*)s]+sa1[-x1β+x2(1-p)β+x3pβ]+sa2[-x1γ+x2(1-p)γ+x3pγ]+a1[x1βs*+x3α-x2(α+μ)]+a2[x1rs*-x3(π+μ)+x4π]+r[x3η-x4(μ+η)].
The positive constantsx1,  x2,  x3, and x4 are chosen such that the coefficients of sa1,  sa2,  a1,  a2, and r are equal to zero, that is,
(3.12)-x1+x2(1-p)+x3p=0,x1βs*-x2(α+μ)+x3α=0,x1γs*-x3(π+μ)+x4π=0,x3η-(η+μ)x4=0.
From (3.12), we have
(3.13)x1=1,x2=μ(π+μ+η)βs*+α(η+μ)γs*μ(α+μ)(η+μ+π),x3=(η+μ)γs*μ(π+η+μ),x4=ηγs*μ(π+η+μ).
We regroup terms in V˙ such that V˙=V˙1+V˙2+V˙3, where
(3.14)V˙1=μs*(2-ss*-s*s),V˙2=x1s*(βa1*+γa2*)+x2(α+μ)a1*+x3(π+μ)a2*+x4(η+μ)r*,V˙3=-x1(s*)2(βa1*+γa2*)s-x2(1-p)βsa1*-x3pγsa2*-x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-x3αa1a2*a2-x3ηra2*a2-x4πa2r*r.
Using the values for x1,  x2,  x3, and x4 in (3.13), relations in (3.12), and the equilibrium relations,
(3.15)(α+pμ)(η+μ)a1*=(1-p)μ(π+η+μ)a2*,x3ηr*=x4πa2*.
We can rewrite V˙2 as
(3.16)V˙2=2x2(1-p)βs*a1*+2x3pγs*a2*+4x3pβs*a1*+3(1-p)αγs*a2*α+pμ+2x4πa2*.
Similarly, we can rewrite V˙3 as
(3.17)V˙3=[-x2(1-p)βsa1*-x2(1-p)β(s*)2a1*s]+[-x3pγsa2*-x3pγ(s*)2a2*s]+[-yx2(1-p)γa2sa1*a1-x3αa1a2*a2-yx2(1-p)γa2*(s*)2s]+[-(1-y)x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-(1-y)x2(1-p)γa2*(s*)2s-x3pβa1*(s*)2s]+[-x3ηra2*a2-x4πa2r*r],
where
(3.18)y=α(α+pμ)x2,1-y=pβs*(α+pμ)x2.
Let V˙3=V˙a+V˙b+V˙c+V˙d+V˙e, with each term representing the expression enclosed in a pair of big square brackets. Using the arithmetic-geometric mean inequality, we obtain
(3.19)V˙a=-x2(1-p)βsa1*-x2(1-p)β(s*)2a1*s≤-2[x22(1-p)2β2(s*)2(a1*)2]1/2=-2x2(1-p)βs*a1*,V˙b=-x3pγsa2*-x3pγ(s*)2a2*s≤-2[x32p2γ2(s*)2(a2*)2]1/2=-2x3pγs*a2*,V˙c=-yx2(1-p)γa2sa1*a1-x3αa1a2*a2-yx2(1-p)γa2*(s*)2s≤-3[x22(1-p)2x3αa1*y2(s*)2γ2(a2*)2]1/3=-3(1-p)αs*γa2*α+pμ,V˙d=-(1-y)x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-(1-y)x2(1-p)γa2*(s*)2s-x3pβa1*(s*)2s≤-4[x22(1-p)2p2x32(1-y)2γ2(a1*)2β2(a2*)2(s*)4]1/4=-4x3pβs*a1*,V˙e=-x3ηra2*a2-x4πa2r*r≤-2[x3ηr*x4πa2*]1/2=-2x4πa2*.
From (3.19), we have
(3.20)V˙3≤-2x2(1-p)βs*a1*-2x3pγs*a2*-4x3pβs*a1*-3(1-p)αγs*a2*α+pμ-2x4πa2*.
From (3.16) and (3.20), we have V˙2+V˙3≤0. Therefore, V˙≤0. Furthermore, V˙=0 if and only if V˙1=0 and V˙2+V˙3=0. We can show that V˙=0⇔(s,a1,a2,r)=(s*,a1*,a2*,r*). Thus, V˙ is a negative definite with respect to P*. So the alcohol-present equilibrium E* is globally asymptotically stable.


## 3.1. Alcohol-Free Equilibrium and the Reproduction Number

The model has an alcohol-free equilibrium given by(3.1)E0=(1,0,0,0).
In the following, the basic reproduction number of system (2.4) will be obtained by the next-generation matrix method formulated in [18].Letx=(a1,a2,r,s)T, then system (2.4) can be written as
(3.2)dxdt=ℱ(x)-𝒱(x),
where
(3.3)ℱ(x)=((1-p)(βa1+γa2)sp(βa1+γa2)s00),𝒱(x)=((α+μ)a1(π+μ)a2-αa1-ηr(η+μ)r-πa2μs+(βa1+γa2)s-μ).
The Jacobian matrices of ℱ(x) and 𝒱(x) at the alcohol-free equilibrium E0 are, respectively,
(3.4)Dℱ(E0)=(F3×3000),D𝒱(E0)=(V3×30ββ0μ),
where
(3.5)F=((1-p)β(1-p)γ0pβpγ0000),V=(α+μ00-απ+μ-η0-πη+μ).
The model reproduction number, denoted by R0, is thus given by
(3.6)R0=ρ(FV-1)=(1-p)βμ(π+η+μ)+(1-p)γα(η+μ)+pγ(α+μ)(η+μ)μ(α+μ)(π+η+μ).
Following Theorem 2 of [18], we have the following result on the local stability of E0.Theorem 3.1.
The alcohol-free equilibriumE0 is locally asymptotically stable for R0<1 and unstable otherwise.


## 3.2. Global Stability ofE0

Theorem 3.2.
For system (2.4), the alcohol-free equilibrium E0 is globally asymptotically stable if R0<1.Proof.
We use the comparison theorem to prove the global stability of the alcohol-free equilibrium. The rate of change of the variables(a1,a2,r) of system (2.4) can be rewritten as
(3.7)(a˙1a˙2r˙)=(F-V)(a1a2r)-(1-s)((1-p)β(1-p)γ0pβpγ0000)(a1a2r),
where F and V are defined in (3.5). Since s≤1 for all t≥0 in Ω, then
(3.8)(a˙1a˙2r˙)≤(F-V)(a1a2r).
Since the eigenvalues of the matrix F-V all have negative real parts (this comes from the local stability results in Lemma 1 in [18]), then system (2.4) is stable whenever R0<1. So, (a1,a2,r)→(0,0,0) as t→∞. By the comparison theorem [19], it follows that (a1,a2,r)→(0,0,0) and s→1 as t→∞. The (s,a1,a2,r)→E0 as t→∞. So, E0 is globally asymptotically stable for R0<1.


## 3.3. Alcohol-Present Equilibrium

### 3.3.1. Existence of the Alcohol-Present Equilibrium

IfR0>1, system (2.4) has a unique alcohol-present equilibrium E*(s*,a1*,a2*,r*), where
(3.9)s*=1R0,a1*=μ(1-p)α+μ(1-1R0),a2*=p(α+μ)(η+μ)+α(1-p)(η+μ)(1-p)μ(π+η+μ)a1*,r*=πη+μa2*.


### 3.3.2. Global Stability of the Alcohol-Present Equilibrium

Theorem 3.3.
IfR0>1, the alcohol-present equilibrium E* is globally asymptotically stable.Proof.
To study the global stability of the alcohol-present equilibrium, motivated by [20–22], we use a Lyapunov function V as follows:
(3.10)V=x1(s-s*lns)+x2(a1-a1*lna1)+x3(a2-a2*lna2)+x4(r-r*lnr).
Applying the identity μ=s*(βa1*+γa2*)+μs*, the derivative of V is given by
(3.11)V˙=x1(1-s*s)s˙+x2(1-a1*a1)a˙1+x3(1-a2*a2)a˙2+x4(1-r*r)r˙=x1[μ-s(βa1+γa2)-μs-μs*s+s*(βa1+γa2)+μs*]+x2[(1-p)s(βa1+γa2)-(α+μ)a1-(1-p)βa1sa1*a1-(1-p)γa2sa1*a1+(α+μ)a1*]+x3[ps(βa1+γa2)+ηr+αa1-(π+μ)a2-pβa1sa2*a2-pγa2sa2*a2-ηra2*a2-αa1a2*a2+(π+μ)a2*]+x4[πa2-(η+μ)r-πa2r*r+(η+μ)r*]=x1μs*(2-ss*-s*s)+[x1s*(βa1*+γa2*)+x2(α+μ)a1*+x3(π+μ)a2*+x4(η+μ)r*]-[x1(s*)2(βa1*+γa2*)s+x2(1-p)βsa1*+x2(1-p)γa2sa1*a1+x3pβa1sa2*a2+x3pγsa2*+x3ηra2*a2+x3αa1a2*a2+x4πa2r*r(s*)2(βa1*+γa2*)s]+sa1[-x1β+x2(1-p)β+x3pβ]+sa2[-x1γ+x2(1-p)γ+x3pγ]+a1[x1βs*+x3α-x2(α+μ)]+a2[x1rs*-x3(π+μ)+x4π]+r[x3η-x4(μ+η)].
The positive constantsx1,  x2,  x3, and x4 are chosen such that the coefficients of sa1,  sa2,  a1,  a2, and r are equal to zero, that is,
(3.12)-x1+x2(1-p)+x3p=0,x1βs*-x2(α+μ)+x3α=0,x1γs*-x3(π+μ)+x4π=0,x3η-(η+μ)x4=0.
From (3.12), we have
(3.13)x1=1,x2=μ(π+μ+η)βs*+α(η+μ)γs*μ(α+μ)(η+μ+π),x3=(η+μ)γs*μ(π+η+μ),x4=ηγs*μ(π+η+μ).
We regroup terms in V˙ such that V˙=V˙1+V˙2+V˙3, where
(3.14)V˙1=μs*(2-ss*-s*s),V˙2=x1s*(βa1*+γa2*)+x2(α+μ)a1*+x3(π+μ)a2*+x4(η+μ)r*,V˙3=-x1(s*)2(βa1*+γa2*)s-x2(1-p)βsa1*-x3pγsa2*-x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-x3αa1a2*a2-x3ηra2*a2-x4πa2r*r.
Using the values for x1,  x2,  x3, and x4 in (3.13), relations in (3.12), and the equilibrium relations,
(3.15)(α+pμ)(η+μ)a1*=(1-p)μ(π+η+μ)a2*,x3ηr*=x4πa2*.
We can rewrite V˙2 as
(3.16)V˙2=2x2(1-p)βs*a1*+2x3pγs*a2*+4x3pβs*a1*+3(1-p)αγs*a2*α+pμ+2x4πa2*.
Similarly, we can rewrite V˙3 as
(3.17)V˙3=[-x2(1-p)βsa1*-x2(1-p)β(s*)2a1*s]+[-x3pγsa2*-x3pγ(s*)2a2*s]+[-yx2(1-p)γa2sa1*a1-x3αa1a2*a2-yx2(1-p)γa2*(s*)2s]+[-(1-y)x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-(1-y)x2(1-p)γa2*(s*)2s-x3pβa1*(s*)2s]+[-x3ηra2*a2-x4πa2r*r],
where
(3.18)y=α(α+pμ)x2,1-y=pβs*(α+pμ)x2.
Let V˙3=V˙a+V˙b+V˙c+V˙d+V˙e, with each term representing the expression enclosed in a pair of big square brackets. Using the arithmetic-geometric mean inequality, we obtain
(3.19)V˙a=-x2(1-p)βsa1*-x2(1-p)β(s*)2a1*s≤-2[x22(1-p)2β2(s*)2(a1*)2]1/2=-2x2(1-p)βs*a1*,V˙b=-x3pγsa2*-x3pγ(s*)2a2*s≤-2[x32p2γ2(s*)2(a2*)2]1/2=-2x3pγs*a2*,V˙c=-yx2(1-p)γa2sa1*a1-x3αa1a2*a2-yx2(1-p)γa2*(s*)2s≤-3[x22(1-p)2x3αa1*y2(s*)2γ2(a2*)2]1/3=-3(1-p)αs*γa2*α+pμ,V˙d=-(1-y)x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-(1-y)x2(1-p)γa2*(s*)2s-x3pβa1*(s*)2s≤-4[x22(1-p)2p2x32(1-y)2γ2(a1*)2β2(a2*)2(s*)4]1/4=-4x3pβs*a1*,V˙e=-x3ηra2*a2-x4πa2r*r≤-2[x3ηr*x4πa2*]1/2=-2x4πa2*.
From (3.19), we have
(3.20)V˙3≤-2x2(1-p)βs*a1*-2x3pγs*a2*-4x3pβs*a1*-3(1-p)αγs*a2*α+pμ-2x4πa2*.
From (3.16) and (3.20), we have V˙2+V˙3≤0. Therefore, V˙≤0. Furthermore, V˙=0 if and only if V˙1=0 and V˙2+V˙3=0. We can show that V˙=0⇔(s,a1,a2,r)=(s*,a1*,a2*,r*). Thus, V˙ is a negative definite with respect to P*. So the alcohol-present equilibrium E* is globally asymptotically stable.


## 3.3.1. Existence of the Alcohol-Present Equilibrium

IfR0>1, system (2.4) has a unique alcohol-present equilibrium E*(s*,a1*,a2*,r*), where
(3.9)s*=1R0,a1*=μ(1-p)α+μ(1-1R0),a2*=p(α+μ)(η+μ)+α(1-p)(η+μ)(1-p)μ(π+η+μ)a1*,r*=πη+μa2*.


## 3.3.2. Global Stability of the Alcohol-Present Equilibrium

Theorem 3.3.
IfR0>1, the alcohol-present equilibrium E* is globally asymptotically stable.Proof.
To study the global stability of the alcohol-present equilibrium, motivated by [20–22], we use a Lyapunov function V as follows:
(3.10)V=x1(s-s*lns)+x2(a1-a1*lna1)+x3(a2-a2*lna2)+x4(r-r*lnr).
Applying the identity μ=s*(βa1*+γa2*)+μs*, the derivative of V is given by
(3.11)V˙=x1(1-s*s)s˙+x2(1-a1*a1)a˙1+x3(1-a2*a2)a˙2+x4(1-r*r)r˙=x1[μ-s(βa1+γa2)-μs-μs*s+s*(βa1+γa2)+μs*]+x2[(1-p)s(βa1+γa2)-(α+μ)a1-(1-p)βa1sa1*a1-(1-p)γa2sa1*a1+(α+μ)a1*]+x3[ps(βa1+γa2)+ηr+αa1-(π+μ)a2-pβa1sa2*a2-pγa2sa2*a2-ηra2*a2-αa1a2*a2+(π+μ)a2*]+x4[πa2-(η+μ)r-πa2r*r+(η+μ)r*]=x1μs*(2-ss*-s*s)+[x1s*(βa1*+γa2*)+x2(α+μ)a1*+x3(π+μ)a2*+x4(η+μ)r*]-[x1(s*)2(βa1*+γa2*)s+x2(1-p)βsa1*+x2(1-p)γa2sa1*a1+x3pβa1sa2*a2+x3pγsa2*+x3ηra2*a2+x3αa1a2*a2+x4πa2r*r(s*)2(βa1*+γa2*)s]+sa1[-x1β+x2(1-p)β+x3pβ]+sa2[-x1γ+x2(1-p)γ+x3pγ]+a1[x1βs*+x3α-x2(α+μ)]+a2[x1rs*-x3(π+μ)+x4π]+r[x3η-x4(μ+η)].
The positive constantsx1,  x2,  x3, and x4 are chosen such that the coefficients of sa1,  sa2,  a1,  a2, and r are equal to zero, that is,
(3.12)-x1+x2(1-p)+x3p=0,x1βs*-x2(α+μ)+x3α=0,x1γs*-x3(π+μ)+x4π=0,x3η-(η+μ)x4=0.
From (3.12), we have
(3.13)x1=1,x2=μ(π+μ+η)βs*+α(η+μ)γs*μ(α+μ)(η+μ+π),x3=(η+μ)γs*μ(π+η+μ),x4=ηγs*μ(π+η+μ).
We regroup terms in V˙ such that V˙=V˙1+V˙2+V˙3, where
(3.14)V˙1=μs*(2-ss*-s*s),V˙2=x1s*(βa1*+γa2*)+x2(α+μ)a1*+x3(π+μ)a2*+x4(η+μ)r*,V˙3=-x1(s*)2(βa1*+γa2*)s-x2(1-p)βsa1*-x3pγsa2*-x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-x3αa1a2*a2-x3ηra2*a2-x4πa2r*r.
Using the values for x1,  x2,  x3, and x4 in (3.13), relations in (3.12), and the equilibrium relations,
(3.15)(α+pμ)(η+μ)a1*=(1-p)μ(π+η+μ)a2*,x3ηr*=x4πa2*.
We can rewrite V˙2 as
(3.16)V˙2=2x2(1-p)βs*a1*+2x3pγs*a2*+4x3pβs*a1*+3(1-p)αγs*a2*α+pμ+2x4πa2*.
Similarly, we can rewrite V˙3 as
(3.17)V˙3=[-x2(1-p)βsa1*-x2(1-p)β(s*)2a1*s]+[-x3pγsa2*-x3pγ(s*)2a2*s]+[-yx2(1-p)γa2sa1*a1-x3αa1a2*a2-yx2(1-p)γa2*(s*)2s]+[-(1-y)x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-(1-y)x2(1-p)γa2*(s*)2s-x3pβa1*(s*)2s]+[-x3ηra2*a2-x4πa2r*r],
where
(3.18)y=α(α+pμ)x2,1-y=pβs*(α+pμ)x2.
Let V˙3=V˙a+V˙b+V˙c+V˙d+V˙e, with each term representing the expression enclosed in a pair of big square brackets. Using the arithmetic-geometric mean inequality, we obtain
(3.19)V˙a=-x2(1-p)βsa1*-x2(1-p)β(s*)2a1*s≤-2[x22(1-p)2β2(s*)2(a1*)2]1/2=-2x2(1-p)βs*a1*,V˙b=-x3pγsa2*-x3pγ(s*)2a2*s≤-2[x32p2γ2(s*)2(a2*)2]1/2=-2x3pγs*a2*,V˙c=-yx2(1-p)γa2sa1*a1-x3αa1a2*a2-yx2(1-p)γa2*(s*)2s≤-3[x22(1-p)2x3αa1*y2(s*)2γ2(a2*)2]1/3=-3(1-p)αs*γa2*α+pμ,V˙d=-(1-y)x2(1-p)γa2sa1*a1-x3pβa1sa2*a2-(1-y)x2(1-p)γa2*(s*)2s-x3pβa1*(s*)2s≤-4[x22(1-p)2p2x32(1-y)2γ2(a1*)2β2(a2*)2(s*)4]1/4=-4x3pβs*a1*,V˙e=-x3ηra2*a2-x4πa2r*r≤-2[x3ηr*x4πa2*]1/2=-2x4πa2*.
From (3.19), we have
(3.20)V˙3≤-2x2(1-p)βs*a1*-2x3pγs*a2*-4x3pβs*a1*-3(1-p)αγs*a2*α+pμ-2x4πa2*.
From (3.16) and (3.20), we have V˙2+V˙3≤0. Therefore, V˙≤0. Furthermore, V˙=0 if and only if V˙1=0 and V˙2+V˙3=0. We can show that V˙=0⇔(s,a1,a2,r)=(s*,a1*,a2*,r*). Thus, V˙ is a negative definite with respect to P*. So the alcohol-present equilibrium E* is globally asymptotically stable.


## 4. Numerical Simulation

To illustrate the analytic results obtained above, we give some simulations using the parameter values in Table1. Numerical results are displayed in the following figures. First, we choose p=0.8, α=0.2; numerical simulation gives R0=0.4091<1, then Theorem 3.2 indicates that youths do not have alcohol problems. Figure 2 confirms this conclusion. Second, we choose p=0.1,  α=0.4; numerical simulation gives R0=1.8031>1; then Theorem 3.3 shows that alcohol problems persist. Figure 3 further validates the conclusion. Finally, we choose α=0.1; numerical simulation gives the relation between p and R0. Figure 4 confirms that R0 shows a decline, while the number of admitting drinkers increases. From the figures above, we find that making more people admit having alcohol problems can reduce the alcohol problems. Then, the treatment of admitting drinkers is an effective measure in alcohol problems.Table 1
Description and estimation of parameters.
Parameter
Description
Estimated value
Date source
μ
Recruitment rate of the population
0.25 year−1
[15]
β
Transmission coefficient of the nonadmitting compartment
0.3 year−1
[15]
γ
Transmission coefficient of the admitting compartment
0.25 year−1
Estimate
α
The fraction ofS being infected by A1 and entering A2
Variable
π
The fraction ofA2 going into treatment
0.05 year−1
[15]
η
The fraction ofR who relapse into A2
0.8 year−1
[15]
p
Probability of a newly infected individual who is admitting the problem
Variable
1
-
p
Probability of a newly infected individual who is not admitting the problem
VariableWhenp=0.8, α=0.2, and  R0=0.4091<1, the alcohol-free equilibrium E0 is globally asymptotically stable.Whenp=0.1,  α=0.4, and R0=1.8031>1, the alcohol-present equilibrium E* is globally asymptotically stable.Figure 4
Whenp increased, R0 decline.


## 5. Discussion

We have formulated a binge drinking model with two stages and investigated their dynamical behaviors. Depending on the basic reproduction numberR0, the steady state is either the alcohol free or the alcohol-present. By using the comparison theorem, we prove that all solutions converge to E0 when the basic reproduction number is less than one, that is, the alcohol problems disappear eventually. By constructing the Lyapunov function, we prove that the unique alcohol-present equilibrium is globally stable, that is, the alcohol problems will persist in the population, and the number of binge drinking individuals tends to be a positive constant when the basic reproduction number exceeds one.To better understand the binge drinking model with two stages, we consider the more realistic model and investigate the global stability of equilibria of the model for binge drinking problem. The global stability is not discussed in the literature [15].The basic reproductive number in this paper is given by(5.1)R0=(1-p)βμ(π+η+μ)+(1-p)γα(η+μ)+pγ(α+μ)(η+μ)μ(α+μ)(π+η+μ).If the model has only one transfer from the class of susceptible individuals towards the class of nonadmitting drinkers(p=0), then the basic reproductive number is
(5.2)R1=βμ(π+η+μ)+γα(η+μ)μ(α+μ)(π+η+μ).
It is similar to the basic reproductive number in [15]. Comparing with the two expressions, we see that R0 is a function on p. In Figure 4, we find that the basic reproductive number R0 decreases when the number of admitting drinkers increases. In general, the nonadmitting drinkers are also unlikely to seek help. Making more people admit to have alcohol problems through the awareness programs, such as education and media programs is an effective measure to control the alcohol problems.In deriving our model (2.2), we consider the direct transfer from the class of susceptible individuals towards the class of admitting drinkers and assume the relapse from the class of treatment individuals towards the class of admitting drinkers. However, if we hope to include another relapse from the class of treatment individuals towards the class of nonadmitting drinkers, we could modify (2.2) to the following model:
(5.3)S˙=μN-S(βA1+γA2)N-μS,A˙1=(1-p)S(βA1+γA2)N+ξR-(α+μ)A1,A˙2=pS(βA1+γA2)N+ηR+αA1-(π+μ)A2,R˙=πA2-(μ+η+ξ)R,
where ξ represents relapse from the class of treatment individuals towards the class of nonadmitting drinkers.In another way, we might consider the class of quit drinkers and also consider the following model:(5.4)S˙=μN-S(βA1+γA2)N-μS,A˙1=(1-p)S(βA1+γA2)N+ξR-(α+μ)A1,A˙2=pS(βA1+γA2)N+ηR+αA1-(π+μ)A2,R˙=πA2-(μ+η+ξ+δ)R,Q˙=δR-μQ,
where Q(t) represents the population that has stopped drinking permanently at time t, say quit drinkers. The δ term is the rate from the class of treatment individuals towards the class of quit drinkers.  We leave these works for the future.



---

*Source: 829386-2012-11-11.xml*