
==== Front
J Med Chem
J Med Chem
jm
jmcmar
Journal of Medicinal Chemistry
0022-2623
1520-4804
American Chemical Society

37697621
10.1021/acs.jmedchem.3c00107
Article
Inclusion of Control Data in Fits to Concentration–Response Curves Improves Estimates of Half-Maximal Concentrations
https://orcid.org/0000-0001-5302-3271
La Van Ngoc Thuy †
Nicholson Stanley ‡
Haneef Amna †
https://orcid.org/0000-0002-6000-3436
Kang Lulu ‡
https://orcid.org/0000-0002-4802-2618
Minh David D. L. *§
† Department of Biology, Illinois Institute of Technology, Chicago, Illinois 60616, United States
‡ Department of Applied Mathematics, Illinois Institute of Technology, Chicago, Illinois 60616, United States
§ Department of Chemistry, Illinois Institute of Technology, Chicago, Illinois 60616, United States
* Email: dminh@iit.edu.
12 09 2023
28 09 2023
12 09 2024
66 18 1275112761
18 01 2023
© 2023 The Authors. Published by American Chemical Society
2023
The Authors
https://creativecommons.org/licenses/by/4.0/ Permits the broadest form of re-use including for commercial purposes, provided that author attribution and integrity are maintained (https://creativecommons.org/licenses/by/4.0/).

Concentration–response curves, in which the effect of varying the concentration on the response of an assay is measured, are widely used to evaluate biological effects of chemical compounds. While National Center for Advancing Translational Sciences guidelines specify that readouts should be normalized by the controls, recommended statistical analyses do not explicitly fit to the control data. Here, we introduce a nonlinear regression procedure based on maximum likelihood estimation that determines parameters for the classical Hill equation by fitting the model to both the curve and the control data. Simulations show that the proposed procedure provides more precise parameters compared with previously prescribed practices. Analysis of enzymatic inhibition data from the COVID Moonshot demonstrates that the proposed procedure yields a lower asymptotic standard error for estimated parameters. Benefits are most evident in the analysis of the incomplete curves. We also find that Lenth’s outlier detection method appears to determine parameters more precisely.

National Science Foundation 10.13039/100000001 1905324 Pritzker Institute of Biomedical Science and Engineering 10.13039/100012795 NA document-id-old-9jm3c00107
document-id-new-14jm3c00107
ccc-price
==== Body
pmcIntroduction

Some data are just underappreciated. Maybe they look different or come from a different background from most other data. Maybe they do not fit neatly into common notions of what data on a “curve” should look like. Whatever the case, they are pigeonholed into a restricted role that limits their contributions. But if people would only give them the opportunity, they could improve the entire fit.

This is the case for control data in concentration–response curves (CRCs). CRCs describe the response of a biological assay to different concentrations of a chemical compound. When the assay measures the response of an organism, the data are termed a dose–response curve. CRCs are an important step in evaluating the biological effects of a chemical compound. Widely used in pharmacology and toxicology, CRCs have become even more ubiquitous with the emergence of high-throughput screening.

CRC data are often summarized by the concentration required to induce a response halfway between the minimum and maximum. For inhibition assays, the IC50 specifies the half-maximal inhibitory concentration. For an agonist/stimulator assay, the EC50 is the concentration of substance that elicits half of the maximal response. IC50/EC50 values are often used to classify the potency of a compound against a target. For example, drugs may be categorized as high potency (below 1 μ M), medium potency (between 1 and 10 μ M), or low potency (above 10 μ M).1

While it is possible to analyze CRC data and determine IC50/EC50 values in many ways, standard approaches that are considered best practices in the scientific community are codified in the Assay Guidance Manual (AGM).2 The AGM is a free online resource written by over 100 authors and managed by the National Center for Advancing Translational Sciences with input from experts working in industry, academia, and government. Its scope includes “appropriate statistical ways to analyze assay results and accommodate minor changes to assay protocols to ensure robustness.”3

According to the AGM,2 CRC data should first be normalized by control data. As described in the section “Data Standardization for Results Management”, the specific type of control data depends on the type of assay. In enzymatic inhibition assays, negative controls contain the enzyme including cofactors and substrate(s) but no putative inhibitor. Positive controls also include substrate(s) and cofactors, but the enzyme is either absent or fully inhibited. Regardless of the assay, normalization is performed by subtracting the response by the minimum, dividing by the range (difference between the minimum and maximum), and multiplying by 100 to obtain a percentage. For most assays, the minimum is the negative control and maximum the positive control. (The exceptions are for in vitro functional assays of agonist and inverse agonists, in which the AGM recommends fitting to the CRC of a reference agonist to obtain the minimum and maximum.)

After normalization, the model most commonly fit to the CRC data is the classical Hill equation. As CRC data are usually monotonic with a sigmoid shape on a semilog plot, most curves are reasonably modeled by a logistic function in which the dependent variable x is the logarithm of the concentration c, x = log10(c/cθ), where cθ = 1 M is the standard concentration. The classical Hill equation scales the logistic function to interpolate between the minimum response at the bottom of the curve Rb and maximum response at the top of the curve Rt as1

where x50 is the base 10 logarithm of the concentration that induces the half-maximal response, H is the Hill slope, and θ = {Rt,Rb,x50,H} denotes the full set of parameters for the equation. The x50 corresponds to the base 10 logarithm of IC50 or EC50, depending on the type of assay.

While the control data are used for normalization, statistical analyses recommended by the AGM do not explicitly include control data in curve fitting. The AGM goes into considerable detail about fitting the classical Hill equation to the CRC data. Fitting to competitive binding and functional assay data is discussed in the section “Assay Operations for SAR Support”. Fitting to in vivo assay data is described in section 3.5 of the chapter on “In Vivo Assay Guidelines”. In both chapters, the AGM suggests fitting the four-parameter logistic (4PL) model, in which all of θ is allowed to vary, to most CRC data. The control data are not included in the fit and therefore have no effect on the estimated x50 or H. For some data, the AGM recommends the three-parameter logistic fixed bottom (3PLFB) model, where Rb is fixed at zero, or the three-parameter logistic fixed top (3PLFT) model, where Rt is fixed at 100. Conceptually, three-parameter models make the most sense when control data can be reasonably construed as measurements of the minimal response, Rb or maximal response Rt. For example, in an enzymatic inhibition assay, if the response is proportional to the concentration of a product, then the negative control without an inhibitor is a measurement of Rt. Fixing Rb or Rt can strongly influence the estimated values of both x50 and H, especially for incomplete curves.4 While three-parameter fits do not explicitly incorporate control data, these procedures actually give the control data a privileged position; these fits essentially set the value of a parameter based on the control as opposed to multiple concentrations along the curve.

Explicitly including control data in curve fitting, neither unduly ignoring nor privileging them, is a simple and intuitive idea. However, we found only two papers that evaluated whether it is a good idea. Weimer et al.5 compared various fitting procedures to androgen receptor binding assay data and to simulated data with similar parameters. While negative controls were included in all fits, positive controls were either included or excluded. They found that including positive controls in fits to experimental data most significantly changes the estimate of Rb. In fits to simulated curves, including positive control data improved the precision of the x50 estimates (Table 2 of their paper). Thus, they recommended fits including both positive and negative controls. Kappenberg et al.6 considered the problem of deviating controls - when the response of the negative control deviates from the response at low concentrations. In a literature review, they found deviating controls to be a common occurrence in the toxicological literature. They found that if deviations are small, the influence of fitting with or without controls depends on whether the curves are complete, with response values near both Rt and Rb. For complete curves, inclusion of controls had little influence on the fraction of estimates deemed acceptable. On the other hand, including control data was clearly beneficial for the analysis of incomplete curves.

Curve fitting including control data is only rarely mentioned in software documentation and the affiliated literature. Among commercial packages, the documentation of Origin7 does not mention control data whereas GraphPad documentation8 mentions a suggestion from Weimer et al.:5 Even if software is unable to accept concentrations of zero and infinity, control data may be incorporated by assigning them extremely high (for positive control data) or low (for negative control data) concentrations. The oft-cited book “Fitting models to biological data using linear and nonlinear regression: A practical guide to curve fitting”,9 coauthored by the founder and chief executive officer of GraphPad, spans 351 pages. As with the AGM, the use of control data is to fix curve parameters. It also describes how CRCs from control compounds may be used for global fitting and to evaluate whether parameters for a compound are different from those of the control compound. However, it does not discuss the inclusion of control conditions in curve fits. Beyond commercial software, the past decade has seen the publication of a number of papers describing software tools for fitting the Hill equation and other models to CRC data.10−15 These papers have focused on how the tools make curve fitting freely available and easier to use,11,12,15 fit data with more complex multiphasic models,10,14 and perform Bayesian uncertainty quantification.13 Other papers since 2010 have described how to make curve fitting more robust via an evolutionary algorithm,16 a systematic search of parameter space,17 or outlier detection.18 In only a few of these papers is fitting to control data mentioned at all. One of these mentions is in the paper about the popular drc package12 within the statistical computing environment R,19 which as of April 17, 2023 has received 1555 citations according to the Web of Science. It explains that model functions are implemented in a way such that they are defined at a concentration of zero. eeFit also formulates the Hill equation to avoid a singularity at zero concentration.15 However, these papers do not explain the benefits of including control data in curve fits.

Here, we describe and evaluate a method for incorporating control data into CRC fits. Based on the assumption that measurement error is Gaussian distributed, our method is based on maximizing the likelihood of observing all of the data given the full set of parameters θ. We also assume homoscedasticity, that is, that all the data on the curve have the same variance, but the approach may be readily generalized to heteroscedastic data. Given that they may contain a different number of species and require a different number of dilution operations, controls are assumed to have a variance different from the curve. Without control data, the maximum likelihood estimate of model parameters is the standard procedure of nonlinear least-squares regression; it minimizes the sum of square deviation between the model and the observed data. On the other hand, when control data are included, the estimated model parameters are distinct (Figure 1) and, as we show, usually more accurate. We demonstrate the method on simulations of complete and incomplete CRCs. We also apply the proposed procedure to 1830 CRCs for two types of SARS-CoV-2 main protease enzymatic inhibition assays that were collected for the COVID Moonshot, an open-source effort to develop a COVID-19 antiviral.20 Lastly, we describe the results of an outlier detection method.

Figure 1 Example fits to CRC data with (blue line) or without (orange line) controls. Curve data are shown with black circles and control data as green triangles. The x position of the control data is not meaningful; data on the left correspond to x → –∞ and data on the right correspond to x → ∞. The vertical lines indicate the true x50 (black dotted) or the estimated x50 from CRC fitting with (blue dashed) or without (orange dashed) controls.

Results and Discussion

As detailed in the Experimental Section, curve fitting was performed with five methods: our proposed statistical approach and up to four established procedures. Our proposed approach may be summarized as a four-parameter fitting to a logistic function with controls (4PL+C). The established procedures are a four-parameter fit without controls (4PL); a three-parameter fit with fixed bottom such that Rb = 0 (3PLFB); a three-parameter fit with a fixed top such that Rt = 100 (3PLFT); and two parameter fit with a fixed bottom and fixed top such that Rb = 0 and Rt = 100 (2PL). Simulations were analyzed with all five procedures, while biochemical assay data were analyzed with the 4PL+C, 4PL, and 3PLFB procedures. In analyses with the four established procedures, data were normalized based on the sample mean of the controls. With the proposed approach, data were normalized based on and obtained from the fit.

The 4PL+C Procedure Improves the Accuracy of All Parameter Estimates Based on Simulated Data

The proposed 4PL+C procedure results in parameter estimates that have a lower error than all other tested statistical analysis procedures. For all tested procedures, the distribution of parameter estimates appears to be Gaussian, symmetric about the true value (Figures S1 and S2 of Supporting Information 1). For a set of N estimates of a parameter θ, with n ∈ {1,2,3, ...,N}, the mean square error (MSE) is . This metric accounts for both the variance and the bias in an estimate. However, as mean values of all parameters are indistinguishable from the true values, none of the estimators are biased, and distinctions between are due to variance. The 4PL+C procedure leads to a reduction in MSE that is especially pronounced for Rb (∼0.17 of 4PL) and Rt (∼0.48 of 4PL) and smaller but still statistically significant for x50 and H (∼0.75 of 4PL for both) (Table 1). The change in MSE is larger than the standard deviation across 5 independent sets of curves.

Table 1 Comparison of the Error in Parameter Estimates from Simulated CRCsa

MSE	4PL+C	4PL	3PLFB	3PLFT	2PL	
Low variance	
Rb	3.8 × 10–1 (1.1 × 10–2)	2.3 (4.6 × 10–2)	 	2.3 (4.6 × 10–2)	 	
Rt	6.2 × 10–1 (1.3 × 10–2)	1.3 (2.7 × 10–2)	1.3 (2.5 × 10–2)	 	 	
x50	4.4 × 10–4 (9.1 × 10–6)	5.8 × 10–4 (8.7 × 10–6)	5.0 × 10–4 (1.3 × 10–5)	6.1 × 10–4 (1.0 × 10–5)	6.0 × 10–4 (1.2 × 10–5)	
H	1.8 × 10–3 (2.1 × 10–5)	2.4 × 10–3 (2.7 × 10–5)	2.1 × 10–3 (3.1 × 10–5)	2.8 × 10–3 (3.9 × 10–5)	2.2 × 10–3 (2.5 × 10–5)	
High variance	
Rb	1.1 (3.5 × 10–2)	6.9 (1.4 × 10–1)	 	7.1 (1.4 × 10–1)	 	
Rt	1.9 (4.0 × 10–2)	3.9 (8.2 × 10–2)	3.8 (7.5 × 10–2)	 	 	
x50	1.3 × 10–3 (2.7 × 10–5)	1.7 × 10–3 (2.7 × 10–5)	1.5 × 10–3 (4.0 × 10–5)	1.8 × 10–3 (3.2 × 10–5)	1.8 × 10–3 (3.9 × 10–5)	
H	5.6 × 10–3 (7.8 × 10–5)	7.5 × 10–3 (9.6 × 10–5)	6.4 × 10–3 (1.0 × 10–4)	8.8 × 10–3 (1.3 × 10–4)	6.9 × 10–3 (7.5 × 10–5)	
a Mean and standard deviation (in parentheses) of the MSE of parameter estimates from 5000 simulations across 10 independent sets. Note that the parameters Rb and Rt correspond to normalized data, opposed to parameters Rb° and Rt° which relate to the original data.

While there are statistically significant benefits to using 4PL+C to analyze complete curves, the improvement is even more evident for incomplete curves. We obtained incomplete curves by removing data corresponding to the two highest concentrations of the simulated completed curves (Figure 2). Removing these data increases the MSE of all parameter estimates (Table 2). While the increase in MSE is approximately 4-fold for the standard 4PL procedure, the increase is subtle for 4PL+C; using the positive control appears to mostly compensate for losing data points at the highest concentrations. Because 3PLFB also uses the controls to set Rb, this procedure also compensates for the loss of these data points, and the precision is similar to 4PL+C.

Table 2 Comparison of Error in Parameter Estimates from Simulated Incomplete CRCsa

MSE	4PL+C	4PL	3PLFB	
Rb	1.0 × 10–1 (3.5 × 10–3)	6.6 × 101 (1.0)	 	
Rt	1.9 × 100 (3.9 × 10–2)	3.9 (8.2 × 10–2)	3.8 (7.9 × 10–2)	
x50	1.4 × 10–3 (3.1 × 10–5)	8.3 × 10–3 (1.8 × 10–4)	1.5 × 10–3 (3.8 × 10–5)	
H	5.8 × 10–3 (8.5 × 10–5)	2.1 × 10–2 (2.8 × 10–4)	6.1 × 10–3 (9.1 × 10–5)	
a Mean and standard deviation (in parentheses) of the MSE of parameters estimated from 5000 curves across 10 independent sets.

Figure 2 Example fits to incomplete data with (blue line) and without (orange line) controls. Curve data are shown with black circles and control data as green triangles. In the analysis of the incomplete curves, the points represented by red crosses are removed. The x position of the control data is not meaningful; data on the left correspond to x → –∞ and on the right correspond to x → ∞. The vertical lines indicate the true x50 (black dotted) or estimated x50 from CRC fitting with (blue dashed) or without (orange dashed) controls.

While it is clear that including control data will improve the estimation of parameters, the extent of improvement will depend on a large number of factors. These factors include the variance of the control and curve data. They also include the set of concentrations at which responses are measured, including the number of points along the curve and whether the curve is complete or incomplete. Hence the aforementioned reduction in the MSE may not be representative of all situations.

Our result that including control data is especially beneficial for the analysis of incomplete curves is consistent with the previous observations. Kappenberg et al.6 found that in the “difficult” situation of an incomplete curve, estimation with controls outperforms estimation without controls unless the controls are strongly deviating. Sebaugh4 found that fixing Rb or Rt can strongly influence estimates of x50 and H in incomplete curves. As mentioned in the introduction, fixing Rb or Rt privileges the control data by using it (and not the curve data) to set values for these parameters. The 4PL+C approach provides similar benefits for the analysis of incomplete curves without the subjective choice of which parameters to fix and without giving unwarranted privilege to the control data.

The Precision of Estimates from Repeated Experiments Is Consistent with Trends Observed in Simulations

The COVID Moonshot data include some repeated experiments. The fluorescence assay data include 9 independent CRCs for inhibition of MPro by a compound with database identifier CVD-0002707. The mass spectroscopy assay data include 25 independent CRCs for ebselen. The parameters estimated by applying 4PL versus 4PL+C are shown in Figure S3 of Supporting Information 1. In this figure, we compared 4PL+C to 4PL instead of 3PLFB because the former is more standard and is more widely used. For both data sets, 4PL+C reduces the standard deviation of estimated Rb and Rt compared to 4PL (Table S1 of Supporting Information 1). There was not a statistically significant difference between the standard deviation of pIC50 and Hill slope. These results are inconclusive but consistent with the trends observed in the simulated data.

There is a large variation in the estimated Hill slope for ebselen. This is likely because ebselen is a covalent inhibitor. For a covalent inhibitor, the results are dependent on the time between incubation and measurement, which could be highly variable. This variability does not affect the comparison between fitting procedures, as the Hill slope reported for each curve is relatively insensitive to the fitting procedure (Figure S3 of Supporting Information 1).

4PL+C Is Reasonably Fast and a Simplification Is Widely Available

Although curve fitting with the 4PL+C procedure is slower than that with 4PL, both are still quite fast. Performing fits for all 1668 fluorescence curves took 6 min and 2 s for 4PL+C but only 18 s for 4PL. For all 238 mass spectrometry curves, the fitting took 36 s for 4PL+C versus 8 s for 4PL. As described in the methods section, the 4PL+C procedure alternates between optimizing the Hill equation parameters θ and the estimated variances. In contrast, the 4PL procedure requires optimizing only the Hill equation parameters. Nonetheless, 4PL+C is quite fast, and the additional computing cost required for the curve fitting should not be a barrier to using the procedure.

While existing data analysis software do not, to our knowledge, implement 4PL+C, all nonlinear regression software implement a simplified version of 4PL+C. Some packages like drc(12) use model functions that are defined at a concentration of zero. Even in software that do not, Weimer et al.5 noted that it is straightforward to include control data by assigning them extreme concentration values such that the model response is a close approximation to the asymptotes. This procedure leads to a simplication of 4PL+C in which the controls are assumed to have the same variance as the curve. However, at least for the COVID Moonshot data, we see that this appears to be an unreasonable assumption.

Variances Differ between the Curve and Controls

In the fluorescence and mass spectroscopy CRC data collected for the COVID Moonshot, the estimated variance of the curve and the controls differ (Figure 3). In both assay types, the bottom control has the least variance. For the fluorescence measurements, the variance of the top controls is comparable to that of the curve. For mass spectroscopy experiments, the variance of the curve is often higher than the top control. In both experiments, there is no correlation between the estimated variance of the top control and the curve (data not shown).

Figure 3 Estimated standard deviations of the controls and of the curves. Histogram of estimated standard deviations of the bottom controls (green dashed dotted line), the top controls (blue dashed line), and the curve (red line) for (a) fluorescence and (b) mass spectroscopy experiments.

Differences in variance can be explained by experimental settings, which have been described elsewhere.20 Bottom controls have buffer but no enzyme and are subject to less variation due to concentration errors that may occur, for example, due to pipetting error or enzyme degradation. In the fluorescence measurements, the top controls and curves have comparable sources of error. In the mass spectroscopy experiments, the top control contains only DMSO (no inhibitor) and does not have error due to the inhibitor concentration.

The 4PL+C Procedure Reduces the Estimated ASE of Experimental CRCs

For the vast majority of the fluorescence and mass spectroscopy CRC data collected for the COVID Moonshot, 4PL+C yields a smaller asymptotic standard error (ASE) than 4PL. Examples of fitting curves using 4PL+C and 4PL can be found in Supporting Information 2. We fit 1589 of 1668 fluorescence curves and 205 of 238 mass spectrometry curves with a high coefficient of determination, R2 > 50%, using both procedures. In over 95% of the curves of both types, the ASE of Rt and Rb was lower using 4PL+C than 4PL (Figures 4 and 5). The 4PL+C procedure provides a lower ASE than 4PL in over 75% of the x50 and H estimates. The reduction in ASE suggests that the 4PL+C determines parameters more accurately than 4PL.

Figure 4 Comparison of ASE estimates for 4PL+C and 4PL analysis of fluorescence CRC from the COVID Moonshot. For each parameter, the fraction of data above the diagonal is Rb (99.69%), Rt (98.16%), x50 (89.51%), and H (98.49%).

Figure 5 Comparison of ASE estimates for 4PL+C and 4PL analysis of mass spectroscopy CRC from the COVID Moonshot. For each parameter, the fraction of data above the diagonal is Rb (100.00%), Rt (96.57%), x50 (76.96%), and H (88.24%).

Outlier Detection and Curve Refitting Can Improve the Precision of Parameter Estimates in Both 4PL and 4PL+C

For simulated CRCs, outlier detection and curve refitting did not improve the precision of the parameter estimates. Lenth’s method21 detected an outlier in about 40% of curves. However, removing the outliers and refitting had a negligible effect on the MSE of estimated parameters, yielding results identical with those in Table 1. The normally distributed error in the simulated curves does not lead to significant perturbations of parameter estimates.

For data from both the fluorescence and mass spectroscopy assays, removing outliers using Lenth’s method21 reduces the ASE for the vast majority of estimated parameters. Of the 1668 fluorescence assay curves, 590 had at least one outlier. To show the effect of outlier removal on the precision of each parameter, we plotted a histogram of the difference in the ASE before versus after outlier removal. For the majority of these curves, removing the outlier(s) and refitting the data led to a lower ASE. (Figure 6). Of the 238 mass spectroscopy curves, 94 had at least one outlier. For most of these, removing the outlier(s) and refitting led to lower ASE (Figure 7). Similar results were observed when refitting with the standard procedure (Figure S4 in the Supporting Information 1).

Figure 6 Histogram of the change in ASE after outlier detection and refitting using the 4PL+C procedure for data from the fluorescence assay. In most parameter estimates (Rb: 86.86%, Rt: 97.44%, x50: 98.46%, H: 94.37%), the ASE is reduced (blue), but in some, the ASE increases (orange).

Figure 7 Histogram of the change in ASE after outlier detection and refitting using the 4PL+C procedure for data from the mass spectroscopy assay. In most parameter estimates (Rb: 78.72%, Rt: 96.81%, x50: 94.68%, H: 89.36%), the ASE is reduced (blue), but in some the ASE increases (orange).

Conclusions

As demonstrated by analyses of simulation and experimental data, the proposed 4PL+C procedure results in more accurate estimates of Hill equation parameters, including the half-maximal concentration and Hill slope, than the established regression methods. Benefits are especially evident in the analysis of incomplete curves. Moreover, we have shown that an outlier detection and curve refitting leads to a lower ASE of the pIC50 and Hill slope for both 4PL+C and 4PL.

The proposed method provides clear benefits without a significant cost. No more data need to be collected than is required for standard statistical approaches. Implementation is only minimally more complex. Computational expense increases are not prohibitive. While most curve fitting software can perform a variant of 4PL+C in which all control and curve data are assumed to have the same variance, this appears to be, at least with the COVID Moonshot data, a bad assumption. Thus, we recommend that the described 4PL+C procedure, with different variances for the curve and top and bottom controls, be implemented by data analysis software developers. We recommend that it be applied by scientists as a new standard, replacing 4PL.

Experimental Section

Statistical Approaches

Maximum Likelihood Estimation without Control Data

Our procedure is based on the maximum likelihood estimation of parameters that fit the Hill equation to CRC data. We use to denote the experimental data on the curve . Here xi’s are base 10 logarithms of concentrations of the compound where response measurements are recorded, and we assume that there are m unique xi values. At each xi, there are ni repeated measurements of response ri,j for j = 1, ..., ni. Therefore, there are Nc=∑i = 1mni observations that are part of the concentration–response curve. The most common practice is to use the same number of replicates at each xi, such that ni = n for i = 1, ..., m and Nc = m × n.

To obtain the likelihood function of the CRC, we assume that the data follow the classical Hill equation in eq 1 with an additive measurement error:2

The measurement error ϵi,j independently and identically follows the normal distribution with constant variance σc2. Unknown parameters include θ = (Rt, Rb, x50, H) in the Hill equation and the nuisance parameter σc2. Based on these model assumptions, the log-likelihood of the CRC data is3

The maximum likelihood estimator (MLE) is the solution of the following maximization problem:4

Note that for the parameters that maximize the likelihood, the gradient of the log likelihood with respect to is zero. This condition may be expressed in two equations: and , where . Due to the structure of , the two equations can be solved separately. Since the Hill equation is nonlinear in θ, it is suitable to use a nonlinear equation solver to solve the first equation; this is the typical nonlinear regression approach. Thus, the MLE of is . However, the MLE of σc2 is actually biased. We instead use an unbiased estimator of σc2:

where p = 4 is the number of parameters in θ.

Maximum Likelihood Estimation Using Control Data

In many cases, it is reasonable to use control data as measurements of extreme responses. Consider an enzyme inhibition assay. For a negative control without any inhibitor, because it is the base 10 logarithm of the concentration, x is undefined, and the data technically are not part of the curve. However, it is reasonable to treat the negative control data as corresponding to an extremely low concentration such that x → –∞. For a positive control with substrate(s) and cofactor(s) but no enzyme, the response should be the same as if the enzyme was present but fully inhibited by an extremely high concentration of inhibitor, such that x → ∞.

We obtain the likelihood of the control data based on the assumptions that they correspond to extreme responses and that measurement error has a similar structure to the curve. The bottom control data include the observations rb,j for j = 1, ..., Nb and are considered as measurements in which x is extremely large, such that x → ∞. In this limit, the Hill equation yields Rb and measured responses are rb,j = Rb + ϵb,j, where ϵb,j’s are the measurement errors. Because the controls may be measurements of solutions containing substances different from the curve, we cautiously assume that the measurement errors follow different normal distributions for the control and curve data. Specifically, ϵb,j’s independently and identically follow the normal distribution , where σb2 is a constant, the variance of the measurement error. Like σc2, σb2 is an unknown nuisance parameter. At the other extreme, the top control data include the observations rt,j for j = 1, ..., Nt and correspond to x → –∞. Based on the Hill equation in this limit, rt,j = Rt + ϵt,j, where and are independently and identically distributed. Based on these assumptions, the log-likelihood of a bottom controls is,5

and of the top controls is6

The log-likelihood of curve and control data is simply the sum of the log-likelihood of each group of data:7

The MLE of the parameters based on both the curve and control data is the solution of8

The two estimators and are different as the former is only based on the data set Dc, whereas the latter is based on Dc, Dt, and Db. Moreover, the estimators maximize different likelihood functions (eq 4 versus eq 8).

Solving eq 8 is slightly more complicated than solving eq 4. Maximizing is equivalent to minimizing R(θ,σc2,σt2,σb2), where9

At the minimal R(θ,σc2,σt2,σb2), the gradients are zero, such that , , and . Compared to the MLE without control data, the equations for θ and (σc2,σt2,σb2) are more intertwined. Therefore, we iteratively solve these equations in terms of θ for fixed (σc2,σt2,σb2) values and then update (σc2,σt2) values using the latest solution of θ. Iteration are stopped when some convergence conditions are achieved: when the changes in R function values and parameter values are less than a prespecified tolerance value. Denoting as the optimal solution for θ, the estimates for the three variances are10

11

12

Compared to the fitting without control data, it is more complicated to adjust the estimates to correct the bias for . For larger sample sizes, such as in this study, the bias is not substantial.

Inference of Asymptotic Standard Error

Next, we discuss the asymptotic standard error of the two MLEs. In general, the asymptotic distribution of an MLE θ̂ is normal with zero mean and covariance given by the inverse of the Fisher information matrix, i.e.,13

where N and θ are the generic notation for sample size and parameters, and θ0 is the vector of the true parameter values. The matrix I(θ) is the Fisher information matrix evaluated at θ and defined as

Since θ0 is unknown, I(θ0) is usually approximated by the sample mean,14

where Li is the likelihood of the ith observation. Consequently, the asymptotic covariance of θ̂ is15

The ASE of each parameter is simply the diagonal of the asymptotic covariance matrix. Therefore, a confidence interval of (1 – α) × 100% is given by , where zα/2 is the number of standard deviations (z-score) for the given value of α/2 such that the probability of a variable being contained in the interval is Pr(Z ≤ zα/2) = 1 – α/2.

To infer the ASE for the MLE based on data excluding controls (eq 4), we used eq 3 in the general ASE expression (eq 15). This leads to,16

The gradient and Hessian of Q are given by

Detailed formulas for gradient and Hessian of r(xi, θ) are given in Appendix S1 of Supporting Information 1. To compute , we can plug in either the MLE or unbiased estimator for .

To infer the ASE for the MLE based on data including controls (eq 8), we use eq 7 in the general ASE expression. In this case, the covariance of θctb is approximately,17

where Ntotal = Nc + Nt + Nb is the total number of observations from all the data sets. The gradient and Hessian of R are given in Appendix S2 of Supporting Information 1. We plug in the estimates eqs 10, 11, and 12 to compute .

Normalization

Thus, far in the Statistical Approaches, we have not discussed the effect of normalization. It is a common practice in the analysis of the CRC experiments to normalize the data prior to fitting the Hill equation model. Normalization helps “accommodate minor changes to assay protocols to ensure robustness”, a key goal of the AGM.2 Minor changes in the assay protocols could include variations in the amount of time between the initiation of the enzymatic activity and measurement of the product concentration. Here we address the question of how normalization affects the MLE estimators and .

Normalization is a linear transformation of the original data ri,j° for j = 1, ..., ni and i = 1, ..., m,18

to yield a percentage response. As explained in the chapter “Data Standardization for Results Management” in the AGM, in most assays, the bottom and top of the curve are based on the controls. Usually, and are sample means of the two sides’ control data. It is also possible to use their medians, which are less sensitive to outliers.

Based on the same model as in eq 2, the original data are denoted by ri,j° and19

where r◦(x, θ) is,

where and . Compared to the normalized measurement error, the original measurement error ϵi,j° is scaled, . Hence, the variance of ϵi,j° is . Thus, if we use the original data {ri,j°,xi} for j = 1, ..., ni and i = 1, ..., m, the corresponding parameters are θ°=(Rb°, Rt°, x50, H) with Rb° and Rt° defined above and x50 and H remaining the same as in the normalized formulation.

The procedure for obtaining the MLE is exactly the same as for normalized data (see the sections entitled Maximum Likelihood Estimation without Control Data and Maximum Likelihood Estimation Using Control Data) except the data are replaced with original data and the parameters are changed to θ◦. The log-likelihood of the original data (excluding control data) is20

It is easy to see that the MLE for θ◦, denoted by , has the same relationship with , as θ′ with θ, and the MLE for x50 and H remains the same whether the data are normalized or not.

If the control data are included in the estimation, would the normalization affect the estimator ? No, normalization does not affect the estimator either. Note the relationships:

where and , , and Rb° and Rt° are previously defined. Therefore, the models for the normalized control data follow the same model as the unnormalized ones as long as the parameters are transformed accordingly. Using the same logic as for , we can directly conclude that the MLE using the original data have the same estimates for x50 and H as using the normalized data. For Rb and Rt, using either one of the two approaches, it is always the case that and .

Data Analysis

Data

We analyzed simulated and experimental data. Simulations are advantageous because we know the true values and can repeat calculations many times. For experimental data, even without true values and as many repetitions, it is helpful to compare the estimated ASE for different statistical analyses. We analyzed experimental data from the COVID Moonshot project.20 The COVID Moonshot was a worldwide collaboration to develop an open-source antiviral that treats COVID-19 by targeting the SARS-CoV-2 main protease. We analyzed CRC in which protease activity was measured by fluorescence (1668 curves) and mass spectrometry (238 curves). These data sets can be found in an Excel spreadsheet as Supporting Information 2.

CRC data were simulated by adding random Gaussian noise to the classical Hill equation. Parameters for the Hill equation were selected to be comparable to mass spectrometry measurements of enzymatic inhibition assays from the COVID Moonshot:20Rb° and Rt° were set at 0.227 and 54.074, while H was fixed at 1 and x50 at −6.0759. For the simulated curves, 10 concentrations were used, with the highest concentration at 100 μM and the others diluted by a factor of 5 each. Random Gaussian noise with mean 0 and variance was added. The number of replicates for each concentration on the curve was n = 4 and the number of observations of each control was Nb = Nt = 16. Simulations were performed in low-variance (, ) and high-variance (, ) conditions. In each condition, 10 sets of 5000 curves were simulated. Incomplete curves were obtained by removing data corresponding to the two highest concentrations from simulations in the high-variance condition.4

Curve Fitting

All of the statistical approaches, 4PL+C, 4PL, 3PLFB, 3PLFT, and 2PL, were implemented in Python. The code is freely available at https://github.com/vanngocthuyla/crc/tree/main/nls_crc. Maximization of the log likelihood was performed using functions curve_fit and minimize from scipy.optimize.22 Initial values of Rb and Rt were set at the minimum and maximum values of the percentage activity. During the fitting procedure, these parameters were restricted to their initial values, plus/minus a quarter of the data range. For x50, the initial value was set to be observed concentration with the response closest to the mean of the minimum and maximum response. The initial Hill slope was set to 1. x50 and H were restricted to the ranges from −50 to 0 and from 0 to 50, respectively. Optimization was performed until the absolute difference in the log likelihood was less than a tolerance of 1 × 104. If this tolerance was not achieved before 100 iterations, no fitted values were reported.

Fitting was performed on a computer cluster with the following key specifications: a dual AMD Opteron 6376 processor with a clock speed of 2.30 GHz and 16MB cache, 32 cores, and 64 GB RAM. The job required one core and 2 GB of RAM.

Outlier Detection and Curve Refitting

In addition to analyses without preprocessing, we also performed some analyses of experimental data with automated outlier detection based on Lenth’s method.21 Starting with the residual between the observed data and model values from the initial fit, , s0 = 1.5median|ϵi,j|. The pseudostandard error (PSE) is defined as21

For every data point, the coefficient of variance is the ratio of the residual, ϵi,j, over the PSE22

A data point with CV larger than 1.5 times the interquartile range—the distance between 25th (Q1) and 75th (Q3) quantiles of all CVi,j—is considered to be an outlier and is removed.23 If a curve has at least one outlier, the outliers are removed, and the curve is refit.

Supporting Information Available

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jmedchem.3c00107.Appendices for the gradients and Hessians of functions reference in the main text, histograms of parameter estimates from simulated CRCs, parameter estimates and precision thereof from repeated experiments, and a histogram summarizing the effects of outlier detection in the ASE of the 4PL procedure (PDF)

CRCs from the COVID Moonshot and example fits with and without control data (XLSX)

Supplementary Material

jm3c00107_si_001.pdf

jm3c00107_si_002.xlsx

The authors declare no competing financial interest.

Acknowledgments

We thank John Chodera (ORCID 0000-0003-0542-119X) for helpful discussions. Chodera asked us to analyze CRC data from the COVID Moonshot, inspiring us to start this project. Later, he made a suggestion that we look to the AGM as the standard in the field. We also thank members of the COVID Moonshot consortium for providing access to their data. This project was supported by National Science Foundation Award 1905324 to J. Chodera, D. Minh, and L. Kang. S. Nicholson was supported by the RES-MATCH undergraduate research program of the Pritzker Institute of Biomedical Science and Engineering at Illinois Tech.

Abbreviations Used

2PL two parameter logistic

3PLFB three-parameter logistic fixed bottom

3PLFT three-parameter logistic fixed top

4PL four parameter logistic

4PL+C four parameter logistic with controls

AGM Assay Guidance Manual

ASE asymptotic standard error

COVID coronavirus disease

CRC concentration response curve

MLE maximum likelihood estimator

MSE mean square error

SARS-CoV-2 severe acute respiratory syndrome coronavirus 2
==== Refs
References

Krippendorff B.-F. ; Lienau P. ; Reichel A. ; Huisinga W. Optimizing Classification of Drug-Drug Interaction Potential for CYP450 Isoenzyme Inhibition Assays in Early Drug Discovery. Journal of Biomolecular Screening 2007, 12 , 92–99. 10.1177/1087057106295897.17130250
Assay Guidance Manual eBook; Markossian S. , Eds.; Eli Lilly & Company and the National Center for Advancing Translational Sciences: Bethesda, MD, 2004. https://www.ncbi.nlm.nih.gov/books/NBK53196/ (accessed 2022-07-01).
Assay Guidance Manual; National Center for Advancing Translational Sciences: Bethesda, MD, 2016. https://ncats.nih.gov/expertise/preclinical/agm (accessed 2022-07-01).
Sebaugh J. L. Guidelines for Accurate EC50/IC50 Estimation. Pharmaceutical Statistics 2011, 10 , 128–134. 10.1002/pst.426.22328315
Weimer M. ; Jiang X. ; Ponta O. ; Stanzel S. ; Freyberger A. ; Kopp-Schneider A. The Impact of Data Transformations on Concentration–Response Modeling. Toxicol. Lett. 2012, 213 , 292–298. 10.1016/j.toxlet.2012.07.012.22828011
Kappenberg F. ; Brecklinghaus T. ; Albrecht W. ; Blum J. ; van der Wurp C. ; Leist M. ; Hengstler J. G. ; Rahnenfuhrer J. Handling Deviating Control Values in Concentration-Response Curves. Arch. Toxicol. 2020, 94 , 3787–3798. 10.1007/s00204-020-02913-0.32965549
Dose Response Relationships/Sigmoidal Curve Fitting/Enzyme Kinetics. https://www.originlab.com/index.aspx?go=Solutions%2FCaseStudies&pid=106 (accessed 2023-04-18).
GraphPad Prism 9 Curve Fitting Guide - Pros and Cons of Normalizing the Data. https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_pros_and_cons_of_normalizing.htm (accessed 2023-04-08).
Motulsky H. ; Christopoulos A. Fitting Models to Biological Data Using Linear and Nonlinear Regression. A Practical Guide to Curve Fitting; Oxford University Press: New York, NY, 2004.
Di Veroli G. Y. ; Fornari C. ; Goldlust I. ; Mills G. ; Koh S. B. ; Bramhall J. L. ; Richards F. M. ; Jodrell D. I. An Automated Fitting Procedure and Software for Dose-Response Curves with Multiphasic Features. Sci. Rep. 2015, 5 , 14701 10.1038/srep14701.26424192
Gadagkar S. R. ; Call G. B. Computational Tools for Fitting the Hill Equation to Dose–Response Curves. Journal of Pharmacological and Toxicological Methods 2015, 71 , 68–76. 10.1016/j.vascn.2014.08.006.25157754
Ritz C. ; Baty F. ; Streibig J. C. ; Gerhard D. Dose-Response Analysis Using R. PLoS One 2015, 10 , e0146021 10.1371/journal.pone.0146021.26717316
Johnstone R. H. ; Bardenet R. ; Gavaghan D. J. ; Mirams G. R. Hierarchical Bayesian Inference for Ion Channel Screening Dose-Response Data. Wellcome Open Research 2016, 1 , 6 10.12688/wellcomeopenres.9945.2.27918599
Wang Z.-J. ; Liu S.-S. ; Qu R. JSFit: A Method for the Fitting and Prediction of J- and S-shaped Concentration–Response Curves. RSC Adv. 2018, 8 , 6572–6580. 10.1039/C7RA13220D.35540430
Vivaudou M. eeFit: A Microsoft Excel-embedded Program for Interactive Analysis and Fitting of Experimental Dose–Response Data. BioTechniques 2019, 66 , 186–193. 10.2144/btn-2018-0136.30987445
Beam A. L. ; Motsinger-Reif A. A. Optimization of Nonlinear Dose- and Concentration-Response Models Utilizing Evolutionary Computation. Dose-Response 2011, 9 , 387–409. 10.2203/dose-response.09-030.Beam.22013401
Wang Y. ; Jadhav A. ; Southal N. ; Huang R. ; Nguyen D.-T. A Grid Algorithm for High Throughput Fitting of Dose-Response Curve Data. Current Chemical Genomics 2010, 4 , 57–66. 10.2174/1875397301004010057.21331310
Nguyen T. T. ; Song K. ; Tsoy Y. ; Kim J. Y. ; Kwon Y.-J. ; Kang M. ; Edberg Hansen M. A. Robust Dose-Response Curve Estimation Applied to High Content Screening Data Analysis. Source Code for Biology and Medicine 2014, 9 , 27 10.1186/s13029-014-0027-x.25614758
R Core Team R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021.
The COVID Moonshot Consortium, Open Science Discovery of Oral Non-Covalent SARS-CoV-2 Main Protease Inhibitor Therapeutics. bioRxiv, January 20, 2022, ver. 1.10.1101/2020.10.29.339317 (accessed 2020-12-21).
Lenth R. V. Quick and Easy Analysis of Unreplicated Factorials. Technometrics 1989, 31 , 469–473. 10.1080/00401706.1989.10488595.
Virtanen P. ; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods 2020, 17 , 261–272. 10.1038/s41592-019-0686-2.32015543
Tukey J. W. Exploratory Data Analysis; Addison-Wesley, 1977.
