Stochastic Dynamics of Hepatitis B Virus Infection: Analysis, Stability, and Numerical Simulation
August 10, 2023
Abstract
This study develops and analyzes a stochastic differential equation (SDE) model for the dynamics of hepatitis B virus (HBV) infection. While deterministic frameworks have yielded important insights into viral behavior, they cannot adequately describe the intrinsic randomness and fluctuations present in biological processes. To address this limitation, we construct a stochastic model incorporating multiplicative environmental noise to account for variability in infection rates, cellular mortality, and viral replication. We establish a rigorous theoretical foundation by proving the existence, uniqueness, and global positivity of solutions for all biologically relevant initial conditions. Stability properties are investigated in detail, including stability in probability and almost sure exponential stability, with particular emphasis on conditions under which random perturbations stabilize the infection-free state. Furthermore, we demonstrate the existence of a unique ergodic stationary distribution and derive convergence properties of the uninfected hepatocyte population. Numerical simulations, performed via the Euler–Maruyama method with sufficiently small time step to ensure positivity and accuracy, validate the analytical results and illustrate the impact of stochastic fluctuations on system dynamics. The simulations confirm that environmental noise can induce viral extinction even in parameter regimes where deterministic analysis predicts persistence. These findings enhance the mathematical understanding of HBV infection dynamics and underscore the significant role of stochastic effects in shaping long-term disease outcomes.
Keywords: Stochastic differential equations, Hepatitis B virus, Almost sure exponential stability, Ergodic stationary distribution, Environmental noise, Viral dynamics
1 Introduction↩︎
Hepatitis B virus (HBV) infection remains a significant global health concern. According to recent estimates by the World Health Organization, approximately 296 million individuals are chronically infected, with 1.5 million new infections reported annually and over 820,000 deaths attributed each year, mainly due to cirrhosis and hepatocellular carcinoma [1]. Although safe and effective vaccines, along with potent antiviral regimens, are available, the biological processes governing HBV persistence, clearance, and treatment response remain highly complex. The intrinsic nonlinearities of viral replication and immune regulation, combined with patient heterogeneity, necessitate quantitative frameworks that can explain the observed variability in clinical outcomes.
Mathematical modeling has long provided such frameworks. The seminal contributions of Nowak and colleagues established the paradigm for within-host viral dynamics by coupling hepatocyte infection, viral replication, and immune responses in deterministic systems of ordinary differential equations (ODEs) [2]. Subsequent studies further refined these approaches to describe persistence, immune-mediated clearance, and therapeutic response in HBV and other viral infections [3]–[6]. For HBV, deterministic analyses linked observed multiphasic declines in viral load under therapy to parameter regimes that describe infection rates, viral production, and loss of infected hepatocytes [7]–[10]. These results, reinforced by the theory of threshold parameters such as the basic reproduction number [11], [12], have formed the backbone of mathematical virology.
However, deterministic models rest on the assumption that biological processes evolve with perfect regularity, an assumption that is rarely satisfied in practice. Viral replication, immune activity, and hepatocyte turnover are subject to stochastic variability at multiple scales. Fluctuations arise from the discrete nature of infection events, variability in host immune responses, and environmental perturbations. These stochastic effects can be decisive in determining whether infection is cleared, maintained at chronic levels, or reactivated. Lloyd [13] demonstrated that realistic infectious-period distributions alter persistence properties in epidemic systems, highlighting the need to account for randomness explicitly. At the cellular level, stochastic fluctuations in gene expression and viral production can create heterogeneity even in genetically identical populations [14]. In hepatitis C virus modeling, Guedj and Perelson [15] showed that stochastic analyses provided sharper estimates of viral half-life under therapy. Such findings underscore that the omission of randomness can lead to misleading conclusions about disease dynamics and therapeutic strategies.
Stochastic differential equations (SDEs) provide a natural formalism to capture these random fluctuations while retaining continuous-time dynamics. SDE-based epidemic models have revealed phenomena that are inaccessible to deterministic models, including noise-induced extinction, stochastic resonance, and altered thresholds for persistence [16]–[20]. For HBV, preliminary stochastic models have examined positivity and extinction under nonlinear incidence structures, demonstrating the analytical challenges and the potential of noise to reshape qualitative outcomes [21]. The mathematical theory of SDEs has developed powerful tools for addressing these questions, particularly Lyapunov techniques for establishing global existence, nonnegativity, and long-term stability [22]. Numerical methods, such as the Euler–Maruyama and Milstein schemes, extend these theoretical insights to simulation, although they require careful implementation to avoid artifacts, including spurious negativity [23], [24].
The present work develops a rigorous and comprehensive stochastic framework for HBV within-host dynamics. Building on earlier deterministic formulations [7]–[9], [25]–[27], we introduce multiplicative noise terms to represent variability in infection, viral production, and cellular death rates. We establish global existence and positivity of solutions, derive verifiable conditions for almost sure extinction of infected compartments, and prove convergence properties for the uninfected hepatocyte population. Long-term behavior is characterized through Foster–Lyapunov criteria, which guarantee the existence and uniqueness of invariant measures and ergodicity under biologically plausible conditions. In addition, we examine numerical approximations via Euler–Maruyama and Milstein methods, aligning theoretical stability conditions with computational practice.
Our findings demonstrate that stochastic perturbations can stabilize the infection-free state even when deterministic theory predicts persistence. This phenomenon of noise-induced extinction has profound implications, suggesting that intrinsic variability may act as a natural aid to immune clearance and that therapies exploiting stochastic effects could enhance viral control. By integrating rigorous analysis with biologically motivated modeling, this study provides both a methodological framework and novel insights into the stochastic nature of HBV infection.
2 Mathematical Model Formulation↩︎
2.1 Deterministic Foundation↩︎
We start with the deterministic framework of [25], which describes the interaction between hepatocytes and free virions through the nonlinear system of ordinary differential equations: \[\label{DM} \begin{align} \frac{dx}{dt} &= \Lambda - \mu_1 x - (1-\eta)\beta xz + qy, \\ \frac{dy}{dt} &= (1-\eta)\beta xz - \mu_2 y - qy, \\ \frac{dz}{dt} &= (1-\epsilon)py - \mu_3 z, \end{align}\tag{1}\] where \(x(t)\) denotes the population of uninfected hepatocytes, \(y(t)\) the population of infected hepatocytes, and \(z(t)\) the concentration of free viral particles. The biological interpretation of the parameters is given in Table 1.
| Parameter | Interpretation |
|---|---|
| \(\Lambda\) | Influx rate of uninfected hepatocytes |
| \(\mu_1\) | Death rate of uninfected hepatocytes |
| \(\mu_2\) | Death rate of infected hepatocytes |
| \(\mu_3\) | Clearance rate of free viral particles |
| \(\beta\) | Infection rate constant |
| \(p\) | Viral production rate per infected hepatocyte |
| \(q\) | Non-cytolytic cure rate of infected hepatocytes |
| \(\eta\) | Efficacy of treatment blocking infection \((0 \leq \eta \leq 1)\) |
| \(\epsilon\) | Efficacy of treatment blocking viral production \((0 \leq \epsilon \leq 1)\) |
The system 1 admits the infection–free equilibrium \[E_0 = \left(\frac{\Lambda}{\mu_1}, 0, 0\right),\] and, provided \(R_0>1\), a positive endemic equilibrium \(E_1=(x^*,y^*,z^*)\). The basic reproduction number, derived using the next–generation operator approach [11], takes the form \[\label{R0} R_0 = \frac{(1-\eta)(1-\epsilon)\beta p \Lambda}{(\mu_2+q)\mu_3 \mu_1}.\tag{2}\] The threshold \(R_0=1\) separates regimes of viral extinction and persistence in the deterministic setting.
2.2 Stochastic Extension↩︎
Deterministic equations capture only mean behavior. At the cellular and molecular scales, infection events, immune responses, and treatment effects fluctuate randomly. To incorporate such variability, we perturb the deterministic dynamics by an Itô noise of multiplicative type. The resulting system of stochastic differential equations is \[\label{SM} \begin{align} dx(t) &= \big[\Lambda - \mu_1 x - (1-\eta)\beta xz + qy\big]\,dt + \sigma_1 x \, dW_1(t), \\ dy(t) &= \big[(1-\eta)\beta xz - \mu_2 y - qy\big]\,dt + \sigma_2 y \, dW_2(t), \\ dz(t) &= \big[(1-\epsilon)py - \mu_3 z\big]\,dt + \sigma_3 z \, dW_3(t), \end{align}\tag{3}\] where \((W_1(t),W_2(t),W_3(t))\) are independent standard Brownian motions on a filtered probability space \((\Omega,\mathcal{F},\{\mathcal{F}_t\}_{t\geq 0},\mathbb{P})\), and \(\sigma_i>0\) \((i=1,2,3)\) denote the intensities of the random fluctuations.
The choice of multiplicative noise reflects the empirical observation that larger populations exhibit larger absolute fluctuations. Meanwhile, the vanishing of noise at the boundary \(x=0\), \(y=0\), \(z=0\) ensures biological feasibility by preventing the emergence of negative solutions. Since the drift coefficients in 3 are locally Lipschitz and satisfy linear growth bounds, the classical existence–uniqueness theorem for SDEs ensures the well–posedness of strong solutions [16], [22].
2.3 Vector Formulation↩︎
For analytical convenience, system 3 can be written in vector form: \[\label{SM1} du(t) = f(u(t))\,dt + B(u(t))\,dW(t), \qquad u(0)=u_0 \in \mathbb{R}^3_{+},\tag{4}\] with \[u(t) = \begin{pmatrix}x(t) \\ y(t) \\ z(t)\end{pmatrix}, \quad dW(t) = \begin{pmatrix}dW_1(t) \\ dW_2(t) \\ dW_3(t)\end{pmatrix},\] and \[f(u) = \begin{pmatrix} \Lambda - \mu_1 x - (1-\eta)\beta xz + qy \\ (1-\eta)\beta xz - \mu_2 y - qy \\ (1-\epsilon)py - \mu_3 z \end{pmatrix}, \qquad B(u) = \begin{pmatrix} \sigma_1 x & 0 & 0 \\ 0 & \sigma_2 y & 0 \\ 0 & 0 & \sigma_3 z \end{pmatrix}.\] This compact representation allows the direct use of Itô calculus and Lyapunov techniques in the subsequent analysis of stability, extinction, and ergodicity.
3 Mathematical Preliminaries↩︎
This section establishes the analytic framework required for the subsequent results. We recall standard facts from the theory of stochastic differential equations (SDEs), Itô calculus, and auxiliary inequalities that will be used in the stability and ergodicity analysis.
3.1 Stochastic Differential Equations↩︎
Let \((\Omega,\mathcal{F},\{\mathcal{F}_t\}_{t\geq t_0},\mathbb{P})\) be a complete filtered probability space satisfying the usual conditions, and let \(W(t) \in \mathbb{R}^m\) denote an \(m\)-dimensional standard Brownian motion adapted to \(\{\mathcal{F}_t\}\). We consider a \(d\)-dimensional Itô SDE of the form \[\label{gSM1} du(t) = f(u(t),t)\,dt + B(u(t),t)\,dW(t), \qquad t \geq t_0,\tag{5}\] where the drift \(f:\mathbb{R}^d \times [t_0,T]\to \mathbb{R}^d\) and the diffusion coefficient \(B:\mathbb{R}^d \times [t_0,T]\to \mathbb{R}^{d\times m}\) are Borel measurable.
Definition 1 (Strong solution of an SDE). A stochastic process \(\{u(t)\}_{t_0 \leq t \leq T}\) is called a strong solution of 5 if:
\(u(t)\) is continuous, \(\mathcal{F}_t\)-adapted, and takes values in \(\mathbb{R}^d\);
\(f(u(t),t) \in L^1([t_0,T];\mathbb{R}^d)\) almost surely and \(B(u(t),t) \in L^2([t_0,T];\mathbb{R}^{d\times m})\) almost surely;
For all \(t\in [t_0,T]\), \[\label{eq:integral95form} u(t) = u(t_0) + \int_{t_0}^t f(u(s),s)\,ds + \int_{t_0}^t B(u(s),s)\,dW(s) \quad \text{a.s.}\qquad{(1)}\]
The solution is unique if whenever \(\bar{u}(t)\) is another solution of 5 , one has \[\mathbb{P}\{u(t)=\bar{u}(t)\;\;\forall\,t\in[t_0,T]\}=1.\]
The classical existence and uniqueness theorem states that if \(f\) and \(B\) satisfy the local Lipschitz condition and the linear growth condition, then 5 admits a unique strong solution on \([t_0,T]\) (see [16], [22]).
3.2 Itô’s Formula and Infinitesimal Generator↩︎
Let \(u(t)\) satisfy 5 , and let \(V:\mathbb{R}^d\to \mathbb{R}\) be twice continuously differentiable. By Itô’s formula, \[\label{ito95formula} dV(u(t)) = \mathcal{L}V(u(t))\,dt + \nabla V(u(t)) \cdot B(u(t),t)\,dW(t),\tag{6}\] where the infinitesimal generator \(\mathcal{L}\) is given by \[\label{generator} \mathcal{L}V(u) = \sum_{i=1}^d f_i(u)\,\frac{\partial V}{\partial u_i}(u) + \tfrac{1}{2}\sum_{i,j=1}^d \big[B(u)B^T(u)\big]_{ij}\,\frac{\partial^2 V}{\partial u_i \partial u_j}(u).\tag{7}\] This operator plays a central role in establishing stability, extinction, and ergodic properties of stochastic epidemic models.
3.3 Auxiliary Inequality↩︎
The following inequality will be used repeatedly in Lyapunov-based estimates.
Lemma 1. For all \(v>0\), one has \[v \;\leq\; 2\,(v+1-\ln v) - (4-2\ln 2),\] With equality if and only if \(v=2\).
Proof. Define \[g(v) = 2(v+1-\ln v) - v - (4-2\ln 2).\] Then \[g'(v) = 1 - \frac{2}{v}, \qquad g''(v) = \frac{2}{v^2} > 0.\] Hence \(g\) is strictly convex on \((0,\infty)\). The critical point occurs when \(g'(v)=0\), i.e. \(v=2\). Since \(g''(v)>0\), this is the global minimizer. Evaluating gives \(g(2)=0\), so \(g(v)\geq 0\) for all \(v>0\), with equality only at \(v=2\). The inequality follows immediately. ◻
4 Existence, Uniqueness, and Positivity of Solutions↩︎
The first step in analyzing the stochastic system 3 is to establish that solutions exist for all time, are unique, and remain positive. These properties are essential in biological modeling, since negative states have no physical meaning.
Theorem 1 (Global Existence, Uniqueness, and Positivity). Let \(u_0=(x_0,y_0,z_0)\in \mathbb{R}^3_{+}\). Then the stochastic system 3 admits a unique global strong solution \(u(t)=(x(t),y(t),z(t))\) defined for all \(t\geq 0\) almost surely, and the solution satisfies \(u(t)\in \mathbb{R}^3_{+}\) for all \(t\geq 0\) almost surely.
Proof. Step 1. Local existence and uniqueness. The drift term \(f(u)\) and diffusion coefficient \(B(u)\) of system 3 are locally Lipschitz on \(\mathbb{R}^3_{+}\) and satisfy a linear growth condition. By the standard existence–uniqueness theorem for SDEs (see [16], [22]), there exists a unique local strong solution \(u(t)\) defined on \([0,\tau_e)\), where \(\tau_e\) is the explosion time.
Step 2. Construction of stopping times. To show that \(\tau_e=\infty\) almost surely, define the sequence of stopping times \[\tau_n = \inf\big\{t\geq 0: \min\{x(t),y(t),z(t)\}\leq \tfrac{1}{n} \;\; \text{or} \;\; \max\{x(t),y(t),z(t)\}\geq n \big\},\] with \(\tau_n\) increasing in \(n\) and \(\tau_\infty = \lim_{n\to\infty}\tau_n\). Clearly \(\tau_n\leq \tau_e\).
Step 3. Lyapunov function. Consider the function \[V(u) = \big(x+1-\ln x\big) + \big(y+1-\ln y\big) + \big(z+1-\ln z\big).\] This function is \(C^2\) on \(\mathbb{R}^3_{+}\), satisfies \(V(u)\geq 3\), and \(V(u)\to \infty\) as any component tends to \(0\) or \(\infty\).
Step 4. Application of Itô’s formula. Applying Itô’s formula to \(V(u(t))\) for \(t<\tau_n\), we obtain \[dV(u(t)) = \mathcal{L}V(u(t))\,dt + dM(t),\] where \(M(t)\) is a local martingale and \[\mathcal{L}V(u) = \sum_{i=1}^3 f_i(u)\frac{\partial V}{\partial u_i}(u) + \tfrac{1}{2}\sum_{i=1}^3 \sigma_i^2 u_i^2 \frac{\partial^2 V}{\partial u_i^2}(u).\]
Since \(\partial V/\partial x = 1 - \tfrac{1}{x}\), \(\partial^2 V/\partial x^2 = \tfrac{1}{x^2}\), and similarly for \(y\) and \(z\), direct substitution yields \[\begin{align} \mathcal{L}V(u) &= \left(1-\tfrac{1}{x}\right)(\Lambda - \mu_1 x - (1-\eta)\beta xz + qy) \\ &\quad + \left(1-\tfrac{1}{y}\right)((1-\eta)\beta xz - \mu_2 y - qy) \\ &\quad + \left(1-\tfrac{1}{z}\right)((1-\epsilon)py - \mu_3 z) \\ &\quad + \tfrac{1}{2}\big(\sigma_1^2 + \sigma_2^2 + \sigma_3^2\big). \end{align}\]
Step 5. Bounding the generator. Rearranging and using the inequality of Lemma 1, one finds constants \(K_1,K_2>0\) such that \[\mathcal{L}V(u) \leq K_1 + K_2 V(u), \qquad u\in \mathbb{R}^3_{+}.\]
Step 6. Expectation estimates. Taking expectation and applying Dynkin’s formula gives \[\mathbb{E}[V(u(t\wedge \tau_n))] \leq V(u_0) + K_1 t + K_2\int_0^t \mathbb{E}[V(u(s\wedge \tau_n))]\,ds.\] By Gronwall’s inequality, \[\label{boundV} \mathbb{E}[V(u(t\wedge \tau_n))] \leq \big(V(u_0) + K_1 t\big)e^{K_2 t}, \qquad t\geq 0.\tag{8}\] Thus \(\sup_{n}\mathbb{E}[V(u(t\wedge \tau_n))]<\infty\) for each fixed \(t\).
Step 7. Non-explosion. Suppose, for contradiction, that \(\mathbb{P}(\tau_\infty\leq T)>\varepsilon\) for some \(T>0\) and \(\varepsilon>0\). Then on \(\{\tau_n\leq T\}\) one has \[V(u(\tau_n)) \geq \min\big\{n+1-\ln n, \tfrac{1}{n}+1+\ln n\big\}\to \infty \quad \text{as } n\to\infty.\] This contradicts the uniform bound 8 . Hence \(\tau_\infty=\infty\) almost surely, implying \(\tau_e=\infty\) almost surely. Therefore, the solution exists globally.
Step 8. Positivity. Positivity follows from two facts: (i) the Lyapunov function \(V(u)\) blows up as any coordinate approaches \(0\), so the solution cannot reach the boundary in finite time; (ii) the diffusion coefficients vanish at \(x=0\), \(y=0\), or \(z=0\), so no component can cross into the negative axis. Thus \(u(t)\in\mathbb{R}^3_{+}\) for all \(t\geq 0\) almost surely.
This completes the proof. ◻
5 Stability Analysis↩︎
The long-term dynamics of stochastic epidemic systems require a more refined analysis than their deterministic counterparts. Random attractors or stationary distributions replace deterministic equilibria, and stability concepts must account for probabilistic fluctuations. We therefore begin by recalling the precise definitions of stochastic stability and then proceed to establish extinction conditions for the infected compartments, together with the limiting behavior of the uninfected hepatocytes.
5.1 Stability Concepts↩︎
Definition 2 (Stochastic Stability [16], [22]). Let \(u(t)\) be the solution of an SDE 5 with equilibrium point \(u^*=0\).
The equilibrium is stable in probability* if for every \(\varepsilon\in(0,1)\) and \(r>0\), there exists \(\delta>0\) such that \[|u_0|<\delta \quad \Longrightarrow \quad \mathbb{P}\{|u(t)|<r \text{ for all } t\geq 0\} \geq 1-\varepsilon.\]*
The equilibrium is asymptotically stable in probability* if it is stable in probability and \[\lim_{t\to\infty} \mathbb{P}\{u(t)\to 0\}=1 \qquad \text{for all sufficiently small } |u_0|.\]*
The equilibrium is almost surely exponentially stable* if \[\limsup_{t\to\infty} \frac{1}{t}\ln |u(t)| < 0 \qquad \text{a.s.}\]*
5.2 Extinction of Infected Compartments↩︎
The compartments \(y(t)\) (infected hepatocytes) and \(z(t)\) (free virions) do not admit positive equilibria when subject to strong stochastic fluctuations. Their asymptotic decay can be established by Lyapunov analysis.
Theorem 2 (Exponential Extinction of \(y(t)\) and \(z(t)\)). Let \(u(t)=(x(t),y(t),z(t))\) satisfy 3 with \(u_0\in\mathbb{R}^3_{+}\). Suppose the following conditions hold: \[\begin{align} &(1-\epsilon)p - q - \mu_2 + \tfrac{1}{2}\sigma_2^2 < 0, \label{stab95cond1}\\ &\Delta_{\text{stab}} \leq 0, \label{stab95cond2} \end{align}\] {#eq: sublabel=eq:stab95cond1,eq:stab95cond2} where \[\begin{align} \Delta_{\text{stab}} &= \big[(1-\eta)\beta \bar{x} - \mu_3\big]\big[(1-\epsilon)p - q - \mu_2\big] \\ &\quad - \big[(1-\epsilon)p - q - \mu_2 + \tfrac{1}{2}\sigma_2^2\big]\big[(1-\eta)\beta \bar{x} - \mu_3 + \tfrac{1}{2}\sigma_3^2\big], \end{align}\] with \(\bar{x}=\Lambda/\mu_1\). Then \(y(t)\) and \(z(t)\) satisfy \[\limsup_{t\to\infty} \frac{1}{t}\ln(y(t)+z(t)) < 0 \quad \text{a.s.},\] i.e., both infected compartments decay to zero exponentially fast almost surely.
Proof. Define the Lyapunov function \[V(y,z)=\ln(y+z), \qquad (y,z)\in\mathbb{R}^2_{+}.\] Applying Itô’s formula gives \[\begin{align} dV &= \frac{1}{y+z}\Big( [(1-\epsilon)p - q - \mu_2]y + [(1-\eta)\beta x - \mu_3]z \Big)\,dt \\ &\quad + \frac{1}{2}\frac{\sigma_2^2 y^2 + \sigma_3^2 z^2}{(y+z)^2}\,dt + \frac{\sigma_2 y}{y+z}\,dW_2 + \frac{\sigma_3 z}{y+z}\,dW_3. \end{align}\]
Introduce the vector \(v=(y,z)^T\) and observe that \[\mathcal{L}V = \frac{v^T A v}{(y+z)^2},\] where \[A = \begin{pmatrix} (1-\epsilon)p - q - \mu_2 + \tfrac{1}{2}\sigma_2^2 & (1-\eta)\beta \bar{x} - \mu_3 \\ (1-\epsilon)p - q - \mu_2 & (1-\eta)\beta \bar{x} - \mu_3 + \tfrac{1}{2}\sigma_3^2 \end{pmatrix}.\]
Note that \(x(t)\) remains bounded by \(\bar{x}+\delta\) for arbitrarily small \(\delta>0\) almost surely, since \(x(t)\) converges in probability to the linear process ?? (see Theorem 3). Thus, we can replace \(x(t)\) by \(\bar{x}\) in the asymptotic analysis.
By conditions ?? –?? , the matrix \(A\) is negative definite. Let \(\lambda_{\max}<0\) denote its largest eigenvalue. Then \[\mathcal{L}V \leq \frac{\lambda_{\max}}{(y+z)^2}(y^2+z^2).\] Using \((y^2+z^2)/(y+z)^2 \geq \tfrac{1}{2}\), we obtain \[\mathcal{L}V \leq \tfrac{1}{2}\lambda_{\max}.\]
Thus \[V(y(t),z(t)) \leq V(y(0),z(0)) + \tfrac{1}{2}\lambda_{\max} t + M(t),\] where \(M(t)\) is a martingale with \(\lim_{t\to\infty}M(t)/t=0\) almost surely. Dividing by \(t\) and letting \(t\to\infty\) yields \[\limsup_{t\to\infty}\frac{1}{t}\ln(y(t)+z(t)) \leq \tfrac{1}{2}\lambda_{\max} < 0.\] Hence, both \(y(t)\) and \(z(t)\) converge to zero exponentially fast almost surely. ◻
Remark 1. The negative definiteness of \(A\) reflects the combined damping effect of biological death rates and environmental noise. The inequalities ?? –?? demonstrate that noise intensities \(\sigma_2,\sigma_3\) can shift the system from persistence to extinction, even in regimes where the deterministic basic reproduction number \(R_0>1\).
5.3 Dynamics of Uninfected Hepatocytes↩︎
The uninfected hepatocyte population \(x(t)\) does not approach zero but fluctuates around the deterministic equilibrium \(\Lambda/\mu_1\). To capture this behavior, we study an auxiliary linear SDE.
Lemma 2 (Analysis of Linear Reference Process). Consider \[\label{x1-eq} dx_1(t) = \big(\Lambda - \mu_1 x_1(t)\big)\,dt + \sigma_1 x_1(t)\,dW_1(t), \qquad x_1(0)>0.\qquad{(2)}\]
The explicit solution is \[\begin{align} x_1(t) &= x_1(0)\exp\!\left(-\Big(\mu_1+\tfrac{1}{2}\sigma_1^2\Big)t + \sigma_1 W_1(t)\right) \\ &\quad + \Lambda\int_0^t \exp\!\left(-\Big(\mu_1+\tfrac{1}{2}\sigma_1^2\Big)(t-s) + \sigma_1(W_1(t)-W_1(s))\right)\,ds. \end{align}\]
Its expectation satisfies \[\mathbb{E}[x_1(t)] = x_1(0)e^{-\mu_1 t} + \frac{\Lambda}{\mu_1}(1-e^{-\mu_1 t}),\] so that \(\lim_{t\to\infty}\mathbb{E}[x_1(t)]=\Lambda/\mu_1\).
(iii) If \(2\mu_1>\sigma_1^2\), the process \(\{x_1(t)\}\) admits a unique invariant probability measure \(\pi\) on \((0,\infty)\) with mean \(\mathbb{E}_\pi[x_1]=\Lambda/\mu_1\). .
Proof. Part (i) follows from the integrating factor method for linear SDEs. For (ii), take expectations in (i). Using the independence of increments and \[\mathbb{E}\!\left[\exp\big(\sigma_1(W_1(t)-W_1(s))\big)\right] = \exp\!\Big(\tfrac{1}{2}\sigma_1^2(t-s)\Big),\] the integral simplifies to the deterministic term \(\tfrac{\Lambda}{\mu_1}(1-e^{-\mu_1 t})\). Finally, (iii) is a standard result: linear mean-reverting SDEs with multiplicative noise are ergodic with stationary Gamma-type distributions (see [22]). The mean of this distribution coincides with the deterministic equilibrium \(\Lambda/\mu_1\). ◻
Theorem 3 (Asymptotics of \(x(t)\)). Let \((x(t),y(t),z(t))\) satisfy 3 and let \(x_1(t)\) satisfy ?? . Under the assumptions of Theorem 2, \[\lim_{t\to\infty}[x(t)-x_1(t)] = 0 \quad \text{in probability}.\]
Proof. Define for \(r>0\) the auxiliary process \[dx_r(t) = \big(\Lambda-(\mu_1+r)x_r(t)\big)\,dt + \sigma_1 x_r(t)\,dW_1(t).\] Comparison theorems for SDEs imply \(x_r(t)\leq x(t)\leq x_1(t)\). Subtracting gives \[\begin{align} x(t)-x_r(t) &= \phi(t)\int_0^t \phi(s)^{-1}\Big(rx_r(s) - (1-\eta)\beta x(s)z(s) + qy(s)\Big)\,ds, \end{align}\] where \[\phi(t) = \exp\!\left(-\Big(\mu_1+r+\tfrac{1}{2}\sigma_1^2\Big)t + \sigma_1 W_1(t)\right).\]
By continuity of eigenvalues with respect to parameters and the boundedness in probability of \(x(t)\), there exists \(\delta>0\) such that the principal minors remain negative for \(x\in[\bar{x}-\delta,\bar{x}+\delta]\), ensuring \(A\) is negative definite along almost all sample paths for large \(t\).
By Theorem 2, \(y(t),z(t)\to 0\) exponentially almost surely. Consequently, the integral term converges to \(0\) as \(t\to\infty\), and since \(r>0\) is arbitrary, the sandwich estimate \(x_r(t)\leq x(t)\leq x_1(t)\) yields the desired convergence in probability. ◻
6 Ergodic Stationary Distributions↩︎
We establish the existence and uniqueness of an invariant probability measure and the ergodic property for the HBV SDE 3 . The argument proceeds via a Foster–Lyapunov drift condition together with a local nondegeneracy condition on the diffusion, following standard criteria for positive recurrence and ergodicity of continuous-time Markov processes.
6.1 Ergodicity Criterion↩︎
Let \(U(t)\) be a time-homogeneous Markov process in \(\mathbb{R}^d\) governed by the Itô SDE \[dU(t)=f(U(t))\,dt+\sum_{k=1}^m B_k(U(t))\,dW_k(t),\] with generator \[\mathcal{L}g(u)=\sum_{i=1}^d f_i(u)\,\partial_{u_i} g(u)+\tfrac12\sum_{i,j=1}^d a_{ij}(u)\,\partial^2_{u_i u_j} g(u),\qquad A(u)=[a_{ij}(u)]=\sum_{k=1}^m B_k(u)B_k(u)^{\top}.\]
Lemma 3 (Positive recurrence and ergodicity [18], [19]). Suppose there exists a bounded open set \(D\subset\mathbb{R}^d\) and a function \(V\in C^2(\mathbb{R}^d)\) with \(V\ge 1\), \(V(u)\to\infty\) as \(|u|\to\infty\), such that \[\begin{align} &\text{\emph{(A1) local nondegeneracy:}}\quad \exists\,\underline{\lambda}>0\;\text{with}\; \xi^{\top}A(u)\,\xi\ge \underline{\lambda}|\xi|^2,\quad \forall\,u\in D,\;\forall\,\xi\in\mathbb{R}^d; \label{A1}\\ &\text{\emph{(A2) Foster--Lyapunov drift:}}\quad \exists\,\theta>0,\;C<\infty\;\text{with}\; \mathcal{L}V(u)\le -\theta \;\;\text{for all } u\in D^{c}. \label{A2} \end{align}\] {#eq: sublabel=eq:A1,eq:A2} Then \(U(t)\) admits a unique invariant probability measure \(\pi\) and is positive recurrent. Moreover, for any \(g\) with \(\int |g|\,d\pi<\infty\), \[\lim_{T\to\infty}\frac{1}{T}\int_0^T g\big(U(t)\big)\,dt=\int_{\mathbb{R}^d} g(u)\,\pi(du)\quad\text{almost surely.}\]
6.2 Application to the HBV SDE↩︎
We apply Lemma 3 to \(U(t)=(x(t),y(t),z(t))\) solving 3 . Throughout, write \(\bar{x}=\Lambda/\mu_1\).
Theorem 4 (Existence and uniqueness of an ergodic stationary distribution). Assume the extinction conditions of Theorem 2 hold. Then the Markov process \(U(t)=(x(t),y(t),z(t))\) generated by 3 admits a unique invariant probability measure \(\pi\) on \(\mathbb{R}^3_{+}\) and is ergodic in the sense of Lemma 3. In particular, \[\lim_{T\to\infty}\frac{1}{T}\int_0^T g\big(x(t),y(t),z(t)\big)\,dt=\int_{\mathbb{R}^3_{+}} g(u)\,\pi(du)\quad\text{a.s.}\] for any \(\pi\)-integrable \(g\).
Proof. We verify ?? –?? .
Step 1: Local nondegeneracy on a bounded set. The diffusion matrix of 3 is \[A(x,y,z)=\operatorname{diag}\big(\sigma_1^2 x^2,\;\sigma_2^2 y^2,\;\sigma_3^2 z^2\big).\] Fix \(0<\rho<R<\infty\) and set the bounded open set \[D:=\big\{(x,y,z)\in\mathbb{R}^3_{+}:\;\rho<x<R,\;\rho<y<R,\;\rho<z<R\big\}.\] For \(u=(x,y,z)\in D\) and any \(\xi\in\mathbb{R}^3\), \[\xi^{\top}A(u)\xi = \sigma_1^2 x^2 \xi_1^2+\sigma_2^2 y^2 \xi_2^2+\sigma_3^2 z^2 \xi_3^2 \;\ge\;\underline{\lambda}|\xi|^2,\quad \underline{\lambda}:=\rho^2 \min\{\sigma_1^2,\sigma_2^2,\sigma_3^2\}>0.\] Hence ?? holds.
Step 2: Foster–Lyapunov drift outside \(D\). Consider the coercive Lyapunov function \[V(x,y,z)= (x-\bar{x})^2 + c_2\,y + c_3\,z,\] with positive constants \(c_2,c_3\) to be chosen. Note \(V\in C^2\), \(V\ge 1\) on \(\mathbb{R}^3_{+}\) (after increasing by \(1\) if needed), and \(V(u)\to\infty\) as \(|u|\to\infty\).
Compute \(\mathcal{L}V\) using the generator of 3 . Since \(V_{xx}=2\), \(V_{yy}=V_{zz}=0\), we obtain \[\begin{align} \mathcal{L}V &= 2(x-\bar{x})\big(\Lambda-\mu_1 x-(1-\eta)\beta x z + q y\big) + c_2\big((1-\eta)\beta x z-(\mu_2+q)y\big) \\ &\quad + c_3\big((1-\epsilon)p y - \mu_3 z\big) \;+\;\sigma_1^2 x^2. \label{eq:LVstar} \end{align}\tag{9}\] Split terms and bound each contribution.
(i) \(x\)–drift and \(x\)–diffusion: \[2(x-\bar{x})(\Lambda-\mu_1 x)=2(x-\bar{x})(\mu_1\bar{x}-\mu_1 x)=-2\mu_1(x-\bar{x})^2.\] Moreover \(x^2\le 2(x-\bar{x})^2+2\bar{x}^2\), hence \[\sigma_1^2 x^2\le 2\sigma_1^2 (x-\bar{x})^2 + 2\sigma_1^2 \bar{x}^2.\] Combining, \[2(x-\bar{x})(\Lambda-\mu_1 x)+\sigma_1^2 x^2 \le -(2\mu_1-2\sigma_1^2)(x-\bar{x})^2 + 2\sigma_1^2\bar{x}^2.\]
(ii) \(xz\)–coupling terms: \[-2(1-\eta)\beta (x-\bar{x})\,x z + c_2(1-\eta)\beta x z.\] For any \(\delta>0\), \[-2(1-\eta)\beta (x-\bar{x})\,x z \le \delta (x-\bar{x})^2 + \frac{(1-\eta)^2\beta^2}{\delta} x^2 z^2.\] We control \(x^2 z^2\) outside \(D\) by the linear terms \(y\) and \(z\) via the drift of \(y,z\); this is handled below using Young’s inequalities and the negativity of linear terms \(-(\mu_2+q)c_2 y\) and \(-\mu_3 c_3 z\). Precisely, for any \(\kappa_1,\kappa_2>0\), \[x^2 z^2 \le \kappa_1 z + C_1(\kappa_1)\,x^4 z^3,\qquad x^2 z^2 \le \kappa_2 y + C_2(\kappa_2)\,x^4 z^2 y^{-1},\] and since outside \(D\) at least one coordinate is either \(<\rho\) or \(>R\), one can choose the domain \(D\) large and \(\rho\) small enough so that the linear dissipative terms in \(y\) and \(z\) dominate the above remainders. A simpler and classical way that avoids high-order remainders is to choose a linear Lyapunov weight in \(z\) and absorb \(xz\) by \(z\): \[c_2(1-\eta)\beta x z \le \frac{\mu_3 c_3}{2} z + \frac{c_2^2 (1-\eta)^2 \beta^2}{2\mu_3 c_3} x^2.\] Using again \(x^2\le 2(x-\bar{x})^2+2\bar{x}^2\), we obtain \[\begin{align} &-2(1-\eta)\beta (x-\bar{x})\,x z + c_2(1-\eta)\beta x z \\ &\qquad \le \delta (x-\bar{x})^2 + \frac{c_2^2 (1-\eta)^2 \beta^2}{\mu_3 c_3}\,(x-\bar{x})^2 + \frac{c_2^2 (1-\eta)^2 \beta^2}{\mu_3 c_3}\,\bar{x}^2 + \frac{\mu_3 c_3}{2} z. \end{align}\]
(iii) \(y\)– and \(z\)–linear drifts and the \(qy\) coupling: \[2(x-\bar{x})\,q y - c_2(\mu_2+q) y + c_3(1-\epsilon)p\, y - \mu_3 c_3 z.\] Bound \(2(x-\bar{x})\,q y\) via Young’s inequality with constants \(\varepsilon_y>0\): \[2q |x-\bar{x}|\, y \le \varepsilon_y (x-\bar{x})^2 + \frac{q^2}{\varepsilon_y}\, y^2.\] Since \(V\) is linear in \(y\) and \(z\), we use the fact that outside \(D\) either \(y\) or \(z\) is large. To maintain a linear Foster drift, we apply \(y^2\le K_y y + K'_y\) outside \(D\) for suitable \(K_y,K'_y\) depending on \((\rho,R)\); thus the \(y^2\) term can be absorbed into \(-c_2(\mu_2+q)y\) by choosing \(c_2\) sufficiently large. Moreover, choose \(c_2,c_3\) so that \[(\mu_2+q)c_2 > c_3(1-\epsilon)p \quad \text{and} \quad \mu_3 c_3 > \frac{\mu_3 c_3}{2},\] which ensures that the net linear drift in \(y\) and \(z\) is negative after including the \(\frac{\mu_3 c_3}{2}z\) term coming from part (ii).
(iv) Collecting terms. Combining (i)–(iii) in 9 yields, for any fixed small \(\delta,\varepsilon_y>0\) and appropriate large \(c_2,c_3\), \[\begin{align} \mathcal{L}V(x,y,z) &\le -\Big[\,2\mu_1-2\sigma_1^2 - \delta - \varepsilon_y - \frac{c_2^2 (1-\eta)^2 \beta^2}{\mu_3 c_3}\,\Big](x-\bar{x})^2 \\ &\quad - \Big[\,(\mu_2+q)c_2 - c_3(1-\epsilon)p - C_y\,\Big]\, y \\ &\quad - \Big[\,\mu_3 c_3 - \tfrac{\mu_3 c_3}{2}\,\Big]\, z \;+\;C_0, \end{align}\] for constants \(C_y,C_0<\infty\) depending on \((\rho,R,\bar{x},\beta,p,q,\sigma_1)\) but independent of \((x,y,z)\).
Choose \(\delta,\varepsilon_y>0\) small and \(c_2,c_3>0\) large such that the three square brackets are strictly positive. Then there exist \(\theta>0\) and \(C<\infty\) with \[\mathcal{L}V(x,y,z) \le -\theta\big[ (x-\bar{x})^2 + y + z \big] + C.\] Consequently, there is \(R_0>0\) such that \[\mathcal{L}V(x,y,z)\le -\tfrac{\theta}{2}\qquad \text{whenever}\quad (x,y,z)\notin D,\] with \(D=\{V\le R_0\}\) intersected with the box \(\{\rho<x,y,z<R\}\). This verifies ?? .
Step 3: Conclusion. By Lemma 3 the process is positive recurrent and admits a unique invariant distribution \(\pi\); ergodicity follows from the same lemma. ◻
Remark 2. Under the extinction conditions of Theorem 2, \(y(t)\to 0\) and \(z(t)\to 0\) exponentially almost surely, while \(x(t)\) fluctuates about \(\bar{x}\) with the same long-time average as the linear reference process ?? . In this regime, the invariant measure \(\pi\) is supported near the set \(\{(x,0,0):x>0\}\) and integrates polynomials in \(x\) of any fixed order; the Foster drift above yields polynomial moment bounds under \(\pi\).
7 Numerical Methods and Simulations↩︎
The analytical results developed above establish fundamental properties of the HBV stochastic system. To complement these theoretical findings, we implement numerical simulations. The stochastic system 3 lacks a closed-form solution; therefore, a discretization scheme is required. We adopt the Euler–Maruyama method, which is well-suited for systems with multiplicative noise and preserves positivity with sufficiently small time step \(\Delta t\).
7.1 Euler–Maruyama Numerical Method↩︎
To approximate the stochastic system 3 , we apply the Euler–Maruyama method. Let \(t_n=n\Delta t\), \(n=0,1,\dots,N\), with total time \(T=N\Delta t\). For each step, the Wiener increments \(\Delta W_{i,n}=W_i(t_{n+1})-W_i(t_n)\) \((i=1,2,3)\) are independent Gaussian random variables with distribution \(\mathcal{N}(0,\Delta t)\). The discretized update equations are \[\begin{align} x_{n+1} &= x_n + \big(\Lambda - \mu_1 x_n - (1-\eta)\beta x_n z_n + qy_n\big)\Delta t + \sigma_1 x_n \Delta W_{1,n}, \\ y_{n+1} &= y_n + \big((1-\eta)\beta x_n z_n - \mu_2 y_n - qy_n\big)\Delta t + \sigma_2 y_n \Delta W_{2,n}, \\ z_{n+1} &= z_n + \big((1-\epsilon)py_n - \mu_3 z_n\big)\Delta t + \sigma_3 z_n \Delta W_{3,n}. \end{align}\]
The simulation procedure is implemented as follows.
7.1.0.1 Remark (positivity).
We apply a nonnegativity projection after each Euler–Maruyama update. While such a projection can introduce a slight bias, we verified that halving \(\Delta t\) did not alter the qualitative conclusions. Positivity-preserving variants (e.g., truncated/tamed EM or logarithmic transforms) could also be used [23], [24].
7.2 Parameter Values and Setup↩︎
The baseline parameters used in simulations are reported in Table 2. Unless otherwise noted, we fix \(T=200\) days, \(\Delta t=10^{-3}\), and use \(N=2\times 10^5\) iterations. Initial conditions are \((x_0,y_0,z_0)=(10^7,10^5,10^4)\).
| Symbol | Description | Value (units) | Source |
|---|---|---|---|
| \(\Lambda\) | Influx of uninfected hepatocytes | \(1.0\times 10^{7}\;\text{cells day}^{-1}\) | [7] |
| \(\mu_1\) | Death rate of uninfected hepatocytes | \(0.014\;\text{day}^{-1}\) | [8] |
| \(\mu_2\) | Death rate of infected hepatocytes | \(0.24\;\text{day}^{-1}\) | [9] |
| \(\mu_3\) | Viral clearance rate | \(0.24\;\text{day}^{-1}\) | [7] |
| \(\beta\) | Infection rate constant | \(2.4\times 10^{-8}\;\text{(cell}\cdot\text{day)}^{-1}\) | [7] |
| \(p\) | Virions produced per infected cell per day | \(2.4\;\text{virions}\,(\text{cell}\cdot\text{day})^{-1}\) | [7] |
| \(q\) | Non-cytolytic cure rate | \(0.67\;\text{day}^{-1}\) | [8] |
| \(\eta\) | Efficacy blocking infection | baseline: \(0.36\) | [9] |
| \(\epsilon\) | Efficacy blocking production | baseline: \(0.24\) | [9] |
| \(\sigma_1,\sigma_2,\sigma_3\) | Noise intensities | \(\sigma_1=0.01,\;\sigma_2=0.05,\;\sigma_3=0.05\) | [16], [24] |
7.3 Numerical Observations↩︎
The simulations confirm the following properties:
All trajectories remain positive, consistent with Theorem 1.
The infected states \(y(t),z(t)\) decay exponentially under the stability conditions of Theorem 2.
The uninfected hepatocytes converge in distribution to a stationary regime around \(\bar{x}=\Lambda/\mu_1\).
For large noise intensities \((\sigma_2,\sigma_3)\), infection clearance occurs even when the deterministic \(R_0>1\), verifying noise-induced stabilization.
7.4 Disease–free versus endemic regimes↩︎
Figure 2 displays the qualitative regimes induced by treatment efficacy. For \(\eta=0.9\) the infected compartments decay exponentially, whereas for \(\eta=0.1\) the infection persists. The \((x,y)\) and \((y,z)\) portraits separate cleanly across regimes, while the 3D trajectories illustrate attraction toward the disease–free manifold or toward a stochastic neighborhood of an endemic balance.
7.5 Effect of noise intensity on dynamics↩︎
Figure 3 shows a monotone reduction of infected loads as the standard noise amplitude increases. For sufficiently large \(\sigma\), both \(y\) and \(z\) cross below numerical resolution and remain negligible thereafter. This observation is consistent with the negative–definiteness condition used in the proof of Theorem 2, where diffusion compensates deterministic growth.
7.6 Stochastic ensemble versus deterministic path↩︎
In Figure 4, the stochastic mean path of \(x\) tracks the deterministic trend and settles near \(\bar{x}=\Lambda/\mu_1\), while \(y\) and \(z\) decline more rapidly on average when noise is present. This aligns with the explicit solution and expectation identity for the linear surrogate ?? and with the comparison argument in Theorem 3.
7.7 Stability indicators and reproduction numbers↩︎
Here, we visualize a heuristic stochastic reproduction index \(R_0^{\mathrm{stoch}}(\sigma)\) obtained by replacing the linearized drift terms with their Itô-corrected counterparts; in our parameter range, this decreases monotonically with \(\sigma\) and correlates with the extinction behavior observed in simulations.
The top panels of Figure 5 visualize the dynamical counterpart of the analytical thresholds. The bottom panels show the deterministic reproduction number as a function of \(\eta\) and a commonly used stochastic correction that decreases with \(\sigma\), reflecting the role of diffusion in the Lyapunov analysis. The simulated decays in the top panels are in line with these indicators.
7.8 Long–time statistics and stationary laws↩︎
Figure 6 uses a long run to probe the terminal distribution of the process. The mean of \(x\) over the terminal window matches \(\Lambda/\mu_1\) within a small relative error, and both infected components are concentrated near the absorbing edge. These observations support the ergodic stationary regime under the drift–diffusion inequalities used in Theorem 4.
7.9 Sensitivity of \(R_0\) and two–parameter exploration↩︎
These parametric views in Figure 7 complement the time–series evidence by illustrating how control parameters shift the system across the epidemic threshold. The joint surface highlights a region where moderate treatment, combined with moderate noise, is sufficient to reduce the effective reproduction below one.
7.10 Phase–space geometry and vector field inspection↩︎
Figure 8 confirms the geometric picture: the flow field favors decay of \(y\) and \(z\) under strong treatment or strong noise, while for weaker control the trajectories organize around a small enduring set, in accordance with the analytical drift–diffusion tradeoffs.
8 Discussion and Biological Implications↩︎
8.1 Mathematical Contributions↩︎
The present study advances the mathematical analysis of stochastic models for viral infections in several significant ways. First, we established global existence, uniqueness, and positivity of solutions for the HBV system. These properties are indispensable for biological interpretability and provide a rigorous foundation for all subsequent analysis. Second, we developed a detailed stability theory for the infected compartments. Our results demonstrate that stochastic perturbations can alter the qualitative behavior of the system and, under suitable parameter regimes, drive the extinction of infection even when the corresponding deterministic threshold parameter exceeds unity. This phenomenon of noise-induced stabilization has no analogue in deterministic models, illustrating the importance of stochastic methods in epidemiology. Third, by employing Foster–Lyapunov criteria, we obtained explicit conditions for the existence of an ergodic stationary distribution. This provides a mathematically rigorous description of the long-term behavior of the process and clarifies the sense in which the stochastic model equilibrates over time. Taken together, these contributions enrich the general theory of stochastic epidemic systems and extend earlier frameworks such as those developed in [16], [18], [19].
8.2 Biological Significance↩︎
The mathematical findings carry several important biological implications. Environmental fluctuations are intrinsic to biological processes, arising from variability in host responses, heterogeneity in viral replication, and external factors such as co-infections or non-adherence to treatment. The analysis shows that such fluctuations may promote viral clearance in situations where deterministic models predict persistence. This suggests that the immune system or therapeutic interventions may benefit from variability, an observation consistent with emerging studies that emphasize the role of stochasticity in infection outcomes [4]. Moreover, the results on ergodicity indicate that long-term infection dynamics do not converge to fixed values but rather to stationary distributions. This aligns with clinical observations of chronic HBV carriers who exhibit fluctuating but bounded viral loads over extended periods. Finally, the demonstration that uninfected hepatocyte populations converge to a stationary regime around their natural equilibrium highlights the resilience of healthy cells, even under persistent environmental variability.
8.3 Limitations and Future Directions↩︎
Several limitations of the present analysis should be acknowledged. The formulation assumes white noise perturbations, which are temporally uncorrelated. In reality, biological fluctuations often exhibit temporal correlation or burst-like behavior. Incorporating colored noise or Lévy jump processes may capture these features more accurately. Parameter estimation for stochastic epidemic models remains a substantial challenge, particularly in quantifying noise intensities from limited clinical data. Addressing this problem requires integration of statistical inference methods with stochastic modeling [13]. Furthermore, the current model does not explicitly account for adaptive immune responses, spatial heterogeneity within the liver, or interactions with other viruses. Each of these factors can substantially influence HBV dynamics.
Future research should therefore pursue several directions. These include the development of jump-diffusion and hybrid models to describe discrete random events such as immune cell activation, multi-scale approaches that connect intracellular variability with population-level effects, and systematic methods for parameter inference from longitudinal patient data. Another promising line of work lies in optimal control under uncertainty, where treatment protocols are designed with explicit consideration of stochastic fluctuations. Such investigations may provide a theoretical basis for exploiting randomness in the design of antiviral therapies.
9 Conclusion↩︎
The analysis carried out in this work provides a detailed mathematical investigation of hepatitis B virus infection dynamics under the influence of random environmental fluctuations. We established the global well-posedness of the stochastic system, proving that solutions exist for all time, remain unique, and preserve positivity, a property of central importance in any biologically grounded model.
The study of stability revealed that the infected hepatocyte and viral compartments may undergo almost sure exponential extinction, even in parameter regimes where deterministic models predict persistence. This phenomenon, driven by noise-induced stabilization, highlights the potential for stochastic fluctuations to suppress infection in ways that are inaccessible to deterministic dynamics. The behavior of uninfected hepatocytes was analyzed separately, showing convergence to a stochastic equilibrium in agreement with the deterministic average, and the existence of ergodic stationary distributions was demonstrated under appropriate conditions using Foster–Lyapunov techniques.
Numerical simulations based on the Euler–Maruyama scheme confirmed and illustrated the theoretical results. The trajectories consistently remained positive, the extinction of infection was observed in the regimes predicted by the analytical criteria, and the uninfected hepatocyte population approached its stochastic equilibrium. In addition, simulations revealed the existence of stationary distributions for long-term dynamics and provided clear evidence of the capacity of stochasticity to induce viral clearance even when the basic reproduction number is greater than one.
Taken together, these results highlight the crucial role of randomness in determining the outcome of infection. The framework developed here is not limited to HBV but can be applied to other viral infections and epidemic systems where environmental variability exerts a decisive influence. Beyond its theoretical contributions, this work suggests that stochastic effects could be leveraged in the design of intervention strategies, offering new perspectives on treatment and control of viral diseases.
The interplay between rigorous stochastic analysis and biological interpretation underscores the importance of incorporating randomness into models of infectious disease. Future research should further explore the quantitative implications of noise-induced extinction and ergodicity, as well as the development of computational approaches that allow for efficient simulation of increasingly complex stochastic epidemic systems.