Analysis of hepatitis C viral dynamics using Latin hypercube sampling

Commun Nonlinear Sci Numer Simulat 17 (2012) 5125–5130

Contents lists available at SciVerse ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Analysis of hepatitis C viral dynamics using Latin hypercube sampling Gaurav Pachpute, Siddhartha P. Chakrabarty ⇑ Department of Mathematics, Indian Institute of Technology Guwahati, Guwahati 781 039, Assam, India

a r t i c l e

i n f o

Article history: Received 8 November 2011 Accepted 25 March 2012 Available online 7 April 2012 Keywords: Hepatitis C Latin hypercube sampling Steady states Critical efficacy Statistical tests

a b s t r a c t We consider a mathematical model comprising four coupled ordinary differential equations (ODEs) to study hepatitis C viral dynamics. The model includes the efficacies of a combination therapy of interferon and ribavirin. There are two main objectives of this paper. The first one is to approximate the percentage of cases in which there is a viral clearance in absence of treatment as well as percentage of response to treatment for various efficacy levels. The other is to better understand and identify the parameters that play a key role in the decline of viral load and can be estimated in a clinical setting. A condition for the stability of the uninfected and the infected steady states is presented. A large number of sample points for the model parameters (which are physiologically feasible) are generated using Latin hypercube sampling. An analysis of the simulated values identifies that, approximately 29.85% cases result in clearance of the virus during the early phase of the infection. Results from the v2 and the Spearman’s tests done on the samples, indicate a distinctly different distribution for certain parameters for the cases exhibiting viral clearance under the combination therapy. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Hepatitis C is an infectious disease caused by the hepatitis C virus (HCV) and has a widespread global prevalence. HCV is estimated to have infected 170 million individuals worldwide [1,2]. HCV is usually transmitted through blood contact with an infected person. In places like the United States and Europe, the primary mode of transmission of HCV is through injecting drug use [1]. In India, however, the lack of effective and reliable anti-HCV screening amongst blood donors is a major source of this infection [3,4]. New HCV infections (through intravenous drug usage or otherwise) can be primarily classified as acute and chronic [1]. While the acute cases in general do not have major detrimental implications and the virus is cleared, the ramifications of chronic cases are far reaching [1]. It could lead to liver cirrhosis and may eventually develop into hepatocellular carcinoma (HCC). About 50–80% of HCV infections are chronic in nature [5]. Of these 10–20% develop liver cirrhosis out of which about 5% are likely to be afflicted with HCC [5]. The current treatment for HCV infection involves the combination therapy of pegylated interferon (IFN) and ribavirin [2], which yields long term response in only 50% of the cases. There is very little therapeutic alternative for the cases of nonresponders. Feld and Hoofnagle [6] dwell on the mechanism of action of combination treatment of IFN and ribavirin in HCV infected patients. They report limited success (6–12% for a six-month treatment and 16–20% for a one-year treatment) in case of an IFN based monotreatment. The combination therapy, however, results in significant improvement (more than 50%) in response rates. Perelson et al. [7] in their review article, discuss the role of this combination therapy for HCV infection, taking into account the fall in efficacy levels of the drugs between dosing intervals.

⇑ Corresponding author. Tel.: +91 361 2582606; fax: +91 361 2582649. E-mail addresses: [email protected] (G. Pachpute), [email protected] (S.P. Chakrabarty). 1007-5704/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cnsns.2012.03.035

5126

G. Pachpute, S.P. Chakrabarty / Commun Nonlinear Sci Numer Simulat 17 (2012) 5125–5130

Several mathematical models have been proposed for the study of hepatitis C viral dynamics. One of the earliest models was proposed by Neumann et al. [8], wherein the dynamics of HCV and the effect of interferon-a-2b were studied in vivo. The quantitative representation of the dynamics involved the incorporation of some aspects of the earlier successful models for HIV and HBV [2,8]. The model involves three coupled ODEs, where the key factors are identified as uninfected hepatocytes, productively infected hepatocytes and free HCV virions. The model assumes the growth of the uninfected hepatocytes at a constant rate accompanied by a natural death rate. In addition, HCV is assumed to infect the hepatocytes, thereby creating infected hepatocytes, which has a natural death rate. The infected hepatocytes in turn abets the growth of HCV. A term for the efficacy of IFN is incorporated to reduce the rate of infection of hepatocytes and to block the production of HCV from infected hepatocytes. The model, in presence of therapy, exhibits a decline in the levels of infected hepatocytes and HCV accompanied by an increase in the level of uninfected hepatocytes. A single phasic decline in HCV levels is observed if the efficacy of IFN in blocking the production of virions is taken to be zero, which is inconsistent with clinical observations [2,8]. Otherwise it shows a biphasic decline in viral load which is more realistic from the biomedical point of view. The model indicates that the key role of IFN is in blocking the production of virions from infected cells, as compared to blocking the infection of hepatocytes, which is minimal [2,8]. Dixit et al. [9] have extended the work of Neumann et al. [8] by including the action of ribavirin. In this model [2,9], it is assumed that, ribavirin acts (either on its own or in conjunction with IFN) by rendering a fraction of the newly produced HCV non-infectious, thereby dividing the virion population into infectious and non-infectious virions. Another assumption of the model (in contrast to Neumann et al. [8]) is that, the role of IFN in blocking the production of infected hepatocytes is not significant. The model predicts that, ribavirin does not have an impact on the IFN induced first phase decline [2]. Furthermore, if the efficacy of IFN is large enough, then the second phase decline is not significantly affected by ribavirin. However, if the IFN efficacy is much smaller than 1, then the impact of ribavirin in the second phase decline is more profound. A plausible explanation for this can be attributed to a high IFN efficacy which results in low HCV levels. This, in turn, reduces the role of ribavirin in rendering the virions non-infectious. The model is able to explain why ribavirin enhances the second phase decline in some cases and not in others [2]. Dahari et al. [10] incorporate the proliferation of both infected and uninfected hepatocytes which is not considered in the earlier models. They include density dependent proliferation for both infected and uninfected hepatocytes, which restricts the growth of liver to a maximum possible size. This model is able to explain the limitations of the model of Dixit et al. [9], which cannot account for the non-response and triphasic decline patterns in some cases. The model defines a critical threshold efficacy below which there cannot be sustained virological response [9,10]. If the efficacy is above the threshold one observes a decline in the viral load, which could be biphasic or triphasic. The triphasic decay can be attributed to homeostatic liver regeneration. 2. Mathematical model The HCV model under consideration follows after the earlier models [8–11] and involves the dynamics of uninfected (T) and infected (I) hepatocytes as well as HCV. The model assumes that the uninfected hepatocytes are being produced at a constant density independent rate s and have a natural death rate d. The uninfected hepatocytes are infected by the virions at rate b. The infected hepatocytes thereby produced have a death rate of d. Hepatocytes, both uninfected and infected, proliferate at rate r with T max being the maximum hepatocyte density. In the absence of any kind of treatment, the infected hepatocytes produce infectious virions at a rate p. The administration of IFN lowers this production by a factor of ð1  p Þ, where p is the efficacy of IFN. Ribavirin acts by rendering a fraction q of the newly produced infectious virions (V I ), noninfectious (V NI ), with a clearance rate c for both of these virion populations. The mathematical model under the above assumptions is given by the following system of four coupled ODEs,

  dT T þI  dT  bTV I ; ¼ s þ rT 1  dt T max   dI T þI  dI; ¼ bTV I þ rI 1  dt T max dV I ¼ ð1  qÞð1  p ÞpI  cV I ; dt dV NI ¼ qð1  p ÞpI  cV NI : dt

ð2:1Þ

The model under consideration admits two steady states [10], 1. Uninfected steady state,

T ðuÞ ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi# " T max 4rs ; r  d þ ðr  dÞ2 þ T max 2r

2. Infected steady state,

IðuÞ ¼ 0;

ðuÞ

VI

¼ 0;

ðuÞ

V NI ¼ 0:

G. Pachpute, S.P. Chakrabarty / Commun Nonlinear Sci Numer Simulat 17 (2012) 5125–5130

T

ðiÞ

5127

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi# " 1 4sT max ; D þ D2 þ ¼ 2 rA2

IðiÞ ¼ T ðiÞ ðA  1Þ þ T max  B; ð1  p Þð1  qÞpIðiÞ ; c

V NI ¼

ð1  p Þð1  qÞpbT max ; cr

B¼

ðiÞ

VI ¼

ðiÞ

qV ðiÞ I ; ð1  qÞ

where,

A¼

dT max ; r

D¼

   1 dB 1 T max þ B þ1 : A dA A

Under the physiological conditions of r > d and s 6 dT max it can be shown that the condition for stability of the uninfected steady state is,

ð1  p Þð1  qÞ <

cðdT max þ rT ðuÞ  rT max Þ pbT ðuÞ T max

;

while that of the infected steady state is,

ð1  p Þð1  qÞ >

cðdT max þ rT ðuÞ  rT max Þ pbT ðuÞ T max

:

Thus there is a transcritical bifurcation at,

Cr  ¼

cðdT max þ rT ðuÞ  rT max Þ pbT ðuÞ T max

¼ ð1  p Þð1  qÞ:

3. Latin hypercube sampling In this article, we use the method of Latin hypercube sampling to generate a collection of parameter values from a multivariate normal distribution with a specified mean vector ~ l and covariance matrix R. The theoretical and computational aspects of this method can be found in the pioneering work of McKay et al. [13] which was further advanced by Stein [14]. Latin hypercube sampling is a multidimensional version of one-dimensional stratification [15]. In this method, random variates are generated from a l-dimensional uniform distribution over the hypercube [0, 1)l. Random variates from other ldimensional distributions can be generated from the uniform variates using the inverse transform method [15]. Stratification in one-dimensional uniform distribution can be achieved by dividing the interval into K strata. The extension of this idea to multidimensional stratification is not trivial, as the sample size of K l is computationally expensive even for moderate values of l (since K has to be large enough). We present a brief outline of this method (after Glasserman [15]). ðjÞ To begin with, we generate random variates V i from a uniform distribution over ½ðj  1Þ=K; j=KÞ for i ¼ 1; 2; . . . ; l and a set ðjÞ ðjÞ ðjÞ of independent permutations r1 ; . . . ; rl over the set f1; 2; . . . ; Kg (from K! possibilities). The vectors ðV 1 ; V 2 ; . . . ; V l Þ> (for ðjÞ ri ðjÞ l j ¼ 1; 2; . . . ; K) represent uniformly distributed points over the hypercube ½0; 1Þ . We then set V i V i , which gives us a stratified sample over the hypercube. An obvious way to achieve this is to set, ðjÞ

Vi ¼





ðjÞ ri ðjÞ  1 þ U ðjÞ =K; U i  Uniform½0; 1Þ; i ¼ 1; . . . ; l; j ¼ 1; . . . ; K: i

In generating this construction, one of the crucial assumptions is the independence of marginals, for introducing a correlation between parameter distributions might alter the stratification properties of Latin hypercube sampling. This is very much evident in generating a sample from a normal distribution which does not have a diagonal covariance matrix. The values so generated will not in general be stratified. In order to generate from a normal distribution with mean vector ~ l ¼ ðli Þli¼1 and diagonal covariance matrix R ¼ ðRij Þli;j¼1 we set,

  ðjÞ ðjÞ Z i ¼ U1 li ;Rii V i ;

i ¼ 1; . . . ; l; j ¼ 1; . . . ; K;

where Uli ;Rii is the one-dimensional cumulative normal distribution with mean

li and variance Rii .

4. Results and discussion We implement Latin hypercube sampling procedure described in the previous section using MatLab™ and generate K ¼ 214 sample points for all the parameters. The distribution of the parameters is assumed to be normal and independent of each other (on the lines of [16]). The ranges for the parameter values as given by Dahari et al. [12] are given in Table 1. The mean and standard deviation of the normal distribution used for sampling are chosen so that 95% of the values generated lie

5128

G. Pachpute, S.P. Chakrabarty / Commun Nonlinear Sci Numer Simulat 17 (2012) 5125–5130 Table 1 Range of parameter values used for the Latin hypercube sampling [12]. Parameter

Range of values

Units

s d p b

1—dT max 0.001–0.014 0.1–15.0

cell ml1 day1 day1 virions cell1 day1 ml virions1 day1

c d r T max

108 —106 0.8–22.2 d–3.0 0.1–3.0 0:4  107 —1:3  107

day1 day1 day1 cells ml1

Table 2 Mean and standard deviation of parameter values used for the Latin hypercube sampling. Parameter

Mean

Standard deviation

s

2:005  10

d

7:5  103 7.55

3:3163  103 3.801

5:05  107 11.5 1.507 1.55

2:5255  107 5.4592 0.76173 0.7398

8:5  106

2:2959  106

p b c d r T max

3

1:0202  103

in the parameter range specified, and are given in Table 2. The maximum for the confidence interval of s is taken to be 4  103 . Similarly, the minimum for the range of d is taken to be 0.014. The stratified samples generated from this distribution are tested for positivity as well as the physiological conditions (r > d and s 6 dT max ) as given in Section 2. They were also tested for the d 6 d condition [12] as given in Table 2 and discarded accordingly. These tests left us with about 86% of the original samples. A number of tests are performed to understand how the sample parameter values correlate. v2 -tests are performed on the accepted sample points against the frequency of positive variates from their respective intended distributions proposed earlier. The goodness of fit null hypothesis (set of sample points which are biomedically feasible do not deviate from the intended distribution) is not rejected for any of the sample parameters. On the other hand, the tests reject the independence null hypothesis for all parameters with respect to the significance level chosen (P < 0.05). Since, physiological conditions add ‘‘restrictions’’ to the intended distribution, we can surmise from these tests that the accepted values do not deviate from this distribution. The accepted sample points are categorized into three parts depending on whether the value

Cr  ¼

cðdT max þ rT ðuÞ  rT max Þ pbT ðuÞ T max

is greater than 1, less than 0 or in the interval ½0; 1. It can be proven that the condition Cr  < 0 is equivalent to d > d. Biologically, these cases exhibit a higher death rate of uninfected hepatocytes than that of infected hepatocytes. These are also the cases for which a sustained virological response cannot be achieved for all possible values of p and q. Such values are already eliminated from the samples as a result of the condition d 6 d specified in the first paragraph. Therefore, the parameter samples exclusively contain sample sets that correspond to either 0 6 Cr  6 1 or Cr  > 1. Notice that, when 0 6 Cr  6 1, both steady states can be stable on two disjoint sets of feasible p and q values (due to the transcritical bifurcation). The percentage of such sample points, which emerged from our simulations, is 70.15%. On the other hand, the sample parameter sets with Cr  > 1, can only lead to the stability of the uninfected steady state since, in this case the stability condition for the uninfected steady state is satisfied for all possible values of p and q. The percentage of cases with Cr  > 1 is 29.85%, which is in line with biomedical observations [5,10]. We consider values of p and q in the range [0, 1] and in increments of 0.01. For these 101  101 pairs of values, we keep a count of the sample points which satisfy the condition ð1  p Þð1  qÞ < Cr  . Recall that, this condition corresponds to the stability of the uninfected steady state. Thus, this gives us the cumulative distribution of the percentage of cases, which give a stable uninfected steady state for the corresponding drug efficacies. The results are presented in a contour plot in Fig. 1 and in a surface plot in Fig. 2. We plot the results for the various values of p against q (Fig. 3) and observe that, for smaller values of IFN efficacy p , the response (in terms of percentage of cases with stable uninfected steady state) is very sensitive to change in the efficacy q of ribavirin. The impact of q, however, is much less evident in case of IFN being highly effective. This observation is consistent with the findings of Dixit et al. [9].

G. Pachpute, S.P. Chakrabarty / Commun Nonlinear Sci Numer Simulat 17 (2012) 5125–5130

5129

Contour plot of percentage of cases with stable uninfected steady state 1

100

0.9

90

0.8

80

0.7

70

0.6

ρ

60 0.5 50 0.4 40

0.3

30

0.2

20

0.1 0

0

0.2

0.4

ε

0.6

0.8

1

10

p

Fig. 1. Contour plot of percentage of cases with stable uninfected steady state.

100 80 60 40 20 0 1

1 0.5

0.5

ρ

0

0

εp

Fig. 2. Surface plot of percentage of cases with stable uninfected steady state.

v2 -tests done between these grouped parameters and their original distribution under the independence null hypothesis show that for the samples which have Cr  values greater than 1 or between 0 and 1, the variables p; b; c; d and T max reject the null hypothesis. It implies that, the sample points from these two sets tend to deviate more from their original distribution which leads them to be more restricted in nature. The cases presenting with higher values of d achieve virological response for smaller efficacies as opposed to smaller values of d. Next, we choose a threshold value 0 6 Cr thr 6 1 and analyze the independence of different parameters which have Cr  values on both sides of Cr thr . If Cr  > Cr thr , the uninfected steady state can be reached for smaller values of drug efficacies as compared to Cr  < Cr thr . By varying the value of Cr thr from 0.2 to 0.4, a Spearman’s test shows a modest negative monotonic relationship between p and b for the samples with Cr  > Cr thr and not in the cases with Cr  < Cr thr . Biologically, b represents the rate of infection and p represents the rate with which infected hepatocytes are converted into virions. The results highlight the fact that for the uninfected steady state to be easily reachable, these two rates cannot be simultaneously high. A similar relation is seen between c and d for Cr  < Cr thr , indicating that cases having a high clearance rate of virions (c) as well as a high death rate of infected hepatocytes (d) reach the uninfected steady state for smaller values of efficacies. However, one of the two values being low would make it harder to achieve sustained virological response for low efficacies.

5130

G. Pachpute, S.P. Chakrabarty / Commun Nonlinear Sci Numer Simulat 17 (2012) 5125–5130

Percentage of cases with stable uninfected steady state for various εp 100 90 80 70 60 50

εp=0.25

40

ε =0.75

εp=0.50 p

εp=0.80

30

ε =0.85 p

ε =0.90

20 10

p

εp=0.99

0

0.2

0.4

ρ

0.6

0.8

1

Fig. 3. Percentage of cases with stable uninfected steady state for various ep.

5. Conclusion Latin hypercube sampling enables one to choose appropriate distributions which best fit the observations (such as previous findings based on observed data and physiological constraints of the mathematical model). In this paper, we apply this sampling technique to generate a large number of sample parameter values for an HCV model. Once the sampling is done, we choose values which are positive and satisfy the physiological conditions. A mathematical justification is stated for the death rate of infected hepatocytes being greater than that of the uninfected hepatocytes. Goodness of fit and independence tests can be performed at every stage of adding constraints to ascertain how much and in what way they affect the distribution of the generated samples. In our case, the values found feasible are subjected to the v2 -test and the Spearman’s test. A v2 -test shows that the accepted values do not deviate from their intended distribution. In the cases where viral clearance is achieved (with or without the therapy), the parameters are somewhat restricted. Spearman’s test reveals a monotonic relationship between the rate of infection of hepatocytes and the rate of viral replication for the cases exhibiting viral clearance under the combination therapy. This suggests that, a higher rate of infection needs to be countered by a higher IFN efficacy to lower the rate of viral replication. The methods and tests mentioned in this paper are also performed on a sample generated from the log-normal distribution, which gave similar results. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

Razali K et al. Modelling the hepatitis C virus epidemic in Australia. Drug Alcohol Depen 2007;91:228–35. Dixit NM. Advances in the mathematical modelling of hepatitis C virus dynamics. J Indian Inst Sci 2008;88(1):37–43. Pal SK, Chalamalasetty BS, Choudhuri G. Hepatitis C: a major health problem of India. Curr Sci 2002;83(9):1058–9. Mukhopadhya A. Hepatitis C in India. J Biosci 2008;33(4):465–73. Roe B, Hall WW. Cellular and molecular interactions in coinfection with hepatitis C virus and human immunodeficiency virus. Exp Rev Mol Med 2008;10:1–18. Feld JJ, Hoofnagle JH. Mechanism of action of interferon and ribavirin in treatment of hepatitis C. Nature 2005;436:967–72. Perelson AS, Herrmann E, Micol F, Zeuzem S. New kinetic models for the hepatitis C virus. Hepatology 2005;42(4):749–54. Neumann AU et al. Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-a therapy. Science 1998;282:103–7. Dixit NM, Layden-Almer JE, Layden TJ, Perelson AS. Modelling how ribavirin improves interferon response rates in hepatitis C virus infection. Nature 2004;432:922–4. Dahari H, Lo A, Ribeiro RM, Perelson AS. Modeling hepatitis C virus dynamics: liver regeneration and critical drug efficacy. J Theor Biol 2007;247:371–81. Dahari H, Ribeiro RM, Perelson AS. Triphasic decline of hepatitis C virus RNA during antiviral therapy. Hepatology 2007;46(1):16–21. Dahari H, Shudo E, Cotler SJ, Layden TJ, Perelson AS. Modeling hepatitis C virus kinetics: the relationship between the infected cell loss rate and the final slope of viral decay. Antiviral Ther 2009;14(3):459–64. McKay MD, Beckman RJ, Conover WJ. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 1979;21(2):239–45. Stein M. Large sample properties of simulations using Latin hypercube sampling. Technometrics 1987;29(2):143–51. Glasserman P. Monte Carlo methods in financial engineering. Springer; 2004. Guedj J, Bazzoli C, Neumann AU, Mentre F. Design evaluation and optimization for models of hepatitis C viral dynamics. Stat Med 2011;30:1045–56.