QSAR Models for Predicting Cathepsin B Inhibition by Small Molecules: Continuous and Binary QSAR Models to Classify Cathepsin B Inhibition Activities of Small Molecules

doi:10.1016/j.jmgm.2010.01.009

J Mol Graph Model. Author manuscript; available in PMC 2011 June 1.

Published in final edited form as:

J Mol Graph Model. 2010 June; 28(8): 714–727.

Published online 2010 February 1. doi: 10.1016/j.jmgm.2010.01.009

PMCID: PMC2873115

NIHMSID: NIHMS175992

QSAR Models for Predicting Cathepsin B Inhibition by Small Molecules

Continuous and Binary QSAR Models to Classify Cathepsin B Inhibition Activities of Small Molecules

Zhigang Zhou, Yanli Wang, and Stephen H. Bryant^*

National Center for Biotechnology Information, National Library of Medicine, National Institutes of Health Bethesda, MD, 20894, USA

^* Corresponding author. Tel: (301)-435–7792; Fax: (301)-480–9241; Email: zhougeor/at/ncbi.nlm.nih.gov, Email: bryant/at/ncbi.nlm.nih.gov. NCBI/NIH, Bethesda, MD 20894

The publisher's final edited version of this article is available at J Mol Graph Model

See other articles in PMC that cite the published article.

Abstract

Cathepsin B is a potential target for the development of drugs to treat several important human diseases. A number of inhibitors targeting this protein have been developed in the past several years. Recently, a group of small molecules were identified to have inhibitory activity against cathepsin B through high throughput screening (HTS) tests. In this study, traditional continuous and binary QSAR models were built to classify the biological activities of previously identified compounds and to distinguish active compounds from inactive compounds for drug development based on the calculated molecular and physicochemical properties. Strong correlations were obtained for the continuous QSAR models with regression correlation coefficients (r²) and cross-validated correlation coefficients (q²) of 0.77 and 0.61 for all compounds, and 0.82 and 0.68 for the compound set excluding 3 outliers, respectively. The models were further validated through the leave-one-out (LOO) method and the training-test set method. The binary models demonstrated a strong level of predictability in distinguishing the active compounds from inactive compounds with accuracies of 0.89 and 0.94 for active and inactive compounds, respectively, in non-cross-validated models. Similar results were obtained for the cross-validated models. Collectively, these results demonstrate the models’ ability to discriminate between active and inactive compounds, suggesting that the models may be used to pre-screen compounds to facilitate compound optimization and to design novel inhibitors for drug development.

Keywords: Binary QSAR, regression, partial least squares, docking, cathepsin B protein, screening, PubChem Bioassay

Introduction

As members of the papain superfamily,¹ cathepsins are involved in many biological processes related to human diseases and disorders.¹^–⁴ Previously identified cathepsins include cathepsin B, D, H, K, L, and S. Several of these proteins have been selected as biological targets to develop therapeutic treatments, and a number of inhibitors have been identified and developed for many of these enzymes. Cathepsin B inhibitors are highly sought after chemical agents since many diseases, such as neurodegenerative disorder, cardiovascular disease, cancer, inflammation, rheumatoid arthritis, and Alzheimer’s disease⁵^–¹², have been connected with unusual levels or abnormal function of cathepsin B. As an ubiquitous lysosomal cysteine proteases, cathepsin B has been found to be responsible for intracellular as well as extracellular proteolysis in mammalian cells and can facilitate cell migration by dissolving the extracellular barriers, which result in tumor metastasis and angiogenesis.¹³^–¹⁵ The biological activity and function of cathepsin B is also important during viral infection and replication for several viruses, such as Ebola, SARS (Severe Acute Respiratory Syndrome) in human cells.¹⁶^,¹⁷ Due to its important biological functions, which are directly related to several important human diseases, cathepsin B has been chosen as a drug development target in many efforts.

A number of compounds have been found to inhibit cathepsin B activity, and some of these compounds have been tested and are effective in animal experiments.⁵^–¹²^,¹⁸^–²⁰ Most of these cathepsin inhibitors disable the biological activity of cathepsin B through forming irreversible covalent chemical bond in the catalytic site of the enzyme. These irreversible inhibitors include dipeptidyl nitriles²¹, vinyl sulfones, expoxysuccinates, acyloxymethyl ketones, fluoromethyl ketones, hydrazides, and bis-α-amidoketones.²² Structural studies have provided detailed insights into the biological mechanisms of these inhibitors. Inhibitors bind to the catalytic active site of cathepsin B, and form an irreversible covalent bond with the protein in the active site.²¹^,²³ Meanwhile, computational studies, such as docking and virtual screening, were also used to explore the binding and inhibition mechanism of these inhibitors.²⁴^–²⁶ These computational approaches have proved efficient and essential for optimizing and designing chemical agents with improved biological activities.²⁷^–²⁹

In an effort to identify chemical probes through high throughput screening technology (HTS), a number of chemicals have been found to be active against cathepsin B in the screening campaign for inhibitor identification, with the NIH Molecular Libraries Program (MLP). The screening results were deposited into the PubChem (http://pubchem.ncbi.nlm.nih.gov), a database for molecular structures and associated biological activities which is available to general public. The database contains the 2D structures and the biological activity (IC₅₀) derived from dose-response HTS experiments for the active compounds tested.

In a previous study,³⁰^,³¹ docking simulation was used to model the binding structures of the active compounds identified by the high throughput screening (HTS) tests to the binding site of the protein. A relative binding affinity was calculated for each compound studied based on the respective modeled bound structure using the linear response molecular mechanics Poisson Boltzmann-surface area method (LR-MM-PBSA). Strong correlations between the calculated binding affinities and experimental biological activities were obtained.³⁰^,³¹ Three-dimension (3D) Comparative Molecular Field Analysis (CoMFA) quantitative structure-activity relationships (QSAR) models were also established based on the multi-conformation method with high correlation coefficients for these active compounds (to be published). Given the importance of the enzyme in disease treatment and the broad research interest to screen and identify potent inhibitors, efforts have been taken further to build the QSAR models to correlate the molecular and physicochemical properties with the biological activities of the compounds tested in the Cathepsin B inhibition screening experiments. Efforts were also undertaken to build statistical models to classify inhibition activities against Cathepsin B for specific compounds. This work reports the results of the continuous (continuous in bioactivity space) and binary (two discrete bioactivity status, e.g. ‘active’ vs ‘inactive’) QSAR models in order to obtain insight into the relationship of chemical structure and physicochemical properties/biological activity. Additionally, the results of this work may provide means to pre-screen compounds to facilitate in silico design of de novo cathepsin B inhibitors and the development of drugs for the treatment of varies human diseases using some other molecular targets.

Materials and Methods

Data Set and Molecule Preparation

The 3D chemical structures of all small molecule compounds were first converted from the 2D molecular structures obtained from the PubChem database (PubChem bioAssay accession: 820 and 523),³² which were subsequently subjected to energy minimization calculation until the root mean square deviation (RMSD) of potential energy was smaller than 0.001, using the Optimized Potentials for Liquid Simulations (OPLS) force field³³^–³⁶ by the Molecular Operating Environment (MOE) program (Version 2007.09, developed by Chemical Computing Group, Montreal, Canada). In Bioassay AID (Pubchem bioassay accession) 820, 75 compounds were tested by a dose response confirmatory screening experiment, whereas 37 compounds were confirmed to have inhibitory activity, 35 compounds were confirmed to have no inhibitory activity, and the bioactivity outcome for the remaining 3 compounds were reported as “inconclusive”. By further examining these active compounds using different experimental protocol, assay depositors suggested that some of them were likely artifacts. Further investigation using several independent bioassay tests (different biological experiments), which may be affected by similar causes of artifact, suggested that five compounds bear greater chance as being false positives, and they were removed from the current analysis. So, 32 active (37–5) and 35 inactive compounds from bioassay 820 were included in the work. Bioassay 523 employed a similar experimental protocol to the Bioassay 820, where 27 compounds were tested, 10 compounds were confirmed having inhibitory activity, 16 compounds were confirmed having no inhibitory activity, and one (1) compound was reported as “inconclusive”. Four out of the 10 active compounds are unique structures and were added to the active compound list. Total 36 active compounds (32 from bioassay 820 and 4 from bioassay 523) and 51 inactive compounds (35 from bioassay 820 and 16 from bioassay 523) were used in the modeling. The average IC₅₀ reported in the bioassay depositions was used for the bioactivity of the active compounds. The PubChem compound accession (CID), molecular formula, and biological activity (IC₅₀) for the 86 compounds are listed in Table 1 (one active compound was not included in QSAR modeling as no 3D conformer was obtained in docking simulation). Molecular structure and other information can be obtained from PubChem website based on CIDs (http://pubchem.ncbi.nlm.nih.gov) and the bioassay protocol information can be obtained from PubChem website by BioAssay ID (AID) 820 and 523 (http://pubchem.ncbi.nlm.nih.gov/assay/assay.cgi?aid=820 and http://pubchem.ncbi.nlm.nih.gov/assay/assay.cgi?aid=523).³² The compounds tested in the confirmatory assays were cherry picked based on the screening results for over 60,000 compounds in a single dose HTS assay for potential cathepsin B inhibition activity (PubChem AID:453) by measuring the release of the fluorophore aminomethyl coumarin (AMC) from the hydrolysis of an AMC-labeled dipeptide, The large number of inactive compounds from the primary HTS assay were not included in the model since the statistical methods typically do not perform well on unbalanced data sets. Therefore claiming all compounds as inactive will yield an accuracy of 99.94% as pointed out by Weis and coworkers.³⁷

Table 1

The compound structure and biological activity of the active and inactive compounds from PubChem bioassay database.¹

The three dimensional structure coordinate of cathepsin protein with bound Dipeptidyl Nitrile (DPN) was obtained from Protein Databank (PDB code: 1GMY²¹). The PDB file for the crystallographic structure complex of cathepsin B contains three chains A, B and C. The first two chains form a dimer and were used in docking simulation (Figure 1). Hydrogen atoms and partial charges were added to the protein and then short steps of minimization were performed to relax the newly-added hydrogen atoms and the potential steric contacts in the original PDB coordinates. The minimizations were performed using an OPLS force field following the standard protocol of the Glide program (version 8 in FirstDiscovery suite).³⁸^,³⁹ The minimized protein structure was then used to generate a docking grid which was used to docked all active compounds into the DPN binding site. A docking box (exterior box) of 16 × 18 × 22 Å for the placement of all ligand atoms and a restraining box (interior box) of 8 × 10 × 14 Å for the placement of ligand geometric center were used to generate the docking grid. The boxes were centered at the geometric center of DPN ligand. Default values were used for all other settings within Glide (version 8).

Figure 1

The active site located in the dimer interface of cathepsin protein complexed with ligand DNP (stick-ball mode) in crystal complex (PDB code: 1GMY). A: The dimer proteins are shown with surface mode with the primary protein in red and the second protein (more ...)

Conformation Determination and Molecule Alignment

All active compounds were flexibly docked into the active site of the cathepsin protein based on the generated docking grid. All docking simulations were carried out using Glide program. Details of the docking simulation can be obtained from previously reported work.³⁰^,³¹^,⁴⁰ Docking simulation was used to determine the “active” conformation of each compound according to the structural and physicochemical properties of the binding site. The 3D conformation for each inactive compound was modeled by flexibly aligning it to the 3D conformer of DPN, the ligand in the crystal complex. The alignment was carried out based on molecular shape (volume) and pharmacophore elements using the flexible alignment module implemented in MOE, Default parameters were used for the flexible alignment. The conformation for each active compound was also obtained using the same alignment for the comparison with the conformation obtained from docking simulations in QSAR modeling.

Statistical Analysis

Partial Least-Squares analysis (PLS) method was used to conduct statistical analysis and to derive a continuous QSAR model based on the descriptors calculated for each inhibitor. The PLS method has been proved to be a priority method for such statistical analysis especially in the case where the number of descriptors is large. All PLS calculations were carried out with the module implemented in MOE. The log₁₀IC₅₀ value was used as a dependent variable and the calculated molecular and physicochemical properties were used as independent variables for the analysis. Over 50 descriptors were calculated and preliminary factor analysis was conducted to obtain the optimum correlation between each individual descriptor and the biological activity (logIC₅₀). The nine descriptors listed in Table 2 were selected and used to build QSAR models. An optimal number of components (9 components for continuous models and 8 for binary models) were used to build final models.

Table 2

The calculated molecular and physicochemical properties (descriptors) of 86 compounds used to build QSAR models

Model validation was carried out by the Leave-one-out (LOO) cross-validation procedure and training-test set approach in which the compounds were split into training sets for model building and test sets for which activities were predicted by the built model. Both continuous and binary models were validated by LOO and three (3) training-test experiments in which nine compounds were randomly selected as test compounds and all other compounds were used as training set for QSAR model building. The LOO technique provides a good way to quantitatively evaluate the predictive ability and robustness of a model by predicting each compound s activity using a QSAR model built based on information of the remaining compounds, which avoids the effect of a compound on its own activity prediction. This approach tries to eliminate the over-fitting problem normally existing in conventional regression method which includes all compounds in model construction. A cross-validated correlation coefficient (q²) was used to measure the predictability of a model and the conventional correlation coefficient (r²) was used to measure the quality of a model.

Results and Discussion

Conformation determination and alignment construction

In order to calculate the physicochemical and molecular property-based descriptors to build QSAR models, 3D confirmation for each studied compound was needed. The 3D structures of these compounds were first converted from their 2D structures. As stated in the method section, totally 36 active compounds and 51 inactive compounds with unique structures were obtained from two confirmatory bioassay data entries (PubChem AID 820 and 523). Docking simulations and molecular alignment methods were applied to obtain active 3D confirmers for the 36 active compounds. In the docking simulations, the active compounds were docked into the active site of cathepsin B using the same methodology and confirmation determination procedure as described previously.³⁰ Reliable docked poses were obtained for all but one compound. Detail results and discussion on conformer selection and comparison for most compounds can be found in previous report.³⁰ As a result, the 35 compounds with docking results were included in the QSAR modeling. The docked cluster of the active compounds is depicted in Figure 2. It is seen that all these compounds were docked in the binding site with comparable binding conformers and roughly the similar locations.

Figure 2

The cluster of docked poses of the active compounds in the protein binding site.

In the molecular alignment, all 86 (active and inactive) compounds were superposed onto DPN using flexible molecular alignment method. The best fitted conformation of each compound was selected. The fitting was scored based on molecular shape (occupied volume) and pharmacophore matching which includes hydrogen donor, hydrogen acceptor, hydrophobic atoms, polar hydrogen atoms, and aromaticity of the two molecules. The 3D conformation for each compound was selected based on the best fit from MOE. Two separate sets of molecular structures were prepared, one based on docked structures for active compounds and molecular alignment for inactive compounds and the other based on the molecular alignment for all compounds. The two sets of 3D structures were used to build QSAR models for comparison and they yielded the very similar results. The information, including PubChem compound accession, e.g. CID, molecular formula, SMILES (Simplified Molecular Input Line Entry Specification) string, molecular weight (MW), and the activity (IC₅₀), as obtained from PubChem BioAssay database (AID: 820 and 523) of all 86 compounds used in the work are listed in Table 1. To build a binary QSAR model, the activities of the active and the inactive compounds were assigned as of 1 and 0, respectively.

Descriptor calculation

Over 250 descriptors are available with MOE QSAR module. Initial selection for the descriptors was attempted empirically based on the nature of the descriptors, the features of the studied molecular set, and the problem studied. Those which are not relevant to the modeling, such as the total energy, heat of formation, the numbers of bromine, boron, fluorine, or phosphorus atoms, were eliminated from further consideration in the work. Over 50 descriptors of physicochemical molecular properties were selected and calculated for all studied compounds based on the constructed 3D conformations using the methods implemented in MOE (version 2007.09). Factor analyses were performed to examine the correlation between each individual descriptor and the biological activity (logIC₅₀) of these compounds. The nine descriptors that demonstrated an apparent correlation with the observed biological activity were chosen to build QSAR models, which are listed in Table 2. Among these descriptors, N_acc+don is the sum of the counts of hydrogen acceptor and donor. Apol, MR, and Dipole are the sum of atomic polarizabilities,⁴¹ molar refractivity, and the dipole moment of a molecule calculated from the partial charges of the molecule, respectively. SA and SA_pol are the water accessible surface areas (WASA) of whole molecule and the hydrophilic part of WASA calculated using a radius of 1.4 Å for the probe, respectively. Vol, E_sol, and logP_(o/w) are the molecular volume, the empirical solvation energy calculated based on OPLS force field, and the log value of octanol/water partition coefficient of a compound, respectively. Most of these descriptors are related to the molecular solvation properties to some degree. Calculation of molecular solvation or partition of a compound between hydrophobic and hydrophilic compartments still remains challenging. Combination of various computational methods would be a better choice to estimate the chemical s partition property which is well regarded to have direct and large effect on its biological activities.

Continuous QSAR models

Conventional QSAR models were constructed for the 35 active compounds based on the measured biological activities (logIC₅₀) using the PLS regression method. A model for the entire compound set was first built to examine model coherence of all compounds and to identify the potential outliers. The model was then validated using the LOO method and three (3) training-test sets to examine the predictability and robustness of the QSAR models. To extend the prediction spectrum to include inactive compounds, seven (7) inactive compounds were randomly chosen and included in the models. The inhibitory activity (reading) of these compounds at the max tested concentration (50 μM) was close to the reading of the control experiments. It was assumed that their IC₅₀ (if exists) was much high than 50 μM. To enable the numerical analysis, a value of −3.00 was assigned as the IC₅₀ of these five inactive compounds in the QSAR analysis.

The predicted results and prediction errors for the model based on all 42 compounds (All Model) are listed in Table 3. Two thresholds, one and half units, respectively were used to validate the calculated results. The former corresponds to a ± 10 fold variation ranges or with a range of a single order of magnitude on each direction from the experimental value of IC₅₀ and the latter corresponds to ± 3.33 fold variation ranges or within a window of a single order of magnitude on IC₅₀. The non-validated model based on 9 descriptors (properties) showed strong ability to predict the biological activities with errors (difference between the predicted value and experimental value) of smaller than one unit for all 42 compounds. Twenty six compounds (62%) have predicted errors smaller than or equal to 0.5 unit by the model. Model evaluation using the LOO validation method generated satisfactory results. The observed activities vs. the predicted activities by the non-validated fitting model and cross-validated model are plotted in Figure 3A. By further checking the predicted results of the cross-validation model, it was noted that the three compounds (compounds 24, 26 and 34) had predicted errors close to one unit. When these three compounds were treated as outliers and excluded from the model construction, the new model (3 Outlier Model) consisting of the remaining 39 compounds demonstrated an improved quality. The predicted activities for these 39 compounds are listed in Table 4. All compounds have predicted errors smaller than one unit in the non-validated model and 37 compounds have predicted errors smaller than on unit in the cross-validated model. Twenty seven and more than half compounds have predicted errors smaller than half (0.5) unit according to the non-validated and cross-validated results, respectively. These cross-validation results demonstrated the strong predictability of the constructed models for the application of cathepsin B inhibition activity prediction.

Table 3

Continuous All Model: The predicted activities and errors (difference between predicted and observed activity (logIC₅₀)) of 42 compounds predicted by two statistic methods based on the QSAR model.

Figure 3

The plots of the predicted vs. observed bioactivity (logIC₅₀) for the continuous models. A: The Continuous All Model. B: The training-test Model 1. C: The training-test Model 2. D: The training-test Model 3. For B–D: The predicted activity for (more ...)

Table 4

3 Outlier Model: The predicted activities and residues of 37 compounds predicted by two statistic methods based on the QSAR model excluding three outliers (compounds 24, 26, and 34) from the All Model.

To further evaluate the model s predictability and robustness, three additional models were built using the training-test-set method, in which the compounds were randomly split into a training set and test set. In this work, six compounds were randomly selected for the test set and this process was repeated three times. The predicted results for the three models are listed in Table 5–7. The plots of the observed activities vs. the predicted activities for the training compounds and testing compounds are shown in Figure 3B--2D.2D. It is seen that the predicted results for the testing compounds are close to their experimental activities with reasonable deviations. All three models based on the training set of compounds demonstrated similar predictability results.

Table 5

The Training-test Model 1: The predicted activities and residues of 36 compounds in the training set and 6 compounds in the test set.

Table 7

The Training-test Model 3: The predicted activities and residues of 36 compounds in the training set and 6 compounds in the test set.

As a summary, the non-validated regression coefficient, cross-validated coefficient, and root mean square error (RMSE) for the All Model are 0.771, 0.608, and 0.483, respectively (Table 8). The non-validated, cross-validated coefficients, and RMSE for the 3 Outlier Model are 0.821, 0.683, and 0.424, respectively (Table 8). The statistics results for the three training-test models are also listed in Table 8. All three models based on the training set of compounds have similar statistical properties as the All Model described above, with non-validated correlation coefficients of 0.804–0.807 and cross-validated coefficients of 0.637–0.675. The predicted errors for all test compounds but one (compound 34) in the three training sets are smaller than one unit. The Model 1 and 3 (Table 5 and and7)7) has two test compounds predicted with errors larger than half unit and the Model 2 has the most compounds predicted with errors larger than half unit. These results show that the predictability of these models is acceptable. Considering the fact that the compound structures and bioactivities studied are reasonably diverse, it is reasonable to conclude that the models demonstrated relatively strong prediction power for estimating the activity of the “unknown” compounds based on the known” compounds.

Table 8

Summary of statistics and relative importance of descriptors for three fitting and 3 training-test QSAR models.

Binary QSAR models to classify active and inactive compounds

The availability of confirmed active and inactive compounds provides an opportunity to develop binary models for classifying active cathepsin B inhibitors from inactive compounds. Although modern HTS technology allows rapid screening of tens of thousands of compounds in a matter of a few days, such experiments are expensive, and can hardly cover all chemical space. In silico screening helps to eliminate potentially inactive compounds and to build compound libraries with enriched promising lead candidates. Thus in silico screening is a cost efficient strategy to facilitate the identification of potential drug candidates. Such prescreening provides a way to eliminate potentially inactive molecules in novel compound design for synthesis and make the limited experimental resources available for screening only the potentially active molecules.

The same nine descriptors used in the continuous models were also used to build binary QSAR models relying on the method developed by the Chemical Computing Group which was implemented as a module in MOE.⁴²^,⁴³ The binary QSAR model was developed as an economic QSAR model to be used distinguish active vs. inactive compounds in HTS experiments to assist model drug development. All binary models were built using the binary method module of MOE.⁴²^,⁴³ Similarly to the process of continuous model construction, the “All Binary Model” was first built with all 86 compounds to examine the predictability of the model for the “inside” compounds, the compounds used to build a model, and followed by validations using the LOO method and three training-test sets to examine the predictability and robustness of the QSAR models for “outside” compounds. The predicted results for the All Binary Model are listed in Table 9 and the statistical results are listed in Table 13. All 86 compounds, except for four active and three inactive ones, were predicted correctly with this non-validated model (Table 9). The prediction accuracies (Table 13) for all compounds are 0.919, and for active and inactive compounds are 0.886 and 0.941 respectively. This demonstrates the strong classification power of this model for discriminating inactive compounds from active compounds. In addition, four cross-validated models were built, one from LOO and three from the training-test-set method. In the latter three models, the compounds were randomly split into training set (80 compounds) and test set (6 compounds). The four models were designated as “Binary Model”, “Bin. Train-test 1”, “Bin. Train-test 2”, and “Bin. Train-test 3”, and the predicted results from these four cross-validated models are listed in Tables 9–12. The three training-test models produced consistent results that are similar to those from the LOO-based Binary Model. For the 80 training set compounds, the non-validated prediction accuracies for all compounds are 0.906–9.25, and for active and inactive compounds are 0.906 and 0.938, respectively. The cross-validated prediction accuracies for all compounds are 0.900, active and inactive compounds are 0.844 and 0.938, respectively. For the six test compounds, five compounds were predicted correctly, the prediction accuracies for the randomly selected test compounds in all three training-test models are 0.833 (Table 13). These validated models clearly demonstrated their abilities to distinguish active vs. inactive compounds with high efficiency and accuracy, which indicates their strong potential for in silico screening of small molecules.

Table 9

Binary Model: The predicted active vs. inactive for all active and inactive compounds predicted by two statistic methods based on the binary QSAR model.

Table 13

Summary of statistics and importance of descriptors for the binary model and three training-test binary models.

Table 12

Binary Training-test Model 3: The predicted active vs. inactive for training set compounds and the predicted results for the tested compounds.

It is noted that the prediction accuracies for active compounds were slightly, but consistently, lower than those for the inactive compounds. Such phenomenon could arise from the unbalanced nature of the data set regarding to bioactivity classes (active and inactive in this work). The larger population of inactive compounds over that of the active ones shifts the splitting point used to determine the active vs. inactive to the active side. Such shift results in more active compounds at the border line of being classified as “inactive”. To solve the problem, one could design a weighting schema to “overestimate” the type of compounds that are underrepresented in the data set. In many cases, however, such unbalance would not cause a severe effect on the application of the method when the focus is to select a molecule with a given property, such as eliminating inactive molecules to further develop active compounds in designing cathepsin B inhibitors.

Conclusion

The availability of molecular structures and their inhibition activity against cathepsin B identified through HTS experiments provided an opportunity to conduct a structure-activity relationship study to facilitate chemical probe development, and potential drug development for related disease treatment. The “active” conformers of active compounds were modeled using docking simulations, and the conformations of the inactive compounds were determined based on flexible alignment of each compound to the 3D conformer of the original ligand (DPN) in the crystal structure complex. Conventional continuous models were built for predicting cathepsin B inhibition activity of active small molecules. Binary models were constructed for the “active” vs. “inactive” classification. Nine (9) molecular and physicochemical properties were calculated based on 3D conformational information, and were selected and used as descriptors to build the QSAR models.

The continuous models demonstrated reasonable correlations between the predicted and the observed activities with non-validated (r²) and cross-validated (q²) regression correlation coefficients of 0.771 and 0.608 for all compounds, and 0.821 and 0.683 for the compounds excluding 3 outliers, respectively. The model showed strong predictability with reasonable predicted error rates over the entire test set and across several models. Almost all the 42 compounds used in the model have predicted errors smaller than one unit, while the majority of the compounds have predicted errors less than a half unit from both the non-validated and cross-validated results. The predictability of the models were further validated using the three training-test methods with the results that the predicted errors for all test compounds are smaller than one unit, while a significant fraction of them are less than a half unit in all three models together.

The consistent performance of the four cross-validated models demonstrated that the models are robust in predicting the inhibition activities of small molecules for the biological system of cathepsin B. The success of the QSAR models proves that the modeling approach based on docking conformation determination and the nine physicochemical properties may be a promising method to predict biological activity of a small molecule for drug development based on the knowledge derived from the identified cathepsin B inhibitors.

Furthermore, binary models were built for classifying cathepsin B inhibition activities, e.g. active vs. inactive categories. The statistical results showed that the models have demonstrated excellent performance with the prediction accuracies for bioactivity classes of 0.919, 0.886, and 0.941 for all 86 tested, 35 active, and 51 inactive compounds in the non-validated binary models, respectively. All compounds, except for four active and three inactive compounds, were correctly predicted by the non-validated fitting model. Similar prediction results were obtained with the cross-validated model. With the three training-test models, the bioactivity categories for the five out of six tested compounds were predicted correctly, with a prediction accuracy of 0.833. The factor that all of these models yielded high accuracies upon the prediction of cathepsin B inhibition activity demonstrates that such models enable the classification of active vs. inactive inhibitory compounds with high accuracies. The robustness and feasibility of this binary method indicates its potential in several applications to facilitate molecular design and drug development. In particular, this approach may prove highly efficient for pre-screening in silico designed compounds to optimize the biological activities based on the information of known bioactive compounds. Such in silico pre-screening may be performed to filter out the ones that are unlikely to have the desired biological activities, thus to allow the limited efforts to be focused on the highly promising candidates. With such advantages, we anticipate the utilization of the combination of the continuous and binary classification models will potentially benefit modern drug development for therapeutic purposes.

Table 6

The Training-test Model 2: The predicted activities and residues of 36 compounds in the training set and 6 compounds in the test set.

Table 10

Binary Training-test Model 1: The predicted active vs. inactive for training set compounds and the predicted results for the tested compounds.

Table 11

Binary Training-test Model 2: The predicted active vs. inactive for training set compounds and the predicted results for the tested compounds.

Acknowledgments

This research was supported by the Intramural Research Program of the NIH, NLM. The authors would like to thank the NIH Fellows Editorial Board (FEB) for reviewing the manuscript and the developers of Pymol software for sharing the program to prepare the molecular figures used in this paper.

References

1. Rooprai HK, McCormick D. Proteases and their inhibitors in human brain tumours: a review. Anticancer Res. 1997;17:4151–4162. [PubMed]

2. Barrett AJ, Kirschke H. Cathepsin B, Cathepsin H, and cathepsin L. Methods Enzymol. 1981;80 Pt C:535–561. [PubMed]

3. Chapman HA, Jr, Munger JS, Shi GP. The role of thiol proteases in tissue injury and remodeling. Am J Respir Crit Care Med. 1994;150:S155–159. [PubMed]

4. Roberts LR, Adjei PN, Gores GJ. Cathepsins as effector proteases in hepatocyte apoptosis. Cell Biochem Biophys. 1999;30:71–88. [PubMed]

5. Giusti I, D’Ascenzo S, Millimaggi D, Taraboletti G, Carta G, Franceschini N, Pavan A, Dolo V. Cathepsin B mediates the pH-dependent proinvasive activity of tumor-shed microvesicles. Neoplasia. 2008;10:481–488. [PMC free article] [PubMed]

6. Gounaris E, Tung CH, Restaino C, Maehr R, Kohler R, Joyce JA, Plough HL, Barrett TA, Weissleder R, Khazaie K. Live imaging of cysteine-cathepsin activity reveals dynamics of focal inflammation, angiogenesis, and polyp growth. PLoS ONE. 2008;3:e2916. [PMC free article] [PubMed]

7. Ha SD, Martins A, Khazaie K, Han J, Chan BM, Kim SO. Cathepsin B is involved in the trafficking of TNF-alpha-containing vesicles to the plasma membrane in macrophages. J Immunol. 2008;181:690–697. [PubMed]

8. Haque A, Banik NL, Ray SK. New insights into the roles of endolysosomal cathepsins in the pathogenesis of Alzheimer’s disease: cathepsin inhibitors as potential therapeutics. CNS Neurol Disord Drug Targets. 2008;7:270–277. [PubMed]

9. Sandes E, Lodillinsky C, Cwirenbaum R, Arguelles C, Casabe A, Eijan AM. Cathepsin B is involved in the apoptosis intrinsic pathway induced by Bacillus Calmette-Guerin in transitional cancer cell lines. Int J Mol Med. 2007;20:823–828. [PubMed]

10. Lutgens SP, Cleutjens KB, Daemen MJ, Heeneman S. Cathepsin cysteine proteases in cardiovascular disease. Faseb J. 2007;21:3029–3041. [PubMed]

11. Hook V, Kindy M, Hook G. Cysteine protease inhibitors effectively reduce in vivo levels of brain beta-amyloid related to Alzheimer’s disease. Biol Chem. 2007;388:247–252. [PubMed]

12. Downs LS, Jr, Lima PH, Bliss RL, Blomquist CH. Cathepsins B and D activity and activity ratios in normal ovaries, benign ovarian neoplasms, and epithelial ovarian cancer. J Soc Gynecol Investig. 2005;12:539–544. [PubMed]

13. Vasiljeva O, Korovin M, Gajda M, Brodoefel H, Bojic L, Kruger A, Schurigt U, Sevenich L, Turk B, Peters C, Reinheckel T. Reduced tumour cell proliferation and delayed development of high-grade mammary carcinomas in cathepsin B-deficient mice. Oncogene. 2008;27:4191–4199. [PubMed]

14. Vasiljeva O, Papazoglou A, Kruger A, Brodoefel H, Korovin M, Deussing J, Augustin N, Nielsen BS, Almholt K, Bogyo M, Peters C, Reinheckel T. Tumor cell-derived and macrophage-derived cathepsin B promotes progression and lung metastasis of mammary cancer. Cancer Res. 2006;66:5242–5250. [PubMed]

15. Jane DT, Morvay L, Dasilva L, Cavallo-Medved D, Sloane BF, Dufresne MJ. Cathepsin B localizes to plasma membrane caveolae of differentiating myoblasts and is secreted in an active form at physiological pH. Biol Chem. 2006;387:223–234. [PubMed]

16. Chandran K, Sullivan NJ, Felbor U, Whelan SP, Cunningham JM. Endosomal Proteolysis of the Ebola Virus Glycoprotein Is Necessary for Infection. Science. 2005;308:1643–1645. [PubMed]

17. Simmons G, Gosalia DN, Rennekamp AJ, Reeves JD, Diamond SL, Bates P. Inhibitors of cathepsin L prevent severe acute respiratory syndrome coronavirus entry. Proc Natl Acad Sci USA. 2005;102:11876–11881. [PubMed]

18. Hook VY, Kindy M, Hook G. Inhibitors of cathepsin B improve memory and reduce beta-amyloid in transgenic Alzheimer disease mice expressing the wild-type, but not the Swedish mutant, beta-secretase site of the amyloid precursor protein. J Biol Chem. 2008;283:7745–7753. [PubMed]

19. Hosokawa M, Kashiwaya K, Eguchi H, Ohigashi H, Ishikawa O, Furihata M, Shinomura Y, Imai K, Nakamura Y, Nakagawa H. Over-expression of cysteine proteinase inhibitor cystatin 6 promotes pancreatic cancer growth. Cancer Sci. 2008;99:1626–1632. [PubMed]

20. Parker BS, Ciocca DR, Bidwell BN, Gago FE, Fanelli MA, George J, Slavin JL, Moller A, Steel R, Pouliot N, Eckhardt BL, Henderson MA, Anderson RL. Primary tumour expression of the cysteine cathepsin inhibitor Stefin A inhibits distant metastasis in breast cancer. J Pathol. 2008;214:337–346. [PubMed]

21. Greenspan PD, Clark KL, Tommasi RA, Cowen SD, McQuire LW, Farley DL, van Duzer JH, Goldberg RL, Zhou H, Du Z, Fitt JJ, Coppa DE, Fang Z, Macchia W, Zhu L, Capparelli MP, Goldstein R, Wigg AM, Doughty JR, Bohacek RS, Knap AK. Identification of dipeptidyl nitriles as potent and selective inhibitors of cathepsin B through structure-based drug design. J Med Chem. 2001;44:4524–4534. [PubMed]

22. Otto HH, Schirmeister T. Cysteine Proteases and Their Inhibitors. Chem Rev. 1997;97:133–172. [PubMed]

23. Yamamoto A, Tomoo K, Hara T, Murata M, Kitamura K, Ishida T. Substrate specificity of bovine cathepsin B and its inhibition by CA074, based on crystal structure refinement of the complex. J Biochem. 2000;127:635–643. [PubMed]

24. Mladenovic M, Ansorg K, Fink RF, Thiel W, Schirmeister T, Engels B. Atomistic Insights into the Inhibition of Cysteine Proteases: First QM/MM Calculations Clarifying the Stereoselectivity of Epoxide-Based Inhibitors. J Phys Chem B. 2008;112:11798–11808. [PubMed]

25. Redzynia I, Ljunggren A, Abrahamson M, Mort JS, Krupa JC, Jaskolski M, Bujacz G. Displacement of the occluding loop by the parasite protein, chagasin, results in efficient inhibition of human cathepsin B. J Biol Chem. 2008;283:22815–22825. [PubMed]

26. Watanabe D, Yamamoto A, Tomoo K, Matsumoto K, Murata M, Kitamura K, Ishida T. Quantitative evaluation of each catalytic subsite of cathepsin B for inhibitory activity based on inhibitory activity-binding mode relationship of epoxysuccinyl inhibitors by X-ray crystal structure analyses of complexes. J Mol Biol. 2006;362:979–993. [PubMed]

27. Markt P, McGoohan C, Walker B, Kirchmair J, Feldmann C, Martino GD, Spitzer G, Distinto S, Schuster D, Wolber G, Laggner C, Langer T. Discovery of Novel Cathepsin S Inhibitors by Pharmacophore-Based Virtual High-Throughput Screening. J Chem Inf Model. 2008;48:1693–1705. [PubMed]

28. Beavers MP, Myers MC, Shah PP, Purvis JE, Diamond SL, Cooperman BS, Huryn DM, Smith AB., 3rd Molecular docking of cathepsin L inhibitors in the binding site of papain. J Chem Inf Model. 2008;48:1464–1472. [PMC free article] [PubMed]

29. Shah PP, Myers MC, Beavers MP, Purvis JE, Jing H, Grieser HJ, Sharlow ER, Napper AD, Huryn DM, Cooperman BS, Amos I, Smith B, Diamond SL. Kinetic Characterization and Molecular Docking of a Novel, Potent, and Selective Slow-Binding Inhibitor of Human Cathepsin L. Mol Pharmacol. 2008;74:34–41. [PMC free article] [PubMed]

30. Zhou Z, Wang Y, Bryant SH. Computational Analysis of the Cathepsin B Inhibitors Activities Through LR-MMPBSA Binding Affinity Calculation Based on Docked Complex. Journal of Computational Chemistry. 2009;30:2165–2175. [PMC free article] [PubMed]

31. Zhou Z, Bates M, Madura JD. Structure modeling, ligand binding, and binding affinity calculation (LR-MM-PBSA) of human heparanase for inhibition and drug design. Proteins: Structure, Function, and Bioinformatics. 2006;65:580–592. [PubMed]

32. Myers MC, Napper AD, Motlekar N, Shah PP, Chiu C-H, Beavers MP, Diamond SL, Huryn DM, ABS Identification and characterization of 3-substituted pyrazolyl esters as alternate substrates for cathepsin B: the confounding effects of DTT and cysteine in biological assays. Bioorg Med Chem Lett. 2007;17:4761–4766. [PMC free article] [PubMed]

33. Jorgensen WL, Tirado-Rives J. The OPLS Potential function for proteins Energy minimizations for crystals of cyclic peptides and crambin. J Am Chem Soc. 1988;110:1657–1666.

34. Damm W, Frontera A, Tirado-Rives J, Jorgensen WL. OPLS all-atom force field for carbohydrates. Journal of Computational Chemistry. 1997;18:1955–1970.

35. Rizzo RC, Jorgensen WL. OPLS All-Atom Model for Amines: Resolution of the Amine Hydration Problem. Journal of the American Chemical Society. 1999;121:4827–4836.

36. Jorgensen WL, Maxwell DS, Tirado-Rives J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. Journal of the American Chemical Society. 1996;118:11225–11236.

37. Weis DC, DPV, Faulon J-L. Data mining PubChem using a support vector machine with the Signature molecular descriptor: Classification of factor XIa inhibitors. Journal of Molecular Graphics and Modelling. 2008;27:466–475. [PubMed]

38. Friesner RA, Banks JL, Murphy RB, Halgren TA, Klicic JJ, Mainz DT, Repasky MP, Knoll EH, Shelley M, Perry JK, Shaw DE, Francis P, Shenkin PS. Glide: A new approach for rapid, accurate docking and scoring. 1. method and assessment of docking accuracy. Journal of Medicinal Chemistry. 2004;47:1739–1749. [PubMed]

39. Halgren TA, Murphy RB, Friesner RA, Beard HS, Frye LL, Pollard WT, Banks JL. Glide: A new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening. Journal of Medicinal Chemistry. 2004;47:1750–1759. [PubMed]

40. Zhou Z, Khaliq M, Suk JE, Patkar C, Li L, Kuhn RJ, Post CB. Antiviral Compounds Discovered by Virtual Screening of Small Molecule Libraries against Dengue Virus E Protein. ACS Chemical Biology. 2008;3:765–775. [PMC free article] [PubMed]

41. Lide DR. CRC Handbook of Chemistry and Physics. CRC Press; Boca Raton, FL: 1994.

42. Labute P. Binary QSAR: a new method for the determination of quantitative structure activity relationships. Pacific Symposium on Biocomputing. 1999;4:444–455. [PubMed]

43. Gao H, Williams C, Labute P, Bajorath J. Binary Quantitative Structure-Activity Relationship (QSAR) Analysis of Estrogen Receptor Ligands. Journal of Chemical Information and Computer Science. 1999;39:164–168. [PubMed]