This paper presents a deterministic SATQ-type mathematical model (including susceptible, alcoholism, treating, and quitting compartments) for the spread of alcoholism with two control strategies to gain insights into this increasingly concerned about health and social phenomenon. Some properties of the solutions to the model including positivity, existence and stability are analyzed. The optimal control strategies are derived by proposing an objective functional and using Pontryagin’s Maximum Principle. Numerical simulations are also conducted in the analytic results.

1. Introduction

Alcoholism, also known as alcohol dependence, is a disease that includes the desire for alcohol and continuing to drink it despite its negative effect on individual’s health, relationships, and social status [1]. Similar to all other drug addictions, alcoholism can be regarded as a treatable disease. The World Health Organization estimates that about 140 million people throughout the world suffer from alcohol dependence with related problems, such as being sick, losing a job, among a host of other things [2]. Particularly, young people’s alcoholism problem is a major concern to public health. US surveys indicate that approximately 90% of college students have consumed alcohol at least once [3], and more than 40% of college students have engaged in binge drinking [4, 5]. Unfortunately, the biological mechanisms underpinning alcoholism are not known; however, risk factors include social environment, stress, mental health, genetic sensitivity, age, ethnic group, and sex [6, 7]. Long-term alcohol abuse will produce negative changes in the brain such as tolerance and physical dependence. The subtle changes make it difficult for the alcoholics to stop drinking and result in alcohol withdrawal symptoms upon discontinuation of alcohol consumption. Alcohol damages almost all parts of the body and contribute to a number of human diseases including but not limited to liver cirrhosis, pancreatitis, heart disease, and sexual dysfunction and can eventually be deadly [8]. Damage to the central and peripheral nervous systems can take place from sustained alcohol consumption [913].

Although alcoholism is becoming more and more dangerous and serious as well as a widespread social phenomenon, only much less work has been done in the mathematical modelling of alcoholism as a growing health problem, including a few studies which offered some mathematical approaches to understand the growing burden of alcoholism [10, 1419]. In [10], a SIR-type model was proposed; the authors used standard contact rate between susceptibles and alcoholism, getting alcoholism reproductive number and discussing the existence and stability of two equilibria. In [14], a framework where drinking was modeled as a socially contagious process in low- and high-risk connected environments was introduced; they found that high levels of social interaction between light and moderate drinkers in low-risk environments can diminish the importance of the distribution of relative drinking times on the prevalence of heavy drinking. In [15], neurophysiological examinations of 100 long-term alcohol dependent patients, who were having neuropsychiatric treatment, showed symptoms of polytopic damage of the peripheral and central nervous system. The results showed that for recognition of the damage an extensive diagnostic programme must be used. In [16], the authors considered a kind of binge drinking model with two equal infectivity drunk states; mathematical analyses established that the global dynamics of the model were determined by the basic reproduction number. In [17], the authors modified the model from [16]; that is, they considered different infectivity of two drunk states, and a SEIR-type model of alcoholism was thus presented, in which two alcohol related states were involved, namely, no alcohol dependent consumers and alcohol dependent consumers . In [18], the authors formulated a deterministic model for evaluating the impact of heavy alcohol drinking on the reemerging gonorrhea epidemic, and both analytical and numerical results were provided to ascertain whether heavy alcohol drinking had an impact on the transmission dynamics of gonorrhea. The approach of the literature [18] was very meaningful, since it provided a new direction of thinking when the cross-infection between alcoholism and other pathological diseases occurs. In recent monograph [19], the authors also proposed a SIR-type model to investigate alcohol abuse phenomenon and generated some useful insights; for example, the basic reproductive number was not always the key to controlling drinking within the population. For other papers that study the model of giving up smoking or quitting drinking, please see [20, 21] and references cited therein.

As living standard and health awareness get improved, more and more people who fall into binge drinking state are actively seeking the quitting alcoholism measures and treatment methods [1, 11, 22]. In [22], treatment strategy was introduced into a simple SIR-type alcoholics quitting model, in which the authors used the bilinear incidence to depict the “infection" between the occasional drinkers and problem drinkers . Motivated by some aforementioned documents [10, 19, 22], in this paper, we will formulate a more reasonable alcoholics quitting model. The fact that our model is more reasonable is embodied from the following three aspects.(1)Taking into account that alcoholism is a widespread social phenomenon, so the standard incidence is superior to bilinear incidence when we portray the relationship between the alcoholism and the susceptibles during the course of infection. While in [22], the authors adopted bilinear incidence, we will adopt standard incidence in this paper.(2)Since alcohol is harmful to health, moreover, as we all know, alcoholism is treatable if we can take approximate measure in time, for example, artificial isolation from alcoholisms, medications, persuasion, and education programing on alcoholism. So it is necessary to take effective measures to avoid alcohol or to treat after alcoholism. Documents [10, 19] have not considered these aspects.(3)Since there is effective prevention and treatment in describing the phenomenon of alcoholism, there are some people who will never drink due to successful prevention or some people who no longer drink after successful treatment. Therefore, when we formulate the model in this paper, it's reasonable to introduce a new compartment Q, the people in which will never drink for ever. Obviously, the models of [10, 19, 22] are not involving the quitting compartment .

Based on the above considerations, we will premeditate two treating methods, namely, prevention of susceptibles from alcoholism and treatment on alcoholism as control variables; hence, we will derive a SATQ-type model. We note the fact that many authors are interested in solving optimal control problems, such as cost minimization and optimal control of various disease, especially with biological background and various mathematical models [2224]. In this paper, we will propose an objective functional which considers not only alcohol quitting effects but also the cost of controlling alcohol. Then, we consider a range of issues related to the optimal control with the method of Pontryagin’s Maximum Principle, including optimal control existence, uniqueness, and characterization.

The organization of this paper is as follows. In the next section, the alcoholism model with prevention for the susceptibles and treatment for alcoholism is formulated. In Section 3, the basic reproduction number and the existence of equilibria are investigated. The stability of the disease free and endemic equilibria is proved in Section 4. Optimal control strategies by the classic method of PMP (Pontryagin’s Maximum Principle) are discussed in Section 5. In Section 6, we give some numerical simulations. We give some discussions and conclusions in the last section.

2. The Model Formulation and Some Fundamental Properties

In this section, we introduce a mathematical model with prevention and treatment for the alcoholism and then study some important properties such as the boundness and positivity of its solutions.

2.1. Model Formulation and Parameter Explanation

The total population is partitioned into four compartments: the susceptible compartment which refers to the persons who never drink or drink moderately without affecting the physical health, the alcoholism compartment which refers to the persons who binge drink and affect the physical health seriously, the treatment compartment which refers to the persons who have been receiving treatments by taking pills or other medical interventions after alcoholism, and the quitting compartment which refers to the persons who recover from alcoholism after treatment and stay off alcohol hereafter. In this paper, we focus on a closed environment, such as a community, a university, or a village. So the total number of population to be considered is a constant; we denote it as . The population flow among those compartments is shown in the following diagram (Figure 1).

The schematic diagram leads to the following system of ordinary differential equations: Here, is the birth number of the population; is the natural death rate of the population; is the fraction of the susceptible individuals who successfully avoid to stay off the alcoholism; is the fraction of the alcoholics who take part in treatments; here, , , and they will be considered as two control variables in Section 5; is the transmission coefficient of the “infection” for the susceptible individuals from the alcoholic individuals; is the rate coefficient of the person who fail to be treated and return to the alcoholism compartment mostly due to their own weak will; is the rate coefficient of the person who have received effective treatment and recovered from alcoholism forever.

2.2. Boundedness of Solutions to System and Positively Invariant Region

It is important to show positivity and boundedness for the system (1) as they represent populations. Firstly, we present the positivity of the solutions. System (1) can be put into the matrix form where and is given by It is easy to check that Due to Lemma   in [25], any solution of (1) is   for all .

We denote ; summing equations in (1) yields so (denoted as ), and the set is a positively invariant region for (1). Therefore, we will consider the global stability of (1) on the set .

3. The Basic Reproduction Number and Existence of Alcoholism Equilibria

3.1. The Basic Reproduction Number

In epidemiology, the basic reproduction number (sometimes called basic reproductive rate or basic reproductive ratio) of an infection is the number of infectious cases that one infectious case generates on average over the course of its infectious period. While in this context, it means the number of persons that an alcoholic will “infect” during his “infectious” period in the pure susceptible environment so that the infected persons will enter the alcoholism compartment. It is easy to see that the model has an alcohol free equilibrium . In the following, the basic reproduction number of system (1) will be obtained by the next generation matrix method formulated in [26].

Let , then system (1) can be written as where The Jacobian matrices of and at the alcohol free equilibrium are, respectively, where The basic reproduction number, denoted by , is thus given by It is easy to see that both of the control parameters contributed to reducing the alcoholism. From this point, the control measures are meaningful.

3.2. Existence of Alcoholism Equilibrium

The endemic equilibrium of system (1) is determined by equations The third equation in (12) leads to From the last equation in (12), we have From the first equation of (12), and together with (13), we can get Substituting (13)–(15) into the second equation of (12) gives By simplifying (16), we can get where

Hence, we get two explicit solutions to (17); one is , which is corresponding to the alcohol free equilibria, and the other is which should be corresponding to the alcoholism equilibria on condition that ; otherwise, the alcoholism equilibria are nonexistent. It is enough to show the positivity of to make sure the existence of alcoholism equilibria on the condition . By some simple calculations, we simplify the expression of to be Since is equivalent to the right side of this inequality is exactly equal to . Hence, we have proved the existence of , so are the alcoholism equilibria. We summarize this result in Theorem 1.

Theorem 1. For system (1), there is always an alcohol free equilibrium . When , besides alcohol free equilibrium , system (1) also has a unique alcoholism equilibrium , where

4. Stability Analysis of Equilibria

For the convenience of subsequent proof, we denote a vector and So the Jacobian matrix of about vector is as the following:

Theorem 2. For system (1), the alcohol free equilibrium is locally asymptotically stable if .

Proof. Since we can easily get that two of the eigenvalues are , while satisfy Thus, Since is equivalent to so and then while Similarly from , we can derive the inequality so It reduces to Hence, , . The proof is complete.

Next, we will turn to investigate the global stability of .

Theorem 3. For system (1), the alcohol free equilibrium is globally asymptotically stable if .

Proof. Consider the subsystem of (1) as follows: Equation (35) can be rewritten as Since and , then for all , we can get According to Lemma 1 in [26], all the eigenvalues of matrix have negative real parts, so the solutions of this subsystem are stable whenever . So as . By the comparison theorem [27], and based on the fact that the total population is constant , it follows that and as . So the alcohol free equilibrium is globally asymptotically stable; the proof is complete.

Theorem 4. For system (1), the alcoholism equilibrium is globally asymptotically stable if .

Proof. Since the total population in model (1) is a constant number , in order to prove the global stability of system (1), it is sufficed to prove the corresponding stability of subsystem (38): We make normalization transform and still use the same symbols to denote the variables; then (38) can be transformed into From (39), we can easily know that the equilibria satisfy the following three equalities to be used later:
Let ; then To eliminate the cross-term and two single-variable terms and , we let By solving them, we can get Next, we let and then Due to so Hence, if and only if , , . According to LaSalle’s invariance principle [28], we can derive the conclusion that the alcoholism equilibria are globally asymptotically stable; the proof is complete.

5. Optimal Control Problem

5.1. The Existence of Optimal Control

In order to investigate an effective campaign to control alcoholism in a community which pursue the goals of the minimized alcoholisms and more recovered individuals, we will reconsider the system (1) and use two control variables to reduce the numbers of alcoholics. The difference is that we will change the parameters into control variable . Their aforementioned definitions allow us to do so. is used to limit the proportion of the susceptible individual to contact with alcoholism, usually by propaganda and education, so that the susceptible individual can stay off alcoholism consciously and be free of “infection,” we can understand the effect of is to prevent the the susceptible from contacting with the alcoholism. The control variable is used to control the alcoholism to take appropriate treatment measures, such as taking pills or seeking other medical help. However, just as a coin has two sides, there will be a lot of costs generated during the control process. So it is advisable to balance between the costs and the alcohol effects. In view of this, our optimal control problem to minimize the objective functional is given by which subjects to system with initial conditions Here, , is the end time to be controlled, is an admissible control set, , and , are weight factors (positive constants) that adjust the intensity of two different control measures.

Next, we will investigate the existence of the optimal control of the above-mentioned problem.

Theorem 5. There exists an optimal control pair such that subjects to the control system (1) with initial conditions (50).

Proof. To prove the existence of an optimal control, according to the classic literature [29], we have to show the following.(1)The control and state variables are nonnegative values.(2)The control set is convex and closed.(3)The right side of the state system is bounded by linear function in the state and control variables.(4)The integrand of the objective functional is concave on .(5)There exist constants and such that the integrand of the objective functional satisfies statements (1), (2) and (3) are obvious satisfied, we only need to test and verify the latter two ones. Since the four state variables have been all proved to be up bounded by , we will get the following equalities: so the fourth condition is set up. As for the last condition, is also true, when we choose , and for all . The proof is complete.

We next come to the core of this section.

5.2. The Characterization of the Optimal Control

With the existence of the optimal control pairs established, we now present the optimality system and use a result from [30]; we can easily know the existence of the solutions to the optimality system (71) which will be gotten later. Firstly, we come to discuss the theorem that relates to the characterization of the optimal control. The optimality system can be used to compute candidates for optimal control pairs. To do this, we begin by defining an augmented Hamiltonian with penalty terms for the control constraints as follows: where are the penalty multipliers satisfying

Theorem 6. Given optimal control pairs and solutions of the corresponding state system (50), there exist adjoint variables , , satisfying with the terminal conditions Furthermore, are represented by

Proof. According to Pontryagin Maximum Principle [2931], we first differentiate the Hamiltonian operator , with respect to states. Then the adjoint system can be written as The terminal condition (56) of adjoint equations is given by , .
To obtain the necessary conditions of optimality (59), we also differentiate the Hamiltonian operator , with respect to and set them equal to zero; then
By solving the optimal control, we obtain
To determine an explicit expression for the optimal control without and , a standard optimality technique is utilized [29]. We consider the following three cases.(i) On the set , we have . Hence, the optimal control is (ii) On the set , we have . Hence, This implies that (iii) On the set , we have . Hence, This implies that
Combining these results, the optimal control is characterized as Using the similar arguments, we can also obtain the other optimal control function The proof is complete.

We point out that the optimality system consists of the state system (50) with the initial conditions , the adjoint (or costate) system (58) with the terminal conditions (59), and the optimality condition (60). Any optimal control pairs must satisfy this optimality system. For the convenience of subsequent numerical simulation in Section 6, we give the optimality system as follows:

5.3. The Uniqueness of Optimal Control

Due to the a priori boundedness of the state, adjoint functions, and the resulting Lipschitz structure of the ODEs, we can obtain the uniqueness of the optimal control.

Lemma 7 (see [23]). The function is Lipschitz continuous in , where are some fixed positive constants.

Theorem 8. For all  , the solution to the optimality system (71) is unique.

Proof. Suppose and are two different solutions of our optimality system (71). Let where is to be chosen.
Accordingly, we have Now we substitute , , , , ,  ,  ,   and ,  ,  ,  ,  ,  ,  ,    into the first ODE of (71), respectively; then we can obtain for and , respectively. Similarly, we can derive for and , respectively; for and , respectively; for and , respectively; for and , respectively; for and , respectively; for and , respectively; for and , respectively.
By Lemma 7, we can obtain
The equations for and the equations for are subtracted, respectively; then we multiply each equation by appropriate difference of functions and integrate from 0 to . Next, we add all eight integral equations and some inequality techniques to obtain uniqueness. The following calculation is similar; for the sake of simplicity, we only take and for an example:
Multiplying both sides of (83) by and integrating from 0 to gives
In the above derivation, we use many scaling techniques for inequality or absolute inequality. Particularly, what should be noted is that to get the first inequality of above derivation, we use the estimation of which has been given before; besides, for the sake of convenience, we note . Furthermore, we notice that the coefficients of all the eight terms in the last formula: , , , , , , , , namely, , , , , , , , are nonnegative and bounded. So there exists a positive constant such that
Combining eight of these inequalities gives Thus, from the above inequality we can conclude that where depends on the coefficients and the bounds depend on . If we choose such that , then , , , , , , , and . Hence, the solution to the optimality system is unique. The proof is complete.

6. Numerical Simulation

6.1. The Simulation of State System (1) without Control Parameters

For the sake of simplicity but without loss of generality, we will perform the numerical simulation of state system (1) with parameters , . Before illustrating the analytic properties of the alcoholism model (1), we will target the populations in the environment of a community or a university, for example, the school of material science and engineering in our university, that is, Lan zhou University of Technology (LUT for short), owing to the accurate and available information we can obtain. Referring to the information provided by the admissions office of LUT, this school will enroll almost 1200 undergraduates and almost 300 various postgraduates at the beginning of fall semester; at the same time, there will be almost 1500 various students graduated and left this school, so the scale of students in school remained almost 6000; we can take the total population . In this simulation, we will take September as the initial time and units in one week, period in one year. According to the investigations of the student union implemented in September every year, we can take initial values as , , , and . It seems that the alcoholism is a little bit more, but it is rather natural because many freshmen feel confused when they are faced with the new environment and a new lifestyle; many of them have no better choice but gather together to drink in small groups to mediate the anxiety and get to know each other; over time, some of them develop the habit of drinking. To a certain extent, for example, the frequent drinking badly affects their study; we can classify them into the alcoholism compartment. Other initial values seem more reasonable, so we need no more explanation. As we know, alcoholism death is seldom happen within one year, so we omit mortality from alcoholism; then how to understand the recruitment rate as well as natural death rate ? We can treat the freshmen admission as the recruitment population and graduation students as the natural “death” parts. So we can take , which is exactly consistent with the value in [16]. As for the infection rate and recovery rate , we will let them be variables, since the drinking behaviors are related to many factors such as the season and the pressure.

According to the data we get from the student union, we choose . To summarize, we list the values of the parameters in Table 1. Using the values of parameters in Table 1, we can plot Figures 2 and 3 which are on condition and , respectively. From Figure 2, we easily know when holds; the solution of system (1) tends to the alcohol free equilibrium and verifies the global stability of . While seen from Figure 3, we also know that if holds, the solution of system (1) tends to the alcoholism equilibrium and verifies the global stability of .

6.2. The Sensitivity of about Two Control Parameters

Although from the expression of the model reproduction number , we can easily find out the fact that the two control variables, that is, and , attribute to reducing the severity of alcoholism; we will still depict the graph between and the two control variables to see more intuitiveness see Figure 4. It seems from the figure that is a monotonically decreasing function about two control parameters, so it is advisable to take two approaches simultaneously to control the alcoholics.

6.3. The Simulation of Optimality System (71)

In this subsection, we will investigate numerically the optimal solution to optimality system (71) by numerical method from [32]; the optimality system is solved with a fourth-order Runge-Kutta scheme. Beginning with a guess for the control variables, the state system is solved forward in time and then those values of state are used to solve the adjoint equations backward in time. The controls are updated at the end of each iteration using the values of optimal controls obtained lastly. The iterations continue until convergence takes place.

In the simulations, we choose the available variable values as Table 1 shows; besides, . The initial value of model (1) is assumed to be , , , and as before.

The ideal weights in objective functional are very difficult to obtain in reality; it needs much work on data mining and fitting. Hence, the acquisition of appropriate practical weights is still a difficult problem and remains for further investigations. The cost associated with and mainly includes the cost of dangerous behaviours during the alcoholism time and educating the public, while the cost associated with mainly comes from health professional and the medical resource including medicines and nursing care. In view of this and taking the expressions of into account, after many numerical simulations, we finally give weighting coefficients as . It should be pointed out that the weights here are of only theoretical interest to reveal the control strategies proposed in this paper. Another point to note is that the maximum control is very difficult to achieve in reality, so we will omit the situation of the maximum control during the series of simulations.

Next, we will make some necessary instructions and explanations to the above simulation graphs. Figures 5, 6, 7, and 8 depict the number of four compartments under different control levels when we choose the weight coefficients in objective function to be . From the four simulation graphs, we can observe the following simple facts, in reducing the total number of alcoholisms and increasing the number of susceptibles; the effectiveness of various control measures is as follows: optimal control is evidently better than middle control, and middle control is better than single control , single control is better than single control , while single control is much better than no control.