A non-standard numerical scheme for an alcohol-abuse model with induced-complications

The prevalence of alcohol-related fatalities worldwide is on the ascendancy not only Ghana, but worldwide. Although the ramifications of alcohol consumption have been the subject of several studies, alcoholism remains a serious concern in public health. This study investigates the dynamics of alcoholism in a population with consumption-induced complications using a deterministic Modelling framework. Using a novel technique, we determined a threshold parameter R0 which we call the basic alcohol-abuse initiation number which is similar to the basic reproduction number for infectious diseases. The model has two mutually-exclusive fixed points whose existence depend on whether or not the R0 is less or greater than unity. Global asymptotic stability of the alcohol-abuse-free fixed point is shown to be associated with R0≤1. Further, forward bifurcation is observed to occur at R0=1, indicating the possibility of eradication of the phenomenon of alcoholism if R0 can be kept below unity over a sufficiently long period of time. Sensitivity analysis also revealed that the probability of initiation into alcohol-abuse by moderate drinkers (β1), followed by the probability of initiation into alcohol-abuse by heavy drinkers (β2) are the most the parameters with the most influence on R0 and consequently on alcohol-abuse persistence. A non-standard finite difference scheme is also developed to numerically simulate the model so as to demonstrate the findings derived from the analysis and also to observe the impact of some epidemiological factors on the dynamics of alcohol-abuse.


Introduction
Alcohol drinking is one of the most dangerous activities that young people erroneously regard as fun or enjoyment.Besides, it was falsely recommended in 2019 that, alcohol was very effective in the treatment of corona virus  and this led to a lot of alcohol poisoning [1].The initiation into alcohol-abuse should be considered as a public concern because there is no such thing as a safe consumption level [2].The rate of alcohol consumption during the outbreak of COVID-19 was a little higher for males as compared to the previous 2-3 years and significantly higher in females [3].It is estimated that alcohol consumption is responsible for the death of about 2.5 million people globally each year [4].The World Health Organization (WHO) indicates that, alcoholabuse is expected to cause approximately 60 types of disorders including liver cancer, liver cirrhosis, heart disease, esophageal cancer, homicide, pancreatitis, stroke, epileptic, neuropathy disorders, brain disorder, fatty liver problems, hepatic, type 2 diabetes, hypertension and dyspepsia.

E.A.B. Sandow, B. Seidu and S. Abagna
Mathematics has played an important role in researching into many disciplines.Among these are its application to disease control.Researchers have developed many mathematical models to predict and explain the effects of epidemics.However, the world still battles with a lot of diseases which are caused by either natural or artificial factors.All these explain why mathematicians have to do more research in building mathematical models to help reduce the risk of being exposed to epidemic diseases.Mathematical models have demonstrated high potential in helping to improve our knowledge of the dynamics of infectious agents and diseases.With alcoholism being a public health concern, many researchers have concentrated on building models to study the dynamics of alcoholism.Sa'nchez et al. [5] proposed an alcohol drinking model with recovery relapse and described the transmission dynamics of alcohol within the context of the classic SIR model.Huo and Song [6] presented an alcohol-abuse model to study binge drinking among youth.Another deterministic alcohol-abuse model was in [7] in order to gain more insight into this health concern and social phenomenon and to determine the most cost-effective way of addressing the canker.Xiang et al. [8] studied the global stability of an alcohol-use model with public health awareness campaign.Khajji et al. [9] developed a model with liver complications that specifically deals with educating alcohol-abusers and giving treatment to those with disorders.Pe'rez [10] presented a non-linear differential equation model, describing the difficulties and problems associated with alcohol consumption habits among college students in Colombia.Khajji et al. [11] proposed a model in order to look at the interaction between different groups of people with different levels of alcohol consumption.They focused on treatment of heavy alcoholics drinkers.In [12] a discrete mathematical model was presented to elucidate the interaction between the various compartments of drinkers.They also sought to determine optimal strategies for the minimization of the population of drinkers and maximization of acceptance by heavy drinkers to seek treatment.On the basis of the above research works and several others, the purpose of this study is to put forth a mathematical model, suitable for describing the dynamics of alcohol abuse to include the quitting and recovery rate moderate alcohol users, which is often neglected by most research on alcohol-use models.In this study, the effect of moderate drinkers, heavy drinkers, and individuals with complications recovering is discussed and also, the likelihood of the recovered individuals rejoining the chain of drinking is considered.Most researches involving infectious diseases and substance-abuse often use explicit methods to numerically solve the models.However, it is known that explicit methods such as the Euler and Runge-Kutta methods have inherent challenges including poor stability and efficiencies, which may give wrong dynamical behavior which are not inherent in the studied dynamical systems [13].For instance, use of the explicit methods may indicate bifurcations and chaos which the system may not really posses.The implicit methods have the ability to avoid these contrived dynamics while maintaining their efficiency and stability.This makes the implicit methods a better choice over their explicit counterpart.The non-standard finite difference methods developed in [14,15] have motivated several other schemes and the use of same in solving numerical schemes.Following these works, a non-standard finite difference scheme is developed to use in the numerical experimentation on the model.
We structure the remainder of the paper in the following manner: Section 2 presents the development of the Mathematical model in focus.Qualitative analysis including stability of equilibrium points, bifurcation analysis, and sensitivity analysis of the model are presented in Section 3. In section 4, we develop a non-standard numerical system to be used for simulation of the model.In Section 5, the model is solved using the developed numerical scheme and the discussion of findings of the study are presented.Finally, the conclusions of the study are captured in section 6.

Development of the alcohol-abuse model
In this section, the deterministic mathematical model of interest is developed to elucidate the dynamics of alcoholism with induced-complications.We categorize the population into five distinct and non-overlapping compartments: The compartment of Potential drinkers (P) denotes people who do not consume alcohol but are prone to adopting alcohol consumption behavior; Moderate alcohol abusers (M) represents individuals who slightly consume alcohol and have the ability to control the intake of alcohol; Heavy alcohol abusers (H) denotes people who are obsessive to alcohol consumption (unable to control the intake of alcohol) and; Alcohol-abusers who have developed alcohol-induced complications (C) representing people who develop disorders such as liver cirrhosis, brain disorders, neuropathy complications, fibrosis, pancreatitis, and alcoholic hepatitis.We assume that all inflows into the population are Potential drinkers at the rate of .We define  = (similar to force of infection) as the rate at which potential drinkers adopt alcohol-consumption behavior due to contact with alcohol-consumers.The parameter,  1 is the probability that a potential drinker starts to abuse alcohol as a result of interaction with a moderate drinker, and  2 is the probability that potential drinker begins to abuse alcohol due to interaction with a heavy drinker.The moderate drinkers are assumed to progress into heavy drinking at a rate of  and they quit drinking alcohol at the rate of  3 .Heavy drinkers are assumed to develop complications at the rate of  1 and quit alcohol consumption at the rate of  2 .Individuals with induced-complications are assumed to recover at the rate of  1 .Persons who recovered from alcohol-abuse are assumed to rejoin the compartment of potential drinkers at the rate of  2 per unit time.The parameters ,  1 , and  2 respectively represent death rate due to nature, death rate do to alcohol-abuse among heavy drinkers, and death rate associated with complications caused by alcohol-abuse.The general dynamics of alcohol abuse with induced-complications are depicted in the schematic diagram in Fig. 1.
The alcohol-abuse model described so far is thus presented by the following coupled ordinary differential equations.In subsequent discussions, where necessary, the following conventions will be used.
The description of the model parameters is summarized in Table 1.The baseline values of the parameters used later for simulation purposes are also presented in Table 1.Some qualitative properties of the model are presented in the following section.

Qualitative properties
It is easy to observe that the functions in the right-hand-side of (1) are Schlitz continuous, which is a necessary condition for Picard's existence theorem [16].Therefore, we can conclude that, the alcohol-abuse model (1) has a unique solution.
The epidemiological and mathematical well-posedness of the model are presented in Theorem 3.1.
Theorem 3.1.With non-negative initial conditions, the solution of the alcohol-abuse model (1) are always non-negatives for all time,  > 0.
Proof.Assume that for some time  1 with  1 >  > 0 we have  ( 1 ) = 0, d d Then from the first equation of (1), we have d (  1 ) d =  > 0, which contradicts our assumption.Therefore  () > 0. Similar arguments show that the remainder of the state variables are also non-negative, concluding the proof.□ The set of epidemiologically-reasonable solutions of the alcohol-abuse model ( 1) is named as the feasible region, wherein all solutions of the model with reasonable initial conditions are contained.We define this set in Lemma 3.1.

Lemma 3.1. Future solution of the alcohol-abuse model (1) with non-negative initial conditions are all non-negative and contained in
, which is positively invariant.
Proof.Adding all equations in (1) gives whose solution is given by Therefore, lim sup →∞  () ≤   .Therefore, the solutions of the alcohol-abuse model ( 1) are uniformly bounded in Ω, completing the proof.□

Equilibrium points and the basic reproduction number of the alcohol-abuse model (1)
The model ( 1) has two fixed points; the alcohol-abuse-free equilibrium ( 0 ) and the alcohol-abuse-persistent equilibrium  * .
The basic reproduction number denoted as  0 is an epidemiological quantity used to describe the mean number of secondary infections that arises from the introduction of a single infectious agent into an otherwise completely Susceptible population.In this study,  0 is named basic alcohol-abuse initiation number.We define  0 as the average number of alcohol-consumers that are initiated by a single alcohol-abuser introduced into an otherwise Potential drinkers population.In this paper, we adopt the Jacobian-Determinant method of [17] to determine  0 of the alcohol-abuse model (1).Following [17], the infected sub-system is given as follows: The Jacobian of the infected sub-system (2) evaluated at  0 is given by , whose determinant is given by Following the formulation in [17], the basic alcohol-abuse initiation number is given by The alcohol-abuse persistent equilibrium point can be shown to be given by  * = ( * ,  * ,  * ,  * ,  * ), where

Local stability of fixed points of the alcohol-abuse model (1)
Theorems 3.2 and 3.3 present the results on local stability of the fixed points of the model (1).The Lyapunov second technique is employed to study the local stability of the equilibria.
The Jacobian matrix of the alcohol-abuse model ( 1) is given by The characteristic polynomial corresponding to the alcohol-abuse-free critical point  0 is given by Now, the coefficients of equation ( 4) are positive whenever  0 < 1, so that by the Routh Hurwitz criterion, all the zeros of the polynomial in equation (4) will have negative real parts.This establishes the following result.
The characteristic polynomial of  evaluated at  * is given by; where, The following results follows from application of the Routh-Hurwitz criteria on equation ( 5).
Theorem 3.3.The following conditions guarantee the local asymptotic stability of the alcohol-abuse-persistent fixed point  * of model (1): > 0.

Global stability of  𝐴0
In this brief, the global stability of the alcohol-abuse-free fixed point  0 using the first technique of Lyapunov.To determine the Lyapunov function for the model, we adopted the technique described in [17] as follows: Since the infected sub-system consists of  ,  and , we use a Lyapunov function candidate give by Recall that, .
Following [17], the infectivity vector for this model is given by  = (  1 ,  2 , 0 ) , so that the Lyapunov coefficients are obtained as follows.
1 is obtained by replacing the first row of   (  0 ) with  to get.
2 is obtained by replacing the second row of   (  0 ) with  to get.
) with  to get.
Therefore, we define the Lyapunov function as Then From the Lasalle's invariance principle [18], every solution of (1) with initial conditions in Ω converges to  0 as  → ∞ whenever  0 ≤ 1.This establishes the following result.

Bifurcation analysis
(from  0 = 1) be the bifurcation parameter.Then  evaluated  0 has a simple eigenvalue associated the following right and left eigenvectors: where ,  4 = 0,  5 = 0.
Note that  2 and  2 can be found using ).
The bifurcation coefficients are thus obtained (using Clearly,  < 0 and  > 0 and hence, from Theorem 4.1 of [19], the model exhibits forward bifurcation at  0 = 1 as illustrated in Fig. 2.

Sensitivity analysis
We conducted sensitivity analysis to determine the impact of small changes in model parameter values on ( 0 ) and the state variables at the endemic equilibrium,  * .The sensitivity index of variable  with respect to a model parameter  is computed using the Normalized Forward Sensitivity Index defined as: The sensitivity indices of  * (presented in ( 3)) and the  0 were calculated and evaluated using parameter values in Table 1 and presented in Table 2.It is observed from Table 2 that, all parameters but two ( 1 and  2 ) have an inverse relationship with the basic alcohol-abuse initiation number ( 0 ).Thus, increasing (or decreasing) these parameters leads to a decrease (or an increase) in  0 .Similarly, increasing (or decreasing)  1 and  2 leads to an increase (or a decrease) in  0 .The parameters with the most affect on  0 are  followed by  1 and then  2 .Thus, in order to reduce  0 , efforts should be made to increase  and reduce  1 and  2 .As expected, increased  0 is associated with disease progression and hence the relationship between  0 and model parameters are expected to be the same for the infected compartments and model parameters.This is illustrated by the positive signs of the sensitivity indices of  0 ,  * and  * with respect to  1 and  2 .A counter-intuitive result is however observed for  which is observed to have an inverse relationship with  0 but a positive relationship with  * and  * .Also the sensitivity index with respect to  1 seems to be counter-intuitive in the case of  * .However, this can be explained as follows: An increase in alcohol-induced deaths ( 1 ) among heavy-drinkers reduces the population of heavy drinkers and consequently the basic initiation number.This is likely to lead to more and more people changing their behavior towards heavy drinking to moderate drinking, causing the number of moderate drinkers to increase.
The relationship between some model parameters and  0 using surface plots are depicted in Fig. 3. Fig. 3 graphically shows how the model parameters (in pairs) influence the basic alcohol-consumption-initiation number.It is observed in Fig. 3a that  0 linearly depends on the probabilities of adopting alcohol consumption behaviors due to contact with moderate drinkers ( 1 ) and heavy drinkers ( 2 ).As these parameters increase (decrease), the average number of new alcohol consumers will increase (respectively decrease).The other parameters are observed to have a nonlinear impact on  0 (see Figs. 3b, 3c, 3d, 3e, 3f, 3g, 3h, and 3i).This suggests that, attempts to control the dynamics of alcoholism through vary these parameters needs a closer look than doing so by varying  1 and  2 .

Development of the non-standard finite difference scheme (NFDS)
Numerical schemes used to solve differential equations are many and are often chosen based on some conditions which may include convergence, speed and ease of application.The most used method for solving ordinary differential equations models is the fourth-order Runge-Kutta method.This method is part of a family of explicit methods which have been shown to have inherent challenges such as poor stability and inefficiency.To circumvent these challenges, implicit numerical schemes may be used.In this section, a non-standard finite difference scheme that will be used to solve the system is developed (see [15] for further details on the method).
Since, all equations in (1) are in  1 (ℝ), the time derivatives are approximated using the forward difference scheme given by In the equation above, the function  (ℎ) is called the denominator function, which must satisfy  (ℎ) = 1 −  −ℎ .The model (1) can be approximated as follows: Rearranging the equations in ( 6) gives; ) , Now, since the equations in (7) have factors of the form 1 +  ℎ, we choose denominator functions of the form (ℎ, ) = ( 1 −  − ℎ ) ∕ so that the NFDS above can be written as follows [14]:

Numerical results and discussions
In this section, the non-standard finite difference scheme (8) developed in section 4 is used to solve the model (1).The essence of the numerical solution is to illustrate the analytical findings and also to conduct a deeper examination of how model parameters influence the prevalence of alcohol abuse within the population.The model parameter values in Table 1 were used to carry out the simulations.
To illustrate the superiority of the non-standard finite difference scheme, we solved the model using three different methods; The forward Euler method, the Runge-Kutta fourth order scheme, and the developed NSFDS.The results are presented in Fig. 4. It is observed that the developed non-standard finite difference scheme outperforms the Euler's scheme both for large step sizes (see Fig. 4b) and very small step sizes (see Fig. 4a).The results for our method and the RK method are not so different.However, we note that our method is computationally cheaper than the RK method.This gives us a good reason to adopt our developed method for the simulation.
Fig. 5 illustrates the impact of alcohol abuse behavior adoption caused by the interaction between potential and moderate drinkers.To study this, we varied the parameter  1 and plotted the results for the various values of the transmission parameters.From Fig. 5a, it is observed that increasing  1 leads to decrease in the potential drinkers population.This is as a result of the fact that increasing  1 increases the chance of people adopting alcohol consumption behaviors and subsequently the population of the potential drinkers.It is also observed that, increasing the probability of acquiring drinking habits through interaction with drinkers leads to increase in the number of moderate drinker (see Fig. 5b) and heavy drinkers (see Fig. 5c) and subsequently the number of individuals who develop alcohol-consumption-related complications (see Fig. 5d).The dynamics are however observed to switch for a while after some time where the number of drinkers is reduced for increased values of  1 .This switch can be attributed to the fact that as more individuals develop complications the number of drinkers available to introduce potential drinkers into alcohol consumption will reduce and consequently the number of drinkers will reduce.This explains why it may be impossible for alcoholism to be adopted by the whole population.
The impact of probability of alcohol-abuse adoption caused by heavy drinkers is presented in Fig. 6.We observe that increasing the probability of alcohol-consumption transference from heavy drinkers increases the moderate drinkers (see Fig. 6b), heavy drinkers (see Fig. 5c) and individuals with complications (see Fig. 6d).These increases in the alcohol consumers populations lead to a corresponding decrease in the potential drinkers population (see Fig. 6a) as  2 increases.We therefore note that, efforts should be made to keep the chances of individuals involved in alcoholism influencing others from engaging in the habit of drinking.This is more especially important in the case of minors in certain parts of the developing world who may often be sent on errands to buy alcoholic beverages.It is found in Table 2 that when the rates of recoveries ( 2 and  3 ) are increased, the basic initiation number decreased.This implies that, moderate and heavy drinkers should be given the necessary support in order to get them fully recovered.This will go a long way to curb the spread of alcohol abuse infections in the society.
The impact of the rate of transition from moderate to heavy drinkers was also studied by solving the model for varying values of the associated parameter, .It can be observed from Fig. 7 that, when () is increased, the number of potential drinker and moderate drinkers decrease (see Fig. 7a and Fig. 7b) while the number of heavy drinkers and those with complications increase (see Fig. 7c and Fig. 7d).Thus, measures need to be taken to ensure that moderate drinking does not progress into heavy drinking.This can be achieved by ensuring that the frequency of drinking be reduced to the barest minimum.Also, the phenomenon of resorting to drinking to overcome problems or stress need to be curbed through counseling as that has the potential of resulting into heavy drinking if the problems persist.

Conclusions
In this paper, we put forth a nonlinear differential equations model to elucidate the dynamics of alcoholism in a population.The population is split into five mutually-exclusive compartments of Potential drinker's moderate drinkers, heavy drinkers and individuals who have developed complications due to excessive consumption of alcohol.The model is developed by treating alcoholism as a contagion and following the theory of disease transmission from infected persons to Susceptible persons, consider the process of drinkers influencing non-drinkers as a transmission using standard incidence function.We adopted a recently developed technique [17] to determine basic alcohol-abuse initiation number which is similar to the basic reproduction number to describe the average number of secondary alcohol consumers that are initiated into alcoholism by a single alcoholic that is introduced into an otherwise completely non-drinkers population.The alcohol-abuse-free equilibrium is shown to be globally asymptotically stable whenever the basic alcohol-abuse initiation number is not more that 1.This global stability results was further confirmed by the existence of forward bifurcation in the model at  0 = 1.The sensitivity analysis revealed that, the most sensitive parameter to the  0 is rate of transition into heavy drinking  followed by probability of initiation per contact with a moderate drinker  1 and then the probability of initiation per contact with a heavy drinker  2 .In order to conduct simulation on the model, we developed and adopted an implicit non-standard numerical scheme which overcomes the inherent challenges of standard explicit schemes.The numerical simulations confirmed that, the spread of alcohol infections and alcohol-induced complications can be minimized by lowering the probabilities of initiation per contact between potential drinkers and alcohol abusers.
There are multiple potential extensions of this study: The incidence functions have been shown to have interesting effect on the dynamics of models.While this work adopted a standard incidence function to describe initiation into alcoholism, future research could adopt saturated incidence functions as in [20,21, and related works].Recently, fractional and fractional-fractal operators have been adopted in the development epidemic models.These models have been said to be more realistic than their integer-order counterparts due to the inherent memory property the fractional differential operator.The current can be extended by incorporating fractional and fractal-fractional operators as in [22,23, and similar works].Also, a discrete-time model versions of the model can be

with d d𝑡 = 0
if and only if  =  = 0. Furthermore, the alcohol-abuse-free equilibrium point of the model (1),  0 is the largest invariant compact set in { ( , , , , )

Fig. 3 .
Fig. 3. Surface plots showing the relationship between  0 and some model parameters.

Fig. 4 .
Fig. 4. Plots of heavy drinkers population using the NSFDS and two other numerical schemes for different steps sizes.

Fig. 5 .
Fig. 5. Time series plots showing the impact of alcohol-abuse transference ( 1 ) on population dynamics.

Table 1
Description and baseline values of the model parameters.

Table 2
Sensitivity indexes of  0 and state variables of at the endemic equilibrium.