Mathematical Modeling of the Spread of Alcoholism Among Colombian College Students

In this paper, we present a nonlinear mathematical model, describing the spread of high-risk alcohol consumption behavior among college students in Colombia. We proved the existence and stability of the alcohol-free and drinking state equilibrium by means of Lyapunov function and LaSalle’s invariance principle. Also, we apply optimal control to study the impact of a preventive measure on the spread of drinking behavior among college students. Finally, we use numerical simulations and available data provided by the United Nations Office on Drugs and Crime (UNODC) and the Colombian Ministry of Justice to validate the obtained mathematical model.


Introduction
Drinking alcohol is widely socially accepted and associated with relaxation and pleasure, and some people drink alcohol without experiencing harmful effects. However, according the World Health Organization (WHO) in 2016, the harmful use of alcohol resulted in some 3 million deaths (5.3% of all deaths) worldwide and 132.6 million disability-adjusted life years (DALYs)-i.e. 5.1% of all DALYs in that year. Mortality resulting from alcohol consumption is higher that caused by diseases such as tuberculosis, HIV/AIDS and diabetes. Among men in 2016, an estimated 2.3 million deaths and 106.5 million DALYs were attribute to the consumption of alcohol. Women experienced 0.7 million deaths and 26.1 million DALYs attributable to alcohol consumption.
The report, which comes out every four years, reveals that all deaths attributable to alcohol consumption worldwide, 28.7% were due to injuries, 21.3% due to digestive diseases, 19% due to cardiovascular diseases, 12.9% due to infectious diseases and 12.6% due to cancers. About 49% of alcohol-attributable DALYs are due to noncommunicable and mental health conditions, and about 40% are due to injuries. Harmful use of alcohol caused some 1.7 million deaths from noncommunicable diseases in 2016, including some 1.2 million deaths from digestive and cardiovascular diseases (0.6 million for each condition) and 0.4 million deaths from cancers [1]. Globally an estimated 0.9 million injury deaths were attributable to alcohol, including around 370000 deaths due to road injuries 150000 due to self-harm and around 90000 due to interpersonal violence. Of the road traffic injuries, 187000 alcohol-attributable deaths were among people other than drivers.
Worldwide, alcohol was responsible for 7.2% of all premature (among persons 69 years of age and younger) mortality in 2016. People of younger ages were disproportionately affected by alcohol compared to older persons, and 13.5% of all deaths among those who are 20 − 39 years of age are attributed to alcohol.
Young people often are in situations where drinking is tolerated or even reinforced [2], [3]. Specifically, problem drinking during the college years is a significant public health concern with a variety of negative consequences like: missing class, falling behind in school work, engaging in unplanned and or unprotected sexual activity, getting hurt or injured, damaging property, getting in trouble with law enforcement, or requiring medical treatment for an alcohol overdose [4], [5]. Heavy drinking also increases probability of sexual victimization; alcohol use increased the risk of being victim to sexual assault as well as the severity of the sexual assault in a sample of college women [6], [7]. There are several variables linked to college drinking in the literature, these influences include: demographic variables, personality, drinking history, alcohol expectancies, drinking motives, stress and coping and peer and family influence [4].
The situation in Colombia is no better than the rest of the world. It has been estimated that Colombian youth tend to start drinking alcohol at age 12 [8], [9]. The alcohol is the most common psychoactive substance used in the country. According to the Colombian Ministry of Health and Social Protection around seven million of people between the ages of 12 and 65 years are current drinkers (35 percent of this age group), around 2.4 million reported binge drinking (15 percent of this age group). The most prevalence of drinking is reported between ages 18-24 (46%), followed ing.cienc., vol. 16, no. 32, pp. 195-223, julio-diciembre. 2020. by the ages 25-34 (43%). Among Colombian college students the data provided by the United Nations Office on Drugs and Crime (UNODC) and the Colombian Ministry of Justice [8] show that the percentage of drinkers has decreased from 87.76% in 2012 to 84.28% in 2016 among men, and from 82.51% in 2012 to 79.59% in 2016 among women.
To respond to problems related to drug use, Colombia is implementing a national Policy on Drug Abuse 2014-2021 that is based on four pillars: i)prevention; ii) treatment; iii) risk and harm reduction; and iv) health promotion. The main programme on drug prevention is Strengthening Families: Love and Limits, initiated in 2012 under the leadership of the Ministries of Justice and Law, and Health and Social Protection, and in close cooperation with the World Health Organization and UNODC. Among the new partners to the programme, there are governorships, mayorships, and family welfare institutions reaching 14.000 families, and impacting over 50.000 persons, in 100 municipalities of 24 provinces in the country.
In the present work, our goal is to analyze and explain the dynamic of the prevalence of alcohol consumption among Colombian college students through the analysis of a mathematical model given by a system of nonlinear differential equations.
Mathematical modeling of infectious disease has its origin in the early 20th with an important work of Kermack and McKendrick [10] where it is describes the relationship between susceptible, infected and immune individuals in a population. The Kermack-McKendrick epidemic model was successful in predicting the behavior of outbreaks very similar to that observed in many recorded epidemics [11]. Mathematical models provide useful tools to study the mechanism by which diseases spread and to predict the future course of an outbreak and to evaluate strategies to control an epidemic [12], [13].
In this paper, we consider drinking behavior as a "disease" which can be spread through peer-influence. We use a system of nonlinear differential equations to represent the relationships between non-drinkers, occasional drinkers and alcohol dependent individuals. Our model is a variant of the classical susceptible, exposed and infectious (SEI) model, here we consider two levels of infectives, namely occasional drinkers and alcohol dependent.
Several different mathematical models for drinking had been formulated and studied [14], [15], [16], [17], [18]. In this paper, we use a modified model of [19] to model alcoholism as epidemic and to the best of our knowledge such model has not been used to analyze the problem of alcohol consumption among Colombian college students. One main difference between our model and the one considered in [19] is that a person with alcohol dependence cannot recover unless they pass through a disease control strategy like a treatment. Also, in our model is not possible for a susceptible individual to become an alcohol dependent without passing through the moderate alcohol consumption class. Furthermore, we analyze separately the prevalence of alcohol consumption in both women and men.
In order to study the stability of the model proposed we found the basic reproduction number R 0 [20], [21], [22]. Then the stability analysis of the model is made using the basic reproduction number. Next we apply an optimal control strategy to reduce the number of infected individuals. Finally, numerical simulations confirm our analytic results.
The paper is organized as follows: In section 2, we present the model description and the stability analysis of drinking-free and endemic equilibria. Section 3 presents an optimal control problem, which we use to study the impact of a treatment on the spread of drinking behavior. After the mathematical modelling, we use the obtained model to simulate it in Section 4 with the parameters estimated from recent statistical data based on the UNODC and the Colombian Ministry of Justice report of the 2017 [8]. We end with Section 5 of conclusions and discussions of our results.

Peer influence
According to B. Bradford Brown, adolescents are influenced by the world through four different means: family, school, the media, and their peers. Adolescence initiates a process in which individuals search for their own personality and a sense of identity. In this period of life, the most important role in the development of values and attitudes is played by peers [23]. Brown points that as adolescents grow older, friends become more influential with regard to substance use. When children start to drink at an (nonnormative) early age, alcohol use is mainly affected by familiar factors such as parental alcohol use and parenting practices, and ing.cienc., vol. 16, no. 32, pp. 195-223, julio-diciembre. 2020. by genetic susceptibility. However, when adolescent's drinking goes beyond occasionally trying out alcohol, the role of friends becomes stronger. This is partly due to the increasing importance of friendships and intimate relationships during the course of adolescence and the fact that when teenagers go out to parties and public drinking places, their drinking is concentrated in these peer settings [24]. Therefore, adolescents' alcohol consumption is closely associated with the drinking behavior of their peers. Thus, the peer influence exerts a role in explaining the willingness to drink alcohol [25], [26]. For instance, students who associate with more friends who drink tended to consume more alcohol than those students who fewer friends who drank according to the student's self-report [27], [28].
Youth in Colombia as in other parts of the world congregate in more or less recognized peer groups and engage in collective behavior. They provide referential parameters and behavioral codes [23]. The university environment promotes the creation of natural environments, around different activities, such as: academic, social and recreational activities. For example, in the case of males, football (soccer) represents an element of identification among Colombian adolescents and young adults living in the cities. After or during the game, depending on whether one is a player or only a spectator, the consumption of alcohol constitutes the main activity for many youths and one of the reasons to be part of the group [23]. That is, how most of the Colombian colleges students reported consuming alcohol with their group of friends (76%), followed by the family (24.9%), partner (15.3%), and coworkers (6.6%) [29].
Our model assumes the effects of peer pressure in the dynamics of alcohol abuse disorders. As presented, it only considers peer influence in the recruitment of new alcohol users (the nonlinear terms).

Model formulation and equilibrium discussion
We propose a nonlinear system of three ordinary differential equations, which is analyzed to find the equilibrium points and their stability, including the threshold parameter R 0 known as the basic reproduction number. The basic reproduction number, is defined as the expected number of secondary infections produced by an index case in a completely susceptible population. This number is a measure of the potential for disease spread within a population. Usually R 0 is defined as the dominant eigenvalue of a matrix K, called the next-generation matrix (NGM) that relates the numbers of newly infected individuals in the various categories in consecutive generations [20], [21].
We consider three classes of individuals in a constant population size of N . The class S(t) for those individuals that are susceptible to becoming regular alcohol users (people who did not drink any alcohol in the previous 12 month period, this includes former drinkers and lifetime abstainers), E(t) who consumed alcohol in the previous 12 month period, but have not become alcohol dependent and D(t) alcohol consumers now dependent on alcohol (to identify the individuals in this class we use the Alcohol Use Disorders Identification Test, which is a simple and effective method of screening for unhealthy alcohol use). We describe the dynamics of drinking behavior by the following three nonlinear differential equations: where µ is the combined birth/death rate, β 1 is the effective influence rate between E and S, β 2 is the effective influence rate between D and S, α is the rate at which occasional drinkers become alcohol dependent.  Figure 1: Dynamics of the model. The boxes represent the subpopulation and the arrows the transition between the subpopulations.
All the parameters mentioned above are positive constants. Notice thaṫ N =Ṡ +Ė +Ḋ = 0 thus N (t) is constant for all t. Thus, the solutions of (1) with positive initial conditions are bounded in Ω 1 given by It's straightforward verify that Ω 1 is positively invariant with respect to (1). (1) with positive initial conditions is positive in Ω 1 .
Proof. Considering that the right-hand of (1) is continuous with continuous partial derivatives, then (1) has a unique solution in Ω 1 . And sincė we have This completes the proof.
Since the sum of the right hand sides of (1) adds to zero, so S(t)+E(t)+ D(t) is constant over time and equal to the size of the total population N . Consequently, the system (1) can be reduced to the following: Now, we show some numerical simulations of the dE dt = 0 and dD dt = 0 isoclines for different values of the parameters α, β 1 , β 2 , µ. Here the vector field is defined by the system (2) an horizontal and vertical, respectively, in the E − D plane. The flow toward en equilibrium is indicated by arrows in the regions between these isoclines. Notice that fromḊ = αE − µD, theḊ = 0 isocline is a line through the origin with the slope α µ . ing.cienc., vol. 16, no. 32, pp. 195-223, julio-diciembre. 2020.

Equilibrium points and stability of the mathematical model
Setting the right-hand side of equations (1) equal to zero we can find the equilibrium points. The system (1) has a unique disease-free equilibrium, P 0 (N, 0, 0). To analyze the stability of P 0 we compute the basic reproductive number R 0 for the model (1). We apply the next generation technique as proposed in [21,22]. The infectious classes in this model are E(t) and D(t). Thus, at P 0 we obtain F = β 1 β 2 0 0 and V = α + µ 0 α −µ where the matrix F is related to the rate of increase of secondary infections and V to the rate of the disease progression, death and recovery. The next generation matrix, K = F V −1 , is non-negative and therefore it has a non-negative eigenvalue, R 0 = ρ(F V −1 ) and non-negative eigenventor ω associated with R 0 . Thus, applying the next generation matrix technique we obtain the basic reproduction number Let the right side of each of the two differential equations equal to zero in system (2), giving the equations: and From (4) we obtain D = αE µ and substituting in (5) we have Solving equation (6) for E we obtain the solutions E * = 0 or Thus, the system (1) has an unique positive endemic equilibrium provided R 0 > 1.

Sensitivity analysis of R 0
Now, we examine the sensitivity of R 0 to each of its parameters, following [30], the normalized forward sensitivity index with respect to each of the parameters is calculated. Indeed: Thus, we conclude that the basic reproduction number R 0 is most sensitive to changes in β 1 and β 2 . Therefore, an increase (or decrease) in the value of β 1 or β 2 leads to a corresponding increase (or decrease) in R 0 . Additionally, the parameter α may have directly or inversely proportional relationship with the reproduction number depending of sign(β 2 − β 1 ). Thus, given R 0 's sensitivity to β 1 and β 2 , it seems sensible to focus efforts on the reduction of β 1 and β 2 .

Parameter estimation
We consider the population of college students from 18 years of age to 28 years of age. According to [9] among Colombian college students the age range of highest alcohol consumption is 21 − 22 years. Now, in the simplest case Then, Thus, 1 − E(3) E(0) is the proportion of Colombian college students that leave E to become alcohol dependent in three year. It is estimated [9] that 12.3% of men and 6.6% of women leave E. Also, 11% of men and 6.2% of women leave S to become occasional drinkers. Therefore, 1 − E(3) E(0) = 0.123, 1 − e −3α = 0.123, thus α 0.04 for men, and 0.02 for women. Then, the average time as occasional drinkers E for men is When R 0 > 1 the system has a unique equilibrium, Thus, near the equilibrium dS dt = −βS, where β = β 1 E * +β 2 D * N . Then, And, using that 11% of men leave S, we get 1 − e −3β = 0.11, so β 0.03. Now, notice that Therefore, assuming the mortality rate µ = 0.1, we can obtain (R 0 − 1)µ = 0.03, then R 0 = 1.3 for men, and similarly, we can get R 0 = 1.2 for women.
Theorem 2.1. If R 0 ≤ 1, the disease-free equilibrium P 0 is globally asymptotically stable in Ω 1 . If R 0 > 1, the disease-free equilibrium P 0 is an unstable saddle point.

Proof. Constructing a suitable Lyapunov function
Its derivative along the solution to the system (1) iṡ Clearly the last expression is negative provided R 0 ≤ 1. Furthermore, L = 0, if and only if E = 0. Therefore, the maximum invariant set in {(S, E, D) ∈ Ω 1 :L = 0} is the singleton {P 0 }. LaSalle's invariance principle implies that P 0 is globally asymptotically stable in Ω 1 .
When R 0 > 1, to analyze the stability of P 0 , we use the method of first approximation. Indeed, the Jacobian matrix of the system (1) at a point P 0 (N, 0, 0) is Its characteristic equation is det(J(P 0 ) − λI) = 0, where I is the unit matrix. Thus, the characteristic equation becomes to Here (7) has one positive root and two negative roots provided αµ + µ 2 − (αβ 2 + β 1 µ) < 0 which is equivalent to R 0 > 1. Consequently P 0 is an unstable saddle point.
Theorem 2.2. When R 0 > 1, the unique positive endemic equilibrium P * is locally asymptotically stable in Ω o 1 .

Optimal control
Recently, different control strategies have been used to prevent the spread of different diseases [31], [32], [33], [34], [35], [36], [37]. Now, our objective is to better understand and to predict the effect of a preventive measure on the alcohol dependent individuals over time. Thus, we introduce some preventive control strategy to the system (1) using a parameter T ∈ [0, 1], where T = 0 indicates that no control strategy is applied, while T = 1 means that preventive control strategies are 100 percent effective. Preventive control strategies are any measures implemented to reduce the spread of the disease, in our particular case of the spread of risky alcohol consumption, some of them could be therapy through a rehabilitation program. Thus, we consider the following optimal control problem: subject to given initial conditions S(0), E(0), D(0) > 0, where T denotes a preventive measure (T is the control variable). Note that in the absence of control measures, that is, when T (t) ≡ 0, system (8) reduces to (1). We consider the set of admissible control functions where t f denote the duration of the treatment program, and T max ≤ 1.

For the system (8) the basic reproduction number with control is given by
Based on the stability results of the previous section it can be deduced that the disease disappears when R 0 (T ) ≤ 1 (see figure 11), that is, when Note that the last expression holds only if R 0 > 1 and α + µ − β 1 > 0; that is, in the presence of the disease. Which is reasonable, since it is only in the presence of the disease that a control strategy should be implemented. Now, our goal is to minimize the number of alcohol dependent individuals and the cost required to control the disease by treating the infected.
Thus, the optimal control problem consists in subject toṠ with 0 < k 1 , k 2 < ∞. To solve the problem, we apply the Pontryagin maximun principle [38]: the Hamiltonian is given by From the conditions ∂H ∂T = 0 and T ∈ Ω 2 the extremal control is given by while the adjoint system asserts that the co-state variables p i (t), 1 ≤ i ≤ 3, satisfẏ and given that the minimization problem stops at t f , we have a boundary condition (called transversality condition) for the terminal values of the co-state variables (see [35] for more details about the necessary conditions for finite-horizon optimal control), that is: Therefore, in order to solve the optimal control problem, we solve the following boundary value problem: ing.cienc., vol. 16, no. 32, pp. 195-223, julio-diciembre. 2020.
with given initial conditions and transversality conditions (16), where T (t) is given by (14). We solve numerically the system (17) in Section 4.

Numerical simulation
In this section, we use numerical simulations to show the dynamical behavior of our results. We solve the system (1) using an adaptive Runge-Kutta-Fehlberg method of order four. The optimal control problem is solved numerically by using a direct multiple shooting method and MATLAB bvp4c routine.
In order to establish the initial conditions for the system (1) we present some available data provided by the United Nations Office on Drugs and Crime (UNODC) and the Colombian Ministry of Justice [8] about alcohol consumption among college students in Colombia.  In Table 2 and Table 3 we can observe the statistics of alcohol consumption in the university population during the years 2009, 2012 and 2016. Table 2 shows the percentage of occasional drinkers who consume alcohol but have not become alcohol dependent [8]. In Table 3 we can see the number of students at high risk of alcohol dependence. For instance, according to the criteria AUDIT(Alcohol Use Disorders Identification Test) 13.9% of men presented alcohol dependence compared with 5.5% of women during the year 2009 [8]. Next, we perform numerical simulations of the college students scenario, for both female and male population linked to the Colombian higher education system. We use the model (1) to study the dynamics of the different subpopulations S(t), E(t) and D(t). We simulate several years in order to observe the dynamics of these populations. Table 4 shows numerical simulation for the model (1), there we take the information obtained from alcohol consumption in 2009.    Table 5 shows the variation of β 2 and its effects on the population with alcohol dependence (see Figure 7). We observe that for R 0 = 4.25 the percentage of individuals with alcohol dependence is 57.3% and for R 0 = 1.85 is 34.5% for both men and women.  The solution for the optimal control problem (17) is illustrated in Figures 8 and 9. Figure 8 shows the evolution of the number of susceptible and occasional drinkers individuals. Figure 9 shows the evolution of the number of alcohol dependent and the variation of co-state variables.
In what follows, we assume that t f = 5 because the World Health Organization goals for most diseases are usually fixed for 5 years periods.   Figure 10 illustrates the importance of the control strategy in the reduction of the number of alcohol dependent individuals. However, graphs still show a disease prevalence over time.  Figure 11 shows that if T max is greater than T c (see equation (10)) there is a decrease in the number of alcohol dependent population while if T max is lower than T c there is an increase of alcohol dependent population.

Conclutions
We present a simple mathematical model of spread of the alcohol abuse among college students in Colombia, with peer pressure as the main element in the recruitment mechanism of new drinkers. We introduced three population classes, namely non-drinkers, occasional drinkers and alcohol dependent. Also, we divided our target population into men and women enrolled in higher education in Colombia.
Here we have found R 0 = αβ 2 + β 1 µ µ(α + µ) as the basic reproduction number of the model system (1). We prove that dynamics of system (1) are completely determined by the threshold parameter R 0 . Sensitivity analysis of R 0 identifies β 1 and β 2 as the most important parameters for the reduction of R 0 , the third most important parameter for the threshold R 0 is α, which is related to the transition from occasional drinkers to alcohol dependent. Notice, that α could have directly or inversely proportional relationship with R 0 depending of sign(β 2 − β 1 ).
Consequently, sensitivity analysis tells us that, some possible strategies for reducing drinking problems on college students are to minimize the ability of occasional drinkers and alcohol dependent to directly recruit non-drinkers, and to increase alcohol prevention programs for youth.
Our model has two steady states, alcohol-free state P 0 and drinking state P * . Using Lyapunov function theory we show that if R 0 ≤ 1, then P 0 is globally asymptotically stable, and if R 0 > 1, P 0 is an unstable saddle point. Also, we prove that is R 0 > 1, the unique endemic equilibrium P * is locally asymptotically stable.
An important conclusion is about the prevalence of alcohol abuse among higher education students in Colombia. An analysis based on parameter estimation shows that the basic reproduction number R 0 is greater than one, it indicates that, there is a greater risk for the spread of alcohol dependence among Colombian college students.
The numerical simulation shows that an optimal control strategy indeed helps to reduce the number of alcohol dependent individuals. Also, we can observe numerically that when the applied control has an effectiveness T > T c the disease is controlled. Finally, all our important mathematical findings are numerically verified using MATLAB. We have numerically verified that, P 0 is stable when R 0 ≤ 1 and when R 0 > 1, endemic equilibrium P * becomes stable. Future research directions include extending this study to incorporate more complex transitions, and more classes such as treatment populations, among others.