Deterministic stability and random behavior of a Hepatitis C model

The deterministic stability of a model of Hepatitis C which includes a term defining the effect of immune system is studied on both local and global scales. Random effect is added to the model to investigate the random behavior of the model. The numerical characteristics such as the expectation, variance and confidence interval are calculated for random effects with two different distributions from the results of numerical simulations. In addition, the compliance of the random behavior of the model and the deterministic stability results is examined.


Introduction
Hepatitis C is an infectious liver disease. The virus which causes this disease was identified in 1989 but the worldwide presence of the virus shows that it has been active for a much longer period. It is estimated that around 150 million people are chronically infected by the World Health Organization (WHO) reports [1]. Hepatitis C Virus (HCV) causes 3-4 million new infections per year. HCV infections occur in two basic stages, acute and chronic infections. The terms 'acute' and 'chronic' refer to the duration of the disease and not the severity. The illness can range from a mild illness which lasts less than a month to serious infections which can last several months, and even a lifetime in some cases [1]. The acute stage of the disease is largely asymptomatic and about one fifth of these cases resolve spontaneously due to adequate response to the HCV by the immune system. Less than one fifth of acute infections show mild symptoms like fatigue and jaundice. Infections that last up to 6 months are called 'acute' and acute infections have 1% mortality rate [2,3]. Those that last longer are called 'chronic'. Nearly 80% of HCV infections develop into the chronic stage which can last asymptomatic for more than 20 years [2,4], [1]. The long period of asymptomatic infections makes the diagnosis of the disease difficult therefore hepatitis C is sometimes called the 'silent epidemic' [2,3]. In about 30 years, more than one fifth of the infected develop cirrhosis and 1%-3% develop lung cancer. More than 300,000 people die yearly from diseases that are related with HCV [1,2].
Stability analysis of the equilibrium points of the system provides a better understanding of the behavior of the system in a long period of time without the need to find the solutions of the model. Local stability of an equilibrium point suggests that if the system is close to this point, it PLOS  will eventually reach the equilibrium state at this point. The global stability of an equilibrium point suggests that the system will reach equilibrium state at this point, whether it is close to this point or not [5]. Considering models in medicine, biology and etc., global stability of equilibrium points can indicate extinction or persistence of the disease according to the equilibrium point under consideration [6]. Local stability analysis examines the effects of small variations in each of the variables on the results of the model. Hence, a motivation for the random analysis of the model containing random effects in the parameters is to visualize the random behavior of the model, which can be linked with the stability of the equilibrium points under small changes in the conditions of the system. The behavior of the solution of Hepatitis C virus model is examined by Ahmed and El-Saka [7]. The model in [7] is a fractional model which compares the results of the system for various powers of differentiation. Fractional calculus, with the use of its memory effect property, may provide useful information on the results of the model. However, we concentrate on the system of ordinary differential equations (α = 1), since we want to add random effects to the parameters and investigate the randomness of the event. Similar dynamical modeling studies can also be reviewed for an investigation of the framework of disease transmission models consisting of ordinary differential equations [8][9][10][11][12]. It should also be noted that the use of spatial effects may also provide useful results for analyzing Hepatitis C transmission dynamics [13][14][15][16]. These models use mathematical tools to guide and enhance studies in medicine, biology and etc. and with the addition of random effects, we intend to extend these analyses with the addition of a statistical point of view.
The components of the basic three-component model are uninfected hepatocytes, infected hepatocytes and the virus, which are denoted by T(t), I(t) and V(t), respectively. The flowchart of this model can be visualized as Fig 1, which has been obtained by a modification of the flowchart in [17].
The parameters of the model describe the rates of change in the uninfected hepatocytes, infected hepatocytes and virus during treatment. s describes the constant production rate of uninfected hepatocytes per cell, while d and β describe their constant death and infection rates per cell. δ is the constant rate of death per cell for infected hepatocytes. p is the constant rate of production for viral particles per infected hepatocyte and c is their constant clearance rate per virion. The treatment effects are included in the model with two parameters p and η, which describe the virion production blockage and new infection reduction respectively. For example, p = 0.95 indicates 95% efficacy in blocking of virion production. c 2 describes the rate of immunity response.
These parameters must be assigned values for the simulations of the Model (1). Values and descriptions of the parameters of model, which were obtained from [7] are explained in Table 1. Simulations are done with these parameter values above and the initial conditions T (0) = 2.4 × 10 6 , I(0) = 2.0 × 10 6 and V(0) = 4.0 × 10 5 . Study motivation of this work is covered by the earlier study of Merdan and Khaniyev [18].

Basic properties of the model
Firstly, the basic reproduction number and some properties of the model are examined. The disease-free equilibrium (DFE) for the System (1) The basic reproduction number is used for the analysis of the spread and control of the disease mathematically. It indicates whether the disease will spread through the community or be taken under control, which is very useful information. We will use the formula from [17,[19][20][21][22] to calculate this number for System (1). A similar study [23] can also be reffered to for the calculation of the equilibrium points of a similar model for HCV transmission. Let X = [I, V, T] T , then for System (1) ; WðX Þ ¼ dI À ð1 À p ÞpI þ cV À s þ dT þ ð1 À ZÞbVT : Jacobian of F at DFE and Jacobian of W are obtained respectively as follows: A The inverse of V is used with F(X) to obtain; The spectral Radius R 0 is , as found in [20]. [24] can also be referred to for the calculation of the basc reproduction number. Proposition 1. If R 0 > 1, then there exists positive equilibria for the System (1) and one of these is the endemic equilibrium Let all the equations of System (1) equal to zero, then Substituting Eq (3) into the first equation of Eq (2), we get is a quadratic function, and notice that where R 0 > 1. Thus, f(I) = 0 has positive form and two roots. Therefore, the endemic equilibriums of the System (1) are given by not have any positive roots. However it can be seen that I Ã 2 becomes negative for every value of R 0 , thus the equilibrium point E 2 is biologically irrelevant and should be ignored.

Local stability
In this part of the study, the local stability of the model is investigated for the equilibrium points.
Theorem 1. For the disease-free equilibrium point E 0 of the System (1), we have the following: Proof. When evaluated at the point E 0 , System (1) has the Jacobian matrix The characteristic equation of J 0 is given by The Routh-Hurwitz criterion for the cubic equation is as follows: Thus, the Routh-Hurwitz criterion is satisfied. So, the System (1) is locally asymptotically stable in the neighborhood of the disease-free equilibrium (DFE) E 0 . If R 0 > 1, then Q 2 < 0, implying that E 0 is unstable. The referred study [20] can be reviewed for further notes on the stability of the disease-free equilibrium.
Theorem 2. The endemic equilibrium point (1) is locally asymptotically stable when R 0 > 1 and unstable otherwise.
Proof. When evaluated at the point E 1 , System (1) has the Jacobian matrix The characteristic equation of J 1 is given by The Routh-Hurwitz criterion for the cubic equation is as follows: Thus, by Routh-Hurwitz criterion, the endemic equilibrium point

Global stability
The global stabilities of the disease-free equilibrium E 0 and the endemic equilibrium E 1 of the Hepatitis C Virus transmission model are investigated in this section [25][26][27]. Theorem 3. The disease-free equilibrium E 0 is globally asymptotically stable in the positively variant set O for R 0 < 1.
Proof. Consider the following Lyapunov function U(t) = (1 − p )pI + δV [28]. Calculating the derivative dUðtÞ dt for the solutions of System (1), we get Provided that R 0 < 1, we find that dUðtÞ dt 0, which, considering LaSalle's invariance principle [29], indicates that the disease-free equilibrium E 0 is globally asymptotically stable.
In the following, we deal with the geometric approach developed in [27] for the proof of the global stability of endemic equilibrium point E 1 . A simple sufficient condition guarantees the global asymptotical stability of the epidemic equilibrium E 1 . We begin with a summary of the geometric approach.
Consider the autonomous dynamical system: where f:  (4). If x Ã is stable then it is globally stable in D.
Thus, p f is the matrix given by, and J [2] is the second additive compound matrix of the Jacobian matrix (J(x) = Df(x)), while μ (B) is the Lozinski measure of B with respect to a vector norm |.| in R n ; N ¼ C 2 n and is defined by The following result for global stability is proved by [27]. Lemma 2. ( [30]) Suppose that D is simply connected and that (H 2 ) and (H 3 ) are satisfied. Then the unique equilibrium The following result can be found with the lemmas above.
Proof. Firstly, we consider the Jacobian matrix of System (1) Its second additive compound matrix is  : Let (a 1 , a 2 , a 3 ) be a vector in R 3 . Its norm L 1 k . k is defined as k ða 1 ; a 2 ; a 3 Þ k¼ maxfja 1 j; ja 2 j þ ja 3 jg: Denote the Lozinski measure with respect to this norm by μ(B). It follows from the notation in [31] that we have μ(B) sup{g 1 , g 2 }, where |B 12 |, |B 21 | are matrix norms with respect to the L 1 vector norm, and μ 1 denotes the Lozinski measure with respect to the L 1 norm. Then where h ¼ minfd þ ð1 À ZÞbV; d À 2dI c 2 g > 0. Therefore, we have From the equation System (1), we get Then we have Furthermore, we obtain mðBÞ supfg 1 ; As a result, endemic equilibrium

Random behavior of the model
We investigate the behavior of the randomized components to visualize the effects of small variations on the model output. The parameters in the deterministic System (1) are added random effects to investigate the random characteristics of the model. A random analysis of models using random differential equations with random parameters has been used before for a model of avian-influenza by Merdan and Khaniyev in 2008 and for bacterial resistance by Merdan et al., which was the motivation for the method used in this work [18,32]. The random model is analyzed to investigate the numerical characteristics of the event and thus comment on the random behavior of the model components. The randomness in the parameters of the model can be linked to the stability of the deterministic model since the stability of the equilibrium points are essentially the ability of the system of maintain its position on these points under small variations. We use normally and symmetrically distributed random effects to model real lifer variations in the parameters of the model. Normal distribution is commonly used for random variables with unknown distributions since the central limit theorem states that a sufficiently large number of independent random variables will be approximately normally distributed under certain conditions. A triangularly distributed random variable has a high probability around its mean and the probability decreases for values that are far away from the mean. A symmetrical triangular distribution and a normal distribution were used for the random effects since they are similar in the above mentioned sense and hence a comparison can be made for the randomness of the results.

Investigation of the model under normally varying random effects
The parameters s, d, η, β, δ, c 2 , p , p and c are considered to be random variables with normal distribution in order to investigate the model under normally distributed random effects.
Using the "randn" command in MATLAB, which generates random variables with standard normal distribution, we can generate random variables s, d, η, β, δ, c 2 , p , p and c which have normal distribution. These generated random variables will have the following forms: c = c 0 + σ 1 η 1 , d = d 0 + σ 2 η 2 , β = β 0 + σ 3 η 3 , δ = δ 0 + σ 4 η 4 , s = s 0 + σ 5 η 5 , p = p 0 + σ 6 η 6 , η = η 0 + σ 7 η 7 , p = p 0 + σ 8 η 8 , c 2 = c 2 0 + σ 9 η 9 , where the random variables Z i ; i ¼ 1; 9 are independent random variables with standard normal distribution and s i ; i ¼ 1; 9 are the corresponding deviations used for each of the parameters s, d, η, β, δ, c 2 , p , p and c. The following values will be used for the deviations of parameters: As it can be seen in the list above, normally distributed random variables, which are denoted as Z i ; i ¼ 1; 9, are added to the initial values of the parameters s, d, η, β, δ, c 2 , p , p and c, which are denoted by the zero-indexed parameter, therefore forming the new parameters which have normal distribution. The deviations of the random effects are determined to be powers of 10, so that the random effect added to the initial values of the parameters is around 1% to 4.4% of the initial value.
The random system produces deterministic differential equations which are assigned different coefficients for every trial of the event.
Note that a random variable is a measureable function from the set of possible events to real numbers R meaning it is a real valued function. Hence for every trial of the event, the random variables produce different real numbers according to their random distribution. For the random model, this would mean that we would get a different set of differential equations for various trials which would all be deterministic differential equations with different valued coefficients. Since the random model produces deterministic equations with variations in the set of parameters, the random analysis of the behavior of model components will be based on the statistical properties of the numerical solutions of these deterministic equations. Using 10 5 simulations of the numerical solutions of the equations, comments are made on the numerical characteristics of the model with random coefficients.

Expectations. T(t), I(t)
and V(t) under random effects are random variables. Hence, their moments can be evaluated using the law of large numbers.  Fig 2. Expectations of T(t), I(t) and V(t).

., N are the results obtained from the simulation of the process T(t). The graphs of E(T(t)), E(I(t)) and E(V(t)) are given in Fig 2.
https://doi.org/10.1371/journal.pone.0181571.g002

Confidence intervals. The confidence intervals of T(t), I(t) and V(t) are found in the form of [E(T(t)) − Kσ(T(t)), E(T(t)) + Kσ(T(t))
] by using the expected values and standard deviations which were previously calculated. For K = 3, confidence intervals are given at approximately 99%, meaning that there is 99% probability that the given interval includes the real value of T(t). The confidence intervals can be seen in

Investigation of the model under triangularly varying random effects
The parameters are considered to be random variables with symmetrical triangular distribution in the interval (−1, 1) in order to investigate the model. Using the property above and the 'rand' command in MATLAB, which generates uniformly distributed random variables from (0, 1), we can generate random variables which have symmetrical triangular distribution in the interval (−1, 1). A similar modeling approach will produce the system under triangular effects as: þð1 À Z 0 þ 0:01ðZ 71 À Z 72 ÞÞðb 0 þ 0:00000001ðZ 31 À Z 32 ÞÞVT; dI dt ¼ ð1 À Z 0 þ 0:01ðZ 71 À Z 72 ÞÞðb 0 þ 0:00000001ðZ 31 À Z 32 ÞÞVT À ðd 0 þ 0:01ðZ 41 À Z 42 ÞÞI 1 À I c 2 0 þ 100000 Â ðZ 91 À Z 92 Þ ! ; dV dt ¼ ð1 À p 0 þ 0:01ðZ 61 À Z 62 ÞÞðp 0 þ 0:1ðZ 81 À Z 82 ÞÞI À ðc 0 þ 0:1ðZ 11 À Z 12 ÞÞV: Here, the random variables Z ij ; i ¼ 1; 9; j ¼ 1; 2 are uniformly distributed independent random variables from (0, 1), so that their difference can produce independent symmetrical A similarity between the results of the expectations for the normally and triangularly distributed effects can be seen (Figs 2 and 5). The only difference between the results in these two cases are the extreme values of the expectations. While the shapes of the graphs for var(T(t)), var(I(t)) and var(V(t)) are the same for both normally and triangularly distributed effects (Figs 3 and 6), there is considerable difference in the values of these variances. It can be said that the variances show similar behavior, but with different values. The minimum and maximum values for the variances are as follows: max(var (T(t))) = 6.7858 × 10 7 is obtained at t = 20, while min(var(T(t))) = 0 is obtained at t = 0. max (var(I(t))) = 1.8583 × 10 8 is obtained at t = 5.4, while min(var(I(t))) = 0 is obtained at t = 0. max(var(V(t))) = 1.2817 × 10 7 is obtained at t = 0.6, while min(var(V(t))) = 0 is obtained at t = 20. Hence, it can be said that Fig 6 suggests that the results for the random effects with symmetrical triangular distribution have a smaller variance for T(t), I(t) and V(t).

Comparison of results
More than 10 5 simulations were made for the model under both normally and symmetrically triangularly varying random effects. Results for solution curves, expectations, variances and confidence intervals were given above. However, the results for standard deviation, second moments, third and fourth central moments, skewness and kurtosis were also calculated from the simulations. When the results for the random effect with normal distribution and symmetrical triangular distribution are compared, some differences can be seen. The higher variance in the results for the normally distributed random effects is the first to be noticed ( Table 2).
The results show that there is about 6 times more variance in the normal results compared to the symmetrical triangular results, meaning that the random results for the normal results are 6 times more scattered around the mean values compared to the symmetrically triangularly varying results. The cause of this difference is the characteristics of the distributions used for the random effect. Standard normal distribution has a variance of 1 while the variance of symmetrical triangular distribution in the interval (−1, 1) is 1 6 . The difference in variances causes a difference in standard deviations and also a difference in confidence intervals, as expected. The standard deviation of results for the normally varying random effects are about 2.45 (which is the square-root of 6) times bigger than the standard deviation of the results for the symmetrically triangularly varying effects. As a consequence of the bigger standard deviation, the confidence intervals for the normally distributed effects are larger than the confidence intervals for the symmetrically triangularly distributed effects. Again, these differences can be traced back to the characteristics of the distributions. Further investigation of other characteristics such as higher central moments, skewness and kurtosis yields that while there is not much difference in fourth central moments and kurtosis, there is a noticeable variation Fig 6. Variances of T(t), I(t) and V(t). https://doi.org/10.1371/journal.pone.0181571.g006

Fig 7. Confidence intervals of T(t), I(t) and V(t).
https://doi.org/10.1371/journal.pone.0181571.g007 between the third central moments and kurtosis, which can be interpreted as the consequence of the properties of the distributions used. Finally, it can be said that the symmetrical triangular distribution could be accepted as more favorable for this study since the scattering of the values would mean too much change in the variables of the model, which could affect the stability of the model.

Conclusion
In this study, the deterministic stability and the random behavior of the model from [7] were investigated. The disease-free equilibrium E 0 and the endemic equilibrium point E 1 were found. In addition to the equilibrium points, the spectral radius of the system R 0 was found. Once the spectral radius was calculated, the local stability of the model was examined. The results show that the disease-free equilibrium point of the model is locally asymptotically stable if R 0 1 and unstable if R 0 > 1. It is also shown that the endemic equilibrium E 1 is locally asymptotically stable if R 0 > 1 and unstable otherwise. A Lyapunov function was constructed for the global stability analysis of the disease free equilibrium. The results show that the disease free equilibrium is globally asymptotically stable provided that R 0 < 1. A geometric approach is used for the global stability analysis of the endemic equilibrium. Thus, it is shown that the endemic equilibrium point E 1 is globally asymptotically stable if R 0 > 1. In the last part, the random behavior of the model is examined. Computer simulations are performed for the model with both normally and triangularly varying random effects. Numerical characteristics of the model such as, expectation, variance and confidence intervals are calculated from the simulations and the distributions for the random effect are compared. It is seen that the random behavior of the model is in compliance with the global and local deterministic stability analysis results. The spectral radius calculated with the parameter values obtained from [7] matches the outcomes of the random simulations. The stability analysis performed in this study can be used for a wide number of models in different research areas both on local and global scales. Various other probability distributions such as bilateral exponential (Laplace) and generalized beta distributions could be used for the random effect. Furthermore, random behavior analysis and comparison with the deterministic results can be improved by using Brownian motion to form stochastic differential equations for the models. A stochastic model formed by using stochastic differential equations with Brownian motion could provide better results for the accuracy of the equation system in modeling the real process. The methods of this study could be used on models for other diseases and provide a better understanding of the disease dynamics hence making way for better treatment.