# Stability Analysis of an Alcoholism Model with Public Health Education and NSFD Scheme

Authors: Zhaofeng An; Suxia Zhang; Jinhu Xu
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In this paper, an alcoholism model of SEAR type with different susceptibilities due to public health education is investigated, with the form of continuous differential equations as well as discrete differential equations by applying the Mickens nonstandard finite difference (NSFD) scheme to the continuous equations. Threshold dynamics of the continuous model are performed by constructing Lyapunov functions. The analysis of a discrete model indicates that the alcohol-free equilibrium is globally asymptotically stable if the basic reproductive numberR0<1, and conversely, the alcohol-present equilibrium is globally asymptotically stable if R0>1, revealing the consistency and efficiency of the discrete model to preserve the dynamical properties of the corresponding continuous model. In addition, stability preserving and the impact of the parameters related with public health education are conducted by numerical simulations.


---

## Body

## 1. Introduction

Alcohol abuse is a major public health concern due to its potential negative health and social effect, especially among college students. There was evidence showing that more than40% of college students have experienced alcoholism [1, 2], and it has been observed that binge drinking and drinking behaviors have caused negative effects in different levels on students’ academic study, social relations, risk-taking, and health [3, 4]. Particularly, alcohol-related injuries and deaths per 100,000 college students increased by 6% from 1998 to 2001 [5]. Studies have shown that alcoholics are more likely to engage in risky behaviors that may result in serious consequences, such as sexually transmitted infections (STIs) [6]. Although much attention has been paid to the reduction in initiation and maintenance of alcohol consumption, the management of such problems in youths still faces a challenge in many countries.In an attempt to address some of the issues encountered in alcohol control and gain insights into drinking dynamics, many theoretical studies have been carried out by formulating mathematical models that focus on different issues [6–11]. According to the fact that social factors play an important role in modelling the spread of diseases, a number of studies have captured the impact of social interactions on the dynamics of alcohol drinking [8, 12–17]. Benedict presented an SAR-type model with standard incidence to describe alcoholism and studied the existence and stability of the equilibria [12]. In consideration of low-risk and high-risk social transmission process, the authors found that the social interaction between light drinkers and moderate drinkers in a low-risk environment may weaken the significance of the influence of heavy drinking on the disease [16]. Mulone and Straughan in [8] investigated a model for binge drinking in particular with two types of alcoholics and studied the threshold dynamics, which was extended by including the possibility of complete recovery from alcohol dependence in [9]. Different from the assumption of equal infectivity in [8], Huo and Song made a modification of taking into account the different infectious risks of alcoholics in the modelling [14]. In view of alcohol abuse among students on and around college campuses, there have been compartmental models established to study the role of environmental factors on the dynamics of this specific drinking group [13, 15]. Other social behaviors of individuals in the epidemic population, such as treatment, migration, and public health education campaign, have been also incorporated in models to analyze epidemics [18, 19] and alcohol drinking [10, 20–22].In the literature, different approaches to modelling alcoholism have been used, including continuous ordinary differential equation models [7, 10, 20], partial differential equation models [11, 23], stochastic differential equation model [13], and network models [24, 25]. Meanwhile, the discrete model, with its obvious advantage of describing epidemiology [26–28], can also be employed to understand the process of drinking, which is rarely seen in the abovementioned references. Hence, based on the framework of the SAR-type model, we construct two models with the form of continuous differential equation and discrete differential equation, respectively. In terms of the discretization, the Mickens method is adopted for its efficiency in preserving the global asymptotic stability of equilibria to the corresponding continuous model [29, 30], which has attracted much attention [31–34]. On the contrary, it is known that the function form of the incidence rate in the epidemic model has a crucial role in determining the qualitative behavior of the proposed model and in giving a reasonable description of the dynamics. To account for the influence of public health education on the transmission of alcohol consumption, a general nonlinear incidence is utilized in our models, which is different from that in [20].The organization of this paper is as follows. In Section2, we construct and analyze the continuous model. The corresponding discrete model by applying Mickens nonstandard finite difference method is established and conducted in Section 3. In Section 4, numerical simulations are carried out to illustrate the stability preserving and the effect of educational campaigns, and a brief conclusion is given in Section 5.


## 2. Continuous Drinking Model

### 2.1. Model Formulation

According to the drinking status, we divide the host population into four categories, i.e., susceptible drinkersSt, who drink moderately and are susceptible to becoming alcoholic without public health education; educated drinkers Et, who are susceptible drinkers with public health education; alcoholics At, who have drinking problems or addictions; and nondrinkers Rt, who quit drinking or never have drinking problems. Let Nt denote the total number of all individuals at time t, i.e., Nt=St+Et+At+Rt. In the modelling, we make some assumptions as follows (see [20]):(i)
A1: the recruitment is given by a constant fluxΛ, which can enter into the nondrinker class R with proportion of q and the rest 1−qΛ into the susceptible class S(ii)
A2: the susceptible individuals inS can move into the educated class E due to their knowledge by public health education, and individuals in these two classes can be both affected by alcoholics, but with different rates(iii)
A3: alcoholics inA can become nondrinkers because of some factors, such as counselling, treatment, and prohibitionWith these assumptions, we have the following system of continuous ordinary differential equations:(1)dStdt=1−qΛ−βStfAt−μS+pSt,dEtdt=pSt−σβEtfAt−μE+εEt,dAtdt=St+σEtβfAt−μA+γAt,dRtdt=qΛ+γAt+εEt−μRRt.In this model, parameterσ is used to denote the infectiousness of educated individuals relative to susceptible ones, which reveals the effectiveness of public health education. σ=0 implies that the education is completely effective and σ=1 accounts for its useless. β is the transmission coefficient. p denotes the rate moving from the susceptible to the educated class. The individuals in the educated class move to the nondrinker class at a rate ε, and the alcoholics can recover at a rate γ to the nondrinker class. μS, μE, μA, and μR denote the death rates in the susceptible, educated, alcoholic, and nondrinker class, respectively. In model (1), to give the incidence rate, we employ a general function fA, which is sufficiently smooth and is assumed to satisfy the following conditions:(iv)
A4:f0=0 for A=0 and f′A>0,f″A≤0 for A≥0.The assumptions forfA are fundamental and biologically motivated, and it is easy to check that the function satisfying the above conditions includes the bilinear and the saturation incidences. Note that the variable R in model (1) does not appear in other equations, and thus, the R equation will be neglected in the remaining analysis. The dynamics of the following system without the R equation would be the same:(2)dStdt=1−qΛ−βStfAt−μS+pSt,dEtdt=pSt−σβEtfAt−μE+εEt,dAtdt=St+σEtβfAt−μA+γAt.All the parameters in system (2) are nonnegative, and we are interested in those nonnegative solutions. It is easy to show that the solution of system (2) with nonnegative initial values is nonnegative and exists for all t≥0. From the first equation of (2), we have(3)limt⟶∞supSt≤1−qΛμS+p.Furthermore,(4)S+E+A′≤1−qΛ−μSt+Et+At,with μ=minμS,μE+ε,μA+γ, and we can study the dynamical behavior of system (2) in the positively invariant set:(5)Ω≔S,E,A∈R+3:0≤S≤1−qΛμS+p,0≤S+E+A≤1−qΛμ.It is easy to see that model (2) always has a alcohol-free equilibrium P0=S0,E0,0 with S0=1−qΛ/μS+p and E0=p1−qΛ/μE+εμS+p. Define the basic reproductive number as(6)R0=βf′0S0+σE0μA+γ,which gives the expected number of secondary cases produced in a completely susceptible and educated population by a typical alcoholics during the time he spends in the infectious period. Then, the alcohol-present equilibrium P∗=S∗,E∗,A∗ exists if R0>1, where(7)S∗=1−qΛβfA∗+μS+p,E∗=p1−qΛβfA∗+μS+pσβfA∗+μE+ε.Indeed, the alcohol-present equilibrium satisfies the following equations:(8)1−qΛ=βS∗fA∗+μS+pS∗,pS∗=σβE∗fA∗+μE+εE∗,μA+γA∗=S∗+σE∗βfA∗.The expression ofS∗ can be directly obtained by the first equality and E∗ by the second in (8), and then it follows from the last equality that A∗ remains positive for all R0>1.


### 2.2. Stability Analysis

In this subsection, the globally asymptotic stabilities ofP0 and P∗ will be discussed by constructing suitable Lyapunov functions.Theorem 1.
IfR0<1, then the alcohol-free equilibrium P0 is globally asymptotically stable.Proof.
Firstly, we use the Hurwitz criterion to discuss the local stability ofP0. The characteristic equation at P0 is(9)λ+μS+p0βS0f′0−pλ+μE+εσβE0f′000λ+μA+γ−S0+σE0βf′0=0,which has three eigenvalues(10)λ1=−μS+p,λ2=−μE+ε,λ3=S0+σE0βf′0−μA+γ.It is clear that λ3<0 when R0<1, then P0 is locally stable.
Define a Lyapunov function:(11)LP0=1μA+γS−S0−S0lnSS0+E−E0−E0lnEE0+A.
Calculating the derivative ofLP0 along the solution of system (2), we can obtain(12)LP0′=1μA+a+γ1−S0S1−qΛ−βSfA−μS+pS+1−E0EpS−σβEfA−μE+εE+S+σEβfA−μA+γA.
Note that1−qΛ=μS+pS0 and pS0=μE+εE0, then LP0′ can be written as(13)LP0′=1μA+γμS+pS02−S0S−SS0+μE+εE02−E0E−EE0+pS0SS0−1−E0SES0+EE0+βS0fA+σβE0fA−μA+γA=1μA+γμSS02−S0S−SS0+μE+εE03−S0S−EE0−SE0S0E+βSfA1−S0S+βS0fA+σβE0fA−μA+γA≤1μA+γμSS02−S0S−SS0+μE+εE03−S0S−EE0−SE0S0E+1μA+γβS0f′0A+σβE0f′0A−μA+γA=1μA+γμSS02−S0S−SS0+μE+εE03−S0S−EE0−SE0S0E+R0−1A.
Notice the fact(14)2−S0S−SS0≤0,3−S0S−EE0−SE0S0E≤0,for all S,E>0 and the equalities hold if and only if S=S0 and E=E0. Hence, we have LP0′≤0 when R0<1 and LP0′=0 only at P0, which leads to the global stability of P0. This completes the proof.Theorem 2.
IfR0>1, the alcohol-present equilibrium P∗ is globally asymptotically stable.Proof.
Construct a Lyapunov function as follows:(15)LP∗=S−S∗−S∗lnSS∗+E−E∗−E∗lnEE∗+A−A∗−A∗lnAA∗.Then, the time derivative of LP∗ along the solution of (2) satisfies(16)LP∗′=1−S∗S1−qΛ−βSfA−μS+pS+1−E∗EpS−σβEfA−μE+εE+1−A∗AS+σEβfA−μA+γA.
Noticing that(17)1−qΛ=βS∗fA∗+μS+pS∗,pS∗=μE+εE∗+σβE∗fA∗,we have(18)LP∗′=μS+pS∗2−S∗S−SS∗+1−S∗SβS∗fA∗−βSfA+μE+εE∗2−E∗E−EE∗+pS∗SS∗−1−E∗SES∗+E∗E+1−E∗EσβE∗fA∗−σβEfA+1−A∗AβSfA+σβEfA−βS∗fA∗+σβE∗fA∗AA∗=μSS∗2−S∗S−SS∗+μE+εE∗3−S∗S−EE∗−SE∗S∗E+βS∗fA∗2−S∗S−SfAA∗S∗fA∗A∗+fAfA∗−AA∗+σβE∗fA∗3−S∗S−SE∗S∗E−EfAA∗E∗fA∗A+fAfA∗−AA∗=μSS∗2−S∗S−SS∗+μE+εE∗3−S∗S−EE∗−SE∗S∗E+βS∗fA∗−gS∗S−gSfAA∗S∗fA∗A+fAfA∗−AA∗+σβE∗fA∗−gS∗S−gSE∗S∗E−gEfAA∗E∗fA∗A+fAfA∗−AA∗=μSS∗2−S∗S−SS∗+μE+εE∗3−S∗S−EE∗−SE∗S∗E+βS∗fA∗−gS∗S−gSfAA∗S∗fA∗A−gAfA∗A∗fA+fAfA∗−AA∗1−fA∗fA+σβE∗fA∗−gS∗S−gSE∗S∗E−gEfAA∗E∗fA∗A−gAfA∗A∗fA+fAfA∗−AA∗1−fA∗fA,where gx=x−1−lnx≥0 for all x>0. Furthermore, the assumption of fA implies that(19)fAfA∗−AA∗1−fA∗fA≤0,which leads to LP∗′≤0 and the largest compact invariant set of LP∗′≤0 is the singleton P∗. Therefore, by LaSalle’s invariance principle, P∗ is globally asymptotically stable when R0>1.


## 2.1. Model Formulation

According to the drinking status, we divide the host population into four categories, i.e., susceptible drinkersSt, who drink moderately and are susceptible to becoming alcoholic without public health education; educated drinkers Et, who are susceptible drinkers with public health education; alcoholics At, who have drinking problems or addictions; and nondrinkers Rt, who quit drinking or never have drinking problems. Let Nt denote the total number of all individuals at time t, i.e., Nt=St+Et+At+Rt. In the modelling, we make some assumptions as follows (see [20]):(i)
A1: the recruitment is given by a constant fluxΛ, which can enter into the nondrinker class R with proportion of q and the rest 1−qΛ into the susceptible class S(ii)
A2: the susceptible individuals inS can move into the educated class E due to their knowledge by public health education, and individuals in these two classes can be both affected by alcoholics, but with different rates(iii)
A3: alcoholics inA can become nondrinkers because of some factors, such as counselling, treatment, and prohibitionWith these assumptions, we have the following system of continuous ordinary differential equations:(1)dStdt=1−qΛ−βStfAt−μS+pSt,dEtdt=pSt−σβEtfAt−μE+εEt,dAtdt=St+σEtβfAt−μA+γAt,dRtdt=qΛ+γAt+εEt−μRRt.In this model, parameterσ is used to denote the infectiousness of educated individuals relative to susceptible ones, which reveals the effectiveness of public health education. σ=0 implies that the education is completely effective and σ=1 accounts for its useless. β is the transmission coefficient. p denotes the rate moving from the susceptible to the educated class. The individuals in the educated class move to the nondrinker class at a rate ε, and the alcoholics can recover at a rate γ to the nondrinker class. μS, μE, μA, and μR denote the death rates in the susceptible, educated, alcoholic, and nondrinker class, respectively. In model (1), to give the incidence rate, we employ a general function fA, which is sufficiently smooth and is assumed to satisfy the following conditions:(iv)
A4:f0=0 for A=0 and f′A>0,f″A≤0 for A≥0.The assumptions forfA are fundamental and biologically motivated, and it is easy to check that the function satisfying the above conditions includes the bilinear and the saturation incidences. Note that the variable R in model (1) does not appear in other equations, and thus, the R equation will be neglected in the remaining analysis. The dynamics of the following system without the R equation would be the same:(2)dStdt=1−qΛ−βStfAt−μS+pSt,dEtdt=pSt−σβEtfAt−μE+εEt,dAtdt=St+σEtβfAt−μA+γAt.All the parameters in system (2) are nonnegative, and we are interested in those nonnegative solutions. It is easy to show that the solution of system (2) with nonnegative initial values is nonnegative and exists for all t≥0. From the first equation of (2), we have(3)limt⟶∞supSt≤1−qΛμS+p.Furthermore,(4)S+E+A′≤1−qΛ−μSt+Et+At,with μ=minμS,μE+ε,μA+γ, and we can study the dynamical behavior of system (2) in the positively invariant set:(5)Ω≔S,E,A∈R+3:0≤S≤1−qΛμS+p,0≤S+E+A≤1−qΛμ.It is easy to see that model (2) always has a alcohol-free equilibrium P0=S0,E0,0 with S0=1−qΛ/μS+p and E0=p1−qΛ/μE+εμS+p. Define the basic reproductive number as(6)R0=βf′0S0+σE0μA+γ,which gives the expected number of secondary cases produced in a completely susceptible and educated population by a typical alcoholics during the time he spends in the infectious period. Then, the alcohol-present equilibrium P∗=S∗,E∗,A∗ exists if R0>1, where(7)S∗=1−qΛβfA∗+μS+p,E∗=p1−qΛβfA∗+μS+pσβfA∗+μE+ε.Indeed, the alcohol-present equilibrium satisfies the following equations:(8)1−qΛ=βS∗fA∗+μS+pS∗,pS∗=σβE∗fA∗+μE+εE∗,μA+γA∗=S∗+σE∗βfA∗.The expression ofS∗ can be directly obtained by the first equality and E∗ by the second in (8), and then it follows from the last equality that A∗ remains positive for all R0>1.


## 2.2. Stability Analysis

In this subsection, the globally asymptotic stabilities ofP0 and P∗ will be discussed by constructing suitable Lyapunov functions.Theorem 1.
IfR0<1, then the alcohol-free equilibrium P0 is globally asymptotically stable.Proof.
Firstly, we use the Hurwitz criterion to discuss the local stability ofP0. The characteristic equation at P0 is(9)λ+μS+p0βS0f′0−pλ+μE+εσβE0f′000λ+μA+γ−S0+σE0βf′0=0,which has three eigenvalues(10)λ1=−μS+p,λ2=−μE+ε,λ3=S0+σE0βf′0−μA+γ.It is clear that λ3<0 when R0<1, then P0 is locally stable.
Define a Lyapunov function:(11)LP0=1μA+γS−S0−S0lnSS0+E−E0−E0lnEE0+A.
Calculating the derivative ofLP0 along the solution of system (2), we can obtain(12)LP0′=1μA+a+γ1−S0S1−qΛ−βSfA−μS+pS+1−E0EpS−σβEfA−μE+εE+S+σEβfA−μA+γA.
Note that1−qΛ=μS+pS0 and pS0=μE+εE0, then LP0′ can be written as(13)LP0′=1μA+γμS+pS02−S0S−SS0+μE+εE02−E0E−EE0+pS0SS0−1−E0SES0+EE0+βS0fA+σβE0fA−μA+γA=1μA+γμSS02−S0S−SS0+μE+εE03−S0S−EE0−SE0S0E+βSfA1−S0S+βS0fA+σβE0fA−μA+γA≤1μA+γμSS02−S0S−SS0+μE+εE03−S0S−EE0−SE0S0E+1μA+γβS0f′0A+σβE0f′0A−μA+γA=1μA+γμSS02−S0S−SS0+μE+εE03−S0S−EE0−SE0S0E+R0−1A.
Notice the fact(14)2−S0S−SS0≤0,3−S0S−EE0−SE0S0E≤0,for all S,E>0 and the equalities hold if and only if S=S0 and E=E0. Hence, we have LP0′≤0 when R0<1 and LP0′=0 only at P0, which leads to the global stability of P0. This completes the proof.Theorem 2.
IfR0>1, the alcohol-present equilibrium P∗ is globally asymptotically stable.Proof.
Construct a Lyapunov function as follows:(15)LP∗=S−S∗−S∗lnSS∗+E−E∗−E∗lnEE∗+A−A∗−A∗lnAA∗.Then, the time derivative of LP∗ along the solution of (2) satisfies(16)LP∗′=1−S∗S1−qΛ−βSfA−μS+pS+1−E∗EpS−σβEfA−μE+εE+1−A∗AS+σEβfA−μA+γA.
Noticing that(17)1−qΛ=βS∗fA∗+μS+pS∗,pS∗=μE+εE∗+σβE∗fA∗,we have(18)LP∗′=μS+pS∗2−S∗S−SS∗+1−S∗SβS∗fA∗−βSfA+μE+εE∗2−E∗E−EE∗+pS∗SS∗−1−E∗SES∗+E∗E+1−E∗EσβE∗fA∗−σβEfA+1−A∗AβSfA+σβEfA−βS∗fA∗+σβE∗fA∗AA∗=μSS∗2−S∗S−SS∗+μE+εE∗3−S∗S−EE∗−SE∗S∗E+βS∗fA∗2−S∗S−SfAA∗S∗fA∗A∗+fAfA∗−AA∗+σβE∗fA∗3−S∗S−SE∗S∗E−EfAA∗E∗fA∗A+fAfA∗−AA∗=μSS∗2−S∗S−SS∗+μE+εE∗3−S∗S−EE∗−SE∗S∗E+βS∗fA∗−gS∗S−gSfAA∗S∗fA∗A+fAfA∗−AA∗+σβE∗fA∗−gS∗S−gSE∗S∗E−gEfAA∗E∗fA∗A+fAfA∗−AA∗=μSS∗2−S∗S−SS∗+μE+εE∗3−S∗S−EE∗−SE∗S∗E+βS∗fA∗−gS∗S−gSfAA∗S∗fA∗A−gAfA∗A∗fA+fAfA∗−AA∗1−fA∗fA+σβE∗fA∗−gS∗S−gSE∗S∗E−gEfAA∗E∗fA∗A−gAfA∗A∗fA+fAfA∗−AA∗1−fA∗fA,where gx=x−1−lnx≥0 for all x>0. Furthermore, the assumption of fA implies that(19)fAfA∗−AA∗1−fA∗fA≤0,which leads to LP∗′≤0 and the largest compact invariant set of LP∗′≤0 is the singleton P∗. Therefore, by LaSalle’s invariance principle, P∗ is globally asymptotically stable when R0>1.


## 3. Discrete Drinking Model

### 3.1. Model Formulation

By applying the Mickens nonstandard finite difference scheme, we discretize the continuous model (1) and construct the following discrete system:(20)Sn+1−Snh=1−qΛ−βSn+1fAn−μS+pSn+1,En+1−Enh=pSn+1−σβEn+1fAn−μE+εEn+1,An+1−Anh=Sn+1+σEn+1βfAn−μA+γAn+1,Rn+1−Rnh=qΛ+γAn+1+εEn+1−μRRn+1,where n∈N, h>0, denoting the time step size. Sn, En, An, and Rn are the approximations of variables St, Et, At, and Rt in (1) at the discrete time t=nh, respectively, and the parameters in (20) are the same as in (1).Similar to system (1), the solution of (20) with nonnegative initial values is nonnegative and exists for all t≥0, and it is suffice to consider the following reduced system in Ω:(21)Sn+1−Snh=1−qΛ−βSn+1fAn−μS+pSn+1,En+1−Enh=pSn+1−σβEn+1fAn−μE+εEn+1,An+1−Anh=Sn+1+σEn+1βfAn−μA+γAn+1.Furthermore, it is easy to show that system (20) may have exactly the same equilibria as P0 and P∗ in the corresponding continuous system (1). The positive alcohol-present equilibrium P∗ exists for all R0>1 with the alcoholics determined by the following equation:(22)1βfA∗+μS+p+σpβfA∗+μS+pσβfA∗+μE+ε=μA+γA∗1−qΛβfA∗,where R0 is the basic reproductive number defined in (6).


### 3.2. Stability Analysis

Since the discretization of the continuous model generally performs as a method of approximation for the analytic solution, it is supposed for the discrete model to preserve the dynamical properties of the corresponding continuous model as much as possible. In the following, we will conduct the global stability of the equilibria of system (20) to examine its consistency and efficiency.Theorem 3.
For any time steph>0, the alcohol-free equilibrium P0 of system (20) is globally asymptotically stable when R0<1.Proof.
To use the Jury criterion, we rearrange the discrete system (20) and get the Jacobian matrix at P0 as(23)11+hμS+p0−S0hβf′01+hμS+php1+hμE+ε11+hμE+ε−E0hσβf′01+hμE+ε001+hS0+σE0βf′01+hμA+γ,noting that(24)S0+h1−qΛ=S01+hμS+p,E0+hpS0=E01+hμE+ε.
Therefore, the characteristic equation has eigenvalues:(25)λ1=11+hμS+p,λ2=11+hμE+ε,λ3=1+hS0+σE0βf′01+hμA+γ,with λi<1i=1,2,3 for any h>0 due to the fact S0+σE0βf′0<μA+γ, which implies that P0 is locally stable.
Then, define the following Lyapunov function:(26)Ln=cSn−S0−S0lnSnS0+En−E0−E0lnEnE0+1+hμA+γAn,where c=1/μA+γ. Noticing that the inequality lnx≤x−1 is valid for any x>0, then we have(27)Ln+1−Ln=cSn+1−Sn+S0lnSnSn+1+1+hμA+γAn+1−An+En+1−En+E0lnEnEn+1≤c1−S0Sn+1Sn+1−Sn+1+hμA+γAn+1−An+1−E0En+1En+1−En=hc1−S0Sn+11−qΛ−βSn+1fAn−μS+pSn+1+1−E0En+1pSn+1−σβEn+1fAn−μE+εEn+1+hcSn+1+σEn+1βfAn−μA+γAn.
In view of the equalities1−qΛ=μS+pS0 and pS0=μE+εE0, we obtain(28)Ln+1−Ln≤hcμS+pS02−S0Sn+1−Sn+1S0+μE+εE02−E0En+1−En+1E0+βS0fAn+σβE0fAn−μA+γAn+pS0Sn+1S0−1−E0Sn+1En+1S0+E0En+1=hcμSS02−S0Sn+1−Sn+1S0+μE+εE03−S0Sn+1−En+1E0−Sn+1E0S0En+1+βS0fAn+σβE0fAn−μA+γAn+βSn+1fAn1−S0Sn+1≤hcμSS02−S0Sn+1−Sn+1S0+μE+εE03−S0Sn+1−En+1E0−Sn+1E0S0En+1+hcβS0f′0An+σβE0f′0An−μA+γAn=hcμSS02−S0Sn+1−Sn+1S0+μE+εE03−S0Sn+1−En+1E0−Sn+1E0S0En+1+R0−1An.
Therefore,R0<1 ensures that Ln+1−Ln≤0, which implies that Ln is a monotonic decreasing sequence. Thus, the limit of this nonnegative sequence exists and limn⟶∞Ln+1−Ln=0, which is equivalent to(29)limn⟶∞Sn=S0,limn⟶∞En=E0,limn⟶∞An=A0,when R0<1. It is observed from (20) that limn⟶∞An=0, then it yields that P0 is globally asymptotically stable when R0<1.Theorem 4.
For any time steph>0, if R0>1, the unique alcohol-present equilibrium P∗ is globally asymptotically stable.Proof.
Define a Lyapunov function as(30)Un=Sn−S∗−S∗lnSnS∗+En−E∗−E∗lnEnE∗+An−A∗−A∗lnAnA∗,and then the following holds:(31)Un+1−Un=Sn+1−Sn+S∗lnSnSn+1+En+1−En+E∗lnEnEn+1+An+1−An+A∗lnAnAn+1≤1−S∗Sn+1Sn+1−Sn+1−E∗En+1En+1−En+1−A∗An+1An+1−An=h1−S∗Sn+11−qΛ−βSn+1fAn−μS+pSn+1+h1−E∗En+1pSn+1−σβEn+1fAn−μE+εEn+1+h1−A∗An+1Sn+1+σEn+1βfAn−μA+γAn+1.
For convenience, letgx=x−1−lnx, then gx is nonnegative for all x>0, and it has the minimum value gx=0 only at x=1. Noting that(32)1−qΛ=βS∗fA∗+μS+pS∗,pS∗=μE+εE∗+σβE∗fA∗,then substituting them into the above relation of Un+1−Un leads to(33)Un+1−Un≤μS+pS∗2−S∗Sn+1−Sn+1S∗+1−S∗Sn+1βS∗fA∗+βSn+1fAn+μE+εE∗2−E∗En+1−En+1E∗+pS∗Sn+1S∗−1−E∗Sn+1En+1S∗+E∗En+1+1−E∗En+1σβE∗fA∗−σβEn+1fAn+1−A∗An+1βSn+1fAn+σβEn+1fAn−βS∗fA∗+σβE∗fA∗An+1A∗=hμSS∗2−S∗Sn+1−Sn+1S∗+μE+εE∗3−S∗Sn+1−En+1E∗−Sn+1E∗S∗En+1+βS∗fA∗−gS∗Sn+1−gSn+1fAnA∗S∗fA∗An+1+fAnfA∗−AnA∗+σβE∗fA∗−gS∗Sn+1−gSn+1E∗S∗En+1−gEn+1fAnA∗E∗fA∗An+fAnfA∗−AnA∗=hμSS∗2−S∗Sn+1−Sn+1S∗+μE+εE∗3−S∗Sn+1−En+1E∗−Sn+1E∗S∗En+1+βS∗fA∗−gS∗Sn+1−gSn+1fAnA∗S∗fA∗An+1−gAn+1fA∗A∗fAn+fAnfA∗−An+1A∗1−fA∗fAn+σβE∗fA∗−gS∗Sn+1−gSn+1E∗S∗En+1−gEn+1fAnA∗E∗fA∗An+1−gAn+1fA∗A∗fAn+fAnfA∗−AnA∗1−fA∗fAn.
From the assumption offA, we obtain that fAn/fA∗−An/A∗≥0 for An≤A∗ or, conversely, fAn/fA∗−An/A∗≤0 for An≥A∗, which gives(34)fAnfA∗−AnA∗1−fA∗fAn≤0,and the equality holds only if An=A∗. For all S,E>0, it is easy to verify that(35)2−S∗Sn+1−Sn+1S∗≤0,3−S∗Sn+1−En+1E∗−Sn+1E∗S∗En+1≤0,and the equalities are valid only for Sn=S∗ and En=E∗. Therefore, from the above deduction, we have Un+1−Un≤0 and Un+1−Un=0 if and only if Sn=S∗, En=E∗, and An=A∗. As a result, P∗ is globally asymptotically stable when R0>1.
From the theoretical results in this section, it can be seen that the basic reproductive number is a crucial threshold which determines the dynamics of the model. From epidemiological viewpoint, alcohol abuse will eventually die out whenR0<1. Otherwise, it will persist at the alcohol-present equilibrium when R0>1.


## 3.1. Model Formulation

By applying the Mickens nonstandard finite difference scheme, we discretize the continuous model (1) and construct the following discrete system:(20)Sn+1−Snh=1−qΛ−βSn+1fAn−μS+pSn+1,En+1−Enh=pSn+1−σβEn+1fAn−μE+εEn+1,An+1−Anh=Sn+1+σEn+1βfAn−μA+γAn+1,Rn+1−Rnh=qΛ+γAn+1+εEn+1−μRRn+1,where n∈N, h>0, denoting the time step size. Sn, En, An, and Rn are the approximations of variables St, Et, At, and Rt in (1) at the discrete time t=nh, respectively, and the parameters in (20) are the same as in (1).Similar to system (1), the solution of (20) with nonnegative initial values is nonnegative and exists for all t≥0, and it is suffice to consider the following reduced system in Ω:(21)Sn+1−Snh=1−qΛ−βSn+1fAn−μS+pSn+1,En+1−Enh=pSn+1−σβEn+1fAn−μE+εEn+1,An+1−Anh=Sn+1+σEn+1βfAn−μA+γAn+1.Furthermore, it is easy to show that system (20) may have exactly the same equilibria as P0 and P∗ in the corresponding continuous system (1). The positive alcohol-present equilibrium P∗ exists for all R0>1 with the alcoholics determined by the following equation:(22)1βfA∗+μS+p+σpβfA∗+μS+pσβfA∗+μE+ε=μA+γA∗1−qΛβfA∗,where R0 is the basic reproductive number defined in (6).


## 3.2. Stability Analysis

Since the discretization of the continuous model generally performs as a method of approximation for the analytic solution, it is supposed for the discrete model to preserve the dynamical properties of the corresponding continuous model as much as possible. In the following, we will conduct the global stability of the equilibria of system (20) to examine its consistency and efficiency.Theorem 3.
For any time steph>0, the alcohol-free equilibrium P0 of system (20) is globally asymptotically stable when R0<1.Proof.
To use the Jury criterion, we rearrange the discrete system (20) and get the Jacobian matrix at P0 as(23)11+hμS+p0−S0hβf′01+hμS+php1+hμE+ε11+hμE+ε−E0hσβf′01+hμE+ε001+hS0+σE0βf′01+hμA+γ,noting that(24)S0+h1−qΛ=S01+hμS+p,E0+hpS0=E01+hμE+ε.
Therefore, the characteristic equation has eigenvalues:(25)λ1=11+hμS+p,λ2=11+hμE+ε,λ3=1+hS0+σE0βf′01+hμA+γ,with λi<1i=1,2,3 for any h>0 due to the fact S0+σE0βf′0<μA+γ, which implies that P0 is locally stable.
Then, define the following Lyapunov function:(26)Ln=cSn−S0−S0lnSnS0+En−E0−E0lnEnE0+1+hμA+γAn,where c=1/μA+γ. Noticing that the inequality lnx≤x−1 is valid for any x>0, then we have(27)Ln+1−Ln=cSn+1−Sn+S0lnSnSn+1+1+hμA+γAn+1−An+En+1−En+E0lnEnEn+1≤c1−S0Sn+1Sn+1−Sn+1+hμA+γAn+1−An+1−E0En+1En+1−En=hc1−S0Sn+11−qΛ−βSn+1fAn−μS+pSn+1+1−E0En+1pSn+1−σβEn+1fAn−μE+εEn+1+hcSn+1+σEn+1βfAn−μA+γAn.
In view of the equalities1−qΛ=μS+pS0 and pS0=μE+εE0, we obtain(28)Ln+1−Ln≤hcμS+pS02−S0Sn+1−Sn+1S0+μE+εE02−E0En+1−En+1E0+βS0fAn+σβE0fAn−μA+γAn+pS0Sn+1S0−1−E0Sn+1En+1S0+E0En+1=hcμSS02−S0Sn+1−Sn+1S0+μE+εE03−S0Sn+1−En+1E0−Sn+1E0S0En+1+βS0fAn+σβE0fAn−μA+γAn+βSn+1fAn1−S0Sn+1≤hcμSS02−S0Sn+1−Sn+1S0+μE+εE03−S0Sn+1−En+1E0−Sn+1E0S0En+1+hcβS0f′0An+σβE0f′0An−μA+γAn=hcμSS02−S0Sn+1−Sn+1S0+μE+εE03−S0Sn+1−En+1E0−Sn+1E0S0En+1+R0−1An.
Therefore,R0<1 ensures that Ln+1−Ln≤0, which implies that Ln is a monotonic decreasing sequence. Thus, the limit of this nonnegative sequence exists and limn⟶∞Ln+1−Ln=0, which is equivalent to(29)limn⟶∞Sn=S0,limn⟶∞En=E0,limn⟶∞An=A0,when R0<1. It is observed from (20) that limn⟶∞An=0, then it yields that P0 is globally asymptotically stable when R0<1.Theorem 4.
For any time steph>0, if R0>1, the unique alcohol-present equilibrium P∗ is globally asymptotically stable.Proof.
Define a Lyapunov function as(30)Un=Sn−S∗−S∗lnSnS∗+En−E∗−E∗lnEnE∗+An−A∗−A∗lnAnA∗,and then the following holds:(31)Un+1−Un=Sn+1−Sn+S∗lnSnSn+1+En+1−En+E∗lnEnEn+1+An+1−An+A∗lnAnAn+1≤1−S∗Sn+1Sn+1−Sn+1−E∗En+1En+1−En+1−A∗An+1An+1−An=h1−S∗Sn+11−qΛ−βSn+1fAn−μS+pSn+1+h1−E∗En+1pSn+1−σβEn+1fAn−μE+εEn+1+h1−A∗An+1Sn+1+σEn+1βfAn−μA+γAn+1.
For convenience, letgx=x−1−lnx, then gx is nonnegative for all x>0, and it has the minimum value gx=0 only at x=1. Noting that(32)1−qΛ=βS∗fA∗+μS+pS∗,pS∗=μE+εE∗+σβE∗fA∗,then substituting them into the above relation of Un+1−Un leads to(33)Un+1−Un≤μS+pS∗2−S∗Sn+1−Sn+1S∗+1−S∗Sn+1βS∗fA∗+βSn+1fAn+μE+εE∗2−E∗En+1−En+1E∗+pS∗Sn+1S∗−1−E∗Sn+1En+1S∗+E∗En+1+1−E∗En+1σβE∗fA∗−σβEn+1fAn+1−A∗An+1βSn+1fAn+σβEn+1fAn−βS∗fA∗+σβE∗fA∗An+1A∗=hμSS∗2−S∗Sn+1−Sn+1S∗+μE+εE∗3−S∗Sn+1−En+1E∗−Sn+1E∗S∗En+1+βS∗fA∗−gS∗Sn+1−gSn+1fAnA∗S∗fA∗An+1+fAnfA∗−AnA∗+σβE∗fA∗−gS∗Sn+1−gSn+1E∗S∗En+1−gEn+1fAnA∗E∗fA∗An+fAnfA∗−AnA∗=hμSS∗2−S∗Sn+1−Sn+1S∗+μE+εE∗3−S∗Sn+1−En+1E∗−Sn+1E∗S∗En+1+βS∗fA∗−gS∗Sn+1−gSn+1fAnA∗S∗fA∗An+1−gAn+1fA∗A∗fAn+fAnfA∗−An+1A∗1−fA∗fAn+σβE∗fA∗−gS∗Sn+1−gSn+1E∗S∗En+1−gEn+1fAnA∗E∗fA∗An+1−gAn+1fA∗A∗fAn+fAnfA∗−AnA∗1−fA∗fAn.
From the assumption offA, we obtain that fAn/fA∗−An/A∗≥0 for An≤A∗ or, conversely, fAn/fA∗−An/A∗≤0 for An≥A∗, which gives(34)fAnfA∗−AnA∗1−fA∗fAn≤0,and the equality holds only if An=A∗. For all S,E>0, it is easy to verify that(35)2−S∗Sn+1−Sn+1S∗≤0,3−S∗Sn+1−En+1E∗−Sn+1E∗S∗En+1≤0,and the equalities are valid only for Sn=S∗ and En=E∗. Therefore, from the above deduction, we have Un+1−Un≤0 and Un+1−Un=0 if and only if Sn=S∗, En=E∗, and An=A∗. As a result, P∗ is globally asymptotically stable when R0>1.
From the theoretical results in this section, it can be seen that the basic reproductive number is a crucial threshold which determines the dynamics of the model. From epidemiological viewpoint, alcohol abuse will eventually die out whenR0<1. Otherwise, it will persist at the alcohol-present equilibrium when R0>1.


## 4. Numerical Simulations

In this section, numerical simulations of model (1) and model (20) are provided to illustrate the theoretical results, including the global stabilities of the equilibria when the basic reproductive number R0<1 and R0>1, respectively, with the aim to demonstrate the consistency and efficiency of the discrete model in preserving the dynamic behaviors of the corresponding continuous model. In addition, the influence of some parameters on the basic reproductive number is examined, such as p, σ, and ε, which are related with the public health education. In models (1) and (20), the incidence rate is taken as the general form. In the following, we conduct numerical simulations by choosing two specific types of incidence, i.e., fA=A and fA=A/1+αA.


### 4.1. Stability Preserving

(i)
Case 1:fA=A.
In this case, the parameters are chosen asΛ=3.469, q=0.75, μS=0.25, p=0.9, μE=0.25, μA=0.3, ε=0.2, μR=0.25, γ=0.2, β=0.3, and σ=0.2. With these parameters, the basic reproductive number is computed as R0=0.6335<1. Figure 1(a) shows that the alcohol-free equilibrium P0=0.7542,1.5082,0,11.6136 is globally asymptotically stable (variable Rt is omitted for clarity in plots). When p and ε are reduced to p=0.265 and ε=0.1, and the other parameters remain the same, we have R0=1.1634>1, which leads to the global stability of P∗=1.4560,1.0536,0.2689,11.0437, as plotted in Figure 1(b).
For comparison, we conduct simulations for the discrete model (20) with fA=A. The parameters are taken the same as in Figure 1. Firstly, we have R0=0.6335<1, corresponding to the case in Figure 1(a). Let the time step be h=0.4,2,5,10, respectively. As shown in Figure 2, the alcohol-free equilibrium P0 is globally asymptotically stable, which is consistent with Theorem 3. To illustrate the dynamical behaviors of model (20) when R0 being greater than 1, we compare it with Figure 1(b), i.e., R0=1.1634>1. Figure 3 shows that the alcohol-present equilibrium P∗ is globally stable, which is consistent with Theorem 4. From Figures 2 and 3, we can see that the discrete model (20) preserves exactly the dynamical properties of the corresponding continuous model (1), and the increasing time step size from h=0.4 to h=10 has no effect on the stability of the equilibria, indicating the efficiency and consistency of the proposed method for discretization.
(ii)
Case 2:fA=A/1+αA.Figure 1
Global stabilities of the equilibriaP0 and P∗ for model (1) with fA=A. (a) P0 is globally asymptotically stable when R0=0.6335<1. (b) P∗ is globally asymptotically stable when R0=1.1634>1.
(a)
(b)Figure 2
Global stability ofP0 for model (20) with fA=A when R0=0.6335<1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)Figure 3
Global stability ofP∗ for model (20) with fA=A when R0=1.1634>1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)When the general incidence function is specified as this type of saturated form, numerical solutions of model (1) and model (20) are simulated when the basic reproductive number is less or greater than 1, respectively. If the parameters are chosen as Λ=3.8, q=0.75, μS=0.25, p=0.85, μE=0.25, μA=0.28, ε=0.15, μR=0.25, γ=0.26, β=0.3, σ=0.22, and α=0.5, then for this set of parameter values, R0=0.7041<1 and the alcohol-free equilibrium P0=0.8637,1.8352,0,12.5012 is globally asymptotically stable, as shown in Figure 4(a). If p is reduced from 0.85 to 0.24 and the other parameters remain the same,then R0=1.2193>1, and all the solutions of model (1) converge to the alcohol-present equilibrium P∗=1.7498,1.0201,0.1935,12.2133, as plotted in Figure 4(b).Figure 4
Global stabilities of the equilibriaP0 and P∗ for model (1) with fA=A/1+αA. (a) P0 is globally asymptotically stable when R0=0.7041<1. (b) P∗ is globally asymptotically stable when R0=1.2193>1.
(a)
(b)For the corresponding discrete model (20), the numerical simulations are also conducted to check its stability preserving in this case. All the parameter values are taken the same as in Figure 4 and set the time step size h as 0.4,2,5,and10, respectively. With these parameters, we have R0=0.7041<1, corresponding to Figure 4(a). As desired, the alcohol-free equilibrium P0 is asymptotically stable, as shown in Figure 5. When R0=1.2193>1, the simulation result is presented in Figure 6, which is corresponding to Figure 4(b). From the plot, we can see that the numerical solutions of model (20) converge to P∗, showing the expected consistency of the obtained discrete model (20) with the original continuous model (1). Furthermore, it is depicted in Figures 5 and 6 that the dynamical properties will not be influenced by increasing the time step size; i.e., whether alcohol abuse persists or not is completely determined by the basic reproductive number, no matter what the time step is taken.Figure 5
Global stability ofP0 for model (20) with fA=A/1+αA when R0=0.7041<1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)Figure 6
Global stability ofP∗ for model (20) with fA=A/1+αA when R0=1.2193<1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)Remark.
When the general incidence in models (1) and (20) is taken as other forms, such as Beddington–DeAngelis response (βSA/1+aS+bA) and Crowley–Martin response (βSA/1+aS+bA+abSA), numerical simulations about stability preserving present the similar results as in Figures 1–6.


### 4.2. Effect of Public Health Education

In the modelling, we introduce two different susceptibilities due to consideration of public health education, which is perceived as an effective strategy to reduce alcohol abuse. When people are influenced by these educational campaigns, some of them may consequently have a relatively low probability of being affected by alcoholics as well as a high rate of moving to the nondrinker class. In the following, the parameters that reflect the effort of public health education are highlighted, i.e.,p (rate moving from the susceptible to the educated class), σ (reduced infectiousness of educated individuals relative to susceptible ones), and ε (rate moving from the educated class to the no-drinker class).To evaluate the educational campaign and demonstrate the impact of these parameters, we investigate the sensitivity of parameter variations on the basic reproductive numberR0, which plays an important role in determining the dynamics of the model. By (6), we can get the derivatives of R0 about p, σ, and ε, respectively, as follows:(36)∂R0∂p=βf′01−qΛσμS−μE−εμA+γμE+εμS+p2,∂R0∂ε=−βf′01−qΛσpμA+γμE+ε2μS+p<0,∂R0∂σ=βf′01−qΛpμA+γμE+εμS+p>0.For the derivative aboutp, it is negative if only σμS−μE−ε<0. In fact, noting that μS and μE are the death rates in the susceptible and the educated class, the difference between them is typically small enough in most actual situations, implying that ∂R0/∂p<0 can generally hold for all positive parameter values. Therefore, R0 is negatively related to p and ε but positively related to σ, revealing that the more effectively the public health education performs, the better results the alcoholic consumption control will produce.The plots in Figure7 illustrate the dependence of R0 on the parameter pairs of p,ε, p,σ, and σ,ε, using the function of fA=A. The intersecting curves of R0=1 by the two surfaces R0 and the constant 1 reveal the parameter ranges in which the basic reproductive number is less than 1 and then the alcohol abuse can be controlled eventually. In Figure 8, the percent of S, E, and A in the whole population is presented when the three focused parameters p, ε, and σ vary from 0 to 1. As depicted in Figure 8, the percent of alcoholics can be held to a very low level (probably even to zero) as the effect of educational campaign rises, reflected by the increasing p and ε. However, it is also shown by the percent-σ plot that the lowered education strength allows for an obvious increase in the percent of alcoholics, which is induced by the waning knowledge of alcoholism in the public.Figure 7
Dependence ofR0 on the parameter pairs of p,ε, p,σ, and σ,ε.
(a)
(b)
(c)Figure 8
Percent ofS, E, and A in the whole population with varying (a) p, (b) ε, and (c) σ, respectively.
(a)
(b)
(c)


## 4.1. Stability Preserving

(i)
Case 1:fA=A.
In this case, the parameters are chosen asΛ=3.469, q=0.75, μS=0.25, p=0.9, μE=0.25, μA=0.3, ε=0.2, μR=0.25, γ=0.2, β=0.3, and σ=0.2. With these parameters, the basic reproductive number is computed as R0=0.6335<1. Figure 1(a) shows that the alcohol-free equilibrium P0=0.7542,1.5082,0,11.6136 is globally asymptotically stable (variable Rt is omitted for clarity in plots). When p and ε are reduced to p=0.265 and ε=0.1, and the other parameters remain the same, we have R0=1.1634>1, which leads to the global stability of P∗=1.4560,1.0536,0.2689,11.0437, as plotted in Figure 1(b).
For comparison, we conduct simulations for the discrete model (20) with fA=A. The parameters are taken the same as in Figure 1. Firstly, we have R0=0.6335<1, corresponding to the case in Figure 1(a). Let the time step be h=0.4,2,5,10, respectively. As shown in Figure 2, the alcohol-free equilibrium P0 is globally asymptotically stable, which is consistent with Theorem 3. To illustrate the dynamical behaviors of model (20) when R0 being greater than 1, we compare it with Figure 1(b), i.e., R0=1.1634>1. Figure 3 shows that the alcohol-present equilibrium P∗ is globally stable, which is consistent with Theorem 4. From Figures 2 and 3, we can see that the discrete model (20) preserves exactly the dynamical properties of the corresponding continuous model (1), and the increasing time step size from h=0.4 to h=10 has no effect on the stability of the equilibria, indicating the efficiency and consistency of the proposed method for discretization.
(ii)
Case 2:fA=A/1+αA.Figure 1
Global stabilities of the equilibriaP0 and P∗ for model (1) with fA=A. (a) P0 is globally asymptotically stable when R0=0.6335<1. (b) P∗ is globally asymptotically stable when R0=1.1634>1.
(a)
(b)Figure 2
Global stability ofP0 for model (20) with fA=A when R0=0.6335<1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)Figure 3
Global stability ofP∗ for model (20) with fA=A when R0=1.1634>1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)When the general incidence function is specified as this type of saturated form, numerical solutions of model (1) and model (20) are simulated when the basic reproductive number is less or greater than 1, respectively. If the parameters are chosen as Λ=3.8, q=0.75, μS=0.25, p=0.85, μE=0.25, μA=0.28, ε=0.15, μR=0.25, γ=0.26, β=0.3, σ=0.22, and α=0.5, then for this set of parameter values, R0=0.7041<1 and the alcohol-free equilibrium P0=0.8637,1.8352,0,12.5012 is globally asymptotically stable, as shown in Figure 4(a). If p is reduced from 0.85 to 0.24 and the other parameters remain the same,then R0=1.2193>1, and all the solutions of model (1) converge to the alcohol-present equilibrium P∗=1.7498,1.0201,0.1935,12.2133, as plotted in Figure 4(b).Figure 4
Global stabilities of the equilibriaP0 and P∗ for model (1) with fA=A/1+αA. (a) P0 is globally asymptotically stable when R0=0.7041<1. (b) P∗ is globally asymptotically stable when R0=1.2193>1.
(a)
(b)For the corresponding discrete model (20), the numerical simulations are also conducted to check its stability preserving in this case. All the parameter values are taken the same as in Figure 4 and set the time step size h as 0.4,2,5,and10, respectively. With these parameters, we have R0=0.7041<1, corresponding to Figure 4(a). As desired, the alcohol-free equilibrium P0 is asymptotically stable, as shown in Figure 5. When R0=1.2193>1, the simulation result is presented in Figure 6, which is corresponding to Figure 4(b). From the plot, we can see that the numerical solutions of model (20) converge to P∗, showing the expected consistency of the obtained discrete model (20) with the original continuous model (1). Furthermore, it is depicted in Figures 5 and 6 that the dynamical properties will not be influenced by increasing the time step size; i.e., whether alcohol abuse persists or not is completely determined by the basic reproductive number, no matter what the time step is taken.Figure 5
Global stability ofP0 for model (20) with fA=A/1+αA when R0=0.7041<1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)Figure 6
Global stability ofP∗ for model (20) with fA=A/1+αA when R0=1.2193<1: (a) h = 0.4; (b) h = 2; (c) h = 5; (d) h = 10.
(a)
(b)
(c)
(d)Remark.
When the general incidence in models (1) and (20) is taken as other forms, such as Beddington–DeAngelis response (βSA/1+aS+bA) and Crowley–Martin response (βSA/1+aS+bA+abSA), numerical simulations about stability preserving present the similar results as in Figures 1–6.


## 4.2. Effect of Public Health Education

In the modelling, we introduce two different susceptibilities due to consideration of public health education, which is perceived as an effective strategy to reduce alcohol abuse. When people are influenced by these educational campaigns, some of them may consequently have a relatively low probability of being affected by alcoholics as well as a high rate of moving to the nondrinker class. In the following, the parameters that reflect the effort of public health education are highlighted, i.e.,p (rate moving from the susceptible to the educated class), σ (reduced infectiousness of educated individuals relative to susceptible ones), and ε (rate moving from the educated class to the no-drinker class).To evaluate the educational campaign and demonstrate the impact of these parameters, we investigate the sensitivity of parameter variations on the basic reproductive numberR0, which plays an important role in determining the dynamics of the model. By (6), we can get the derivatives of R0 about p, σ, and ε, respectively, as follows:(36)∂R0∂p=βf′01−qΛσμS−μE−εμA+γμE+εμS+p2,∂R0∂ε=−βf′01−qΛσpμA+γμE+ε2μS+p<0,∂R0∂σ=βf′01−qΛpμA+γμE+εμS+p>0.For the derivative aboutp, it is negative if only σμS−μE−ε<0. In fact, noting that μS and μE are the death rates in the susceptible and the educated class, the difference between them is typically small enough in most actual situations, implying that ∂R0/∂p<0 can generally hold for all positive parameter values. Therefore, R0 is negatively related to p and ε but positively related to σ, revealing that the more effectively the public health education performs, the better results the alcoholic consumption control will produce.The plots in Figure7 illustrate the dependence of R0 on the parameter pairs of p,ε, p,σ, and σ,ε, using the function of fA=A. The intersecting curves of R0=1 by the two surfaces R0 and the constant 1 reveal the parameter ranges in which the basic reproductive number is less than 1 and then the alcohol abuse can be controlled eventually. In Figure 8, the percent of S, E, and A in the whole population is presented when the three focused parameters p, ε, and σ vary from 0 to 1. As depicted in Figure 8, the percent of alcoholics can be held to a very low level (probably even to zero) as the effect of educational campaign rises, reflected by the increasing p and ε. However, it is also shown by the percent-σ plot that the lowered education strength allows for an obvious increase in the percent of alcoholics, which is induced by the waning knowledge of alcoholism in the public.Figure 7
Dependence ofR0 on the parameter pairs of p,ε, p,σ, and σ,ε.
(a)
(b)
(c)Figure 8
Percent ofS, E, and A in the whole population with varying (a) p, (b) ε, and (c) σ, respectively.
(a)
(b)
(c)


## 5. Conclusions

In this paper, we investigated the threshold dynamics of an SEAR alcoholism model with the form of continuous equations and discrete equations, noting that the discrete model is formulated by applying the Michens nonstandard finite difference scheme to the corresponding continuous model. With this method of discretization, as demonstrated by the theoretical and numerical results obtained in this paper, the dynamical properties of the continuous system can be preserved with much high efficiency when compared with the traditional schemes such as forward Euler and Runge–Kutta, which sometimes fail generating oscillations, bifurcations, and chaos, and even result in false steady states [35, 36]. Specifically, not only the expressions of equilibria and the basic reproductive number but also the threshold dynamics are identical in both models (1) and (20), which confirm its consistency of the discrete model with the original continuous model in the qualitative properties.For some troubling contagions, such as alcohol drinking, improving public knowledge forms part of the perceived strategy to prevent the spread. Public health education, which may substantially contribute to an effective management of this health problem, is incorporated in our model by introducing two different susceptibilities. To analyze the impact of parameters related with public health education, numerical simulations are conducted to demonstrate the effect of them on the basic reproductive number and on the percentage of a different drinking status in whole population. The results show that educational campaign plays an important role in controlling alcoholism (seen in Figures7 and 8), which can be held to a very low level (eventually dies out) when some parameters related with public health education remain in certain intervals. Meanwhile, educational campaign is also expected to perform efficiently in reducing the initiation and maintenance of alcohol consumption.



---

*Source: 4272914-2020-01-13.xml*