$\displaystyle\leq~{}d(\mu_{a},\mu_{1})+O(N^{1-3\alpha/8})=O(N^{1-3\alpha/8}),$
(98)

where the last step follows from the fact that, the partial derivatives of
$I_{a}(\cdot)$ with respect to $N_{a}$ is $d(\mu_{a},x_{1,a}(N))\leq
d(\mu_{a},\mu_{1})$. As a result, since arm $a$ was pulled in iteration $N$,
the increment $I_{a}(N)-I_{a}(N-1)$ in $I_{a}(\cdot)$ is upper bounded by
$d(\mu_{a},\mu_{1})$.

Therefore (F.6) implies,

$\displaystyle
S_{a,b}(\mathbold{\widetilde{N}}(N+R),\mathbold{\widetilde{\mu}}(N+R))~{}$
$\displaystyle\leq~{}\Delta\widetilde{N}_{1}\cdot(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-\Delta\widetilde{N}_{b}\cdot
d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})$
$\displaystyle~{}+S_{a,b}(\mathbold{\widetilde{N}}(N),\mathbold{\widetilde{\mu}}(N))+O\left(\left(N^{-3\alpha/8}+RN^{-1}\right)R\right)$
$\displaystyle\leq~{}\Delta\widetilde{N}_{1}\cdot(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-\Delta\widetilde{N}_{b}\cdot
d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})$
$\displaystyle+O\left(N^{1-3\alpha/8}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right).$
(99)

Now we consider two possibilities:

Case I
$d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b})<0$
: In this case, by (F.6),
$S_{a,b}(\mathbold{\widetilde{N}}(N+R),\mathbold{\widetilde{\mu}}(N+R))$ is
upper bounded by,

$-\Delta\widetilde{N}_{b}\cdot
d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})+O\left(N^{1-3\alpha/8}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right).$

By Lemma F.13 and (33), since both $\widetilde{N}_{1}(N)$ and
$\widetilde{N}_{b}(N)$ are $\Theta(N)$, we have
$d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})=\Theta(1)$. As a result, (F.16)
follows.

Case II
$d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b})\geq
0$ : In this case, the RHS of (F.6) can be rewritten as,

$\Delta\widetilde{N}_{b}\left(\frac{\Delta\widetilde{N}_{1}}{\Delta\widetilde{N}_{b}}(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})\right)+O\left(N^{1-3\alpha/8}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right).$
(100)

Using (88) in Lemma F.14, the upper bound becomes

$\displaystyle\Delta\widetilde{N}_{b}\left(\frac{\widetilde{N}_{1}}{\widetilde{N}_{b}}(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})\right)+O\left(N^{1-3\alpha/8}(1+RN^{-1})^{3}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right)$
$\displaystyle=\Delta\widetilde{N}_{b}\left(\frac{\widetilde{N}_{1}}{\widetilde{N}_{b}}(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})\right)+O\left(N^{1-3\alpha/8}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right)\quad\text{(since
$R\leq N$).}$ (101)

By (F.6), we have,

$\displaystyle\mathcal{I}_{a}(N)~{}$
$\displaystyle\leq~{}\mathcal{I}_{b}(N)+O(N^{1-3\alpha/8}).$ (102)

Expanding the empirical index terms, we get,

$\widetilde{N}_{1}d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})+\widetilde{N}_{a}d(\widetilde{\mu}_{a},\widetilde{x}_{1,a})~{}\leq~{}\widetilde{N}_{1}d(\widetilde{\mu}_{1},\widetilde{x}_{1,b})+\widetilde{N}_{b}d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})+O(N^{1-3\alpha/8}).$

By Lemma F.13 we have $\widetilde{N}_{b}=\Theta(N)$. As a result, upon
dividing both sides of the above inequality by $\widetilde{N}_{b}$ and after
some re-arrangement of terms, we obtain,

$\frac{\widetilde{N}_{1}}{\widetilde{N}_{b}}(d(\widetilde{\mu}_{1},\widetilde{x}_{1,j})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})~{}\leq~{}-\frac{\widetilde{N_{j}}}{\widetilde{N_{b}}}d(\widetilde{\mu}_{j},\widetilde{x}_{1,j})+O(N^{-3\alpha/8}).$

Using the above inequality in (F.6), we get the upper bound,

$\displaystyle-\Delta\widetilde{N}_{b}\left(\frac{\widetilde{N}_{j}}{\widetilde{N}_{b}}d(\widetilde{\mu}_{a},\widetilde{x}_{1,a})\right)+$
$\displaystyle
O\left(N^{1-3\alpha/8}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right)$ (103)

By Lemma F.13 and (33), we have,
$\frac{\widetilde{N}_{j}}{\widetilde{N}_{b}}d(\widetilde{\mu}_{j},\widetilde{x}_{1,j})=\Theta(1)$.
Therefore (F.16) follows.

IAT2 Algorithm: The proof of (F.16) for IAT2 is very similar to the proof of
(F.16) for AT2. We first consider the difference of the modified empirical
indexes between the two arms $a$ and $b$,

$S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N),\mathbold{\widetilde{\mu}}(N))~{}=~{}\mathcal{I}_{a}^{(m)}(N)-\mathcal{I}_{b}^{(m)}(N).$

Following a similar procedure as the AT2 algorithm, we apply the mean value
theorem to obtain,

$\displaystyle
S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N+R),\mathbold{\widetilde{\mu}}(N+R))-S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N),\mathbold{\widetilde{\mu}}(N))~{}$
$\displaystyle=~{}\sum_{j=1,a,b}\frac{\partial
S_{a,b}^{(m)}}{\partial\mu_{j}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})\cdot\Delta\widetilde{\mu}_{j}$
$\displaystyle+\frac{\partial S_{a,b}^{(m)}}{\partial
N_{1}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})\cdot\Delta\widetilde{N}_{1}+\frac{\partial
S_{a,b}^{(m)}}{\partial
N_{b}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})\cdot\Delta\widetilde{N}_{b},$
(104)

where
$\Delta\widetilde{\mu}_{j}=\widetilde{\mu}_{j}(N+R)-\widetilde{\mu}_{j}(N)$,
and
$(\mathbold{\hat{N}},\mathbold{\hat{\mu}})=(\hat{N}_{1},\hat{N}_{a},\hat{N}_{b},\hat{\mu}_{1},\hat{\mu}_{a},\hat{\mu}_{b})$,
such that, $\hat{\mu}_{j}$ lies between $\widetilde{\mu}_{j}(N)$ and
$\widetilde{\mu}_{j}(N+R)$, and
$\hat{N}_{j}\in\left[\widetilde{N}_{j}(N),\widetilde{N}_{j}(N+R)\right]$, for
every $j=1,a,b$.

The contribution to (F.6) due to noise in $\mathbold{\widetilde{\mu}}$ is
$O(N^{1-3\alpha/8})$ by the same argument we proved (F.6) for AT2.

Now we consider the partial derivatives of $S_{a,b}^{(m)}$ with respect to
$N_{1},N_{a},N_{b}$. Following the same steps as used for AT2 algorithm, we
have,

$\displaystyle\frac{\partial S_{a,b}^{(m)}}{\partial
N_{1}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})~{}$
$\displaystyle=~{}\frac{\partial S_{a,b}}{\partial
N_{1}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})~{}=~{}d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b})+O(N^{-3\alpha/8}+RN^{-1}),~{}\text{and}$
$\displaystyle\frac{\partial S_{a,b}^{(m)}}{\partial
N_{b}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})~{}$
$\displaystyle=~{}\frac{\partial S_{a,b}}{\partial
N_{b}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})-\frac{1}{\hat{N}_{b}}~{}\leq~{}-d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})+O(N^{-3\alpha/8}+RN^{-1}).$

As a result, the contribution to (F.6) due to $\widetilde{N}_{1}$ and
$\widetilde{N}_{b}$ is,

$\displaystyle\frac{\partial S_{a,b}^{(m)}}{\partial
N_{1}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})\cdot\Delta\widetilde{N}_{1}+\frac{\partial
S_{a,b}^{(m)}}{\partial
N_{b}}(\mathbold{\hat{N}},\mathbold{\hat{\mu}})\cdot\Delta\widetilde{N}_{b}$
$\displaystyle\leq~{}\Delta\widetilde{N}_{1}\cdot(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-\Delta\widetilde{N}_{b}\cdot
d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})+O\left(\left(N^{-3\alpha/8}+RN^{-1}\right)R\right).$

Therefore, adding the contributions of the noise in
$(\widetilde{\mu}_{j})_{j=1,a,b}$, (F.6) can be further upper bounded by,

$\displaystyle
S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N+R),\mathbold{\widetilde{\mu}}(N+R))-S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N),\mathbold{\widetilde{\mu}}(N))$
$\displaystyle\leq~{}\Delta\widetilde{N}_{1}\cdot(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-\Delta\widetilde{N}_{b}\cdot
d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})+O\left(N^{1-3\alpha/8}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right).$
(105)

Now, we find an upper bound to
$S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N),\mathbold{\widetilde{\mu}}(N))$.
Since the algorithm pulls arm $a$ and doesn’t pull arm $b$ at iteration $N$,
we must have

$\mathcal{I}_{a}^{(m)}(N-1)~{}\leq~{}\mathcal{I}_{b}^{(m)}(N-1)~{}=~{}\mathcal{I}_{b}^{(m)}(N).$

Therefore,

$\displaystyle
S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N),\mathbold{\widetilde{\mu}}(N))~{}$
$\displaystyle=~{}\mathcal{I}_{a}^{(m)}(N)-\mathcal{I}_{b}^{(m)}(N)$
$\displaystyle=~{}(\mathcal{I}_{a}^{(m)}(N)-\mathcal{I}_{a}^{(m)}(N-1))+(\mathcal{I}_{a}^{(m)}(N-1)-\mathcal{I}_{b}^{(m)}(N))$
$\displaystyle\leq~{}\mathcal{I}_{a}^{(m)}(N)-\mathcal{I}_{a}^{(m)}(N-1)$
$\displaystyle=~{}\mathcal{I}_{a}(N)-\mathcal{I}_{a}(N-1)+\log\left(\frac{\widetilde{N}_{a}(N)}{\widetilde{N}_{a}(N-1)}\right)$
$\displaystyle\overset{(1)}{=}~{}\mathcal{I}_{a}(N)-\mathcal{I}_{a}(N-1)+O(1)$
$\displaystyle\overset{(2)}{=}~{}O(N^{1-3\alpha/8}),$ (106)

where (1) follows from the fact that $\widetilde{N}_{a}(N)=\Theta(N)$ for
$N\geq T_{6}$, and (2) follows using the arguments used while bounding
$\mathcal{I}_{a}(N)-\mathcal{I}_{a}(N-1)$ in (F.6).

Putting this in (F.6), we get,

$\displaystyle
S_{a,b}^{(m)}(\mathbold{\widetilde{N}}(N+R),\mathbold{\widetilde{\mu}}(N+R))$
$\displaystyle\leq\Delta\widetilde{N}_{1}\cdot(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-\Delta\widetilde{N}_{b}\cdot
d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})$
$\displaystyle+O\left(N^{1-3\alpha/8}+R(N^{-3\alpha/8}+RN^{-1})\right).$ (107)

Now there can be two cases.

Case I
$d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b})<0$
: Then, using the same argument as was used for Case I of AT2 algorithm, we
get (F.16).

Case II
$d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b})\geq
0$ : Using (88) of Lemma F.12, the RHS in (F.6) is upper bounded by,

$\Delta\widetilde{N}_{b}\left(\frac{\widetilde{N}_{1}}{\widetilde{N}_{b}}(d(\widetilde{\mu}_{1},\widetilde{x}_{1,a})-d(\widetilde{\mu}_{1},\widetilde{x}_{1,b}))-d(\widetilde{\mu}_{b},\widetilde{x}_{1,b})\right)+O\left(N^{1-3\alpha/8}+\left(N^{-3\alpha/8}+RN^{-1}\right)R\right).$

Using (F.6), we have

$\mathcal{I}_{a}^{(m)}(N)~{}\leq~{}\mathcal{I}_{b}^{(m)}(N)+O\left(N^{1-3\alpha/8}\right),$

which implies,

$\mathcal{I}_{a}(N)~{}\leq~{}\mathcal{I}_{b}(N)+\log\left(\frac{\widetilde{N}_{b}}{\widetilde{N}_{a}}\right)+O(N^{1-3\alpha/8})~{}\overset{(1)}{=}~{}\mathcal{I}_{b}(N)+O(N^{1-3\alpha/8}),$

where (1) follows from the fact that $\widetilde{N}_{a}(N)$ and
$\widetilde{N}_{b}(N)$ are both $\Theta(N)$, causing
$\log\left(\frac{\widetilde{N}_{a}(N)}{\widetilde{N}_{b}(N)}\right)=O(1)$ for
$N\geq T_{6}$. Note that the above inequality is same as (102) obtained in
Case II of AT2. Now (F.16) follows using exactly the same argument as in the
Case II of AT2 after (102). ∎

For every $a\in[K]/\\{1\\}$, $N\geq T_{6}$ and $R\in\\{1,2,\ldots,N\\}$, we
have,

$\displaystyle\left|\mathcal{I}_{a}(N+R)-\mathcal{I}_{a}(N+R-1)\right|~{}$
$\displaystyle=~{}\left|I_{a}(N+R)-I_{a}(N+R-1)\right|+O((N+R)^{1-3\alpha/8})\quad\text{(using
Lemma \ref{lem:noise_in_index})}$
$\displaystyle=~{}\left|I_{a}(N+R)-I_{a}(N+R-1)\right|+O(N^{1-3\alpha/8})\quad\text{(since
$R\leq N$)}$
$\displaystyle\overset{(1)}{\leq}~{}\max\\{d(\mu_{1},\mu_{a}),d(\mu_{a},\mu_{1})\\}+O(N^{1-3\alpha/8})~{}=~{}O(N^{1-3\alpha/8}),$
(108)

where (1) follows from the fact that $I_{a}(\cdot)$ can increment by atmost
$\max\\{d(\mu_{1},\mu_{a}),d(\mu_{a},\mu_{1})\\}$ in one iteration. Since, at
every iteration $N\geq 1$ and for every arm, the modified empirical index
differ from the empirical index by atmost $\log(N)$, we have,

$\displaystyle\left|\mathcal{I}^{(m)}_{a}(N+R)-\mathcal{I}^{(m)}_{a}(N+R-1)\right|~{}$
$\displaystyle\leq~{}\left|\mathcal{I}_{a}(N+R)-\mathcal{I}_{a}(N+R-1)\right|+2\log(N+R)$
$\displaystyle\leq~{}\left|\mathcal{I}_{a}(N+R)-\mathcal{I}_{a}(N+R-1)\right|+O(\log(N))\quad\text{(since
$R\leq N$)}$
$\displaystyle=~{}O(N^{1-3\alpha/8})+O(\log(N))~{}=~{}O(N^{1-3\alpha/8}).$
(109)

Using (F.6) and (F.6), the following corollary follows from Lemma F.16,

######  Corollary F.2.

For AT2 (1) and IAT2 (2) algorithms, for $N\geq T_{6}$ and
$R\in\\{1,2,\ldots,N\\}$, if the algorithm pulls some arm $a\in[K]/\\{1\\}$ at
iteration $N$ and doesn’t pull $a$ for the next $R$ iterations, then,

$\displaystyle\text{{AT2:}}\quad\mathcal{I}_{a}(N+R-1)-\mathcal{I}_{b(N,R)}(N+R-1)~{}$
$\displaystyle\leq~{}-C_{3}\Delta\widetilde{N}_{b(N,R)}(N,R)+C_{5}N^{1-3\alpha/8}$
$\displaystyle\quad+O\left(\left(N^{-3\alpha/8}+RN^{-1}\right)R\right),$ (110)

$\displaystyle\text{{IAT2:}}\quad\mathcal{I}^{(m)}_{a}(N+R-1)-\mathcal{I}^{(m)}_{b(N,R)}(N+R-1)~{}$
$\displaystyle\leq~{}-C_{3}\Delta\widetilde{N}_{b(N,R)}(N,R)+C_{5}N^{1-3\alpha/8}$
$\displaystyle\quad+O\left(\left(N^{-3\alpha/8}+RN^{-1}\right)R\right),$ (111)

where $C_{3},C_{5}>0$ are constants independent of the sample paths.

Proof of Lemma F.2: We prove the proposition for the AT2 algorithm. The proof
for IAT2 algorithm follows the exact same argument by replacing $\mathcal{I}$
with $\mathcal{I}^{(m)}$.

In the proof, we argue via contradiction. We show that, there exists constants
$M_{5}\geq 1$ and $C_{1}>0$ independent of the sample paths, such that for
$N\geq\max\\{M_{5},T_{6}\\}$, if the algorithm pulls some arm
$a\in[K]/\\{1\\}$ at iteration $N$ and doesn’t pull it for the next
$R(N)\overset{\text{def.}}{=}\lceil C_{1}N^{1-3\alpha/8}\rceil$ iterations,
then at iteration $\tau(N,R(N))$, the algorithm pulls arm $b(N,R(N))$, even
though arm $a$ has empirical index strictly less than that of arm $b(N,R(N))$.

Using Corollary F.2, there exists constants $C_{3},C_{5},C_{6}>0$ independent
of the sample paths, such that, for $N\geq T_{6}$ and
$R\in\\{1,2,\ldots,N\\}$,

$\displaystyle\mathcal{I}_{a}(N+R-1)-\mathcal{I}_{b(N,R)}(N+R-1)~{}$
$\displaystyle\leq~{}-C_{3}\Delta\widetilde{N}_{b(N,R)}(N,R)+C_{4}N^{1-3\alpha/8}+C_{5}R(N^{-3\alpha/8}+RN^{-1}).$
(112)

We define

$C_{1}=\max\left\\{D_{1},\frac{2C_{4}}{D_{2}C_{3}}\right\\}+1,$

where $D_{1}$ and $D_{2}$ are the constants introduced in Lemma F.15. Let
$M_{52}\geq 1$ to be the smallest number such that, every $N\geq M_{52}$
satisfies $\lceil C_{1}N^{1-3\alpha/8}\rceil<N$. Since $C_{1}>D_{1}$, by the
definition of $M_{51}$ in the proof of Lemma F.15, we have $M_{52}>M_{51}$.

Now let $R(N)=\lceil C_{1}N^{1-3\alpha/8}\rceil$. Then by Lemma F.15, for
$N\geq\max\\{M_{52},T_{6}\\}$, since $R(N)>\lceil D_{1}N^{1-3\alpha/8}\rceil$,
we have,

$\displaystyle\Delta\widetilde{N}_{b(N,R(N))}(N,R(N))~{}$
$\displaystyle\geq~{}D_{2}R(N)~{}\geq~{}\frac{2C_{4}}{C_{3}}N^{1-3\alpha/8}.$
(113)

Since arm $b(N,R)$ is not pulled between the iterations $\tau(N,R)$ and $N+R$,
we have

$\Delta\widetilde{N}_{b(N,R)}(N,t)~{}=~{}\Delta\widetilde{N}_{b(N,R)}(N,R)$

for all
$t\in\\{~{}\tau(N,R)-N,~{}~{}\tau(N,R)-N+1,~{}~{}\ldots,~{}~{}R~{}\\}$.

As a result, by definition of $b(N,R)$ we have,

$b(N,R(N))=b(N,\tau(N,R(N))-N)\quad\text{and}\quad\Delta\widetilde{N}_{b(N,R(N))}(N,R(N))=\Delta\widetilde{N}_{b(N,R)}(N,\tau(N,R(N))-N).$

In the rest of the proof, we denote $b(N,R(N))$ and $\tau(N,R(N))-N$,
respectively, using $b(N)$ and $\tau_{b}(N)$.

Therefore, using Corollary F.2, we have, for $N\geq\max\\{M_{52},T_{6}\\}$,

$\displaystyle\mathcal{I}_{a}(N+\tau_{b}(N)-1)-\mathcal{I}_{b(N)}(N+\tau_{b}(N)-1)$
$\displaystyle\leq~{}-C_{3}\Delta\widetilde{N}_{b(N)}(N,\tau_{b}(N))+C_{4}N^{1-3\alpha/8}+C_{5}\tau_{b}(N)\times(N^{-3\alpha/8}+\tau_{b}(N)N^{-1})$
$\displaystyle\overset{(1)}{\leq}~{}-C_{3}\times\frac{2C_{4}}{C_{3}}N^{1-3\alpha/8}+C_{4}N^{1-3\alpha/8}+C_{5}R(N)\times(N^{-3\alpha/8}+R(N)N^{-1})$
$\displaystyle\overset{(2)}{\leq}~{}-C_{4}N^{1-\alpha/8}+C_{5}(C_{1}+1)(C_{1}+2)N^{1-3\alpha/8}\times
N^{-3\alpha/8}$
$\displaystyle\overset{(3)}{=}~{}-(C_{4}-C_{6}N^{-3\alpha/8})N^{1-3\alpha/8},$
(114)

where: (1) follows using (113) and $\tau_{b}(N)\leq R(N)$, (2) follows since
$R(N)\leq(C_{1}+1)N^{-3\alpha/8}$, and (3) follows by letting
$C_{6}=C_{5}(C_{1}+1)(C_{1}+2)$.

We now take $M_{53}\geq 1$ to be large enough, such that,
$C_{6}M_{53}^{-3\alpha/8}<C_{4}$. Let $M_{5}=\max\\{M_{52},M_{53}\\}$. Then
(F.6) implies, for $N\geq\max\\{M_{5},T_{6}\\}$ and $R(N)=\lceil
C_{1}N^{1-3\alpha/8}\rceil$, if the algorithm picks some arm $a\in[K]/\\{1\\}$
at iteration $N$ and doesn’t pick $a$ for the next $R(N)$ iterations, then, at
iteration $N+\tau_{b}(N)$, the algorithm picks arm $b(N)$, even though,
$\mathcal{I}_{a}(N+\tau_{b}(N)-1)-\mathcal{I}_{b(N)}(N+\tau_{b}(N)-1)<0$. Thus
we arrive at a contradiction. ∎

##  Appendix G Proof of Theorem 3.1

By Proposition 3.1, we know that, for every $a\in[K]$ and $N\geq T_{stable}$,

$|\widetilde{\omega_{a}}(N)-\omega_{a}^{\star}|~{}\leq~{}C_{1}N^{-3\alpha/8}.$

Recall the functions $W_{a}(\cdot,\cdot)$ defined in Appendix F.2.2. Note that
for every allocation $\mathbold{\omega}=(\omega_{a})_{a\in[K]}$, and
$a\in[K]/\\{1\\}$, we have,

$\displaystyle\frac{\partial
W_{a}}{\partial\omega_{1}}(\omega_{1},\omega_{a})~{}$
$\displaystyle=~{}d(\mu_{1},x_{1,a}(\omega_{1},\omega_{a})),\quad\text{and}$
$\displaystyle\frac{\partial
W_{a}}{\partial\omega_{a}}(\omega_{1},\omega_{a})~{}$
$\displaystyle=~{}d(\mu_{a},x_{1,a}(\omega_{1},\omega_{a})),$
$\displaystyle\text{where}\quad x_{1,a}(\omega_{1},\omega_{a})$
$\displaystyle=\frac{\omega_{1}\mu_{1}+\omega_{a}\mu_{a}}{\omega_{1}+\omega_{a}}.$

As a result, for every $a\in[K]/\\{1\\}$, the partial derivatives of
$W_{a}(\omega_{1},\omega_{a})$ with respect to $\omega_{1}$ and $\omega_{a}$
are both $O(1)$. Therefore, using the mean value theorem, for every
$a\in[K]/\\{1\\}$ and $N\geq T_{stable}$, we have,

$\displaystyle|W_{a}(\widetilde{\omega}_{1}(N),\widetilde{\omega}_{a}(N))-W_{a}(\omega_{1}^{\star},\omega_{a}^{\star})|~{}$
$\displaystyle=~{}O(N^{-3\alpha/8}).$ (115)

Define the normalized index of every arm as

$H_{a}(N)~{}=~{}\frac{I_{a}(N)}{N}~{}=~{}W_{a}(\widetilde{\omega}_{1}(N),\widetilde{\omega}_{a}(N)).$

Also by Corollary 4.1,
$W_{a}(\omega_{1}^{\star},\omega_{a}^{\star})=I^{\star}=T^{\star}(\mathbold{\mu})^{-1}$
for every alternative arm $a\in[K]/\\{1\\}$.

Therefore (115) gives us,

$\displaystyle|H_{a}(N)-T^{\star}(\mathbold{\mu})^{-1}|~{}$
$\displaystyle=~{}O(N^{-3\alpha/8}),\quad\text{for}~{}a\in[K]/\\{1\\}~{}\text{and}~{}N\geq
T_{stable}.$ (116)

By Lemma F.3, we also know,

$\displaystyle|\mathcal{I}_{a}(N)-NH_{a}(N)|~{}$
$\displaystyle=~{}|\mathcal{I}_{a}(N)-I_{a}(N)|~{}=~{}O(N^{1-3\alpha/8}).$
(117)

Combining (116) and (117), we get

$\displaystyle|\mathcal{I}_{a}(N)-NT^{\star}(\mathbold{\mu})^{-1}|~{}$
$\displaystyle\leq~{}|\mathcal{I}_{a}(N)-NH_{a}(N)|+N|H_{a}(N)-T^{\star}(\mathbold{\mu})^{-1}|$
$\displaystyle=~{}O(N^{1-3\alpha/8}),$

for every $a\in[K]/\\{1\\}$. Hence, we can find a constant $C_{2}>0$,
independent of the sample paths, such that,

$\displaystyle\mathcal{I}_{a}(N)~{}\geq~{}\frac{N}{T^{\star}(\mathbold{\mu})}-C_{2}N^{1-3\alpha/8},$

for every $N\geq T_{stable}$ and $a\in[K]/\\{1\\}$. As a result, for $N\geq
T_{stable}$, we have,

$\displaystyle\min_{a\in[K]/\\{1\\}}\mathcal{I}_{a}(N)~{}$
$\displaystyle\geq~{}\frac{N}{T^{\star}(\mathbold{\mu})}-C_{2}N^{1-3\alpha/8}.$
(118)

The threshold function $\beta(N,\delta)$ deciding out stopping condition is of
the form,

$\beta(N,\delta)~{}=~{}\log(1/\delta)+C_{3}\log\log(1/\delta)+C_{4}\log\log(N)+C_{5}$

for constants $C_{3},C_{4},C_{5}>0$. For every $\delta>0$, we define the
deterministic quantity,

$\displaystyle t_{\max,\delta}~{}$ $\displaystyle=~{}\min\left\\{~{}N\geq
T_{stable}~{}~{}\Big{|}~{}~{}\frac{N}{T^{\star}(\mathbold{\mu})}-C_{2}N^{1-3\alpha/8}>\beta(N,\delta)~{}\right\\}.$
(119)

Now we make the following observations about $t_{\max,\delta}$,

  1. Observation 1:

Note that $\frac{N}{T^{\star}(\mathbold{\mu})}-C_{2}N^{1-3\alpha/8}$ increases
linearly in $N$ and $\beta(N,\delta)$ is $O(\log\log(N)+\log(1/\delta))$, for
a fixed $\delta>0$. Hence, $t_{\max,\delta}$ is finite for every $\delta>0$.

  2. Observation 2:

We have $\beta(N,\delta)\geq\log(1/\delta)$ and
$\frac{N}{T^{\star}(\mathbold{\mu})}-C_{2}N^{1-3\alpha/8}<\frac{N}{T^{\star}(\mathbold{\mu})}$.
As a result, $t_{\max,\delta}$ is atleast the iteration at which
$\frac{N}{T^{\star}(\mathbold{\mu})}$ exceeds $\log(1/\delta)$, which is
atleast $T^{\star}(\mathbold{\mu})\log(1/\delta)$. This implies
$t_{\max,\delta}\geq T^{\star}(\mathbold{\mu})\log(1/\delta)$. As a result,
$t_{\max,\delta}\to\infty$ as $\delta\to 0$.

  3. Observation 3:

If $\tau_{\delta}>T_{stable}$, then $\min_{a\in[K]/\\{1\\}}\mathcal{I}_{a}(N)$
exceed $\beta(N,\delta)$ before the lower bound of
$\min_{a\in[K]/\\{1\\}}\mathcal{I}_{a}(N)$ in (118) exceeds $\beta(N,\delta)$.
As a result, $\tau_{\delta}\leq t_{\max,\delta}$, whenever $\tau_{\delta}\geq
T_{stable}$. This gives us the upper bound,

$\displaystyle\tau_{\delta}~{}$
$\displaystyle\leq~{}\max\\{T_{stable},t_{\max,\delta}\\}\quad\text{a.s.
in}~{}\mathbb{P}_{\mathbold{\mu}}.$ (120)

We have $\mathbb{E}_{\mathbold{\mu}}[T_{stable}]<\infty$, which also implies
$T_{stable}<\infty$ a.s. in $\mathbb{P}_{\mathbold{\mu}}$. Now using (120), we
get,

$\limsup_{\delta\to
0}\frac{\mathbb{E}_{\mathbold{\mu}}[\tau_{\delta}]}{\log(1/\delta)}~{}\leq~{}\limsup_{\delta\to
0}\frac{\mathbb{E}_{\mathbold{\mu}}[T_{stable}]+t_{\max,\delta}}{\log(1/\delta)}~{}=~{}\limsup_{\delta\to
0}\frac{t_{\max,\delta}}{\log(1/\delta)}.$

Similarly, we have,

$\limsup_{\delta\to
0}\frac{\tau_{\delta}}{\log(1/\delta)}~{}\leq~{}\limsup_{\delta\to
0}\frac{T_{stable}+t_{\max,\delta}}{\log(1/\delta)}~{}=~{}\limsup_{\delta\to
0}\frac{t_{\max,\delta}}{\log(1/\delta)}\quad\text{a.s.
in}~{}\mathbb{P}_{\mathbold{\mu}}.$

Therefore to prove asymptotic optimality, it is sufficient to prove
$\limsup_{\delta\to 0}\frac{t_{\max,\delta}}{\log(1/\delta)}\leq
T^{\star}(\mathbold{\mu})$.

Let $s_{\max,\delta}=t_{\max,\delta}-1$. Note that $\limsup_{\delta\to
0}\frac{t_{\max,\delta}}{\log(1/\delta)}=\limsup_{\delta\to
0}\frac{s_{\max,\delta}}{\log(1/\delta)}$. By definition of $t_{\max,\delta}$,
we have

$\frac{s_{\max,\delta}}{T^{\star}(\mathbold{\mu})}-C_{2}s_{\max,\delta}^{1-3\alpha/8}~{}\leq~{}\beta(s_{\max,\delta},\delta)=\log(1/\delta)+C_{3}\log\log(1/\delta)+C_{4}\log\log(s_{\max,\delta})+C_{5}.$

After some rearrangement of terms, upon dividing both sides by
$\log(1/\delta)$, we get,

$\displaystyle\frac{1}{T^{\star}(\mathbold{\mu})}\frac{s_{\max,\delta}}{\log(1/\delta)}\left(1-C_{2}s_{\max,\delta}^{-3\alpha/8}-C_{3}s_{\max,\delta}^{-1}\log\log(s_{\max,\delta})-C_{5}s_{\max,\delta}^{-1}\right)~{}$
$\displaystyle\leq~{}1+C_{3}\frac{\log\log(1/\delta)}{\log(1/\delta)}.$

By Observation 2, $s_{\max,\delta}=t_{\max,\delta}-1\to\infty$ as $\delta\to
0$. As a result, the above inequality implies,

$\frac{1}{T^{\star}(\mathbold{\mu})}\limsup_{\delta\to
0}\frac{s_{\max,\delta}}{\log(1/\delta)}~{}\leq~{}1.$

Hence the theorem stands proved.


# Stochastic modeling and computational simulations of HBV infection dynamics

Abdallah Alsammani
Department of Mathematics, Jacksonville University
2800 University Blvd N, Jacksonville, FL 32211, Jacksonville, FL32211, USA
<EMAIL_ADDRESS>
###### 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. The results of this
research provide a robust tool to understand the complex behavior of HBV
dynamics, which offers a significant contribution to the ongoing quest for
more effective HBV control and prevention strategies.

Keywords: Stochastic dynamics, HBV, stability in probability, Euler-Maruyama,
Milstein

##  1 Introduction

Hepatitis B is a serious viral infection that affects the liver, leading to
acute and chronic conditions. The virus is commonly transmitted through
childbirth, contact with infected blood or body fluids, and unsafe injections.
In 2019, 296 million people were living 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 [18]. 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 [15, 14, 4, 16, 1, 3].

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 [12]. 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 [17, 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) [14, 4, 16]. 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 [12, 2]. 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 [17, 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 Model Formulation

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

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

$\begin{split}\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{split}$ (1)

The parameters in these equations are defined in Table 1.

Parameter | Description
---|---
$\Lambda$ | Production rate of uninfected cells $x$.
$\mu_{1}$ | Death rate of $x$-cells.
$\mu_{2}$ | Death rate of $y$-cells.
$\mu_{3}$ | Free virus cleared rate.
$\eta$ | Fraction that reduced infected rate after treatment with an antiviral drug.
$\epsilon$ | Fraction that reduced free virus rate after treatment with an antiviral drug.
$p$ | Free virus production rate from $y$-cells.
$\beta$ | Infection rate of $x$-cells by free virus $z$.
$q$ | Spontaneous cure rate of $y$-cells by non-cytolytic process.
Table 1: Parameter descriptions for the deterministic model.

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

$\begin{split}dx(t)&=(\Lambda-\mu_{1}x-(1-\eta)\beta
xz+qy)dt+\sigma_{1}xdW_{1}(t),\\\ dy(t)&=((1-\eta)\beta
xz-\mu_{2}y-qy)dt+\sigma_{2}ydW_{2}(t),\\\
dz(t)&=((1-\epsilon)py-\mu_{3}z)dt+\sigma_{3}zdW_{3}(t).\end{split}$ (2)

Here, $W_{1}(t)$, $W_{2}(t)$, and $W_{3}(t)$ represent standard Wiener
processes, and $\sigma_{1}$, $\sigma_{2}$, and $\sigma_{3}$ represent the
noise coefficients that control the amplitude of the stochastic fluctuations.

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

$du_{t}=f(u(t),t)dt+B(u,t)dW\,,\quad\quad t\geq 0$ (3)

Here, the initial condition is $u_{0}\in\mathbb{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 $dW$, are
defined as follows:

$\displaystyle f(u,t)$
$\displaystyle=\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},$
$\displaystyle B(u,t)$ $\displaystyle=\begin{pmatrix}\sigma_{1}x&0&0\\\
0&\sigma_{2}y&0\\\ 0&0&\sigma_{3}z\end{pmatrix},$ $\displaystyle dW$
$\displaystyle=\begin{pmatrix}dW_{1}\\\ dW_{2}\\\ dW_{3}\end{pmatrix}.$

###  2.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 [13], lay the
groundwork for the forthcoming analysis.

In general, equation 3 can be reformulated as a d-dimensional stochastic
equation within the context of a complete probability space
$(\Omega,\mathcal{F},\mathbb{P})$, paired with a filtration
$\\{\mathcal{F}_{t}\\}_{t\geq 0}$. This reformulation can be expressed as

$du_{t}=f(u,t)dt+B(u,t)dW\,,\quad\quad t\geq 0$ (4)

Here, $f(u,t):\mathbb{R}^{d}\times[t_{0},T]\rightarrow\mathbb{R}^{d}$ and
$B(u,t):\mathbb{R}^{d}\times[t_{0},T]\rightarrow\mathbb{R}^{d\times n}$ are
Borel measurable, and the white noise
$W(t)=(w_{1}(t),w_{2}(t),\cdots,w_{n}(t))\in\mathbb{R}^{n}$, for $t\geq 0$.
Assuming that $0\leq t_{0}\leq T<\infty$ and that the initial value $u_{0}$ is
$\mathcal{F}_{t_{0}}$-measurable $\mathbb{R}^{d}$ random variable such that
$E|u_{0}|^{2}<\infty$, equation 4 can be recognized as an Itô type stochastic
differential equation.

######  Definition 2.1.

A stochastic process $\\{u(t)\\}_{0\leq t_{0}\leq T}$ in $\mathbb{R}^{d}$ is
identified as a solution of equation 4 if it satisfies the following
conditions:

  1. i.

$\\{u(t)\\}$ is continuous and $\mathcal{F}_{t}$-adapted;

  2. ii.

$\\{f(u(t),t)\\}\in\mathcal{L}^{1}([t_{0},T];\mathbb{R}^{d})$ and
$\\{B(u(t),t)\\}\in\mathcal{L}^{2}([t_{0},T];\mathbb{R}^{3}{d\times n})$;

  3. 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:

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

######  Lemma 2.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

$du(t)=f(u(t),t)dt+B(u(t),t)dW(t),$ (5)

Where $u(t)=(x_{1}(t),x_{2}(t),\cdots,x_{d}(t))$ with the initial condition
$u(t_{0})=u_{0}\in\mathbb{R}^{d}$, and W(t) is the m-dimensional white noise
defined on a complete probability space
$(\Omega,\mathcal{F},\\{\mathcal{F}_{t}\\}_{t\geq 0},\mathbb{P})$.

###  2.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)]_{ij}\frac{\partial^{2}}{\partial
u_{i}\partial u_{j}}.$ (6)

Consider $V(u,t)$, a nonnegative twice differentiable function defined on
$\mathbb{R}^{d}\times[t_{0},\infty)$, applying operator $\mathcal{L}$ on $V$,
we get:

$LV=V_{t}+V_{u}f_{i}(u,t)+\frac{1}{2}\text{trace}[W^{T}(u,t)V_{uu}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}},\cdots\frac{\partial
V}{\partial x_{d}})$, and $V_{uu}=\left(\frac{\partial^{2}V}{\partial
u_{i}\partial u_{j}}\right)_{d\times d}$.

Therefore, the Itô formula can be defined as:

$dV=LV(u(t),t)dt+V_{u}(u(t),t)B(u(t),t)dW(t).$ (8)

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

###  2.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 [8].

######  Theorem 2.1.

For any initial value $u_{0}\in\mathbb{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\mathbb{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\mathbb{R}^{3}_{+}$ for any
initial $u_{0}\in\mathbb{R}^{3}_{+}$, where $t\in[0,\tau)$. To establish that
the solution is global, we must show that $\tau=\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

$\tau_{n}=\inf\\{t\in[0,\tau):u(t)\notin(1/n,n)\\}.$

We must show that $\tau_{n}$ is an empty set, assuming that
$\inf{\Phi}=\infty$. Since $\tau_{n}$ is an increasing sequence, let
$\tau_{\infty}=\lim_{n\rightarrow\infty}\tau_{n}$. By definition,
$\tau_{\infty}<\tau$ a.s.

To complete the proof, we must show that $\tau_{\infty}=\infty$ a.s., which
implies $\tau=\infty$. Assume by contradiction that this is not true. Then
there exists a pair of constants $T>0$ and $\epsilon\in(0,1)$ such that
$\mathbb{P}\\{\tau_{\infty}>T\\}>\epsilon$, which implies the existence of an
integer $n_{1}>n_{0}$ such that

$\mathbb{P}\\{\tau_{n}<T\\}\geq\epsilon,\quad\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 [6]). By
applying Itô’s formula on $G(u)$, we get the following:

$\displaystyle dV$
$\displaystyle=[(1-\frac{1}{x})(\lambda-\mu_{1}x-(1-\eta)\beta xz+qy)$
$\displaystyle\quad+(1-\frac{1}{y})((1-\eta)\beta xz-\mu_{2}y-qy)$
$\displaystyle\quad+(1-\frac{1}{z})((1-\epsilon)py-\mu_{3}z)+\frac{1}{2}(\sigma_{1}^{2}+\sigma_{2}^{2}+\sigma_{3}^{2})]dt$
$\displaystyle\quad+\sigma_{1}(x-1)dW_{1}+\sigma_{2}(y-1)dW_{2}+\sigma_{3}(z-1)dW_{3},$

which simplifies to

$\displaystyle dV$
$\displaystyle\leq[a+bV(t)]dt+\sigma_{1}(x-1)dW_{1}+\sigma_{2}(y-1)dW_{2}+\sigma_{3}(z-1)dW_{3},$

where
$a=\lambda+\mu_{1}+\mu_{2}+\mu_{3}+q+\frac{1}{2}(\sigma_{1}^{2}+\sigma_{2}^{2}+\sigma_{3}^{2})$,
and $b=\max\\{(1-\epsilon)p,(1-\eta)\beta\\}$.

Let $c=\max\\{a,b\\}$. Integrating the above inequality, if $t_{1}\leq T$ we
get

$\displaystyle\int_{0}^{\tau_{n}\wedge t_{1}}dV(u(t))$
$\displaystyle\leq\int_{0}^{\tau_{n}\wedge
t_{1}}c(1+V(u(t)))dt+\int_{0}^{\tau_{n}\wedge t_{1}}\sigma_{1}(x-1)dW_{1}$
$\displaystyle+\int_{0}^{\tau_{n}\wedge
t_{1}}\sigma_{2}(y-1)dW_{2}+\int_{0}^{\tau_{n}\wedge
t_{1}}\sigma_{3}(z-1)dW_{3}$

By definition, this implies

$\displaystyle EV(u(\tau_{n}\wedge t_{1}))$ $\displaystyle\leq
V(u_{0})+E\int_{0}^{\tau_{n}\wedge t_{1}}c(1+V(u(t)))dt,$ $\displaystyle\leq
V(u_{0})+ct_{1}+cE\int_{0}^{\tau_{n}\wedge t_{1}}V(u(t))dt,$
$\displaystyle\leq V(u_{0})+cT+cE\int_{0}^{t_{1}}V(\tau_{n}\wedge t_{1})dt,$
$\displaystyle=V(u_{0})+cT+c\int_{0}^{t_{1}}EV(\tau_{n}\wedge t_{1})dt.$

From the Gronwall inequality, we obtain

$EV(\tau_{n}\wedge t_{1})\leq c_{1}=(V(u_{0})+cT)e^{cT}.$ (10)

We then set $\Omega_{n}=\\{\tau_{n}\leq T\\}$ for $n\leq n_{1}$ and, by
inequality 9, $P(\Omega_{n})\geq\epsilon$. Note that there exists some $x,y,$
or $z$ such that $u(\tau_{n},\omega)=n$ or $1/n$, for every $\omega\in\Omega$.
Therefore, $V(u(\tau_{n},\omega))$ 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(\tau_{n},\omega))\geq\min\left(n-1-\ln(n),\frac{1}{n}+1+\ln(n)\right).$

From the inequalities 9 and 10, we obtain

$\displaystyle c_{1}$ $\displaystyle\geq
E[1_{\Omega_{n}}(\omega)V(u(\tau_{n},\omega))]$ $\displaystyle\geq
E\left[\min\left(n-1-\ln(n),\frac{1}{n}+1+\ln(n)\right)\right],$

where $1_{\Omega_{n}}$ is the indicator function of $\Omega_{n}$. Taking the
limit as $n$ approaches $\infty$ yields $c_{1}=\infty$, which is a
contradiction. Therefore, $\tau_{\infty}=\infty$ a.s., which completes the
proof. ∎

##  3 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, please refer to [13] and [5].

######  Definition 3.1 ([13], pages 110, 119).

  * (i)

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

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

whenever $|u_{0}|<\delta$. If not, it is termed stochastically unstable.

  * (ii)

The trivial solution is considered stochastically asymptotically stable if it
is stochastically stable and, for each $\epsilon\in(0,1)$ and $r>0$, there
exists a $\delta=\delta(\epsilon,r,t_{0})>0$ such that

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

whenever $|u_{0}|<\delta$.

  * (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\mathbb{R}^{d}$ and $r>0$, there is a
$\delta=\delta(\epsilon,r,t_{0})>0$ such that

$\mathbb{P}\\{\lim_{t\rightarrow\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\rightarrow\infty}\sup\frac{1}{t}\ln|u(t;t_{0},u_{0})|<0$

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

Observe that the first equation of our system, 2, does not possess a direct
equilibrium point in $\mathbb{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
$x(t)$ equation 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 3.1.

In the system 2, $y(t)$ and $z(t)$ are almost surely exponentially stable if
and only if the following conditions are met

  1. a.

$(1-\epsilon)p-q-\mu_{2}+\frac{1}{2}\sigma_{2}^{2}<0$

  2. b.

$[(1-\eta)\beta\gamma-\mu_{3}][(1-\epsilon)p-q-\mu_{2}]\leq[1-\epsilon)p-q-\mu_{2}+\frac{1}{2}\sigma_{2}^{2}][(1-\eta)\beta\gamma-\mu_{3}+\frac{1}{2}\sigma_{3}^{2}]$

where $\gamma=\max\\{x\\}$.

The proof proceeds as follows

###### Proof.

By adding equations $y(t)$ and $z(t)$ equations, we get

$d(y+z)=\left[((1-\epsilon)p-\mu_{2}-q)y+((1-\eta)\beta
x-\mu_{3})z\right]dt+\sigma_{2}ydW_{2}+\sigma_{3}zdW_{3},$

Let $V(y,z)=\ln(y(t)+z(t))$ for $y,z\in\mathbb{R}_{+}$. Applying Itô’s
formula, we have

$\displaystyle dV$
$\displaystyle=\left[\frac{1}{y+z}\left((1-\epsilon)p-q-\mu_{2}\right)y\right]+\left[\frac{1}{y+z}\left((1-\eta)\beta
x-\mu_{3}\right)z\right]$
$\displaystyle\quad+\frac{1}{2}\frac{\sigma_{2}^{2}y^{2}}{(y+z)^{2}}+\frac{1}{2}\frac{\sigma_{3}^{2}z^{2}}{(y+z)^{2}}\,dt$
$\displaystyle\quad+\frac{y}{y+z}\sigma_{2}dW_{2}+\frac{z}{y+z}\sigma_{3}dW_{3}$
$\displaystyle=\frac{1}{(y+z)^{2}}\left[(y+z)\left((1-\epsilon)p-q-\mu_{2}\right)y+\left((1-\eta)\beta
x-\mu_{3}\right)z\right]$
$\displaystyle\quad+\frac{1}{2}\sigma_{2}^{2}y^{2}+\frac{1}{2}\sigma_{3}^{2}z^{2}$
$\displaystyle\quad+\frac{y}{y+z}\sigma_{2}dW_{2}+\frac{z}{y+z}\sigma_{3}dW_{3}$
$\displaystyle\leq\frac{1}{(y+z)^{2}}\left[(y+z)\left((1-\epsilon)p-q-\mu_{2}\right)y+\left((1-\eta)\beta\gamma-\mu_{3}\right)z\right]$
$\displaystyle\quad+\frac{1}{2}\sigma_{2}^{2}y^{2}+\frac{1}{2}\sigma_{3}^{2}z^{2}$
$\displaystyle\quad+\frac{y}{y+z}\sigma_{2}dW_{2}+\frac{z}{y+z}\sigma_{3}dW_{3},$
$\displaystyle=\frac{1}{(y+z)^{2}}\left\\{\begin{pmatrix}y&z\end{pmatrix}\begin{pmatrix}(1-\epsilon)p-q-\mu_{2}+\frac{1}{2}\sigma_{2}^{2}&(1-\eta)\beta\gamma-\mu_{3}\\\
(1-\epsilon)p-q-\mu_{2}&(1-\eta)\beta\gamma-\mu_{3}+\frac{1}{2}\sigma_{3}^{2}\end{pmatrix}\begin{pmatrix}y\\\
z\end{pmatrix}\right\\}$
$\displaystyle\quad+\frac{y}{y+z}\sigma_{2}dW_{2}+\frac{z}{y+z}\sigma_{3}dW_{3}.$

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

$\begin{pmatrix}(1-\epsilon)p-q-\mu_{2}+\frac{1}{2}\sigma_{2}^{2}&(1-\eta)\beta\gamma-\mu_{3}\\\
(1-\epsilon)p-q-\mu_{2}&(1-\eta)\beta\gamma-\mu_{3}+\frac{1}{2}\sigma_{3}^{2}\end{pmatrix}$

is negative-definite, which implies it has at least one negative eigenvalue.
Let $\lambda_{\max}$ denote the largest eigenvalue. With this, the inequality
above can be expressed as:

$\displaystyle dV(y,z)$
$\displaystyle\leq-|\lambda_{\max}|\frac{1}{(y+z)^{2}}(y^{2}+z^{2})dt+\frac{\sigma_{2}y}{y+z}dW_{2}+\frac{\sigma_{3}z}{y+z}dW_{3}.$

By utilizing the inequality $y^{2}+z^{2}\geq 2yz$, it follows that
$(y^{2}+z^{2})/(y+z)^{2}\geq 1/2$, yielding:

$dV=d\ln(y(t)+z(t))\leq-\frac{1}{2}|\lambda_{\max}|dt+\frac{\sigma_{2}y}{y+z}dW_{2}+\frac{\sigma_{3}z}{y+z}dW_{3}.$

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

$\limsup_{t\rightarrow\infty}\frac{1}{t}|W_{i}(t)|=0\quad\text{for }i=2,3,$

we derive

$\limsup_{t\rightarrow\infty}\frac{1}{t}\ln(y(t)+z(t))\leq-\frac{1}{2}|\lambda_{\max}|<0\quad\text{almost
surely}.$

This concludes the proof. ∎

######  Remark 3.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$ or
$R_{0}>1$. Moreover, it is crucial to note that the conditions within theorem
3.1 cannot hold in the deterministic case when both $\sigma_{2}=\sigma_{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 $\lambda/\mu_{1}$. However, before proceeding, let’s first
introduce some necessary lemmas.

######  Lemma 3.1.

Assume $W_{1}(t)$ is a one-dimensional standard Brownian motion. Then, the
expectation is given by

$E\\{e^{\sigma_{1}(W_{1}(t)-W_{1}(s))}\\}=e^{\frac{\sigma_{1}^{2}}{2}(t-s)},\quad\text{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,

$\displaystyle E\\{e^{\sigma_{1}W}\\}$
$\displaystyle=E\\{e^{\sigma_{1}(W_{1}(t)-W_{1}(s))}\\}$
$\displaystyle=\frac{1}{\sqrt{2\pi(t-s)}}\int_{-\infty}^{\infty}e^{\sigma_{1}w}\cdot
e^{-\frac{w^{2}}{2(t-s)}}\,dw$
$\displaystyle=\frac{1}{\sqrt{2\pi(t-s)}}\int_{-\infty}^{\infty}e^{-\frac{(w+\sigma_{1}(t-s))^{2}}{2(t-s)}}\cdot
e^{\frac{\sigma_{1}^{2}}{2}(t-s)}\,dw$
$\displaystyle=e^{\frac{\sigma_{1}^{2}}{2}(t-s)}.$

∎

######  Lemma 3.2.

Assume that $x_{1}(t)$ is a solution to

$\mathrm{d}x_{1}(t)=(\lambda-\mu_{1}x_{1}(t))\mathrm{d}t+\sigma_{1}x_{1}(t)\mathrm{d}W_{1}(t)$
(11)

Then, we have $\lim_{t\rightarrow\infty}\mathbb{E}[x_{1}(t)]=\lambda/\mu_{1}$
for any initial $x_{1}(0)\in\mathbb{R}_{+}$.

###### Proof.

Given an initial value $x_{1}(0)\in\mathbb{R}_{+}$, there is a unique solution
$x_{1}(t)$ to equation (11) with the explicit form

$x_{1}(t)=x_{1}(0)e^{-(\mu_{1}+\frac{1}{2}\sigma_{1}^{2})t}+\lambda\int_{0}^{t}e^{-(\mu_{1}+\frac{1}{2}\sigma_{1}^{2})(t-s)}\cdot
e^{\sigma_{1}(W_{1}(t)-W_{1}(s))}\,ds$

Applying the expectation to the above equation and using the fact that
$W_{1}(0)=0$, we get

$\displaystyle\mathbb{E}[x_{1}(t)]$
$\displaystyle=\mathbb{E}\left[x_{1}(0)e^{-(\mu_{1}+\frac{1}{2}\sigma_{1}^{2})t}+\lambda\int_{0}^{t}e^{-(\mu_{1}+\frac{1}{2}\sigma_{1}^{2})(t-s)}\cdot
e^{\sigma_{1}(W_{1}(t)-W_{1}(s))}\,ds\right]$
$\displaystyle=x_{1}(0)e^{-\mu_{1}t}+\frac{\lambda}{\mu_{1}}(1-e^{-\mu_{1}t})$

By applying lemma 3.1, we find

$\lim_{t\rightarrow\infty}\mathbb{E}[x_{1}(t)]=\frac{\lambda}{\mu_{1}}$

∎

######  Lemma 3.3.

Assume that $x_{1}(t)$ is a solution of equation 11. Then, for any initial
value $x_{1}(0)\in\mathbb{R}_{+}$, we obtain the following results:

  * i.

$x_{1}(t)$ admits a unique stationary distribution $\pi$, where
$\Gamma(\cdot)$ represents the Gamma function.

  * ii.

$x_{1}(t)$ satisfies the equation

$\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{t}x_{1}(s)ds=\int_{0}^{\infty}x_{1}\pi
dx_{1}=\int_{0}^{\infty}f(x_{1})dx_{1}=\frac{\lambda}{\mu_{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

$LV(x_{1})=(1-\frac{1}{x_{1}})(\lambda-\mu_{1}x_{1})+\frac{1}{2}\sigma_{1}^{2}=-\mu-1-\frac{\lambda}{x_{1}}+\lambda+\mu_{1}+\frac{1}{2}\sigma_{1}^{2}.$

By choosing a sufficiently small $\epsilon$ and defining
$\mathbb{D}=(\epsilon,1/\epsilon)$, we find that

$LV(x_{1})\leq-1,\quad\text{for any}\quad x_{1}\in\mathbb{D}^{c}.$

This completes the first part of the proof.

  * ii.

Employing the ergodicity of $x_{1}$, we obtain

$\mathbb{P}\left\\{\lim_{t\rightarrow-\infty}\frac{1}{t}\int_{0}^{t}x_{1}(s)ds=\int_{0}^{\infty}x_{1}\pi(dx_{1})\right\\}=1.$

Therefore,

$\lim_{t\rightarrow-\infty}\frac{1}{t}\int_{0}^{t}x_{1}(s)ds=\int_{0}^{\infty}x_{1}\pi(dx_{1})=\int_{0}^{\infty}f(x_{1})x_{1}dx_{1}=\frac{\lambda}{\mu_{1}},\quad\text{almost
surely.}$

where we used the fact $\lim_{t\rightarrow\infty}E[x_{1}(t)]=\lambda/\mu_{1}.$

∎

######  Theorem 3.2.

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 3.1, we can
confirm that

$\lim_{t\rightarrow\infty}[x(t)-x_{1}(t)]=0\quad\text{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
$\liminf_{t\rightarrow\infty}[x(t)-x_{1}(t)]\geq 0$ a.s.

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

$dx_{r}(t)=[\lambda-(\mu_{1}+r)x_{r}(t)]dt+\sigma_{1}x_{r}(t)dW_{1}(t),$ (14)

where the initial condition is $x_{r}(0)=x(0)$.

Recall $x(t)$ in the system 2,

$\displaystyle dx(t)$ $\displaystyle=[\lambda-\mu_{1}x-(1-\eta)\beta
xz+qy]dt+\sigma_{1}xdW_{1}(t),$

and equation 11,

$\displaystyle dx_{1}(t)$
$\displaystyle=(\lambda-\mu_{1}x_{1}(t))dt+\sigma_{1}x_{1}(t)dW_{1}(t).$

From fact that

$\displaystyle\liminf_{t\rightarrow\infty}[x(t)-x_{1}(t)]$
$\displaystyle=\liminf_{t\rightarrow\infty}(x(t)-x_{r}(t))+\liminf_{t\rightarrow\infty}(x_{r}(t)-x_{1}(t)]$
$\displaystyle\geq\liminf_{t\rightarrow\infty}[x(t)-x_{r}(t)]+\liminf_{t\rightarrow\infty}[x_{r}(t)-x_{1}(t)],$

We will proceed to prove the following claims.

#### claim 1:

$\liminf_{t\rightarrow\infty}[x(t)-x_{r}(t)]\geq 0$ almost surely.

###### Proof.

Subtracting the given equations, we find

$\displaystyle d(x(t)-x_{r}(t))$
$\displaystyle=[-\mu_{1}(x-x_{r})+rx_{r}-(1-\eta)\beta
xz+qy]dt+\sigma-1(x-x_{r})dW_{1}$
$\displaystyle=\left[-(\mu_{1}+r)(x-x_{r})+(r-(1-\eta)\beta
z)x+qy\right]dt+\sigma_{1}(x-x_{r})dW_{1}.$

This yields a solution.

$x(t)-x_{r}(t)=\phi(t)\int_{0}^{t}\phi^{-1}(s)((r-(1-\eta)\beta z)x+qy)dx(s),$

where

$\phi(t)=e^{-(\mu_{1}+r+\frac{1}{2}\sigma_{1}^{2})t+\sigma_{1}W_{1}(t)}.$

By Theorem 3.1, we know that $y(t)\rightarrow 0$ and $z(t)\rightarrow 0$
almost surely as $t\rightarrow\infty$. Therefore,

$\displaystyle x(t)-x_{r}(t)$
$\displaystyle=\phi(t)(\int_{0}^{T}\phi^{-1}(s)((r-(1-\eta)\beta z)x+qy)dx(s)$
$\displaystyle+\int_{T}^{t}\phi^{-1}(s)((r-(1-\eta)\beta z)x+qy)dx(s)),$

for all $\omega\in\Omega$ and $t>T$. Therefore,

$x(t)-x_{r}(t)\geq\phi(t)\kappa(T),$

where

$\kappa(T)=\int_{0}^{T}\phi^{-1}(s)((\epsilon-(1-\eta)\beta
z)x(s)+qy(s))dx(s).$

Given $|\kappa(T)|<\infty$ and $\phi(t)\rightarrow 0$ almost surely, we
conclude

$\liminf_{t\rightarrow\infty}[x(t)-x_{r}(t)]\geq 0\quad\quad\text{a.s.}$ (15)

∎

#### claim 2

$\liminf_{t\rightarrow\infty}[x_{r}(t)-x_{1}(t)]\geq 0\quad$ a.s.

###### Proof.

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

$\displaystyle d(x_{r}(t)-x_{1}(t))$
$\displaystyle=[-\mu_{1}(x_{r}(t)-x_{1}(t))-rx_{r}]dt+\sigma_{1}(x_{r}(t)-x_{1}(t))dW(t),$

which leads to the solution

$\displaystyle x_{r}(t)-x_{1}(t)$
$\displaystyle=-r\int_{0}^{t}x_{r}e^{-(\mu_{1}+\sigma_{1}^{2}/2)(t-s)-\sigma_{1}(W_{1}(t)-W_{1}(s))}ds.$

analog cognize $x_{r}$ as the solution for equation 14, which can be expressed
explicitly as

$\displaystyle
x_{r}=\lambda\int_{0}^{t}e^{-(\mu_{1}+r+\sigma_{1}^{2}/2)(t-s)+\sigma_{1}(W_{1}(t)-W_{1}(s))}ds.$

Hence,

$\displaystyle|x_{r}(t)-x_{1}(t)|$
$\displaystyle=r\int_{0}^{t}x_{r}e^{-(\mu_{1}+\sigma_{1}^{2}/2)(t-s)-\sigma_{1}(W_{1}(t)-W_{1}(s))}ds.$

Applying the expectation and invoking lemma 3.1, we find

$\displaystyle E|x_{r}(t)-x_{1}(t)|$
$\displaystyle=rE\left[\int_{0}^{t}x_{r}e^{-(\mu_{1}+\sigma_{1}^{2}/2)(t-s)-\sigma_{1}(W_{1}(t)-W_{1}(s))}ds\right]$
$\displaystyle=r\int_{0}^{t}E\left[x_{r}e^{-(\mu_{1}+\sigma_{1}^{2}/2)(t-s)}\right]\cdot
E\left[e^{-\sigma_{1}(W_{1}(t)-W_{1}(s))}ds\right]$
$\displaystyle=r\int_{0}^{t}e^{-\mu_{1}(t-s)}E[x_{r}]ds,\quad\text{(by lemma
}\ref{lem4}).$

Additionally, we know

$\displaystyle E[x_{r}]$
$\displaystyle=\lambda\int_{0}^{t}e^{-(\mu_{1}+r)(t-s)}ds\leq\frac{\lambda}{\mu_{1}+r},$

which implies

$\displaystyle E|x_{r}(t)-x_{1}(t)|$ $\displaystyle\leq\frac{\lambda
re^{-\mu_{1}t}}{\mu_{1}+r}(e^{\mu_{1}t}-1)$ $\displaystyle\lim_{r}\rightarrow
0\lim_{t}\rightarrow E|x_{r}(t)-x_{1}(t)|$ $\displaystyle=0.$

Therefore,

$\lim_{r\rightarrow
0}\lim_{t\rightarrow\infty}|x_{r}(t)-x_{1}(t)|=0\quad\text{(In probability).}$
(16)

∎

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

###  3.1 Existence of Ergodic Stationary Distribution

Let $U(t)$ be a Markov process in $\mathbb{R}^{d}$ represented by the
following stochastic differential equations:

$dU(t)=f(U(t))dt+\sum_{k=1}^{n}B_{k}(U(t))dW_{k}(t).$ (17)

The diffusion matrix is defined as

$A(u)=\left(a_{ij}(u)\right),\quad\text{where}\quad
a_{ij}(u)=\sum_{k=1}^{n}B_{k}^{i}(u)B_{k}^{j}(u).$ (18)

######  Lemma 3.4.

(See [11, 20]) The model in Equation 17 is positive recurrent if there exists
a boundary open subset $D\subset\mathbb{R}^{d}$ with a regular boundary and
the following conditions hold:

  1. (A1)

There is a positive number $M$ such that

$\sum_{i,j=1}^{d}a_{ij}(u)\xi_{i}\xi_{j}\geq|\xi|^{2},\quad u\in D\text{ and
}\xi\in\mathbb{R}^{d}.$ (19)

  2. (A2)

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

$\mathbb{P}\left\\{\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}f(U(t))dt=\int_{\mathbb{R}^{d}}f(u)\mu(du)\right\\}=1,$
(20)

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

######  Theorem 3.3.

Assume that $R_{0}>1$. Under the conditions
$\sigma_{1}^{2}\geq\mu_{1}-\mu^{*}$,
$\sigma_{2}^{2}\geq\mu_{2}-\mu^{*}+\frac{1}{2}(1-\epsilon)p$, and
$\mu_{3}\geq\mu_{3}-\frac{1}{2}(1-\epsilon)p$, the system 2 has a unique
ergodic stationary distribution $\pi(\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\mathbb{R}_{+}^{3}$.

The diffusion matrix of the system is given by

$A=\left(\begin{matrix}\sigma_{1}^{2}x^{2}&0&0\\\ 0&\sigma_{2}^{2}y^{2}&0\\\
0&0&\sigma_{3}^{2}z^{2}\end{matrix}\right)$

Then

$\sum_{i,j=1}^{3}a_{ij}(u)\xi_{i}\xi_{j}=\sigma_{1}^{2}x^{2}\xi_{1}^{2}+\sigma_{2}^{2}y^{2}\xi_{2}^{2}+\sigma_{3}^{2}z^{2}\xi_{3}^{2}\geq
M|\xi|^{2},$ (21)

where,
$M=\min\\{\sigma_{1}^{2}x^{2},\sigma_{2}^{2}y^{2},\sigma_{3}^{2}z^{2}\\}$,
thus condition $(A1)$ satisfied.

Now, we want to show that condition $(A2)$ is also satisfied by constructing a
nonnnegative Lyapunov function $V$ such that $LV<0$.

Consider the positive functions

$V_{1}(x,y)=\frac{1}{2}(x-\bar{x}+y-\bar{y})^{2},\quad\text{and}\quad
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

$LV_{1}(x,y)=(x-\bar{x}+y-\bar{y})(\lambda-\mu_{1}x-\mu_{2}y)+\frac{1}{2}\sigma_{1}^{2}x^{2}+\frac{1}{2}\sigma_{2}^{2}y^{2}$

since we have
$\lambda-\mu_{1}\bar{x}-\mu_{2}\bar{y}=0,\quad\Rightarrow\quad\lambda=\mu_{1}\bar{x}+\mu_{2}\bar{y}$,
then we get

$\displaystyle LV_{1}(x,y)$
$\displaystyle=(x-\bar{x}+y-\bar{y})[-\mu_{1}(x-\bar{x})-\mu_{2}(y-\bar{y})]+\frac{1}{2}\sigma_{1}^{2}x^{2}+\frac{1}{2}\sigma_{2}^{2}y^{2},$
$\displaystyle=-\mu_{1}(x-\bar{x})^{2}-\mu_{2}(y-\bar{y})^{2}-\mu_{1}(x-\bar{x})(y-\bar{y})-\mu_{2}(x-\bar{x})(y-\bar{y})$
$\displaystyle+\frac{1}{2}\sigma_{1}^{2}x^{2}+\frac{1}{2}\sigma_{2}^{2}y^{2}$

By using the fact that, $x^{2}\leq 2(x-a)^{2}+2a^{2}$, we get

$\frac{1}{2}\sigma_{1}^{2}x^{2}\leq\sigma_{1}^{2}(x-\bar{x})^{2}+\sigma_{1}^{2}\bar{x}^{2}\quad\text{and}\quad\frac{1}{2}\sigma_{2}^{2}y^{2}\leq\sigma_{2}^{2}(y-\bar{y})^{2}+\sigma_{2}^{2}\bar{y}^{2}$

thus

$\displaystyle LV_{1}(x,y)$
$\displaystyle\leq-\mu_{1}(x-\bar{x})^{2}-\mu_{2}(y-\bar{y})^{2}-2\mu^{*}(x-\bar{x})(y-\bar{y})+\sigma_{1}^{2}(x-\bar{x})^{2}+\sigma_{1}^{2}\bar{x}^{2}$
$\displaystyle+\sigma_{2}^{2}(y-\bar{y})^{2}+\sigma_{2}^{2}\bar{y}^{2}$

where $\mu^{*}=\min\\{\mu_{1},\mu_{2}\\}$, and since
$-2\mu^{*}(x-\bar{x})(y-\bar{y})\leq\mu*(x-\bar{x})^{2}+\mu^{*}(y-\bar{y})^{2}$,
then

$\displaystyle LV_{1}(x,y)$ $\displaystyle\leq$
$\displaystyle-(\mu_{1}-\mu^{*}-\sigma_{1}^{2})(x-\bar{x})^{2}-(\mu_{2}-\mu^{*}-\sigma_{2}^{2})(y-\bar{y})^{2}$
$\displaystyle+\sigma_{1}^{2}\bar{x}+\sigma_{2}^{2}\bar{y}$

Similarly,

$LV_{2}(z)\leq-(\mu_{3}-\frac{1}{2}(1-\epsilon)p-\sigma_{3}^{2})(z-\bar{z})^{2}+\frac{1}{2}(1-\epsilon)p(y-\bar{y})^{2}+\sigma_{3}^{2}\bar{z}^{2}$

Thus,

$\displaystyle LV(x,y,z)$ $\displaystyle=LV_{1}(x,y)+LV_{2}(z)$
$\displaystyle\leq-(\mu_{1}-\mu^{*}-\sigma_{1}^{2})(x-\bar{x})^{2}-(\mu_{2}-\mu^{*}-\sigma_{2}^{2})(y-\bar{y})^{2}+\sigma_{1}^{2}\bar{x}+\sigma_{2}^{2}\bar{y}$
$\displaystyle-(\mu_{3}-\frac{1}{2}(1-\epsilon)p-\sigma_{3}^{2})(z-\bar{z})^{2}+\frac{1}{2}(1-\epsilon)p(y-\bar{y})^{2}+\sigma_{3}^{2}\bar{z}^{2}$
$\displaystyle=-(\mu_{1}-\mu^{*}-\sigma_{1}^{2})(x-\bar{x})^{2}-(\mu_{2}-\mu^{*}+\frac{1}{2}(1-\epsilon)p-\sigma_{2}^{2})(y-\bar{y})^{2}$
$\displaystyle-(\mu_{3}-\frac{1}{2}(1-\epsilon)p-\sigma_{3}^{2})(z-\bar{z})^{2}+\sigma_{1}^{2}\bar{x}+\sigma_{2}^{2}\bar{y}+\sigma_{3}^{2}\bar{z}^{2}$
$\displaystyle=-k_{1}(x-\bar{x})^{2}-k_{2}(y-\bar{y})^{2}-k_{3}(z-\bar{z})^{2}+\omega$

where $k_{1}=\mu_{1}-\mu^{*}-\sigma_{1}^{2}$,
$k_{2}=\mu_{2}-\mu^{*}+\frac{1}{2}(1-\epsilon)p-\sigma_{2}^{2}$,
$k_{3}=\mu_{3}-\frac{1}{2}(1-\epsilon)p-\sigma_{3}^{2}$, and
$\omega=\sigma_{1}^{2}\bar{x}+\sigma_{2}^{2}\bar{y}+\sigma_{3}^{2}\bar{z}^{2}$

Since $k_{1},k_{2},k_{3}$ are positive and by the same computation as in [19],
we obtained

$\limsup_{t\rightarrow\infty}\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}]ds\leq\omega.$

then the ellipsoid

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

lies entirely in $\mathbb{R}_{+}^{3}$, so we can take any neighborhood $D$ of
this ellipsoid, such that

$LV(u)<-\theta,\quad\text{for all}\quad u\in D^{c}.$

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

##  4 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 [10]. 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 [9].

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.

###  4.1 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:

  1. 1.

Discretize the time domain $[0,T]$ into $n$ intervals of size $\Delta t=T/n$,
where $t_{i}=i\Delta t$. Denote $x(t_{i})$, $y(t_{i})$, and $z(t_{i})$ by
$x_{i}$, $y_{i}$, and $z_{i}$ respectively.

  2. 2.

Initialize $x_{0}$, $y_{0}$ and $z_{0}$ as the initial conditions.

  3. 3.

For each time step $i=0,1,...,n-1$, compute the increments $\Delta
W_{1}=W_{1}(t_{i+1})-W_{1}(t_{i})$ , $\Delta
W_{2}=W_{2}(t_{i+1})-W_{2}(t_{i})$, and $\Delta
W_{3}=W_{3}(t_{i+1})-W_{3}(t_{i})$

  4. 4.

Update the solution from time $t_{i}$ to $t_{i+1}$ as follows

$\displaystyle x_{i+1}$
$\displaystyle=x_{i}+(\Lambda-\mu_{1}x_{i}-(1-\eta)\beta
x_{i}z_{i}+qy_{i})\Delta t+\sigma_{1}x_{i}\Delta W_{1},$ $\displaystyle
y_{i+1}$ $\displaystyle=y_{i}+((1-\eta)\beta
x_{i}z_{i}-\mu_{2}y_{i}-qy_{i})\Delta t+\sigma_{2}y_{i}\Delta W_{2},$
$\displaystyle z_{i+1}$
$\displaystyle=z_{i}+((1-\epsilon)py_{i}-\mu_{3}z_{i})\Delta
t+\sigma_{3}z_{i}\Delta W_{3}.$

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 $\Delta W_{1}$, $\Delta W_{2}$, and $\Delta
W_{3}$, normally distributed random numbers with mean 0 and variance $\Delta
t$.

###  4.2 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:

  1. 1.

As in the Euler-Maruyama method, discretize the time domain $[0,T]$ into $n$
intervals of size $\Delta t=T/n$ and initialize $x_{0}$, $y_{0}$, and $z_{0}$
as the initial conditions.

  2. 2.

For each time step $i=0,1,...,n-1$, calculate the increments $\Delta W_{1}$,
$\Delta W_{2}$, and $\Delta W_{3}$ in the same way as the Euler-Maruyama
method.

  3. 3.

Update the solution from time $t_{i}$ to $t_{i+1}$ by

$\displaystyle x_{i+1}$
$\displaystyle=x_{i}+(\Lambda-\mu_{1}x_{i}-(1-\eta)\beta
x_{i}z_{i}+qy_{i})\Delta t+\sigma_{1}x_{i}\Delta
W_{1}+\frac{1}{2}\sigma_{1}^{2}x_{i}(\Delta W_{1}^{2}-\Delta t),$
$\displaystyle y_{i+1}$ $\displaystyle=y_{i}+((1-\eta)\beta
x_{i}z_{i}-\mu_{2}y_{i}-qy_{i})\Delta t+\sigma_{2}y_{i}\Delta
W_{2}+\frac{1}{2}\sigma_{2}^{2}y_{i}(\Delta W_{2}^{2}-\Delta t),$
$\displaystyle z_{i+1}$
$\displaystyle=z_{i}+((1-\epsilon)py_{i}-\mu_{3}z_{i})\Delta
t+\sigma_{3}z_{i}\Delta W_{3}+\frac{1}{2}\sigma_{3}^{2}z_{i}(\Delta
W_{3}^{2}-\Delta t).$

This approach takes into account the second moment of the diffusion term,
thereby providing more accurate approximations to 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

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 2: Comparison between Euler-Maruyama and Milstein Methods Parameter | $\Lambda$ | $\mu_{1}$ | $\mu_{2}$ | $\mu_{3}$ | $\beta$ | $\eta$ | $\epsilon$ | $p$ | $q$ | $\sigma_{1}$ | $\sigma_{2}$ | $\sigma_{3}$
---|---|---|---|---|---|---|---|---|---|---|---|---
Value | 100 | 20 | 5 | 7 | 0.6 | 0.6 | 0.2 | 2 | 5 | 0.5 | 0.6 | 0.8
Table 3: Parameter values for the stochastic and deterministic models Figure
1: Figure describes the solution of the model2 using Euler-Maruyama Method.
Parameter values are defined in the table3 Figure 2: Figure describes the
solution of the model 2 using Milstein Method. Parameter values are defined in
the table3 Figure 3: Figure describes the deterministic solution of the model
2. Parameter values are defined in the table3 Figure 4: Figure describes the
comparison of the stochastic ”Euler-Maruyama” and deterministic solutions.
Parameter values are defined in the table3 Figure 5: Figure describes the
comparison of the stochastic ” Milstein” and deterministic solutions.
Parameter values are defined in the table3 Figure 6: Figure describes the
comparison of the stochastic ” Milstein” and deterministic solutions.
Parameter values are defined in the table3 Figure 7: Figure describes the
comparison of the stochastic ” Milstein” and deterministic solutions.
Parameter values are defined in the table3 Figure 8: Figure describes the
comparison of the stochastic ” Milstein” and deterministic solutions.
Parameter values are defined in the table3

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.

##  5 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 nature 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 engages 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 advancement in
the field of HBV infection dynamics. It offers a deeper comprehension of the
subject and fosters 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.

## References

  * [1] I. Abdulrashid, A. Alsammani, and X. Han.  Stability analysis of chemotherapy model with delays.  Discrete and Continuous Dynamical Systems Series B, 24, 2019.
  * [2] Abdallah Alsammani.  Mathematical analysis of autonomous and nonautonomous hepatitis b virus transmission models.  In Computational Science and Its Applications – ICCSA 2023 Workshops, pages 327–343, Cham, 2023. Springer Nature Switzerland.
  * [3] Abdallah Alhadi Mahadi Alsammani.  Dynamical Behavior of Nonautonomous and Stochastic HBV Infection Model.  PhD thesis, Auburn University, 2020.
  * [4] Charles R Bangham.  Population dynamics of immune responses to persistent viruses.  Science, 272(5258):74–79, 1996.
  * [5] Tomas Caraballo and Xiaoying Han.  Applied Nonautonomous and Random Dynamical Systems.  SpringerBriefs, 2016.
  * [6] N. Dalal, D. Greenhalgh, and X. Mao.  A stochastic model for internal hiv dynamics.  Journal of Mathematical Analysis and Applications, 2008.
  * [7] Jeremie Guedj and Alan S Perelson.  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, 110(10):3991–3996, 2013.
  * [8] R. Z. Has’minskii.  Stochastic Stability of Differential Equations.  Sijthoff Noordhoff, Alphen aan den Rijn, The Netherlands, 1980.
  * [9] Desmond J Higham.  Algorithmic introduction to numerical simulation of stochastic differential equations.  SIAM review, 43(3):525–546, 2001.
  * [10] Peter E Kloeden and Eckhard Platen.  Numerical solution of stochastic differential equations.  Springer Science & Business Media, 2013.
  * [11] Dan Li, Jing’an Cui, Meng Liu, and Shengqiang Liu.  The evolutionary dynamics of stochastic epidemic model with nonlinear incidence rate.  Bulletin of mathematical biology, 77:1705–1743, 2015.
  * [12] Alun L Lloyd.  Realistic distributions of infectious periods in epidemic models: Changing patterns of persistence and dynamics.  Theoretical Population Biology, 60(1):59–71, 2001.
  * [13] X. Mao.  Stochastic Differential Equations and Applications.  Horwood Publishing Limited, 2 edition, 2007.
  * [14] Martin A Nowak et al.  Viral dynamics of primary viremia and antiretroviral therapy in simian immunodeficiency virus infection.  Journal of virology, 70(10):6735–6741, 1996.
  * [15] R.J.H. Payne, M.A. Nowak, and B. Blumberg.  The dynamics of hepatitis b virus infection.  Proc. Natl. Acad. Sci. USA, 93:6542–6546, 1996.
  * [16] Alan S Perelson and Patrick W Nelson.  Modelling viral and immune system dynamics.  Nature Reviews Immunology, 2(1):28–36, 2002.
  * [17] Andre S Ribeiro.  Effects of stochastic population dynamics on the quantification of gene expression.  Physical Review E, 71(1):011912, 2005.
  * [18] World Health Organization.  Hepatitis b, 2023.  Accessed: 2023-07-23.
  * [19] Qingshan Yang, Daqing Jiang, Ningzhong Shi, and Chunyan Ji.  The ergodicity and extinction of stochastically perturbed sir and seir epidemic models with saturated incidence.  Journal of Mathematical Analysis and Applications, 388(1):248–271, 2012.
  * [20] Chao Zhu and George Yin.  Asymptotic properties of hybrid diffusion systems.  SIAM Journal on Control and Optimization, 46(4):1155–1179, 2007\.