Abstract
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 [9–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–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 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 [22–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.
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 [29–31], 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.
Figures 9 and 10 depict the optimal control law of , respectively. In the beginning of the simulation, the control effort of 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 , 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 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 ; , and we will list the corresponding numerical simulation results as Figures 11–16 show.
When we change the weight coefficients in objective function into , , 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.
7. 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. From the expression of the basic reproduction number and 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 real-time curve of two optimal controls, we point out the specific implementation methods of optimal control which can be achieved in practice.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was partially supported by the NSF of Gansu Province of China (2013GS09485), the NNSF of China (10961018), the NSF of Gansu Province of China (1107RJZA088), the NSF for Distinguished Young Scholars of Gansu Province of China (1111RJDA003), the Special Fund for the Basic Requirements in the Research of University of Gansu Province of China and the Development Program for Hong Liu Distinguished Young Scholars in Lanzhou University of Technology.