Stochastic modeling and computational simulations of HBV infection dynamics

Preprint

Article

Stochastic modeling and computational simulations of HBV infection dynamics

Altmetrics

Downloads

Views

Comments

Abdallah Alsammani^*

This version is not peer-reviewed

Submitted:

02 September 2023

Posted:

05 September 2023

You are already at the latest version

Alerts

Abstract

This study investigates the stochastic dynamics of hepatitis B virus (HBV) infection using a newly proposed stochastic model. Contrary to deterministic models that fail to encapsulate the inherent randomness and fluctuations in biological processes, our stochastic model provides a more realistic representation of HBV infection dynamics. It incorporates random variability, thereby acknowledging the changes in viral and cellular populations and uncertainties in parameters such as infection rates and immune responses. We examine the solution's existence, uniqueness, and positivity for the proposed model, followed by a comprehensive stability analysis. We provide the necessary and sufficient conditions for local and global stability, offering deep insight into the infection dynamics. Furthermore, we utilize numerical simulations to corroborate our theoretical results. This research provides a robust tool for understanding the complex behavior of HBV dynamics, contributing significantly to the ongoing quest for more effective HBV control and prevention strategies.

Keywords:

Subject: Computer Science and Mathematics - Mathematical and Computational Biology

1. Introduction

Hepatitis B is a severe viral infection that affects the liver, leading to acute and chronic conditions. The virus is commonly transmitted through pregnancy, contact with infected blood or body fluids, and unsafe injections. In 2019, 296 million people lived with chronic hepatitis B, with 1.5 million new infections and an estimated 820,000 deaths, primarily from liver-related complications. Effective vaccines are available that offer essential prevention [1]. Mathematical modeling is a powerful tool that extensively enriches our comprehension of the complex HBV infection dynamics and the consequential effects of antiviral therapies. Classical models of HBV infection dynamics are typically anchored in ordinary differential equations (ODEs), reflecting the intricate interactions between the virus and the host immune response [2,3,4,15,17,26].

Although these deterministic models provide pivotal insights into HBV’s pathogenesis and therapeutic interventions’ impacts, they tend to neglect the inherent randomness and variability intrinsic to biological processes. This aspect can be critically important, particularly in shaping the infection dynamics [5]. Therefore, stochastic models have been suggested to offer a more nuanced and accurate representation of HBV dynamics. These models embed random variability into the equations, thereby encapsulating fluctuations in viral and cellular populations and the uncertainty associated with parameters such as infection rates and immune responses [6,7].

One of the first works to provide valuable insight into the nature of the infection and treatment effects are deterministic models, mainly based on ordinary differential equations (ODEs) [2,3,4]. These models represent the virus-host immune response interaction, capturing the core dynamics of the infection.

However, deterministic models’ primary limitation is their inability to consider the inherent variability in biological processes [5,42]. This variability can be integral to understanding the dynamics of infection, including fluctuations in viral and cellular populations and uncertainty in parameters such as infection rates and immune responses.

Recently, stochastic models have been suggested to tackle this limitation and provide a more realistic representation of the HBV infection dynamics [6,7]. These models incorporate random variations into the system, thus accounting for the variability inherent in biological processes.

This research presents a mathematical analysis of the stochastics HBV infection dynamics model. The existence, uniqueness, and probability of the solution are discussed then we studied the local and global stability, where we show the necessary and sufficient conditions of stability in probability. In addition, we show the existence of ergodic stationary distributions. To validate our theoretical findings, we substantiate the results with numerical simulations, offering a robust and practical tool for understanding the complex behavior of HBV dynamics.

2. Introduction

3. Model Formulation

Our research employs a modified stochastic version of a deterministic model that represents the dynamics of HBV infection.

The deterministic model [42] is presented by the following system of ordinary differential equations (ODEs):

\begin{matrix} \frac{d x}{d t} & = Λ - μ_{1} x - (1 - η) β x z + q y, \\ \frac{d y}{d t} & = (1 - η) β x z - μ_{2} y - q y, \\ \frac{d z}{d t} & = (1 - ϵ) p y - μ_{3} z . \end{matrix}

(1)

The parameters in these equations are defined in Table 1.

The stochastic model, built upon the deterministic model, adds stochastic noise to each equation to capture the inherent variability in biological systems:

\begin{matrix} d x (t) & = (Λ - μ_{1} x - (1 - η) β x z + q y) d t + σ_{1} x d W_{1} (t), \\ d y (t) & = ((1 - η) β x z - μ_{2} y - q y) d t + σ_{2} y d W_{2} (t), \\ d z (t) & = ((1 - ϵ) p y - μ_{3} z) d t + σ_{3} z d W_{3} (t) . \end{matrix}

(2)

Here,

W_{1} (t)

W_{2} (t)

, and

W_{3} (t)

represent standard Wiener processes, and

σ_{1}

σ_{2}

, and

σ_{3}

represent the noise coefficients that control the amplitude of the stochastic fluctuations.

The system (2) can be reformulated in the following way:

d u_{t} = f (u (t), t) d t + B (u, t) d W, t \geq 0

(3)

Here, the initial condition is

u_{0} \in R^{3}

, and

u (t)

is represented as a vector

(x (t), y (t), z (t))

The functions

f (u, t)

and

B (u, t)

, along with the differential

d W

, are defined as follows:

\begin{matrix} f (u, t) & = (\begin{matrix} λ - μ_{1} x - (1 - η) β x z + q y \\ (1 - η) β x z - μ_{2} y - q y \\ (1 - ϵ) p y - μ_{3} z \end{matrix}), \\ B (u, t) & = (\begin{matrix} σ_{1} x & 0 & 0 \\ 0 & σ_{2} y & 0 \\ 0 & 0 & σ_{3} z \end{matrix}), \\ d W & = (\begin{matrix} d W_{1} \\ d W_{2} \\ d W_{3} \end{matrix}) . \end{matrix}

3.1. Preliminaries

Before exploring the solution properties and the stability analysis of the stochastic system 2, it is essential to establish a foundation by defining some pertinent terms and theories. These concepts, drawn from [18], lay the groundwork for the forthcoming analysis.

In general, equation 3 can be reformulated as a stochastic equation of d dimension within the context of a complete probability space

(Ω, F, P)

, paired with a filter

{F_{t}}_{t \geq 0}

. This reformulation can be expressed as

d u_{t} = f (u, t) d t + B (u, t) d W, t \geq 0

(4)

Here,

f (u, t) : R^{d} \times [t_{0}, T] \to R^{d}

and

B (u, t) : R^{d} \times [t_{0}, T] \to R^{d \times n}

are Borel measurable, and the white noise

W (t) = (w_{1} (t), w_{2} (t), \dots, w_{n} (t)) \in R^{n}

, for

t \geq 0

. Assuming that

0 \leq t_{0} \leq T < \infty

and that the initial value

u_{0}

F_{t_{0}}

-measurable

R^{d}

random variable such that

E | u_{0} |^{2} < \infty

, equation 4 can be recognized as an Itô type stochastic differential equation.

Definition 1.

A stochastic process

{u (t)}_{0 \leq t_{0} \leq T}

R^{d}

is identified as a solution of equation 4 if it satisfies the following conditions:

i.: ${u (t)}$ is continuous and $F_{t}$ -adapted;
ii.: ${f (u (t), t)} \in L^{1} ([t_{0}, T]; R^{d})$ and ${B (u (t), t)} \in L^{2} ([t_{0}, T]; R^{3} d \times n)$ ;
iii.: equation 4 is valid for every $t \in [t_{0}, T]$ with probability 1.

A solution

{u (t)}

is deemed unique if any other solution

{\bar{u} (t)}

fulfills the following condition:

P {u (t) = \bar{u} (t) for all t_{0} \leq t \leq T} = 1

Lemma 1.

For any

v > 0

, the following inequality holds.

v \leq 2 (v + 1 - ln (v)) - (4 - 2 ln 2) .

Proof.

The proof is straightforward since the function

f (v) = v + 2 - 2 ln (v)

has a minimum at

v = 2

. □

A d-dimensional stochastic equation can, in general, be expressed as

d u (t) = f (u (t), t) d t + B (u (t), t) d W (t),

(5)

Where

u (t) = (x_{1} (t), x_{2} (t), \dots, x_{d} (t))

with the initial condition

u (t_{0}) = u_{0} \in R^{d}

, and W(t) is the m-dimensional white noise defined on a complete probability space

(Ω, F, {F_{t}}_{t \geq 0}, P)

3.2. Itô’s Formula

We define a differential operator L as follows:

L = \frac{\partial}{\partial t} + \sum_{i = 1}^{d} f_{i} (u, t) \frac{\partial}{\partial x_{i}} + \frac{1}{2} \sum_{i, j = 1}^{d} {[W^{T} (u, t) W (u, t)]}_{i j} \frac{\partial^{2}}{\partial u_{i} \partial u_{j}} .

(6)

Consider

V (u, t)

, a nonnegative twice differentiable function defined on

R^{d} \times [t_{0}, \infty)

, applying operator

L

on V, we get:

L V = V_{t} + V_{u} f_{i} (u, t) + \frac{1}{2} trace [W^{T} (u, t) V_{u u} W (u, t)],

(7)

where

V_{t} = \frac{\partial V}{\partial t}

V_{u} = (\frac{\partial V}{\partial x_{1}}, \frac{\partial V}{\partial x_{2}}, \dots \frac{\partial V}{\partial x_{d}})

, and

V_{u u} = {(\frac{\partial^{2} V}{\partial u_{i} \partial u_{j}})}_{d \times d}

Therefore, the Itô formula can be defined as:

d V = L V (u (t), t) d t + V_{u} (u (t), t) B (u (t), t) d W (t) .

(8)

In the next section, we will discuss the properties of the solution to the system given by Eq. (2).

3.3. Properties of Solution

We discuss the solution properties of the system 2, such as existence, uniqueness, and positivity. The dynamics of a population model can be understood by demonstrating that the solution is global and positive for all time

t \geq 0

. The coefficients of the above stochastic system are locally Lipschitz and satisfy the linear growth condition [12].

Theorem 1.

For any initial value

u_{0} \in R_{+}^{3}

, there exists a unique solution to the system (2) for all

t \geq 0

. Furthermore, this solution remains positive for all

t \geq 0

with probability 1, i.e.,

u (t) \in R_{+}^{3}

for all

t \geq 0

almost surely.

Proof.

The coefficients of the equations (2) are continuous and locally Lipschitz. Thus, there is a unique local solution

u (t) \in R_{+}^{3}

for any initial

u_{0} \in R_{+}^{3}

, where

t \in [0, τ)

. To establish that the solution is global, we must show that

τ = \infty

almost surely.

Choose

n_{0} \geq 0

large enough such that

u_{0} \in [1 / n_{0}, n_{0}]

, and let

n > n_{0}

. Define

τ_{n} = inf {t \in [0, τ) : u (t) \notin (1 / n, n)} .

We must show that

τ_{n}

is an empty set, assuming that

inf Φ = \infty

. Since

τ_{n}

is an increasing sequence, let

τ_{\infty} = {lim}_{n \to \infty} τ_{n}

. By definition,

τ_{\infty} < τ

a.s.

To complete the proof, we must show that

τ_{\infty} = \infty

a.s., which implies

τ = \infty

. Assume by contradiction that this is not true. Then there exists a pair of constants

T > 0

and

ϵ \in (0, 1)

such that

P {τ_{\infty} > T} > ϵ

, which implies the existence of an integer

n_{1} > n_{0}

such that

P {τ_{n} < T} \geq ϵ, \forall n \geq n_{1} .

(9)

Now, let’s define a function

G (u)

V (u (t)) = V (x (t), y (t), z (t)) = x + 1 - ln x + y + 1 - ln y + z + 1 - ln z,

which is non-negative due to the inequality

v + 1 - ln v \geq 0

(see [19]). By applying Itô’s formula on

G (u)

, we get the following:

\begin{matrix} d V & = [(1 - \frac{1}{x}) (λ - μ_{1} x - (1 - η) β x z + q y) \\ + (1 - \frac{1}{y}) ((1 - η) β x z - μ_{2} y - q y) \\ + (1 - \frac{1}{z}) ((1 - ϵ) p y - μ_{3} z) + \frac{1}{2} (σ_{1}^{2} + σ_{2}^{2} + σ_{3}^{2})] d t \\ + σ_{1} (x - 1) d W_{1} + σ_{2} (y - 1) d W_{2} + σ_{3} (z - 1) d W_{3}, \end{matrix}

which simplifies to

\begin{matrix} d V & \leq [a + b V (t)] d t + σ_{1} (x - 1) d W_{1} + σ_{2} (y - 1) d W_{2} + σ_{3} (z - 1) d W_{3}, \end{matrix}

where

a = λ + μ_{1} + μ_{2} + μ_{3} + q + \frac{1}{2} (σ_{1}^{2} + σ_{2}^{2} + σ_{3}^{2})

, and

b = max {(1 - ϵ) p, (1 - η) β}

Let

c = max {a, b}

. Integrating the above inequality, if

t_{1} \leq T

we get

\begin{matrix} \int_{0}^{τ_{n} \land t_{1}} d V (u (t)) & \leq \int_{0}^{τ_{n} \land t_{1}} c (1 + V (u (t))) d t + \int_{0}^{τ_{n} \land t_{1}} σ_{1} (x - 1) d W_{1} \\ + \int_{0}^{τ_{n} \land t_{1}} σ_{2} (y - 1) d W_{2} + \int_{0}^{τ_{n} \land t_{1}} σ_{3} (z - 1) d W_{3} \end{matrix}

By definition, this implies

\begin{matrix} E V (u (τ_{n} \land t_{1})) & \leq V (u_{0}) + E \int_{0}^{τ_{n} \land t_{1}} c (1 + V (u (t))) d t, \\ \leq V (u_{0}) + c t_{1} + c E \int_{0}^{τ_{n} \land t_{1}} V (u (t)) d t, \\ \leq V (u_{0}) + c T + c E \int_{0}^{t_{1}} V (τ_{n} \land t_{1}) d t, \\ = V (u_{0}) + c T + c \int_{0}^{t_{1}} E V (τ_{n} \land t_{1}) d t . \end{matrix}

From the Gronwall inequality, we obtain the following.

E V (τ_{n} \land t_{1}) \leq c_{1} = (V (u_{0}) + c T) e^{c T} .

(10)

We then set

Ω_{n} = {τ_{n} \leq T}

for

n \leq n_{1}

and, by inequality 9,

P (Ω_{n}) \geq ϵ

. Note that there exist some

x, y,

or z such that

u (τ_{n}, ω) = n

1 / n

, for every

ω \in Ω

. Therefore,

V (u (τ_{n}, ω))

is greater than

n - 1 - ln (n)

and

\frac{1}{n} + 1 - ln (1 / n) = \frac{1}{n} + 1 + ln (n)

, i.e.,

V (u (τ_{n}, ω)) \geq min (n - 1 - ln (n), \frac{1}{n} + 1 + ln (n)) .

From the inequalities 9 and 10, we obtain

\begin{matrix} c_{1} & \geq E [1_{Ω_{n}} (ω) V (u (τ_{n}, ω))] \\ \geq E [min (n - 1 - ln (n), \frac{1}{n} + 1 + ln (n))], \end{matrix}

where

1_{Ω_{n}}

is the indicator function of

Ω_{n}

. Taking the limit as n approaches ∞ yields

c_{1} = \infty

, which is a contradiction. Therefore,

τ_{\infty} = \infty

a.s., which completes the proof. □

4. Stability Analysis

This section is dedicated to the exploration of the stability analysis of the system (2). For this discussion, we will revisit some requisite definitions and theorems. For further details, see [18] and [13].

Definition 2

([18], pages 110, 119).

(i): The trivial solution of equation 4 is deemed stable in probability if, given any $ϵ \in (0, 1)$ and $r > 0$ , there exists a $δ = δ (ϵ, r, t_{0}) > 0$ such that

$P {| u (t; t_{0}, u_{0}) | < r for all t \geq t_{0}} \geq 1 - ϵ$

whenever $| u_{0} | < δ$ . If not, it is termed stochastically unstable.
(ii): The trivial solution is considered stochastically asymptotically stable if it is stochastically stable and, for each $ϵ \in (0, 1)$ and $r > 0$ , there exists a $δ = δ (ϵ, r, t_{0}) > 0$ such that

$P {lim_{t \to \infty} u (t; t_{0}, u_{0}) = 0} \geq 1 - ϵ$

whenever $| u_{0} | < δ$ .
(iii): The trivial solution is designated as stochastically asymptotically stable in the large if it is stochastically stable and, additionally, for all $u_{0} \in R^{d}$ and $r > 0$ , there is a $δ = δ (ϵ, r, t_{0}) > 0$ such that

$P {lim_{t \to \infty} u (t; t_{0}, u_{0}) = 0} = 1 .$
(iv): The trivial solution of equation 4 is characterized as almost surely exponentially stable if

$lim_{t \to \infty} sup \frac{1}{t} ln | u (t; t_{0}, u_{0}) | < 0$

for all $u_{0} \in R^{d}$ .

Observe that the first equation of our system, 2, does not possess a direct equilibrium point in

R

. However, the second and the third equations of the system 2 have a common equilibrium point at

(y, z) = (0, 0)

The forthcoming theorem will demonstrate that the trivial solution

(0, 0)

is stable. For now, let us concentrate on

y, z

variables. We will return to the equation

x (t)

later to demonstrate that it is stable in distribution. The upcoming theorem asserts that

y (t)

and

z (t)

are exponentially stable under certain conditions.

Theorem 2.

In the system 2,

y (t)

and

z (t)

are almost surely exponentially stable if and only if the following conditions are met

a.: $(1 - ϵ) p - q - μ_{2} + \frac{1}{2} σ_{2}^{2} < 0$
b.: $[(1 - η) β γ - μ_{3}] [(1 - ϵ) p - q - μ_{2}] \leq [1 - ϵ) p - q - μ_{2} + \frac{1}{2} σ_{2}^{2}] [(1 - η) β γ - μ_{3} + \frac{1}{2} σ_{3}^{2}]$

where

γ = max {x}

The proof proceeds as follows.

Proof.

By adding equations

y (t)

and

z (t)

equations, we get

d (y + z) = [((1 - ϵ) p - μ_{2} - q) y + ((1 - η) β x - μ_{3}) z] d t + σ_{2} y d W_{2} + σ_{3} z d W_{3},

Let

V (y, z) = ln (y (t) + z (t))

for

y, z \in R_{+}

. Applying Itô’s formula, we have

\begin{matrix} d V & = [\frac{1}{y + z} ((1 - ϵ) p - q - μ_{2}) y] + [\frac{1}{y + z} ((1 - η) β x - μ_{3}) z] \\ + \frac{1}{2} \frac{σ_{2}^{2} y^{2}}{{(y + z)}^{2}} + \frac{1}{2} \frac{σ_{3}^{2} z^{2}}{{(y + z)}^{2}} d t \\ + \frac{y}{y + z} σ_{2} d W_{2} + \frac{z}{y + z} σ_{3} d W_{3} \\ = \frac{1}{{(y + z)}^{2}} [(y + z) ((1 - ϵ) p - q - μ_{2}) y + ((1 - η) β x - μ_{3}) z] \\ + \frac{1}{2} σ_{2}^{2} y^{2} + \frac{1}{2} σ_{3}^{2} z^{2} \\ + \frac{y}{y + z} σ_{2} d W_{2} + \frac{z}{y + z} σ_{3} d W_{3} \\ \leq \frac{1}{{(y + z)}^{2}} [(y + z) ((1 - ϵ) p - q - μ_{2}) y + ((1 - η) β γ - μ_{3}) z] \\ + \frac{1}{2} σ_{2}^{2} y^{2} + \frac{1}{2} σ_{3}^{2} z^{2} \\ + \frac{y}{y + z} σ_{2} d W_{2} + \frac{z}{y + z} σ_{3} d W_{3}, \\ = \frac{1}{{(y + z)}^{2}} \{(\begin{matrix} y & z \end{matrix}) (\begin{matrix} (1 - ϵ) p - q - μ_{2} + \frac{1}{2} σ_{2}^{2} & (1 - η) β γ - μ_{3} \\ (1 - ϵ) p - q - μ_{2} & (1 - η) β γ - μ_{3} + \frac{1}{2} σ_{3}^{2} \end{matrix}) (\begin{matrix} y \\ z \end{matrix})\} \\ + \frac{y}{y + z} σ_{2} d W_{2} + \frac{z}{y + z} σ_{3} d W_{3} . \end{matrix}

Assuming the conditions of the theorem are met, then the matrix

(\begin{matrix} (1 - ϵ) p - q - μ_{2} + \frac{1}{2} σ_{2}^{2} & (1 - η) β γ - μ_{3} \\ (1 - ϵ) p - q - μ_{2} & (1 - η) β γ - μ_{3} + \frac{1}{2} σ_{3}^{2} \end{matrix})

is negative-definite, which implies it has at least one negative eigenvalue. Let

λ_{max}

denote the largest eigenvalue. With this, the inequality above can be expressed as:

\begin{matrix} d V (y, z) & \leq - | λ_{max} | \frac{1}{{(y + z)}^{2}} (y^{2} + z^{2}) d t + \frac{σ_{2} y}{y + z} d W_{2} + \frac{σ_{3} z}{y + z} d W_{3} . \end{matrix}

By utilizing the inequality

y^{2} + z^{2} \geq 2 y z

, it follows that

(y^{2} + z^{2}) / {(y + z)}^{2} \geq 1 / 2

, yielding:

d V = d ln (y (t) + z (t)) \leq - \frac{1}{2} | λ_{max} | d t + \frac{σ_{2} y}{y + z} d W_{2} + \frac{σ_{3} z}{y + z} d W_{3} .

By integrating the inequality above and applying the fact from [18] that

\underset{t \to \infty}{lim sup} \frac{1}{t} | W_{i} (t) | = 0 for i = 2, 3,

we derive

\underset{t \to \infty}{lim sup} \frac{1}{t} ln (y (t) + z (t)) \leq - \frac{1}{2} | λ_{max} | < 0 almost surely .

This concludes the proof. □

Remark 1.

The stability of components y and z has been achieved without reliance on the reproductive number

R_{0}

, irrespective of whether

R_{0} < 1

R_{0} > 1

. Moreover, it is crucial to note that the conditions within theorem 2 cannot hold in the deterministic case when both

σ_{2} = σ_{3} = 0

We aim to demonstrate the initial component

x (t)

stability. We aim to prove that

x (t)

is stable in distribution, which implies that it is stable around the mean value

λ / μ_{1}

. However, before proceeding, let’s first introduce some necessary lemmas.

Lemma 2.

Assume

W_{1} (t)

is a one-dimensional standard Brownian motion. Then, the expectation is given by

E {e^{σ_{1} (W_{1} (t) - W_{1} (s))}} = e^{\frac{σ_{1}^{2}}{2} (t - s)}, for s \leq t .

Proof.

Define

W = W_{1} (t) - W_{1} (s)

. By the definition of Brownian motion, we know that

W \sim N (0, t - s)

. Hence,

\begin{matrix} E {e^{σ_{1} W}} & = E {e^{σ_{1} (W_{1} (t) - W_{1} (s))}} \\ = \frac{1}{\sqrt{2 π (t - s)}} \int_{- \infty}^{\infty} e^{σ_{1} w} \cdot e^{- \frac{w^{2}}{2 (t - s)}} d w \\ = \frac{1}{\sqrt{2 π (t - s)}} \int_{- \infty}^{\infty} e^{- \frac{{(w + σ_{1} (t - s))}^{2}}{2 (t - s)}} \cdot e^{\frac{σ_{1}^{2}}{2} (t - s)} d w \\ = e^{\frac{σ_{1}^{2}}{2} (t - s)} . \end{matrix}

□

Lemma 3.

Assume that

x_{1} (t)

is a solution to

d x_{1} (t) = (λ - μ_{1} x_{1} (t)) d t + σ_{1} x_{1} (t) d W_{1} (t)

(11)

Then, we have

{lim}_{t \to \infty} E [x_{1} (t)] = λ / μ_{1}

for any initial

x_{1} (0) \in R_{+}

Proof.

Given an initial value

x_{1} (0) \in R_{+}

, there is a unique solution

x_{1} (t)

to equation (11) with the explicit form

x_{1} (t) = x_{1} (0) e^{- (μ_{1} + \frac{1}{2} σ_{1}^{2}) t} + λ \int_{0}^{t} e^{- (μ_{1} + \frac{1}{2} σ_{1}^{2}) (t - s)} \cdot e^{σ_{1} (W_{1} (t) - W_{1} (s))} d s

Applying the expectation to the above equation and using the fact that

W_{1} (0) = 0

, we get

\begin{matrix} E [x_{1} (t)] & = E [x_{1} (0) e^{- (μ_{1} + \frac{1}{2} σ_{1}^{2}) t} + λ \int_{0}^{t} e^{- (μ_{1} + \frac{1}{2} σ_{1}^{2}) (t - s)} \cdot e^{σ_{1} (W_{1} (t) - W_{1} (s))} d s] \\ = x_{1} (0) e^{- μ_{1} t} + \frac{λ}{μ_{1}} (1 - e^{- μ_{1} t}) \end{matrix}

By applying lemma 2, we find

lim_{t \to \infty} E [x_{1} (t)] = \frac{λ}{μ_{1}}

□

Lemma 4.

Assume that

x_{1} (t)

is a solution of equation 11. Then, for any initial value

x_{1} (0) \in R_{+}

, we obtain the following results:

i.: $x_{1} (t)$ admits a unique stationary distribution π, where $Γ (\cdot)$ represents the Gamma function.
ii.: $x_{1} (t)$ satisfies the equation

$lim_{t \to \infty} \frac{1}{t} \int_{0}^{t} x_{1} (s) d s = \int_{0}^{\infty} x_{1} π d x_{1} = \int_{0}^{\infty} f (x_{1}) d x_{1} = \frac{λ}{μ_{1}} .$

Proof.

i.: Define $V (x_{1}) = x_{1} - 1 - ln x_{1}$ , a twice continuously differentiable function. Applying the It $\hat{o}$ formula, we have

$L V (x_{1}) = (1 - \frac{1}{x_{1}}) (λ - μ_{1} x_{1}) + \frac{1}{2} σ_{1}^{2} = - μ - 1 - \frac{λ}{x_{1}} + λ + μ_{1} + \frac{1}{2} σ_{1}^{2} .$

By choosing a sufficiently small $ϵ$ and defining $D = (ϵ, 1 / ϵ)$ , we find that

$L V (x_{1}) \leq - 1, for any x_{1} \in D^{c} .$

This completes the first part of the proof.
ii.: Employing the ergodicity of $x_{1}$ , we obtain

$P \{lim_{t \to - \infty} \frac{1}{t} \int_{0}^{t} x_{1} (s) d s = \int_{0}^{\infty} x_{1} π (d x_{1})\} = 1 .$

Therefore,

$lim_{t \to - \infty} \frac{1}{t} \int_{0}^{t} x_{1} (s) d s = \int_{0}^{\infty} x_{1} π (d x_{1}) = \int_{0}^{\infty} f (x_{1}) x_{1} d x_{1} = \frac{λ}{μ_{1}}, almost surely .$

where we used the fact ${lim}_{t \to \infty} E [x_{1} (t)] = λ / μ_{1} .$

□

Theorem 3.

Assume

(x (t), y (t), z (t))

to be a solution of the system 2, and

x_{1} (t)

to be a solution of equation 11. Under the condition of theorem 2, we can confirm that

lim_{t \to \infty} [x (t) - x_{1} (t)] = 0 in probability .

(12)

Proof.

By leveraging the comparison theorem of stochastic differential equations, we can assert that

x (t) \leq x_{1} (t)

, i.e.,

x (t) - x_{1} (t) \leq 0

(13)

To conclude this proof, it is crucial to demonstrate that

{lim inf}_{t \to \infty} [x (t) - x_{1} (t)] \geq 0

a.s.

We now introduce a stochastic differential equation, which will assist us in this proof:

d x_{r} (t) = [λ - (μ_{1} + r) x_{r} (t)] d t + σ_{1} x_{r} (t) d W_{1} (t),

(14)

where the initial condition is

x_{r} (0) = x (0)

Recall

x (t)

in the system 2,

\begin{matrix} d x (t) & = [λ - μ_{1} x - (1 - η) β x z + q y] d t + σ_{1} x d W_{1} (t), \end{matrix}

and equation 11,

\begin{matrix} d x_{1} (t) & = (λ - μ_{1} x_{1} (t)) d t + σ_{1} x_{1} (t) d W_{1} (t) . \end{matrix}

From fact that

\begin{matrix} \underset{t \to \infty}{lim inf} [x (t) - x_{1} (t)] & = \underset{t \to \infty}{lim inf} (x (t) - x_{r} (t)) + \underset{t \to \infty}{lim inf} (x_{r} (t) - x_{1} (t)] \\ \geq \underset{t \to \infty}{lim inf} [x (t) - x_{r} (t)] + \underset{t \to \infty}{lim inf} [x_{r} (t) - x_{1} (t)], \end{matrix}

We will proceed to prove the following claims.

claim 1:

{lim inf}_{t \to \infty} [x (t) - x_{r} (t)] \geq 0

almost surely.

Proof. Subtracting the given equations, we find

\begin{matrix} d (x (t) - x_{r} (t)) & = [- μ_{1} (x - x_{r}) + r x_{r} - (1 - η) β x z + q y] d t + σ - 1 (x - x_{r}) d W_{1} \\ = [- (μ_{1} + r) (x - x_{r}) + (r - (1 - η) β z) x + q y] d t + σ_{1} (x - x_{r}) d W_{1} . \end{matrix}

This yields a solution.

x (t) - x_{r} (t) = ϕ (t) \int_{0}^{t} ϕ^{- 1} (s) ((r - (1 - η) β z) x + q y) d x (s),

where

ϕ (t) = e^{- (μ_{1} + r + \frac{1}{2} σ_{1}^{2}) t + σ_{1} W_{1} (t)} .

By Theorem 2, we know that

y (t) \to 0

and

z (t) \to 0

almost surely as

t \to \infty

. Therefore,

\begin{matrix} x (t) - x_{r} (t) & = ϕ (t) (\int_{0}^{T} ϕ^{- 1} (s) ((r - (1 - η) β z) x + q y) d x (s) \\ + \int_{T}^{t} ϕ^{- 1} (s) ((r - (1 - η) β z) x + q y) d x (s)), \end{matrix}

for all

ω \in Ω

and

t > T

. Therefore,

x (t) - x_{r} (t) \geq ϕ (t) κ (T),

where

κ (T) = \int_{0}^{T} ϕ^{- 1} (s) ((ϵ - (1 - η) β z) x (s) + q y (s)) d x (s) .

Given

| κ (T) | < \infty

and

ϕ (t) \to 0

almost surely, we conclude

\underset{t \to \infty}{lim inf} [x (t) - x_{r} (t)] \geq 0 a . s .

(15)

□

claim 2

{lim inf}_{t \to \infty} [x_{r} (t) - x_{1} (t)] \geq 0

a.s.

Proof. Beginning with the first equation in the system 2 and 11, we can write

\begin{matrix} d (x_{r} (t) - x_{1} (t)) & = [- μ_{1} (x_{r} (t) - x_{1} (t)) - r x_{r}] d t + σ_{1} (x_{r} (t) - x_{1} (t)) d W (t), \end{matrix}

which leads to the solution

\begin{matrix} x_{r} (t) - x_{1} (t) & = - r \int_{0}^{t} x_{r} e^{- (μ_{1} + σ_{1}^{2} / 2) (t - s) - σ_{1} (W_{1} (t) - W_{1} (s))} d s . \end{matrix}

analog cognize

x_{r}

as the solution for equation 14, which can be expressed explicitly as

\begin{matrix} x_{r} = λ \int_{0}^{t} e^{- (μ_{1} + r + σ_{1}^{2} / 2) (t - s) + σ_{1} (W_{1} (t) - W_{1} (s))} d s . \end{matrix}

Hence,

\begin{matrix} | x_{r} (t) - x_{1} (t) | & = r \int_{0}^{t} x_{r} e^{- (μ_{1} + σ_{1}^{2} / 2) (t - s) - σ_{1} (W_{1} (t) - W_{1} (s))} d s . \end{matrix}

Applying the expectation and invoking lemma 2, we find

\begin{matrix} E | x_{r} (t) - x_{1} (t) | & = r E [\int_{0}^{t} x_{r} e^{- (μ_{1} + σ_{1}^{2} / 2) (t - s) - σ_{1} (W_{1} (t) - W_{1} (s))} d s] \\ = r \int_{0}^{t} E [x_{r} e^{- (μ_{1} + σ_{1}^{2} / 2) (t - s)}] \cdot E [e^{- σ_{1} (W_{1} (t) - W_{1} (s))} d s] \\ = r \int_{0}^{t} e^{- μ_{1} (t - s)} E [x_{r}] d s, (by lemma) . \end{matrix}

Additionally, we know

\begin{matrix} E [x_{r}] & = λ \int_{0}^{t} e^{- (μ_{1} + r) (t - s)} d s \leq \frac{λ}{μ_{1} + r}, \end{matrix}

which implies

\begin{matrix} E | x_{r} (t) - x_{1} (t) | & \leq \frac{λ r e^{- μ_{1} t}}{μ_{1} + r} (e^{μ_{1} t} - 1) \\ lim_{r} \to 0 lim_{t} \to E | x_{r} (t) - x_{1} (t) | & = 0 . \end{matrix}

Therefore,

lim_{r \to 0} lim_{t \to \infty} | x_{r} (t) - x_{1} (t) | = 0 (In probability) .

(16)

□

With 13, 15, and 16, we conclude the proof of the theorem.

4.1. Existence of Ergodic Stationary Distribution

Let

U (t)

be a Markov process in

R^{d}

represented by the following stochastic differential equations:

d U (t) = f (U (t)) d t + \sum_{k = 1}^{n} B_{k} (U (t)) d W_{k} (t) .

(17)

The diffusion matrix is defined as

A (u) = (a_{i j} (u)), where a_{i j} (u) = \sum_{k = 1}^{n} B_{k}^{i} (u) B_{k}^{j} (u) .

(18)

Lemma 5.(See [38,39]) The model in Equation 17 is positive recurrent if there exists a boundary open subset

D \subset R^{d}

with a regular boundary and the following conditions hold:

(A1): There is a positive number M such that

$\sum_{i, j = 1}^{d} a_{i j} (u) ξ_{i} ξ_{j} \geq {| ξ |}^{2}, u \in D and ξ \in R^{d} .$

(19)
(A2): There exists a nonnegative $C^{2}$ -function $V : D^{c} \to R$ such that $L V (u) < - θ$ for some $θ > 0$ , and any $u \in D^{c}$ . Moreover, the positive recurrent process $u (t)$ has a unique stationary distribution $π (\cdot)$ , and

$P \{lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} f (U (t)) d t = \int_{R^{d}} f (u) μ (d u)\} = 1,$

(20)

for all $u \in R^{d}$ , where $f (\cdot)$ is an integrable function with respect to the measure $π (\cdot)$ .

Theorem 4. Assume that

R_{0} > 1

. Under the conditions

σ_{1}^{2} \geq μ_{1} - μ^{*}

σ_{2}^{2} \geq μ_{2} - μ^{*} + \frac{1}{2} (1 - ϵ) p

, and

μ_{3} \geq μ_{3} - \frac{1}{2} (1 - ϵ) p

, the system 2 has a unique ergodic stationary distribution

π (\cdot)

Proof. We have seen that the system 2 has a unique positive solution

x (t), y (t), z (t)

for any initial value

(x_{0}, y_{0}, z_{0}) \in R_{+}^{3}

The diffusion matrix of the system is given by

A = (\begin{matrix} σ_{1}^{2} x^{2} & 0 & 0 \\ 0 & σ_{2}^{2} y^{2} & 0 \\ 0 & 0 & σ_{3}^{2} z^{2} \end{matrix})

Then

\sum_{i, j = 1}^{3} a_{i j} (u) ξ_{i} ξ_{j} = σ_{1}^{2} x^{2} ξ_{1}^{2} + σ_{2}^{2} y^{2} ξ_{2}^{2} + σ_{3}^{2} z^{2} ξ_{3}^{2} \geq M {| ξ |}^{2},

(21)

where,

M = min {σ_{1}^{2} x^{2}, σ_{2}^{2} y^{2}, σ_{3}^{2} z^{2}}

, thus condition

(A 1)

satisfied.

We want to show that condition

(A 2)

is also satisfied by constructing a nonnegative Lyapunov function V such as

L V < 0

Consider the positive functions

V_{1} (x, y) = \frac{1}{2} {(x - \bar{x} + y - \bar{y})}^{2}, and V_{2} (z) = \frac{1}{2} {(z - \bar{z})}^{2}

Now let

V (x, y, z) = V_{1} (x, y) + V_{2} (z) = \frac{1}{2} {(x - \bar{x} + y - \bar{y})}^{2} + \frac{1}{2} {(z - \bar{z})}^{2}

By applying It

\hat{o}

formula we get

L V_{1} (x, y) = (x - \bar{x} + y - \bar{y}) (λ - μ_{1} x - μ_{2} y) + \frac{1}{2} σ_{1}^{2} x^{2} + \frac{1}{2} σ_{2}^{2} y^{2}

since we have

λ - μ_{1} \bar{x} - μ_{2} \bar{y} = 0, \Rightarrow λ = μ_{1} \bar{x} + μ_{2} \bar{y}

, then we get

\begin{matrix} L V_{1} (x, y) & = (x - \bar{x} + y - \bar{y}) [- μ_{1} (x - \bar{x}) - μ_{2} (y - \bar{y})] + \frac{1}{2} σ_{1}^{2} x^{2} + \frac{1}{2} σ_{2}^{2} y^{2}, \\ = - μ_{1} {(x - \bar{x})}^{2} - μ_{2} {(y - \bar{y})}^{2} - μ_{1} (x - \bar{x}) (y - \bar{y}) - μ_{2} (x - \bar{x}) (y - \bar{y}) \\ + \frac{1}{2} σ_{1}^{2} x^{2} + \frac{1}{2} σ_{2}^{2} y^{2} \end{matrix}

By using the fact that,

x^{2} \leq 2 {(x - a)}^{2} + 2 a^{2}

, we get

\frac{1}{2} σ_{1}^{2} x^{2} \leq σ_{1}^{2} {(x - \bar{x})}^{2} + σ_{1}^{2} {\bar{x}}^{2} and \frac{1}{2} σ_{2}^{2} y^{2} \leq σ_{2}^{2} {(y - \bar{y})}^{2} + σ_{2}^{2} {\bar{y}}^{2}

thus

\begin{matrix} L V_{1} (x, y) & \leq - μ_{1} {(x - \bar{x})}^{2} - μ_{2} {(y - \bar{y})}^{2} - 2 μ^{*} (x - \bar{x}) (y - \bar{y}) + σ_{1}^{2} {(x - \bar{x})}^{2} + σ_{1}^{2} {\bar{x}}^{2} \\ + σ_{2}^{2} {(y - \bar{y})}^{2} + σ_{2}^{2} {\bar{y}}^{2} \end{matrix}

where

μ^{*} = min {μ_{1}, μ_{2}}

, and since

- 2 μ^{*} (x - \bar{x}) (y - \bar{y}) \leq μ * {(x - \bar{x})}^{2} + μ^{*} {(y - \bar{y})}^{2}

, then

\begin{matrix} L V_{1} (x, y) & \leq & - (μ_{1} - μ^{*} - σ_{1}^{2}) {(x - \bar{x})}^{2} - (μ_{2} - μ^{*} - σ_{2}^{2}) {(y - \bar{y})}^{2} \\ + σ_{1}^{2} \bar{x} + σ_{2}^{2} \bar{y} \end{matrix}

Similarly,

L V_{2} (z) \leq - (μ_{3} - \frac{1}{2} (1 - ϵ) p - σ_{3}^{2}) {(z - \bar{z})}^{2} + \frac{1}{2} (1 - ϵ) p {(y - \bar{y})}^{2} + σ_{3}^{2} {\bar{z}}^{2}

Thus,

\begin{matrix} L V (x, y, z) & = L V_{1} (x, y) + L V_{2} (z) \\ \leq - (μ_{1} - μ^{*} - σ_{1}^{2}) {(x - \bar{x})}^{2} - (μ_{2} - μ^{*} - σ_{2}^{2}) {(y - \bar{y})}^{2} + σ_{1}^{2} \bar{x} + σ_{2}^{2} \bar{y} \\ - (μ_{3} - \frac{1}{2} (1 - ϵ) p - σ_{3}^{2}) {(z - \bar{z})}^{2} + \frac{1}{2} (1 - ϵ) p {(y - \bar{y})}^{2} + σ_{3}^{2} {\bar{z}}^{2} \\ = - (μ_{1} - μ^{*} - σ_{1}^{2}) {(x - \bar{x})}^{2} - (μ_{2} - μ^{*} + \frac{1}{2} (1 - ϵ) p - σ_{2}^{2}) {(y - \bar{y})}^{2} \\ - (μ_{3} - \frac{1}{2} (1 - ϵ) p - σ_{3}^{2}) {(z - \bar{z})}^{2} + σ_{1}^{2} \bar{x} + σ_{2}^{2} \bar{y} + σ_{3}^{2} {\bar{z}}^{2} \\ = - k_{1} {(x - \bar{x})}^{2} - k_{2} {(y - \bar{y})}^{2} - k_{3} {(z - \bar{z})}^{2} + ω \end{matrix}

where

k_{1} = μ_{1} - μ^{*} - σ_{1}^{2}

k_{2} = μ_{2} - μ^{*} + \frac{1}{2} (1 - ϵ) p - σ_{2}^{2}

k_{3} = μ_{3} - \frac{1}{2} (1 - ϵ) p - σ_{3}^{2}

, and

ω = σ_{1}^{2} \bar{x} + σ_{2}^{2} \bar{y} + σ_{3}^{2} {\bar{z}}^{2}

Since

k_{1}, k_{2}, k_{3}

are positive and by the same computation as in [40], we obtained

\underset{t \to \infty}{lim sup} \frac{1}{t} \int_{0}^{t} [k_{1} {(x (s) - \bar{x})}^{2} + k_{2} {(y (s) - \bar{y})}^{2} + k_{3} {(z (s) - \bar{z})}^{2}] d s \leq ω .

then the ellipsoid

k_{1} {(x - \bar{x})}^{2} + k_{2} {(y - \bar{y})}^{2} + k_{3} {(z - \bar{z})}^{2} + ω = 0

lies entirely in

R_{+}^{3}

, so we can take any neighborhood D of this ellipsoid, such that

L V (u) < - θ, for all u \in D^{c} .

Therefore, condition (A2) also holds and completes the proof. □

5. Numerical Solution

Numerous numerical methods exist for solving stochastic dynamical systems. Some prominent techniques include the Euler-Maruyama Method, Milstein Method, Stochastic Runge-Kutta Methods, Strong and Weak Taylor Methods, Split-Step Methods, Stochastic Theta Method, Multilevel Monte Carlo Methods, Gaussian Process Emulators, Stochastic Collocation and Galerkin Methods, Particle Filters, and Hybrid Methods. This study considers only the Euler-Maruyama and Milstein methods.

Euler-Maruyama method: As one of the simplest numerical methods for SDEs, the Euler-Maruyama method extends the deterministic Euler method to the stochastic context [9]. Despite its simplicity and ease of implementation, the method only provides strong convergence of order $0.5$ .
Milstein method: The Milstein method improves upon the Euler-Maruyama method by including additional terms in the Taylor series expansion of the SDE. This method provides a strong convergence of order $1.0$ under suitable conditions, doubling the convergence rate of the Euler-Maruyama method [10].

These methods and many others provide a rich toolbox for the numerical solution of SDEs. The choice of method depends on several factors, such as the specific form of the SDE, the regularity of its coefficients, and the desired balance between accuracy and computational cost. For our model 3, we will solve it numerically using the Euler-Maruyama and Milstein methods as follows.

Euler-Maruyama Method

To solve the stochastic model numerically, we use the Euler-Maruyama scheme, a stochastic analog of the well-known Euler method for ordinary differential equations. The system gives the stochastic model to be solved (2). The Euler-Maruyama scheme for system (2) can be described as follows:

Discretize the time domain $[0, T]$ into n intervals of size $Δ t = T / n$ , where $t_{i} = i Δ t$ . Denote $x (t_{i})$ , $y (t_{i})$ , and $z (t_{i})$ by $x_{i}$ , $y_{i}$ , and $z_{i}$ respectively.
Initialize $x_{0}$ , $y_{0}$ and $z_{0}$ as the initial conditions.
For each time step $i = 0, 1, . . ., n - 1$ , compute the increments $Δ W_{1} = W_{1} (t_{i + 1}) - W_{1} (t_{i})$ , $Δ W_{2} = W_{2} (t_{i + 1}) - W_{2} (t_{i})$ , and $Δ W_{3} = W_{3} (t_{i + 1}) - W_{3} (t_{i})$
Update the solution from time $t_{i}$ to $t_{i + 1}$ as follows

$\begin{matrix} x_{i + 1} & = x_{i} + (Λ - μ_{1} x_{i} - (1 - η) β x_{i} z_{i} + q y_{i}) Δ t + σ_{1} x_{i} Δ W_{1}, \\ y_{i + 1} & = y_{i} + ((1 - η) β x_{i} z_{i} - μ_{2} y_{i} - q y_{i}) Δ t + σ_{2} y_{i} Δ W_{2}, \\ z_{i + 1} & = z_{i} + ((1 - ϵ) p y_{i} - μ_{3} z_{i}) Δ t + σ_{3} z_{i} Δ W_{3} . \end{matrix}$

This procedure gives a numerical approximation of the solution to the stochastic differential equations in the system (2). To implement this scheme, one needs to generate increments

Δ W_{1}

Δ W_{2}

, and

Δ W_{3}

, normally distributed random numbers with mean 0 and variance

Δ t

5.1. Milstein Method

The Milstein method is another widely used numerical approach to solving stochastic differential equations, which extends the Euler-Maruyama method by including an additional term to account for the diffusion term in the stochastic system.

The scheme for system (2) using the Milstein method can be detailed as follows:

As in the Euler-Maruyama method, discretize the time domain $[0, T]$ into n intervals of size $Δ t = T / n$ and initialize $x_{0}$ , $y_{0}$ , and $z_{0}$ as the initial conditions.
For each time step $i = 0, 1, . . ., n - 1$ , calculate the increments $Δ W_{1}$ , $Δ W_{2}$ , and $Δ W_{3}$ in the same way as the Euler-Maruyama method.
Update the solution from time $t_{i}$ to $t_{i + 1}$ by

$\begin{matrix} x_{i + 1} & = x_{i} + (Λ - μ_{1} x_{i} - (1 - η) β x_{i} z_{i} + q y_{i}) Δ t + σ_{1} x_{i} Δ W_{1} + \frac{1}{2} σ_{1}^{2} x_{i} (Δ W_{1}^{2} - Δ t), \\ y_{i + 1} & = y_{i} + ((1 - η) β x_{i} z_{i} - μ_{2} y_{i} - q y_{i}) Δ t + σ_{2} y_{i} Δ W_{2} + \frac{1}{2} σ_{2}^{2} y_{i} (Δ W_{2}^{2} - Δ t), \\ z_{i + 1} & = z_{i} + ((1 - ϵ) p y_{i} - μ_{3} z_{i}) Δ t + σ_{3} z_{i} Δ W_{3} + \frac{1}{2} σ_{3}^{2} z_{i} (Δ W_{3}^{2} - Δ t) . \end{matrix}$

This approach considers the second moment of the diffusion term, thereby providing more accurate approximations of the solution of the stochastic system (2). The implementation involves a similar procedure to the Euler-Maruyama method but with additional terms to be computed at each step.

Omparison between Euler-Maruyama and Milstein Methods

Figure 1. Figure describes the solution of the model 2 using Euler-Maruyama Method. Parameter values are defined in the Table 3

Figure 2. Figure describes the solution of the model 2 using Milstein Method. Parameter values are defined in the Table 3

Figure 3. Figure describes the deterministic solution of the model 2. Parameter values are defined in the Table 3

Figure 4. Figure compares the stochastic "Euler-Maruyama" and deterministic solutions. Parameter values are defined in the Table 3

Figure 5. Figure compares the stochastic " Milstein" and deterministic solutions. Parameter values are defined in the Table 3

Figure 6. Figure compares the stochastic " Milstein" and deterministic solutions. Parameter values are defined in the Table 3

Figure 7. Figure compares the stochastic " Milstein" and deterministic solutions. Parameter values are defined in the Table 3

Figure 8. Figure compares the stochastic " Milstein" and deterministic solutions. Parameter values are defined in the Table 3

Table 2 summarizes the main differences between the two methods, highlighting the trade-offs between simplicity, computational cost, accuracy, and stability. Both methods provide valuable tools for numerically solving stochastic models and can be chosen according to the specific needs and constraints of the problem.

6. Conclusion

This paper presents a refined stochastic model for Hepatitis B Virus (HBV) infection dynamics, capturing biological processes’ inherent variability and randomness. By introducing preliminary concepts, the study fosters easy comprehension, paving the way for an analysis that ensures the solution’s existence, uniqueness, and positive for all positive initial values.

The stability analysis of the model is explored, delineating necessary and sufficient conditions for stability in probability. This contributes to a more profound understanding of the dynamics of HBV infection and is a robust foundation for further investigation.

The study employs two distinguished numerical methods, Euler-Maruyama and Milstein, to solve the stochastic differential equations in the system (2). While the Euler-Maruyama method provides a simple solution with a lower convergence rate, the Milstein method enhances accuracy and stability by accounting for the second moment of the diffusion term. These numerical simulations corroborate the theoretical findings and illustrate the practical applicability of the model.

Despite certain limitations, this research marks a significant advance in the field of HBV infection dynamics. It offers a deeper understanding of the subject and encourages progress toward devising effective strategies for control and prevention. The insights gained from this study serve as a vital stepping stone for future explorations in this critical area of public health.

Conflicts of Interest:

Declare conflicts of interest or state “The authors declare no conflict of interest.” Authors must identify and declare any personal circumstances or interests that may be perceived as inappropriately influencing the representation or interpretation of reported research results. Any role of the funders in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results must be declared in this section. If there is no role, please state “The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results”.

References

World Health Organization, Hepatitis B, 2023, https://www.who.int/news-room/fact-sheets/detail/hepatitis-b, Accessed: 2023-07-23.
Nowak, Martin A. and others, Viral dynamics of primary viremia and antiretroviral therapy in simian immunodeficiency virus infection, Journal of Virology, vol. 70, no. 10, pp. 6735–6741, 1996, Am Soc Microbiol. [CrossRef]
Bangham, Charles R., Population dynamics of immune responses to persistent viruses, Science, vol. 272, no. 5258, pp. 74–79, 1996, American Association for the Advancement of Science. [CrossRef]
Perelson, Alan S. and Nelson, Patrick W., Modelling viral and immune system dynamics, Nature Reviews Immunology, vol. 2, no. 1, pp. 28–36, 2002, Nature Publishing Group.
Lloyd, Alun L., Realistic distributions of infectious periods in epidemic models: Changing patterns of persistence and dynamics, Theoretical Population Biology, vol. 60, no. 1, pp. 59–71, 2001, Elsevier. [CrossRef]
Ribeiro, Andre S., Effects of stochastic population dynamics on the quantification of gene expression, Physical Review E, vol. 71, no. 1, pp. 011912, 2005, APS.
Guedj, Jeremie and Perelson, Alan S., Modeling shows that the NS5A inhibitor daclatasvir has two modes of action and yields a shorter estimate of the hepatitis C virus half-life, Proceedings of the National Academy of Sciences, vol. 110, no. 10, pp. 3991–3996, 2013, National Academy of Sciences. [CrossRef]
Gardiner, Crispin, Handbook of stochastic methods: for physics, chemistry and the natural sciences, 2004, Springer Science & Business Media.
Kloeden, Peter E. and Platen, Eckhard, Numerical solution of stochastic differential equations, 2013, Springer Science & Business Media.
Higham, Desmond J., Algorithmic introduction to numerical simulation of stochastic differential equations, SIAM review, vol. 43, no. 3, pp. 525–546, 2001, SIAM. [CrossRef]
Rackauckas, Christopher and Nie, Qing, Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling with memory, Discrete and continuous dynamical systems. Series B, vol. 22, no. 7, pp. 2731–2761, 2017, American Institute of Mathematical Sciences. [CrossRef]
Has’minskii, R. Z., Stochastic Stability of Differential Equations, 1980, Sijthoff Noordhoff, Alphen aan den Rijn, The Netherlands. [CrossRef]
Caraballo, Tomas and Han, Xiaoying, Applied Nonautonomous and Random Dynamical Systems, 2016, SpringerBriefs. [CrossRef]
arXiv:2305.08210, 2023.Alsammani, Abdallah, Mathematical Analysis of Autonomous and Nonautonomous Hepatitis B Virus Transmission Models, arXiv preprint arXiv:2305.08210, 2023.
Alsammani, Abdallah Alhadi Mahadi, Dynamical Behavior of Nonautonomous and Stochastic HBV Infection Model, Ph.D. thesis, Auburn University, 2020.
Nowak, Martin A. and others, Viral dynamics in hepatitis B virus infection, Proceedings of the National Academy of Sciences, vol. 93, no. 9, pp. 4398–4402, 1996.
Payne, R.J.H. and Nowak, M.A. and Blumberg, B., The dynamics of hepatitis B virus infection, Proc. Natl. Acad. Sci. USA, vol. 93, pp. 6542–6546, 1996.
Mao, X., Stochastic Differential Equations and Applications, 2nd ed., Horwood Publishing Limited, 2007.
Dalal, N. and Greenhalgh, D. and Mao, X., A Stochastic Model for Internal HIV Dynamics, Journal of Mathematical Analysis and Applications, 2008.
Hui, H. and Nie, L., Analysis of a Stochastic HBV Infection Model with Nonlinear Incidence Rate, Journal Of Biological Systems, vol. 27, no. 3, 2019. [CrossRef]
Smith, Hall L., Monotone Dynamical Systems, AMS, 1995. [CrossRef]
Leenheer, Patric DE and Smith, Hall L., Virus Dynamics: a Global Analysis, SIAM J. Appl Math, vol. 63, no. 4, 2003. [CrossRef]
Nowak, Martin A. and Bangham, Charles R.M., Population Dynamics of Immune Responses to Persistent Viruses, Science, vol. 271, pp. 74–79, 1996. [CrossRef]
Nowak, Martin A. and Bangham, Charles R.M., Virus Dynamics, Oxford University Press, 2000.
Perelson, A.S. and Nelson, P.W., Mathematical analysis of HIV-1 dynamics in vivo, SIAM Rev., vol. 41, pp. 3–44, 1999. [CrossRef]
Abdulrashid, I. and Alsammani, A. and Han, X., Stability Analysis of Chemotherapy Model with Delays, Discrete and Continuous Dynamical Systems Series B, vol. 24, 2019. [CrossRef]
Perelson, A.S. and Neumann, A.U. and Markowitz, M. and Leonard, J.M. and Ho, D.D., HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time, Science, vol. 271, pp. 1582–1586, 196. [CrossRef]
Gourley, S.A. and Kuang, Y. and Nagy, J.D., Dynamics of a delay differential model of hepatitis B virus infection, J. Biological Dynamics, vol. 2, pp. 140–153, 2008. [CrossRef]
Lewin, S.R. and Riberiro, R.M. and Walters, T., and others, Analysis of Hepatitis B viral load decline under potent therapy: Complex decay profiles observed, Hepatology, vol. 34, pp. 1012–1020, 2001. [CrossRef]
Dahari, H. and Lo, A. and Ribeiro, M. and Pereslon, A.S., Modeling Hepatitis C Virus Dynamics: Liver regeneration and critical drug efficacy, J. Theor. Biol., vol. 47, pp. 371–381, 2007. [CrossRef]
Reluga, T.C. and Dahari, H. and Perelson, A.S., Analysis of Hepatitis C Virus infection Models with Hepatocyte Homeostasis, SIAM J. APPL. MATH, vol. 69, no. 4, pp. 999–1023, 2008.
Alsammani, Abdallah Alhadi Mahadi, Dynamical Behavior of Nonautonomous and Stochastic HBV Infection Model, Ph.D. thesis, 2020.
Perko, Lawrence, Differential Equations and Dynamical Systems, Springer-Verlag New York, 1991.
Stafford, M.A. and Corey, L. and Cao, Y. and Daar, E.S. and Ho, D.D. and Perelson, A.S., Modeling plasma virus concentration during primary HIV infection, J. Theor. Biol., vol. 203, pp. 285–301, 2000. [CrossRef]
Hale, Jack K., Ordinary Differential Equations, New York, 2009.
Kloeden, Peter E. and Rasmussen, Martin, Nonautonomous Dynamical Systems, American Mathematical Society, Providence, RI, 2011.
P. van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., vol. 180, pp. 29-48, 2002. [CrossRef]
Li, Dan and Cui, Jing’an and Liu, Meng and Liu, Shengqiang, The evolutionary dynamics of a stochastic epidemic model with nonlinear incidence rate, Bulletin of mathematical biology, vol. 77, pp. 1705–1743, 2015. [CrossRef]
Zhu, Chao and Yin, George, Asymptotic properties of hybrid diffusion systems, SIAM Journal on Control and Optimization, vol. 46, no. 4, pp. 1155–1179, 2007. [CrossRef]
Yang, Qingshan and Jiang, Daqing and Shi, Ningzhong and Ji, Chunyan, The ergodicity and extinction of stochastically perturbed SIR and SEIR epidemic models with saturated incidence, Journal of Mathematical Analysis and Applications, vol. 388, no. 1, pp. 248–271, 2012. [CrossRef]
Alsammani, Abdallah, Mathematical Analysis of Autonomous and Nonautonomous Hepatitis B Virus Transmission Models, In Computational Science and Its Applications – ICCSA 2023 Workshops, Springer Nature Switzerland, pp. 327–343, 2023. [CrossRef]
Alsammani, Abdallah, Mathematical Analysis of Autonomous and Nonautonomous Hepatitis B Virus Transmission Models, In Computational Science and Its Applications – ICCSA 2023 Workshops, Springer Nature Switzerland, pp. 327–343, 2023. [CrossRef]

Table 1. Parameter descriptions for the deterministic model.

Parameter	Description
$Λ$	Production rate of uninfected cells x.
$μ_{1}$	Death rate of x-cells.
$μ_{2}$	Death rate of y-cells.
$μ_{3}$	Free virus cleared rate.
$η$	Fraction that reduced infected rate after treatment with an antiviral drug.
$ϵ$	Fraction that reduced free virus rate after treatment with an antiviral drug.
p	Free virus production rate from y-cells.
$β$	Infection rate of x-cells by free virus z.
q	Spontaneous cure rate of y-cells by non-cytolytic process.

Table 2. Comparison between Euler-Maruyama and Milstein Methods

Aspect	Euler-Maruyama	Milstein
Simplicity	High	Medium
Accuracy	First-order	Second-order (diffusion)
Computational Cost	Low	Medium
Stability	May vary	Generally better
Implementation	Easier	More complex

Table 3. Parameter values for the stochastic and deterministic models

Parameter	$Λ$	$μ_{1}$	$μ_{2}$	$μ_{3}$	$β$	$η$	$ϵ$	p	q	$σ_{1}$	$σ_{2}$	$σ_{3}$
Value	100	20	5	7	0.6	0.6	0.2	2	5	0.5	0.6	0.8

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Stochastic modeling and computational simulations of HBV infection dynamics

Abstract

1. Introduction

2. Introduction

3. Model Formulation

3.1. Preliminaries

3.2. Itô’s Formula

3.3. Properties of Solution

4. Stability Analysis

4.1. Existence of Ergodic Stationary Distribution

5. Numerical Solution

Euler-Maruyama Method

5.1. Milstein Method

Omparison between Euler-Maruyama and Milstein Methods

6. Conclusion

Conflicts of Interest:

References

MDPI Initiatives

Important Links

Subscribe