Modeling Hepatitis B and Alcohol Effects on Liver Cirrhosis Progression

1  Introduction

The late stage of liver disease, known as liver cirrhosis, occurs when the liver organ fails to function correctly due to the possible transformation of healthy liver tissue into scar tissue caused by alcohol processing stress or viral infections. In the world, cirrhosis is responsible for one million fatalities annually, indicating its significant morbidity and mortality rates [1]. Ultimately, a chronic or continual injury causes the liver to slowly languish and become incapable of functioning properly. A distinctive feature of the liver is its ability to regenerate; however, this process is also inhibited by liver cirrhosis. Hepatitis B and C virus infections, prolonged alcohol consumption, other liver illnesses and disorders can all lead to cirrhosis. The risk of developing cirrhosis is higher in those with chronic viral infections [2]. Although before the extensive liver damage, the liver cirrhosis has no apparent symptoms. Hepatitis B is a major worldwide health concern, caused by a virus called hepatitis B. Fever, fatigue, joint pain, loss of appetite, nausea, vomiting, and jaundice may occur two weeks later as symptoms. An incubation period of 45−180 days (average 60−90 days) meant that many of the infected would not experience symptoms. Mother-to-child transmission and contact with contaminated bodily fluids such as blood, saliva, vaginal secretions, and semen are potential transmission modes [3,4]. In 2022, the World Health Organization estimates that 254 million individuals have a chronic hepatitis B infection, with 1.2 million new cases reported annually. An estimated 1.1 million people died mostly from cirrhosis and hepatocellular carcinoma (HCC) as a result of hepatitis B in 2022 [5]. Another significant risk factor for the development of chronic liver cirrhosis is heavy alcohol use, with 38% of adults over 15 years, consuming more than 17 liters of pure alcohol per year [6]. The prevalence of cirrhosis is approximately 10 to 20 percent of heavy drinkers after ten or more years, especially in individuals taking at least 50 grams of pure alcohol daily over 10 to 20 years, and there is a high degree of risk of fast liver destruction when alcohol intake is combined with HBV infection [7]. However, a lot of studies have demonstrated that the risk of developing cirrhosis is not increased if people do not consume more than 50 grams of pure alcohol daily. Consuming more than 50 grams of pure alcohol per day, however, speeds up the development towards liver cirrhosis in chronic hepatitis B carriers for both men as well as women. It is poisonous and fatal if a person weighing 60 kg takes in more than 300 grams of alcohol daily [8–10]. Hepatitis B effective vaccine available with 95% effective antibodies [11]. Chronic infection treatment is crucial to reduce the risk of severe problems like liver cancer or cirrhosis. Treatment duration varies depending on the genotype and medication used, ranging from six months to a year [12].

In studying the dynamical behavior of epidemic models, the incidence rates play a vital role. Capasso and Serio [13], for the first time, proposed the idea of saturation incidence rate β(A)SI1+α2I. In a population that is completely susceptible to the illness, β(A)I1+α2I achieves the saturation level when I grows and α2I measures the force of infection after entering the disease. In this incidence rate 11+α2I measures the inhibitory impact of behavioral changes in susceptible people when the number of infected people rises or the crowding effect of infective individuals. Compared to bilinear, the saturated incidence rate is more general.

The dynamical behavior of infectious illnesses, such as cholera, liver cirrhosis, Hepatitis B infection, and others, is commonly studied using mathematical models. Apart from modelling of infectious diseases, various research has also been done for anticancer medication and locating numbers [14,15]. Mathematical models are considered to be the tool helping to discuss hypotheses, verify the experiments, and simulate the dynamics of complex objects. Zhou and Fan [16] created an SIR (Susceptible-Infective-Recovered) model for susceptibility, infection, and recovery analysis to evaluate human behavior during medical resource shortages. Ilhan and Sahin have used the Morgan-Voyce collocation method (M-VCM) to determine the approximate solution of the SIR model with the vaccination effect. Their paper provides an alternative criterion of the certainty of the approximate solutions. They primarily tried to determine the exact solutions of the SIR model of vaccination [17]. Researchers have developed multiple models to study relationships between alcohol abuse patterns alongside cirrhosis evolution and hepatitis B virus infection processes. Park et al. [18] investigated the factors linked to alcohol consumption among Hepatitis B carriers in Korea and concluded that hazardous alcohol consumption is defined as consuming more than sixty grams of alcohol on one occasion for males or forty grams of alcohol on several occasions for females. A continuous and a discrete mathematical model were developed by Khajji et al. [19] to investigate the dynamical behavior of alcohol use and the effects of both public and private addiction treatment centers. By utilizing a mathematical model and an optimum control technique [20,21] elucidated the dynamics of infectious diseases and came to the conclusion that the illness may be controlled by vaccination and treatment. Zhou et al. [22] compared moderate and excessive drinking and came to the conclusion that moderate alcohol consumption has no discernible impact on the progression of liver cirrhosis, whereas excessive drinking causes liver inflammation in patients with chronic Hepatitis B infection, ultimately accelerating the progression of the disease. A mathematical model developed by Dano et al. [23], consisting of four classes, namely susceptible, acutely infected, cirrhotic, and recovered, with a logistic model to study the combined effect of Hepatitis B infection and heavy alcohol consumption on liver cirrhosis progression dynamics.

The novelty of our research is considering the effect of exposed and treated classes. Since Hepatitis B has a latent phase, and treatment effects. With the inclusion of the exposed class, a latent period is introduced, which delays the transmission of individuals into the acutely infected class. We assumed that a fraction of individuals move from the exposed to the recovered compartment directly. The treated class is included to model the impact of medical interventions, which helps in analyzing how treatment reduces the infectious (both acutely infected and cirrhotic) population and reduces the progression to liver cirrhosis. Moreover, the existence of a unique solution for the proposed model is examined, which is ignored in most of the papers. Furthermore, we used a saturated incidence rate to study the saturation effects. Two basic reproductive numbers are calculated when the minimum and maximum amount of pure alcohol is consumed, respectively. The 3D analysis of these basic reproductive numbers are also carried out. To the author’s utmost knowledge, no one has yet considered the same.

We organized this research study as follows. The model formulation, description, as well as associated assumptions are presented in Section 2. In Section 3, the qualitative analysis of the proposed model is presented. That is, the solution is of existence, uniqueness, positivity, as well as boundedness. Moreover, the feasible region is defined, equilibrium points, the basic reproductive number with 3D-type simulation, and local and global asymptotic stability at equilibria are calculated. In Section 4, the bifurcation phenomena are analyzed, where we proved the existence of forward bifurcation of the proposed model. The sensitivity analysis is carried out in Section 5, while the numerical simulations are explored in Section 7, and the conclusion is presented in Section 8.

2  Model Formulation

The total human population is denoted by N(t) and partitioned into six different compartments namely the susceptible, exposed, acutely infected, liver cirrhotic, treated and recovered denoted by S(t), E(t), I(t), C(t), T(t) and R(t), respectively. The susceptible individuals S(t) include those who are at risk of becoming infected but have not yet been infected at time t. People in the susceptible group are not immune to a possible infection. The exposed compartment E(t) includes those individuals who are infected but not yet infectious, meaning that they are in the incubation period, which implies that in their bodies the pathogen is present and replicating, but they are not capable of transmitting the disease to others. The population with acute infection I(t) is those who are currently infected with the disease, have developed a robust enough immunity to eradicate the disease from their bodies, and have the ability to transmit the disease to others. People with liver cirrhosis C(t) are asymptomatic infection individuals with end-stage liver disease who are capable of spreading the disease at any time t. Treated compartment T(t) includes those individuals who do not have enough immunity to clear the disease from their bodies; therefore, they are treated successfully, and they are capable of transmitting the disease to other people at any time t. The recovered class R(t) includes those individuals who gained enough immunity and eradicated the disease. The sub-model of alcohol consumption in our Hepatitis B virus transmission model is nonlinear, and this is because of the effects of behavioral saturation. This type of modeling has been followed in the recent works on hemodynamic modelling [24], where the blood flow through the arteries’ dynamics is described by the non-linear wave equations with the Bernoulli-type terms. These models make clear the fact that, in both fluid flow and behavioral epidemiology, nonlinearities play a critical role in relation to adequately describing biological phenomena.

We impose the following assumptions on the proposed model:

•   a1. The recruitment to the susceptible population is purely due to new births.

•   a2. Saturated incidence rate is considered.

•   a3. The disease transmission coefficient β(A)=β0+β1(A(t)−A0Amax) is dependent upon variations in alcohol intake among Hepatitis B infected individuals.

•   a4. The recovered population has permanent immunity.

•   a5. Individuals with successful vaccination goes to the recovered class.

•   a6. The vaccination may result in temporary immunity.

The model accounts for individual variations in alcohol metabolism through the dynamic transmission coefficient β(A), which adjusts based on alcohol intake levels. By incorporating a logistic growth term for alcohol consumption and a saturation incidence rate, the framework implicitly captures the nonlinear effects of alcohol metabolism on Hepatitis B progression. While the model does not explicitly include genetic or physiological factors, the alcohol-dependent transmission rate and saturation effects provide a proxy for metabolic variability, ensuring that the progression dynamics reflect the heterogeneity observed in real-world populations. Future refinements could integrate explicit metabolic parameters to further enhance the model’s precision.

The proposed model is given by the following system of differential equations, and Fig. 1 represents its schematic diagram, while parameters of the model are listed in Table 1 with their descriptions, values, and sources.

Figure 1: Flow diagram of the formulated model

The transmission function β(A) incorporates alcohol consumption due to its significant biological and behavioral impact on Hepatitis B dynamics, despite Hepatitis B being primarily blood-borne. Chronic alcohol consumption exacerbates liver damage, increases viral replication, and weakens immune responses, leading to higher viral loads in transmissible bodily fluids such as blood and semen. This biological mechanism is supported by studies such as Ganesan et al. [7], who demonstrated that alcohol accelerates Hepatitis B progression by promoting viral persistence and liver fibrosis. Furthermore, alcohol use is associated with riskier behaviors (e.g., unsafe injections or sexual practices), indirectly amplifying transmission rates. The function β(A)=β0+β1(A(t)−A0Amax) quantifies this relationship, where β0 represents the baseline transmission rate and β1 captures the incremental risk from alcohol. This formulation aligns with clinical evidence and ensures the model reflects the synergistic effects of alcohol and Hepatitis B on liver cirrhosis progression.

The model accounts for the threshold alcohol consumption required to significantly impact Hepatitis B transmission rates through the transmission coefficient β(A), defined as β(A)=β0+β1(A(t)−A0Amax). Here, β0 represents the baseline transmission rate without alcohol, while β1 scales the additional risk due to alcohol consumption. The logistic growth of alcohol consumption, governed by dAdt=r(A(t)−A0)(1−A(t)Amax), ensures that A(t) remains bounded between A0 and Amax. This formulation captures the saturation effect of alcohol on transmission rates, reflecting empirical observations that excessive alcohol intake accelerates Hepatitis B progression. The sensitivity analysis further confirms the pronounced impact of β1 on ℛ0, particularly when A(t) exceeds A0, aligning with clinical evidence on the synergistic effects of alcohol and Hepatitis B.

3  Qualitative Analysis

3.1 Existence and Uniqueness of the Solution

In this subsection, we prove that the solution to system (1) exists and is unique.

Theorem 1: The solution of the model (1) equations together with the initial conditions exists in R+7.

Proof: The RHS of model (1) can be written as follows,

F1(S,E,I,C,T,R,A)=b−β(A)S(t)I(t)1+α2I(t)−(ν+μ0)S(t),F2(S,E,I,C,T,R,A)=β(A)S(t)I(t)1+α2I(t)−(ψ+μ0)E(t),F3(S,E,I,C,T,R,A)=ψγ1E(t)−(ρ+τ1+μ0+μ1)I(t),F4(S,E,I,C,T,R,A)=ργ2I(t)−(τ2+γ3+μ0+μ2)C(t),F5(S,E,I,C,T,R,A)=τ1I(t)+τ2C(t)−(γ4+μ0)T(t),F6(S,E,I,C,T,R,A)=ψ(1−γ1)E(t)+ρ(1−γ2)I(t)+νS(t)+γ3C(t)+γ4T(t)−μ0R(t),F7(S,E,I,C,T,R,A)=r(A(t)−A0)(1−A(t)Amax).

Let Ω denote the region which is given by

Ω={(S(t),E(t),I(t),C(t),T(t),R(t),A(t))∈R+7:N(t)≤bμ0}

and let U⊆Ω represent the region |t−t0|≤δ, ||x−x0||≤ϵ, where x=(x1,x2,...,xn) and x0=(x10,x20,...xn0) also suppose that a(t,x) satisfies the Lipchitz condition:

||a(t,x1)−a(t,x2)||≤k||x1−x2||

whenever the pairs (t,x1), (t,x2) belongs to U⊆Ω where k is a positive constant, then there exist a positive constant δ≥0, such that there exists a unique and continuous vector solution x(t) of the system (1) in interval |t−t0|<δ. This condition is satisfied, if ∂Fi∂xj∀i,j are bounded in U⊆Ω, where x1=S, x2=E, x3=I, x4=C, x5=T, x6=R, x7=A.

For F1;

|∂F1∂S|=|−β(A)I(t)1+α2I(t)−(ν+μ0)|<∞,|∂F1∂E|=0<∞,|∂F1∂I|=|−β(A)S(t)(1+α2I(t))2|<∞,|∂F1∂C|=0<∞,|∂F1∂T|=0<∞,|∂F1∂R|=0<∞,|∂F1∂A|=0<∞

The Partial derivative exist, similarly for F2, F3, F4, F5, F6 and F7. This shows that all the partial derivatives ∂Fi∂xj ∀ i,j exists, and are continuous as well as bounded in U⊆Ω. This demonstrates that all of the partial derivatives ∂Fi∂xj ∀ i,j exist, are continuous and are bounded in U⊆Ω. Hence, by Lipchitz condition, the model (1) has a unique solution. ◻

3.2 Positivity of the Solution

In this sub-section, we establish positivity of the solutions of model (1) equations.

Theorem 2: For all given positive initial values, solutions S(t), E(t), I(t), C(t), T(t), R(t) and A(t) of system (1) are non-negative ∀ t>0.

Proof: Consider first equation of system (1), which implies that,

dS(t)dt+(β(A)I(t)1+α2I(t)+(ν+μ0))S(t)≥0,(2)

multiplying both sides of inequality (2) by integrating factor and integrating, gives:

S(t)≥C e−((ν+μ0)t+∫β(A)I(t)1+α2I(t)dt),(3)

where C is a constant of integration. By applying initial condition S(0)=S0 and solving gives C=S0 then substituting in inequality (3), the final solution yields,

S(t)≥S0 e−((ν+μ0)t+∫β(A)I(t)1+α2I(t)dt).(4)

In inequality (4), S0>0 and the exponentials are always non-negative. Hence S(t)≥0, ∀ t>0.

The positivity of the remaining state variables E(t), I(t), C(t), T(t), R(t) and A(t) for t>0 can be proved by similar approach. ◻

3.3 Feasible Region

In this subsection, we establish the boundedness of the state variables and define the feasible region.

Theorem 3 [27]: The solution for the system (1) is bounded and the closed set Ω is biologically feasible region of the system (1) such that Ω={(S,E,I,C,T,R)∈R+6:S>0,(E,I,C,T,R)≥0, N(t)≤bμ0}.

Proof: Since, the total population is represented by N(t), where:

N(t)=S(t)+E(t)+I(t)+C(t)+T(t)+R(t),(5)

differentiating Eq. (5) w.r.tt, substituting the values from system (1) and simplifying, we get:

dN(t)dt=b−μ0(S+E+I+C+T+R)−μ1I(t)−μ2C(t),(6)

now, substituting Eq. (5) into Eq. (6), implies that

dN(t)dt=b−μ0N−μ1I−μ2C,

this implies that,

dN(t)dt+μ0N(t)≤b,(7)

multiplying both sides of inequality (7) by the integrating factor and integrating, we obtain:

eμ0tN(t)≤bμ0eμ0t+C,(8)

where C is a constant of integration. Further, using the initial condition N(0)=N0, and solving for C gives C=N0−bμ0, then inequality (8) implies that:

N(t)≤bμ0+(N0−bμ0)e−μ0t.(9)

In inequality (9) when t→∞, then N(t)→bμ0 implying that 0≤N(t)≤bμ0. As a result, we establish that each solution of system (1) associated with an initial starting point x0∈R+6 remains in Ω. Hence, the closed region Ω is invariant positively at all times t. Therefore, it is sufficient to examine the dynamics of this model in the region Ω as it is both mathematically well-posed and epidemiologically significant. ◻

3.4 Equilibrium Points

In this subsection, we calculate two disease-free and disease-endemic equilibrium points when A(t)=A0 and A(t)=Amax.

3.4.1 Disease-Free Equilibrium Point

The disease-free equilibrium point is given by the set 𝒫DFE=[𝒫A0DFE,𝒫AmaxDFE] such that,

𝒫A0DFE=(S0,E0,I0,C0,T0,R0,A0)=(bν+μ0,0,0,0,0,νbμ0(ν+μ0),A0),

𝒫AmaxDFE=(S0,E0,I0,C0,T0,R0,A0)=(bν+μ0,0,0,0,0,νbμ0(ν+μ0),Amax).

3.4.2 Disease-Endemic Equilibrium Point

The disease-endemic equilibrium point is given by the set 𝒫DEE=[𝒫A0DEE,𝒫AmaxDEE] such that,

𝒫A0DEE={S∗=(ρ+τ1+μ0+μ1)(ψ+μ0)+bα2ψγ1β0bψγ1+α2ψγ1(ν+μ0)E∗=β0bψγ1−(ρ+τ1+μ0+μ1)(ψ+μ0)(ν+μ0)(ψ+μ0)ψγ1(β0+(ν+μ0)α2)I∗=β0bψγ1−(ρ+τ1+μ0+μ1)(ψ+μ0)(ν+μ0)(ρ+τ1+μ0+μ1)(ψ+μ0)(β0+(ν+μ0)α2)C∗=ργ2(β0bψγ1−(ρ+τ1+μ0+μ1)(ψ+μ0)(ν+μ0))(ρ+τ1+μ0+μ1)(ψ+μ0)(β0+(ν+μ0)α2)(τ2+γ3+μ0+μ2)T∗=(β0bψγ1−(ρ+τ1+μ0+μ1)(ν+μ0)(ψ+μ0))(τ1(τ2+γ3+μ0+μ2)−τ2ργ2)(ρ+τ1+μ0+μ1)(τ2+γ3+μ0+μ2)(ψ+μ0)(γ4+μ0)(β0+(ν+μ0)α2)R∗=ψ(1−γ1)E∗+ρ(1−γ2)I∗+νS∗+γ3C∗+γ4T∗μ0(10)

where β(A0)=β0, and

𝒫AmaxDEE={S∗=(ρ+τ1+μ0+μ1)(ψ+μ0)+bα2ψγ1(β0+zβ1)bψγ1+α2ψγ1(ν+μ0)E∗=(β0+zβ1)bψγ1−(ρ+τ1+μ0+μ1)(ψ+μ0)(ν+μ0)(ψ+μ0)ψγ1((β0+zβ1)+(ν+μ0)α2)I∗=(β0+zβ1)bψγ1−(ρ+τ1+μ0+μ1)(ψ+μ0)(ν+μ0)(ρ+τ1+μ0+μ1)(ψ+μ0)((β0+zβ1)+(ν+μ0)α2)C∗=ργ2((β0+zβ1)bψγ1−(ρ+τ1+μ0+μ1)(ψ+μ0)(ν+μ0))(ρ+τ1+μ0+μ1)(ψ+μ0)((β0+zβ1)+(ν+μ0)α2)(τ2+γ3+μ0+μ2)T∗=((β0+zβ1)bψγ1−(ρ+τ1+μ0+μ1)(ν+μ0)(ψ+μ0))(τ1(τ2+γ3+μ0+μ2)−τ2ργ2)(ρ+τ1+μ0+μ1)(τ2+γ3+μ0+μ2)(ψ+μ0)(γ4+μ0)((β0+zβ1)+(ν+μ0)α2)R∗=ψ(1−γ1)E∗+ρ(1−γ2)I∗+νS∗+γ3C∗+γ4T∗μ0(11)

where, β(Amax)=(β0+zβ1) and z=1−A0Amax.

3.5 The Basic Reproductive Number ℛ0

The fundamental reproductive number, ℛ0, represents the number of secondary infections caused by an infected individual in an entirely vulnerable population. To obtain the basic reproductive number ℛ0 for model (1), the work of Watmough and Driessche [28] is being followed; accordingly, the non-infected and infected classes must be separated. From model (1), we take the infected classes E(t),I(t),C(t) and T(t) and let 𝒳=(E(t),I(t),C(t),T(t)), we have,

d𝒹𝒳dt=ℱ−𝒱

the Jacobian matrices of ℱ and 𝒱 evaluated at the disease-free equilibrium point are,

ℱ∗=(0β(A)S000000000000000)and𝒱∗=(ψ+μ0000−ψγ1ρ+τ1+μ0+μ1000−ργ2τ2+γ3+μ0+μ200−τ1−τ2γ4+μ0),

the next-generation matrix is given by,

ℱ∗.𝒱∗−1=(β(A)ψγ1S0Z1Z2β(A)S0Z200000000000000),

whereZ1=ψ+μ0andZ2=ρ+τ1+μ0+μ1.

The basic reproductive number ℛ0 is determined as the maximum eigenvalue (the spectral radius, see [29]) of the matrix ℱ∗.𝒱∗−1. Hence,

ℛ0=β(A)bψγ1(ψ+μ0)(ν+μ0)(ρ+τ1+μ0+μ1).(12)

From Eq. (12), two different expressions are calculated for ℛ0. We put A(t)=A0 in Eq. (12) to compute the first basic reproduction number denoted by ℛA0. That is, when a patient with chronic hepatitis B infection consumes even a minimal amount of alcohol, thus,

ℛA0=β0bψγ1(ψ+μ0)(ν+μ0)(ρ+τ1+μ0+μ1),(13)

by substituting A(t)=Amax in Eq. (12), we compute the other basic reproductive number denoted by ℛAmax. That is, when a patient with chronic hepatitis B infection ingests the maximum amount of alcohol, thus,

ℛAmax=(β0+zβ1)bψγ1(ψ+μ0)(ν+μ0)(ρ+τ1+μ0+μ1),(14)

assuming z=1−A0Amax. We can write ℛ0=[ℛA0,ℛAmax] in compact form by observing that 0<ℛA0<ℛAmax.

3.6 3D-Type Simulation of ℛ0

Fig. 2 shows the effect of different model parameters on ℛA0 and ℛAmax, respectively.

Figure 2: 3D type simulation of ℛA0 and ℛAmax

3.7 Stability Analysis

In this subsection, we present both local and global asymptotic stabilities of the equilibria of system (1).

3.7.1 Local Stability

Theorem 4: The disease-free equilibrium (DFE) point 𝒫DFE=[𝒫A0DFE,𝒫AmaxDFE] is locally asymptotically stable (LAS) when ℛA0<ℛAmax<1, otherwise unstable.

Proof: The Jacobian matrix of system (1), at disease-free equilibrium point 𝒫DFE when A∗∈[A0,Amax] is represented by J(𝒫DFE) and given as:

J(𝒫DFE)=(−Z00−β(A∗)S000000−Z1β(A∗)S000000ψγ1−Z2000000ργ2−Z300000τ1τ2−Z400νZ5Z6γ3γ4−μ00000000−Z)

where,

Z0=ν+μ0,Z1=ψ+μ0,Z2=ρ+τ1+μ0+μ1,Z3=τ2+γ3+μ0+μ2,Z4=γ4+μ0,Z5=ψ(1−γ1),Z6=ρ(1−γ2),Z=r(2A∗−(A0+Amax))Amax.

The matrix J(𝒫DFE) can be written in terms of block matrices as,

J(𝒫DFE)=(𝒟1𝒪𝒟2𝒟3)

where, 𝒪3×4 is a void matrix and,

𝒟1=(−Z00−β∗(A∗)S00−Z1β∗(A∗)S00ψγ1−Z2),𝒟2=(00ργ200τ1νZ5Z6000),𝒟3=(−Z3000τ2−Z400γ3γ4−μ00000−Z),

hence, the matrix J(𝒫DFE) is a block diagonal matrix, and thus its eigenvalues are given by:

{Eigenvalues of 𝒟1}∪{Eigenvalues of 𝒟3}.

For matrix 𝒟3, the eigenvalues are obtained as, λ1,2,3,4=−Z3, −Z4, −μ0, −Z while, matrix 𝒟1 has eigenvalues with negative real parts iff its trace is negative and its determinant is positive [30]. Obviously,

Tr 𝒟1=−Z0−Z1−Z2⇒Tr 𝒟1<0.

Now, to verify that the determinant must be positive, we have

|𝒟1|=|−Z00−β∗(A∗)S00−Z1β∗(A∗)S00ψγ1−Z2|>0,

expanding by column 1 and simplifying yields,

ψγ1β(A∗)S0Z1Z2<1,ℛ0<1.

Thus, det(𝒟1)>0iffℛ0<1. Therefore, this implies that the disease-free equilibrium point 𝒫DFE is locally asymptotically stable when ℛA0<ℛAmax<1, otherwise unstable. ◻

Theorem 5: The disease-endemic equilibrium point 𝒫DEE=[𝒫A0DEE,𝒫AmaxDEE] is locally asymptotically stable when ℛAmax>ℛA0>1, otherwise unstable.

Proof: The Jacobian matrix of system (1), at disease-endemic equilibrium point 𝒫DEE is represented by J(𝒫DEE) and given as,

J(𝒫DEE)=(−Q00−Q10000Q2−Z1Q100000ψγ1−Z2000000ργ2−Z300000τ1τ2−Z400νZ5Z6γ3γ4−μ00000000−Z)

where,

Q0=β(A)I∗1+α2I∗+(ν+μ0),Q1=β(A)S∗(1+α2I∗)2,Q2=β(A)I∗1+α2I∗Z1=ψ+μ0,Z2=ρ+τ1+μ0+μ1,Z3=τ2+γ3+μ0+μ2,Z4=γ4+μ0,Z=r(2A∗−(A0+Amax))Amax,

since, the eigenvalues of the Jacobian matrix J(𝒫DEE) are obtained by |J(𝒫DEE)−λI|=0. Therefore, the four eigenvalues are λ1,2,3,4=−Z, −μ0, −Z4, −Z3 and the remaining are eigenvalues of the following matrix:

ℬ=(−Q00−Q1Q2−Z1Q10ψγ1−Z2),

using the characteristics equation |ℬ−λI|=0, solving we get the following characteristics polynomial:

λ3+(Q0+Z1+Z2)λ2+(Q0Z1+Q0Z2+Z1Z2−Q1ψγ1)λ+Q0Z1Z2+Q1ψγ1(Q2−Q0)=0,

this implies that, f0λ3+f1λ2+f2λ+f3=0, where, f0=1,f1=Q0+Z1+Z2,f2=Q0Z1+Q0Z2+Z1Z2−Q1ψγ1,f3=Q0Z1Z2+Q1ψγ1(Q2−Q0)

Now, we need to verify the following two conditions:

(a) f1,f2,f3>0 (b) f1f2−f0f3>0

Condition (a) is satisfied as f1>0 and f2,f3>0 iff Q0(Z1+Z2)+Z1Z2>Q1ψγ1 Q1>Q2, respectively. Also condition (b) holds if and only if f1f2>f0f3.

Thus, the Routh-Hurwitz criterion [31] as well as conditions (a) and (b) imply that the characteristic equation has all roots with negative real parts. Thus λ5,6,7<0, which guarantees the local stability of disease-endemic equilibrium point 𝒫DEE=[𝒫A0DEE,𝒫AmaxDEE]. ◻

3.7.2 Global Stability

Theorem 6 [32]: The disease-free equilibrium point 𝒫DFE=[𝒫A0DFE,𝒫AmaxDFE], of system (1) is globally asymptotically stable in Ω if ℛA0<ℛAmax<1.

Proof: To examine global asymptotic stability of the disease-free equilibrium point of system (1), let us first construct a suitable candidate Lyapunov function 𝒢:Δ⊂R+6→R where Δ=(S,E,I,C) such that,

𝒢(S,E,I,C)=12[(S(t)−S0)+E(t)+I(t)+C(t)]2(15)

clearly, 𝒢:Δ⊂R+6→R is strictly positive definite and equal to zero at the disease-free equilibrium point 𝒫DFE=[𝒫A0DFE,𝒫AmaxDFE]. Differentiating Eq. (15) with respect to time t, we have,

d𝒢dt=[(S(t)−S0)+E(t)+I(t)+C(t)][dSdt+dEdt+dIdt+dCdt](16)

substitute values from system (1) to Eq. (16), and simplifying we get,

d𝒢dt=[(S−S0)+E+I+C]×[b−(ν+μ0)S−μ0E−ψE+ψγ1E−τ1I−μ1I−μ0I−ρI+ργ2I−τ2C−γ3C−μ2C−μ0C],d𝒢dt=−[(S−S0)+E+I+C]×[(ν+μ0)(S+S0)+μ0(E+I+C+T)+(1−γ1)ψE+(1−γ2)ρI+(τ1+μ1)I+(τ2+γ3+μ2)C].

Thus, d𝒢dt<0ifℛA0<ℛAmax<1, and d𝒢dt=0iff S(t)=S0 and E(t)=I(t)=C(t)=0. The condition ℛAmax<1 is the cornerstone for global eradication. It ensures that even under the worst-case scenario of maximum alcohol consumption (Amax), which amplifies transmission, the infection cannot sustain itself. When this condition holds, the disease-free state acts as a global attractor, meaning the population will inevitably recover from any initial outbreak and converge to a disease-free state over time. This mathematical guarantee is robust within the model’s framework, though its real-world feasibility depends on the system parameters remaining within biologically plausible ranges, as utilized in our numerical simulations.

Hence, the only largest compact positively invariant set is the singleton set 𝒫DFE in {(S,E,I,C,T,R)∈Ω:(d𝒢/dt)=0}. ◻