Qualitative behaviour of a stochastic hepatitis C epidemic model in cellular level

1.   Introduction

Hepatitis C is a liver inflammatory disease caused by a viral infection. Hepatitis C virus (HCV) is an RNA-type virus from the Flaviviridae family (genus Hepacivirus) which has a high replication rate. Hepatitis can develop into liver cancer ^[1] and can cause serious complications even death ^[2,3]. Around 60 – 85% of acute hepatitis C may develop into a chronic condition, 10–15% may develop cirrhosis, and 25% may develop liver cancer ^[1,2,4]. It is estimated that 58 million people have a chronic hepatitis C, with 1.5 million new infections occurring per year over the world^[5]. There is no vaccine available for HCV ^[3,4,6,7], thus it is important to gain a better insight into the nature of this disease.

One of the accurate tools to understand the dynamics of hepatitis C is using a mathematical model. There are several mathematical models developed for HCV for example in ^[8,9] studying the HCV model in population level. The HCV model incorporating treatment, therapy, or isolation were investigated in ^{[10,11,12,13,14]}. In reference ^[15], the optimal control strategies for HCV epidemics considering the uncertainty of the model was discussed.

The spread of HCV at the cellular level was also investigated by several researchers e.g., ^{[16,17,18,19,20]}. In 1998, Neumann et al. studied the dynamic of HCV and the effect of antiviral. They have found that HCV is highly dynamic. An existing deterministic model cannot be applied to describe randomness in many biological factors, for instance, the random occurrence of cell infection, mutation, and apoptosis ^[21,22]. Furthermore, during the spread of HCV, different cells and infectious virus particles reacting in the same environment can give different effects. In other words, some uncertain factors influence the spread of HCV in body cells such as lifestyle (alcohol consumption or smoking), patient compliance level, lipid metabolism, metabolic syndrome, and body weight ^[23]. Motivated by this phenomenon, we are interested in extending the mathematical model given in ^[24] by considering uncertainty factors. A stochastic model will be derived from the deterministic model. This model provides better insight into the uncertainty and variability of the disease dynamic. Moreover, the solution of the stochastic model is in the form of distribution^[25], while the solution of the deterministic model only produces one predictive value ^[26]. Several papers also discussed the stochastic HCV model, e.g., ^{[27,28,29,30,31,32]}. In contrast, in this paper, the noise parameter is added in the transmission rate representing the characteristic of variability and the treatment in the model is considered.

The outline of this paper is organized as follows. In Section 2, the HCV model involving stochastic disturbances in the transmission rate are developed. In addition, the existence, uniqueness, and boundedness of the solution are established. In Section 3, we derive the extinction and persistence in mean condition. In Section 4, the numerical simulations are carried out to illustrate the analytical results. The conclusions are presented in the last section.

2.   Materials and methods

2.1.   Mathematical model

Recall a deterministic mathematical model of HCV given in ^[24] as

where $T$ is the number of uninfected cells, $I$ is the number of infected cells, and $V$ is the number of free viruses. The uninfected cells are produced at rate $\Lambda$ and die naturally at a constant rate $\delta_1$. Cells become infected when they are interacting with a virus with a constant rate $\beta$, and once infected, they will die at a constant rate $\delta_2$. HCV is produced by infected cells at a constant rate $k$ and cleared at a constant rate $c$.

For the deterministic model $(2.1)$, the disease-free equilibrium point is $E_{0} = \left(\frac{\Lambda}{\delta_{1}}, 0, 0\right)$ and the endemic equilibrium point is $E_{1} = \left(\frac{\delta_{2}c}{(1-\eta)(1-\epsilon)k\beta}, \frac{\Lambda}{\delta_{2}}-\frac{\delta_{1}c}{(1-\eta)(1-\epsilon)k\beta}, \frac{(1-\epsilon)k\Lambda}{\delta_{2}c}-\frac{\delta_{1}}{(1-\eta)\beta}\right).$ Using the next-generation matrix method ^[33], the basic reproduction number of the system $(2.1)$ is $R_{0}^{d} = \frac{(1- \epsilon)(1-\eta)\beta k \Lambda}{c \delta_{1} \delta_{2}}$.

In this paper, we generalize system (2.1) by incorporating a stochastic noise parameter in the transmission rate as

where $B(t)$ is a standard Brownian motion and $\sigma$ is a real constant which is known as the intensity of noise. The value $\sigma$ is the standard deviation of transmission rate data that represents the variability of the transmission rate. The other description of parameters model (2.2) is given in Table 1. Clearly, the system (2.1) is the special case of the system (2.2), where $\sigma = 0.$

2.2.   Existence and uniqueness of the solution

In this section, the theorem of the existence, uniqueness, and boundedness of the solution system $(2.2)$ are established. Let $(\Omega, \mathcal{F}, P)$ be a complete probability space with a filtration $\left\{\mathcal{F}_t\right\}_{t\geq 0}$ satisfying the conditions that it is right continuous and $\mathcal{F}_0$ contains all P-null sets.

Definition 2.1. ^[34]

Let $T > 0$, $F(., .):[0, T]\times \mathbb{R}^n \rightarrow \mathbb{R}^n, G(., .):[0, T]\times \mathbb{R}^n \rightarrow \mathbb{R}^{n\times m}$ be measurable function, and $X(t)$ satisfy

where F(., .) and G(., .) are the coefficients of (2.3). Then the coefficients (2.3) are locally Lipschitz, if

for some constant $C_{1}\geq 0$.

Definition 2.2. ^[34]

The coefficients of (2.3) satisfy linear growth condition, if

for some constant $C_{2}$.

Lemma 2.3. The system (2.2) with the initial condition in ${\mathbb{R}}^3 _+$ is uniformly ultimately bounded and belongs to the following closed and bounded positively invariant set

for every $t\geq 0.$

Proof. Let $a(t)$ be a function such that $a(t)\le \Lambda$. We know that $N(t) = T(t)+I(t)$ is the concentration of uninfected and infected cells at time $t$. Take $\mu = \min \{\delta _1, \delta _2 \}$, then

Let $dN = (a(t) - \mu N) dt$, it follows that

thus,

Since $a(t) \leq \Lambda$, we obtain

By taking a limit $t \to \infty$ of Eq ($2.9$), we get

Analog with the technique for solving Eq (2.7), for $V(0) = V_0$ and $I(t) < N(t) \leq \frac{\Lambda}{\mu }$, we get

then

If we take $\mathop {\lim }\limits_{t \to \infty } V(t)$, then

Since the right-hand side of Eq (2.2) satisfies the Lipschitz condition, the solution exists and unique on $[0, b)$ for some $b > 0$. Assume that there exists $t_{1} \in (0, b)$ such that $V(t_{1}) = 0$ and all other variables are positive on $t_{1} \in (0, b)$. Therefore, for all $t \in [0, t_{1}]$

thus

This is contradiction with $V(t_{1}) = 0$. Then $V(t) > 0.$ Analoguely, we get $I(t) > 0$ and $T(t) > 0.$ Therefore, we obtain that $(T(t), I(t) \in \left(0, \frac{\Lambda}{\mu}\right)$ for all $t \in [0, T]$ and $0 < V(t)\leq (1 - \varepsilon)\left({\frac{{k \Lambda}}{{c\mu }}} \right).$

Next, we define the necessary condition that guarantees the existence and uniqueness of time-global solution of Eq (2.2).

Theorem 2.4. Let the coefficients of the system ($2.2$) satisfy the Lipschitz condition. Then for any initial value $(T(0), I(0), V(0)) \in \Gamma$, there exists a unique time-global solution $(T(t), I(t), V(t)) \in \Gamma \subset {\mathbb{R}}^3_+, t \ge 0$ with probability $1$.

Proof. According to Definition 2.1, Definition 2.2, and Theorem 5.2.1 in ^[34], there exists a unique global solution. However, we have that the coefficients of the system ($2.2$) only satisfy the Lipschitz condition ^[35], thus system ($2.2$) has a unique local solution on $t\in [0, \tau_e)$ for any initial value $(T(0), I(0), V(0))\in \Gamma$ where $\tau_e$ is the explosion time (i.e. the time when the solution tends to infinity). To guarantee the solution of system ($2.2$) is a unique global solution, it is necessary to show that $\tau_{e} = \infty$ ^[36,37].

Let $k_0 > 0$ be sufficiently large such that every component of $(T(0), I(0), V(0))$ is in the interval $\left[ {\frac{1}{{{k_0}}}, {k_0}} \right]$. For every $k \ge k_0$, we define the stopping time or the first passage time (i.e. first period when the stochastic process penetrates the barrier) as

Throughout this paper, we set $\inf \emptyset = \infty$. It is known that the lower bound of $\mathbb{R}$ is an empty set and the largest infimum of the empty set is infinity. Since $\tau_k$ is increasing as $k \to \infty$, and

then it follows that

thus $\tau_\infty \le \tau_e$ almost sure. Next, it is necessary to show

We will prove by contradiction. Suppose that $P(\tau_\infty < \infty) < 1$. If $\tau_\infty < \infty$, then there exist $T^* > 0$ and $\varepsilon \in (0, 1)$ such that $P\left\{ {{\tau _k} \le T^*} \right\} > \varepsilon, \forall k \ge k_0$. In this case, the proof technique analogue is to ^[37,38,39]. Define a function ${C^{2, 1}}, Q:{[0, \infty) \times \mathbb{R}}_ + ^3 \to {{\mathbb{R}}_ + }$ where

Since

then Q is a non-negative function. By using It$\hat{o}$'s formula, we get

where $M$ is a positive constant. Integrating Eq (2.12) from $0$ to ${{\tau _k} \wedge T^*}$ yields

Taking the expectation of both sides in Eq (2.13) leads to

Define ${\Omega _k} = \left\{ {{\tau _k} \le T^*} \right\}$ for any $k \ge k_0$. Then for every $w\in \Omega_k$ there are components of $Q\left({T({\tau _k} \wedge T^*), I({\tau _k} \wedge T^*), V({\tau _k} \wedge T^*)} \right)$ which are equal to either $k$ or $\frac{1}{k}$, thus

Combining Eqs (2.14) and (2.15) yields

where ${1_{{\Omega _k}}}$ is the indicator function. By taking a limit $k \to \infty$ of Eq (2.16), we obtain

This is contradiction. We get $\tau _\infty = \infty$ or $P(\tau _\infty = \infty) = 1.$ Therefore, there exists a unique time-global solution $(T(t), I(t), V(t)) \in\Gamma \subset {\mathbb{R}}^3_+, t \ge 0$ with probability $1$.

3.   The equilibria and their stabilities

In this section, the theorem of the extinction and persistence in the mean condition system (2.2) is derived.

3.1.   Extinction

In this part, we investigate the almost surely exponential stability of the disease-free equilibrium point by using the suitable Lyapunov function and another technique of stochastic analysis.

Lemma 3.1. If

the disease-free equilibrium point $(T, I, V) = (\Lambda/\delta_{1}, 0, 0)$ is almost surely exponentially stable in $\Gamma$.

Proof. The technical proof of this theorem follows from ^[38,40]. Define a function $Q = \ln (I+V)$. Using It$\hat{o}$'s formula, we obtain

where $Z = \frac{VT}{I+V}.$ Futhermore, we have

Taking the integral of ($3.2$) and dividing both sides by t and then computing the limit superior $t\rightarrow \infty$ yield

By the strong law of large number for Martingales ^[41], we have

Thus, from Eqs (3.3) and (3.4), we obtain

By Lemma 3.1, if the noise is increasing and satisfies the condition (2.10) HCV will die out. However, in this study, we also get that a bounded variation of infection rate could also lead to extinction. This is presented in the following theorem.

Theorem 3.2. Let $(T(0), I(0), V(0)) \in \Gamma$ be the initial value of system (2.2) and $(T(t), I(t), V(t)) \in \Gamma$ be the corresponding solution. If

and $\sigma^{2} < \frac{4 \delta_{1} \delta_{2} \beta }{(1-\epsilon) (1- \eta) k \Lambda}$, then HCV will be extinct almost surely, i.e.,

Proof. The proof technique follows from ^{[38,40,42,43]}. Consider the first equation of system (2.2), then we get

or

For simplicity, we define an integrable function $\langle X \rangle = \frac{1}{t} \int_0^t X(s) ds$, $X(t)\in [0, +\infty)$, thus

Taking the limit of Eq ($3.5$) as t tend to infinity, we obtain

Next, define a function $Q = \ln ((1- \epsilon)kI+ \delta_{2} V)$. Applying Ito's formula, we obtain

Integrating both sides of Eq ($3.6$) from 0 to t yields

Dividing both sides of Eq ($3.7$) by $t$ and taking a limit superior for $t \rightarrow 0$, we get

By the strong law of large number for Martingales ^[41],

we obtain

Thus, $\mathop {\lim }\limits_{t \to \infty } {I(t)} = 0$ and $\mathop {\lim }\limits_{t \to \infty } {V(t)} = 0.$ In other words, HCV will be extinct.

3.2.   Persistence in mean

In this subsection, we give a sufficient condition to guarantee the persistence in the mean condition.

Definition 3.3. ^[31,40,44]

System (2.2) is said to be persistent in mean if

almost surely.

Theorem 3.4. Let $(T(0), I(0), V(0)) \in \Gamma$ be the initial value of system (2.2) and $(T(t), I(t), V(t)) \in \Gamma$ be the corresponding solution. If

then HCV will be persistent in mean, i.e.,

Proof. Analogue to the proof technique in ^[40] and ^[44]. Integrating the system ($2.2$) yields

Dividing both sides Eq (3.8) by $t$, we get

where

Define a function $Q = \ln [(1-\epsilon)k I(t) + (1+\delta_{2})V(t)].$ Applying It$\hat{o}$'s formula, we obtain

Integrating Eq (3.10) from 0 to t and dividing both sides by t, we have

According to Eqs (3.9) and (3.11), we get

where

Taking the limit inferior of Eq (3.12) as $t \rightarrow 0$ yields

This completes the proof.

4.   Numerical simulation

In this section, we carry out a numerical simulation in order to provide an interpretation of the solution. We use the Euler-Maruyama method to determine the solution of system (2.2). The discretization of system (2.2) is given as follows

where $t\in [t_{0}, t_{N}], \Delta t = \frac{t_{N}-t_{0}}{N}$, and $\zeta$ is normally distributed N(0, 1).

Three simulations are conducted. The first simulation is deterministic system. The second simulation is to represent the solution of a stochastic model in the extinction condition. The third simulation is illustrating the solution of a stochastic model when Lemma 3.1 and Theorem 3.2 are violated. The parameters and initial values of the model in each simulation are given in Table 2.

Figure 1 shows the solution of a deterministic model ($\sigma = 0$) with $R_{0}^{d} > 1$. It can be seen that the solution is tend to the endemic equilibrium point $E_{1} = (6.6\times 10^6, 2.6\times 10^6, 7\times 10^5)$ which means that HCV disease is persistent.

Figure 2 illustrates the stochastic model with $\sigma = 4\times 10^{-7}$ the condition

and $R_{0}^{s} < 1$ where $\sigma^{2} = 1.6\times10^{-13} < \frac{4 \delta_{1} \delta_{2} \beta }{(1-\epsilon) (1- \eta) k \Lambda} = 5.4\times 10^{-13}$ satisfies Lemma 3.1 and Theorem 3.2 respectively. More precisely, when the choosen parameters satisfy Lemma 3.1 and Theorem 3.2, the point $E_{0} = (8\times 10^6, 0, 0)$ is almost exponentially stable. It can be seen from Figure 2 where the concentration of infected cells and free viruses go to zero. Therefore, HCV will go extinct almost surely.

Figure 3 describes the condition when basic reproduction number $R_{0}^{s} = 0.0214 < 1$ and $\sigma^{2} = 4.9\times10^{-15} < \frac{4 \delta_{1} \delta_{2} \beta }{(1-\epsilon) (1- \eta) k \Lambda} = 5.4\times 10^{-13}$. This is satisfying Theorem 3.2 and violating Lemma 3.1. That figure shows a stochastic noise impact the disease extinction.

Figure 4 shows that if the intensity of the disturbance is smaller where $\sigma^{2} = 8.1\times10^{-16} < 1.9\times10^{-12}$, thus Lemma 3.1 is violated. Taking the time interval [0,400] with 4000 points and 10 paths, it appears in Figure 4(a) that the concentration of uninfected cells rises to a peak around at t = 30 days, then drops sharply to around t = 45 days and then they rise slowly. The graphics in Figures 4(b), (c) show the concentration of infected cells and free virus after t = 50 days which tend to certain interval values. This shows that hepatitis C remains exist.

Figure 5 describes the condition of the HCV extinction when Lemma 3.1 is satisfied while Theorem 3.2 is violated. It can be seen that figures show the concentration of infected cells and free virus from initial value that tends to zero for a long period of time.

Figure 6(a) shows the distribution value of uninfected cells in simulation 2 with a confidence interval of 95 %. After time $t$ about 50 days the concentration of uninfected cells will go to a positive constant value ($8\times 10^6$). From Figure 6(b), the simulation at t = 5 is more diverse than t = 10. After time $t$ 20 days, the concentration of infected cells will go to zero. Figure 6(c) shows the free virus concentration at the time $t$ from 3 to 15 days is in the certain interval value. After time $t$ 20 days, the virus concentration will go to zero. Therefore, from simulation 2, it can be concluded that stochastic disturbances can affect the behavior of the spread of hepatitis C. When conditions are met according to Lemma 3.1 and Theorem 3.2, at a certain time t hepatitis C will disappear. Mathematically, in this case, the basic reproduction number in the extinction condition is less than one.

Figure 7 shows the confidence intervals for simulation 3 where we reduced the sigma values while preserving the other parameter values. Figure 7(a) shows the concentration of uninfected cells will reach peaks around 23 to 30 days and then decline. After 75 days the concentration will move to a certain interval value. The green lines show the mean value of the simulations. Figures 7(b), (c) show that after 75 days, the concentration of infected cells and free virus lead to a certain positive number interval value with a confidence interval of approximately 95 %. In other words, hepatitis C disease will persist.

Lastly, we conduct a simulation for persistence in mean condition where $\beta = 5.4\times 10^{-7}, \delta_{1} = 0.0047, k = 1, c = 0.6, \sigma = 9 \times 10^{-8}$, and the rest of parameters and initial values state in Table 2. Figure 8 illustrates the solution of the stochastic model compared to the deterministic. Figure 9 shows the number of concentrations of healthy cells, infected cells, and free virus at a certain value which are more than zero. When $t \geq 50$, it appears that the concentrations of uninfected cells, infected cells and free viruses go to an equilibrium state. According to the numerical calculations, we obtain $R_{1}^{s} > 1$ and the endemic equilibrium point $E_{1} = (T, I, V) = (4.4\times 10^6, 2.6\times10^6, 1.08\times 10^6)$. According to Theorem 3.4 hepatitis C remains to exist which are also depicted in Figure 9. Figure 10 describes the condition if the noise intensity is increasing, then the solution of model $(2.2)$ will be strongly oscillating around the endemic equilibrium point of the model $(2.1)$.

Furthermore, in Figure 11, we generate a 95% confidence interval for $t \in [0,400]$ with 4000 paths. Figure 11(a) describes the number of uninfected cells that tends to a consistent interval value after 75 days. Figure 11(b), (c) show the confidence intervals of infected liver cells and free viruses. In the time t span of 0 to 400 days the number of concentrations in a row increases to a peak around day 27 for infected cells and free viruses, then decreases. After 30 days, leading to a consistent value for infected cells from $1.9 \times 10 ^ 6$ to $3.9 \times 10 ^ 6$, and free viruses from $0.75 \times 10 ^ 6$ to $1.5 \times 10 ^ 6$.

5.   Conclusions

In this paper, we propose an epidemic model for HCV transmission at the cellular level incorporating statistical noise. The results extend the paper of Neumann, et.al.^[24] in understanding the dynamics of a deterministic hepatitis C virus. Various types of Lyapunov functions are designed to study the persistence in mean conditions as well as the extinction of the stochastic system. Based on this model, there exists a unique time-global solution for any given positive initial value. In addition, numerical simulations are carried out to describe the solution behaviour of the model. We analyze that if the noise intensity is increasing, then the disease will go to extinction. If the basic reproduction number of the extinction condition is less than one, then hepatitis C will extinct. If the basic reproduction number of the persistence in the mean condition is more than one, then hepatitis C remains to exist. When the noise intensity increases, the solution will give a strong oscillation under the condition of persistence in the mean. It can provide some useful strategies for controlling the dynamics of that disease. In future research, we can consider an optimal control of the stochastic model for HCV, research on memory effect, and fractional derivatives ^[47]. Furthermore, it needs to investigate the stochastics noise for other parameters apart from infection rate. It is also important to consider a stochastic model for HCV transmission involving the response immune system.

Acknowledgments

The first author thanks LPDP Indonesia for the financial support under the Doctoral Program Scholarship. The authors thanks to Department of Mathematics, Universitas Gadjah Mada (UGM) and Directorate for the Higher Education, Ministry of Research, Technology, and Higher Education of Indonesia for the Research Grant "Penelitian Disertasi Doktor (PDD)" UGM 2021 under grant no. 2261/UN1/DITLIT/DIT-LIT/PT/2021. The special thanks also for some colleagues for the discussions during the research.

Conflict of interest

The authors declare that they have no competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.