Deterministic and stochastic model for the hepatitis C with different types of virus genome

1.   Introduction

Hepatitis C is a disease caused by a virus that infects the liver. The virus, called hepatitis C virus or HCV, is just one of the hepatitis viruses. Global statistics estimate that nearly 350 million people worldwide are infected with HCV, which significantly affects human health, especially the liver ^[16]. Mathematical models are currently a useful way to assess the transmissibility of many infectious diseases, predict future morbidity, and assess the efficacy of prevention and therapy ^{[1,11,14,17,21,34,35,36,42,46]}. Many successful mathematical models investigating fundamental hepatitis C viral dynamics have been produced in recent decades, including Feng et al.^[12], Lestari et al. ^[23] and Li et al. ^[24] to cite only a few. Moneim and Mosa ^[31,32], constructed a mathematical model to study the spread of HCV-subtype 4a amongst the Egyptian population. They divided the population into three groups, susceptible class, $S(t)$, the hepatitis C subtype 4a infective class, $I_1$ and the HCV from all other subtype classes, $I_2$. They assumed that all types of the class $I_2$ could mutate to $I_1$ at a constant rate $\mu > 0$. A few years after the appearance of ^[31,32], effective treatment for the hepatitis C virus appeared. For this, we will add to the mathematical model another class $(R(t))$, which represents those who have been cured of this virus. Following ^{[8,18,19,39,47]}, the saturated incidence rate will replace the bilinear incidence rate. The hepatitis C virus becomes as follows.

where $A$ denotes the recruitment rate of susceptible individuals, $a$ expresses the half-saturation constant for susceptible individuals with HCV. $\gamma$ and $\delta$ are the recovery rates, $k_1$ is the transmission rate of virus C when susceptible individuals $S(t)$ contact with corresponding infected $I_1(t)$ individuals, and $k_2$ is the transmission rate when $S(t)$ contact with corresponding $I_2(t)$ individuals. Following ^[31,32], we assume that the population is mixing in a homogenous manner, i.e., every person has the same chance of coming in contact with an infected person also we assume that the birth and death rates are equal and positive constant rate $b$.

The deterministic HCV system (1.1) ignores the possible importance of environmental noise. In reality, stochastic effects can be important during the transmission of hepatitis C viruses, because various cells and infective virus particles react differently in the same environment. Consequently, the deterministic models do not provide an adequate understanding of viral dynamics due to the stochastic behavior of viruses and the complexity of the immune system. The primary purpose of this paper is to propose and analyze a deterministic and stochastic model of HCV. The paper is arranged as follows: In Section 2, the dynamics of the deterministic system are verified. The stochastic extension of the deterministic model is performed in Section 3, and the existence of a unique global positive solution for the stochastic model is investigated, and the sufficient conditions for the extinction of the hepatitis C virus from the stochastic system are obtained. In Section 4, some numerical simulations are presented to verify the obtained theoretical results. Finally, Section 5 contains the conclusion.

2.   Dynamics of deterministic system

In the following, two critical parameters $R_{01}$, and $R_{02}$, can be used to classify the dynamics of the hepatitis C virus model (1.1). The threshold parameter $R_{01}$ defined by $R_{01} = \frac{ Ak_1}{\rho(A+ a b)},$ $\rho = b+\gamma$. While the threshold parameter $R_{02}$ defined by $R_{02} = \frac{ A k_2}{\theta(A+ a b)}$, $\theta = b+\delta$. Using the next generation method, one can obtain the basic reproduction number $R_{0}$ defined by $R_{0} = {\rm{max}} \left\lbrace R_{01}, R_{02} \right\rbrace$.

The hepatitis C virus model (1.1) has the following four equilibrium points:

(1) The disease-free equilibrium point $E_0 = \left(\frac{A}{b}, 0, 0, 0 \right)$, which always exists. At this point, all individuals are susceptible, and there is no infection present in the population.

(2) The equilibrium point $E_{1} = (S_1, I_{11}, 0, R_1)$, where

$E_1$ exists if $k_1 > \rho$ and $R_{01} > 1$.

(3) The equilibrium point $E_{2} = (S_2, 0, I_{22}, R_2)$, where

$E_2$ exists if $k_2 > \theta$ and $R_{02} > 1$.

(4) The coexistence equilibrium point $E_3 = (S_3, I_{13}, I_{23}, R_3)$, where

$\beta = a b-A+\frac{\rho k_2}{\mu}-\frac{\theta k_1}{\mu}$. The coexistence equilibrium point $E_3$ exists if $\frac{k_2 S_3}{\theta(a+S_3)} > 1$ and $\frac{k_1 S_3}{\rho(a+S_3)} < 1.$

The boundedness of the solutions of model (1.1) is given as follows. Let $(S(t), I_1(t), I_2(t), R(t))$ be any solution of system (1.1) with non-negative initial conditions and assume that the total population size $N(t) = S(t)+I_1(t)+I_2(t)+R(t),$ then $\frac{d N}{dt}+b N = A,$ consequentially $N(t) = \frac{A}{b}+N(0) e^{-b t}$, thus it follows that $0\le N(t)\le\frac{A}{b}, \; \rm{as}\; t\to\infty.$

The locally asymptotically stable equilibrium points of system (1.1) are now investigated. The Jacobian matrix is given as follows:

The eigenvalues of $J(E_{0})$ are $-b, \, \, \, -b, \, \, \, \, \rho (R_{01}-1)$, and $\theta (R_{02}-1).$ Thus, $E_0$ is locally asymptotically stable if $R_0 < 1$. The eigenvalues of $J(E_1)$ are $\lambda_1 = -b$ and $\lambda_2 = \frac{k_2 S_1}{a+S_1}-\mu I_{11}-\theta$. The other roots are determined by

The eigenvalues $\lambda_3$ and $\lambda_4$ have negative real parts. Thus, if $\frac{k_2 S_1}{(\mu I_{11}+\theta)(a+S_1)} < 1$ then the equilibrium point $E_1$ is locally asymptotically stable. The eigenvalues of $J(E_2)$ are $\lambda_1 = -b$ and $\lambda_2 = \frac{k_1 S_2}{(a+S_2)}+\mu I_{22}-\rho$. The other roots are determined by

The eigenvalues $\lambda_3$ and $\lambda_4$ have negative real parts. Thus, if $\frac{k_1 S_2}{(a+S_2)}+\mu I_{22} < \rho$ then the equilibrium point $E_2$ is locally asymptotically stable. The stability of the coexistence equilibrium point $E_{3} = (S_3, I_{13}, I_{23}, R_3)$ is investigated as follows

The first eigenvalues of $J(E_3)$ is $\lambda_1 = -b$. The other roots are determined by

where

It is clear that $c_1, c_2, c_3 > 0$ and $c_1 c_2 -c_3 > 0$. It follows from the Routh-Hurwitz criterion that all roots of (2.1) have negative real parts. Thus, the equilibrium point $E_3$ is locally asymptotically stable.

3.   Dynamics of stochastic model

In this section, we will perform a stochastic extension of model (1.1). The hepatitis virus C deterministic model (1.1) will be extended to include the environmental noise as follows.

where $W = \left\lbrace W_1, W_2, W_3, W_4, t\geq0\right\rbrace$ represents the four-dimensional standard Brownian motions with $W_i(0) = 0$ and $\sigma_i^2(i = 1, 2, 3, 4)$ denote the intensities of the white noise. The white noise is defined in a complete probability space $(\Omega, \mathcal F_{t\geq0}, \mathbb P)$ with a filtration $\mathcal F_{t\geq0}$ satisfying the usual conditions. In many applications, the solution of the $\rm{Itô}$ stochastic differential equation must preserve the positivity of the solutions ^[9,10,37]. According to Theorem 2.2 and Corollary 1 in ^[9], the solutions of (3.1) emanating from nonnegative initial data (almost surely) remain nonnegative as long as they exist. The next theorem gives another approach according to ^[30] to prove the existence and uniqueness of a positive global solution of the system (3.1). This approach has recently been used in many papers, and one can highlight ^{[4,5,6,7,13,17,20,22,27,28,33,38,40,41,42,43,44,45,48]}.

Theorem 1. For any given initial value $(S(0), I_1 (0), I_2(0), R(0)) \in \mathbb{R}_+^4$, there exists a unique solution $(S(t), I_1 (t), I_2(t), R(t))$ of system (3.1) for $t\geq 0$ and the global positive solution remains in $\mathbb{R}_+^4$ with probability one.

Proof. Firstly, one can consider the local solution $(S(t), I_1 (t), I_2(t), R(t))$ of system (3.1) for $t\in \left[ 0, \tau_e \right)$, where $\tau_e$ is the explosion time ^[30]. By making the transformation of variables

Using $\rm{Itô}$ formula, one can change system (3.1) as follows

For $X_0 = (x_0, y_0, z_0, w_0) \in \mathbb{R}_+^4$, the coefficients of system (3.2) satisfy the local Lipschitz conditions, consequently, there exists a unique local solution $(S(t), I_1 (t), I_2(t), R(t)) = (e^{x(t)}, e^{y(t)}, e^{z(t)}, e^{w(t)})$ on $\left[ 0, \tau_e \right)$. To ensure that this solution is global, one needs to prove that $\tau_e = \infty$ a.s. Let $s_0 > 0$ be sufficiently large for every coordinate $(S(0), I_1 (0), I_2(0), R(0))$ in the interval $[\frac{1}{s_0}, s_0]$. For each integer $s > s_0$, we define the stopping time

From (3.3), one can note that $\tau_s$ is increasing as $s\to \infty$. Assume $\tau_{\infty} = \lim_{s \to \infty} \tau_s$, then $\tau_\infty\leq \tau_e$ almost sure. In the next, one needs to verify that $\tau_\infty = \infty$. If this is not true, then there exists a constant $T > 0$ and $\epsilon \in (0, 1)$ such that $\mathbb P(\tau_\infty \leq T)\geq \epsilon$. As a result, there exists an integer $s_1\geq s_0$ such that $\mathbb P(\tau_s \leq T)\geq \epsilon, \, \, \, s\geq s_1$. Define the following $C^2$ positive definite function $V_1 (S, I_1, I_2, R)$ as

Using $\rm{Itô}$'s formula, one obtains

Using the inequality $x\leq 2(x+1-{\rm{ln}} x)$, for any $x > 0$, one obtains

which means that

where $K = {\rm{max}} \left\lbrace K_1, K_2\right\rbrace$, $K_1 = A+2b+\rho+\theta+ \frac{1}{2}\sum_{i = 1}^{4} \sigma_i^2$ and $K_2 = {\rm{max}} \left\lbrace 2, 2(\frac{k_1}{a}+\mu+\gamma), 2(\frac{k_2}{a}+\delta)\right\rbrace$. For $t_1 \leq T$, integrating both sides of (3.4) from $0$ to $t_1 \wedge \tau_s$ and then taking the expectation leads to

Following ^{[25,26,29,44]}, applying Grownwall's inequality, one gets

The rest of the proof is similar to ^[25,26,29] and hence is omited here. This complete the proof.

The above theorem shows that the stochastic HCV system (3.1) have positive global solution remain in $\mathbb{R}_+^4$ with probability one. The non-explosion characteristic in an epidemic model is essential but often insufficient. Hence, one needs to indicate that the solution will be ultimately bounded with a large probability. In the following, we will establish the stochastically ultimate boundedness property of the HCV system (3.1).

Lemma 2. Let $N(t) = S(t)+I_1(t)+I_2(t)+R(t)$, then for any given initial value $(S(0), I_1 (0), I_2(0), R(0)) \in \mathbb{R}_+^4$, the following inequality holds:

Proof. It follows from HCV model (3.1) that

Then, the solution of Eq (3.5) has the following form

where

which is a contiuous local martingal with $M(0) = 0$, almost surely. Define $\Lambda(t) = \frac{A}{b}[1-e^{-b t}]$ and $U(t) = N(0)[1-e^{-b t}],$ then we have $N(t) = N(0)+\Lambda(t)-U(t)+M(t).$ Clearly, $\Lambda(t)$ and $U(t)$ are continuous adapted increasing process on $t\geq 0$ with $\Lambda(0) = U(0) = 0$. Then by Theorem 3.9 in ^[30], we have $\lim_{t \to \infty} N(t) < \infty$ a.s.

Remark 3. It follows from the first equation of HCV system (3.1) that

and it follows that

Theorem 4. If $\frac{\sigma_1^2}{2} +\frac{1}{2} < b, \, \, \, k_1+\frac{\sigma_2^2}{2} +\frac{1}{2} < \rho, \, \, \, k_2+\frac{\sigma_3^2}{2} +\frac{1}{2} < \theta, \, \, \, \sigma_4^2 +\frac{1}{2} < b$, then the solutions of (3.1) are stochastically ultimate bounded.

Proof. For $(S(t), I_1 (t), I_2(t), R(t))\in \mathbb{R}_+^4$, define the following function

By $\rm{Itô}$ formula, one has

where

Assume $f(S(t), I_1 (t), I_2(t), R(t)) = (\sigma_1^2 -2b+1)S^2+\left(2 k_1+\sigma_2^2 -2\rho+1\right) I_1^2+$$\left(2k_2+\sigma_3^2 -2\theta+1\right) I_2^2+\left(\sigma_4^2 -2b+1\right) R^2+2\left(A S+\mu I_1^2 I_2+\gamma I_1 R+\delta I_2 R\right).$ According to Lemma 2 and Remark 3, one can find that the function $f(S(t), I_1 (t), I_2(t), R(t))$ has an upper bound.

Let $N_1 = \sup f(S(t), I_1 (t), I_2(t), R(t))+1$. As a result

By $\rm{Itô}$ formula, one obtains

Integrating both sides of the above equation from 0 to t and then taking expectations, one gets

hence

Using Chebyshev's inequality, one can obtains

where $H = \frac{\sqrt{N_1}}{\sqrt{\epsilon}}, \, \, \epsilon > 0.$ Then

This completes the proof.

In the following, we will establish the conditions for extinction of hepatitis C virus from the stochastic system (3.1). The threshold parameter $R_{01}^s$ defined by $R_{01}^s = \frac{ k_1}{\rho+\frac{\sigma_2^2}{2}},$ while the threshold parameter $R_{02}^s$ defined by $R_{02}^s = \frac{ k_2}{\theta+\frac{\sigma_3^2}{2}}$ and the threshold parameter $R_{03}^s$ defined by $R_{03}^s = \frac{k_1+k_2}{\rho+\theta}$. The stochastic basic reproduction number $R_{0}^s$ defined by $R_{0}^s = {\rm{max}} \left\lbrace R_{01}^s, R_{02}^s, R_{03}^s\right\rbrace$.

Theorem 5. The disease die out exponentially with probability one if $R_{0}^s < 1.$

Proof. Let $V_3(I_1, I_2) = {\rm{ln}} (I_1+I_2)$. Using $\rm{Itô}$ formula, one obtains

using integration from $0$ to $t$, one obtains

Applying strong law of large numbers for local martingales one gets

as a result, diseases $I_1(t)$ and $I_2(t)$ die out and tend to zero exponentially a.s., if $R_{0}^s < 1.$

Remark 6. Theorem 3 indicates that the disease die out exponentially a.s., if $R_{0}^s < 1$ with the consequence that the recover class $R(t)$ also goes to extinction a.s.

4.   Numerical simulations

In this part, the numerical results will be compared with the theorems formulated in the previous sections. The interactions between populations classes will be simulated by the following parameters ^[32]: $N = 1000000, \mu = 0.02, a = 1, b = 0.02, k_1 = 0.00001, k_2 = 0.00002, \gamma = 0.001, \delta = 0.001$. To give some numerical finding to the HCV system (3.1), we use the Milstein method mentioned in ^[2,3,15]. The stochastic hepatitis virus C system (3.1) reduces to the following discrete system

where $h$ is a positive time increment and $\epsilon_{ij}, (i = 1, 2, 3, 4)$ are independent random Gaussian variables $N(0, 1)$. One can note that for the given parameters, the value of the stochastic basic reproduction number $R_0^s = 0.00095$. As a result, the conditions of Theorem 5 are verified and the disease die out exponentially with probability one if $R_{0}^s < 1$. Figure 1 represents the dynamical behavior of model (3.1) when the noise strength law ($\sigma_i = 0.01$). For $k_1 = k_2 = 0.0308$, the conditions of Theorem 5 are verified as $R_{0}^s = 1.130769 > 1$ as shown in Figure 2. It is shown that the trajectories of the stochastic system (3.1) oscillates around the coexistence equilibrium of the deterministic system (1.1).

5.   Conclusions

The novelty of this study is that we introduced a new compartment and saturated incidence rate into the classical hepatitis C virus genome model. A deterministic and stochastic model for hepatitis C with different types of virus genomes is proposed and analyzed. We perform a stochastic extension of the deterministic model to study the fluctuation between environmental factors. Firstly, the existence of a unique global positive solution for the stochastic model is investigated. Secondly, sufficient conditions for the extinction of hepatitis C virus from the stochastic system are obtained. Theoretical results are illustrated using numerical simulations. Theoretical and numerical results show that the smaller white noise can ensure the persistence of susceptible and infected populations while the larger white noise can lead to the extinction of populations. For the deterministic model, some sufficient conditions are obtained to ensure the stability of the equilibrium points. The results show that when the basic reproduction number $R_0 < 1$, the deterministic model is asymptotic stable around the free disease equilibrium point $E_0$. By introducing the basic reproduction number $R_0$ and the stochastic basic reproduction number $R_0^s$, the conditions that cause the disease to die out are indicated. Biologically, it can be noted that by social distancing with hepatitis C patients, this can lead to a fluctuation in the value of $k_1$ and $k_2$, and thus $R_0^s$ becomes less than one, which leads to the disappearance of the epidemic.

Acknowledgments

The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research. The authors thank Dr. Sanling Yuan of University of Shanhai for Science and Technology for his kind cooperation.

The researcher would like to thank the Deanship of Scientific Research, Qassim University, for funding the publication of this project.

Conflict of interest

The authors declare that there is no conflict of interests.