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 Λ and die naturally at a constant rate δ1. Cells become infected when they are interacting with a virus with a constant rate β, and once infected, they will die at a constant rate δ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 E0=(Λδ1,0,0) and the endemic equilibrium point is E1=(δ2c(1−η)(1−ϵ)kβ,Λδ2−δ1c(1−η)(1−ϵ)kβ,(1−ϵ)kΛδ2c−δ1(1−η)β). Using the next-generation matrix method [33], the basic reproduction number of the system (2.1) is Rd0=(1−ϵ)(1−η)βkΛcδ1δ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 σ is a real constant which is known as the intensity of noise. The value σ 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 σ=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 (Ω,F,P) be a complete probability space with a filtration {Ft}t≥0 satisfying the conditions that it is right continuous and F0 contains all P-null sets.
Definition 2.1. [34]
Let T>0, F(.,.):[0,T]×Rn→Rn,G(.,.):[0,T]×Rn→Rn×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 C1≥0.
Definition 2.2. [34]
The coefficients of (2.3) satisfy linear growth condition, if
for some constant C2.
Lemma 2.3. The system (2.2) with the initial condition in R3+ is uniformly ultimately bounded and belongs to the following closed and bounded positively invariant set
for every t≥0.
Proof. Let a(t) be a function such that a(t)≤Λ. We know that N(t)=T(t)+I(t) is the concentration of uninfected and infected cells at time t. Take μ=min{δ1,δ2}, then
Let dN=(a(t)−μN)dt, it follows that
thus,
Since a(t)≤Λ, we obtain
By taking a limit t→∞ of Eq (2.9), we get
Analog with the technique for solving Eq (2.7), for V(0)=V0 and I(t)<N(t)≤Λμ, we get
then
If we take limt→∞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 t1∈(0,b) such that V(t1)=0 and all other variables are positive on t1∈(0,b). Therefore, for all t∈[0,t1]
thus
This is contradiction with V(t1)=0. Then V(t)>0. Analoguely, we get I(t)>0 and T(t)>0. Therefore, we obtain that (T(t),I(t)∈(0,Λμ) for all t∈[0,T] and 0<V(t)≤(1−ε)(kΛcμ).
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))∈Γ, there exists a unique time-global solution (T(t),I(t),V(t))∈Γ⊂R3+,t≥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∈[0,τe) for any initial value (T(0),I(0),V(0))∈Γ where τ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 τe=∞ [36,37].
Let k0>0 be sufficiently large such that every component of (T(0),I(0),V(0)) is in the interval [1k0,k0]. For every k≥k0, 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∅=∞. It is known that the lower bound of R is an empty set and the largest infimum of the empty set is infinity. Since τk is increasing as k→∞, and
then it follows that
thus τ∞≤τe almost sure. Next, it is necessary to show
We will prove by contradiction. Suppose that P(τ∞<∞)<1. If τ∞<∞, then there exist T∗>0 and ε∈(0,1) such that P{τk≤T∗}>ε,∀k≥k0. In this case, the proof technique analogue is to [37,38,39]. Define a function C2,1,Q:[0,∞)×R3+→R+ where
Since
then Q is a non-negative function. By using Itˆo's formula, we get
where M is a positive constant. Integrating Eq (2.12) from 0 to τk∧T∗ yields
Taking the expectation of both sides in Eq (2.13) leads to
Define Ωk={τk≤T∗} for any k≥k0. Then for every w∈Ωk there are components of Q(T(τk∧T∗),I(τk∧T∗),V(τk∧T∗)) which are equal to either k or 1k, thus
Combining Eqs (2.14) and (2.15) yields
where 1Ωk is the indicator function. By taking a limit k→∞ of Eq (2.16), we obtain
This is contradiction. We get τ∞=∞ or P(τ∞=∞)=1. Therefore, there exists a unique time-global solution (T(t),I(t),V(t))∈Γ⊂R3+,t≥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)=(Λ/δ1,0,0) is almost surely exponentially stable in Γ.
Proof. The technical proof of this theorem follows from [38,40]. Define a function Q=ln(I+V). Using Itˆo's formula, we obtain
where Z=VTI+V. Futhermore, we have
Taking the integral of (3.2) and dividing both sides by t and then computing the limit superior t→∞ 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))∈Γ be the initial value of system (2.2) and (T(t),I(t),V(t))∈Γ be the corresponding solution. If
and σ2<4δ1δ2β(1−ϵ)(1−η)kΛ, 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 ⟨X⟩=1t∫t0X(s)ds, X(t)∈[0,+∞), thus
Taking the limit of Eq (3.5) as t tend to infinity, we obtain
Next, define a function Q=ln((1−ϵ)kI+δ2V). 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→0, we get
By the strong law of large number for Martingales [41],
we obtain
Thus, limt→∞I(t)=0 and limt→∞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))∈Γ be the initial value of system (2.2) and (T(t),I(t),V(t))∈Γ 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−ϵ)kI(t)+(1+δ2)V(t)]. Applying Itˆ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→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∈[t0,tN],Δt=tN−t0N, and ζ 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 (σ=0) with Rd0>1. It can be seen that the solution is tend to the endemic equilibrium point E1=(6.6×106,2.6×106,7×105) which means that HCV disease is persistent.
Figure 2 illustrates the stochastic model with σ=4×10−7 the condition
and Rs0<1 where σ2=1.6×10−13<4δ1δ2β(1−ϵ)(1−η)kΛ=5.4×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 E0=(8×106,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 Rs0=0.0214<1 and σ2=4.9×10−15<4δ1δ2β(1−ϵ)(1−η)kΛ=5.4×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 σ2=8.1×10−16<1.9×10−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×106). 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 β=5.4×10−7,δ1=0.0047,k=1,c=0.6,σ=9×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≥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 Rs1>1 and the endemic equilibrium point E1=(T,I,V)=(4.4×106,2.6×106,1.08×106). 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∈[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×106 to 3.9×106, and free viruses from 0.75×106 to 1.5×106.
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.