Global dynamics of a delayed alcoholism model with the effect of health education

1.   Introduction

Alcoholism can be viewed as a social epidemic ^[1], which can bring great harm to individuals, families and society. The study of alcoholism has become an important aspect of social epidemic. Alcohol abuse can also lead to a range of negative social effects such as violence, antisocial and criminal behavior. Alcohol consumption has been identified as a major contributor to the global burden of chronic disease, injury and economic cost ^[2,3,4]. The World Health Organization reports the harmful use of alcohol causing approximately 3.3 million deaths every year (or 5.9% of all the global deaths), and 5.1% of the global burden of disease is attributable to alcohol consumption ^[5]. Therefore, alcoholism is the main targeted healthy risky behavior due to the high relevance of negative health and social effects. The spread of healthy risky behavior within a community can be viewed as a diffusion process with its own incidence rate. In this situation, the social interaction is considered to be the key factor in spreading the alcohol behavior which can result in adverse health effects. The prevention and control of alcoholism is an urgent problem.

Mathematical modelling is a powerful tool for solving problems in various fields ^[6,7,9,8,10]. In the last few decades, mathematical models for human behaviors related to addictions have been developed from epidemiological models for the spread of infectious diseases, including drinking, smoking and drug use, etc., see ^{[11,12,13,14]} and the references contained therein. In particular, there are several different mathematical models for alcohol problems that have been formulated and studied recently ^{[15,16,17,18,19,20,21,22,23,24,25,26,27]}. Benedict ^[15] modelled alcoholism as a contagious disease and studied how "infected" drinking buddies spread problem drinking. Manthey et al. ^[16] studied campus drinking and suggested that the reproductive numbers are not sufficient to predict whether drinking behavior will persist on campus and that the pattern of recruiting new members play a significant role in the reduction of campus alcohol problems. Santonja et al. ^[17] proposed a mathematical model for alcohol consumption in Spanish population. Predictions about the future behavior of the alcohol consumption in Spain are presented using this model. Environmental and peer influence combinations to create a culture of drinking were studied in ^{[18,19,20,21]}. In addition, the two-stage models: one stage where people who admit having the alcohol problem and other stage where people who do not admit to having the alcohol problem have been developed in ^[22,23]. Bhunu ^[24] studied the co-interaction of alcoholism and smoking in a community. Walters et al. ^[25] also discussed alcohol problems, and their results showed that an increase in the recovery rate decreased the proportion of binge drinkers in the population. Huo et al. ^[26] considered the effect of constant immigration on drinking behavior. Wang et al. ^[27] presented an alcoholism model with two control strategies to gain insights into this increasingly concern about health and social phenomenon. The optimal control strategies are derived by proposing an objective functional and using Pontryagins Maximum Principle.

It is well known that the prevalence of any epidemic is strongly dependent on the social behavior of individuals. Human behavior and social response play a very important role in the transmission of social epidemic. So, the rational way is to make people aware about the alcohol problems through the media, which can not only influence the individuals' behavior but also increase the governmental health care involvement to control the spread of heavy drinking. In recent years, many mathematical models have been used for studying the impact of public health education by media on epidemic outbreaks ^{[28,29,30,31,32,33,34,35,36,37,38,39]} and the references cited therein. These studies suggested that education and media have a huge impact in controlling the spread of infectious diseases. Recently, for drinking problem, Huo et al. ^[40] have studied drinking dynamics and focused on awareness programs and treatment in the modelling process. They have extended the model in ^[39] via including a treatment class and established some sufficient conditions for the stability of the alcohol-free and alcohol-present equilibria. Xiang et al. ^[41] also studied a drinking model with public health educational campaigns. Using Lyapunov function, the global stability of equilibria of the model is derived. Their results showed that the public health educational campaigns is one of the commonly used effective measures in order to prevent and reduce the alcohol problems. On the other hand, recently, Ma et al. ^[42,43] proposed that the multi-group alcoholism model can contribute to the control of alcohol problems in more realistic situations.

At the same time, time delays have been considered into the infections diseases models by many authors ^{[44,45,46,47,48]}. Delay differential equations exhibit much more complicated dynamics than ordinary differential equations since the time delay may affect the stability of the system, even lead to instability, oscillation or bifurcation phenomena. Because time delay is common in development of alcohol consumption habit in population, we should note that there are delays in the "infection" process which the alcohol consumption habit develops in a individual who is "infected" by alcoholics. So, it is more realistic to consider the time delay in the modelling alcoholism process. To our knowledge, the results about alcoholism epidemic model with delays are comparatively scarce. Recently, Huo et al. ^[49] introduced a more realistic binge drinking model with delay. They concluded that, regardless of the time delay length, the alcohol-free equilibrium is globally asymptotically stable. Numerical simulations showed that the alcohol-present equilibrium is globally asymptotically stable. Ma et al. ^[50] formulated a dynamic alcohol consumption model with awareness programs and one delay. The results showed that the time delay in alcohol consumption habit which develops in susceptible population may result in a Hopf bifurcation by increasing the value of time delay.

Motivated by the above works and based on the previous work in ^[41], we formulate an alcoholism model that incorporates both public health education and three delays to study the dynamics and control of drinking. We incorporate the delays $\tau_{1}$, $\tau_{2}$ to describe the time needed for a susceptible individual to become an alcohol user, and the delay $\tau_{3}$ to describe temporarily recovered population that will take a period of time to become a alcohol consumer again. Here, time delay for the "infection" and "recovery" are included together. The goal of this paper is to analyze the effects of public health education and three delays on alcohol control. So, the research work of this paper improves the existing results.

The remaining part of this paper is organized as follows. In next section, we formulate the delayed mathematical model. The positivity and ultimate boundedness of the solutions for system are shown in Section 3. In Section 4, for four different delay cases, the global dynamics of the delayed alcoholism model is discussed in detail. In Section 5, the impacts of public health education and delays on alcohol control will be discussed. The sensitivity analysis of system parameters is given in Section 6. Then, in order to support our theoretical predictions, some numerical simulations are included in Section 7. Finally, a brief conclusion is also given in the last section.

2.   Model formulation

In 2015, Xiang et al. ^[41] investigated a SEARQ drinking model with the public health educational campaigns, as the following drinking dynamics:

where the total population in the model is divided into five types: at time $t$, $S(t)$ denotes the numbers of susceptible individuals who consume alcohol in moderation and do not accept the public health education, but may develop problems with alcohol; $E(t)$ is referred to educated susceptible individuals who consume alcohol in moderation and accept the public health education, in which the average rate of turning to drink is reduced relatively to the tuen rate of uneducated individuals; $A(t)$ is alcoholics class who have drinking problems or addictions; $R(t)$ is temporarily recovered drinkers, that is, former alcoholics who have entered treatment and are abstaining from alcohol; $Q(t)$ represent quit drinkers who permanently quit drink. All parameters are assumed to be positive constants. It is also assumed that the drinking only spreads due to the direct contact between susceptible individuals or educated drinkers and alcoholics. The global stability of equilibria are obtained by constructing suitable Lyapunov functions (see Theorem 4.1 in ^[41]).

As is known to us, there are delays in the "infection" or "recovery" process of drinking. So we propose a mathematical model with health education and three delays to study the dynamics and control of alcoholism. In this paper, we incorporate three delays into the alcoholism model (2.1). Assume that homogenous mixing and each individual has the same chance to get infected. The variables and model structure are described in Figure 1.

Our model can be described by the following ordinary differential equations:

Here the meanings of parameters are the same as system (2.1). Time delay $\tau_{1}\geq0$ to describe the time needed for a uneducated susceptible individual to become an heavy alcohol consumer, $\tau_{2}\geq0$ to describe the time needed for a educated susceptible individual to become an heavy alcohol consumer, and time delay $\tau_{3}\geq0$ to describe the time needed for a temporarily recovered drinkers re-enter the alcoholics groups. There is a time delay during which the alcohol consumption habit develops in his/her body. It is only after such time delay that he/she becomes a alcohol consumer.

Let $\tau = \max\{\tau_{1}, \tau_{2}, \tau_{3}\}$. We denote by $C = C([-\tau, 0], \mathbb{R})$ the Banach space of continuous real-valued functions on the interval $[-\tau, 0]$, with the sup-norm

The nonnegative cone of $C$ is defined as

We suppose that the initial conditions for system (2.2) take the forms:

Remark 2.1. Assume $\tau_{i}\equiv0 (i = 1, 2, 3)$, the system (2.2) becomes the system (2.1) in ^[41]. Therefore, system (2.2) expands system (2.1).

3.   Positivity and ultimate boundedness

Since the variable $Q$ does not appear explicitly in the first four equations in (2.2), then the dynamics of system (2.2) is the same as the following reduced system:

Thus, in the rest of the article, we are concerned with the above system (3.1). For system $\rm(3.1)$, it is important to prove that all the state variables remain nonnegative (for biological reasons). We give the following lemma.

Lemma 3.1. Consider system (3.1) with the initial conditions (2.3), we have

(i) All solutions of system $\rm(3.1)$ are strictly positive;

(ii) All solutions of system $\rm(3.1)$ are ultimately bounded if $\delta < \mu$.

Proof. (i) We prove the system $\rm(3.1)$ with initial conditions (2.3) has strictly positive solutions for all $t\geq 0$.

Assuming the contrary and letting $t_{1} > 0$ be the first time such that $S(t_{1}) = 0$. From Eq $\rm(3.1a)$, we have $S^{'}(t_{1}) = q\mu\Lambda > 0$, hence for sufficiently small $\eta > 0$, $S(t) < 0$ for $t\in(t_{1}-\eta, t_{1})$. But this contradicts $S(t) > 0$ for $t\in[0, t_{1})$. It prove that $S(t) > 0$ for all $t\geq0$. Suppose that $t_{2} > 0$ be the first time such that $E(t_{2}) = 0$. From Eq $\rm(3.1b)$, we have $E^{'}(t_{2}) = (1-q)\mu\Lambda +pS(t_{2}) > 0$. We use the same method to prove $E(t) > 0$ for all $t\geq0$.

From Eq $\rm(3.1c)$, we can obtain

By using the method of variation of constant and the step-to-step integration, we get

which implies that $A(t) > 0$ for all $t > 0$.

Then from Eq $\rm(3.1d)$, we obtain

hence

which shows $R(t) > 0$ for all $t > 0$ because of $R(0) = \phi_{4}(0) > 0$.

(ii) We show that the solutions of system $\rm(3.1)$ are ultimately bounded for all $t\geq0.$

According to Eq $\rm(3.1a)$, we have

So, $S(t)$ is ultimately bounded.

According to Eq $\rm(3.1b)$, we have

and $E(t)$ is ultimately bounded, too. Then, we define

we have $X(t)\geq0$, and yield

From Eq (3.2), we choose $T > 0$ so large such that $t\geq T-\tau_{1}-\tau_{3} > 0$, then

and from the definition of $X(t)$, we obtains

Then, we get from Eq (3.3) that

where $M = \mu\Lambda(1+\frac{q}{\mu+p}).$ Following the method of Kuang ^[45], we assume that $\delta < \mu,$ by considering the auxiliary system

there is a unique positive steady state $z^{*} = \frac{M}{\mu-\delta}$ is globally asymptotically stable. It follows from Eq (3.4) that

Therefore, $S(t), E(t), A(t)$ and $R(t)$ are ultimately bounded. This proof is completed.

As a consequence of Lemma 3.1, we know that the dynamics of system $\rm(3.1)$ can be analyzed in the following bounded feasible region, namely, for sufficiently small $\varepsilon > 0$ such that

Hence, the region $\Omega$ is positive invariant set for system $\rm(3.1)$ and the system is well posed.

Remark 3.1. From the above discussion, we can see the parameter $\delta$ (i.e. the fraction of compartment $R$ who relapse into drinking compartment $A$) could not be too large ($\delta < \mu$). The public health educational campaigns is an effective measure in reducing the value of return proportion $\delta$.

4.   Equilibria and stability

Following the compartment approach, the alcoholism generation number $R_{0}$ of the alcohol consumption model (2.2) can be easily obtained by the next generation matrix method ^[51], which is given by the following expression:

where $S_0 = \frac{q\mu\Lambda}{\mu+p}, E_0 = \frac{\mu\Lambda[p+\mu(1-q)]}{(\mu+\varepsilon)(\mu+p)}$, $m = \mu+a_{1}+\gamma$, $n = \mu+a_{2}+\xi+\delta$, and $mn > \gamma\delta$. It acts as a threshold as is shown in the following results.

For convenience and simplicity, we still set $P_0$ is the alcohol-free equilibrium and $P^{*}$ is the unique alcohol-present equilibrium for system $\rm(3.1)$. In addition, from the system $\rm(3.1)$, we have

where $A^{*}$ is the unique positive root of equation $H(A) = 0$ (see Eq (13) in ^[41]). Then, we have the following results.

Lemma 4.1. Consider system (3.1) with the initial conditions (2.3), we have

(i) System $\rm(3.1)$ always has the alcohol-free equilibrium $P_0(S_0, E_0, 0, 0)$;

(ii) System $\rm(3.1)$ also has a unique alcohol-present equilibrium $P^{*}(S^{*}, E^{*}, A^{*}, R^{*})$ if $R_{0} > 1$.

In the following subsection, we investigate the stability of the equilibria of system $\rm(3.1)$. Considering the impacts of health education on the susceptible individuals and treatment on the temporarily recovered individuals, the time delays of three compartments are different, we will be divided into several cases to discuss.

4.1.   Case I: $\tau_{1}=\tau_{2}=\tilde{\tau}=0, \tau_{3} > 0$

1) Stability analysis of $P_{0}$

We linearize the system $\rm(3.1)$ about $P_0$, and the characteristic equation of system $\rm(3.1)$ at $P_0$ is given by

Obviously, we have

So, our focus is to discuss the solution of the following equation

First, we set $c_{1} = \beta(S_{0}+\sigma E_{0}) > 0$. If $R_{0} < 1$ and $\tilde{\tau} = 0, \tau_{3} > 0$, the Eq (4.2) becomes

where $d_{1} = m+n-c_{1} > 0$, $d_{0} = n(m-c_{1}) > 0$.

Notice that 0 is not a root of Eq (4.3) because of $R_{0} < 1$. Following the method in ^[54], we let $\lambda = i\omega_{1}(\omega_{1} > 0)$ be a purely imaginary root of Eq (4.3). Substituting it into Eq (4.3) and separating the real and imaginary parts, then squaring and adding both equations, we can obtain

Let $u_{1} = \omega_{1}^{2}$, then Eq (4.4) becomes

where $d^{2}_{1}-2d_{0} = (m-c_{1})^{2}+n^{2} > 0$ and $d^{2}_{0}-(\gamma\delta)^{2} > 0$ because of $d_{0}-\gamma\delta > 0$, so $G(u_{1}) = 0$ has no positive root. Hence, we easily obtain all roots of Eq (4.1) have negative real parts. The alcohol-free equilibrium $P_{0}$ is locally asymptotically stable for any $\tau_{3} > 0$ with $\tau_{1} = \tau_{2} = \tilde{\tau} = 0$. In the following, we prove $P_{0}$ is globally attractive in $\Omega$ for any $\tau_{3} > 0$ with $\tau_{1} = \tau_{2} = \tilde{\tau} = 0$.

We define the following Lyapunov function $L^{1}:C\times C\times C\times C\rightarrow \mathbb{R}$:

where the notation $x_{t}(s) = x(t+s)$ for $s\in[-\tau, 0)$, and $x_{t}(0) = x(t)$. The time derivative of $L^{1}$ along the solution of system $\rm(3.1)$ ($\tilde{\tau} = 0, \tau_{3} > 0$) is

Using Lemma 3.1, we know that $S(t)$ and $E(t)$ are bounded, that is, we obtain that $S(t)\leq S_{0}$ and $E(t)\leq E_{0}$. It is easy to see that $\frac{dL^{1}}{dt}\mid_{(3.1)}\leq 0$ if $\beta(S_{0}+\sigma E_{0}) < \mu+a_{1}$. This condition is equivalent to the transmission rate $\beta$ meets the following condition:

(H1): $\beta < \frac{\mu+a_{1}}{S_{0}+\sigma E_{0}}$.

Further, $\frac{dL^{1}}{dt} = 0$ if and only if $A(t) = 0.$ It can be verified that the maximal invariant set in $\Big\{\frac{dL^{1}}{dt}\mid_{(3.1)} = 0\Big\}$ is the singleton $\{P_{0}\}.$ By the LaSalle's invariance principle, we can conclude that $P_{0}$ is globally attractive in $\Omega$.

Second, for $R_{0} > 1$, it is not difficult to verify that Eq (4.3) has a root with a positive real part. In fact, $F(\tilde{\tau}, \tau_{3}, 0) < 0$ if $R_{0} > 1$, and $\underset{\lambda\rightarrow\infty}{\mbox{lim}}F(\tilde{\tau}, \tau_{3}, \lambda) = \infty$. Therefore, the following theorem is obtained.

Theorem 4.1. In Case I, we have

(i) If $R_{0} < 1$ and (H1) holds, then the alcohol-free equilibrium $P_0$ is globally asymptotically stable for any $\tau_{3} > 0$, $\tilde{\tau} = 0$;

(ii) If $R_{0} > 1$, then the alcohol-free equilibrium $P_0$ is unstable for any $\tau_{3} > 0, \tilde{\tau} = 0$.

2) Stability analysis of $P^{*}$

We let $s(t) = S(t)-S^{*}, e(t) = E(t)-E^{*}, a(t) = A(t)-A^{*}, r(t) = R(t)-R^{*}$, the linearization of system $\rm(3.1)$ at $P^{*}$ is given by

Further, the Jacobi matrix $J$ at $P^{*}$ can be written as

The characteristic equation of system $\rm(3.1)$ at $P^{*}$ is $\text{det}(\lambda I-J) = 0$, where $I$ is the identity matrix and $\lambda$ is the eigenvalue of matrix $J$. By computation, the characteristic equation of system $\rm(3.1)$ at $P^{*}$ can be expressed by

where

$n_{1} = 2\mu+p+\varepsilon+m+n(\sigma+1)\beta A^{*}$,

$n_{2} = \sigma\beta^{2}(A^{*})^{2}+[(\mu+n)(1+\sigma)+\sigma p+\varepsilon+m\sigma+1]\beta A^{*}+\sigma\beta(\mu+p)(\mu+\varepsilon)+n(2\mu+p+\varepsilon)+m(2\mu+p+n+\varepsilon)$,

$n_{3} = (m+n)[\sigma\beta^{2}(A^{*})^{2}+(\mu+\varepsilon)\beta A^{*}+\sigma(p+\mu)\beta A^{*}+(\mu+p)(\mu+\varepsilon)]+mn(\sigma+1)\beta A^{*}+n(2\mu+p+\varepsilon)$,

$n_{4} = mn[\sigma\beta^{2}(A^{*})^{2}+\beta A^{*}(\mu+\varepsilon+\sigma p+\sigma\mu)+(\mu+p)(\mu+\varepsilon)]$,

$r_{1} = \beta S^{*}+\sigma\beta E^{*}$,

$r_{2} = \beta^{2}S^{*}A^{*}(\sigma+2)+\sigma\beta^{2}E^{*}A^{*}(2\sigma+1)+\beta(2\mu+p+n+\varepsilon)(S^{*}+\sigma E^{*})$,

$r_{3} = (\beta S^{*}+\sigma\beta E^{*})[\sigma\beta^{2}(A^{*})^{2}+(\sigma(\mu+p)+n\sigma+1)\beta A^{*}]+(\mu+p)(\mu+\varepsilon)+n(2\mu+p+\varepsilon)+\sigma^{2}\beta^{2}E^{*}A^{*}(\beta A^{*}+\mu+p+n)\\ +\beta^{2}S^{*}A^{*}(\sigma\beta A^{*}+\mu+\varepsilon+n)$,

$r_{4} = (\beta S^{*}+\sigma\beta E^{*})n[\sigma\beta^{2}(A^{*})^{2}+\beta A^{*}(\mu+\varepsilon+\sigma p+\sigma\mu)+(\mu+p)(\mu+\varepsilon)]+n\beta^{2}S^{*}A^{*}(\sigma\beta A^{*}+\mu+\varepsilon)+np\sigma\beta^{2}(A^{*})^{2}\\ +n\sigma^{2}\beta^{2}E^{*}A^{*}(\beta A^{*}+\mu+p)$,

$s_{2} = \gamma\delta$,

$s_{3} = \gamma\delta[\beta A^{*}(1+\sigma)+(2\mu+p+\varepsilon)]$,

$s_{4} = \gamma\delta[\sigma\beta^{2}(A^{*})^{2}+\beta A^{*}(\mu+\varepsilon+\sigma\mu+\sigma p)+(\mu+p)(\mu+\varepsilon)]$.

Here $\tilde{\tau} = 0$, $\tau_{3} > 0$, so the characteristic equation Eq (4.7) can be simplified to

where $w_{1} = n_{1}-r_{1}, w_{2} = n_{2}-r_{2}, w_{3} = n_{3}-r_{3}, w_{4} = n_{4}-r_{4}.$

(H2): Suppose that the following conditions holds:

(i) $w_{1} > 0, w_{1}(w_{2}-s_{2}) > w_{3}-s_{3}, w_{4}-s_{4} > 0;$

(ii) $w_{1}[(w_{3}-s_{3})(w_{2}-s_{2})-w_{1}(w_{4}-s_{4})] > (w_{3}-s_{3})^{2}.$

By the Routh-Hurwitz criteria, we will known that all roots of Eq (4.8) with $\tau_{3} = 0$ have negative real parts. We now claim the locally stability of the alcohol-present equilibrium $P^{*}$ in Case I.

Theorem 4.2. In case I, the alcohol-present equilibrium $P^{*}$ is locally asymptotically stable for any $\tau_{3} > 0$ if $R_{0} > 1$ and (H2) holds.

4.2.   Case II: $\tau_{1}=\tau_{2}=\tilde{\tau} > 0, \tau_{3}=0$

1) Stability analysis of $P_{0}$

If $R_{0} < 1$ and $\tilde{\tau} > 0, \tau_{3} = 0$, from the discussion of part (1) in Case I, the characteristic equation $F(\tilde{\tau}, \tau_{3}, \lambda)$ of system $\rm(3.1)$ at $P_0$ becomes

where $b_{1} = m+n > 0$, $b_{0} = mn-\gamma\delta > 0$, $c_{1} = \beta(S_{0}+\sigma E_{0}) > 0$, $c_{0} = \beta(S_{0}+\sigma E_{0})n = nc_{1} > 0$.

Notice that 0 is not a root of Eq (4.9) because of $R_{0} < 1$. Let $\lambda = i\omega_{2}(\omega_{2} > 0)$ be a purely imaginary root of Eq (4.9), we have

Separating the real and imaginary parts, we can get the following equations

Squaring and adding both equations of Eq (4.10), it follows that

Let $u_{2} = \omega_{2}^{2}$, then Eq (4.11) becomes

From $m > c_{1}$, we easily get $b^{2}_{0}-c^{2}_{0} > 0$ and $b^{2}_{1}-2b_{0}-c^{2}_{1} = n^{2}+2\gamma\delta+m^{2}-c^{2}_{1} > 0$. So Eq (4.12) has no positive root. Thus, all roots of Eq (4.9) have negative real parts. The alcohol-free equilibrium $P_{0}$ is locally asymptotically stable for any $\tilde{\tau} > 0$ with $\tau_{3} = 0$.

In the following, we prove $P_{0}$ is globally attractive in $\Omega$ for any $\tilde{\tau} > 0$, $\tau_{3} = 0$. Define Lyapunov function $L^{2}:C\times C\times C\times C\rightarrow \mathbb{R}$:

where the notation $x_{t}(s)$ and $x_{t}(0)$ are the same meaning as in $L^{1}$. The time derivative of $L^{2}$ along the solution of system $\rm(3.1)$ with $\tau_{1} = \tau_{2} = \tilde{\tau} > 0, \tau_{3} = 0$ is

If (H1) holds, it ensures that $\frac{dL^{2}}{dt}\mid_{(3.1)}\leq 0$, and $\frac{dL^{2}}{dt} = 0$ if and only if $A(t) = 0.$ It can be verified that the maximal invariant set in $\Big\{\frac{dL^{2}}{dt}\mid_{(3.1)} = 0\Big\}$ is the singleton $\{P_{0}\}.$ By the LaSalle's invariance principle, we can conclude that $P_{0}$ is globally attractive in $\Omega$. Therefore, we have the following theorem.

Theorem 4.3. In case II, we have

(i) If $R_{0} < 1$ and (H1) holds, then the alcohol-free equilibrium $P_0$ is globally asymptotically stable for any $\tilde{\tau} > 0$, $\tau_{3} = 0$;

(ii) If $R_{0} > 1$, then the alcohol-free equilibrium $P_0$ is unstable for any $\tilde{\tau} > 0$, $\tau_{3} = 0$.

2) Stability analysis of $P^{*}$

Here $\tilde{\tau} > 0$, $\tau_{3} = 0$, similar to the discussion of part (2) in Case I, the characteristic equation Eq (4.7) can be simplified to

where $h_{1} = n_{1}, h_{2} = n_{2}-s_{2}, h_{3} = n_{3}-s_{3}, h_{4} = n_{4}-s_{4}.$

(H3): The coefficients of Eq (4.14) satisfy the Routh-Hurwitz criteria when $\tilde{\tau} = 0$.

Theorem 4.4. In case II, if $R_{0} > 1$ and (H3) holds, then the alcohol-present equilibrium $P^{*}$ is locally asymptotically stable for any $\tilde{\tau} > 0$.

4.3.   Case III: $\tau_{1}=\tau_{2}=\tau_{3}=\tau > 0$

1) Stability analysis of $P_{0}$

For $\tau_{1} = \tau_{2} = \tau_{3} = \tau > 0$, the characteristic equation of system $\rm(3.1)$ at $P_0$ is given by

First, let $\lambda = i\omega(\omega > 0)$ be a purely imaginary root of Eq (4.15). We have

Then, we get

Separating the real and imaginary parts, we can get the following equations

Squaring and adding both equations of Eq (4.16), it follows that

Let $u = \omega^{2}$, then Eq (4.17) becomes

From $R_{0} < 1$, we easily obtain $m^{2}+n^{2}-\beta^{2}(S_{0}+\sigma E_{0})^{2} > 0$ and $(mn)^{2}-\Big[\gamma\delta+\beta(S_{0}+\sigma E_{0})n\Big]^{2} > 0$. So, Eq (4.18) has no positive root. Thus, all roots of Eq (4.15) have negative real parts, and the alcohol-free equilibrium $P_{0}$ is locally asymptotically stable for any $\tau > 0$.

Second, we prove $P_{0}$ is globally attractive in $\Omega$ for any $\tau > 0$. To prove this, we consider the following Lyapunov function $L^{3}:C\times C\times C\times C\rightarrow \mathbb{R}$:

where the notation $x_{t}(s)$ and $x_{t}(0)$ are the same meaning as in $L^{1}$. The time derivative of $L^{3}$ along the solution of system $\rm(3.1)$ with $\tau_{i} = \tau > 0\; (i = 1, 2, 3)$ is

If (H1) holds, then it ensures that $\frac{dL^{3}}{dt}\mid_{(3.1)}\leq 0$, and $\frac{dL^{3}}{dt} = 0$ if and only if $A(t) = 0.$ It can be verified that the maximal invariant set in $\Big\{\frac{dL^{3}}{dt}\mid_{(3.1)} = 0\Big\}$ is the singleton $\{P_{0}\}.$ By the LaSalle's invariance principle, we can conclude that $P_{0}$ is globally attractive in $\Omega$. We claim the following result.

Theorem 4.5. For Case III, we get

(i) If $R_{0} < 1$ and (H1) holds, then the alcohol-free equilibrium $P_0$ is globally asymptotically stable for any $\tau_{i} = \tau > 0\; (i = 1, 2, 3)$;

(ii) If $R_{0} > 1$, then the alcohol-free equilibrium $P_0$ is unstable for any $\tau_{i} = \tau > 0\; (i = 1, 2, 3)$.

2) Stability analysis of $P^{*}$

For $R_{0} > 1$, system $\rm(3.1)$ has a unique positive equilibrium $P^{*}(S^{*}, E^{*}, A^{*}, R^{*}).$ In the following, we will consider the stability of the alcohol-present equilibrium $P^{*}$ and concern about the effect of time delays on the stability of the alcohol-present equilibrium.

Let $s(t) = S(t)-S^{*}, e(t) = E(t)-E^{*}, a(t) = A(t)-A^{*}, r(t) = R(t)-R^{*}$, the linearization of system $\rm(3.1)$ at $P^{*}$ is given by

Further, the Jacobi matrix $J$ at $P^{*}$ can be written as

By computation, the characteristic equation of system $\rm(3.1)$ at $P^{*}$ is given by

where

$l_{1} = 2\mu+p+\varepsilon+m+n(\sigma+1)\beta A^{*}$,

$l_{2} = \sigma\beta^{2}(A^{*})^{2}+[(\mu+n)(1+\sigma)+\sigma p+\varepsilon+m\sigma+1]\beta A^{*}+\sigma\beta(\mu+p)(\mu+\varepsilon)+n(2\mu+p+\varepsilon)+m(2\mu+p+n+\varepsilon)$,

$l_{3} = (m+n)[\sigma\beta^{2}(A^{*})^{2}+(\mu+\varepsilon)\beta A^{*}+\sigma(p+\mu)\beta A^{*}+(\mu+p)(\mu+\varepsilon)]+mn(\sigma+1)\beta A^{*}+n(2\mu+p+\varepsilon)$,

$l_{4} = mn[\sigma\beta^{2}(A^{*})^{2}+\beta A^{*}(\mu+\varepsilon+\sigma p+\sigma\mu)+(\mu+p)(\mu+\varepsilon)]$,

$u_{1} = \beta S^{*}+\sigma\beta E^{*}$,

$u_{2} = \beta^{2}S^{*}A^{*}(\sigma+2)+\sigma\beta^{2}E^{*}A^{*}(2\sigma+1)+\beta(2\mu+p+n+\varepsilon)(S^{*}+\sigma E^{*})+\gamma\delta$,

$u_{3} = (\beta S^{*}+\sigma\beta E^{*})[\sigma\beta^{2}(A^{*})^{2}+(\sigma(\mu+p)+n\sigma+1)\beta A^{*}]+(\mu+p)(\mu+\varepsilon)+n(2\mu+p+\varepsilon)+\sigma^{2}\beta^{2}E^{*}A^{*}(\beta A^{*}+\mu+p+n)\\ +\beta^{2}S^{*}A^{*}(\sigma\beta A^{*}+\mu+\varepsilon+n)+\gamma\delta[\beta A^{*}(1+\sigma)+(2\mu+p+\varepsilon)]$,

$u_{4} = (\beta S^{*}+\sigma\beta E^{*})n[\sigma\beta^{2}(A^{*})^{2}+\beta A^{*}(\mu+\varepsilon+\sigma p+\sigma\mu)+(\mu+p)(\mu+\varepsilon)]+n\beta^{2}S^{*}A^{*}(\sigma\beta A^{*}+\mu+\varepsilon)+np\sigma\beta^{2}(A^{*})^{2}\\ +n\sigma^{2}\beta^{2}E^{*}A^{*}(\beta A^{*}+\mu+p)+\gamma\delta[\sigma\beta^{2}(A^{*})^{2}+\beta A^{*}(\mu+\varepsilon+\sigma\mu+\sigma p)+(\mu+p)(\mu+\varepsilon)]$.

For $\tau > 0$, we assume that Eq (4.21) have root $\lambda(\tau) = \xi(\tau)+i\zeta(\tau)$, where $\xi(\tau), \zeta(\tau)\in \mathbb{R}$. Set ${\tau}^{*}$ be such that $\xi({\tau}^{*}) = 0$, $\zeta({\tau}^{*}) = \zeta_{*}$. Substituting $\lambda({\tau}^{*})$ into Eq (4.21) and separating the real and imaginary parts, we can get the following equations

Taking squares and adding the above equations of Eq (4.22), we get

where $\rho_{1} = l_{1}^{2}-2l_{2}-u_{1}^{2}, \rho_{2} = l_{2}^{2}-u_{2}^{2}-2l_{1}l_{3}+2l_{4}+2u_{1}u_{3}, \rho_{3} = l_{3}^{2}-u_{3}^{2}-2l_{2}l_{4}+2u_{2}u_{4}, \rho_{4} = l_{4}^{2}-u_{4}^{2}.$

(H4): Equation (4.23) has at least one simple positive root.

Theorem 4.6. For Case III, if (H4) does not hold, then the alcohol-present equilibrium $P^{*}$ is locally asymptotically stable for any $\tau_{i} = \tau > 0\; (i = 1, 2, 3)$.

Remark 4.1. If (H4) holds and set $\zeta_{*}$ is the last such root. We will assert that there exists a ${\tau}^{*} > 0$ such that $\xi({\tau}^{*}) = 0, \zeta({\tau}^{*}) = \zeta_{*}$. System $\rm(3.1)$ may exhibit a Hopf bifurcation at $P^{*}$ as ${\tau}$ passes upwards through ${\tau}^{*} = \zeta^{-1}(\zeta_{*}) > 0$ with $\frac{d\xi}{d{\tau}}\mid_{{\tau} = {\tau}^{*}}\neq0.$

Considering the delay of temporarily recovered people is likely to be different from the delay of two classes of susceptible populations, we assume that $\tau_{1} = \tau_{2} = \tilde{\tau}\neq\tau_{3}$ and show the stability in the following case.

4.4.   Case IV: $\tau_{1}=\tau_{2}=\tilde{\tau} > 0, \tau_{3} > 0$

In this case, we only study the stability analysis of $P_{0}$. From the discussion of part (1) in Case I, the characteristic equation Eq (4.1) can be rewritten as

where $f_{1} = m+n > 0, f_{0} = mn > 0$ and $g_{1} = \beta(S_{0}+\sigma E_{0}) > 0, g_{0} = ng_{1} > 0$. Let $\lambda = i\omega_{3}(\omega_{3} > 0)$ be a purely imaginary root of Eq (4.24). It satisfies

Separating the real and imaginary parts, we obtain

where $A = g_{1}\omega_{3}, B = g_{0}$ and $C = f_{0}-\omega_{3}^{2}-h_{0}\cos(\omega_{3}\tau_{3}), D = f_{1}\omega_{3}+h_{0}\sin(\omega_{3}\tau_{3})$. Then, we obtain

Set

where $m_{1} = 0, m_{2} = f_{1}^{2}-g_{1}^{2}-2f_{0}+2h_{0}\cos(\omega_{3}\tau_{3}), m_{3} = 2f_{1}h_{0}\sin(\omega_{3}\tau_{3}), m_{4} = f_{0}^{2}-2f_{0}h_{0}\cos(\omega_{3}\tau_{3})+h_{0}^{2}-g_{0}^{2}$.

We have

Set

where $p_{1} = \frac{m_{2}}{2}, q_{1} = \frac{m_{3}}{4}.$ Notice that there is no square term in the above equation. In order to analyze the roots of Eq (4.27), we cite the results in ^[55] about the existence of positive roots of the fourth degree polynomial equation. Define

Then, we need the following lemma.

Lemma 4.2. (See ^[55]) Suppose that $m_{4}\geq0$.

(i) If $\Delta\geq0$, then Eq (4.27) has positive roots if and only if $\omega_{3}^{(1)} > 0$ and $G(\omega_{3}^{(1)}) < 0$;

(ii) If $\Delta < 0$, then Eq (4.27) has positive root if and only if there exists at least one $\omega^{*}\in \{\omega_{3}^{(1)}, \omega_{3}^{(2)}, \omega_{3}^{(3)}\}$, such that $\omega^{*} > 0$ and $G(\omega^{*})\leq0$.

Notice that $R_{0} = \frac{nc_{1}}{mn-\gamma\delta} < 1$ $(mn > \gamma\delta)$, we can get

Supposing that condition (i) or (ii) of Lemma 4.2 is satisfied, then Eq (4.27) has finite positive roots $\omega_{31}, ..., \omega_{3k}, k\leq4,$ and for every fixed $\omega_{3i}, \; i = 1, ..., k, \; k\leq4,$ there exists a sequence $\{\tau_{3i}^{j}\mid i = 1, ..., k, k\leq4, j = 0, 1, ...\}$ such that Eq (4.25) holds. Let

and $\omega^{*} = \omega_{3i}$, here $i$ corresponding to the minimum value. Then Eq (4.24) has a pair of purely imaginary roots $\pm i\omega^{*}$ when $\tau = \tau_{3}^{*}(\tilde{\tau})$, $\tilde{\tau}\in[0, +\infty)$.

Theorem 4.7. In Case IV, set $\tilde{\tau}\in[0, +\infty)$. If the condition (i) or (ii) of Lemma 4.2 is satisfied, then the alcohol-free equilibrium $P_0$ is asymptotically stable when $\tau_{3}\in\Big[0, \tau_{3}^{*}(\tilde{\tau})\Big].$

Similarly, we can have the following results.

Theorem 4.8. In Case IV, set $\tau_{3}\in[0, +\infty)$. If exists

such that Eq (4.25) holds, then the alcohol-free equilibrium $P_0$ is asymptotically stable when $\tilde{\tau}\in \Big[0, \tilde{\tau}^{*}(\tau_{3})\Big].$

Remark 4.2. We will leave the following studies as open problems:

1) The above research shows that delays do not destroy the globally asymptotical stability of alcohol-free equilibrium $P_{0}$. However, delays can destroy the stability of alcohol-prsent equilibrium $P^{*}$. We assert that there may arise a Hopf bifurcation at $P^{*}$ by increasing the value of one or more delays.

2) The study of other situations is also very valuable. For example: $\tau_{i} > 0, i = 1, 2, 3$ and $\tau_{1}\neq \tau_{2}\neq \tau_{3}$. However, it is well known that the analysis of the characteristic equation is a formidable and challenging task.

5.   Impacts of public health education and delays on alcohol control

In the process of establishing the model, the public health education and three delay factors are introduced. In this section, the impact of these factors on the number of alcoholics will be discussed in detail.

First of all, we choose a set of values of parameters (parameter values are mainly taken from ^[5,41]):

$\mu = 0.25, \; \beta = 0.8, \; a_{1} = 0.01, \; a_{2} = 0.01, \; q = 0.6, \; \Lambda = 30.496, \; \varepsilon = 0.2, \; \delta = 0.1, \; \gamma = 0.7,$

$\xi = 0.04, \; p = 0.1, \; \tau_{1} = \tau_{2} = \tau_{3} = 0.5.$

Due to the greater the positive impact of media reporting measures on the population, the smaller the value of $\sigma$ is. Figure 2 shows the number of alcoholics $A(t)$ decreased with the decrease of the value of the media coverage impact factor $\sigma$. So, it is clear to see that the public health education (media reporting measure) is beneficial to control the spread of alcohol problems.

Next, the simulation results for different time delays are shown in Figure 3. Let $\sigma = 0.2$, $\tau_{i} = 0.01, 1, 2, 3 (i = 1, 2, 3)$. We can see that the change trend of alcoholics $A(t)$ for different values of $\tau_{1}$ in Figure 3(a). That is, as the value of $\tau_{1}$ increases, the value of $A(t)$ increases rapidly in a short time. But the smaller the value of $\tau_{1}$ is, the smaller the corresponding peak value of $A(t)$ is. Moreover, the larger the value of $\tau_{1}$, the longer it takes for the value of $A(t)$ to reach the highest peak. This indicates that the higher the value of $\tau_{1}$, the more unfavorable it is to control the problem of drinking. From Figure 3(b), as the value of $\tau_{2}$ increases, the value of $A(t)$ decreases rapidly in a short time. This shows that the higher the value of $\tau_{2}$, the more favorable the control of alcohol problem. From Figure 3(c), in the long term, as the value of $\tau_{3}$ increases, the value of $A(t)$ decreases. This shows that the higher the value of $\tau_{3}$, the more favorable the control of alcohol problem.

6.   Sensitivity analysis

In this section, we will use the Latin hypercube sampling method in ^[56,57,58] and Partial rank correlation coefficients method to study the sensitivity analysis of each parameter on the basic reproduction number $R_{0}$ and the number of alcoholics $A(t)$. By discussing the sensitivity analysis of each parameter, we can control the key parameters and take corresponding measures to control alcohol problems.

First, we will explore the effect of various parameters on the basic reproduction number $R_{0}$, see Figure 4. Parameter values are mainly taken from ^[5,41] and given by Table 1.

It is assumed that the parameters of the system obey normal distribution, and 2000 samples are taken within the 95% confidence interval (taking the upper and lower as 0.02). In what follows, we perform the numerical study. Namely, the PRCC-values of $R_0$ with respect to each parameter of the system are given. According to the order of parameters $\mu, \beta, a_{1}, a_{2}, q, \Lambda, \sigma, \varepsilon, \delta, \gamma, \xi, p$, the corresponding PRCC-values are:

Similarly, the P-values of $R_0$ with respect to each parameter of the system are as follows. According to the order of parameters $\mu, \beta, a_{1}, a_{2}, q, \Lambda, \sigma, \varepsilon, \delta, \gamma, \xi, p$, the corresponding P-values are:

Parameters with the highest (or lowest) PRCC values have the largest positive (or negative) impact on $R_0$. On the contrary, $R_0$ is more sensitive to the parameter with smaller P-value. Figure 4(a), (b) depict the parameters that have the most positive influence on $R_0$ are $\beta, \; q, \; \Lambda, \; \sigma$ and $\delta$. While the parameters with the most negative impact on $R_0$ are $\mu, \; a_{1}, \; a_{2}, \; \varepsilon, \; \gamma, \; \xi$ and $p$. It can be seen that alcohol transfer rate $\beta$ has the greatest positive correlation with $R_0$, and population transfer rate caused by media reports $p$ has the greatest negative correlation with $R_0$. Moreover, the sensitivity of $a_{2}$ is the weakest among all parameters. Practical significance: by strengthening the news real-time report and cumulative report on alcohol related hazards, we can reduce the value of $\sigma$ and $\delta$, and increase the values of $\varepsilon$, $\xi$ and $p$ to make the value of $R_0$ decreases to $R_0 < 1$. so as to control the number of alcoholics in the population. Second, we study the sensitivity analysis of parameters to $A(t)$. Focusing on the impact of media reporting measures on alcoholism. Therefore, the following mainly considers several parameters that are directly affected by media reports as shown in Figure 5. Figure 5(a) shows the influence of parameters $\mu$ and $a_{1}$ on $A(t)$ is always negatively correlated. And, the influence of parameter $a_{2}$ on $A(t)$ is positively correlated at first, and as time goes by, the influence of $a_{2}$ on $A(t)$ becomes negatively correlated because the relapsed person drinks again. However, the influence of parameter $\varepsilon$ on $A(t)$ is negatively. Because with the passage of time, the influence of media reports on people will weaken or even disappear, which also shows that it is necessary to carry out regular reports on the dangers of alcoholism to strengthen people's health awareness. In particular, Figure 5(b) shows the PRCC-values of parameter $\sigma$ or $p$ will be reduced to zero in the later stage, with slight positive and negative fluctuations. Research shows that the influence of media reports has saturation effect. This is also consistent with the reality. Regular media coverage of the dangers of alcoholism can control alcoholism in the community.

Remark 6.1. Among the abundant literature on sensitivity measures, the Sobols indices have received much attention since they provide accurate information for most models. Polynomial chaos-based Sobols indices are used in ^[60,61,62]. Recently, using deterministic, frequentist and Bayesian approaches, Calatayud et al. ^[63] studied the effect of social behavior on the increase of the obesity epidemic. And several sensitivity analyses are also conducted. These approaches are interesting alternatives for sensitivity analyses. We will leave these interesting work for the future.

7.   Numerical simulations

In this section, some numerical results of system $\rm(3.1)$ are provided to substantiate the analytic results obtained in this paper.

In the following, we fix parameters:

$\mu = 0.25, \; \beta = 0.1, \; a_{1} = 0.01, \; a_{2} = 0.01, \; p = 0.1, \; q = 0.6, \; \Lambda = 3.496,$

$\sigma = 0.2, \; \varepsilon = 0.2, \; \delta = 0.1,$$\gamma = 0.5, \; \xi = 0.4,$

7.1.   For case I: $\tau_{1}=\tau_{2}=0, \tau_{3} > 0$

Choose $\; \tau_{1} = \tau_{2} = 0, \tau_{3} = 0.5$, we have $R_{0} = 0.2478 < 1$ and meet the previous conditions. Thus the alcohol-free equilibrium

is globally asymptotically stable. The solutions of $P_0$ are shown in Figure 6(a).

Other parameters remain unchanged and let $\; \beta = 0.8, \gamma = 0.7,$ we have $R_{0} = 1.7531 > 1$, then the alcohol-free equilibrium $P_0$ is unstable and the alcohol-present equilibrium

is locally asymptotically stable. The solutions of $P^{*}$ are shown in Figure 6(b).

7.2.   For case II: $\tau_{1}=\tau_{2} > 0, \tau_{3}=0$

Choose $\mu = 0.7, \; \tau_{1} = \tau_{2} = 0.7, \tau_{3} = 0$, we have $R_{0} = 0.1792 < 1$ and the alcohol-free equilibrium

is globally asymptotically stable. The solutions of $P_0$ are shown in Figure 7(a).

Choose $\mu = 0.3, \; \beta = 0.8$, we get $R_{0} = 1.9437 > 1$ and the alcohol-present equilibrium

is locally asymptotically stable. The solutions of $P^{*}$ are shown in Figure 7(b).

7.3.   For case III: $\tau_{1}=\tau_{2}=\tau_{3} > 0$

Choose $\mu = 0.25, \; \beta = 0.05, \; \tau_{1} = \tau_{2} = \tau_{3} = 0.5$, we have $R_{0} = 0.1239 < 1$ and the alcohol-free equilibrium

is globally asymptotically stable. The solutions of $P_0$ are shown in Figure 8(a).

Choose $\mu = 0.25, \; \beta = 0.7$, we get $R_{0} = 1.7363 > 1$ and the alcohol-present equilibrium

is locally asymptotically stable. The solutions of $P^{*}$ are shown in Figure 8(b).

7.4.   For case IV: $\tau_{1}=\tau_{2} > 0, \tau_{3} > 0$

Choose $\; \tau_{1} = \tau_{2} = 0.1, \; \tau_{3} = 0.9$, we have $R_{0} = 0.0423 < 1$ and the alcohol-free equilibrium

is asymptotically stable. The solutions are shown in Figure 9(a).

Choose another set of values $\; \tau_{1} = \tau_{2} = 0.2, \; \tau_{3} = 0.8$, we get $R_{0} = 0.1025 < 1$ and the alcohol-free equilibrium

is also asymptotically stable. The solutions are shown in Figure 9(b).

Through the above theoretical analysis and numerical simulation, it can be seen that different time delays have little effect on the global asymptotic behavior of the alcohol-free equilibrium point $P_0$, but have a greater impact on the dynamic behavior of the alcoholic equilibrium point $P^{*}$.

8.   Conclusions and discussions

There is a large literature (see ^{[15,16,18,19,20,21,22,23,24,25,26,27]} and ^{[28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]}) that has been used for studying alcohol problems and describing the effects of awareness programs by media on the infectious disease or alcohol consumption behavior. From ^[41], we know that the global stability of system $\rm(2.1)$ without delays has been proved by using the Lyapunov function. That is, if $R_{0} < 1$, all solutions converge to the alcohol-free equilibrium $P_0$, and the alcohol problems disappear eventually; if $R_{0} > 1$, the alcohol-present equilibrium $P^{*}$ is globally stable, i.e., the alcohol problems will persist in the population and the number of problem drinkers tends to a positive constant. However, most of these models are ODE models and do not incorporate the effects of the time delay. Recently, Huo et al. ^[49] and Ma et al. ^[50] incorporated one delay into alcoholism model and studied it's stability via the delay differential equations. Their results showed that the stability of system be destroyed by the delay and the system occurs a Hopf bifurcation at the alcohol-present equilibrium under certain conditions.

Compared to previous models, this paper considered an alcoholism model with health education and three delays. The main goal is to analyze the global dynamics of the model and discuss the impacts of health education and three delays on the alcohol consumption behavior. A key novelty of our model is that we do not only consider the impact of health education but also introduce three delays which describe the time needed that a susceptible or temporarily individual becomes an infectious alcohol user. It is necessary and practical significance to consider the dynamic behavior about delayed alcoholism model. A question is how the three delays affect the dynamics of the system $\rm(3.1)$. This is the main purpose of this paper.

Specifically, we discussed a FDE system $\rm(3.1)$ in four different cases respectively. If $R_{0} < 1$, our results showed that the alcohol-free $P_0$ is still globally asymptotically stable after incorporating time delays, which dose not change the dynamic of $P_0$. However, if $R_{0} > 1$, our results shows that incorporation of three delays will change the stability of the alcohol-present equilibrium $P^{*}$. This means that difference of the time delay between two classes of susceptible populations and temporarily recovered people can affect the spread and control of alcohol problems greatly. The higher the value of $\tau_{2}$ or $\; \tau_{3}$, the more conducive to the prevention and control of alcoholism. On the contrary, the higher the value of $\tau_{1}$, the more unfavorable the prevention and control of alcoholism. Under the point of view of the alcohol control, our results show that time delays have important influence in controlling alcohol problem. For condition (H1), that is, the transmission rate $\beta$ plays an essential role in our study. Public health education can reduce the transmission rate $\beta$. Our model improves and generalizes the existing results.

In this article, we only considered several special cases. For the alcohol-present equilibrium $P^{*}$, we assert that there may arise a Hopf bifurcation at $P^{*}$ by increasing the value of delay. The study of other situations is more conform to reality the actual phenomenon. However, it is well known that the analysis of the characteristic equation is a formidable and challenging task. We will leave these studies as open problems.

For simplicity, the random variability is normally distributed in this paper. In fact, the value of the model key parameters can be obtained by using statistic data will make our paper more perfect. From Brauer ^[59], we know that a example (The Great Plague in Eyam): the village of Eyam near Sheffield, England suffered an outbreak of bubonic plague in 1665–1666. The actual data for the Eyam epidemic are remarkably close to the predictions of the model. Our model coefficients should be fixed from medical reports and the delays in model may have an estimation. At present, we lack statistical data on alcohol abuse. So that the parameters are given some values to test the theoretical results in Sections 5 and 7.

In Section 6, we use the Latin hypercube sampling method and Partial rank correlation coefficients method to study the sensitivity analysis. We know that the Sobol's index is a sensitivity grading method based on the proportion of parameters contribution to output variance. Among the abundant literature on sensitivity measures, the Sobol's indices have received much attention since they provide accurate information for most models. Polynomial chaos-based Sobol's indices are used in ^[60,61,62]. By deterministic, frequentist and Bayesian approaches, several sensitivity analyses are also conducted in ^[63]. These alternative methods of sensitivity analysis is effective. At the same time, compared with ordinary differential equation model, stochastic functional differential equation model (see ^{[64,65,66,67]}) will be more suitable for the actual situation. We will leave these interesting work for the future.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11861044 and 11661050), and the HongLiu first-class disciplines Development Program of Lanzhou University of Technology. The authors acknowledge the anonymous referees and editors very much for their fruitful comments and valuable suggestions which made this paper much improved.

Conflict of interests

The authors declare that there is no conflict of interests regarding the publication of this paper.