Estimating Potency in High-Throughput Screening Experiments by Maximizing the Rate of Change in Weighted Shannon Entropy

High-throughput in vitro screening experiments can be used to generate concentration-response data for large chemical libraries. It is often desirable to estimate the concentration needed to achieve a particular effect, or potency, for each chemical tested in an assay. Potency estimates can be used to directly compare chemical profiles and prioritize compounds for confirmation studies, or employed as input data for prediction modeling and association mapping. The concentration for half-maximal activity derived from the Hill equation model (i.e., AC50) is the most common potency measure applied in pharmacological research and toxicity testing. However, the AC50 parameter is subject to large uncertainty for many concentration-response relationships. In this study we introduce a new measure of potency based on a weighted Shannon entropy measure termed the weighted entropy score (WES). Our potency estimator (Point of Departure, PODWES) is defined as the concentration producing the maximum rate of change in weighted entropy along a concentration-response profile. This approach provides a new tool for potency estimation that does not depend on the assumption of monotonicity or any other pre-specified concentration-response relationship. PODWES estimates potency with greater precision and less bias compared to the conventional AC50 assessed across a range of simulated conditions.

Quantitative high-throughput screening (qHTS) assays 1 return thousands of concentration-response profiles for large chemical libraries and are currently driving major advancements in drug discovery 2 and toxicity testing 3 . For example, more than 10,000 substances are now being tested in 15-point concentration-response format in phase II of the Tox21 collaboration, involving the U.S. Environmental Protection Agency (EPA), the U.S. Food and Drug Administration (FDA), the National Institutes of Health (NIH) National Center for Advancing Translational Sciences (NCATS) and the National Institute for Environmental Health Sciences (NIEHS)/National Toxicology Program (NTP) 4 . Response profiles can be summarized by a measure of average activity across tested concentrations, such as the area under the curve (AUC) of concentration-response curves 5 , a weighted version of AUC 6 , or a weighted entropy score (WES) 7 . While these measures are useful for ranking compounds, it is often desirable to estimate the concentration at which a chemical induces a particular effect level using automated data analysis processes. Such potency measures can be applied for rapid identification of pharmacoactive hits or toxicological assessment, or used as input data for prediction modeling 8 or association mapping 5 .
The most common approach used to approximate chemical potency in chemical genomics and large-scale toxicity testing is the AC 50 parameter in the Hill Equation model 9 . The AC 50 parameter estimates the concentration at which a chemical produces the half-maximal response along a sigmoidal curve 10 . Incorporating domain knowledge into the curve fitting process can improve agreement between AC 50 estimates for sigmoidal curves 11 . However, it is not possible to know the underlying shape of the concentration-response relationship before conducting an experiment 12 and complex response patterns may reflect real biological responses 13 . Furthermore, linearizing assumptions can render AC 50 parameter estimation from the Hill model very unreliable, even with increased sample sizes 10,14 . Applying individualized curve fitting procedures can be useful for characterizing screening results. However, in the high-throughput setting manual scrutiny can be restrictively laborious and result in extensive data censoring. Also, while outlier removal and parameter constraints may reduce curve fit error, these procedures do not necessarily increase the repeatability of nonlinear parameter estimation. It is not unusual for AC 50 estimates to be accompanied by large standard errors even when one or both asymptotes can be defined 10 .
A point of departure (POD) represents a concentration derived from observed concentration-response data that is associated with a defined effect. In vitro POD estimates have been calculated based on linear interpolation between the two concentrations that lie on either side of the assay detection threshold 6 or establishing a baseline noiseband using the first two tested concentrations 15 . Other POD metrics include an estimate of the concentration producing a predetermined level of an adverse response (i.e., the benchmark dose or BMD) and the highest tested concentration for which there is no observed adverse effect (i.e., the no-observed-adverse-effect-level or NOAEL) 16,17 . With true experimental replicates, BMD modeling or NOAEL determinations could serve as POD estimates describing the concentration at which the assay response begins to deviate from baseline response levels. Unlike the NOAEL approach, the BMD procedure uses mathematical modeling to make use of the entire observed concentration-response profile. Unfortunately, in qHTS studies there is usually very little, if any, replication at each tested concentration and it is often not appropriate to combine data across different experimental runs because conditions can change substantially between trials 4,10 .
We propose a nonparametric approach based on information theory to improve the precision of compound potency estimation in qHTS studies. Information theoretic concepts were originally developed for communication technology 18 , but these approaches have recently been used to summarize patterns in gene expression microarray data 19,20 , find differential methylation sites 21 and rank chemicals in qHTS experiments 7 . Shannon entropy (H) describes the average information content in a probability distribution 22 , and can be used to describe the extent and uniformity of response in a concentration-response profile. Here, H is computed from the probability distribution obtained from the observed responses and naturally accommodates any concentration-response pattern, not just monotonic trends such as the sigmoidal shape of the Hill equation model.
We define compound potency as the concentration producing the maximal rate of change in entropy. This potency is calculated by finding the maximum first derivative of the entropy measure across the concentration range. However, Shannon entropy does not take into account the uncertainty in response measurements when responses are within the noise region, i.e., measurements that are less than the assay detection limit. We therefore employ a weighted version of Shannon entropy (or WES) 7 . WES weights responses found within the noise region so that profiles with larger WES scores have greater probability mass (i.e., greater average activity) in the detectable region of the assay. Accordingly, the point of departure is found at the concentration where the rate of change in weighted entropy is maximized along the tested concentration range. This new potency estimator is termed POD WES . Unlike the AC 50 value, POD WES does not rely on the shape of the profile far removed from the point of departure. Observed concentration-response profiles that lie entirely within the assay noise region are assigned the outcome "undefined". Profiles which have detectable responses and for which the maximum rate of change in weighted entropy is located at the lowest observed concentration C 1 , where POD WES must be less than C 1 but cannot be estimated from the given data, are assigned the outcome "less than C 1 ".

Results
Computing POD WES for illustrative profiles. Figure 1 summarizes the workflow used to calculate POD WES . To begin, WES and its derivatives are calculated at each tested concentration level. Chemicals with larger WES scores have greater average relative responses across concentrations 7 . If the maximum observed response is less than the assay detection limit, POD WES is "undefined", since a detectable response may have occurred if a larger range of test concentrations had been used. If at least one measured response is detectable, a search for a maximal rate of change in WES is conducted within the observed concentration-response space. If a global dWES dC max extremum is located, POD WES is estimated. However, if POD WES cannot be found, the concentration-response data is extrapolated outside of the observed concentration-response region using finite difference calculus. After extrapolating new responses, WES and its derivatives are recalculated and another search for POD WES is conducted. If POD WES still cannot be quantitatively determined, but dWES dC max is located at the lowest concentration in the extrapolated profile, POD WES must be less than the lowest tested concentration (see Supplementary  Information). Figure 2 depicts hypothetical sigmoidal response profiles for three chemicals. Each chemical follows equation (1) in the Methods with no ERROR. The baseline response R0 is set to 0% of positive control, the maximal response |RMAX| is set to 100% of positive control, the h parameter is set to 1 and the AC 50 is set to 0.001, 0.1, and 10 μ M, for Chemical-1, Chemical-2 and Chemical-3, respectively. This figure shows the normalized responses (row 1), WES computed at each concentration level (row 2), the first derivative of WES at each concentration level (row 3) and the second derivative of WES at each concentration level (row 4). The concentration at which the first derivative of WES is maximized is indicated by an open square.
Chemical-1 is the most potent of the three chemicals shown in Fig. 2, where only the upper asymptote is well defined. This chemical has a "true" AC 50 value equal to 1.00 × 10 −3 μ M, which corresponds to an POD WES of 4.19 × 10 −4 μ M. Chemical-2 has two clearly defined asymptotes with an AC 50 value of 0.1 μ M and a calculated POD WES of 0.07 μ M. One data point, indicated by an open triangle, was extrapolated in order to find POD WES for Chemical-3, which had an AC 50 value of 10 μ M and a calculated POD WES of 3.73 μ M. In this case, a single extrapolated data point was used in order to calculate the deviation of the estimate of d WES dC 2 2 from zero within the prespecified tolerance level (see Supplementary Information for more explanation of the computations). Notice that the value of WES at the kth concentration level becomes smaller as the AC 50 of a profile increases, but the potency Scientific RepoRts | 6:27897 | DOI: 10.1038/srep27897 measure POD WES is located at the concentration for which the the rate of change in WES is increasing most rapidly. Fig. S1 shows additional examples of POD WES calculated from curves generated with the "gain-loss" model given in equation (2) in the Methods.
Evaluating the proposed approach using simulated data. To explore precision and bias of the potency estimates derived from sigmoidal models, we generated 15-point concentration-response profiles from equation (1) in the Methods with R0 = 0% and h = 1 for profiles having [1] only an upper asymptote (AC 50 = 0.001 μ M), [2] both asymptotes well defined (AC 50 = 0.1 μ M) and [3] only a lower asymptote (AC 50 = 10 μ M). In the simulations, |RMAX| values were selected as weak (|RMAX| = 25%), moderate (|RMAX| = 50%), and strong (|RMAX| = 100%) activity. A total of 10,000 profiles were generated for each of these nine combinations of AC 50 and |RMAX| and the residual errors were modeled as ERROR ~ N( μ = 0, σ i 2 ) where σ i = 5% or σ i = 10%. In Table 1, the precision of potency estimation differed markedly between the estimators for the lower error of σ i = 5%. POD WES estimates were generally more repeatable with confidence interval widths (CIWs) ranging from 1.03 to 1.53 orders of magnitude (OM). AC 50 at the same error level ranged from 0.27 to 13.80 OM. The precision of POD WES at σ i = 10% was comparable to the levels achieved at σ i = 5% for curves simulated under conditions in which "true" maximum response is greater than the detection limit of the assay. By contrast, precision of the AC 50 estimator was noticeably lower for σ i = 10% compared to σ i = 5%.
As shown in Table 1, for σ i = 5%, the bias in POD WES estimates was less than about 2.0-fold on the natural scale, ranging from log 10 Bias = − 0.03 (< 1.1-fold less than expectation) to log 10 Bias = − 0.27 (< 1.9-fold less than expectation). For the same data sets, the bias in AC 50 estimates ranged from log 10 Bias = − 0.0002 (< 1.1-fold less than expectation) to log 10 Bias = 1.38 (24.0-fold greater than expected). The estimation bias of POD WES at σ i = 10% was similar to the values found at σ i = 5%. By contrast, the estimation bias of AC 50 was about 10-fold greater for σ i = 10% compared to σ i = 5% in some instances.
In addition to investigating the precision and bias of potency estimators based on the sigmoidal Hill model, we also investigated the precision and bias of estimators using the "gain-loss" model from equation (2) in the represents the first derivative of WES with respect to concentration, and R extr is the vector of response values obtained after data extrapolation.
Methods. As shown in Table 2, the precision of POD WES was less than 1.5 OM at σ i = 5% and σ i = 10%. By contrast, the AC 50 measure could be extremely imprecise for these curve shapes, even reaching 19.78 OM in one case. Similar to the evaluation of Hill model curves, the log 10 Bias of POD WES did not exceed − 0.42 (< 2.7-fold less than expectation). The bias of the AC 50 metric was often greater than 2.0 OM (or 100-fold).
Example compound potency estimates across experimental runs. Previously, the in vitro BG1 estrogen receptor alpha (ERα) assay from phase II of Tox21 was used to screen for agonist activity in an ER reporter gene cell line with an endogenous full length ERα. Approximately 10,000 compounds were assayed in three different experimental runs and activity measurements for 15-point concentration-response curves were obtained as luciferase activity readings from the BG1 ERα cell line 23 . This data was normalized to 100% of the activity of estrone positive control compounds. Weighted entropy scores (WES) and POD WES values were calculated as described here. Ranking profiles based on WES is not based on any pre-specified concentration response model form or direction of response 7 .  Figure 3 shows concentration-response profiles for four representative compounds that are tested once in each of the three experimental runs. Estradiol valerate is a synthetic ester of the positive control compound 17β-estradiol and is consistently ranked within the top ten compounds based on WES. The corresponding potency value for this compound (POD WES ) is assigned a value of "less than the lowest tested concentration" in each run. Gestrinone is a synthetic hormone that elicits an agonistic response of 0.05 ± 0.03 μ M (across runs) in  (2) in the Methods was used to simulate data sets with the specified error. Log 10 Precision is shown with log 10 Bias in parentheses. Only data sets with true responses above the assay detection limit were evaluated.

Table 2. Precision and bias of potency metrics in gain-loss models. Equation
Scientific RepoRts | 6:27897 | DOI: 10.1038/srep27897 this experiment, and is ranked in the top 250 compounds based on WES in each run shown here. As shown in Fig. 3, the response profile for this compound is better represented by the "gain-loss" model than the Hill model, perhaps due to cytotoxic effects at the greater concentrations. The next compound, 4-Nonylphenol, has previously been shown to act as an agonist of the estrogen receptor alpha in MCF7 breast cancer cells 24 . This compound is ranked within the top 1,200 profiles based on WES and has a corresponding in vitro potency of 17.7 ± 8.8 μ M. Finally, the biocide 2-Phenylphenol does not show detectable activity in the assay in any experimental run and is therefore ranked very low based on the WES score in each case. The potency of this compound is assigned the value of "undefined. " The reproducibility of the potency estimates in this data set was evaluated by calculating log 10 potency differences between intra-assay duplicates and inter-assay duplicates interrogated across experimental runs. It is expected that duplicated compounds will have a log 10 mean ratio of 0, which corresponds to a mean ratio of 1.00 on the natural scale. All duplicated chemicals that had at least one observed response above the assay detection limit were included in the analysis. A shown in Fig. 4, there is less variation in log 10 potency differences for intra-array duplicates compared to inter-array duplicates as assessed by the median absolute deviation from zero. The dispersion of log 10 potency differences is noticeably greater for the AC 50 value compared to POD WES , indicating that AC 50 values are less reproducibile potency estimates in this experiment.

Discussion
High-throughput screening of compounds for biological activity can play a fundamental role in the advancement of drug discovery 25 and in efforts to transform toxicology from a mostly observational science into a predictive science 26 . Large-scale qHTS data analyses typically proceed by fitting the Hill equation 9 to the data and utilizing the AC 50 value as an estimate for compound potency. However, the uncertainty (e.g., confidence intervals) of these nonlinear parameter estimates can be extremely large and potentially limit the reproducibility of results obtained from qHTS studies 10 . A new procedure is proposed here to estimate compound potency based on locating the maximum rate of change in weighted entropy. This approach (POD WES ) provides more precise estimates of potency than typically obtained by nonlinear parameter estimation from the Hill model.
Regardless of the level of error used to simulate the concentration-response curves, under most circumstances potency measures examined here were subject to empirical confidence interval widths spanning at least one order of magnitude. However, the CIW for AC 50 estimates extended to greater than 13 orders of magnitude for low efficacy compounds at |RMAX| = 25% (see Table 1). Even so, the CIW of POD WES was less than 1.53 orders of magnitude (less than 34-fold) for data simulated from a Hill equation model or a "gain-loss" model. The bias in POD WES estimates was less than 2.7-fold (or |log 10 Bias| ≤ 0.42), and usually less than 1.5-fold. AC 50 estimates showed less bias than POD WES for Hill model curves generated with two clearly defined asymptotes, but bias was much greater when the data was generated under other conditions. Across-run comparisons of potency can be more variable than within-run comparisons (see Fig. 4). However, high-throughput screening responses can be affected by both random and systematic error, and run-to-run variability should be not ignored 27 . Obtaining experimental replicates can increase the precision of the potency estimates 28 and the interpretation of POD WES may be improved by focusing on robust assays with good agreement between compound measurements 29 and using appropriate signal curation 6 . If potencies are only desired from a pre-specified functional form (e.g., the Hill model), a two-step procedure can be used to (1) find response profiles that are active according to a robust analysis framework designed to detect the desired trend 30 and then (2) estimate potencies from the returned profiles.
The repeatability of AC 50 estimates can be extremely small for compounds with low efficacy or for situations in which one of the asymptotes cannot be established 10 . Furthermore, the assay detection limit can impact the precision of potency estimation. Using 3σ of the negative controls as a detection limit is a common practice in qHTS studies 1,6,28 and the 3σ value performed optimally in our simulation study across a range of σ values when considering bias, precision and the number of profiles with estimable potency values according to the Hill model (Fig. S2) and the "gain-loss" model (Fig. S3). Selective elimination of influential observations will not overcome these issues and may introduce bias because the true concentration-response relationship cannot be known in advance. Difficulties in characterizing the uncertainty of potency estimates derived from pre-specified models may be compounded when response profiles deviate from monotonicity or the incorrect model is employed for nonlinear curve fitting.
Each compound in a qHTS assay can be expected to have a distinctive set of parameters governing its response behavior. However, the approach proposed here to estimate potency using POD WES does not rely on a pre-specified concentration-response pattern, can be applied to complex response patterns without respect to the direction of response and naturally accommodates missing data into the estimation framework.

Methods
This section describes the procedures used to estimate compound potency based on maximizing the rate of change in weighted entropy. Data sets are simulated based on the Hill equation in order to evaluate the precision and bias of estimated potencies across a range of parameter space characterizing qHTS studies. In addition, the new potency measure is applied to an experimental data set assaying for estrogen receptor agonist activity from phase II of Tox21. Description of simulated data. Similar to previous studies 7,10,30 , concentration-response data sets were simulated using the logistic form of the four-parameter Hill equation model, where R i is a normalized response representing a percentage of the positive control activity at concentration C i . RMAX is the maximal response, R0 is the minimal response, AC 50 is the concentration of half-maximal response, h affects the shape of the curve and ERROR is residual error of the model. The logarithm in equation (1) ensures that back-calculated estimates of AC 50 obtained from log 10 AC 50 are restricted to positive values 10 . The concentrations (C i ) are based on equivalent log 10 concentration spacing ranging from 0.0001 to 100 μ M for fifteen-point concentration-response curves. The values of RMAX and AC 50 were set to (25,50, 100% of positive control activity) and (10 −3 , 10 −1 , 10 μ M), respectively, for a total of 9 different data sets. The R0 parameter was set to 0 and h was set to 1. Other data sets were simulated using a "gain-loss" model defined as the product of two Hill equation models,  (1) and (2) were modeled as ERROR ~ N(0,σ 2 ) with σ = 5% or 10%, where σ is related to the percent of negative control responses producing variation levels often seen in Tox21 Phase II assays 6 . Unless otherwise noted, the assay detection limit is taken to be 3σ , a typical detection limit in HTS studies 1,6,28 . A total of 10,000 simulated substances (RMAX = 25%, 50%, or 100% of positive control activity) were simulated for each data set.
Description of estrogen receptor agonist data set. We acquired qHTS data involving approximately 10,000 compounds that were screened for estrogen receptor alpha agonist activity 23 . This screen utilized an endogenous full length estrogen receptor alpha (BG1 cell line) with a luciferase reporter gene producing a single-channel readout 23 . A total of 15 concentrations were evaluated with concentrations typically ranging from ~10 −3 μ M to ~78 μ M. As part of phase II of Tox21, the library is screened three times with compounds located in different well positions during each experimental run 4 . The raw plate reads were normalized using the positive and negative control wells and subsequently corrected for row, column, and plate effects using linear interpolation 23 . Hill equation parameter estimates and activity calls were determined as described previously 30 . In order to assess within-run reproducibility, a set of 88 broadly active duplicates were deliberately included in the Tox21 Phase II 1,536-well assay plates. Concentration-response patterns in this experimental data set encompass many different types of patterns which may deviate substantially from sigmoidal profiles.
Weighted entropy score. The weighted entropy score provides a measure of average relative activity across a concentration-response profile 7 . Briefly, the response vector for a given substance R N = (R 1 , … , R N ) contains an observed response R i for each of N concentrations, where R i corresponds to the response at the ith concentration, C i . The relative response at C i is defined as.
The relative responses p i define a probability mass distribution based on the magnitude of R i , where R i may be positive or negative for activation or inhibition, respectively 7 . The entropy of R i , or surprisal of the ith event, is defined as . When every response value is greater than or equal to the assay detection limit, all w i = 1 so that WES is the same as Shannon entropy (i.e., WES = H = −Σ p i log 2 p i ). However, when R i values are less than the assay detection limit of 3σ the weights w i are defined as the ratio of the surprisal frequency for a relative response within the assay noise region (i.e., −p i,noise log 2 p i,noise , where p i,noise = p i /3σ ) divided by the uncorrected surprisal frequency (i.e., -p i log 2 p i ), or across concentration levels 7 . The entropy at the kth tested concentration (H k or WES k ) is computed by considering only the responses R k = (R 1 , … , R k ) that are observed within the full concentration-response profile R N .
The POD Approach to Estimate Potency. We define the profile-specific potency (denoted Point of Departure, POD) as the concentration along the response profile at which the magnitude of the rate of change in WES is greatest. This maximum rate of change defines the potency regardless of the direction of change, i.e., irrespective of whether the chemical is an activator or inhibitor. The rate of change in WES across the profile is computed as the derivative of WES with respect to concentration, or dWES dC , where concentration is based on log 10 units. In mathematical terms, POD WES is located at the concentration with the maximal value of dWES dC where d WES dC 2 2 is equal to zero and either (a) d WES dC 2 2 changes sign from positive to negative (for activation) or (b) d WES dC 2 2 changes sign from negative to positive (for inhibition) according to "The First Derivative Test" 31 . We compute the derivatives of WES using finite difference calculus, a mathematical procedure based on a Taylor series procedure that provides difference formulas for a grid sampled at discrete data points 32 . If there are no detectable responses in the profile, the potency is declared "undefined". However, if potency cannot be estimated within the observed response profile but a detectable response is found within data region, finite difference calculus is used to predict the assay response beyond the tested concentration range. This extrapolation continues until POD WES is quantitatively estimated or designated "less than C 1 " for profiles that have substantial activity at the lowest observed concentration but no quantitative estimate is obtainable (see Fig. 1). No data points located within the detection window are extrapolated outside of the noise region. The estimation of POD WES is described in greater detail in the Supplementary Information. Evaluating potency estimates. AC 50 estimates were determined according to Shockley 30 . The POD WES approach was described above and presented in Fig. 1. The precision of each potency estimator was investigated by calculating the empirical 95% confidence interval widths (2.5th percentile -97.5th percentile) of the log 10 transformed estimates within a generated data set. Bias was calculated by subtracting the "true" value θ of potency estimator U (e.g., AC 50 or POD WES obtained from profiles simulated with ERROR = "0%") from the mean of the estimated values of U according to θ ∑ − = U n k n k 1 1 . Evaluations of potency estimates are performed using the log 10 transformation so that the distributions of potencies better approximate a normal distribution with constant error 33 . All computations were performed in the statistical software R 34 .