Bifurcation and optimal control analysis of a delayed drinking model

Alcoholism is a social phenomenon that affects all social classes and is a chronic disorder that causes the person to drink uncontrollably, which can bring a series of social problems. With this motivation, a delayed drinking model including five subclasses is proposed in this paper. By employing the method of characteristic eigenvalue and taking the temporary immunity delay for alcoholics under treatment as a bifurcation parameter, a threshold value of the time delay for the local stability of drinking-present equilibrium and the existence of Hopf bifurcation are found. Then the length of delay has been estimated to preserve stability using the Nyquist criterion. Moreover, optimal strategies to lower down the number of drinkers are proposed. Numerical simulations are presented to examine the correctness of the obtained results and the effects of some parameters on dynamics of the drinking model.


Introduction
In recent years, with the improvement of our living standards, life styles are becoming more and more diversified, and drinking plays an increasingly important role in people's daily life. However, the current situation of alcoholism around the world is really worrying. Taking China as an example, from the statistical data in the white paper on the status of moderate drinking among Chinese drinkers in 2017 [1], the drinking rate of men and women is 84.1% and 29.3%, respectively. It also points out that 65% of drinkers show unhealthy drinking and the main problem is excessive drinking. Alcoholism has become one of the public health and social problems considered by the world. The harm of alcoholism is very serious. Alcoholism may not only have a serious impact on the brain, nervous system, liver, stomach, and heart, but it also accelerates the spread of AIDS and sexually transmitted diseases. And it can also lead to traffic accidents, social violence, and other serious consequences of criminal acts. It is reported that more than 3 million people died because of alcohol, which accounts for 5.3% of all deaths [2]. Therefore, alcoholism is not only a serious global public health problem, but also a social problem to be solved.
Health-related risky behaviors related to addictions can be regarded as treatable contagious diseases since they can spread through the social interaction [3,4]. It is worth mentioning that health-related risky behaviors such as smoking [5][6][7][8][9][10][11] and drug abuse [12][13][14][15][16] have been explored by mathematical modeling method. In the same way, some drinking models [17][18][19][20][21][22] based on differential equations have been developed to describe the spreading law of drinking dynamics since the pioneering work of [23]. Specially, there are also some models considering the impact of social environment on drinking spread.
Huo and Zhang [24] established a drinking model taking into account the influence of Twitter and showed that Twitter can affect the spread of drinking. Xiang et al. [25] formulated a drinking model considering the effects of public health educational campaigns and demonstrated that the public health educational campaigns can slow down the drinking dynamics. Recently, Xiang et al. [26] proposed the following drinking model considering the effect of immigration on drinking behavior: where P(t), L(t), S(t), and Q(t) denote the number of moderate drinkers, light alcoholics, heavy alcoholics, and alcoholics under treatment at time t, respectively. is the total recruitment rate due to birth; denotes the total number of immigrants; q 1 is the fraction of immigrants entering into light alcoholics, and q 2 is the fraction of immigrants entering into heavy alcoholics; η, α and β are the contact rates; denotes the departure rate of light alcoholics among which a proportion p (0 < p < 1) accept treatment and the rest 1p enter into heavy alcoholics; μ is the natural death rate; d 1 , d 2 , and d 3 are the death rates due to drinking; ε, ρ, and φ are state transition rates.
Evidently, it is instantaneous for alcoholics under treatment to enter into the compartment of heavy alcoholics in system (1). This is not acceptable in real life. Indeed, alcoholics under treatment may enter into the compartment of heavy alcoholics again due to their own reasons or the influence of other external factors. However, for such a process, it is reasonable that there is usually a certain period due to the impact of the accepted treatment and self-control for alcoholics under treatment. In other words, the alcoholics begin to accept treatment at t -τ , then they may become heavy alcoholics again at time t for some reason. From another perspective, in comparison with ordinary differential equations, delayed differential equations present more complex dynamics. In fact, delayed differential equations have been applied in various fields, for example, mathematical biology [27][28][29][30], epidemiology [31][32][33], and computer virus [34][35][36][37]. All the work above revealed that delay can influence dynamics of the systems significantly. Inspired by the above consideration and related research work, we investigate the following delayed drinking model: where τ is the temporary immunity delay for alcoholics under treatment due to the impact of the accepted treatment and self-control.

Numerical simulations of Hopf bifurcation
We choose the same parameters as the ones used in literature [26] Equation (3) turns into 0.3495) loses stability and the specific case of system (2) experiences a Hopf bifurcation, which can be depicted by Fig. 2. In what follows, we are interested in focusing on the effects of some other parameters on dynamics of system (2).
(i) The number of moderate drinkers decreases as the value of increases, whereas the number of light alcoholics, heavy alcoholics, and alcoholics under treatment increases, which can be displayed in Fig. 3. (ii) The number of moderate drinkers and heavy alcoholics decreases; however, the number of light alcoholics and alcoholics under treatment increases when the value of p or φ increases. This phenomenon can be exhibited in Figs. 4 and 5. It is also interesting to note that the system changes from stable state to limit cycle as there is a decrease in p or φ, which can be observed in Fig. 6 and 7. (iii) The number of moderate drinkers and alcoholics under treatment increases due to the increase of ρ; nevertheless, the number of light alcoholics and heavy alcoholics decreases, which can be shown in Fig. 8. It can also be found in Fig. 9 that the system will lose stability and exhibits limit cycle due to the increase of ρ.

Optimal control analysis
We considered an optimal control model for the drinking model. In this model, control is education or media coverage, which changes the behavior of certain individuals in the heavy alcoholics class. This change in behavior leads to subdividing heavy alcoholics into three subclasses, namely S, S 1 , and S 2 . A proportion of the heavy alcoholics populations S decide to change their behavior due to the effect of a successful education campaign and thus enter the S 1 or S 2 class. These two classes S 1 and S 2 have lower transmission rates than the S class and will contribute to lowering the number of moderate drinkers.
We consider optimal control of an ordinary differential equation model, which describes the interaction of education with heavy alcoholics as follows: with the initial conditions where σ ∈ (-τ , 0].
In the above optimal control model, we need three rates α 1 , α 2 , α 3 for S, S 1 , and S 2 respectively for their interactions with the class P. We also state transition rates φ 1 , φ 1 , φ 1 .
Notice that, as a result of interactions of individuals in class S with the control, education u(t), a proportion of heavy alcoholics leave the general heavy alcoholics class S and move to S 1 and S 2 . The rate of moving into class S i for i = 1, 2 is λ i u(t)S(t) .
Now, we will focus on the use of the variable control function u(t), which represents the level of education or media coverage and enables heavy alcoholics to change their behavior.
Our goal is to find the control u(t) and associated state variables P(t), L(t), S(t), S 1 (t), S 2 (t), and Q(t) to minimize the following objective functional: The control set u is where 0 ≤ u min (t) < u max (t) ≤ 1, C is weighted costs associated with the use of the control u(t).
In this subsection, we derive the first order necessary conditions for the existence of optimal control by constructing Hamiltonian H and then applying Pontryagin's maximum principle. To simplify the notations, we write x(t) = [P(t), L(t), S(t), S 1 (t), S 2 (t), Q(t)] T .ρ(t) = [ρ 1 (t), ρ 2 (t), ρ 3 (t), ρ 4 (t), ρ 5 (t), ρ 6 (t)] T . The Hamiltonian is given by Let χ [t 1 ,t 2 ] (t) be the characteristic function defined by Let u * be the optimal control and x * = [P(t), L(t), S(t), S 1 (t), S 2 (t), Q(t)] T be the corresponding optimal trajectory. Then there exists ρ(t) ∈ R 6 such that the first order necessary condition for the existence of optimal control is given by the equation Now we apply the necessary conditions to the Hamiltonian H in (32).
Here we call formulas in Eq. (36) for u * the characterization of the optimal control. The optimal control and the state are found by solving the optimality system, which consists of the state system Eq. (31), the adjoint system Eq. (34), initial conditions at t = 0, boundary conditions Eq. (35), and the characterization of the optimal control Eq. (36). To solve the optimality system, we use the initial and transversality conditions together with the optimal control given by Eq. (36).

Simulations of optimal control strategy
In this section, the drinking model optimal system is solved using Euler's method. The optimal strategy is obtained by solving the state and adjoint systems and the transversality conditions. When viewing the graphs, remember that each of the individuals without control is marked by blue solid lines. The individuals with control are marked by red dash lines. The weight constant value in the objective functional is C = 0.5. Figure 10 represents the behavior time plots of P, L, S, and Q in a system without control and with control. In Fig. 11, we plot the optimal control u(t) in the system. We select the parameters as follows: = 0.8, q 1 = 0.5, q 2 = 0.  The behavior time plots of u * Figure 10 shows that after effective control such as education activities, media coverage, etc., the number of P, L, S, Q is significantly reduced compared with that without control. The results of numerical simulation prove that the effect of media coverage can reduce the number of moderate drinkers, light alcoholics, heavy alcoholics, and alcoholics under treatment in a certain period of time. In particular, the control resulted in a significant reduction in the number of heavy alcoholics.

Concluding remarks
Death rate because of drinking has been higher than that of tuberculosis, HIV/AIDS, and diabetes by WTH in 2018. Mathematical modeling has been used to predict the spread of drinking by scholars. During its spread in a particular population, drinking shows different kinds of delays which essentially impact the dynamics. We proposed a delayed drinking model with the effect of immigration by introducing the temporary immunity delay for alcoholics under treatment due to the impact of the accepted treatment and self-control into the model formulated in [26].
By considering the time delay as a bifurcation parameter, we have seen the effect of time delay on the stability of the model. We observe that if the time delay is suitably small under some certain conditions, then the four populations in system (2) tend to the drinkingpresent equilibrium, which indicates that the spread of drinking can be dominated. On the contrary, it will be out of control when the time delay is large enough. Thus, it gives a conclusion that the time delay is harmful for the control of spread of drinking. The numerical simulation shows that system (2) changes its stable state into limit cycle as p or φ decreases, which makes it difficult to control the spread of drinking. From this point of view, we should strongly appeal to more and more people to join the team of those abstaining from drinking. On the other hand, it is revealed that system (2) becomes stable as we increase the value of ρ, which indicates that alcoholics under treatment should have strong will.
We also find that the number of light alcoholics, heavy alcoholics, and alcoholics under treatment will increase once the amount of immigrants increases, which gives us a suggestion that we should properly control the immigrants in a particular community. Furthermore, since the number of heavy alcoholics decreases when p or φ increases, it is strongly suggested that drinkers should give up drinking and accept treatment as soon as possible, if necessary. However, the number of both light and heavy alcoholics increases with the decrease in ρ. Thus, we can conclude that once the drinkers accept treatment and begin to give up drinking, they should have strong will and self-control, not give up halfway.
Moreover, we study the optimal control of the drinking model. The behavior change of various classes is caused by media coverage and propaganda messages. This work demonstrates the power of public propaganda and is an easy control tool for making relevant propaganda plans and policy decisions.