A model of hepatitis B virus with random interference infection rate

1.   Introduction

In real life, in the face of emerging and re-emerging infectious diseases, human behavior and environmental noise interference will produce many uncertain factors, which may accelerate the spread of the disease or make the disease exist for a long time, or even make the local epidemic situation beyond control. For example, currently popular novel coronavirus pneumonia (NCP), long-term existence of AIDS (AIDS), hepatitis B virus (HBV) and seasonal influenza virus (influenza virus). For these diseases, although there is no specific drug cure, standard treatment and appropriate prevention and control measures can control the trend of viral infection. In the long term treatment, drug tolerance and human activity with uncertain interference will affect the control of the transmission of the virus.

Many scholars find that the resistance of hepatitis B virus is affected by many factors when they study the mechanism of hepatitis B resistance. Zhang Jing et al. studied the changes of drug-resistant mutants in 268 cases of chronic hepatitis B patients and found that with the extension of drug treatment time, the increase probability of drug-resistant strains in different patients was different, but all of them would increase, even up to 100% ^[1]. Deng Jun et al. studied the mechanism of drug resistance of hepatitis B and found that multiple factors could lead to drug resistance, such as the reliability of viral polymerase, the pressure of drug selection, the patient's past drug use history and the patient's genetic factors (such as congenital metabolic defects), and the sensitive virus strains in the patient could mutate at one or more sites ^[2,3].

In order to be more practical, many scholars have combined the influence of environmental fluctuation factors to study the development trend of infectious diseases with stochastic models. Pang Haiyan et al. considered that variables oscillated around the equilibrium point caused by white noise, and the disturbance value was proportional to the amplitude of oscillation, thus establishing the virus dynamics model under immune response damage, and obtained the conclusion that the positive equilibrium point of the model was stable under random disturbance ^[4,5]. Xie Falan et al. took into account that the regeneration ability of liver cells was interfered by environmental factors, so that the infection rate would fluctuate around a certain average value. They established a random hepatitis B virus infection model with Logistic liver cell growth, and proved that the solution of this model had a unique stable distribution ^[6,7]. Hui Hongwen et al. established an HBV model with nonlinear incidence and random interference, obtained sufficient conditions for the exponential stability of the free virus in infected cells, estimated the oscillating behavior of the model near its deterministic model, and verified it by numerical simulation ^[8]. Bao Kangbo et al. analyzed the influence of environmental noise on stochastic hepatitis B virus dynamics model, and observed that high environmental noise intensity can inhibit the outbreak of hepatitis B, indicating the important role of intervention strategies in the control of hepatitis B ^[9,10]. Taking into account the variables directly proportional to the ambient noise, Liu et al. constructed a three-dimensional stochastic virus infection model, and obtained sufficient conditions for the ergonomically stationary distribution of the model solution and conditions for the extinction of the disease ^[11,12,13].

Xia Peiyan established a random HTLV-I model with CTL immune response by considering the random perturbation of exposure rate, and proved the existence of a unique global positive solution to the model, and found the random threshold values of virus extinction, latent and epidemic ^[14,15,16]. By considering the random interference of the environment on the exposure rate, Cui Xiaowei established the dynamic model of virus infection with random immune effect, proved that the unique global positive solution existed in the model solution, obtained the threshold conditions of virus extinction and persistence, and gave the upper and lower bounds of the persistence, and verified the theory through numerical simulation ^[17]. The research with stochastic model changes the law of the original determined model ^[18,19,20], and it is a meaningful work to pay attention to the problems of environment and personal noise. For patients with hepatitis B virus infection, white noise will interfere with the original infection rate in the presence of external environmental interference or drug resistance during long-term treatment. Therefore, the dynamic model of hepatitis B virus with random infection rate can be used to predict the effect of virus control under treatment more accurately. This is worthwhile work.

In the second part of this paper, according to the statistical data rules of drug treatment of hepatitis B virus, it is found that drug-resistant HBV is uncertain due to viral mutation, which resulted in the randomness of the infection rate of healthy cells, and a stochastic hepatitis B virus dynamics model is established. In the third part, Lyapunov function and Ito's formula are used to prove the existence and uniqueness of the global positive solution of the model and the asymptotic property of the model solution. In the fourth part, several cases of the influence of noise intensity on HBV treatment are simulated numerically, and the rationality of the theorem is also verified. By comparing the results of no noise interference model, it is found that strong noise can affect the variation law of hepatitis B virus.

2.   Model establishment

According to the principle of action of anti-hepatitis B virus drugs and the rule of therapeutic test data, it was found that the drug resistance of the virus was the main reason for the decrease of drug efficacy ^[21]. Therefore, the virus can be divided into drug sensitive virus and drug resistant virus. It is assumed that the inhibition rate of drugs to the resistant virus is related to the mutation rate of the virus, and the following kinetic model of hepatitis B virus with drug resistance is established.

here $T$, ${T_s}$, ${V_s}$, ${T_r}$ and ${V_r}$ respectively represent the number of uninfected hepatocytes at time $t$, the number of hepatocytes infected with drug-sensitive HBV, the number of drug-sensitive HBV, the number of hepatocytes infected with drug-resistant HBV, and the number of drug-resistant HBV; $\lambda$ is the growth rate of uninfected hepatocytes; $\alpha$ and $\delta$ are the death rates of hepatocytes and viruses; ${K_s}$ and ${K_r}$ respectively denote the infection rate of drug-sensitive HBV and drug-resistant HBV on uninfected hepatocytes; ${N_s}$ and ${N_r}$ respectively denote the total number of viruses produced by drug-sensitive infected hepatocytes and drug-resistant infected hepatocytes during their life cycle; the mutation rate between drug-sensitive hepatocytes and drug-resistant hepatocytes is given by $u$; $f({C_N})$ is the rate at which drug therapy causes the decrease of HBV; ${C_N}$ denotes the average steady-state plasma concentration in a patient; $\beta (u)$ denotes the inhibition rate of drug therapy on drug-resistant HBV; all parameters are positive. Suppose that $f({C_N})$ is a bounded function, and $\beta = 1 - u$. The basic regeneration number is as follows

The model has four equilibrium points, namely, the disease-free equilibrium point ${E_0}$, the sensitive virus equilibrium point ${E_s}$, the drug-resistant virus equilibrium point ${E_r}$ and the double virus equilibrium point ${E^{\text{*}}}$. When the basic regeneration number meets different conditions, the equilibrium point is globally asymptotically stable.

In drug therapy, because the mechanism of its action on the virus is different from the body's tolerance, the blood drug concentration is inversely proportional to the rate of decline of HBV DNA, which can be used to determine the occurrence of HBV mutation. Taking lamivudine as an example, Su Ying et al. studied its drug-resistant mutation site in the catalytic region of HBV polymerase or YMDD sequence in the C region, which was due to the mutation of rtM204V /I (methionine M at position 204 was replaced by isoleucine I or valine V), and the drug-resistant strains formed by the mutation of this sequence were called YMDD variants ^[2]. It was also found that there were some variations in the B region. The susceptibility of the virus to the drug is slightly reduced when only the B region variation is present, and is greatly reduced when they are accompanied by the RTM204V/I variation. Due to the differences in the variation types of HBV in patients, some patients only have YMDD variant strains in the body, some patients only have mutations in the B region, and some have these two variants. So the more mutated sites on HBV DNA, the more likely the virus is to infect liver cells, which leads to different rates of infection of the mutant strains on healthy liver cells.

According to the statistical analysis of the duration of medication and the number of virus strains ^[1], the length of treatment would also lead to different probabilities of YMDD variant strains appearing in patients, as shown in Figure 1. Even for the same treatment duration, the percentage of patients with the YMDD variant varies, as shown in Figure 2. Due to the influence of individual genetic factors or variation characteristics in B region in vivo, the replication rate of HBV in patients is different, and the infection rate of the virus on liver cells is uncertain.

In ^[22], the hepatitis B virus genome is double-stranded relaxed circular DNA (rcDNA). When HBV DNA enters the cell, the rcDNA becomes covalently closed circular DNA (cccDNA), which is used as a template to start the replication process of the virus. Therefore, cccDNA is stable in infected liver cells, and there are about 5-50 cccDNA copies in each liver cell. When liver cells are damaged, cccDNA in the liver nucleus is released into the blood, so the amount of cccDNA detected in liver tissue and serum can be used to predict the number of infected cells. Therefore, the total amount of cccDNA and HBV DNA can be measured to determine the relationship between infected liver cells and the virus. In ^[23], changes in liver tissue cccDNA, serum cccDNA, and HBV DNA in 88 patients with chronic hepatitis B were recorded before and at 24, 48, and 96 weeks after treatment. The data are shown in Table 1.

The total amount of cccDNA can be considered to be composed of cccDNA in liver tissue and cccDNA in serum, represented by variable $x$; The amount of serum HBV DNA was denoted by variable $y$. Statistical analysis shows that the two are linearly correlated and satisfy the relationship $y = 0.6302x - 0.6822$, and the degree of correlation is $r = 0.9906$. The drawing graph is shown in Figure 3.

According to the image and fitting results, it can be approximately believed that there is a linear proportional relationship between infected liver cells and virus, as shown in the following expression

To control the infection of healthy liver cells, that means to reduce the infection rate. According to the above data analysis, it is found that multiple factors during treatment will interfere with the infection rate, which may lead to some changes in the results given by Model (1). If the influence of random factors on the infection rate is considered, the infection rate ${K_r}$ is rewritten into the following form

Let $(\Omega , \Phi , P)$ be a complete probability space and $\Phi$ is the $\sigma$ -algebra. In this probability space, a $\sigma$ -algebraic stream ${\left\{ {{F_t}} \right\}_{t \geqslant 0}}$ is defined, which satisfies the usual conditions.

1) For all $0 \leqslant s < t < \infty$, ${\Phi }_{s}\subset {\Phi }_{t}\subset \Phi$,

2) Right continuity: for all $t \geqslant 0$, ${\Phi }_{t} = {\cap }_{s > t}{\Phi }_{s}$.

Moreover, the $B(t)$ belongs to the independent Brownian movement in this probability space, and $B(0) = 0$. ${\sigma ^2} > 0$ is noise intensity.

In summary, the following hepatitis B virus model with random interference infection rate is established based on Model (1):

The variables and other parameters in Model (3) are the same as those in Model (1).

3.   Main research results

3.1.   Existence and uniqueness of global positive solutions

Theorem 1. For any given initial value $(T(0), {T_s}(0), {V_s}(0), {T_r}(0), {V_r}(0)) \in R_ + ^5$ , Model (3) has a unique global positive solution $(T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t))(t \geqslant 0)$ . That is, $(T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t)) \in R_ + ^5$ a.s..

Proof. The function at the right end of Model (3) satisfies the local Lipschitz condition, so any given initial value $(T(0), {T_s}(0), {V_s}(0), {T_r}(0), {V_r}(0)) \in R_ + ^5$, there is a unique local solution $y(t) = (T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t))$ on $[0, {\tau _e})$, where ${\tau _e}$ is the blasting time, the solution is proved to be global, only need to prove ${\tau _e} = \infty$ a.s..

Take a sufficiently large positive number ${k_0} \geqslant 0$, so that every component of the initial $y(0)$ is a member of the interval $[\frac{1}{{{k_0}}}, {k_0}]$. For any $k \geqslant {k_0}$, $k \in {N_ + }$, define the stopping time

here $\inf \emptyset = \infty$ ($\emptyset$ stands for the empty set).

From the definition of the stop time, we know that when $k \to \infty$, ${\tau _k}$ monotonically increases. Let ${\tau _\infty } = \mathop {\lim }\limits_{k \to \infty } {\tau _k}$, then ${\tau _\infty } \leqslant {\tau _e}$ a. s..

If ${\tau _\infty } = \infty$ can be proved for all $t \geqslant 0$ and $(T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t)) \in R_ + ^5$, then ${\tau _e} = \infty$. Use contradiction. If not, there are constants $\tilde t > 0$ and $\delta \in (0, 1)$, such that $P{\text{\{ }}{\tau _k} \leqslant \tilde t{\text{\} }} > \delta$, so there is an integer ${k_1} \geqslant {k_0}$, and for any $k \geqslant {k_1}$, there is

when $t \leqslant {\tau _k}$, for each $k$

Thus, it is obtained from Eq (5)

where

Define the function $V$ : $R_ + ^5 \to {R_ + }$

where $a$ is the positive number to be determined.

It is easy to prove that inequality $u + 1 - \ln u > 0(u > 0)$ is true, so there is $V(T, {T_s}, {V_s}, {T_r}, {V_r}) > 0$. From ${\rm{It\hat o}}$ formula, we get

where

Let

Based on Eq (6), Eq (7) can be reduced to

where

Therefore

Integrate Eq (8) from 0 to ${\tau _k} \wedge \tilde t$ and take the expectation

Let ${\Omega _k} = \{ {\tau _k} \leqslant \tilde t\}$, then for any $k \geqslant {k_1}$, as shown in Eq (4), there is $P({\Omega _k}) \geqslant \delta$, and for any $\omega \in {\Omega _k}$, as can be seen from the definition of stop time, at least one of them is equal to $k$ or $\frac{1}{k}$, so there is

It can be obtained from Eq (9) and the above equation

where ${I_{{\Omega _k}}}$ represents an indicator function of ${\Omega _k}$.

When $k \to \infty$, the above equation is in contradiction with the boundedness of function $V$, so there must be ${\tau _\infty } = \infty$ a.s., that is, the solution is positive and globally unique.

The result of Theorem 1 and Eq (3) can be obtained

For any given initial value $(T(0), {T_s}(0), {V_s}(0), {T_r}(0), {V_r}(0))$, the solution of Model (1) is ultimately bounded, thus the region

is the positive invariant set of Model (3). Let's assume $(T(0), {T_s}(0), {V_s}(0), {T_r}(0), {V_r}(0)) \in {\Gamma ^{\text{*}}}$, so all the questions that follow will be discussed in ${\Gamma ^{\text{*}}}$.

3.2.   The asymptotic behavior of the solution of Model (3)

Theorem 2. If ${R_s} < 1 - u$ and ${R_r} < 1$ , the solution $(T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t))$ of Model (3) has

where $m = \min \{ \alpha , \frac{1}{{2N_s^2}}(\delta + f({C_N})), \frac{1}{{2N_r^2}}(\delta + f({C_N})), \frac{{\frac{\alpha }{{{N_s}(\delta + f({C_N}))}}}}{{2(1 - u)\frac{\lambda }{c}{K_s}}}, \frac{1}{2}{C_{\text{2}}}\alpha {\kern 1pt} {\kern 1pt} \} ,$

Proof. Define the ${C^2}$ -function ${V_i}(i = 1, 2, 3, 4, 5)$

Applying ${\rm{It\hat o}}$ formula to ${V_i}(i = 1, 2, 3, 4, 5)$, we get

From ${R_s} < 1 - u$ and ${R_r} < 1$, Eq (10) can be reduced to the following form

We can define the function $V$ From ${V_i}(i = 1, 2, 3, 4, 5)$

where ${C_1}, {C_2}$ is the normal number to be determined.

Substitute Eqs (11)–(15) into Eq (10) to obtain

where

Inequality $2ab \leqslant {a^2} + {b^2}$ and the expressions in the second, third and fourth lines of Eq (17) can be obtained

Choose

meet

Therefore, when Eqs (19) and (20) are substituted into Eq (17), there is

Then substitute Eq (21) into Eq (16), and take the expectation of both ends after the integration from 0 to $t$, we get

Therefore

where $m = \min \{ \alpha , \frac{1}{{2N_s^2}}(\delta + f({C_N})), \frac{1}{{2N_r^2}}(\delta + f({C_N})), \frac{{\frac{\alpha }{{{N_s}(\delta + f({C_N}))}}}}{{2(1 - u)\frac{\lambda }{c}{K_s}}}, \frac{1}{2}{C_{\text{2}}}\alpha {\kern 1pt} {\kern 1pt} \}$

According to Theorem 2, the solution of Model (3) will oscillate around the equilibrium point of Model (1). The amplitude of oscillation is positively correlated with the noise intensity.

In the special case, when $\sigma = 0$, there is

This indicates that its equilibrium point ${E_0}$ with respect to Model (1) is globally asymptotically stable.

Theorem 3. If ${R_s} > {R_r}$ and ${R_s} > 1$ , and satisfy:

where

where

where

Then, the solution $(T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t))$ of Model (3) has

where

Proof. Define the ${C^2}$ -function ${V_i}(i = 1, 2, 3, 4, 5)$

where ${b_1},$ ${b_2}$ is the normal number to be determined.

Applying ${\rm{It\hat o}}$ formula to ${V_i}(i = 1, 2, 3, 4, 5)$, we get

From the inequality $2ab \leqslant {a^2} + {b^2}$, $L{V_i}(i = 1, 2, 3, 4, 5, 6)$ in Eq (21) can be reduced to the following form

We can define the function $V$ From ${V_i}(i = 1, 2, 3, 4, 5)$

Substituting Eqs (23)–(28) into Eq (22), we get

Then, $LV$ in Eq (26) can be obtained from inequality $2ab \leqslant {a^2} + {b^2}$ and condition (H)

where

There are positive numbers ${b_1}$ and ${b_2}$ such that ${m_2}({b_1}) > 0$ and ${m_3}({b_1}, {b_2}) > 0$ are true.

In fact, it's just a matter of choosing ${b_1}$ and ${b_2}$ so that $f({b_1}) > 0$ and $g({b_1}, {b_2}) > 0$ are true. Because $f({b_1})$ is a downward opening quadratic function of ${b_1}$, the discriminant is

We know from (H₁) that $f({b_1}) = 0$ has two positive roots, and the expression is as follows

So, for any ${b_1} \in ({b_{11}}, {b_{12}})$, there is a $f({b_1}) > 0$. Therefore, ${m_2}({b_1}) > 0$.

We know from condition (H₂) that we choose $b_1^* \in ({b_{11}}, {b_{12}})$ such that $A(b_1^*) > 0$, $B(b_1^*) > 0$.

From ${m_1}(b_1^*, {b_2}) > 0$, we can get

Substituting $b_1^*$ into $g({b_1}, {b_2})$ we get that $g(b_1^*, {b_2})$ is a downward opening quadratic function of ${b_2}$, the discriminant is

We know from (H₃) that $g(b_1^*, {b_2}) = 0$ has two positive roots, and the expression is as follows

So, for any ${b_2} \in ({b_{21}}, {b_{22}})$, there is a $g({b_2}) > 0$. Therefore, ${m_3}(b_1^*, {b_2}) > 0$.

We know from (H₄) that there is $b_2^{\text{*}} \in (\max \{ {b_{21}}, \frac{{3K_r^2}}{{2A(b_1^*)(\delta + \beta (u)f({C_N}))}}\} , {b_{22}})$, such that ${m_1}(b_1^*, b_2^*) > 0$ and ${m_3}(b_1^*, b_2^{\text{*}}) > 0$ are true.

Equation (30) is simplified as

Integrate Eq (29) from 0 to $t$ and take the expectation of both ends, and we get

Therefore

where

According to Theorem 3, the solution of Model (3) will oscillate around the equilibrium point of Model (1).The amplitude of oscillation is positively correlated with the noise intensity.

In the special case, when $\sigma = 0$, there is

This indicates that its equilibrium point ${E^{\text{*}}}$ with respect to Model (1) is globally asymptotically stable.

Theorem 4. If ${R_s} > 1$ and ${R_r} < 1$ , and satisfy:

where

where

where $A(b_1^{\text{*}}) = \alpha + \frac{{c(\alpha + {K_s}{V_{s1}})}}{\lambda } - \frac{{5K_s^2}}{{2b_1^{\text{*}}(\delta + f({C_N}))}},$

Then, the solution $(T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t))$ of Model (3) has

where

Proof. Define the ${C^2}$ -function ${V_i}(i = 1, 2, 3, 4, 5)$

where ${b_1}$, ${b_2}$ is the normal number to be determined.

Applying ${\rm{It\hat o}}$ formula to ${V_i}(i = 1, 2, 3, 4, 5)$, we get

Same as the proof method of Eqs (23)–(28) in Theorem 3, $L{V_i}(i = 1, 2, 3, 4, 5, 6)$ in Eq (32) is simplified as

We can define the function $V$ From ${V_i}(i = 1, 2, 3, 4, 5)$

Substituting Equations (33)–(38) into Eq (32), we get

Then, $LV$ in Eq (39) can be obtained from inequality $2ab \leqslant {a^2} + {b^2}$ and condition (H)

where

Same as the proof of Theorem 3, under the conditions (H₁)–(H₄), we can find positive values $b_1^{\text{*}}$ and $b_2^{\text{*}}$ such that ${m_2}(b_1^{\text{*}}) > 0$ and ${m_3}(b_1^*, b_2^*) > 0$. Here, the discriminants are respectively

$b_1^{\text{*}}$ and $b_2^{\text{*}}$ satisfy the following conditions

The selected parameters $b_1^*$ and $b_2^*$ are substituted into Eq (40), and $LV$ is simplified as

Integrate Eq (39) from 0 to $t$ and take the expectation of both ends, and we get

Therefore

where $m = \min \{ {m_1}(b_1^*, b_2^*), {m_2}(b_2^*), {m_3}(b_1^*, b_2^*), \frac{{b_1^*(\delta + f({C_N}))}}{2}, \frac{{b_2^*(\delta + f({C_N}))}}{2}{\kern 1pt} {\kern 1pt} \} ,$

According to Theorem 4, the solution of Model (3) will oscillate around the equilibrium point ${E^*}$ of Model (1). The amplitude of oscillation is positively correlated with the value of noise intensity.

In the special case, when $\sigma = 0$, there is

This indicates that its equilibrium point ${E_s}$ with respect to Model (1) is globally asymptotically stable.

Theorem 5. If ${R_s} < 1$ and ${R_r} > 1$ , and satisfy:

where

where

where $A(b_1^{\text{*}}) = \alpha + \frac{{c(\alpha + {K_r}{V_{r2}})}}{\lambda } - \frac{{5K_s^2}}{{2b_1^{\text{*}}(\delta + f({C_N}))}}.$

Then, the solution $(T(t), {T_s}(t), {V_s}(t), {T_r}(t), {V_r}(t))$ of Model (3) has

where $m = \min \{ {m_1}(b_1^*, b_2^*), {m_2}(b_2^*), {m_3}(b_1^*, b_2^*), \frac{{b_1^*(\delta + f({C_N}))}}{2}, \frac{{b_2^*(\delta + f({C_N}))}}{2}{\kern 1pt} {\kern 1pt} \} ,$

Proof. Define the ${C^2}$ -function ${V_i}(i = 1, 2, 3, 4, 5)$

where ${b_1},$ ${b_2}$ is the normal number to be determined.

Applying ${\rm{It\hat o}}$ formula to ${V_i}(i = 1, 2, 3, 4, 5)$, we get

Same as the proof method of Equations (23)–(28) in Theorem 3, $L{V_i}(i = 1, 2, 3, 4, 5, 6)$ in Eq (42) is simplified as

We can define the function $V$ From ${V_i}(i = 1, 2, 3, 4, 5)$

Substituting Eqs (43)–(48) into Eq (42), we get

Then, $LV$ in Eq (46) can be obtained from inequality $2ab \leqslant {a^2} + {b^2}$

where

Same as the proof of Theorem 3, under the conditions (H₁)–(H₄), we can find positive values $b_1^{\text{*}}$ and $b_2^{\text{*}}$ such that ${m_2}(b_1^{\text{*}}) > 0$ and ${m_3}(b_1^{\text{*}}, b_2^{\text{*}}) > 0$. Here, the discriminants are respectively

$b_1^{\text{*}}$ and $b_2^{\text{*}}$ satisfy the following conditions

Therefore, defined by Eq (41), the selected parameters $b_1^{\text{*}}$ and $b_2^{\text{*}}$ are substituted into Eq (50) to be simplified as

Integrate Eq (49) from 0 to $t$ and take the expectation of both ends, and we get

Therefore

where

According to Theorem 5, the solution of Model (3) will oscillate around the equilibrium point ${E_r}$ of Model (1). The amplitude of oscillation is positively correlated with the noise intensity $\sigma$.

In the special case, when $\sigma = 0$, there is

This indicates that its equilibrium point ${E_r}$ with respect to Model (1) is globally asymptotically stable.

4.   Numerical simulation

The resistance of hepatitis B virus due to individual differences or the interference of factors such as environmental changes may cause the virus to mutate and change the original infection rate. The drug cannot accurately act on the mutant virus, and the infection rate is uncertain, so it will be affected by the clinic. Treatment brings many unknown results. By studying the asymptotic behavior of the mathematical model solution with random noise, understanding the influence of the interference of uncertain factors on the model solution can provide a theoretical basis for determining the treatment plan. The following is numerical simulation of the influence of random noise intensity on the model solution.

4.1.   Parameter selection

Some of the parameters in Model (1), such as the growth rate $\lambda$ of healthy liver cells, the mortality rate $\alpha$ of liver cells, the mortality rate $\delta$ of viruses, the number of viruses ${N_s}$ and ${N_r}$ produced in the life cycle of sensitive and drug-resistant infected liver cells, are shown in Table 2 ^[24], and the average steady-state plasma concentration is ${C_N} = 0.8$ ng/ml.

Choose four situations to simulate the difference caused by noise. Related parameters such as virus mutation rate $u$, drug inhibition rate $\beta$ of virus, sensitive HBV infection rate ${K_s}$, drug-resistant HBV infection rate ${K_r}$ and noise intensity $\sigma$ are shown in Table 3.

4.2.   Influence of noise intensity on the solution of Model (1)

In order to better understand the interference of noise on virus changes, Models (1) and (3) are simulated and compared. Divide Models (1) and (3) into sensitive sub-models and drug-resistant sub-models. The initial values of the model are taken as follows

Case 1: Substituting the parameters in Tables 2 and 3 into Eq (2) can calculate ${R_s} = 0.5778$ and ${R_r} = 0.1825$. It can be seen from paper ^[14] that Model (1) has a disease-free equilibrium point ${E_0}$, which is a global asymptotically stable point. From Theorem 2, we can see that the solution trajectory of Model (1) is asymptotically oscillating around ${E_0}$. The simulation results are shown in Figures 3 and 4.

Comparing Figures 3 and 4, it can be found that when the noise intensity increases, the vibration amplitude of the solution of the Model (3) around ${E_0}$ will also increase.

Case 2: Substituting the parameters in Tables 2 and 3 into Eq (2) can calculate ${R_s} = 1.1920$ and ${R_r} = 1.1853$. It can be seen from paper ^[14] that Model (1) has a drug-sensitive type virus equilibrium point ${E_s}$, which is a global asymptotically stable point. From Theorem 4, we can see that the solution trajectory of Model (1) is asymptotically oscillating around ${E_s}$. The simulation results are shown in Figures 5 and 6.

Comparing Figures 5 and 6, it can be found that when the noise intensity increases, the vibration amplitude of the solution of the Model (3) around ${E_s}$ will also increase.

Case 3: Substituting the parameters in Tables 2 and 3 into Eq (2) can calculate ${R_s} = 2.1286$ and ${R_r} = 1.2733$. It can be seen from paper ^[14] that Model (1) has double virus equilibrium point ${E^*}$, which is a global asymptotically stable point. From Theorem 3, we can see that the solution trajectory of Model (1) is asymptotically oscillating around ${E^*}$. The simulation results are shown in Figures 7 and 8.

Comparing Figures 7 and 8, it can be found that when the noise intensity increases, the vibration amplitude of the solution of the Model (3) around ${E^*}$ will also increase.

Case 4: Substituting the parameters in Tables 2 and 3 into Eq (2) can calculate ${R_s} = 0.0578$ and ${R_r} = 1.7230$. It can be seen from paper ^[14] that Model (1) has a drug-resistance type virus equilibrium point ${E_r}$, which is a global asymptotically stable point. From Theorem 5, we can see that the solution trajectory of Model (1) is asymptotically oscillating around ${E_r}$. The simulation results are shown in Figures 9 and 10.

Numerical simulation results show that:

Case 4 is different from Cases 1–3. Comparing Figures 9 and 10, when the noise intensity is $\sigma = 0.00000008$, the solution of Model (3) will eventually oscillate around the drug-resistance equilibrium point ${E_r}$; When the noise intensity increases to $\sigma = 0.0000002$, it can be seen from H2) and H4 in Figure 10 that, with the increase of noise intensity, the number of drug-resistant infected cells and drug-resistant HBV will gradually decrease and eventually reach a disease-free equilibrium point. In the other three cases, the oscillations only increased. Therefore, it can be inferred that the increase of noise intensity, which affects the infection rate of drug-resistant HBV, can inhibit the increase of drug-resistant infected cells and drug-resistant HBV when the degree of drug resistance of the virus in patients is high. All Cases 1–4 indicate that with the increase of noise intensity, the oscillation amplitude of the solution of Model (3) will increase, which verifies the conclusion of Theorems 2–5. Case 4 also shows that when drug resistance is high, the increase of noise intensity will inhibit the replication of the virus, that is to say, the interference of human factors can be increased in real life, so that the number of viruses in the patient's body is decrease.

4.3.   Influence of noise intensity on antiviral therapy

If drug-resistant virus has been produced in the treatment process, the interference of noise will produce uncertainty phenomenon, so the change trend of virus cannot be accurately judged, and the subsequent treatment cannot achieve the expected goal. Taking the drug-resistance equilibrium point as an example, the initial value in 4.2 and the parameter value in Table 2 are still selected. The average steady state plasma concentration ${C_N}$ is still 0.8 ng/ml. The infection rate ${K_s}$ of sensitive HBV and the infection rate ${K_r}$ drug-resistant HBV is $9 \times {10^{ - 8}}$ and $5.2 \times {10^{ - 8}}$, the mutation rate $u$ is 0.4, and the drug inhibition rate $\beta$ is 0.6. At this time, the Model (1) has a globally asymptotically stable equilibrium point ${E_r}$. The changing law of the solution of Model (3) is numerically simulated, and the simulation result is shown in Figure 11.

It can be seen from Figure 11 that due to the influence of noise factors, the number of viruses may fluctuate around the lower limit of detection (${10^3}$ copies/mL, the critical value for clinical treatment). For example, in Figure11 i1), during the period of 1070 to 1324 days of antiviral treatment, the number of virus will be below the lower limit of detection many times. From i1) to i3) in Figure 11, it can be seen that as the noise intensity increases, the fluctuations near the lower detection limit (${10^3}$ copies/mL) become more frequent. When it is clinically detected that the HBV DNA of a patient is below the lower limit of detection for many times, the drug will be stopped, but the actual patient is still at the equilibrium point of drug resistance. Stopping the drug will lead to a rapid increase in the number of viruses. The intensity of noise provides an illusion for clinical detection of HBV DNA.

The simulation results found that when the drug-resistant HBV in the patient is the main group of viruses, the increase in noise intensity that affects the infection rate of drug-resistant HBV can effectively inhibit virus replication and reduce drug-resistant infected cells and HBV. But it will also interfere with the correct detection of HBV DNA. Therefore, it is recommended that in the process of clinical treatment, when the patient has strong drug resistance, measures should be taken to suppress the infection rate of the mutant strain, such as dressing change or combined drug therapy, so as to achieve the purpose of treating resistant infection cells and eliminating HBV. When the noise intensity is large, the number of detections should be strengthened during the treatment process, and the time of two detections should be lengthened to improve the accuracy of detection. Even when the number of viruses in the patient's body reaches the lower limit of detection, reasonable consolidation treatment should be carried out to make the patient's treatment more thorough. Avoid noise interference to make the detection inaccurate and affect the patient's antiviral treatment.

5.   Conclusions

This paper established a hepatitis B virus dynamic model with random interference infection rate, proved the existence and uniqueness of the global positive solution, and obtained the asymptotic behavior of the model solution near the equilibrium point of the deterministic model. Through numerical simulation, the conclusion of the theorem is verified, and when the random interference is small enough, the solution of the random model vibrates around the equilibrium point of the determined model. It is found that the amplitude of the oscillation is proportional to the noise intensity.

Acknowledgments

This work was supported by the project of Nature Scientific Foundation of Heilongjiang Province (A2016004).

Conflict of interest

The authors declare there is no conflict of interest.