Mathematical models for the dynamics of alcohol related health risks with changing behavior via cultural beliefs in Tanzania

Alcoholism has continually posed health challenges in many communities for decades. In this paper, a more realistic model for health related risks associated with alcoholism is formulated. It considers a population proportion that has social cultural protection from alcohol consumption. In the context of this paper, such protection emanates from religious beliefs. The Next Generation Matrix (NGM) approach is used to compute the basic risk reproduction number. The risk free equilibrium point is proved to be globally asymptotically stable whenever the basic risk reproduction is less than unity and unstable otherwise. The sensitivity analysis of the basic risk reproduction number and numerical simulation results reveal that for effective control of the health risk problem in the community, the deliberate intervention strategies and policies should focus on discouraging alcoholic behaviors on its onset during initiation stage than focusing other population proportions already at risk.


INTRODUCTION
Alcoholic beverages are an integral part of cultures across the globe with wide use in rituals, societal artifacts and festivals [1,2]. In essence alcoholic beverages are important for both social and economic reasons. Precisely, alcohol beverages are a source of income for the livelihood of the rural community and have health benefits such as prevention of thrombosis when consumed at desired levels [3]. However, despite the aforementioned health and economic benefits of alcoholic beverages; when taken to undesired level it can accelerate to alcoholism, a behavior which poses serious health challenges to consumers. Alcoholism may be defined as the state of addiction to the consumption of alcoholic drinks which eventually turns into the state of alcohol dependency. This is a condition in which a person has a physical or psychological dependence on drinks that contain alcohol. It is a precursor to injury and violence and its negative impacts can spread throughout a community or a country, and beyond, influencing levels and patterns of alcohol consumption across borders [4].
The common symptoms of alcoholism include, but not limited to: strong compulsion to drink; inability to limit ones drinking in any given time; physical dependence; increased uptake of alcohol for optimum effects; and problems associated with alcoholism − injuries, receives multiple drunken driving citations, frequent arguments and poor relationships in families.
The drinking limits or threshold to be referred to as an alcoholic is estimated to a maximum of 21 standard bottles per week for men and 14 drinks per week for women [2,3,5]. This quantity defines the health tolerable amount of alcohol consumption at which an individual attains the most vulnerable state of exposure to alcohol associated health risks.
According to the World Health Organization report, at least three (3) million deaths and 132.6 disabilities due to harmful use of alcohol reported annually. Comparatively, over-consumption of alcohol causes more harm to human health than tuberculosis, HIV/AIDS and diabetes [6].
Further, alcoholism is linked to liver cirrhosis [2,7] and accounts for 10% of the total disease burden attributable to harmful use of alcohol [8]. Owing to the available literature, it was observed that alcohol consumption increases risks of cancer infections, which is attributable to 20% of alcohol related deaths. Common types of cancers related to alcohol overconsumption are cancers of the mouth, oropharynx, esophagus, colorectal, liver, and breast [7,8,9]. Some cultural practices have promoted positive drinking with the aim of reducing health and socio-economic risks. Part of renown active cultural practices are religious beliefs which have strong influence in people's lives. It plays an important role in promoting the health and molding individual's behaviors [10,11]. By the fact that religious communities at least discourages irresponsible alcohol drinking by pronouncing it a sinful behavior. It is therefore important to focus on the impacts of religious communities when dealing with the dynamics of health related risks associated with alcoholism because they provide social protection to a proportion of the population from alcohol consumption or alcohol abuse.
In quest to provide insights on effects and the spread of alcoholism, different mathematical alcohol epidemic models have been formulated and analyzed for the purpose of understanding the effects of alcohol consumption on health and socio-economic aspects of the society. In [12], for example, a link between drug addiction and infectious disease is considered and the importance of using dynamic models to predict trends and generate estimates is emphasized especially where data are sparse. Mundt and collaborators in [13] fitted a stochastic actor − based model to Add Health data of 7th through 11th grade U.S students enrolled between 1995 and 1996.
The study sought to determine factors that influence the dynamic interplay between adolescent friendships and substance abuse and revealed that peer selection and friendship making has a greater role to play in alcoholic behavior among adolescent. Similarly, [14] presented SAT Q− type alcoholic model with two control strategies. Two objective functions for alcohol quitting and cost of controlling alcohol were proposed and analyzed with the help of Pontryagin's Maximum Principle approach. Numerical simulations results recommended reduction in the number of alcoholics and the increase in the number of susceptible as the better control option.
Xiang and collaborators, [15] studied the effects of public health educational campaigns on drinking dynamics. The study explored the use of Lyapunov functional to establish the global stability of the model equilibria. Numerical simulations result proved that public health education campaigns of drinking individuals can slow down the drinking dynamics. In [16] a non − linear SHT R mathematical model for drinking epidemics is presented and analyzed. In this model, the drinking population were placed in heavy drinkers (H) and drinkers in treatment (T ) compartments assuming that after treatment a treated heavy drinker will only become a heavy drinker again after passing through a recovery and susceptible compartments in respective order. The analytical results revealed that drinking epidemic can be controlled by the combined efforts of reducing the contact rate between the non-drinkers and heavy drinkers to refrain from drinking; and increasing the number of drinkers that go into treatment.
In [7], an SIRS alcoholism models with relapse on a weighted networks was formulated to study the impact of the fixed weight and adaptive weight on the spread of alcoholism. The study results established that, if the proportion of recuperator to accept treatment is equal to that of susceptible people to refuse alcohol drinking, then preventing the susceptible people from alcohol drinking become more effective. In [9], a deterministic model to study the spread of alcoholism was formulated by dividing the alcohol drinking population into moderate and alcoholic population. The comparison of alcohol control approaches targeting different alcohol drinking populations revealed that encouraging and supporting moderate drinkers to quit alcoholic consumption produces better results than targeting the alcoholic group only.
Mushanyu and Nyabadza in [17] presented a risk − structured model and used it to understand the phenomenon of the spread of drug abuse with in-patients treatment programs. Analysis of the risk abuse and numerical simulations suggested that education and skills to deal with risky situation may better equip individuals to stand against initiating into drug abuse.
Mayengo and collaborators in [18] suggested that for alcohol related health risks modeling to be more relevant, two aspects of the model formulation are proposed: the influence of external motivations with positive and negative influential effects, and a clearly defined alcohol drinking population compartments have to be included. The crisp model developed in [19] have definitely answered this concerns where the fuzzy logic analysis of the model is carried out. A non-drinker acquire alcohol drinking habits through social contacts [9,19,20] at the force of peer influence where β is the probability of initiation to alcohol drinking, θ 1 and θ 2 are the modification parameters (θ 1 ≤ θ 2 ), c is the effective contact to influence one into alcohol drinking, and N is the total population. The protected class is increased by conversion to cultural and religious beliefs of the susceptible, low risk drinkers, moderate risk drinkers and alcoholics. It is decreased by backsliding at constant rate γ 2 . Low and moderate risk drinkers are increased by initiation to alcohol drinking of the susceptible class. Low risk drinkers progress to moderate drinkers at a constant rate σ . Moderate risk drinkers progress to alcohol addiction at a constant rate δ , recover at a constant rate ξ and are protected at a constant rate τ. Alcohol addicts either recover at the rate η or get protected at the rate ψ. The recovered class become susceptible again at the rate ω. All the classes are subjected to reduction due to natural causes at a constant rate µ.
The proportion ρ ∈ (0, 1) of susceptible individuals are recruited via peer influence into low risk drinking class, while the remaining proportion, (1 − ρ) enters the moderate risk class.
Based on the dynamics on the model compartment in Fig. 1, the following set of equations are (2) Let N be the total population size given by; Then the equation describing changes in the total population is given by, Basic properties of the model solution. In this section, the basic results of solutions for model system (2) is discussed. These properties lay down a foundation of proofs of stability analysis results of the model. Following [21], let x ∈ R n + denote the set of state variables x = (x 1 , x 2 , . . . , x n ) with the positive components, the following results may be established.
. . , f n (x)) with x ∈ R n + , be continuous and there exist a continuous partial derivatives Theorem 2.1. Let f : R n + → R n be a locally Lipschitz continuous and for each i = 1, 2, . . . , n Proof. Applying Theorem 2.1, we define and f 6 (x) =Ṙ where x = (S, P, L, M, A, R). By the properties of continuity over operations, we have continuity of f i for all i = 1, 2, ..., 6. Furthermore, where the derivatives are continuous. Thus we have, By Lemma 1, we know that f is locally Lipschitz continuous. Let Similarly, repeating the procedures for the rest of state variables, the following conditions are established 3. Invariant region. Since the above model in (2) represents human population, it is assumed that all the state variables and parameters of the model are non-negative for all t ≥ 0. Now, from the equation (4) where N 0 is the initial population obtained by evaluating equation (3) at the initial conditions of the respective state variables. It follows that, for µ < π as t → ∞, Hence, the region Ω is positively-invariant under the flow induced by system (2), that it is well posed mathematically and epidemiologically, sufficiently for the dynamics of the flow generated by system (2) to be considered in Ω.

Steady state solutions.
To obtain the equilibrium point, we solve a non-linear equations. (14).

2.4.2.
Risk endemic state. The solution N * − R H S * = 0 from equation (17) leads to endemic state. From system (14) we note that The total population is given by Re-writing the first two equations of the model system (2) in terms of λ * and S * with appropriate substitution of other state variables. Eliminating S * from the subsystem, equation (21) is established with an implication of forward bifurcation.
It follows the establishment of a non-trivial solution E 1 for the risk endemic as 2.5. Basic risk reproduction number. In this section, we compute the basic risk reproduction number of the model using Next Generation Matrix (NGM) method by [19,23,24]. The basic risk reproduction number, denoted as R 0 , may be defined as an average number of secondary risk individuals produced by a single risk individual in an entirely susceptible population during his/her risk duration. It serves as an indicator used to predict the possibility of the occurrence of risk epidemic. Based on [24], the basic risk reproduction number is given as the dominant eigenvalue or spectral radius of the next generation matrix. Considering the system of equations (2), this system has three risky states at different risk levels, L(t), M(t) and A(t); and three non-risky states, S(t), P(t) and R(t). Compositing system (2), we state the risk system as and non-risky system as The matrix F corresponds to "transmissions" and the matrix V to "transitions" of risks factors through different states with different risk levels. Following [23,24], referring to the risk states with indices i and j, with i, j ∈ {1, 2, 3}, the entry F i j is the rate at which individuals in a risk state j give rise to individuals in risk state i, in the linearlized system. Thus, F i j = 0 occurs only when there are no new cases in risk state i can be produced by an individual in risk state j. All epidemiological events that lead to new risks are incorporated in the model via F, and all other events via V .
The linearlization of risk system (23) at the risk free state E 0 yields the following F and V matrices, . Thus by direct computation, we have In epidemiological sense, the interpretation of the entry V −1 i j describes the expected time that an individual with risk state j will spend in a risk state i for the rest of his/her life. For instance, in a matrix V −1 above, individuals who are presently in state L will spend, on average, an amount  [24]. Now, we have the following Next Generation Matrix (28) To interpret the entries of the Next Generation Matrix FV −1 and develop a meaningful definition of R 0 we note that the entry FV −1 i j is the expected number of secondary risk individuals in compartment i produced by individuals initially in compartment j assuming that the drinking environment shared by the individuals remain homogeneous [25]. The basic risk reproduction number, R 0 , is given by the dominant eigenvalue of matrix FV −1 . Therefore, It is important to note that
3.1.1. Local stability. For the differential equations presented in (2), following [21,26] the following theorem can be established.
Theorem 3.1. The risk-free equilibrium E 0 is locally asymptotically stable (L.A.S) for R 0 < 0 and unstable otherwise. Proof. Consider the Jacobian matrix evaluated at the risk free equilibrium point bellow It follows that the trace and determinant of J(E 0 ) are respectively, given by which upon simplification we get Using trace determinant approach, the risk free equilibrium point E 0 is locally stable if and only if Tr(J(E 0 )) < 0 and Det(J(E 0 )) > 0. Since µ < π, it follows that, E 0 is asymptotically stable whenever R 0 < 1 and unstable when R 0 > 1 provided that

3.1.2.
Global stability conditions for the risk-free equilibrium. The roles R 0 in global stability analysis can be traced back from [26,27] who established two axioms ((H1) and (H2)), that need to be satisfied to guarantee the global asymptotic stability of the risk free state. Following [26] we present system (2) in the form: where X = (S, P, R) T ∈ R 3 whose components denote the number of individuals at risk free (non alcoholic individuals) and Z = (L, M, A) T ∈ R 3 whose components denote the number of alcoholic individuals at different risk levels. The coordinate (X * 0 , 0) ∈ R 6 denotes the risk free equilibrium for the model and (·) T denotes a vector transpose. Now with the system (2) in terms of X and Z, the following axioms combined guarantee globally asymptotically stability.
Theorem 3.2. The risk free equilibrium E 0 of the model (2) is globally asymptotically stable (g.a.s) of the system (2) provided that axioms (H1) and (H2) are satisfied.
Proof. We re-write the model equation (2) into f (X, Z) and g(X, Z) such that X = (S, P, R) , Z = (L, M, A). We have Evaluating the subsystem (33) at X * 0 , 0 we get which satisfies axiom (H1). Now, from the susbsystem (34) we re-write g(x, z) such that Since B is a Metzler matrix andg(X, Z) ≥ 0 provided that S(t) ≥ 1 and ρ ≤ 1, then axioms (H1) and (H2) are satisfied. It is clear that the equilibrium point E 0 is a globally asymptotically stable equilibrium. Hence, by the above theorem E 0 is globally asymptotically stable.  (2) and their sensitivity indices on R 0

Sensitivity analysis.
To determine the best solution in reduction of health risks and mortality associated to alcoholism, it is important to understand the relative importance of the different parameters responsible for risk transmission dynamics and prevalence [30]. We perform sensitivity analysis of sensitive parameters in order to determine the input parameters with the most contribution to the output variability of R 0 [7, 31, 32]. We do this based on understanding that, the initial spread of health risk associated to alcoholism in the model presented in (2) is directly related to the basic risk reproduction number R 0 [30].
Using normalized forward sensitivity index method -a partial differential technique -we calculate the sensitivity coefficients, ϒ R 0 y i , for each of the input parameters y i ∈ y for R 0 in equation (29) to the output variable R 0 [31,30]. We therefore have a vector y given as a set of parameters in which some of them serve as independent input for the corresponding dependent output R 0 . The sensitivity coefficients, defined as the measure of the relative change in the dependent variable when the independent variables change one at a time [30,31,32], explain the impact of each parameter value in the health risk transmission threshold. Researchers conduct sensitivity analyses for a number of reasons in an attempt to answer their research questions.
However, the robustness of the model predictions to parameter values is the main impetus of performing sensitivity analysis to many modelers [30,32].
Consider an explicit formula for R 0 given in equation (29), to each of the input parameters we derive an analytical expression for the sensitivity coefficient of R 0 , with respect to the parameter y i as where y i is the i th parameter as shown in Table 1. For example, the sensitivity index of R 0 with respect to β is given by Following the similar procedures we calculate the sensitivity coefficients of the rest of parameters and presents their results in Table 1.
Learning from the sensitivity coefficients in Table 1 we know that the most sensitive parameters in making significant changes in R 0 include: the natural mortality rate, µ; recruitment rate, π; measure of influence of risky individuals, β ; and the necessary contact rate between a susceptible member and a drinker required to convince the susceptible member to drink, c. The rest of the parameters have smaller sensitivity indices which may not require as much attention to estimate since small perturbation in those parameters lead to insignificant changes in the output variable [32].  Table 1 The sensitivity coefficients ϒ R 0 µ = −1.6901, and ϒ R 0 π = ϒ R 0 β = ϒ R 0 c = +1.0000 imply that, µ is negatively correlated with R 0 while π, c, and β are positively correlated with R 0 . If we increase (or decrease) µ by 10%, the resultant R 0 is also expected to decrease (or increase) by 16.901% of its original value. Similarly, if we increase (or decrease) either π, β or c by 10%, the value of R 0 also increases (or decreases) by 10%. However, despite the significant sensitivity coefficient of π, and µ, they may not be the suitable decision variables in this case.
Any deliberate efforts to increase the natural death rate or decrease the recruitment rate in order to reduce the intensity of the spread of the health risks in the community defeats the ethical requirements. We therefore remove π and µ from the list of targeted parameters.
According to the sensitivity analysis results for input parameters for R 0 , intervention strategies focusing on discouraging the drinking behaviors at its initiation stage is more effective than  Table 1 where the particular parameter is not considered as independent variable of R 0 .

Numerical analysis.
Where necessary and for the purpose of simulation we use parameter values in Table 1 Table 1. Consequently, it increases the recruitment of susceptible population proportion into low risk compartment (see Fig. 3). The variation of β shows no significant effects on the medium risk population compartment, so is the variation of c.
Since alcohol related risks is a staged process depending on the patterns and frequencies of alcohol consumption among others factors, as β increases it is expected that the population size   Table 1.
of L(t) will increase quantitatively. However, at higher values of β the population proportion of the low risk class converges to a common value within the period of one year. Similarly, we observe that the number of high risk drinkers decreases with the increase of β values for the first three years and later it increases with the increase of β . The effects of variation of the contact rates between risky and non risky individuals giving easy access to alcohol beverages follows similar patterns as that of β variation with β being more efficient in terms of time. That is to say, β based decisions will yield the desired results inn considerably shorter time than is the case for c based decisions.

CONCLUSION
This paper studies the dynamics of health risks associated with alcoholism in a community in which a proportion of susceptible population receives a social cultural protection acquired from the religious beliefs. Except for occasional drinkers, this work considers that some people voluntarily quit alcoholic consumption as their personal efforts to change their drinking status or upon medical grounds. The equilibria of the model system are found and their stability are analyzed. The basic risk reproduction number of the model is computed and the sensitivity analysis for its input parameters is established. The sensitivity analysis of R 0 suggested that discouraging alcohol drinking behavior at its initiation stage is more effective in the control and reduction of the health related risks associated with alcoholism than targeting other populations already at risk.
Numerical simulations are performed to illustrate the effects of the most sensitive parameters to various risky classes. The numerical simulations results are used to confirm the sensitivity analysis findings. For effective control of the health related risks associated with alcoholism in the community, public educational campaigns would do a better job. When the influential characters are involved in the public education campaign, better results are expected since societal influence is the key factor recruiting people into the behavior. Since alcohol drinking compromises the quality of decision making to its consumers, it can be closely associated with other irresponsible behaviors such as irresponsible sexual behaviors, the combination of alcohol related risks and the risk of contacting sexually transmitted diseases in one model will make an interesting study.

CONFLICT OF INTERESTS
The author(s) declare that there is no conflict of interests.