Mathematical Modeling and Analysis of an Alcohol Drinking Model with the Influence of Alcohol Treatment Centers

In this paper, we present a continuous mathematical model PMHTTQ of alcohol drinking with the influence of private and public addiction treatment centers. We study the dynamical behavior of this model and we discuss the basic properties of the system and determine its basic reproduction number R0. We also study the sensitivity analysis of model parameters to know the parameters that have a high impact on the reproduction number R0. *e stability analysis of the model shows that the system is locally as well as globally asymptotically stable at drinking-free equilibrium E0 when R0 ≤ 1. When R0 > 1, drinking present equilibrium E∗ exists and the system is locally as well as globally asymptotically stable at alcohol present equilibrium E∗.


Introduction
Alcoholism is a social phenomenon that affects all social classes and people with different educational levels and age groups and can have an impact on many aspects of life. It is a chronic disorder that causes the person to drink uncontrollably; as a result, it can destroy people's relationships and make them less productive as well as harm their physical and mental health. Long-term abuse can also cause severe damage to the brain and liver. Alcohol has several social, economic, and health effects on the individual and society as a whole including heart disease, waste of money, poverty, crime, family disintegration, and liver disease and also causes ulcers, diabetes complications, sexual problems, birth defects, bone loss, vision problems, increased risk of cancer, and suppressed immune function which need treatment in hospitals and private or public addiction treatment centers. Alcoholism remains one of the most frequent and important topics discussed by the world community due to its dire consequences [1].
According to the World Health Organization (WHO), in its report "World Status Report on Alcohol and Health 2018," alcohol addiction causes the death of three million people worldwide every year, which represents one in every 20 deaths. For example, alcohol consumption kills more than AIDS, tuberculosis, and violence combined [1].
In 2016, alcohol has killed around 3 million people worldwide, compared with 3.3 million in 2012. ree quarters of the deaths are in men. According to these figures, there is a clear reduction in the number of deaths caused by alcohol addiction. In the European region and in the Americas region, the harmful consequences including illness and injury resulted from the use of alcohol are detrimental in comparison with other regions. Several diseases, more than 200, are related to alcohol consumption. 28% of the 3 million deaths ascribed to alcohol were linked to traffic accidents, violence, suicides, and other violent acts, 21% to digestive disorders, and 19% to cardiovascular diseases. e other cases of deaths are resulted from infectious diseases, cancers, mental disorders, and other health problems [1].
Many studies and research in social, medical, and political sciences have focused on this topic and other related topics ( [2][3][4] and the references cited therein). But the mathematical studies and research on this topic are still limited and most of them have focused on the statistical aspect of the phenomenon ( [5][6][7][8]).
Mathematical models can be used to analyze the spread of infectious diseases or the social behavior of individuals [9][10][11][12][13][14][15][16][17][18][19][20]. As regard to drinking, several different mathematical models have been formulated and studied to help in reducing the number of drinkers [17,[21][22][23][24]. Agrawal et al. [22] developed a nonlinear SHTR mathematical model of alcohol abuse with a nonlinear incidence rate. e stability analysis of the model that they proposed shows that the system is locally asymptotically stable at alcohol-free equilibrium E 0 when R 0 ≤ 1. When R 0 > 1, alcohol equilibrium E * exists and the system becomes locally asymptotically stable at E * and E 0 becomes unstable. Manthey et al. [23] proposed a mathematical model to study the dynamics of campus drinking as an epidemiological model. According to their results, the reproductive numbers were not sufficient to predict whether drinking behavior would persist on campus and that the pattern of recruiting new members played a significant role in the reduction of campus alcohol problems. Sharma et al. [17] developed a mathematical model of alcohol abuse and discussed the existence, local, global stability of drinking-free, endemic equilibria, and sensitivity analysis of R 0 . ey demonstrated that backward bifurcation can occur when R 0 � 1. Sharma et al. [24] developed a two-stage (four compartments) model for youths with serious drinking problems and their treatment having admitted and having not admitted. ey also analyzed the stability of all the equilibria. Adu [21] used a nonlinear SHTR mathematical model to study the dynamics of drinking epidemic. He divided the population into four compartments: nondrinkers (S), heavy drinkers (H), drinkers in treatment (T), and recovered drinkers (R). He discussed the existence and stability of drinkingfree and endemic equilibria.
Besides these works, we will study the analysis of a mathematical alcohol model PMHT r T p Q with the following additions: (i) Compartment T p represents the number of the poor heavy drinkers who join public addiction treatment centers, which may not have sophisticated equipment and do not provide high quality treatment services in comparison with the private ones, especially in the developing countries, and those individuals who do not have the financial capacity to go to private centers (ii) Compartment T r represents the number of the rich heavy drinkers who join private addiction treatment centers that encompass special facilities and provide advanced treatment for such individuals with sufficient financial capacity (iii) e death rate induced by the heavy drinkers δ roughout this research, we present a nonlinear mathematical alcohol model PMHT r T p Q with private and public addiction treatment centers that describes the dynamics of drinkers' classes with the influence of the private and public addiction treatment centers. Hence, the awareness programs and the addiction treatment centers of alcohol drinking will alert the susceptible individuals and heavy drinkers so they decline drinking. We assume that the heavy drinkers can recover from heavy alcohol drinking due to media, health reasons, joining addiction treatment centers, prohibition, and hike tax on alcohol beverages. A fraction of heavy drinkers' population will join the private addiction treatment centers, whereas the rest will join the public addiction treatment centers. e goal is to analyze the impact of the addiction treatment centers and awareness programs on the drinkers. We also analyze the effect of addiction treatment centers on system stability. Our results show that addiction treatment centers and awareness programs are an effective measure in controlling heavy drinking. By using Routh-Hurwitz criteria and constructing Lyapunov functions, the local and the global stability of alcohol-free equilibrium and alcohol equilibrium are obtained. We also study the sensitivity analysis of model parameters to know the parameters that have a high impact on the reproduction number R 0 . e drinker classes of this model are divided into six compartments: potential drinkers (P), moderate drinkers (M), heavy drinkers (H), the rich heavy drinkers who join private addiction treatment centers (T r ), the poor heavy drinkers who join public addiction treatment centers (T p ), and quitters of drinking (Q). e rest of the paper is organized as follows. In section 2, the formulation of the model and some basic properties are derived. In section 3, equilibria of the proposed model are obtained and their stability discussed. In section 4, the global stability of equilibrium point is discussed. e problem of sensitivity's parameters is discussed in section 5. Some numerical simulations and discussions are given in section 6. Lastly we give the conclusion of the whole paper.

Model Description.
We propose a continuous model PMHT r T p Q to describe the interaction between drinkers classes. e population is divided into six compartments: potential drinkers P(t), moderate drinkers M(t), heavy drinkers H(t), rich heavy drinkers T r (t), poor heavy drinkers T p (t), and quitters of drinking Q(t).
e graphical representation of the proposed model is shown in Figure 1.
We consider the following system of six nonlinear differential equations: where ≥ 0, and Q(0) ≥ 0 are the given initial states.

Potential Drinkers
e potential drinkers P(t) represents individuals whose age is over adolescence and adulthood and may become drinkers. is compartment is increased by the recruitment rate denoted by b and decreased by an effective contact with the moderate drinkers at β 1 rate and natural death μ. It is assumed that potential drinkers can acquire drinking behavior and can become moderate drinkers through effective contact with moderate drinkers in some social occasions such as weddings, celebrating graduation ceremonies, weekend parties, and end of the year celebration. In other words, it is assumed that the acquisition of a drinking behavior is analogous to acquiring disease infection.

Moderate Drinkers
e compartment M is composed of moderate drinkers who can control their consumption during some events and occasions or in a way that is unapparent to their social environment. is category of drinkers does not face any problems or negative consequences. Friends or family do not complain about their intake of alcohol. A moderate drinker does not think about drinking very often or often feel a need to drink. ey may go out to have a few drinks and are able to handle their alcohol consumption without experiencing a loss of control. Alcohol does not dominate their thoughts and they do not need to set limits when they drink. ey are not prone to extreme mood swings, fighting, or being violent. It is increased by potential drinkers who turn to be moderate drinkers at β 1 rate. is compartment is decreased when moderate drinkers become heavy drinkers at a rate β 2 and also by natural death at rate μ.

Heavy Drinkers
e compartment H is composed of heavy drinkers suffering from addiction to alcohol. When an individual becomes an alcoholic, they face a great difficulty to control or set limits for their consumption. e majority of alcoholics begin as potential drinkers and then turn to moderate drinkers. Alcohol seems to exert a control on the alcoholic's life. eir job, their family, social circle, and health are all endangered. Despite these negative consequences, the alcoholic is unable to quit drinking. e alcoholics may begin to disclaim that they have a problem; this disclaim can make it even more difficult for the person to get help. Alcohol addiction is considered to be a disease; it changes chemicals in the addict's brain and has made alcohol the most important thing in their life. At the time a person is an alcoholic, they will usually need to get help at a rehab to overcome their addiction. is compartment becomes larger as the number of heavy drinkers increases by the rate β 2 and decreases when some of them give up drinking at a rate α 3 as well as when they join treatment centers whether private or public at rates α 1 and α 2 . In addition, this compartment decreases by natural death μ and due to deaths caused by diseases resulted from excessive alcohol intake at a rate δ.

Rich Heavy Drinkers
e compartment T r contains the number of heavy drinkers who take advantage of their financial potentials to join private treatment centers of alcohol addiction that are very often well equipped and provide good and quality services. is compartment is increased by the rate α 1 and decreased by the rates c 1 and μ.

Poor Heavy Drinkers
e compartment T p represents the number of heavy drinkers who join public treatment centers of alcohol addiction which may not provide advanced treatment and that b P International Journal of Mathematics and Mathematical Sciences are marked by a shortage of equipment and low-quality services, especially in the developing countries, and it also contains individuals who do not have the financial capacity to join the private centers. is compartment is increased by the rate α 2 and decreased by the rates c 2 and μ.
e compartment Q encompasses the individuals who quit drinking. It is increased with the recruitment of individuals who have been treated in treatment centers of alcohol addiction at rates α 1 and α 2 . It is also increases at the rate α 3 of those who quit alcohol without resort to treatment centers and decreases at the rate μ due to natural deaths. e total population size at time t is denoted by

Invariant Region.
It is necessary to prove that all solutions of system (1) with positive initial data will remain positive for all times t > 0. is will be established by the following lemma.

Lemma 1. All feasible solution P(t), M(t), H(t), T r (t), T p (t), and Q(t) of system equation (1) are bounded by the region:
Ω � P, M, H, T r , T p , Q ∈ IR 6 + : Proof. From system equation (1), implies that and it follows that where N(0) is the initial value of total number of people; thus, then Hence, for the analysis of model (1), we get the region which is given by the set Ω � P, M, H, T r , T p , Q ∈ IR 6 + : which is a positively invariant set for (1), so we only need to consider the dynamics of system (1) on the set Ω nonnegative of solutions.
It is assumed that We multiply equation (15) by exp( t 0 A(s)ds) and we find which implies that erefore, Taking integral with respect to s from 0 to t, we obtain

International Journal of Mathematics and Mathematical Sciences
Multiplying equation (20) by exp(− t 0 A(s)ds), we obtain en, we can obtain So, the solution P(t) is positive. Similarly, from the second equation of system (1), we have where Similarly, from the third, fourth, and fifth equations of system (1), we have erefore, we can see that P(t) > 0, M(t) > 0, H(t) > 0, T r (t) > 0, T p (t), and Q(t) > 0 ∀t ≥ 0, and this completes the proof.
□ e first three equations in system (1) are independent of the variables T r , T p , and Q. Hence, the dynamics of equation system (1) is equivalent to the dynamics of the equation system:

Equilibrium Point.
e standard method is used to analyze model (26). In this model, there are two equilibrium points, that is, drinking-free equilibrium point and drinking present equilibrium point. e equilibrium points are found by setting the right-hand side of equations (1)-(3) equal to zero. e drinking-free equilibrium E 0 ((b/μ), 0, 0) is achieved in the absence of drinking (M � H � 0). e drinking present equilibrium E * (P * , M * , H * ) is achieved when drinkers exist (M ≠ 0 and H ≠ 0), where R 0 is the basic reproduction number that measures the average number of new drinkers generated by single drinker in a population of potential drinkers. e value of R 0 will indicate whether the epidemic could occur or not. e reproduction basic number can be determined by using the next generation matrix method formulated in [25].

Local Stability
Analysis. Now we proceed to study the stability behavior of equilibria E 0 and E * .

Proof.
e Jacobian matrix at E is given by e Jacobian matrix for the drinking-free equilibrium is given by erefore, all the eigenvalues of the characteristic equation J(E 0 ) are clearly real and negative if R 0 < 1. □ We conclude that the drinking-free equilibrium is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

Drinking Present Equilibrium.
In this section, we analyze the local stability of drinking present equilibrium.
To find the drinking present equilibrium of the system of equation (26) we consider dP(t)/dt � 0, dM(t)/dt � 0, and dH(t)/dt � 0 provided that at least one of the infected compartments is nonzero. We evaluate the equilibrium of system (22) by setting the right-hand side of equation of system (26) equal to zero and then solve for P * , M * , and H * .
From the second equation in system (26), we have From the first equation in system (22), we have Also, equation (3) in system (26) gives Let the following theorem analyzes the local stability of the drinking present equilibrium when R 0 > 1.

Theorem 3.
e drinking present equilibrium E * is locally asymptotically stable if R 0 > 1 and unstable otherwise.
us, the drinking present equilibrium E * of system (26) is locally asymptotically stable. (26) is globally asymptotically stable, we use the Lyapunov function theory for both the drinking-free equilibrium and the drinking present equilibrium. First, we present the global stability of the drinking-free equilibrium E 0 when R 0 ≤ 1.

Theorem 4.
e drinking-free equilibruim E 0 is globally asymptotically stable if R 0 ≤ 1 and unstable otherwise.

International Journal of Mathematics and Mathematical Sciences
So, (dV/dt) ≤ 0 if R 0 ≤ 1. Furthermore, (dV/dt) � 0 if and only if M � 0. Hence, by LaSalle's invariance principle [26], E 0 is globally asymptotically stable. e final result of the global stability of E * in this section is as follows.

Theorem 5.
e drinking-present equilibrium point E * is globally asymptotically stable if R 0 > 1.
Proof. Consider the Lyapunov function V: where c 1 and c 2 are positive constants to be chosen later and Γ � (P, M) ∈ Γ/P > 0, M > 0 { }. en, the time derivative of the Lyapunov function is given by en, the time derivative of the Lyapunov function is given by For c 1 � c 2 � 1, we have Also, we obtain dV(P, M) dt Hence, by LaSalle's invariance principle [26], the drinking present equilibrium point E * is globally asymptotically stable on Γ. □

Sensitivity Analysis of R 0
Sensitivity analysis is commonly used to determine the model robustness to parameter values, that is, to help us know the parameters that have a high impact on the reproduction number R 0 .
Using the approach in Chitnis et al. [27], we calculate the normalized forward sensitivity indices of R 0 . Let denote the sensitivity index of R 0 with respect to the parameter m. We obtain From the abovementioned discussion, we observe that the basic reproduction number R 0 is most sensitive to changes in β 1 . If β 1 increases R 0 will also increase with the same proportion, and if β 1 decreases in the same proportion, μ and β 2 will have an inversely proportional relationship with R 0 . So, an increase in any of them will bring about a decrease in R 0 ; however, the size of the decrease will be proportionally smaller. Recall that μ is the natural death rate of the population. Given R 0 's sensitivity to β 1 , it seems sensible to focus efforts on the reduction of β 1 . In other words, this sensitivity analysis tells us that prevention is better than cure. Efforts to increase prevention are more effective in controlling the spread of habitual drinking than efforts to increase the numbers of individuals accessing treatment (Table 1).

Numerical Simulations
In this section, we illustrate some numerical solutions of model (1) for different values of the parameters. e resolution of system (1) was created using the Gauss-Seidel-like implicit finite-difference method developed by Gumel et al. [28], presented in [29] and denoted the GSS1 method. We use the following different initial values such that P + M + H + T r + T p + Q � 1000.
We use the following parameters in Table 2.
We begin by a graphic representation of the drinkingfree equilibrium E 0 , and we use the same parameters and different initial values given in Table 2, R 0 � 0.37 and R 0 < 1.
From these figures, using the different values of initial variables P 0 , M 0 , H 0 , T r 0 , T p 0 , and Q 0 , we obtained the following remarks ( Figure 2): (i) e number of potential drinkers increases and approaches the number P 0 � 1000 (see Figure 2(a)) (ii) e number of moderate drinkers decreases and approaches zero (see Figure 2(b)) (iii) e number of the heavy drinkers increases at first, after that it decreases and approaches zero (see Figure 2(c)) (iv) e number of the poor heavy drinkers decreases and approaches zero (see Figure 2(d)) (v) e number of the rich heavy drinkers decreases and approaches zero (see Figure 2(e)), (vi) e number of the quitters of drinking decreases and approaches zero (see Figure 2(f )) International Journal of Mathematics and Mathematical Sciences Table 2: e academic parameters used for model (1). 600  800  100  300  500  50  150  200 65   Figure 2: When R 0 < 1, the drinking-free equilibrium E 0 is globally asymptotically stable. Table 3: e academic parameters used for model (1).
We use the following parameters in Table 3. Also, we begin by a graphic representation of the drinking present equilibrium E * and we use the same parameters and different initial values given in Table 3, R 0 � 3.66 and R 0 > 1.
From these figures, using the different values of initial variables P 0 , M 0 , H 0 , T r 0 , T p 0 , and Q 0 , we obtained the following remarks ( Figure 3): (i) e number of potential drinkers increases at first, then it decreases slightly and approaches the value P * � 273 (see Figure 3(a)) (ii) e number of moderate drinkers decreases rapidly at first, then increases slightly and approaches the value M * � 230 (see Figure 3(b)) (iii) e number of heavy drinkers increases and approaches the value H * � 461 (see Figure 3(c)) (iv) e number of poor heavy drinkers decreases and approaches the value T p * � 6.89 (see Figure 3(d)) (v) e number of rich heavy drinkers decreases and approaches the value T r * � 7.02 (see Figure 3(e)) (vi) e number of quitters of drinking decreases and approaches the value Q * � 7.03 (see Figure 3(f )) erefore, the solution curves to the equilibrium E * (P * , M * , H * , T p * , T r * , Q * ) when R 0 > 1. Hence, model (1) is globally asymptotically stable.

Conclusion
In this paper, we have presented a continuous mathematical model PMHT r T p Q of alcohol drinking with the influence of Quitters of drinking (Q) (f ) Figure 3: When R 0 < 1, the drinking-present equilibrium E * is globally asymptotically stable.
private and public addiction treatment centers and the dynamical behavior of the model is studied. We have found R 0 � (β1/μ + β 2 ) as basic reproduction number of system (1), which helps us to determine the dynamical behavior of the system. We also studied the sensitivity analysis of model parameters to know the parameters that have a high impact on the reproduction number R 0 . We used the stability analysis theory for nonlinear systems to analyze the mathematical drinking model and to study both the local and global behavior of drinking dynamics. Local asymptotic stability for the drinking-free equilibrium E 0 can be obtained, if the threshold quantity R 0 ≤ 1. On the other hand, if R 0 > 1, then the alcohol present equilibrium E * is locally asymptotically stable. A Lyapunov function was used to show global stability of E 0 . E 0 is globally asymptotically stable if R 0 ≤ 1. Also a Lyapunov function was used to show global stability of E * . E * is globally asymptotically stable if R 0 > 1.

Data Availability
e disciplinary data used to support the findings of this study have been deposited in the Network Repository (http://www.networkrepository.com).