Stochastic Optimal Control Analysis for HBV Epidemic Model with Vaccination

Abstract

In this study, we explore the concept of symmetry as it applies to the dynamics of the Hepatitis B Virus (HBV) epidemic model. By incorporating symmetric principles in the stochastic model, we ensure that the control strategies derived are not only effective but also consistent across varying conditions, and ensure the reliability of our predictions. This paper presents a stochastic optimal control analysis of an HBV epidemic model, incorporating vaccination as a pivotal control measure. We formulate a stochastic model to capture the complex dynamics of HBV transmission and its progression to acute and chronic stages. By leveraging stochastic differential equations, we examine the model’s stationary distribution and asymptotic behavior, elucidating the impact of random perturbations on disease dynamics. Optimal control theory is employed to derive control strategies aimed at minimizing the disease burden and vaccination costs. Through rigorous numerical simulations using the fourth-order Runge–Kutta method, we demonstrate the efficacy of the proposed control measures. Our findings highlight the critical role of vaccination in controlling HBV spread and provide insights into the optimization of vaccination strategies under stochastic conditions. The symmetry within the proposed model equations allows for a balanced approach to analyzing both acute and chronic stages of HBV.

1. Research Background

Viruses infecting the liver cause hepatitis B. Hepatitis B is a significant global health concern, with millions of people affected worldwide. Understanding this disease is crucial for prevention, diagnosis, and treatment. Hepatitis B is a highly contagious virus transmitted through various means. These include sharing needles or syringes, receiving contaminated blood transfusions or organ transplants, or being stuck with a contaminated needle. Unprotected sexual intercourse with an infected partner can also result in transmission. The virus can also be spread through objects like razors or toothbrushes contaminated with infected blood and through an infected mother’s breast milk. Healthcare workers are at risk of exposure to infected blood, making it essential to take necessary precautions when dealing with infected patients [1,2,3,4]. It is possible for individuals with hepatitis B to not show any symptoms, especially during the initial stages of the infection. However, when symptoms do appear, they may include jaundice, loss of appetite, fatigue, abdominal pain, nausea, joint pain, and dark urine. Hepatitis B infections are classified as either acute or chronic. Acute hepatitis B typically lasts a few months, and the body’s immune system can usually clear the virus within this time. However, if the virus persists in the body for over six months, it is considered a chronic infection [5]. The result can be severe liver damage and an increased risk of liver cancer [6,7]. Therefore, it is essential to diagnose and treat Hepatitis B as early as possible to avoid any potential complications. Blood tests can determine if hepatitis B is acute or chronic, detect viral antigens and antibodies, and assess liver function. Hepatitis B can be prevented through vaccination. Practicing safe sex, avoiding the sharing of needles or items that may have blood on them, and maintaining proper infection control measures in healthcare facilities are essential to preventing the spread of the infection. The symptoms of acute hepatitis B can be managed with supportive care rather than specific treatment. On the other hand, chronic hepatitis B is managed with antiviral medications to suppress viral replication and lower liver inflammation. In severe cases, a liver transplant may be necessary for individuals with significant liver damage. Over 250 million people are chronically infected with hepatitis B worldwide [8]. Additionally, it leads to liver-related deaths and morbidities.

Mathematical modeling is one of the most essential tools for understanding disease transmission dynamics [9,10]. Mathematical models have been instrumental in the study of HBV and have helped researchers to identify crucial mechanisms of the virus, such as its transmission, effects on the body, and available treatments. These models have also provided valuable insights for developing effective vaccines [11], which have saved countless lives in the past decades and will continue to be invaluable in the fight against HBV. By enabling scientists to understand the intricate workings of the virus more deeply than ever, these models have allowed them to identify critical factors such as transmission and prevention methods. This breakthrough has revolutionized the fight against the virus, providing scientists the means to combat it more effectively. It has enabled governments to make more informed decisions about responding to the virus. Mathematicians have developed a number of mathematical models to control the spread of infectious diseases. As a result of these models, public health policies have been informed, resources have been targeted to those who are most at risk, and predictions have been made about the spread of the disease. They have also been used to assess the efficacy of different interventions, such as contact tracing and social distancing [12,13]. Ghulam et al. [14] delved into a stochastic epidemic model for hepatitis B featuring saturated incidence rates. Amir Khan et al. [15] explored the impact of environmental fluctuations on transmission rates using a stochastic epidemic model. Murad et al. [16] introduced a stochastic model for the transmission of the Hepatitis B virus (HBV), incorporating environmental noise to elucidate the dynamics. Tahir Khan et al. [17] developed stochastic models to investigate the dynamics of the hepatitis B epidemic, specifically examining population changes and their long-term impact on the epidemic’s behavior. Jiying Mawe and Shasha Ma et al. [18] investigated a transmission model for stochastic HBV with media coverage featuring saturated incidence rates. Mathematical simulation models are widely utilized to study and control infectious diseases, particularly in the context of preventing the spread of contagious diseases [19,20,21,22,23].

In this study, the authors created and analyzed a probability-based model that determines stability. A control theory was also developed for the stochastic and deterministic models as part of the study. In the control aspect, the main focus was determining the control measure of the infection rate to reduce the number of infected individuals with time. Stochastic models were optimized to minimize their expected values. A deterministic model was approached using the Pontryagin principle, while stochastic models were controlled using the Hamilton–Jacobi–Bellman equation.

The subsequent sections of this paper are organized in the following way. Section 2 presents the model that governs the dynamics of HBV. Section 3 presents the uniqueness and existence of a non-negative solution globally. Section 4 presents an in-depth investigation of the disease’s extinction, from which necessary criteria are determined. The suggested model’s ergodicity and the possibility of a specific stationary distribution are examined in Section 5. Section 6 discusses the optimal control for both deterministic and stochastic models. Deterministic model optimal control is based on Pontryagin’s maximum principle, while for stochastic model is based on the Hamilton–Jacobi–Bellman equation. Both methods aim to maximize the expected value of the system’s cost function. In Section 7, the simulation results are displayed. We summarize our results in Section 8 at the conclusion.

2. Model Formulation

An HBV stochastic model will be formulated in this section with different assumptions. The whole population is represented by $N\left(t\right)$ and then divided into six different classes: $S\left(t\right)$ susceptible, $V\left(t\right)$ vaccinated, $A\left(t\right)$ acutely infected, $Z\left(t\right)$ chronically infected, $H\left(t\right)$ hospitalized, and $R\left(t\right)$ recovered. A community’s population at any given t time is represented by $N\left(t\right)$. In other words, $N\left(t\right)=S\left(t\right)+V\left(t\right)+A\left(t\right)+Z\left(t\right)+H\left(t\right)+R\left(t\right)$. A model’s long-term behaviour may be affected by the changing population environment, which is important to note. The following hypothesis is therefore applied to the model.

L₁: 

We assume that both the parameter and system state are non-negative.

L₂: 

Acute duration is short, and if the treatment fails during this period, the individual will be considered in a chronic class.

L₃: 

Our model will include environmental noise by considering the function $B_i\left(t\right)$ for $i=\{1,2,3,4,5,6\}$ as the Brownian standard motion, and ${\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,$ and ${\delta }_6>0$ represent the intensities of white noises. Furthermore, Brownian motion satisfies $B_i\left(0\right)=0$ for all i.

L₄: 

Once an individual recovers from the disease, they will obtain permanent immunity.

The stated assumptions $(L_1-L_4)$ leads to the mathematical model described by stochastic Equation (1):

$$\begin{array}{cc}\hfill dS\left(t\right)& =[\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi )S\left(t\right)]dt+{\delta }_{1}S\left(t\right)d{B}_1\left(t\right),\hfill \\ \hfill dV\left(t\right)& =[\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right)]dt+{\delta }_{2}V\left(t\right)d{B}_2\left(t\right),\hfill \\ \hfill dA\left(t\right)& =[{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )A\left(t\right)]dt+{\delta }_{3}A\left(t\right)d{B}_3\left(t\right),\hfill \\ \hfill dZ\left(t\right)& =[{\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi )Z\left(t\right)]dt+{\delta }_{4}Z\left(t\right)d{B}_4\left(t\right),\hfill \\ \hfill dH\left(t\right)& =[{\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right)]dt+{\delta }_{5}H\left(t\right)d{B}_5\left(t\right),\hfill \\ \hfill dR\left(t\right)& =[{\gamma }_{1}Z\left(t\right)+{\alpha }_{2}A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)]dt+{\delta }_{6}R\left(t\right)d{B}_6\left(t\right).\hfill \end{array}$$

(1)

where $(B_1,B_2,B_3,B_4,B_5,B_6,)\left(t\right)$ are independent standard Brownian motions, and classic Gaussian white noise intensities are $({\delta }_1$, ${\delta }_2$, ${\delta }_3$, ${\delta }_4$, ${\delta }_5$, ${\delta }_6)$ accordingly. The contact between the surrounding environment and an individual in this model is shown by ${\delta }_{1}S\left(t\right)B_1$, ${\delta }_{2}V\left(t\right)B_2$, ${\delta }_{3}A\left(t\right)B_3$, ${\delta }_{4}Z\left(t\right)B_4$, ${\delta }_{5}H\left(t\right)B_5$, and ${\delta }_{6}R\left(t\right)B_6$.

In model (1), all parameters are positive, reflecting various aspects of HBV dynamics, as detailed. The parameters include $\Pi$, the recruitment rate of new individuals; ${\beta }_1$ and ${\beta }_2$, which describe the infection rates from susceptible and vaccinated individuals, respectively; $\eta$, the rate of waning vaccine-induced immunity; and $\chi$, the frequency of natural deaths. Other parameters include $\gamma$, the contact rate between individuals with weakened immune systems and those who are sick; ${\gamma }_1$, the recovery rate for chronically ill individuals; ${\gamma }_2$, the chronic HBV mortality rate; ${\alpha }_1$ and ${\alpha }_2$, which describe the transition rate from acute to chronic infection and recovery rates, respectively; and $\psi$, the proportion of newborns without effective vaccination. Additionally, ${\psi }_1$, ${\psi }_2$, ${\psi }_3$, and ${\psi }_4$ represent hospitalization rates for acute and chronic patients, recovery rates from hospitalization, and disease-induced mortality rates, respectively. These parameters collectively shape the model’s depiction of HBV transmission and progression within the population.

Model (1) has a deterministic form:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi )S\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dV}{dt}}& =\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dA}{dt}}& ={\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )A\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dZ}{dt}}& ={\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi )Z\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dH}{dt}}& ={\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& ={\gamma }_{1}Z\left(t\right)+{\alpha }_{2}A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right).\hfill \end{array}$$

(2)

According to [24,25], deterministic model (2) contains a basic reproduction number. According to a population of susceptible individuals, the basic reproduction number is the number of infected people infected by infected individuals. This number measures the disease’s spread potential. It is a critical indicator of the severity of an epidemic.

$$R_0^D={\displaystyle \frac{{\beta }_1+{\beta }_2}{({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )}}.$$

(3)

Both models (1) and (2) are explained in graphical form in Figure 1a,b. Figure 1a,b clearly show how the HBV model works, with the infected host continuously producing and transmitting the virus to others.

3. Existence and Uniqueness in the Realm of Positive Solutions

According to a biological study of models (1) and (2), their solution should be positive globally. The initial step is to prove that this property is satisfied by the solution of model (1).

Theorem 1.

Model (1) has a unique solution $\left(S\right(t),V(t),A(t),Z(t),H(t),R(t\left)\right)$ regardless of the initial value of $t\ge 0$ $\left(S\right(0),V\left(0\right),A\left(0\right),Z\left(0\right),H\left(0\right),R\left(0\right))\in R_{+}^6$, and the solution remains in $R_{+}^6$ with probability one, as, $S\left(t\right),V\left(t\right),A\left(t\right),Z\left(t\right),H\left(t\right),R\left(t\right)\in R_{+}^6$ fol all $t\ge 0$ a.s.

Proof. 

Evidently, the coefficients governing model (1) exhibit local Lipschitz continuity irrespective of the chosen initial values $\left(S\left(0\right),V\left(0\right),A\left(0\right),Z\left(0\right),H\left(0\right),R\left(0\right)\right)\in {\mathbb{R}}_{+}^6$. This ensures the existence of a unique local solution $\left(S\right(t),V(t),A(t),Z(t),H(t),R(t\left)\right)$ over the interval $t\in [0,{\tau }_e)$, where ${\tau }_e$ represents the explosion time (refer to [26] for a comprehensive explanation). To establish global validity, it becomes imperative to demonstrate that ${\tau }_e$ is almost surely infinite. Select $k_0\ge 0$ to be sufficiently large such that $\left(S\right(0),V(0),A(0),Z(0),H(0),R(0\left)\right)$ all reside within the interval $[{\displaystyle \frac{1}{k_0}},k_0]$. For every integer $k\ge k_0$, we define the subsequent stopping time.

$${\tau }_k=inf\left\{t\in [0,{\tau }_e):minS\left(t\right),V\left(t\right),A\left(t\right),Z\left(t\right),H\left(t\right),R\left(t\right)\}\le {\displaystyle \frac{1}{k}}\mathrm{or}max\{S\left(t\right),V\left(t\right),A\left(t\right),Z\left(t\right),H\left(t\right),R\left(t\right)\ge k\}\right\}.$$

(4)

The set $inf⌀$ throughout this paper equals ∞ (since we normally denote empty sets as ⌀). Based on this definition, ${\tau }_k$ increases as $k\to \infty$. If ${\tau }_{\infty }=li{m}_{k\to \infty }{\tau }_k$, then ${\tau }_{\infty }$ ≤ ${\tau }_e$ a.s., that is, if ${\tau }_{\infty }$ = ∞ a.s. A pair of constants $T>0$ also exists, $ϵ\in$$(0,1)$, if this assertion is false; thus,

$$ϵ<P\{T\ge {\tau }_{\infty }\}.$$

(5)

with an integer $k_1\ge$ $k_0$. Therefore,

$$P\{T\ge {\tau }_k\}\ge ϵ,\forall k_1\le k.$$

(6)

Now, a $C^2$-function $V:R_{+}^6\to R_{+}$, where $R_{+}=x\in R:x\ge 0$, is as follows:

$$\begin{array}{cc}\hfill V(S,V,A,Z,H,R)& =S+V+A+Z+H+R-6-(lnS+lnV+lnA+lnZ+lnH+lnR).\hfill \end{array}$$

(7)

By applying Ito’s formula, we obtain the following result:

$$\begin{array}{cc}\hfill V(S,V,A,Z,H,R)& =LV(S,V,A,Z,H,R)+{\delta }_1\left(S-1\right)d{B}_1\left(t\right)+{\delta }_2(V-1)d{B}_2\left(t\right)+{\delta }_3(A-1)d{B}_3\left(t\right)\hfill \\ & +{\delta }_4(Z-1)d{B}_4\left(t\right)+{\delta }_5(H-1)d{B}_5\left(t\right)+{\delta }_6(R-1)d{B}_6\left(t\right).\hfill \end{array}$$

(8)

where

$$\begin{array}{cc}\hfill LV& =\left(1-{\displaystyle \frac{1}{S}}\right)\left(\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V-(\gamma +\chi )S\right)+{\displaystyle \frac{{\delta }_1^2}{2}}\hfill \\ \hfill & +\left(1-{\displaystyle \frac{1}{V}}\right)\left(\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right)\right)+{\displaystyle \frac{{\delta }_2^2}{2}}\hfill \\ \hfill & +\left(1-{\displaystyle \frac{1}{A}}\right)\left({\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )A\left(t\right)\right)+{\displaystyle \frac{{\delta }_3^2}{2}}\hfill \\ \hfill & +\left(1-{\displaystyle \frac{1}{Z}}\right)\left({\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi )Z\left(t\right)\right)+{\displaystyle \frac{{\delta }_4^2}{2}}\hfill \\ \hfill & +\left(1-{\displaystyle \frac{1}{H}}\right)\left({\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right)\right)+{\displaystyle \frac{{\delta }_5^2}{2}}\hfill \\ \hfill & +\left(1-{\displaystyle \frac{1}{R}}\right)\left({\gamma }_{1}Z\left(t\right)+{\alpha }_{2}A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)\right)+{\displaystyle \frac{{\delta }_6^2}{2}}.\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill LV& =\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi )S\left(t\right)-{\displaystyle \frac{\psi \Pi }{S\left(t\right)}}+{\displaystyle \frac{{\beta }_{1}A\left(t\right)}{N}}-{\displaystyle \frac{\eta V\left(t\right)}{S\left(t\right)}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}\hfill \\ \hfill & +\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{V}}-(\eta +\chi )V\left(t\right)-{\displaystyle \frac{\Pi (1-\psi )}{V\left(t\right)}}-{\displaystyle \frac{\gamma S\left(t\right)}{V\left(t\right)}}+{\displaystyle \frac{{\beta }_{2}A\left(t\right)}{N}}+\eta +\chi \hfill \\ \hfill & +{\displaystyle \frac{{\delta }_2^2}{2}}+{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )A\left(t\right)-{\displaystyle \frac{{\beta }_{1}S\left(t\right)}{N}}-{\displaystyle \frac{{\beta }_{2}V\left(t\right)}{N}}+{\alpha }_1+{\alpha }_2\hfill \\ \hfill & +{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}}+{\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi )Z\left(t\right)-{\displaystyle \frac{{\alpha }_{1}A\left(t\right)}{Z\left(t\right)}}+{\gamma }_1+{\psi }_2+{\gamma }_2+\chi +{\displaystyle \frac{{\delta }_4^2}{2}}\hfill \\ \hfill & +{\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right)-{\displaystyle \frac{{\psi }_{1}A\left(t\right)}{H}}-{\displaystyle \frac{{\psi }_{2}Z\left(t\right)}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+{\gamma }_{1}Z\left(t\right)\hfill \\ \hfill & ++{\alpha }_{2}A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)-{\displaystyle \frac{{\gamma }_{1}Z\left(t\right)}{R}}-{\displaystyle \frac{{\alpha }_{2}A\left(t\right)}{R}}-{\displaystyle \frac{{\psi }_{3}H\left(t\right)}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}.\hfill \\ \hfill & =\psi \Pi +\eta V\left(t\right)+{\displaystyle \frac{{\beta }_{1}A\left(t\right)}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\Pi +\gamma S\left(t\right)+{\displaystyle \frac{{\beta }_{2}A\left(t\right)}{N}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{2}}+{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}\hfill \\ \hfill & +{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}+{\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}}+{\alpha }_{1}A\left(t\right)+{\gamma }_1+{\psi }_2+{\gamma }_2+\chi +{\displaystyle \frac{{\delta }_4^2}{2}}+{\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)\hfill \\ \hfill & +{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+{\gamma }_{1}Z\left(t\right)+{\alpha }_{2}A\left(t\right)+{\psi }_{3}H\left(t\right)+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}-{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}-\gamma S\left(t\right)-\chi S\left(t\right)-{\displaystyle \frac{\psi \Pi }{S\left(t\right)}}\hfill \\ \hfill & -{\displaystyle \frac{\eta V\left(t\right)}{S\left(t\right)}}-\Pi \psi -{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{V}}-\eta V\left(t\right)-\chi V\left(t\right)-{\displaystyle \frac{\Pi (1-\psi )}{V\left(t\right)}}-{\displaystyle \frac{\gamma S\left(t\right)}{V\left(t\right)}}-{\alpha }_{1}A\left(t\right)-{\alpha }_{2}A\left(t\right)-{\psi }_{1}A\left(t\right)\hfill \\ \hfill & -\chi A\left(t\right)-{\displaystyle \frac{{\beta }_{1}S\left(t\right)}{N}}-{\displaystyle \frac{{\beta }_{2}V\left(t\right)}{N}}-{\gamma }_{1}Z\left(t\right)-{\psi }_{2}Z\left(t\right)-{\gamma }_{2}Z\left(t\right)-\chi Z\left(t\right)-{\displaystyle \frac{{\alpha }_{1}A\left(t\right)}{Z\left(t\right)}}-{\psi }_{3}H\left(t\right)-{\psi }_{4}H\left(t\right)\hfill \\ \hfill & -\chi H\left(t\right)-{\displaystyle \frac{{\psi }_{1}A\left(t\right)}{H}}-{\displaystyle \frac{{\psi }_{2}Z\left(t\right)}{H}}-\chi R\left(t\right)-{\displaystyle \frac{{\gamma }_{1}Z\left(t\right)}{R}}-{\displaystyle \frac{{\alpha }_{2}A\left(t\right)}{R}}-{\displaystyle \frac{{\psi }_{3}H\left(t\right)}{R}}.\hfill \\ \hfill & \hfill \le \Pi +{\beta }_1+{\beta }_2+\gamma +6\chi +\eta +{\alpha }_1+{\alpha }_2+{\psi }_1+{\psi }_2+{\psi }_3+{\psi }_4+{\gamma }_1+{\gamma }_2+{\displaystyle \frac{{\delta }_1^2+{\delta }_2^2+{\delta }_3^2+{\delta }_4^2+{\delta }_5^2+{\delta }_6^2}{2}}:=K.\end{array}$$

(10)

In view of the fact that K is a positive independent constant of $S,V,A,Z,H,R$ and t, we come up with the following:

$$\begin{array}{cc}\hfill \mathbb{V}(S,V,A,Z,H,R)& \le Kdt+{\delta }_1\left(S-1\right)d{B}_1\left(t\right)+{\delta }_2(V-1)d{B}_2\left(t\right)+{\delta }_3(A-1)d{B}_3\left(t\right)\hfill \\ & +{\delta }_4(Z-1)d{B}_4\left(t\right)+{\delta }_5(H-1)d{B}_5\left(t\right)+{\delta }_6(R-1)d{B}_6\left(t\right).\hfill \end{array}$$

(11)

By integrating both sides of Equation (11) and taking expectations from 0 to $T\wedge {\tau }_k$, we obtain

$$\begin{array}{cc}\hfill & \mathbb{EV}(S(T\wedge {\tau }_k),V(T\wedge {\tau }_k),A(T\wedge {\tau }_k),Z(T\wedge {\tau }_k),H(T\wedge {\tau }_k),R(T\wedge {\tau }_k))\hfill \\ & \le V\left(S\right(0),V(0),A(0),Z(0),H(0),R(0\left)\right)+KT<\infty .\hfill \end{array}$$

(12)

For each $K\ge k_1$, consider the event ${\mathsf{\Omega }}_k={\tau }_k\le t$ as defined in Equation (6). We have $P\left({\mathsf{\Omega }}_k\right)\ge \infty$. Upon observing each $\omega \in {\mathsf{\Omega }}_k$, at least one of the variables $S({\tau }_k,\omega )$, $V({\tau }_k,\omega )$, $A({\tau }_k,\omega )$, $Z({\tau }_k,\omega )$, $H({\tau }_k,\omega )$, and $R({\tau }_k,\omega )$ takes a value equal to either k or ${\displaystyle \frac{1}{k}}$. Consequently, the values of $S({\tau }_k,\omega )$, $V({\tau }_k,\omega )$, $A({\tau }_k,\omega )$, $Z({\tau }_k,\omega )$, $H({\tau }_k,\omega )$, and $R({\tau }_k,\omega )$ are at least {$k-1-logk$ or ${\displaystyle \frac{1}{k}}-1-log{\displaystyle \frac{1}{k}}$=${\displaystyle \frac{1}{k}}-1+logk$}

Following this,

$$\begin{array}{c}\hfill \mathbb{V}(S({\tau }_k,\omega ),V({\tau }_k,\omega ),A({\tau }_k,\omega ),Z({\tau }_k,\omega ),H({\tau }_k,\omega ),R({\tau }_k,\omega ),)+KT\ge (k-1-lnk)\wedge ({\displaystyle \frac{1}{k}}-1+lnk),\end{array}$$

(13)

The minimum of a and b is represented by $a\wedge b$, where $a\wedge b$ donates the minimum of a and b. Regarding (12) and (13), as of now,

$$\begin{array}{cc}\hfill & \mathbb{V}\left(S\right(0),V(0),A(0),Z(0),H(0),R(0\left)\right)+KT\ge \hfill \\ & \mathbb{E}\left[I_{\mathsf{\Omega }k}\mathbb{V}(S({\tau }_k,\omega ),V({\tau }_k,\omega ),A({\tau }_k,\omega ),Z({\tau }_k,\omega ),H({\tau }_k,\omega ),R({\tau }_k,\omega ))\right]\ge \hfill \\ & \in [(k-1-lnk)\wedge ({\displaystyle \frac{1}{k}}-1+lnk)].\hfill \end{array}$$

(14)

The indicator function for ${\mathsf{\Omega }}_k$ is $I_{\mathsf{\Omega }k}$. Suppose that $k\to$∞ leads to concentration $\infty >\mathbb{V}\left(S\right(0),V(0),A(0),Z(0),H(0),R(0\left)\right)+KT=\infty$ such that we have ${\tau }_{\infty }$ = ∞ a.s.

At this juncture, the demonstration of Theorem 1 reaches its culmination. □

4. Stochastic Analysis

In epidemiology, predicting an infectious disease’s long-term effects is one of the most important concepts. It introduces stochastic Lyapunov functions and delineates criteria for disease extinction within the framework of model (1). It may be possible to eliminate the infected compartment at the disease-free equilibrium of (1). If the disease-free equilibrium achieves global asymptotic stability, stochastic model (1) guarantees the eradication of the infection. Stochastic models are very interesting in extinction theory.

The Extinction of the Disease

This part proves a few lemmas and theorems and shows how the HBV extinction model works. Next, let us take a look at

$$〈Y\left(t\right)〉={\displaystyle \frac{1}{t}}{\int }_0^{t}y\left(x\right)dx.$$

(15)

Then, we have the following lemmas.

Lemma 1

([27,28] (Strong Law of Large Number)). Suppose that $Y=Y_{t\ge 0}$ is a local real-value continuous martingale that becomes zero at $t=0$; then,

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim}{〈Y,Y〉}_t=\infty ,a.s.\Rightarrow \underset{t\to \infty }{lim}{\displaystyle \frac{Y_t}{{〈Y,Y〉}_t}}=0,a.s.,also\hfill \\ \hfill & \underset{t\to \infty }{lim}sup{\displaystyle \frac{{〈Y,Y〉}_t}{t}}<0,a.s.\Rightarrow \underset{t\to \infty }{lim}{\displaystyle \frac{Y_t}{t}}=0,a.s.\hfill \end{array}$$

(16)

Lemma 2.

The following characteristics will be obtained for any define initial point $(S_0,V_0,A_0,Z_0,H_0,R_0)\in R_{+}^6,$ for the solution $\left(S\right(t),V(t),A(t),Z(t),H(t\left)R\right(t))$ of model (1):

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim}{\displaystyle \frac{S\left(t\right)}{t}}=0,\underset{t\to \infty }{lim}{\displaystyle \frac{V\left(t\right)}{t}}=0,\underset{t\to \infty }{lim}{\displaystyle \frac{A\left(t\right)}{t}}=0,\hfill \\ & \underset{t\to \infty }{lim}{\displaystyle \frac{Z\left(t\right)}{t}}=0,\underset{t\to \infty }{lim}{\displaystyle \frac{H\left(t\right)}{t}}=0,\underset{t\to \infty }{lim}{\displaystyle \frac{R\left(t\right)}{t}}=0,a.s.\hfill \end{array}$$

(17)

Moreover, when $\chi >{\displaystyle \frac{1}{2}}({\delta }_1^2\vee {\delta }_2^2\vee {\delta }_3^2\vee {\delta }_4^2\vee {\delta }_5^2\vee {\delta }_6^2)$ holds,

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim}{\displaystyle \frac{1}{t}}{\int }_0^{t}S\left(r\right)d{B}_1\left(r\right)=0,\underset{t\to \infty }{lim}{\displaystyle \frac{1}{t}}{\int }_0^{t}V\left(r\right)d{B}_2\left(r\right)=0,\underset{t\to \infty }{lim}{\displaystyle \frac{1}{t}}{\int }_0^{t}A\left(r\right)d{B}_3\left(r\right)=0,\hfill \\ & \underset{t\to \infty }{lim}{\displaystyle \frac{1}{t}}{\int }_0^{t}Z\left(r\right)d{B}_4\left(r\right)=0,\underset{t\to \infty }{lim}{\displaystyle \frac{1}{t}}{\int }_0^{t}H\left(r\right)d{B}_5\left(r\right)=0,\underset{t\to \infty }{lim}{\displaystyle \frac{1}{t}}{\int }_0^{t}R\left(r\right)d{B}_6\left(r\right)=0,a.s.\hfill \end{array}$$

(18)

Proof. 

In this case, Lemma 2 can be proved in the same manner as Lemma 4.1 in [29], so the proof is omitted.

Furthermore, upon fulfilling the specified conditions, we can assert that model (1) attains a stochastic disease-free state:

• $\left({\mathsf{ℂ}}_1\right):{\mathsf{R}}_0^E<1.$
• $\left({\mathsf{ℂ}}_2\right):\chi >{\displaystyle \frac{1}{2}}({\delta }_1^2\vee {\delta }_2^2\vee {\delta }_3^2\vee {\delta }_4^2\vee {\delta }_5^2\vee {\delta }_6^2)$.

where

$${\mathsf{R}}_0^E={\displaystyle \frac{{\beta }_1+{\beta }_2}{({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +\frac{{\delta }_3^2}{2})}}.$$

(19)

□

Theorem 2.

If the solution of model (1) satisfies both $C_1$ and $C_2$ with the initial conditions $(S_0,V_0,A_0,Z_0,H_0,R_0)$∈$R^{6+}$ for $\left(\right(S\left(t\right),V\left(t\right),A\left(t\right),Z\left(t\right),R\left(t\right))$, then it possesses the following properties:

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim\; sup}〈S\left(t\right)〉={\displaystyle \frac{\Pi (\eta +\chi \psi )}{\chi (\gamma +\eta +\chi )}},\underset{t\to \infty }{lim\; sup}〈V\left(t\right)〉={\displaystyle \frac{\chi \Pi (1-\psi )+\gamma \psi }{\chi (\gamma +\eta +\chi )}}\underset{t\to +\infty }{lim\; sup}{\displaystyle \frac{〈A\left(t\right)〉}{t}}=0,\hfill \\ & \underset{t\to +\infty }{lim\; sup}{\displaystyle \frac{〈Z\left(t\right)〉}{t}}=0,\underset{t\to +\infty }{lim\; sup}{\displaystyle \frac{〈H\left(t\right)〉}{t}}=0,\underset{t\to +\infty }{lim\; sup}{\displaystyle \frac{〈R\left(t\right)〉}{t}}=0,a.s.\hfill \end{array}$$

(20)

As a result, the disease has indeed disappeared.

Proof. 

We obtain the following set of equations by integrating model (1):

$$\begin{array}{cc}\hfill {\displaystyle \frac{S\left(t\right)-S\left(0\right)}{t}}& =\psi \Pi -{\displaystyle \frac{{\beta }_1〈S\left(t\right)A\left(t\right)〉}{〈N〉}}+\eta 〈V\left(t\right)〉-(\gamma +\chi )〈S\left(t\right)〉+{\displaystyle \frac{{\delta }_1}{t}}{\int }_0^{t}S\left(r\right)d{B}_1\left(r\right),\hfill \\ \hfill {\displaystyle \frac{V\left(t\right)-V\left(0\right)}{t}}& =\Pi (1-\psi )+\gamma 〈S\left(t\right)〉-{\displaystyle \frac{{\beta }_2〈V\left(t\right)A\left(t\right)〉}{〈N〉}}-(\eta +\chi )〈V\left(t\right)〉+{\displaystyle \frac{{\delta }_2}{t}}{\int }_0^{t}V\left(r\right)d{B}_2\left(r\right),\hfill \\ \hfill {\displaystyle \frac{A\left(t\right)-A\left(0\right)}{t}}& ={\displaystyle \frac{{\beta }_1〈S\left(t\right)A\left(t\right)〉}{〈N〉}}+{\displaystyle \frac{{\beta }_2〈V\left(t\right)A(t)〉}{〈N〉}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )〈A\left(t\right)〉+{\displaystyle \frac{{\delta }_3}{t}}{\int }_0^{t}A\left(r\right)d{B}_3\left(r\right),\hfill \\ \hfill {\displaystyle \frac{Z\left(t\right)-Z\left(0\right)}{t}}& ={\alpha }_1〈A\left(t\right)〉-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi )〈Z\left(t\right)〉+{\displaystyle \frac{{\delta }_4}{t}}{\int }_0^{t}Z\left(r\right)d{B}_4\left(r\right),\hfill \\ \hfill {\displaystyle \frac{H\left(t\right)-H\left(0\right)}{t}}& ={\psi }_1<A\left(t\right)+{\psi }_2〈Z\left(t\right)+〉-({\psi }_3+{\psi }_4+\chi )〈H\left(t\right)〉+{\displaystyle \frac{{\delta }_5}{t}}{\int }_0^{t}H\left(r\right)d{B}_5\left(r\right),\hfill \\ \hfill {\displaystyle \frac{R\left(t\right)-R\left(0\right)}{t}}& ={\gamma }_1〈Z\left(t\right)〉+{\alpha }_2〈A\left(t\right)〉+{\psi }_3〈H\left(t\right)〉-\chi 〈R\left(t\right)〉+{\displaystyle \frac{{\delta }_6}{t}}{\int }_0^{t}R\left(r\right)d+B_6\left(r\right).\hfill \end{array}$$

(21)

We apply Itô’s formula and the third equation of model (1), and we obtain the following form:

$$\begin{array}{cc}\hfill d(lnA\left(t\right))& =\left[{\displaystyle \frac{{\beta }_{1}SA}{N}}+{\displaystyle \frac{{\beta }_{2}VA}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )A\right]{\displaystyle \frac{1}{A}}dt-{\displaystyle \frac{{\delta }_3^2}{2}}dt+{\delta }_{3}d{B}_3\left(t\right),\hfill \\ \hfill & =\left[{\displaystyle \frac{{\beta }_{1}S}{N}}+{\displaystyle \frac{{\beta }_{2}V}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )\right]dt-{\displaystyle \frac{{\delta }_3^2}{2}}dt+{\delta }_{3}d{B}_3\left(t\right),\hfill \\ \hfill & =\left[{\displaystyle \frac{{\beta }_{1}S}{N}}+{\displaystyle \frac{{\beta }_{2}V}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}})\right]dt+{\delta }_{3}d{B}_3\left(t\right).\hfill \end{array}$$

(22)

By integrating Equation (22) over the interval $[0,t]$ and subsequently dividing both sides by t, we derive

$$\begin{array}{cc}\hfill lnA\left(t\right)-lnA\left(0\right)& ={\int }_0^t\left[{\displaystyle \frac{{\beta }_{1}S}{N}}+{\displaystyle \frac{{\beta }_{2}V}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}})\right]dt+{\delta }_{3}d{B}_3\left(t\right),\hfill \\ \hfill & \le \left[{\beta }_1+{\beta }_2-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}})\right]t+{\delta }_{3}d{B}_3\left(t\right),\hfill \\ \hfill & =\left({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}}\right)\left[R_0^E-1\right]t+{\delta }_{3}d{B}_3\left(t\right).\hfill \end{array}$$

(23)

In accordance with the strong law of large numbers, we can obtain

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}{\displaystyle \frac{B_3\left(t\right)}{t}}=0a.s.\end{array}$$

(24)

By taking the superior limit on both sides of Equation (23) and combining it with (17), we can obtain

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}sup{\displaystyle \frac{lnA\left(t\right)}{t}}\le \left({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}}\right)\left[R_0^E-1\right]<0,a.s.\end{array}$$

(25)

So, this mean that if $R_0^E<1$, we can obtain

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}A\left(t\right)=0a.s.\end{array}$$

Using Lemmas 1 and 2 and conditions $\left({\mathbb{C}}_1\right)$ and $\left({\mathbb{C}}_2\right)$ makes it simple to obtain

$$\underset{t\to \infty }{lim\; sup}{\displaystyle \frac{log〈A\left(t\right)〉}{t}}<0.$$

(26)

By once more employing Lemmas 1 and 2, along with conditions $\left({\mathbb{C}}_1\right)$ and $\left({\mathbb{C}}_2\right)$, and utilizing Equation (26), we have

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim\; sup}〈Z\left(t\right)〉=0,\mathrm{a}.\mathrm{s},\underset{t\to \infty }{lim\; sup}〈\mathsf{H}\left(\mathsf{t}\right)〉=0,\mathrm{a}.\mathrm{s},\underset{t\to \infty }{lim\; sup}〈\mathsf{R}\left(\mathsf{t}\right)〉=0,\mathrm{a}.\mathrm{s},\hfill \\ & \underset{t\to \infty }{lim\; sup}〈\mathsf{S}\left(\mathsf{t}\right)〉={\displaystyle \frac{\Pi (\eta +\chi \psi )}{\chi (\gamma +\eta +\chi )}},\mathrm{a}.\mathrm{s},\underset{t\to \infty }{lim\; sup}〈\mathsf{V}\left(\mathsf{t}\right)〉={\displaystyle \frac{\chi \Pi (1-\psi )+\gamma \psi }{\chi (\gamma +\eta +\chi )}},\mathrm{a}.\mathrm{s}.\hfill \end{array}$$

At this point, the proof of Theorem 2 has been completed. □

5. Stochastic Analysis of the Endemic State

When analyzing the dynamics of an epidemic, our concern extends beyond the mere possibility of the disease dying out. We also consider the scenario where the disease persists within a community. It is commonly acknowledged that the stochastic model being discussed has stable equilibria, including both disease-free and endemic states. Consequently, our attention in this section is directed toward the presence of a stationary distribution. This distribution serves as an indicator of whether the infection is spreading or not, and we will employ Khasminskii’s theory [30] for this purpose. To proceed, we will introduce several lemmas and definitions crucial for proving our main result.

Assuming a homogeneous Markov process $Y\left(t\right)$ in $R_{+}^n$,

$$dfY\left(t\right)=b\left(Y\right)dt+\sum _r^k{\sigma }_{r}d{\mathsf{B}}_r\left(t\right).$$

The following is the form of the diffusion matrix:

$$A\left(Y\right)=\left[a_{ij}\left(x\right)\right],a_{ij}\left(x\right)=\sum _{r=1}^k{\sigma }_r^i\left(x\right){\sigma }_j^r\left(x\right).$$

Lemma 3

([30]). Markov process $X\left(t\right)$ possesses a singular stationary distribution $m(·)$ under the condition that a bounded domain $U\subset {\mathbb{R}}^d$ with a smooth boundary alongside its closure $\overline{U}\in {\mathbb{R}}^d$ satisfies the following stringent criteria:

• Within the neighborhood $\left\{U\right\}$ and across its domain, the minimum eigenvalue of the matrix $\left\{A\right(t\left)\right\}$ is maintained at a specific distance from zero.
• When x is an element of the space ${\mathbb{R}}^d\setminus U$, and the mean time τ needed to reach U from x is not infinite, with ${sup}_{x\in K}{E}^x\tau$ remaining finite for every compact subset $K\subset {\mathbb{R}}^n$, and for any π-measurable function $f(·)$, we deduce the following:

$$P\left(\underset{T\to \infty }{lim}{\displaystyle \frac{1}{T}}{\int }_0^{T}f\left(Y^x\left(t\right)\right)dt={\int }_{R^d}f\left(x\right)\pi \left(dx\right)\right)=1,$$

where $x\in R^d.$

We must specify the parameter because of its importance.

$${\mathsf{R}}_0^S={\displaystyle \frac{{\beta }_1{\gamma }_2{\alpha }_1\psi }{(\gamma +\chi +\frac{{\delta }_1^2}{2})({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +\frac{{\delta }_3^2}{2})({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +\frac{{\delta }_4^2}{2})}}.$$

(27)

The Stationary Distribution and Ergodicity

Ergodicity and the existence of a stationary distribution will be discussed in this section.

Theorem 3.

If ${\mathsf{R}}_{\mathsf{0}}^{\mathsf{S}}>1$, then model (1)’s solution $\left(\mathsf{S}\right(\mathsf{t}),\mathsf{V}(\mathsf{t}),\mathsf{A}(\mathsf{t}),\mathsf{Z}(\mathsf{t}),\mathsf{H}(\mathsf{t}),\mathsf{R}(\mathsf{t}\left)\right)$ is ergodic. Additionally, there is a singular stationary distribution called $\pi (.)$.

Proof. 

To establish claim (2) in Lemma 3, we need to construct a $C^2$-function that maps from the space ${\mathbb{R}}_{+}^6$ to ${\mathbb{R}}_{+}$. The proposed function takes on the following structure:

$${\mathbb{V}}_{⊮}=S+V+A+Z+H+R-c_{1}lnS-c_{2}lnA-c_{3}lnZ.$$

Here, $c_1$, $c_2$, and $c_3$ represent positive real numbers; their values are currently unknown. We obtain the following relationships by applying the Itô formula to model (1):

$$\begin{array}{cc}\hfill & \mathbb{L}(S+V+A+Z+H+R)=\Pi -\chi (S+V+A+Z+H+R)-{\gamma }_{2}Z\left(t\right).\hfill \\ & (-c_{1}lnS)=-{\displaystyle \frac{c_1\psi \Pi }{S}}+{\displaystyle \frac{c_1{\beta }_{1}A\left(t\right)}{N}}-{\displaystyle \frac{c_1\eta V\left(t\right)}{S}}+c_1(\gamma +\chi )+{\displaystyle \frac{c_1{\delta }_1^2}{2}},\hfill \\ & (-c_{2}lnA)=-{\displaystyle \frac{c_2{\beta }_{1}S\left(t\right)}{N}}-{\displaystyle \frac{c_2{\beta }_{2}V\left(t\right)}{N}}+c_2({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )+{\displaystyle \frac{c_2{\delta }_3^2}{2}},\hfill \\ & (-c_{3}lnZ)=-{\displaystyle \frac{c_3{\alpha }_{1}A}{Z}}+c_3({\gamma }_1+{\psi }_2+{\gamma }_2+\chi )+{\displaystyle \frac{c_3{\delta }_4^2}{2}}.\hfill \end{array}$$

(28)

Therefore, we have

$$\begin{array}{cc}\hfill {\mathbb{LV}}_{⊮}& =\Pi -\chi N-{\gamma }_{2}Z-{\displaystyle \frac{c_1\psi \Pi }{S}}+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{1}S}{N}}-{\displaystyle \frac{c_2{\beta }_{2}V}{V}}-{\displaystyle \frac{c_3{\alpha }_{1}A}{Z}}\hfill \\ & +c_1(\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}})+c_2({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}})+c_3({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +{\displaystyle \frac{{\delta }_4^2}{2}}).\hfill \end{array}$$

We have

$$\begin{array}{cc}\hfill {\mathbb{LV}}_{⊮}& \le -4{\left[\left({\gamma }_{2}Z\right)\times \left({\displaystyle \frac{c_3{\alpha }_{1}A}{Z}}\right)\times \left({\displaystyle \frac{c_1\psi \Pi }{S}}\right)\times \left({\displaystyle \frac{c_2{\beta }_{1}S}{N}}\right)\right]}^{{\displaystyle \frac{1}{4}}}-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\hfill \\ & -{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}+c_1(\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}})+c_2({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}})\hfill \\ & +c_3({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +{\displaystyle \frac{{\delta }_4^2}{2}})+\Pi .\hfill \end{array}$$

Let

$$c_1\left(\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}\right)=c_2\left({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +{\displaystyle \frac{{\delta }_3^2}{2}}\right)=c_3\left({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +{\displaystyle \frac{{\delta }_4^2}{2}}\right)=\Pi .$$

The constants in this case are defined as follows:

$$\begin{array}{c}\hfill c_1={\displaystyle \frac{\Pi }{\left(\gamma +\chi +\frac{{\delta }_1^2}{2}\right)}},c_2={\displaystyle \frac{\Pi }{\left({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +\frac{{\delta }_3^2}{2}\right)}},c_3={\displaystyle \frac{\Pi }{\left({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +\frac{{\delta }_4^2}{2}\right)}}.\end{array}$$

(29)

Consequently,

$$\begin{array}{cc}\hfill {\mathbb{LV}}_{⊮}& \le -4\left[{\left({\displaystyle \frac{{\Pi }^4{\beta }_1{\alpha }_1{\gamma }_2\psi }{\left(\gamma +\chi +\frac{{\delta }_1^2}{2}\right)\left({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +\frac{{\delta }_3^2}{2}\right)\left({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +\frac{{\delta }_4^2}{2}\right)}}\right)}^{{\displaystyle \frac{1}{4}}}-\Pi \right]\hfill \\ \hfill & -\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}.\hfill \end{array}$$

After simplifying the expressions, we arrive at the subsequent inequality:

$$\begin{array}{cc}\hfill {\mathbb{LV}}_{⊮}& \le -4\Pi \left[{\left(R_0^S\right)}^{1/4}-1\right]-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}.\hfill \end{array}$$

(30)

Furthermore, we are able to obtain

$$\begin{array}{cc}\hfill {\mathbb{V}}_2& =c_4(S+V+A+Z+H+R-c_{1}lnS-c_{2}lnA-c_{3}lnZ)-lnS-lnV-lnH-lnR+S+V+A+Z+H+R,\hfill \\ \hfill & =(c_4+1)(S+V+A+Z+H+R)-(c_1{c}_4+1)lnS-lnV-c_2{c}_{4}lnA-c_3{c}_{4}lnZ-lnH-lnR.\hfill \end{array}$$

The constant $c_4$, with a positive value, will be defined at a later stage. This can be demonstrated in a useful way.

$$\underset{(S,V,A,Z,H,R)\in R_{+}^6\setminus U_k}{lim\; inf}{\mathbb{V}}_2(S,V,A,Z,H,R)=+\infty ,\mathrm{as}n\to \infty ,$$

(31)

Here, $U_n=({\displaystyle \frac{1}{n}},n)\times ({\displaystyle \frac{1}{n}},n)\times ({\displaystyle \frac{1}{n}},n)\times ({\displaystyle \frac{1}{n}},n)$. Subsequently, it needs to be demonstrated that ${\mathbb{V}}_2(S,V,A,Z,H,R)$ has at least one minimum value ${\mathbb{V}}_2(S_0,V_0,A_0,Z_0,H_0,R_0).$

For the variables $S,V,A,Z,H$, and R, representing the current state, the partial derivative of the function ${\mathbb{V}}_2(S,V,A,Z,H,R)$ can be expressed as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathbb{V}}_2(S,V,A,Z,H,R)}{\partial S}}& =c_4+1-{\displaystyle \frac{1+c_1{c}_4}{S}},\hfill \\ \hfill {\displaystyle \frac{\partial {\mathbb{V}}_2(S,V,A,Z,H,R)}{\partial V}}& =c_4+1-{\displaystyle \frac{1}{V}},\hfill \\ \hfill {\displaystyle \frac{\partial {\mathbb{V}}_2(S,V,A,Z,H,R)}{\partial A}}& =c_4+1-{\displaystyle \frac{c_2{c}_4}{A}},\hfill \\ \hfill {\displaystyle \frac{\partial {\mathbb{V}}_2(S,V,A,Z,H,R)}{\partial Z}}& =c_4+1-{\displaystyle \frac{c_3{c}_4}{Z}},\hfill \\ \hfill {\displaystyle \frac{\partial {\mathbb{V}}_2(S,V,A,Z,H,R)}{\partial H}}& =c_4+1-{\displaystyle \frac{1}{H}},\hfill \\ \hfill {\displaystyle \frac{\partial {\mathbb{V}}_2(S,V,A,Z,H,R)}{\partial R}}& =c_4+1-{\displaystyle \frac{1}{R}}.\hfill \end{array}$$

It becomes evident that the function ${\mathbb{V}}_2$ exhibits a single stagnation point:

$$(S_0,V_0,A_0,Z_0,H_0,R_0)=\left({\displaystyle \frac{c_1{c}_4+1}{c_4+1}},{\displaystyle \frac{1}{c_4+1}},{\displaystyle \frac{c_4{c}_2}{c_4+1}},{\displaystyle \frac{c_3{c}_4}{c_4+1}},{\displaystyle \frac{1}{c_4+1}},{\displaystyle \frac{1}{c_4+1}}\right).$$

(32)

Furthermore, the matrix of the second partial derivatives (Hessian matrix) for ${\mathbb{V}}_2(S,V,A,Z,H,R)$ at $({\mathsf{S}}_{\mathsf{0}},{\mathsf{V}}_{\mathsf{0}},{\mathsf{A}}_{\mathsf{0}},{\mathsf{Z}}_{\mathsf{0}},{\mathsf{H}}_{\mathsf{0}},{\mathsf{R}}_{\mathsf{0}})$ is

$$B=\left[\begin{array}{cccccc}{\displaystyle \frac{c_1{c}_4+1}{S^2\left(0\right)}}& 0& 0& 0& 0& 0\\ 0& {\displaystyle \frac{1}{V^2\left(0\right)}}& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{c_2{c}_4}{A^2\left(0\right)}}& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{c_3{c}_4}{Z^2\left(0\right)}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{1}{H^2\left(0\right)}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{1}{R^2\left(0\right)}}\end{array}\right].$$

(33)

The Hessian matrix is positive definite. Consequently, ${\mathbb{V}}_2(S,V,A,Z,H,R)$ achieves a minimum value at ${\mathbb{V}}_2(S_0,V_0,A_0,Z_0,H_0,R_0)$. Given Equation (31) and acknowledging the continuity of ${\mathbb{V}}_2(S,V,A,Z,H,R)$, it can be inferred that ${\mathbb{V}}_2(S,V,A,Z,H,R)$ possesses a unique minimum value ${\mathbb{V}}_2(S_0,V_0,A_0,Z_0,H_0,R_0)$ within the domain ${\mathbb{R}}_{+}^6$.

A non-negative $C^2$-function $\mathbb{V}:R_{+}^6\to R_{+}$ shall be introduced as follows:

$$\mathbb{V}(S,V,A,Z,H,R)={\mathbb{V}}_2(S,V,A,Z,H,R)-{\mathbb{V}}_2\left(S\right(0),V(0),A(0),Z(0),H(0),R(0\left)\right).$$

We carried out a precise analysis using the stochastic model and Ito’s formula:

$$\begin{array}{cc}\hfill \mathbb{LV}& \le c_4\left\{-4\Pi \left[{\left(R_0^S\right)}^{1/4}-1\right]-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}\hfill \\ \hfill & -{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}\hfill \\ \hfill & -{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \end{array}$$

(34)

The relation stated above will result in the following assertion:

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}\hfill \\ \hfill & -{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}\hfill \\ \hfill & -{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \end{array}$$

(35)

where

$$c_5=4\Pi \left[{\left(R_0^S\right)}^{1/4}-1\right]$$

Our next step will be to define the set forms:

$$D=\left\{{ϵ}_1<S{\displaystyle \frac{1}{{ϵ}_2}},{ϵ}_3<V<{\displaystyle \frac{1}{{ϵ}_4}},{ϵ}_5<A<{\displaystyle \frac{1}{{ϵ}_6}},{ϵ}_7<Z<{\displaystyle \frac{1}{{ϵ}_8}}{ϵ}_9<H<{\displaystyle \frac{1}{{ϵ}_{10}}}{ϵ}_{11}<R<{\displaystyle \frac{1}{{ϵ}_{12}}}\right\}.$$

Here, ${ϵ}_i$ for $(i=1,2,3,4,\cdots ,12)$ is a small positive real number that is currently undetermined. To maintain simplicity, it is advantageous to divide the region ${\mathbb{R}}_{+}^6\setminus D$ into the following subsections:

$$\begin{array}{cc}\hfill D_1& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,0<S\le {ϵ}_1\right\},\hfill \\ \hfill D_2& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,0<V\le {ϵ}_3,S>{ϵ}_2\right\},\hfill \\ \hfill D_3& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,0<A\le {ϵ}_5,V>{ϵ}_4\right\},\hfill \\ \hfill D_4& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,0<Z\le {ϵ}_7,A>{ϵ}_6\right\},\hfill \\ \hfill D_5& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,0<H\le {ϵ}_9,Z>{ϵ}_8\right\},\hfill \\ \hfill D_6& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,0<R\le {ϵ}_{11},Z>{ϵ}_{12}\right\},\hfill \\ \hfill D_7& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,S\ge {\displaystyle \frac{1}{{ϵ}_2}}\right\},\hfill \\ \hfill D_8& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,V\ge {\displaystyle \frac{1}{{ϵ}_4}}\right\},\hfill \\ \hfill D_9& =\left\{((S,V,A,Z,H,R)\in R_{+}^6,A\ge {\displaystyle \frac{1}{{ϵ}_6}}\right\},\hfill \\ \hfill D_{10}& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,Z\ge {\displaystyle \frac{1}{{ϵ}_8}}\right\},\hfill \\ \hfill D_{11}& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,H\ge {\displaystyle \frac{1}{{ϵ}_{10}}}\right\},\hfill \\ \hfill D_{12}& =\left\{(S,V,A,Z,H,R)\in R_{+}^6,R\ge {\displaystyle \frac{1}{{ϵ}_{12}}}\right\}.\hfill \end{array}$$

Continuing our analysis, it is imperative to affirm the negativity of $\mathbb{LV}(S,V,A,Z,H,R)$ across the domain ${\mathbb{R}}_{+}^6\setminus D$. Put simply, it is necessary for each subregion within this set to exhibit the negativity of the mentioned function.

Case 1. If $(S,V,A,Z,H,R)\in D_1$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\psi \Pi }{S}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\psi \Pi }{{ϵ}_1}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_1>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\psi \Pi }{{ϵ}_1}}<0$. Hence, as a result, we can conclude that $\mathbb{LV}<0$ for any $(\mathsf{S},\mathsf{V},\mathsf{A},\mathsf{Z},\mathsf{H},\mathsf{R})\in D_1.$

Case 2. If $(S,V,A,Z,H,R)\in D_2$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\gamma S}{V}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\gamma {ϵ}_4}{{ϵ}_2}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_2>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\gamma {ϵ}_4}{{ϵ}_2}}<0$. Hence, as a result, we can conclude that $\mathbb{LV}<0$ for any $(\mathsf{S},\mathsf{V},\mathsf{A},\mathsf{Z},\mathsf{H},\mathsf{R})\in D_2.$

Case 3. If $(S,V,A,Z,H,R)\in D_3$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\Pi (1-\psi )}{V}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\Pi (1-\psi )}{{ϵ}_3}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_3>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\Pi (1-\psi )}{{ϵ}_3}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(\mathsf{S},\mathsf{V},\mathsf{A},\mathsf{Z},\mathsf{H},\mathsf{R})\in D_3.$

Case 4. If $(S,V,A,C,H,R)\in D_4$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\eta V}{S}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\eta {ϵ}_2}{{ϵ}_4}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_4>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}C}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\eta {ϵ}_2}{{ϵ}_4}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_4.$

Case 5. If $(S,V,A,Z,H,R)\in D_5$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}C}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_{2}A}{N}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_2{ϵ}_5}{{ϵ}_5}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_5>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_2{ϵ}_5}{{ϵ}_5}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_5.$

Case 6. If $(S,V,A,Z,H,R)\in D_6$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{\psi A}{H}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_1{ϵ}_{10}}{{ϵ}_6}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_6>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_1{ϵ}_{10}}{{ϵ}_6}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_6.$

Case 7. If $(S,V,A,Z,H,R)\in D_7$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_{2}A}{N}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_2{ϵ}_5}{{ϵ}_7}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_5>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_2{ϵ}_5}{{ϵ}_7}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_7.$

Case 8. If $(S,V,A,Z,H,R)\in D_8$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\gamma }_{1}Z}{R}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\gamma }_1{ϵ}_{10}}{{ϵ}_8}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_6>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\gamma }_1{ϵ}_{10}}{{ϵ}_8}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_8.$

Case 9. If $(S,V,A,Z,H,R)\in D_9$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_{1}A}{H}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_1{ϵ}_5}{{ϵ}_9}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_9>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_1{ϵ}_5}{{ϵ}_9}}.<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_9.$

Case 10. If $(S,V,A,Z,H,R)\in D_{10}$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_{3}H}{R}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_3{ϵ}_{12}}{{ϵ}_{10}}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_{10}>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_3{ϵ}_{12}}{{ϵ}_{1}0}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_{10}.$

Case 11. If $(S,V,A,Z,H,R)\in D_{11}$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}C}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_{3}H}{R}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_3{ϵ}_9}{{ϵ}_{11}}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_{11}>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\psi }_3{ϵ}_9}{{ϵ}_{11}}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_{11}.$

Case 12. If $(S,V,A,Z,H,R)\in D_{12}$, using Equation (35), we obtain

$$\begin{array}{cc}\hfill \mathbb{LV}& \le -c_4{c}_5+c_4\left\{-\chi N+{\displaystyle \frac{c_1{\beta }_{1}A}{N}}-{\displaystyle \frac{c_1\eta V}{S}}-{\displaystyle \frac{c_2{\beta }_{2}V}{N}}\right\}-{\displaystyle \frac{\psi \Pi }{S}}+{\displaystyle \frac{{\beta }_{1}A}{N}}-{\displaystyle \frac{\eta V}{S}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}-{\displaystyle \frac{\Pi (1-\psi )}{V}}-{\displaystyle \frac{\gamma S}{V}}-{\displaystyle \frac{{\beta }_{2}A}{N}}\hfill \\ \hfill & +\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}-{\displaystyle \frac{{\psi }_{1}A}{H}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}-{\displaystyle \frac{{\gamma }_{1}Z}{R}}-{\displaystyle \frac{{\alpha }_{2}A}{R}}-{\displaystyle \frac{{\psi }_{3}H}{R}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -\chi N-{\gamma }_{2}Z.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_{2}A}{N}}.\hfill \\ \hfill & \le -c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_2{ϵ}_5}{{ϵ}_{12}}}.\hfill \end{array}$$

It is possible to pick a sufficiently small ${ϵ}_{12}>0$ such that $-c_4{c}_5+c_4\left\{{\displaystyle \frac{c_1{\beta }_{1}A}{N}}\right\}+{\displaystyle \frac{{\beta }_{1}A}{N}}+\gamma +\chi +{\displaystyle \frac{{\delta }_1^2}{2}}+\eta +\chi +{\displaystyle \frac{{\delta }_2^2}{N}}+{\displaystyle \frac{{\psi }_{2}Z}{H}}+{\psi }_3+{\psi }_4+\chi +{\displaystyle \frac{{\delta }_5^2}{2}}+\chi +{\displaystyle \frac{{\delta }_6^2}{2}}+\psi -{\displaystyle \frac{{\beta }_2{ϵ}_5}{{ϵ}_{12}}}<0$. As a result, we can conclude that $\mathbb{LV}<0$ for any $(S,V,A,Z,H,R)\in D_{12}.$

$$\mathbb{LV}(S,V,A,Z,H,R)<-W<0\mathrm{for}\mathrm{all}\left(S,V,A,Z,H,R\right)\in R_{+}^6\setminus D.$$

Hence,

$$\begin{array}{cc}\hfill d\mathbb{V}(S,V,A,Z,H,R)& <-Wdt+\left[(c_4+1)S-(c_1{c}_4+1){\delta }_1\right]d{B}_1\left(t\right)+[(c_4+1)V-{\delta }_2]d{B}_2\left(t\right)\hfill \\ \hfill & +[(c_4+1)A-c_2{c}_4{\delta }_3]d{B}_3\left(t\right)+[(c_4+1)Z-c_3{c}_4{\delta }_4]d{B}_4\left(t\right)\hfill \\ \hfill & +[(c_4+1)H-{\delta }_5]B_5\left(t\right)+[(c_4+1)R-{\delta }_6]d{B}_6\left(t\right).\hfill \end{array}$$

(36)

Suppose that $(S_0,V_0,A_0,Z_0,H_0,R_0)=(x_1,x_2,x_3,x_4,x_5,x_6)=x\in R_{+}^6\setminus D$. Further, if ${\tau }^x$ represents the duration it takes for the curve originating from x to reach D, then

$${\tau }_n=inf\{t:|X\left(t\right)|=n\}\mathrm{and}{\tau }^{\left(n\right)}\left(t\right)=min\{{\tau }^x,t,{\tau }_n\}.$$

By integrating Equation (36) over the interval $[0,{\tau }^{\left(n\right)}\left(t\right)]$, employing Dynkin formula, and subsequently calculating the expectation, we have

$$\begin{array}{cc}\hfill & E\mathbb{V}(S\left({\tau }^{\left(n\right)}\left(t\right)\right),V\left({\tau }^{\left(n\right)}\left(t\right)\right),A\left({\tau }^{\left(n\right)}\left(t\right)\right),Z\left({\tau }^{\left(n\right)}\left(t\right)\right),H\left({\tau }^{\left(n\right)}\left(t\right)\right)),R\left({\tau }^{\left(n\right)}\left(t\right)\right))-V\left(x\right),\hfill \\ \hfill & =E{\int }_0^{\tau \left(n\right)\left(t\right)}\mathbb{LV}\left(S\left(u\right),V\left(u\right),A\left(u\right),Z\left(u\right),H\left(u\right),R\left(u\right)\right)du,\hfill \\ \hfill & \le E{\int }_0^{\tau \left(n\right)\left(t\right)}-Wdu=-W\mathbb{E}{\tau }^{\left(n\right)}\left(t\right).\hfill \end{array}$$

As $\mathbb{V}\left(x\right)$ is non-negative,

$$E{\tau }^{\left(n\right)}\left(t\right)\le {\displaystyle \frac{\mathbb{V}\left(x\right)}{W}}.$$

Referring to the proof of Theorem 3, it has been demonstrated that $P\{{\tau }_e=\infty \}=1$. Hence, we need to ensure the consistency of the structure of the model (1). Furthermore, as we allow both $t,n\to \infty$, it is almost certain that $\tau \left(n\right)\left(t\right)\to {\tau }^x$. Additionally, through the application of Fatou’s Lemma, we obtain that

$$E{\tau }^{\left(n\right)}\left(t\right)\le {\displaystyle \frac{\mathbb{V}\left(x\right)}{W}},$$

is finite, indicating that $su{p}_{x\in K}E{\tau }^x$ is also finite, where K represents a compact subset of $R_{+}^6$. As a result, statement (2) in Lemma 3 is proved. Moreover, model (1) is characterized by the diffusion matrix

$$\mathsf{B}=\left[\begin{array}{cccccc}{\delta }_1^2{S}^2& 0& 0& 0& 0& 0\\ 0& {\delta }_2^2{V}^2& 0& 0& 0& 0\\ 0& 0& {\delta }_3^2{A}^2& 0& 0& 0\\ 0& 0& 0& {\delta }_4^2{Z}^2& 0& 0\\ 0& 0& 0& 0& {\delta }_5^2{H}^2& 0\\ 0& 0& 0& 0& 0& {\delta }_6^2{R}^2\end{array}\right].$$

Picking $M={min}_{(S,V,A,Z,H,R)\in \overline{D}\in R_{+}^6}\{{\delta }_1^2{S}^2,{\delta }_2^2{V}^2,{\delta }_3^2{A}^2,{\delta }_4^2{Z}^2,{\delta }_5^2{H}^2,{\delta }_6^2{R}^2\},$ we obtain

$$\sum _{i,j=1}^6{a}_{ij}(S,V,A,Z,H,R){\xi }_i{\xi }_j={\delta }_1^2{S}^2{\xi }^2+{\delta }_2^2{V}^2{\xi }_2^2+{\delta }_3^2{A}^2{\xi }^2+{\delta }_4^2{\xi }_4^2{Z}^2+H^2{\xi }^2+{\delta }_5^2+{\delta }_6^2{\xi }_6^2{R}^2\ge M{\left|\xi \right|}^2,(S,V,A,Z,H,R)\in \overline{D},$$

This concludes verification of statement (1) of Lemma 3, where $\xi =({\xi }_1,{\xi }_2,{\xi }_3,{\xi }_4,{\xi }_5,{\xi }_6)$ belongs to the domain $R_{+}^6$. Consequently, this demonstrates the system’s ergodic behavior and establishes the existence of a single stationary distribution for the system. With this, the theorem’s proof is complete. □

6. Optimal Control

Mathematical biology is increasingly concerned with optimal control theory. These tools can be used to find an effective control strategy by considering suitable control functions [31,32]. Mathematics, dynamical system theory, and economics [33] broadly apply the control approach. This work [34] deduces several optimality conditions and other control theory features. Infectious diseases can be prevented or minimized using several strategies. The most effective method is to reduce the contagious disease with the least cost [35]. Using this method, real-world problems can be linked to the underlying system mathematically. Since most issues require physical knowledge, this step can be very challenging. In numerous scenarios, our available information is restricted to experimental data and fundamental insights from the literature. To develop an accurate model, reviewing a large amount of bibliographic material and following several trials is necessary. To create the best model for experimental results, one must start with a simple model and work their way up to a more complex one. It has been developed and investigated [36] that stochastic delay models contain n competing species.

Additionally, many authors have developed effective harvesting strategies. Liu and Meng [37] explored optimal harvesting strategies by examining the qualitative behavior of stochastic delay systems. The authors employed Hessian matrices, ergodic theory, and optimal harvesting theory to derive empirically supported conclusions regarding the upper limit of sustainable yield. For more information on stochastic optimal control systems, please see [31,32] and references.

It is discussed in this section how model (1) can be controlled equivalently. We will formulate the corresponding stochastic and deterministic control systems by considering $u_1\left(t\right)$ and $u_2\left(t\right)$ as the control measures.

In models (1) and (2), the control variables have the follopwing specific interpretations:

• The control measure $u_1\left(t\right)$ represents vaccination to lower the vulnerable population over time.
• The variable $u_2\left(t\right)$ represents the treatment of HBV patients, leading to a reduction in infected individuals within the community.

The objective of this study is to reduce the count of infected and susceptible individuals by implementing an effective control strategy, simultaneously enhancing the recovery rate.

6.1. Optimal Control of Deterministic Model (2)

In the same way as in [24,25], we apply control theory tools here to decrease the spread of HBV. With $u_1$, $u_2$, and $u_3$ as inputs, we obtain the optimal strategy for (2).

To eliminate the hepatitis B virus, we use model (2) and implement $u_1$, $u_2$, and $u_3$ (three control measures). The objective functions are minimized as follows:

$$\begin{array}{c}\hfill \mathfrak{J}(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right))={\int }_0^T\left[w_{1}A+w_{2}Z+{\displaystyle \frac{1}{2}}{w}_3{u}_1^2\left(t\right)+{\displaystyle \frac{1}{2}}{w}_4{u}_2^2\left(t\right)+{\displaystyle \frac{1}{2}}{w}_5{u}_3^2\left(t\right)\right]dt.\end{array}$$

(37)

Equations (37) and (38) are called optimal control systems per the following system:

$$\begin{array}{cc}\hfill S^{\prime }\left(t\right)& =\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi +u_1\left(t\right))S\left(t\right),\hfill \\ \hfill V^{\prime }\left(t\right)& =\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right),\hfill \\ \hfill A^{\prime }\left(t\right)& ={\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +u_2\left(t\right))A\left(t\right),\hfill \\ \hfill Z^{\prime }\left(t\right)& ={\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +u_3\left(t\right))Z\left(t\right),\hfill \\ \hfill H^{\prime }\left(t\right)& ={\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right),\hfill \\ \hfill R^{\prime }\left(t\right)& ={\gamma }_{1}Z\left(t\right)+({\alpha }_2+u_2\left(t\right))A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)+u_1\left(t\right)S\left(t\right).\hfill \end{array}$$

(38)

$$\begin{array}{c}\hfill S\left(0\right)>0,V\left(0\right)\ge 0,A\left(0\right)\ge 0,Z\left(0\right)\ge 0,H\left(0\right)\ge 0,R\left(0\right)>0.\end{array}$$

(39)

In Equation (37), the terms $w_1\left(t\right)$ and $w_2\left(t\right)$ are the cost incurred due to the infected population. $w_3\left(t\right)$, $w_4\left(t\right)$, and $w_5\left(t\right)$ are positive constants that describe the cost of the control variables vaccination and treatment. Equation (37) represents the total cost, and the integrand function $L_1=w_{1}A+w_{2}Z+{\displaystyle \frac{1}{2}}{w}_3{u}_1^2\left(t\right)+{\displaystyle \frac{1}{2}}{w}_4{u}_2^2\left(t\right)+{\displaystyle \frac{1}{2}}{w}_5{u}_3^2\left(t\right)$ represents the actual and chronic cost at the time t. The notations for $w_i$ in Equation (37) represent positive real numbers for $i=1,2,...,5$. These symbols represent the weights assigned to infectious agents, vaccination, and treatment costs. Guided by the objective function, the aim is to minimize the number of infected individuals and optimize their recovery, all while managing costs effectively. Specifically, our objective is to define $(u_1,u_2,u_3)$ as control variables in the following manner:

$$\begin{array}{c}\hfill \mathfrak{J}(u_1^{\ast },u_2^{\ast },u_3^{\ast })=min\{\mathfrak{J}(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right));u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)\in U\}.\end{array}$$

(40)

System (2) has a control set U that represents the control set

$$\begin{array}{c}\hfill U:=\{(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right))\mid u_i\left(t\right)\mathrm{are}\mathrm{Lebesgue}\mathrm{measureable}\mathrm{on}[0,T]0\le u_1\left(t\right),u_2\left(t\right)\le 1,u_3\left(t\right)\le 1\}.\end{array}$$

(41)

Defining control measures is the first step toward moving forward.

6.1.1. Existence of Solution

This study aimed to determine whether there is a solution to control problems (38) and (39). A non-negative solution of the state system is found through non-negative initial data and Lebesgue measurable control variables [38].

Let

$$\begin{array}{c}\hfill {\displaystyle \frac{d\phi }{dt}}=L\phi +\chi \left(\phi \right),\end{array}$$

(42)

where

$$\begin{array}{cc}\hfill & \phi =\left[\begin{array}{c}S\left(t\right)\\ V\left(t\right)\\ A\left(t\right)\\ Z\left(t\right)\\ H\left(t\right)\\ R\left(t\right)\\ \end{array}\right],\hfill \\ & L=\left[\begin{array}{cccccc}-(\gamma +\chi +u_1\left(t\right))& 0& 0& 0& 0& 0\\ \gamma & -(\eta +\chi )& 0& 0& 0& 0\\ 0& 0& -({\alpha }_1+{\alpha }_2{\psi }_1+\chi +u_2\left(t\right))& 0& 0& 0\\ 0& 0& {\alpha }_1& -({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +u_3\left(t\right))& 0& 0\\ 0& 0& {\psi }_1& {\psi }_2& -({\psi }_1+{\psi }_4+\chi )& 0\\ u_1\left(t\right)& 0& ({\alpha }_2+u_2\left(t\right))& {\gamma }_1& {\chi }_3& -\chi \\ \end{array}\right],\hfill \\ \hfill & \chi \left(\phi \right)=\left[\begin{array}{c}\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}\\ \Pi (1-\psi )-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}\\ {\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}\\ 0\\ 0\\ 0\\ \end{array}\right].\hfill \end{array}$$

(43)

Nonlinear system (42) has a bounded coefficient. Letting

$$\begin{array}{c}\hfill F\left(\phi \right)=L\left(\phi \right)+G\left(\phi \right),\end{array}$$

(44)

fulfills

$$\begin{array}{cc}\hfill \mid G\left({\phi }_1\right)-G\left({\phi }_2\right)\mid & \le \mid S_1\left(t\right)-S_2\left(t\right)\mid +m_2\mid V_1\left(t\right)-V_2\left(t\right)\mid +m_3\mid A_1\left(t\right)-A_2\left(t\right)\mid +m_4\mid Z_1\left(t\right)-Z_2\left(t\right)\mid \hfill \\ & +m_5\mid H_1\left(t\right)-H_2\left(t\right)\mid +m_6\mid R_1\left(t\right)-R_2\left(t\right)\mid ,\hfill \\ & \le M(\mid S_1\left(t\right)-S_2\left(t\right)\mid +\mid V_1\left(t\right)-V_2\left(t\right)\mid +\mid A_1\left(t\right)-A_2\left(t\right)\mid +\mid Z_1\left(t\right)-Z_2\left(t\right)\mid \hfill \\ & +\mid H_1\left(t\right)-H_2\left(t\right)\mid +\mid R_1\left(t\right)-R_2\left(t\right)\mid ).\hfill \end{array}$$

(45)

where $M=max{m}_i,i=1,2,3,4,5,6$ is independent of the state of (38). Additionally, we observe that $\mid F\left({\phi }_1\right)-F\left({\phi }_2\right)\mid \le \mid {\phi }_1-{\phi }_2\mid N$, where $N=max\left|\mathcal{L}\right|,M<\infty$. This ensures that function G is uniformly continuous in a Lipschitz sense. A solution to model (44) exists according to the definitions of the controls; the state variables are as follows: $S\left(0\right)>0,V\left(0\right)\ge 0,A\left(0\right)\ge 0,Z\left(0\right)\ge 0,H\left(0\right)\ge 0,R\left(0\right)>0$. Regarding the optimal controlling existence, the results below are worthy of stating and proving.

Theorem 4.

Control problems (37)–(40) can be solved using a vector $u^{\ast }\left(t\right)=(u_1^{\ast }\left(t\right),u_2^{\ast }\left(t\right),u_3^{\ast }\left(t\right))\in U$ for control.

Proof. 

We will follow [24] to derive the theorem’s conclusion. By doing this, we obtain the following results:

• For all values of t, the states control’s are non-negative.
• Closed and convex sets are defined in (41).
• The boundedness assures the property of compactness of control model (38).
• In expression (37), the integrand exhibits convexity concerning the control functions, ensuring the existence of optimal control variables $(u_1^{\ast },u_2^{\ast },u_3^{\ast })$.

□

6.1.2. Optimality Condition

We must create a Lagrangian and Hamiltonian to determine the optimality of systems (37)–(39). The state variables are ($x=(S,V,A,Z,H,R)$, and the control variables are $u=(u_1,u_2)$. A Lagrangian L can be defined with these variables:

$$\begin{array}{c}\hfill L(x,u)={\int }_0^T\left[w_{1}A\left(t\right)+w_{2}Z\left(t\right)+{\displaystyle \frac{1}{2}}{w}_3{u}_1^2+{\displaystyle \frac{1}{2}}{w}_4{u}_2^2+{\displaystyle \frac{1}{2}}{w}_5{u}_3^2\right]dt.\end{array}$$

(46)

A Hamiltonian $\mathbb{H}$ function becomes

$$\begin{array}{c}\hfill \mathbb{H}(\lambda ,u,x)=-L(u,x)+\lambda .g(u,x),\mathrm{where}\\ \hfill \lambda =({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\lambda }_6),g(u,x)& =(g_1(u,x),g_2(u,x),g_3(u,x),g_4(u,x),g_5(u,x),g_6(u,x)),\mathrm{having}\hfill \\ \hfill g_1(x,u)& =\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi +u_1\left(t\right))S\left(t\right),\hfill \\ \hfill g_2(x,u)& =\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right),\hfill \\ \hfill g_3(x,u)& ={\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +u_2\left(t\right))A\left(t\right),\hfill \\ \hfill g_4(x,u)& ={\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +u_3\left(t\right))Z\left(t\right),\hfill \\ \hfill g_5(x,u)& ={\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right),\hfill \\ \hfill g_6(x,u)& ={\gamma }_{1}Z\left(t\right)+({\alpha }_2+u_2\left(t\right))A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)+u_1\left(t\right)S\left(t\right).\hfill \end{array}$$

As a result, we will have a Hamiltonian of the following form for our control system:

$$\begin{array}{cc}\hfill H(\lambda ,u,x)& =w_{1}A+w_{2}Z+{\displaystyle \frac{1}{2}}{w}_3{u}_1^2\left(t\right)+{\displaystyle \frac{1}{2}}{w}_4{u}_2^2\left(t\right)+{\displaystyle \frac{1}{2}}{w}_5{u}_3^2\left(t\right)\hfill \\ & +{\lambda }_1\left(t\right)\left(\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi +u_1\left(t\right))S\left(t\right)\right)\hfill \\ & +{\lambda }_2\left(t\right)\left(\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right)\right)\hfill \\ & +{\lambda }_3\left(t\right)\left({\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +u_2\left(t\right))A\left(t\right)\right)\hfill \\ & +{\lambda }_4\left(t\right)\left({\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +u_3\left(t\right))Z\left(t\right)\right)\hfill \\ & +{\lambda }_5\left(t\right)\left({\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right)\right)\hfill \\ & +{\lambda }_6\left(t\right)\left({\gamma }_{1}Z\left(t\right)+({\alpha }_2+u_2\left(t\right))A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)+u_1\left(t\right)S\left(t\right)\right).\hfill \end{array}$$

To determine an optimal control solution for the problem, the Pontryagin maximum principle [39,40] establishes the existence of a Lagrange multiplier vector $\lambda$, such as

$$\begin{array}{c}\hfill \left\{{\displaystyle \frac{d\lambda \left(t\right)}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial x}}(x^{\ast }\left(t\right),u^{\ast }\left(t\right),\lambda \left(t\right)),{\displaystyle \frac{\partial \mathbb{H}(x^{\ast },u^{\ast },\lambda )}{\partial u}}=0\right\},\end{array}$$

(47)

with maximality and transversality conditions

$$\begin{array}{c}\hfill \mathbb{H}(u^{\ast }\left(t\right),\lambda \left(t\right),x^{\ast }\left(t\right))=su{p}_{u\in [0,l_1]\times [0,l_2]}\mathbb{H}(u,\lambda \left(t\right),x^{\ast }\left(t\right)),\end{array}$$

(48)

and

$$\begin{array}{c}\hfill \lambda \left(T\right)=0,\end{array}$$

(49)

Theorem 5.

Let $S^{\ast },V^{\ast },A^{\ast },Z^{\ast },H^{\ast }$, and $R^{\ast }$ denote the optimal state solution with the associated optimal variables $(u_1^{\ast },u_2^{\ast },u_3^{\ast })$ for optimal control problems (37)–(39). In this context, there exists am adjoint model that satisfies

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\lambda }_1\left(t\right)}{dt}}& =({\lambda }_1\left(t\right)-{\lambda }_3\left(t\right)){\displaystyle \frac{{\beta }_1{A}^{\ast }\left(t\right)}{N}}-u_1^{\ast }\left(t\right){\lambda }_6\left(t\right)+(\gamma +\chi +u_1^{\ast }\left(t\right)){\lambda }_1\left(t\right)-\gamma {\lambda }_2\left(t\right),\hfill \\ \hfill {\displaystyle \frac{d{\lambda }_2\left(t\right)}{dt}}& =\eta ({\lambda }_2\left(t\right)-{\lambda }_1\left(t\right))+\chi {\lambda }_2\left(t\right)+({\lambda }_2\left(t\right)-{\lambda }_3\left(t\right)){\displaystyle \frac{{\beta }_2{A}^{\ast }\left(t\right)}{N}},\hfill \\ \hfill {\displaystyle \frac{d{\lambda }_3\left(t\right)}{dt}}& =-w+({\lambda }_1\left(t\right)-{\lambda }_3\left(t\right)){\displaystyle \frac{{\beta }_1{S}^{\ast }\left(t\right)}{N}}+({\lambda }_2\left(t\right)-{\lambda }_3\left(t\right)){\displaystyle \frac{{\beta }_2{V}^{\ast }\left(t\right)}{N}}+({\lambda }_3\left(t\right)-{\lambda }_4\left(t\right)){\alpha }_1\hfill \\ & +({\lambda }_3\left(t\right)-{\lambda }_6\left(t\right)){\alpha }_2+({\lambda }_3\left(t\right)-{\lambda }_5\left(t\right)){\psi }_1+(\chi +u_2^{\ast }\left(t\right)){\lambda }_3\left(t\right),\hfill \\ \hfill {\displaystyle \frac{d{\lambda }_4\left(t\right)}{dt}}& =({\lambda }_4\left(t\right)-{\lambda }_6\left(t\right)){\gamma }_1+({\lambda }_4\left(t\right)-{\lambda }_5\left(t\right)){\psi }_2\left(t\right)+(\chi +u_3^{\ast }\left(t\right)){\lambda }_4\left(t\right),\hfill \\ \hfill {\displaystyle \frac{d{\lambda }_5\left(t\right)}{dt}}& =({\lambda }_5\left(t\right)-{\lambda }_6\left(t\right)){\psi }_3+({\psi }_4+\chi ){\lambda }_5\left(t\right),\hfill \\ \hfill {\displaystyle \frac{{\lambda }_6\left(t\right)}{dt}}& =\chi {\lambda }_6\left(t\right).\hfill \end{array}$$

(50)

In incorporating terminal constraints,

$$\begin{array}{c}\hfill {\lambda }_i\left(T\right)=0,i=\{1,2,3,4,5,6\}.\end{array}$$

(51)

The optimal values of the control measures are determined by

$$\begin{array}{c}\hfill u_1^{\ast }\left(t\right)=max\left\{0,min\{{\displaystyle \frac{[{\lambda }_6\left(t\right)-{\lambda }_1\left(t\right)]S^{\ast }}{w_2}},0\},1\right\},\end{array}$$

(52)

$$\begin{array}{c}\hfill u_2^{\ast }\left(t\right)=max\left\{0,min\{{\displaystyle \frac{[{\lambda }_6\left(t\right)-{\lambda }_3\left(t\right)]A^{\ast }}{w_3}},0\},1\right\},\end{array}$$

(53)

and

$$\begin{array}{c}\hfill u_3^{\ast }\left(t\right)=max\left\{0,min\{{\displaystyle \frac{[{\lambda }_6\left(t\right)-{\lambda }_4\left(t\right)]Z^{\ast }}{w_4}},0\},1\right\}.\end{array}$$

(54)

Proof. 

According to the Pontryagin principle, Equation (50) for the deputy is derived from Section 6.1.2. From Equation (49), we can obtain terminal condition (51). By differentiating the Hamiltonian function with respect to each control variable $u_i$ and subsequently solving the equations ${\displaystyle \frac{\partial H}{\partial u_i}}=0$ for $i=1,2$, we can obtain the optimal values of the control measures. It is easy to prove the theorem in Equations (52) and (53) by assuming the lower and upper limits of the controls.

Leveraging the adjoint Equation (50) along with its corresponding conditions (39) and (51), we successfully address the control problem for both the control and state variables. The control measures are also characterized. To continue the study of control theory, we will implement the same measures of control as we did in model (1). □

6.2. Optimal Control of Stochastic Model (1)

The stochastic counterpart of model (1) is formulated as follows:

$$\begin{array}{cc}\hfill dS\left(t\right)& =\left[\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi +u_1\left(t\right))S\left(t\right)\right]dt+{\delta }_{1}S\left(t\right)B_1\left(t\right),\hfill \\ \hfill dV\left(t\right)& =\left[\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right)\right]dt+{\delta }_{2}V\left(t\right)B_2\left(t\right),\hfill \\ \hfill dA\left(t\right)& =\left[{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +u_2\left(t\right))A\left(t\right)\right]dt+{\delta }_{3}A\left(t\right)B_3\left(t\right),\hfill \\ \hfill dZ\left(t\right)& =\left[{\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +u_3\left(t\right))Z\left(t\right)\right]dt+{\delta }_{4}Z\left(t\right)B_4\left(t\right),\hfill \\ \hfill dH\left(t\right)& =\left[{\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right)\right]dt+{\delta }_{5}H\left(t\right)B_5\left(t\right),\hfill \\ \hfill dR\left(t\right)& =\left[{\gamma }_{1}Z\left(t\right)+({\alpha }_2+u_2\left(t\right))A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)+u_1\left(t\right)S\left(t\right)\right]dt+{\delta }_{6}R\left(t\right)B_6\left(t\right).\hfill \end{array}$$

(55)

alongside initial conditions

$$\begin{array}{c}\hfill S\left(0\right)>0,V\left(0\right)\ge 0,A\left(0\right)\ge 0,Z\left(0\right)\ge 0,H\ge 0,R>0.\end{array}$$

(56)

To simplify things, a vector of the following form is introduced:

$$\begin{array}{c}\hfill x\left(t\right)={[x_1,x_2,x_3,x_4,x_5,x_6]}^{\prime },u\left(t\right)={[u_1,u_2,u_3]}^{\prime },\end{array}$$

(57)

Also,

$$\begin{array}{c}\hfill dx\left(t\right)=f(x,u)dt+g\left(x\right)dw\left(t\right).\end{array}$$

(58)

The temporal dependence of functions $u_i$ and $x_i$ is accompanied by the following expression for the initial data:

$$\begin{array}{c}\hfill x\left(0\right)={[x_1,x_2,x_3,x_4,x_5,x_6]}^{\prime }\left(0\right)=x_0.\end{array}$$

(59)

Below are the vectors that contain the functions f and g:

$$\begin{array}{cc}\hfill f_1(x\left(t\right),u\left(t\right))& =\left[\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi +u_1\left(t\right))S\left(t\right)\right]dt+{\delta }_{1}S\left(t\right)B_1\left(t\right),\hfill \\ \hfill f_2(x\left(t\right),u\left(t\right))& =\left[\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right)\right]dt+{\delta }_{2}V\left(t\right)B_2\left(t\right),\hfill \\ \hfill f_3(x\left(t\right),u\left(t\right))& =\left[{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +u_2\left(t\right))A\left(t\right)\right]dt+{\delta }_{3}A\left(t\right)B_3\left(t\right),\hfill \\ \hfill f_4(x\left(t\right),u\left(t\right))& =\left[{\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +u_3\left(t\right))Z\left(t\right)\right]dt+{\delta }_{4}Z\left(t\right)B_4\left(t\right),\hfill \\ \hfill f_5(x\left(t\right),u\left(t\right))& =\left[{\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right)\right]dt+{\delta }_{5}H\left(t\right)B_5\left(t\right),\hfill \\ \hfill f_6(x\left(t\right),u\left(t\right))& =\left[{\gamma }_{1}Z\left(t\right)+({\alpha }_2+u_2\left(t\right))A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)+u_1\left(t\right)S\left(t\right)\right]dt+{\delta }_{6}R\left(t\right)B_6\left(t\right).\hfill \end{array}$$

(60)

where $g_1={\delta }_1$, $g_2={\delta }_2$, $g_3={\delta }_3$, $g_4={\delta }_4$, $g_5={\delta }_5$, and $g_6={\delta }_6$. In this case, we a suppose quadratic cost function of the form

$$\begin{array}{c}\hfill \mathfrak{J}\left(u\right)={\displaystyle \frac{1}{2}}E\left\{{\int }_0^T\left(A_{1}A+A_{2}Z+{\displaystyle \frac{B_1{u}_1^2}{2}}+{\displaystyle \frac{B_2{u}_2^2}{2}}+{\displaystyle \frac{B_3{u}_3^2}{2}}\right)dt+{\displaystyle \frac{k_1{S}^2}{2}}+{\displaystyle \frac{k_1{V}^2}{2}}+{\displaystyle \frac{k_1{A}^2}{2}}+{\displaystyle \frac{k_1{Z}^2}{2}}+{\displaystyle \frac{k_1{H}^2}{2}}+{\displaystyle \frac{k_1{R}^2}{2}}\right\},\end{array}$$

(61)

where $A_1,A_2,B_1,B_2,B_3,$ and $k_1$ for $i=\{1,2,3,4,5,6\}$ are positive real numbers.

Again, we are focused on obtaining a control vector $u^{\ast }\left(t\right)=(u_1^{\ast }\left(t\right),u_2^{\ast }\left(t\right),u_3^{\ast }\left(t\right))$ as follows:

$$\begin{array}{c}\hfill \mathfrak{J}\left(u\right)\ge \mathfrak{J}\left(u^{\ast }\right),\forall u\in U,\end{array}$$

(62)

A controlling set (an admissible set) is U in this case.

$$\begin{array}{c}\hfill U=\left\{u_i\left(t\right);u_i\left(t\right)\in [0,u_i^{max}],\forall u_i\in L^2[0,T],t\in [0,T],i=1,2,3\right\}.\end{array}$$

(63)

In this context, $u_i^{max}$ for $i=1,2,3$ is assumed to be positive and real. The Hamiltonian $H_m(x,u,p,q)$ for the given system is to be formulated based on the stochastic maximum principle.

$$\begin{array}{c}\hfill H(x,u,p,q)=-l(x,u)+\langle g\left(x\right),q\rangle +\langle f(x,u),p\rangle .\end{array}$$

(64)

where $\langle ·,·\rangle$ denotes the Euclidean inner product space, and $q={[q_1,q_2,q_3,q_4,q_5,q_6]}^{\prime }$ and $p={[p_1,p_2,p_3,p_4,p_5,p_6]}^{\prime }$ are two distinct sets of adjoint vectors. The maximum principle can be expressed directly as

$$\begin{array}{c}\hfill d{x}^{\ast }\left(t\right)={\displaystyle \frac{\partial \mathbb{H}(x^{\ast },u^{\ast },p,q)}{\partial p}}dt+g\left(x^{\ast }\left(t\right)\right)dW\left(t\right),\end{array}$$

(65)

$$\begin{array}{c}\hfill d{p}^{\ast }\left(t\right)=-{\displaystyle \frac{\partial \mathbb{H}(x^{\ast },u^{\ast },p,q)}{\partial x}}dt+q\left(t\right)dW\left(t\right),\end{array}$$

(66)

$$\begin{array}{c}\hfill {\mathbb{H}}_m(x^{\ast },u^{\ast },p,q)=mi{n}_{u\in U}{\mathbb{H}}_m(x^{\ast },u^{\ast },p,q).\end{array}$$

(67)

Equations (65) and (66) have initial and terminal conditions for $x\left(t\right)$, where the optimal path is $x^{\ast }\left(t\right)$, and

$$\begin{array}{c}\hfill x^{\ast }\left(0\right)=x_0.\end{array}$$

(68)

$$\begin{array}{c}\hfill p\left(T\right)=-{\displaystyle \frac{\partial h\left(x^{\ast }\left(T\right)\right)}{\partial x}},\end{array}$$

(69)

accordingly. As shown by Equation (67), the $x^{\ast }\left(t\right)$ optimal control is a function of $q\left(t\right)$, $p\left(t\right)$, and $x^{\ast }\left(t\right)$; thus,

$$\begin{array}{c}\hfill d{x}^{\ast }\left(t\right)={\displaystyle \frac{\partial \mathbb{H}(x^{\ast },u^{\ast },p,q)}{\partial p}}dt+g\left(x^{\ast }\left(t\right)\right)dW\left(t\right),\end{array}$$

(70)

$$\begin{array}{c}\hfill dp\left(t\right)=-{\displaystyle \frac{\partial \mathbb{H}(x^{\ast },u^{\ast },p,q)}{\partial x}}dt+q\left(t\right)dW\left(t\right).\end{array}$$

(71)

Therefore, the given Hamiltonian is

$$\begin{array}{cc}\hfill \mathbb{H}& =\left(A_{1}A+A_{2}Z+{\displaystyle \frac{B_1{u}_1^2}{2}}+{\displaystyle \frac{B_2{u}_2^2}{2}}+{\displaystyle \frac{B_3{u}_3^2}{2}}+{\displaystyle \frac{k_1{S}^2}{2}},{\displaystyle \frac{k_2{V}^2}{2}},{\displaystyle \frac{k_3{A}^2}{2}},{\displaystyle \frac{k_4{Z}^2}{2}},{\displaystyle \frac{k_5{H}^2}{2}},{\displaystyle \frac{k_6{R}^2}{2}}\right)\hfill \\ & +P_1\left(\psi \Pi -{\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+\eta V\left(t\right)-(\gamma +\chi +u_1\left(t\right))S\left(t\right)\right)+{\delta }_{1}S{q}_1\hfill \\ & +P_2\left(\Pi (1-\psi )+\gamma S\left(t\right)-{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-(\eta +\chi )V\left(t\right)\right)+{\delta }_{2}V{q}_2\hfill \\ & +P_3\left({\displaystyle \frac{{\beta }_{1}S\left(t\right)A\left(t\right)}{N}}+{\displaystyle \frac{{\beta }_{2}V\left(t\right)A\left(t\right)}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi +u_2\left(t\right))A\left(t\right)\right)+{\delta }_{3}A{q}_3\hfill \\ & +P_4\left(t\right)\left({\alpha }_{1}A\left(t\right)-({\gamma }_1+{\psi }_2+{\gamma }_2+\chi +u_3\left(t\right))Z\left(t\right)\right)+{\delta }_{4}Z{q}_4\hfill \\ & +P_5\left(t\right)\left({\psi }_{1}A\left(t\right)+{\psi }_{2}Z\left(t\right)-({\psi }_3+{\psi }_4+\chi )H\left(t\right)\right)+{\delta }_{5}H{q}_5\hfill \\ & +P_6\left(t\right)\left(({\gamma }_1+u_3\left(t\right))Z\left(t\right)+({\alpha }_2+u_2\left(t\right))A\left(t\right)+{\psi }_{3}H\left(t\right)-\chi R\left(t\right)+u_1\left(t\right)S\left(t\right)\right)+{\delta }_{6}R{q}_6.\hfill \end{array}$$

(72)

Stochastic maximum theory states that

$$\begin{array}{c}\hfill d{p}^{\ast }\left(t\right)-{\displaystyle \frac{\partial \mathbb{H}(x^{\ast },u^{\ast },p,q)}{\partial x}}dt+q\left(t\right)dW\left(t\right).\end{array}$$

(73)

We obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{P}_1\left(t\right)}{dt}}& =(P_1\left(t\right)-P_3\left(t\right)){\displaystyle \frac{{\beta }_{1}A\left(t\right)}{N}}-u_1\left(t\right)P_6\left(t\right)+(\gamma +\chi +u_1\left(t\right))P_1\left(t\right)-\gamma P_2\left(t\right)+{\delta }_1{q}_1,\hfill \\ \hfill {\displaystyle \frac{d{P}_2\left(t\right)}{dt}}& =\eta (P_2\left(t\right)-P_1\left(t\right))+\chi P_2\left(t\right)+(P_2\left(t\right)-P_3\left(t\right)){\displaystyle \frac{{\beta }_{2}A\left(t\right)}{N}}+{\delta }_2{q}_2,\hfill \\ \hfill {\displaystyle \frac{d{P}_3}{dt}}Z& =-A_1+(P_1\left(t\right)-P_3\left(t\right)){\displaystyle \frac{{\beta }_{1}S\left(t\right)}{N}}+(P_2\left(t\right)-P_3\left(t\right)){\displaystyle \frac{{\beta }_{2}V\left(t\right)}{N}}+(P_3\left(t\right)-P_4\left(t\right)){\alpha }_1\hfill \\ & +(P_3\left(t\right)-P_6\left(t\right)){\alpha }_2+(P_3\left(t\right)-P_5\left(t\right)){\psi }_1+(\chi +u_2\left(t\right))P_3\left(t\right)+{\delta }_3{q}_3,\hfill \\ \hfill {\displaystyle \frac{d{P}_4\left(t\right)}{dt}}& =(P_4\left(t\right)-P_6\left(t\right)){\gamma }_1+(P_4\left(t\right)-P_5\left(t\right)){\psi }_2\left(t\right)+({\gamma }_2+\chi +u_2\left(t\right))P_4\left(t\right)+{\delta }_4{q}_4,\hfill \\ \hfill {\displaystyle \frac{d{P}_5\left(t\right)}{dt}}& =(P_5\left(t\right)-P_6\left(t\right)){\psi }_3+({\psi }_4+\chi )P_5\left(t\right)+{\delta }_5{q}_5,\hfill \\ \hfill {\displaystyle \frac{P_6\left(t\right)}{dt}}& =\chi P_6\left(t\right)+{\delta }_6{q}_6.\hfill \end{array}$$

(74)

Auxiliary starting and terminal conditions are also provided:

$$\begin{array}{c}\hfill S^{\ast }\left(0\right)=S^{\prime },V^{\ast }\left(0\right)=V^{\prime },A^{\ast }\left(0\right)=A^{\prime },Z^{\ast }\left(0\right)=Z^{\prime },H^{\ast }\left(0\right)=H^{\prime },R^{\ast }\left(0\right)=R^{\prime },p\left(T\right)=-{\displaystyle \frac{\partial h\left(x^{\ast }\left(T\right)\right)}{\partial x}},\end{array}$$

(75)

and

$$\begin{array}{c}\hfill h(S,V,A,Z,H,R)={\displaystyle \frac{k_1{S}^2}{2}}+{\displaystyle \frac{k_2{V}^2}{2}}+{\displaystyle \frac{k_3{A}^2}{2}}+{\displaystyle \frac{k_4{Z}^2}{2}}+{\displaystyle \frac{k_5{H}^2}{2}}+{\displaystyle \frac{k_6{R}^2}{2}},\end{array}$$

(76)

where $p_1\left(T\right)=-k_{1}S,$$p_2\left(T\right)=-k_{2}V,$$p_3\left(T\right)=-k_{3}A,$$p_4\left(T\right)=-k_{4}Z,$$p_5\left(T\right)=-k_{5}H,$ and $p_6\left(T\right)=-k_{1}R.$ Now, we differentiate the Hamiltonian equation relative to $u_1$ and $u_2$ to obtain the optimal controls: $u_1^{\ast }$ and $u_2^{\ast }$.

$$\begin{array}{c}\hfill u_1^{\ast }\left(t\right)=max\left\{0,min\{{\displaystyle \frac{[P_6\left(t\right)-P_1\left(t\right)]S^{\ast }}{B_2}},0\},1\right\},\end{array}$$

(77)

$$\begin{array}{c}\hfill u_2^{\ast }\left(t\right)=max\left\{0,min\{{\displaystyle \frac{[P_6\left(t\right)-P_3\left(t\right)]A^{\ast }}{B_3}},0\},1\right\},\end{array}$$

(78)

and

$$\begin{array}{c}\hfill u_3^{\ast }\left(t\right)=max\left\{0,min\{{\displaystyle \frac{[P_6\left(t\right)-P_4\left(t\right)]Z^{\ast }}{B_4}},0\},1\right\}.\end{array}$$

(79)

In control theory, the control variables are typically manipulated to achieve the desired objective(s). The objectives that are preplanned will be completed by governing the dynamics of the problem with differential equations and setting the control limits as indicated in Equation (60). To design a cost function, follow relation (61), which requires careful attention, especially when balancing terms and selecting words at $t^{\prime }s$ final value. Control theory’s primary tool (Pontryagin’s maximum principle [41]) must be demonstrated before being applied. In this step, the compactness argument is commonly used based on the control set and control boundedness. Once the boundedness of the cost function has been established, it becomes imperative to construct a minimizing/maximizing sequence of states and controls (depending on the nature of the problem). A crucial aspect of the analysis involves ensuring the convergence of such sequences in the relevant spaces. Consequently, the problem undergoes transformation into an equivalent optimization problem of the Hamiltonian on a point-by-point basis. Thus, the Hamiltonian can be defined as follows, drawing upon [42]:

$$\begin{array}{c}\hfill \mathrm{Hamiltonian}=(\mathrm{RHS}\mathrm{of}\mathrm{Differential}\mathrm{equations}\left)\right(\mathrm{adjoint})+(\mathrm{integrand}\mathrm{of}\mathrm{the}\mathrm{objective}\mathrm{functional}).\end{array}$$

To find the optimal solution, we must differentiate the Hamiltonian for $u^{\ast }$ at the point where it reaches its optimal value, also known as the optimality condition(s). We must take the derivative of the Hamiltonian for the state variables in order to obtain the deputy system. Transversality supports this process (75).

7. Computer Simulations

This study employed numerical simulations to validate the obtained analytical results. To obtain an estimate for the stochastic system, the initial step involved discretizing model (1) as shown below. This section presents the numerical simulations used to verify the obtained analytical results. As a first step, we will discretize model (1) as follows:

$$\begin{array}{cc}\hfill S_{i+1}& =S_i+\left[\psi \Pi -{\displaystyle \frac{{\beta }_1{S}_i{A}_i}{N}}+\eta V_i-(\gamma +\chi )S_i\right]\Delta t+{\delta }_1{S}_i\sqrt{\Delta t}{\zeta }_{1,i}+{\displaystyle \frac{{\delta }_1^2}{2}}{S}_i({\zeta }_{1,i}^2-1)\Delta t,\hfill \\ \hfill V_{i+1}& =V_i+\left[\Pi (1-\psi )+\gamma S_i-{\displaystyle \frac{{\beta }_2{V}_i{A}_i}{N}}-(\eta +\chi )V_i\right]\Delta t+{\delta }_2{V}_i\sqrt{\Delta t}{\zeta }_{2,i}+{\displaystyle \frac{{\delta }_2^2}{2}}{V}_i({\zeta }_{2,i}^2-1)\Delta t,\hfill \\ \hfill A_{i+1}& =A_i+\left[{\displaystyle \frac{{\beta }_1{S}_i{A}_i}{N}}+{\displaystyle \frac{{\beta }_2{V}_i{A}_i}{N}}-({\alpha }_1+{\alpha }_2+{\psi }_1+\chi )A_i\right]\Delta t+{\delta }_3{A}_i\sqrt{\Delta t}{\zeta }_{3,i}+{\displaystyle \frac{{\delta }_3^2}{2}}{A}_i({\zeta }_{3,i}^2-1)\Delta t,\hfill \\ \hfill Z_{i+1}& =Z_i+\left[{\alpha }_1{A}_i-({\gamma }_1+{\gamma }_2+{\psi }_2+\chi )Z_i\right]\Delta t+{\delta }_4{Z}_i\sqrt{\Delta t}{\zeta }_{4,i}+{\displaystyle \frac{{\delta }_4^2}{2}}{Z}_i({\zeta }_{4,i}^2-1)\Delta t,\hfill \\ \hfill H_{i+1}& =H_i+\left[{\psi }_1{A}_i+{\psi }_2{Z}_i-({\psi }_3+{\psi }_4+\chi )H_i\right]\Delta t+{\delta }_5{H}_i\sqrt{\Delta t}{\zeta }_{5,i}+{\displaystyle \frac{{\delta }_5^2}{2}}{H}_i({\zeta }_{5,i}^2-1)\Delta t,\hfill \\ \hfill R_{i+1}& =R_i+\left[{\gamma }_1{Z}_i+{\alpha }_2{A}_i+{\psi }_3{H}_i-\chi R_i\right]\Delta t+{\delta }_6{R}_i\sqrt{\Delta t}{\zeta }_{6,i}+{\displaystyle \frac{{\delta }_6^2}{2}}{R}_i({\zeta }_{6,i}^2-1)\Delta t.\hfill \end{array}$$

(80)

There are independent Gaussian stochastic variables ${\zeta }_{i,j}$ that follow a normal distribution, namely $N(0,1)$. In addition, ${\delta }_i$ (for $i=1,\cdots ,6$) represents the intensity levels of the white noise. Let us assume that time steps have a size of $\Delta t$. The stochastic model will be simulated using MATLAB using the algorithm outlined in (80). A deterministic model will be simulated using the conventional Runge–Kutta fourth-order (Rk4) algorithm. The following sections of the manuscript present the graphical outcomes of the simulations. The following section will present a graphical illustration of the model’s long-term behavior. Moreover, we will demonstrate how manipulating the values of specific sensitive parameters within the model can affect the disease dynamics.

7.1. An Extinction-Based Numerical Simulation

In this part, we are focused on simulating model (1) using the first-order Milstein stochastic method. Due to the involvement of six Brownian motions in the model, denoted by $d{\mathsf{B}}_{\mathsf{i}}\left(\mathsf{t}\right)$, the approach to approximate the double stochastic integral employs tools from stochastic theory.

To begin to demonstrate the disease’s extinction within the population, we will consider the following parameter values: $\Pi =0.03$, $\beta =0.002$, $\gamma =0.5$, $\eta =0.2$, $\chi =0.05$, ${\alpha }_1=0.01$, ${\alpha }_2=0.01$, $\psi =1$, ${\psi }_1=0.01$, ${\psi }_2=0.02$, ${\psi }_3=0.01$, ${\gamma }_1=0.025$, and ${\gamma }_2=0.01$. Likewise, the initial values for the population are taken as $({\mathsf{S}}_{\mathsf{0}},{\mathsf{V}}_{\mathsf{0}},{\mathsf{A}}_{\mathsf{0}},{\mathsf{Z}}_{\mathsf{0}},{\mathsf{H}}_{\mathsf{0}},{\mathsf{R}}_{\mathsf{0}})=(30,20,15,60,40,50)$, and the noise intensities are provided by $({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)=(20,15,5,20,10,10)$. We have set the time scale in our simulations to weeks. The basic reproduction number ($R_0^D$) plays a crucial role in infectious disease models, serving as a determinant of whether the disease can die out or will persist within the population. As a result, in our initial calculations, we determined that ${\mathsf{R}}_0^E=0.93$, which is evidently below one. Therefore, the conclusion of Theorem 2 is valid as evidenced by visual inspection of Figure 2a. The graph clearly illustrates that when ${\mathsf{R}}_0^E<1$ is maintained, the disease will eventually fade away from the population. The basic reproduction number $R_0^D$ for the deterministic model (Equation (1)) has been computed using the same parameter set. Once more, it was evident that $R_0^D=0.95$, which is unambiguously below unity. Drawing from the implications of Theorem 2, we arrive at a parallel outcome for the deterministic model. Alternatively, it can be stated that if ${\mathsf{R}}_0^E<1$, the disease-free equilibrium demonstrates both local and global asymptotic stability. It is crucial to demonstrate that, over time, the solutions of the deterministic model converge to the disease-free equilibrium.

7.2. Numerical Simulations of the Disease Distribution

The proposed stochastic model should have a unique stationary distribution; we will examine an alternative set of parameters [16]: $\Pi =0.03$, $\beta =0.002$, $\gamma =0.5$, $\eta =0.2$, $\chi =0.05$, ${\alpha }_1=0.001$, ${\alpha }_2=0.0002$, $\psi =0.01$, ${\psi }_1=1$, ${\psi }_2=0.35$, ${\psi }_3=0.001$, ${\gamma }_1=0.0025$, and ${\gamma }_2=0.01$. In this scenario, the white noise intensities are set as $({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)=(10,20,15,20,5,10)$, while the initial data for each compartment remain the same as in Section 7.2. In using these parameter values, the threshold parameter $R_0^D$ for the deterministic model was computed to be $1.75$. The value of $R_0^D$ plays a crucial role in assessing the stability of the model’s endemic equilibrium point $E_{\ast }$. If $R_0^D$ exceeds one, it can be inferred that $E_{\ast }$ is globally asymptotically stable. From the data presented in Figure 3a–f, it can be deduced that the infection will persist and spread throughout the population if ${\mathbb{R}}_0^S$ exceeds 1. Using the parameters in the model (1), we found that $R_0^D=1.7$, which is above the threshold of one. This observation indicates that the simulations are consistent with the conclusions outlined in Theorem 3. The curves show slight variations around their corresponding endemic equilibrium points in this scenario due to the small white noise intensities. This indicates that the disease will persist in the population for an extended period. Theorem 3 provides analytical evidence of the existence of an ergodic stationary distribution, further supported by the numerical results shown in Figure 4. Additionally, we created a frequency distribution histogram using the bootstrap-estimated values for each compartment, revealing the consistent persistence of the disease. Therefore, implementing more effective prevention and control measures can only eliminate the condition.

7.3. The Influence of $\beta$ on Stochastic Model (1)

We conducted simulations to show how $\beta$ affects the dynamics of the infected classes. Specifically, we looked at three values of $\beta$: 1, 1.2, and 1.5. Alongside other parameters outlined in Section 7.2, we used stochastic noises $({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)=(10,20,15,20,5,10).$ The results of the simulations are visually presented in Figure 5, where we compare the stochastic model and its deterministic counterpart. The graphs in Figure 5a–c show that as $\beta$ decreases, the populations converge more quickly toward disease extinction. Therefore, lowering the $\beta$ values is necessary to achieve disease extinction. A better understanding of HBV dynamics requires incorporating nonlinear stochastic noises.

7.4. The Influence of Noise $\delta$ on Stochastic System (1)

Let us change the values of $({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)$ while keeping the other parameters the same as in Section 7.2. Figure 6 depicts the partial solutions of the functions $A\left(t\right)$ and $Z\left(t\right)$ under the influence of this change. In comparing stochastic model (1) with its deterministic counterpart, significant perturbations in the stochastic version can lead to a population decrease, potentially reaching zero. However, these perturbations have only a minor impact on the dynamic behavior of infected individuals. According to our simulation analysis, even minor perturbations can have a notable impact on the magnitudes of $A\left(t\right)$ and $Z\left(t\right)$. The comprehensive examination of epidemic models, from both practical and theoretical viewpoints, necessitates the inclusion of nonlinear stochastic noises. Additionally, the numerical simulations, as shown in Figure 6a–c, emphasize the importance of implementing effective measures (such as quarantine strategies) to limit interactions between the environment and humans while regulating human activities. These steps are necessary to reduce the spread of HBV, especially in endemic areas.

7.5. Optimizing Control Strategies through Numerical Simulations

In this segment of the study, we showcase graphical outcomes of optimal control for models (1) and (2). The simulation was conducted utilizing the fourth-order ‘Runge–Kutta’ method.

To address state system (38) under condition (51) within the time interval $0,80$, we will initiate the solution process with the prescribed method. Subsequently, the adjoint equations will be approximated using the backward $RK-4$ method within the same time frame. The transversal conditions of (51) will also be considered.

The following parameters are used in the numerical plotting: $\Pi$ = 95, $\eta$ = 0.027, $\chi$ = 0.1, $\gamma$ = 0.09, ${\gamma }_1$ = 0.02, ${\gamma }_2$ = 0.003, $\beta$ = 0.4, ${\alpha }_1$ = 0.3, ${\alpha }_2$ = 0.1, $\psi$ = 0.07, ${\psi }_1$ = 0.03, ${\psi }_2$ = 0.3, ${\psi }_3$ = 0.1, and ${\psi }_4$ = 0.01. The values of $\delta$ are ${\delta }_1=0.943$, ${\delta }_2=0.405$, ${\delta }_3=0.345$, ${\delta }_4=0.601$, ${\delta }_5=0.601$, and ${\delta }_6=0.701$. The system initiates with the following initial conditions: $S\left(0\right)=70$, $V\left(0\right)=70$, $A\left(0\right)=100$, $Z\left(0\right)=20$, $H\left(0\right)=50$, and $R\left(0\right)=30$.

The results obtained are presented in Figure 7a–f, and Figure 8a–f showcase the control curve for model (2) and the dynamic curves for vulnerable, acutely, and chronically infected individuals. According to the simulation, implementing control measures decreased infected, susceptible, and hospitalized populations, while the recovered and vaccinated populations increased. The comparison between the two cases with and without controls shows significant differences.

To simulate model (1), the stochastic ‘RK4’ technique was used. The optimal solution was found by evaluating approximate adjoint equations and considering transversality conditions. The first step involved solving state system (74) using the RK4 scheme. Then, for the adjoint equations, we use a backward technique to elevate system (74). The conditions of transversality (75) were also taken into consideration.

The controlling agents were determined by combining the controls and the values from (77) and (78). This process was repeated until the unknown values were sufficiently close to those obtained previously. The plots depict a decline in the epidemic’s infectious compartments, which is concurrent with the expansion of the curative effect. Further details can be found in Figure 7 and Figure 8. Illustrations of optimal control schemes for types 1 and 2 are presented in Figure 7b–f and Figure 8b–f. In the context of the HBV model, these control types lead to reductions in the susceptible, acute infectious, chronic infectious, and hospitalized populations, while concurrently increasing the recovered and vaccinated populations.

In Figure 9a,b, both model (1) and (2) are shown with and without control factors. The difference in the graphs is clearly visible. This shows that the control factors have a significant effect on the spread of the disease.

8. Conclusions

This study provides a detailed stochastic optimal control analysis for an HBV epidemic model, with a particular focus on the impact of vaccination. Our work highlights several key contributions to the understanding and management of HBV transmission dynamics. Firstly, by incorporating stochastic differential equations into the model, we accounted for the inherent randomness and environmental fluctuations that significantly affect disease spread. This stochastic approach allowed us to capture the variability in transmission rates and the progression of HBV infections, providing a more realistic and comprehensive representation of the epidemic’s behavior compared to traditional deterministic models. Secondly, our analysis demonstrated the importance of vaccination as a primary control strategy. Through the application of optimal control theory, we derived strategies that not only minimize the number of infected individuals but also optimize the associated costs of vaccination. This dual focus on effectiveness and efficiency is crucial for public health policies, especially in resource-limited settings where maximizing the impact of available resources is imperative. The numerical simulations performed using the fourth-order Runge–Kutta method reinforced the theoretical findings. These simulations illustrated that optimal vaccination strategies could significantly reduce the prevalence of both acute and chronic HBV infections. The results underscore that timely and adequate vaccination can curb the spread of HBV, thereby lowering the incidence of severe liver diseases and related mortalities. Furthermore, our findings emphasize that the integration of stochastic elements into epidemic models enhances their predictive power and robustness. By considering random perturbations, we were able to better understand the potential range of epidemic outcomes and design control strategies that are resilient to uncertainties. This approach is particularly relevant in real-world scenarios where various unpredictable factors can influence disease dynamics. Our study also suggests several directions for future research. Extending the stochastic model to incorporate additional control measures, such as antiviral treatments and public health interventions like contact tracing and education campaigns, could provide a more holistic understanding of HBV management. Additionally, exploring the impact of demographic changes, migration, and other epidemiological factors could further refine the model and its applicability to diverse populations.

Limitation: Our model assumes a minimal impact from the chronic class on hospitalization rates, as NUC therapy has made hospitalization for most chronic patients unnecessary.