Figures
Abstract
The half-maximal inhibitory concentration IC is an important pharmacodynamic index of drug effectiveness. To estimate this value, the dose response relationship needs to be established, which is generally achieved by fitting monotonic sigmoidal models. However, recent studies on Human Immunodeficiency Virus (HIV) mutants developing resistance to antiviral drugs show that the dose response curve may not be monotonic. Traditional models can fail for nonmonotonic data and ignore observations that may be of biologic significance. Therefore, we propose a nonparametric model to describe the dose response relationship and fit the curve using local polynomial regression. The nonparametric approach is shown to be promising especially for estimating the IC
of some HIV inhibitory drugs, in which there is a dose-dependent stimulation of response for mutant strains. This model strategy may be applicable to general pharmacologic, toxicologic, or other biomedical data that exhibits a nonmonotonic dose response relationship for which traditional parametric models fail.
Citation: Zhang H, Holden-Wiltse J, Wang J, Liang H (2013) A Strategy to Model Nonmonotonic Dose-Response Curve and Estimate IC50. PLoS ONE 8(8): e69301. https://doi.org/10.1371/journal.pone.0069301
Editor: Shan-Lu Liu, University of Missouri, United States of America
Received: February 21, 2013; Accepted: June 6, 2013; Published: August 1, 2013
Copyright: © 2013 Zhang et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: The research of Holden-Wiltse and Liang was partially supported by the National Institute of Allergy and Infectious Diseases/National Institutes of Health grants HHSN272201000055C. Liang's research was also partially supported by National Science Foundation grants DMS-1007167 and DMS-1207444. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Motivation
Drugs that inhibit the reverse transcriptase (RT) activity of Human Immunodeficiency Virus (HIV) are widely used to treat HIV infection. RT is an ideal target for antiviral HIV therapy because it is the key required for HIV replication. Non-nucleoside reverse transcriptase inhibitors (NNRTIs) inhibit RT activity by selectively binding RT at a hydrophobic binding pocket adjacent to the polymerase active site. Efavirenz (EFV) is an commonly used NNRTI to treat HIV infection [1]–[3] but patients can develop resistance to this drug because of the development of mutations in the NNRTI binding cite which in turn inhibits NNRTI binding [4] and can lead to resistance mutations such as K101E, K103N, Y188C, G190S, G190A [5], [6], and L100I [7].
Understanding the pharmacodynamic properties associated with the development of NNRTI resistant mutations is vital for devising treatment strategies for HIV. In pharmacodynamics, the drug-target interaction can be modeled by:where
denotes the drug and
the target (usually enzyme). Drug efficiency is primarily determined by the drug target binding affinity. In pharmacodynamic studies, the drug target affinity is usually assessed by comparing dose response curves: the stronger the drug binds target, the steeper the curve is.
Therefore, a critical index of the dose response curve, the half-maximal inhibitory concentration (IC), is commonly used to compare the binding affinities of drugs to the same target. IC
represents the concentration of a drug that is required for 50% of maximal inhibition in vitro. IC
and IC
are the concentrations corresponding to 25% and 75% inhibition, respectively. The dose response curve usually has the steepest portion in the middle. Thus, using IC
, rather than IC
or IC
, minimizes the random error for estimation, making IC
the preferred measure of drug affinity.
To estimate the IC value, several nonlinear functions have been commonly used, for example,
(1)where
is the drug concentration,
is the percentage of inhibition at this concentration, and
is a shape parameter. Other parametric models include the complementary log-log model for asymmetric quantal response data, and the two-parameter Weibull model for carcinogenic experiments [8].
An important feature of the function (1) is that as the value of increases from
to infinity,
increases from
to
, reflecting that a drug's inhibitory potential changes from none to full inhibition as the drug concentration increases (Fig. 1A). As per (1), a steeper dose response curve corresponds to a smaller IC
; for a given IC
, the curve shape is depicted by
(Fig. 1B). In addition, this function curve has the steepest portion in the middle, which is a characteristic of a sigmoidal dose response relationship.
The advantages of this function are that (a) it is symmetrically about the IC; (b) it is monotonic, which equals that the target protein has an inhibitory binding site only, and the drug has an inhibitory effect only; and (c) IC
and
can be easily estimated under certain conditions. On the other hand, these advantages become restrictive when some conditions fail, for instance, if observations are not monotonic. As a consequence, the models can produce remarkably biased estimates or not even fit the observations. For example, Bliss's beetle data show that symmetry is not a required feature of a dose response curve [9]. Another example is that a viral mutation occurs before the drug concentration reaches a certain level, such as in the following example.
A recent study on HIV mutations conferring resistance to NNRTI found that the monotonicity relation does not always hold [10]. The dose response was determined as proportion reduction in HIV replication at a given NNRTI dose relative to viral replication in the absence of drug. As seen in Fig. 2, the study shows that replication of HIV mutation M230L was promoted when the concentration of EFV is lower than 70 nM. Similarly, our data example shows that increasing the concentration of EFV stimulates the replication of HIV K101E+G190S mutant strain when EFV concentrations are below 2000 nM (Fig. 2). It has also been reported that EFV stimulation of the K101E+G190S double mutant strain can be abolished by the presence of additional M41L+T215Y mutation [11]. A potential explanation for this nonsigmoidal dose response relationship is that dimerization is essential for a fully functional RT. For the double mutant K101E+G190S strain, at low concentrations EFV can enhance the dimerization of the two subunits of RT without interfering with the binding of the incoming nucleotide during DNA polymerization.
The HIV strain names are the same as in previous publication [11]. Dose response was determined as proportion reduction in HIV replication at a given EFV dose relative to HIV replication in the absence of EFV. Viral replication above 100% indicates that suboptimal doses of EFV potentiate the ability of the viral strain to replicate compared to the absence of EFV. The data used to generate figures is available upon request.
As the real dose response relationship is nonmonotonic in our example, our preliminary analysis indicates that traditional estimation methods for the sigmoidal model (1) fitting lead to results that either do not converge or are remarkably biased, which may reach an erroneous conclusion. Thus, appropriate estimation of IC for this type of dose response relationship poses statistical challenges. If we use model (1) to fit data when the dose response pattern is nonmonotonic, the fit is poor, and the estimated IC
values are not reliable because the fact that lower EFV concentrations can enhance replication of an HIV mutant strain is ignored.
Thus, to appropriately estimate the pattern of observations and then estimate IC, we developed a robust modeling strategy to test whether: (i) our model fitting is comparable to monotonic parametric models such as model (1) when the observed data are monotonic; and (ii) our model fitting yields reasonable estimates when the data pattern is nonmonotonic and monotonic parametric models, such as model (1), does not work.
The rest of this paper is organized as follows. Section 2 briefly introduces monotonicity testing, our model, estimation, and test methods. Section 3 gives simulation results, including p-values of the monotonicity test. Section 4 presents extensive analysis of our real data example, including estimated ICs using the proposed model and model (1) when it is appropriate.
Methods
We propose that the inhibition percentage and concentration
are related in the form
(2)where
is the measurement error with mean zero and finite variance,
is a mathematical function, but no restrictions are applied on the form of
(i.e., no specified
as quadratic, parametric, or increasing in
, etc.) Since the structure of the model is not fixed, it is called a nonparametric method. Hence, we perform an empirical analysis of the data to estimate
. We use the observations
to estimate
, denoted as
, by appropriate statistical techniques. The IC
may be estimated as the point
that satisfies
. As noted above, we first need to determine whether the function
is monotonic.
Subsection 1 Monotonicity test
Testing the monotonicity of a dose response relationship is of practical interest and has been studies previously. Several parametric and nonparametric methods have been proposed in the statistical literature. For example, Ramsay [12] studied the use of monotone splines to model a dose response function. Bowman et al. [13] developed a monotonicity test by using local linear estimation of the curve, followed by a critical bandwidth test. Hall and Heckman [14] proposed an alternative approach that focuses on “ running gradient” estimation over very short intervals. The method of Hall and Heckman is more effective in estimating the flat part of the curve and is also more sensitive to small dips in the curve.
For our study, we adopt the method of Hall and Heckman, which is based on the following principle. Let be integers and
be constants. For each pair of
, define the estimators of
and
by
and
, as the arguments of the following objective function:
Let , where
satisfies
. Hall and Heckman [14] stated that a large
indicates that the null hypothesis,
being monotonic, should be rejected. To obtain the
on the basis of
, they suggested the following procedure. First, data should be fit with the nonparametric model
. An estimation of
should be obtained by a consistent estimator of
, such as the local linear estimator. Assuming that a constant function is the most difficult nondecreasing form to be tested, Hall and Heckman used
to obtain the
. Specifically, using the estimated
, they resampled and obtained a new dataset
, by which they obtained
. Repeated sampling n times resulted in a set
of size
. By taking the
th ordered
s as the critical value for
, that is, when
obtained from the real data is greater than this
, we claim that the function is nonmonotonic at the
level.
Subsection 2 Nonparametric fitting
We used local linear regression [15] to fit the dose response curve. Assuming that has bounded, continuous second partial derivatives, by Taylor expansion,
, in a neighborhood of
, can be approximated as:
The estimator of
at
is the solution of
by minimizing
subject to (
), where
is a bandwidth controlling the size of the local neighborhood,
, with
being a kernel function assigning weights to each data point.
Simulation
To investigate how the Hall and Heckman [14] test performed for small and moderate sample sizes, we conducted a simulation study. Let and
or
,
for
. Consider 2 cases: (a)
; and (b)
, where
.
Case (a) reflects that and
have a monotonic relationship, whereas case (b) indicates that the monotonicity is violated. Case (b) is based on the relationship between the inhibition of HIV mutant K101E+G190S strain and EFV concentration from our real data example, and is meant to show the performance power of the test. Figure 3 depicts the patterns of
and
, with
showing a pronounced dip around
.
The error follows a normal distribution
, with
. Our data are generated from the model (2) with
or
, and
. For each case, we considered 3 sample sizes
, and generated
independent datasets for each of the 3 error variances.
-value was determined as the probability that
is greater than
. Table 1 shows the
-values of different simulation cases. As
is a nondecreasing function, the
-values should be greater than
, which was set at
; in contrast,
-values for
should be lower than
. All test results based on
gave
-values greater than
, although there was a slight trend of decreasing
-values as
increases. Thus, on the basis of
simulations, the monotonicity test method did not reject the null hypothesis at
and correctly concluded that the relationship is monotonic. In contrast, when
, all test results based on
showed
-values lower than
, indicating that the null hypothesis (monotonicity) would be correctly rejected even when the sample size is very small. When
was increased to
and the sample size was as small as 20, smaller
values, i.e. higher dip sizes (Fig. 3), still gave
-values lower than
. However, greater
values showed
-values slightly higher than 0.05. When the sample size was increased to 50, all
-values were lower than 0.05. When the sample size was further increased to 100, all
-values were lower than
. These results indicate that even with a high noise level, the monotonicity test is still reliable, particularly for large sample sizes. These results show that the monotonicity test is reliable and robust.
Real Data Analysis
We performed a monotonicity test for all HIV mutation dose response curves for percent viral replication compared to no drug (Fig. 2). For each mutant strain, we repeated the monotonicity test 100 times and have reported the average -value in Table 2.
The null hypothesis of monotonicity in HIV mutants K101E+G190S, M230L, and K101E+G190S(D10) was rejected with (Table 2). This result is consistent with the observed shape of the dose response curves (Fig. 2). Surprisingly, the test failed to reject
the monotonicity hypothesis for M41L+K101E+G190S+T215Y(D10) data, probably because the local polynomial fitting of this dataset still gives a non decreasing curve.
We then used the traditional sigmoidal model to fit the data and found that this method did not converge for HIV mutants K101E+G190S, M230L, and K101E+G190S(D10) because of the nonmonotonicity while the model did converge for datasets G190S, K101E, L74V+K101E+G190S, M41L+K101E+G190S+T215Y, M41L+G190S+.
T215Y, M41L+K101E+G190S+T215Y(D10), and M41L+G190S+T215Y(D10). We then fit all datasets again, using the local polynomial regression method (Section 2). Figure 4 shows the fitted curves for datasets K101E+G190S, L74V+K101E+G190S, M230L, and M41L+K101E+G190S+.
Solid line indicates the nonparametric model and dashed line indicates the sigmoidal model if it is available.
T215Y. Fitted curves for datasets L74V+K101E+G190S and M41L+K101E+G190S+T215Y showed that the parametric and nonparametric methods gave comparable results.
As shown in Fig. 4, lower concentrations of EFV clearly stimulate the replication of HIV K101E+G190S and M230L compared to no EFV or at high EFV concentrations. However, this property cannot be recognized by using the traditional sigmoidal model fitting. To test the efficiency of our proposed nonparametric method when the monotonicity property is satisfied, we applied our method to the L74V+K101E+G190S and M41L+K101E+G190S+T215Y datasets. The dose response curves obtained by using the nonparametric model are very similar to those by the sigmoidal model, which confirms the efficiency of the nonparametric method. Table 3 compares the estimated IC values for all datasets obtained by using both methods. When the sigmoidal model works well, as for HIV mutant strains G190S, K101E, L74V+K101E+G190S, M41L+K101E+G190S+T215Y, M41L+G190S+T215Y, M41L+K101E+G190S.
+T215Y(D10), and M41L+G190S+T215Y(D10) (these datasets also satisfy the monotonicity property, as shown in Table 2), the two estimated IC values for the same dataset are close (Table 3). In contrast, because of the lack of monotonicity, the sigmoidal model fails to fit the curves for HIV strains of K101E+G190S, M230L and K101E+G190S(D10) (Table 2). For these datasets, the nonparametric model becomes a better alternative for IC
estimations (Table 3).
Discussion
When the dose response relationship and associated parameters such as IC are studied, data observations, which make the pattern nonmonotonic, are generally deleted in order to use the model (1) or similar monotonic functions. However, by deleting these “unusual” observations, some important information may be lost. For example, in Fig. 2, the observation that a lower EFV concentration stimulates HIV K101E+G190S replication can be neglected if these data points are deleted. Removing the unusual data points leaves only 2 observations in this dataset, making the fitting procedure impossible.
In this paper, we proposed a nonparametric approach as an alternative to the parametric sigmoidal model to fit the dose response curve and estimate IC, and suggested a monotonic check of the dose response relationship at the first stage. If the monotonicity is satisfied, either the traditional sigmoidal model fitting or our local polynomial regression fitting can be applied. If monotonicity is not satisfied, our model is more suited to estimate the IC
. Using this new approach, important dose response features will not be omitted. A similar approach has been used to quantify protein lysate assays [16], although no monotonicity needs to be satisfied in that case.
Our proposed method can also be used for other dose response modeling scenarios, such as hormesis dose response curves. In toxicology, hormesis is a special dose response feature characterized by low dose stimulation and high dose inhibition [17]–[19], giving a J-shape dose response curve. Our nonparametric model is more suited than traditional monotonic models to fit this J-shaped curve. Also, our method can be used to model other nonparametric curves such as U-shaped dose response relationships frequently observed in toxicology and epidemiology studies [20].
The human trefoil peptide (TFF1), a small cysteine-rich secreted protein, stimulates cell migration by chemotaxis. The dose response curve of TFF1 inducing breast cancer cell movement shows a clear bell shape [21]. Similar curves are also seen in may other biomedical studies [22]–[27]. Our nonparametric model may be more suitable than sigmoidal models to fit such dose response curves and estimate the IC.
Our approach for estimating IC can also be used to estimate the half-maximal effective concentration, which is commonly used when the drug enhances its target's activity, and the lethal dose 50%, or the lethal concentration and time of a toxic substance or radiation represents the dose needed to kill half the tested population. Since the results we obtained are based on large sample theory, a potential limitation of our proposed method is that a moderate sample size may be needed, although a minimum sample size is not determined.
Acknowledgments
The authors are grateful to Dr. Heckman for sharing her codes, and the Editor and two referees for their valuable comments and suggestions.
Author Contributions
Conceived and designed the experiments: HZ HL. Performed the experiments: HZ JW HL. Analyzed the data: HZ HL. Contributed reagents/materials/analysis tools: HZ HL. Wrote the paper: HZ JHW HL.
References
- 1. Staszewski S, Morales-Ramirez J, Tashima KT, Rachlis A, Skiest D, et al. (1999) Efavirenz plus zidovudine and lamivudine, efavirenz plus indinavir, and indinavir plus zidovudine and lamivudine in the treatment of HIV-1 infection in adults. New England Journal of Medicine 341: 1865–1873.
- 2. Robbins GK, De Gruttola V, Shafer RW, Smeaton LM, Snyder SW, et al. (2003) Comparison of sequential three-drug regimens as initial therapy for HIV-1 infection. New England Journal of Medicine 349: 2293–2303.
- 3. Gulick RM, Ribaudo HJ, Shikuma CM, Lustgarten S, Squires KE, et al. (2004) Triple-nucleoside regimens versus efavirenz-containing regimens for the initial treatment of HIV-1 infection. New England Journal of Medicine 350: 1850–1861.
- 4. Domaoal RA, Bambara RA, Demeter LM (2006) HIV-1 reverse transcriptase mutants resistant to nonnucleoside reverse transcriptase inhibitors do not adversely affect DNA synthesis: pre-steady-state and steady-state kinetic studies. Journal of Acquired Immune Deficiency Syndromes 42: 405–411.
- 5. Bacheler L, Jeffrey S, Hanna G, D'Aquila R, Wallace L, et al. (2001) Genotypic correlates of phenotypic resistance to efavirenz in virus isolates from patients failing nonnucleoside reverse tran-scriptase inhibitor therapy. Journal of Virology 75: 4999–5008.
- 6. Bacheler LT, Anton ED, Kudish P, Baker D, Bunville J, et al. (2000) Human immunodeficiency virus type 1 mutations selected in patients failing efavirenz combination therapy. Antimicrobial Agents and Chemotherapy 44: 2475–2484.
- 7. Koval CE, Dykes C, Wang J, Demeter LM (2006) Relative replication fitness of efavirenz-resistant mutants of HIV-1: correlation with frequency during clinical therapy and evidence of compensation for the reduced fitness of K103N+ L100I by the nucleoside resistance mutation L74V. Virology 353: 184–192.
- 8. Prentice RL (1976) A generalization of the probit and logit methods for dose response curves. Biometrics 32: 761–768.
- 9. Bliss C (1941) Statistical problems in estimating populations of Japanese beetle larvae. Journal of Economic Entomology 34: 221–232.
- 10. Huang W, Parkin N, Lie Y, Wrin T, Haubrich R, et al. (2000) A novel HIV-1 RT mutation (M230L) confers NNRTI resistance and dose-dependent stimulation of replication. Antiviral Therapy 5: 24–25.
- 11. Wang J, Liang H, Bacheler L, Wu H, Deriziotis K, et al. (2010) The non-nucleoside reverse transcrip-tase inhibitor efavirenz stimulates replication of human immunodeficiency virus type 1 harboring certain non-nucleoside resistance mutations. Virology 402: 228–237.
- 12. Ramsay J (1988) Monotone regression splines in action. Statistical Science 3: 425–441.
- 13. Bowman A, Jones M, Gijbels I (1998) Testing monotonicity of regression. Journal of Computational and Graphical Statistics 7: 489–500.
- 14. Hall P, Heckman NE (2000) Testing for monotonicity of a regression mean by calibrating for linear functions. Annals of Statistics 28: 20–39.
- 15.
Fan J, Gijbels I (1996) Local Polynomial Modelling And Its Applications, volume 66. Chapman & Hall/CRC.
- 16. Hu J, He X, Baggerly KA, Coombes KR, Hennessy BT, et al. (2007) Non-parametric quanti_cation of protein lysate arrays. Bioinformatics 23: 1986–1994.
- 17. Calabrese EJ, Baldwin LA (1998) Hormesis as a biological hypothesis. Environmental Health Perspectives 106: 357–362.
- 18. Sagan LA (1987) What is hormesis and why haven't we heard about it before. Health Physics 52: 521–525.
- 19. Stebbing A (1982) Hormesisthe stimulation of growth by low levels of inhibitors. Science of the Total Environment 22: 213–234.
- 20. Davis JM, Svendsgaard DJ (1990) U-shaped dose-response curves: Their occurrence and implica-tions for risk assessment. Journal of Toxicology and Environmental Health, Part A 30: 71–83.
- 21. Prest SJ, May FE, Westley BR (2002) The estrogen-regulated protein, TFF1, stimulates migration of human breast cancer cells. The FASEB Journal 16: 592–594.
- 22. Visnjić D, Batinić D, Banfić H (1999) Di_erent roles of protein kinase c alpha and delta isoforms in the regulation of neutral sphingomyelinase activity in HL-60 cells. Biochemical Journal 344: 921–928.
- 23.
Poland GA, Jacobson RM, Koutsky LA, Tamms GM, Railkar R, et al.. (2005) Immunogenicity and reactogenicity of a novel vaccine for human papillomavirus 16: a 2-year randomized controlled clinical trial. In: Mayo Clinic Proceedings. Elsevier, volume 80, 601–610.
- 24. Masuda K, Yokomizo T, Izumi T, Shimizu T (1999) cDNA cloning and characterization of guinea-pig leukotriene B4 receptor. Biochemical Journal 342: 79–85.
- 25. Kohn MC (1995) Biochemical mechanisms and cancer risk assessment models for dioxin. Toxicology 102: 133–138.
- 26. Gibbs B, Rathling A, Zillikens D, Huber M, Haas H (2006) Initial Fc epsilon RI-mediated signal strength plays a key role in regulating basophil signaling and deactivation. Journal of Allergy and Clinical Immunology 118: 1060–1067.
- 27. Gallicchio M, Rosa AC, Benetti E, Collino M, Dianzani C, et al. (2006) Substance p-induced cyclooxygenase-2 expression in human umbilical vein endothelial cells. British Journal of Pharma-cology 147: 681–689.