Optimal Control Strategies in an Alcoholism Model

and Applied Analysis 3


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 [9][10][11][12][13].
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,[14][15][16][17][18][19]. 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 2 Abstract and Applied Analysis 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 crossinfection 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 , 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 [22][23][24]. 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.

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.

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; 1 is the fraction of the susceptible individuals who successfully avoid to stay off the alcoholism; 2 is the fraction of the alcoholics who take part in treatments; here, 0 ≤ ≤ 1, = 1, 2, 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.

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 = ( , , , ) ∈ 4 and ( ) is given by It is easy to check that Due to Lemma 2 in [25], any solution of (1) is ( ) ∈ 4 + for all ≥ 0.

The Basic Reproduction Number and Existence of Alcoholism Equilibria
3.1. The Basic Reproduction Number 0 . 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 0 = ( 0 , 0, 0, 0) = ( , 0, 0, 0). 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 F ( ) = ( The Jacobian matrices of F( ) and V( ) at the alcohol free equilibrium 0 are, respectively, The basic reproduction number, denoted by 0 , 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.

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 Hence, we get two explicit solutions to (17); one is 0 = 0, which is corresponding to the alcohol free equilibria, and the other is which should be corresponding to the alcoholism equilibria on condition that * > 0; 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 0 ≥ 1. By some simple calculations, we simplify the expression of to be Since 0 > 1 is equivalent to the right side of this inequality is exactly equal to ( + 2 )( + ) + . Hence, we have proved the existence of * > 0, so are the alcoholism equilibria. We summarize this result in Theorem 1. (1), there is always an alcohol free equilibrium 0 = ( , 0, 0, 0). When 0 > 1, besides alcohol free equilibrium 0 , system (1) also has a unique alcoholism

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: ) .

Abstract and Applied Analysis 5
Proof. Since we can easily get that two of the eigenvalues are 1 = 2 = − < 0, while 3 , 4 satisfy Thus, Since 0 < 1 is equivalent to so and then Similarly from 0 < 1, we can derive the inequality It reduces to Hence, Re 3 < 0, Re 4 < 0. The proof is complete.
Next, we will turn to investigate the global stability of 0 .
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: Abstract and Applied Analysis 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 Abstract and Applied Analysis Hence, = 0 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.

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 1 , 2 into control variable 1 ( ), 2 ( ). Their aforementioned definitions allow us to do so. 1 ( ) 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 1 ( ) is to prevent the the susceptible from contacting with the alcoholism. The control variable 2 ( ) 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, ( ) ∈ ≜ {( 1 , 2 ) | ( ) is measurable and 0 ≤ ( ) ≤ 1, for all ∈ [0, ]}, is the end time to be controlled, is an admissible control set, , and = 1, 2, 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.
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 .
We next come to the core of this section.

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 ( ) ≥ 0 are the penalty multipliers satisfying = 0 at optimal control * 2 .
To obtain the necessary conditions of optimality (59), we also differentiate the Hamiltonian operator , with respect to = ( 1 , 2 ) and set them equal to zero; then By solving the optimal control, we obtain * To determine an explicit expression for the optimal control without 11 and 12 , a standard optimality technique is utilized [29]. We consider the following three cases.
We point out that the optimality system consists of the state system (50) with the initial conditions (0), (0), (0), (0), 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:

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 * ( ) = min ( , max ( , )) is Lipschitz continuous in , where < are some fixed positive constants.
Abstract and Applied Analysis 11 By Lemma 7, we can obtain * The equations for , , , , , , , V and the equations for , , , , , , , V 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 | * 1 ( ) − * 1 ( )| 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:  (85) 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 , , , , , , , V. If we choose such that + > , then = , = , = , = , = , = , = , and V = V. Hence, the solution to the optimality system is unique. The proof is complete.

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 1 = 0, 2 = 0. 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 = 6000. 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 (0) = 4500, (0) = 1000, (0) = 300, and (0) = 200. 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 = 1500/6000 = 0.25, 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 = 0.4. 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 0 < 1 and 0 ≥ 1, respectively. From Figure 2, we easily know when 0 < 1 holds; the solution of system (1) tends to the alcohol free equilibrium 0 and verifies the global stability of 0 . While seen from Figure 3, we also know that if 0 ≥ 1 holds, the solution of system (1) tends to the alcoholism equilibrium * and verifies the global stability of * .

The Sensitivity of 0 about Two Control Parameters.
Although from the expression of the model reproduction number 0 , we can easily find out the fact that the two control variables, that is, 1 and 2 , attribute to reducing the severity of alcoholism; we will still depict the graph between 0 and the two control variables to see more intuitiveness see Figure 4. It seems from the figure that 0 is a monotonically decreasing function about two control parameters, so it is advisable to take two approaches simultaneously to control the alcoholics. 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 values of optimal controls obtained lastly. The iterations continue until convergence takes place.
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 The total populations to be considered 6000 The rate of populations quitting from alcoholism permanently after treatment Variables The rate of populations failed in treatment and returned to be alcoholic 0.  effort of 1 should be decreased from 0.65 to 0.5 within the first month, and over the next week, it should be increased to the maximum control until 50 weeks, then rapidly decreased to 0 at the end of the simulation. As for the control 2 , it should start from around 0.5 due to the initial alcoholics then increase to 0.55 within one week since the rapid infection and next decrease to almost 0 since the effectiveness of treatment in three weeks, but with the infection going on, the control effort of 2 should gradually increase to the maximum and maintain this level until the tenth week for the purpose of consolidation therapy and preventing rebound; hereafter, it should be gradually decreased to the level of almost 0.43 until the fifty weeks then quickly decreased to 0 in the end. In order to investigate the influence of different weight coefficients in the objective functional on the effect of controlling, at the same time, for a better comparison, we will change the weight coefficients in objective function into 2 = 10 4 ; 1 = 10 2 , and we will list the corresponding numerical simulation results as Figures 11-16 show.
When we change the weight coefficients in objective function into 1 = 10 4 , 2 = 10 2 , we find that the results of simulations derived from the graphs are very similar to the ones before. We speculate that the most likely reasons of this result are due to three respects; one is that the weight coefficients are not too sensitive in the numerical simulation, and another possible reason is that both of the two controls are important, in some sense, equivalently important. The last but not the most unlikely reason is that we have not found the most appropriate weight coefficients in the simulation, which is very difficult to find as previously mentioned.

Conclusions
In this paper, we formulate an alcoholics quitting model and firstly investigate the variation discipline of various populations from the perspective of global stability; then we propose an objective functional to examine two different control measures (i.e., prevention and treatment) on the effect of alcohol. The basic reproduction number of the model was derived and the global stability of the two equilibria is given.  related numerical simulation, we can easily see that the two control strategies are effective in the alcoholics process.
Using Pontryagin's Maximum Principle, we firstly determine the necessary conditions for existence of optimal control pairs. The uniqueness of the solution to the optimality system (71) is derived by the classical method of contradiction. Numerical simulations of the model suggest that the two different groups of weights in the objective function have much similar effects on the transmission of the alcoholism; from this point, the two control measures are almost equally important in controlling the alcoholism, although they will probably have great influences on the cost of the objective function. From the simulation figures, it seems that the effect of optimal control, which is measured by the reduction in the number of alcoholics and the increase in the number of susceptibles, is much better than other control strategies as noted earlier in the simulation section. According to the realtime curve of two optimal controls, we point out the specific implementation methods of optimal control which can be achieved in practice.