Epidemic models have been widely used to study contagion processes such as the spread of infectious diseases [1] and rumors [2] . This kind of model has also been used for the spread of social habits, such as the smoking habit [3] , cocaine [4] and alcohol consumption [5] , obesity [6] , corruption [7] , cooperation [8] , ideological conflicts [9] , and also to other problems like rise/fall of ancient empires [10] , dynamics of tax evasion [11] and radicalization phenomena [12] .
The main reason such social behaviors can be modelled by contagion processes is the response by elements of the ensemble to the social context of the studied subject. Both social or peer pressure and positive reinforcement from other agents, regardless if the behavior brings positive or negative consequences to the individual, can influence each one's way of life. Therefore, models for the epidemics of infectious diseases are also able to describe the spread of such tendencies, like alcoholism [13, 14] .
The standard medical way of categorizing alcohol consumption [15] is in three groups -nonconsumers, moderate (or social) consumers and risk (or excessive) consumers; thus, modeling of the interactions and consequent changes of an individual from one group to another is governed by interaction parameters. One interesting aspect that should be taken into consideration when modeling alcohol consumption is the tendency on some individuals to gradually increase their consumption rate, not due to social susceptibility, but when under stressful or depressing circumstances, since alcohol plays a major role both as cause and consequence of depression, for instance [16] . This means that one can attach a probability of a moderate drinker to become an excessive drinker that is dependent only on the actual moderate drinkers population size, instead of the two population groups involved in the change. If one considers the current world situation with the recent coronavirus disease 2019 outbreak (COVID-19), this self-induced increase in alcohol consumption is not only realistic, but also becomes more prominent -this has been observed in a myriad of studies this year detailing the consequences and dangers of both alcohol withdrawal (in places where it has become harder to legally acquire alcohol during the pandemic) and alcohol consumption increase [17, 18, 19, 20] .
This work is organized as follows. In Section 2, we present the model and define the microscopic rules that lead to its dynamics. The analytical and numerical results are presented in Section 3, including comparisons with Brazil's alcohol consumption data for a range of eleven years, used as a case study in order to evaluate the present model. Finally, our conclusions are presented in section 4.
Our model is based on the proposal of references [5, 13, 14, 21, 22, 23, 24, 25] that treat alcohol consumption as a disease that spreads by social interactions. In such case, we consider an epidemic-like model where the transitions among the compartments are governed by probabilities. In this work we consider population homogeneous mixing, i.e., a fully-connected population of N individuals. This population is divided in 3 compartments, namely:
• S: nonconsumer individuals, individuals that have never consumed alcohol or have consumed in the past and quit. In this case, we will call them Susceptible individuals, i.e., susceptible to become drinkers, either again or for the first time;
• M: nonrisk consumers, individuals with regular low consumption. We will call them Moderated drinkers;
• R: risk consumers, individuals with regular high consumption. We will call them Risk drinkers;
To be precise, a moderated drinker is a man who consumes less than 50 cc of alcohol every day or a woman who consumes less than 30 cc of alcohol every day. On the other hand, a risk drinker is a man who consumes more than 50 cc of alcohol every day or a woman who consumes more than 30 cc of alcohol every day [5] 1 . Since we are considering a contagion model, the probabilities related to changes in agents' compartments represent the possible contagions. The transitions among compartments are as following: In the above rules, β represents an "infection" probability, i.e., the probability that a consumer (M or R) individual turns a nonconsumer one into drinker. The Risk drinkers R can also "infect" the Moderated M agents and turn them into Risk drinkers R, which occurs with probability δ. These two infections occur by contagion, in our model, where individuals belonging to a group with a higher degree of consumption can influence others to drink more via social contact. This transition M → R can also occur spontaneously, with probability α, if a given agent increase his/her alcohol consumption -this is the only migration pathway from one group to another, in this model, that does not depend on the population of the receiving compartment, since it corresponds to a self-induced progression from Moderate (M) to Risk (R) drinking. As stated in the introduction, above, the increase of alcohol consumption has been documented to occur under stressful circumstances (like the COVID-19 pandemic) or clinical depression, regardless of social interaction with Risk drinkers. Finally, the probability γ represents the infection probability that turn Risk drinkers R into Susceptible agents S. In this case, it can represent the pressure of social contacts (family, friends, etc) over individuals that drink excessively. We did not take into account transitions from Risk (R) to Moderate (M), assuming that, as a rule, once an individual reaches a behavior of excessive consumption of alcohol, contact with Moderate drinkers does not imply on a tendency to lower one's consumption -meanwhile, it is assumed that contacts that do not drink at all are able to exert a higher pressure on them to quit drinking. It is not that the from Risk to Moderate transition cannot occur -it is just that for our model this probability, when comparing it with the overall picture, is negligible.
For simplicity, we consider a fixed population, i.e., at each time step t we have the normalization condition s(t) + m(t) + r(t) = 1, where we defined the population densities s(t) = S(t)/N, m(t) = M(t)/N and r(t) = R(t)/N. Since we will only deal with the relative proportions among the three different groups in relation to the total population N, i.e. the population densities, we will not take into account birth-mortality relations and populational increase/decrease effects. So, even if N is not a constant number, for all mod-elling purposes it will not matter due to the fact that we will deal only with the s(t), m(t) and r(t) subpopulations in relation to the total population. One other way of looking at this approximation is to consider only the adult population as relevant to our modelling, and assume that new individuals coming of age correspond to the number of deaths [24, 25] .
Based on the microscopic rules defined in the previous subsection, one can write the master equations that describe the time evolution of the densities s(t), m(t) and r(t) as follows,
and we also have the normalization condition
valid at each time step t.
First of all, one can analyze the early evolution of the population, for small times. Considering the initial conditions s(0) ≈ 1, m(0) ≈ 1/N and r(0) = 0, one can linearize Eq. (2) to obtain
that can be directly integrated to obtain m(t) = m 0 e α(R 0 −1) t , where m 0 = m(t = 0) and one can obtain the expression for the basic reproduction number
As it is usual in epidemic models [1, 27] , the disease (alcoholism) will persist in the population if R 0 > 1, i.e., for β > α. One can start analyzing the time evolution of the three classes of individuals. We numerically integrated Eqs. (1), (2) and (3) the effects of the variation of the model's parameters. As initial conditions, we considered s(0) = 0.99, m(0) = 0.01 and r(0) = 0, and for simplicity we fixed α = 0.03 and δ = 0.07, varying the parameters β and γ. In Fig. 1 (a) , (b) and (c) we exhibit results for fixed β = 0.07 and typical values of γ. One can see that the increase of γ causes the increase of s and the decrease of m and r. Remembering that γ models the persuasion of nonconsumers S in the social interactions with risk drinkers R, i.e., the social pressure of individuals that do not consume alcohol over their contacts (friends, relatives, etc) that consume too much alcohol. On the other hand, in Fig. 1 (d) we considered γ = 0.15 and β = 0.07. For this case, where we have β > γ, we see that the densities evolve in time, and in the steady states we observe the survival of only the risk drinkers, i.e., for t → ∞ we have s = m = 0 and r = 1. This last result will be discussed in more details analytically in the following. As we observed in Fig. 1 , the densities s(t), m(t) and r(t) evolve in time, and after some time they stabilize. In such steady states, the time derivatives of Eqs. (1) -(3) are zero. In the t → ∞ limit, Eq. (1) gives us (−β m − β r + γ r) s = 0, where we denoted the stationary values as s = s(t → ∞), m = m(t → ∞) and r = r(t → ∞). This last equation has two solutions, one of them is s = 0 and from the other solution we can obtain a relation between r and m,
Considering now the limit t → ∞ in Eq. (2), one obtains β s r = (α + δ r − β s) m .
If the obtained solution s = 0 is valid, this relation gives us m = 0 and consequently from (4) we have r = 1. This solution represents an absorbing state [28, 29] , since the dynamics becomes frozen due to the absence of S and M agents. We will discuss this solution in more details in the following. Considering now the relation (7) and the normalization condition (4), one can obtain
Substituting (9) and (7) in (8) 
Considering this result (10) in Eqs. (9) and (7) we obtain, respectively
The obtained eqs. (10) -(12) represent a second possible steady state solution of the model, that is a realistic solution since the three fractions s, m and r coexist in the population.
We can look to Eq. (10) in more details. It can be rewritten in the critical phenomena perspective as [30, 31] 
where the critical points are given by The competition among the contagions cause the occurrence of such three regions in the model. From one side we have drinkers (moderated and risk) influencing nonconsumers to consume alcohol, with probability β. On the other hand, we have the social pressure of nonconsumers over risk drinkers, with probability γ, in order to make such alcoholics to begin treatment and stop drinking. Finally, it is important to mention the parameter α, that drives the only transition of the model that does not depend on a direct social interaction. That parameter models the spontaneous increase of alcohol consumption, and it is also responsible for the first phase transition (together with γ), since we have β (1) c = 0 for α = 0. It means that the alcohol consumption (the "disease") cannot be eliminated of the population after a long time if there is a spontaneous increase of alcohol consumption from individuals that drink moderately, which is a realistic feature of the model.
For clarity, we exhibit in Fig. 3 the phase diagram of the model in the plane β versus γ, separating the three above discussed regions. In Fig. 3 , the absorbing phase with s = 1 and m = r = 0 is located in region I for β < β (1) c , the coexistence phase is denoted by II for β
(where the three densities coexist) and the other absorbing phase where s = m = 0 and r = 1 is located in region III. From this figure we see the mentioned competition among the contagions. Indeed, if β is sufficiently high, many nonconsumers become moderated drinkers. Such moderated drinkers will become risk drinkers (via probabilities α and δ), and in the case of small γ we will observe after a long time the disappearance of nonconsumers and moderated drinkers (region III). In the opposite case, i.e., for high γ and small β, the flux into the compartment S is intense, and in the long-time limit the other two subpopulations M and R disappears (region I). Finally, for intermediate values of β and γ the competition among the social interactions lead to the coexistence of the three subpopulations in the stationary states (region II). It is worthwhile to mention that the sizes of regions I and II are directly dependent on probability α, while region III is always fixed due to Eq. (18) . This means that, if parameter α is increased, region I will become gradually larger, which is an indication that the spontaneous evolution from moderate to risk drinking behavior increases the latter's absorbing state. In consequence, since probability α represents a percentage of moderate drinkers that become risk drinkers without the need for social interaction, it is a crucial factor not only to implement the theoretical model but also to identify a possible percentage of the population that has a natural tendency to present excessive alcohol consumption behavior, regardless of their social interaction network. For Fig. 3 , for instance, this value is 3%. Larger values of α narrow the set of parameters that can be chosen in order to realistically describe a real system.
Finally, we compare the model's results with data of drinking abusive consumption in Brazil [26] . Data were collected from 2009 to 2019, thus in Fig. 4 the initial time t = 0 represents the fraction of abusive drinkers for 2009, t = 1 represents the fraction for 2010, and so on. Since the data is for the fraction of people that consume alcohol abusively, we plot the density of risk drinkers r(t) together with the data. In order to compare them with the model, we considered for the initial density of risk drinkers r(0) = 0.185 and numerically integrated Eqs. (1) -(3). The value 0.185 was chosen since it is the fraction of abusive drinkers for 2009 obtained from the database [26] . In addition, we rescaled the time of the simulation results to match the time of real data: such simulation time was multiplied by 0.12 for a better comparison. We find that the simulated drinking trajectories qualitatively correspond to the data. For the numerical results, we considered the parameters β = 0.06, γ = 0.11, α = 0.047 and δ = 0.2, which indicate that the probability of finding an individual that will spontaneously become a risk drinker in Brazil during the last decade is around 4,7%. Furthermore, looking at Eqs. (17) and (18) , it is easy to see that in order to model Brazil's c < β < β (2) c , showing that the model describes the available data in its most realistic spectrum (region II of Figure 3 ). Naturally, in comparison with actual data, models should always present the three different phases, i.e. coexistence between the three different population groups, since descriptions with only nonconsumers or risk drinkers are unrealistic. This qualitative agreement with Brazil's database in the realistic spectrum of the model points to a good, albeit simplistic, modelling.
In this work, we have studied a compartmental model that aims to describe the evolution of drinking behavior in an adult population. We considered a fully-connected population that is divided in three compartments, namely Susceptible individuals S (nonconsumers), Moderated drinkers M and Risk drinkers R. The transitions among the compartments are ruled by probabilities, representing the social interactions among individuals, as well as spontaneous decisions, in particular from moderate evolving into risk drinkers, and we studied the model through analytical and numerical calculations.
From the theoretical point of view, the model is of interest of Statistical Physics since we observed the occurrence of two distinct nonequilibrium phase transitions. These transitions separate the model in three regions: (I) existence of nonconsumers only; (II) coexistence of the three compartments and (III) existence of risk drinkers only. Regions I and III represent two distinct absorbing phases, since the system becomes frozen due to the existence of only one subpopulation for each case -this means that, in order to describe real populational systems, the parameters must be chosen so that the model falls in region II, since populations consisting solely of nonconsumers or risk drinkers do not represent a realistic entity. The critical points of such transitions were obtained analytically. A comparison with available data for Brazil's extreme alcohol consumption for the past decade shows a good qualitative agreement with the model, with the chosen parameters framed within its realistic boundaries. It will be important in a couple of years time to re-evaluate these results in the light of new data comprising years 2020 and 2021, in order to verify the direct effects of the COVID-19 pandemic in the Brazilian population's alcohol consumption. An hypothesis to be tested is a possible increase in parameter α combined with a corresponding decrease in the other parameters, corresponding to social interactions.
The phase transitions observed in the model are active-absorbing phase transitions, and the predicted critical exponent for the order parameter is 1 (m ∼ (β−β c ) 1 ) as in the mean-field directed percolation, that is the prototype of a phase transition to an absorbing state [30, 31] . It would be interesting to estimate numerically other critical exponents of the model, as well as to simulate it in regular d-dimensional lattices (e.g. square and cubic) in order to obtain all the critical exponents. This is important to define precisely the universality class of the model, as well as its upper critical dimension. This extension is left for a future work. Furthermore, it can also be considered the inclusion of heterogeneities in the population, like agents' conviction [32] , time-dependent transition rates [33] , inflexibility [34] , etc. 