PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Chem Inf Model. Author manuscript; available in PMC 2010 June 1.
Published in final edited form as:
PMCID: PMC2702476
NIHMSID: NIHMS118618

Carborane Clusters in Computational Drug Design: A Comparative Docking Evaluation Using Autodock, Flexx, Glide and Surflex#

Abstract

Compounds containing boron atoms play increasingly important roles in the therapy and diagnosis of various diseases, particularly cancer. However, computational drug design of boron-containing therapeutics and diagnostics is hampered by the fact that many software packages used for this purpose lack parameters for all or part of the various types of boron atoms. In the present paper, we describe simple and efficient strategies to overcome this problem, which are based on the replacement of the boron atom types with carbon atom types. The developed methods were validated by docking closo- and nido-carboranyl antifolates into the active site of a human dihydrofolate reductase (hDHFR) using AutoDock, Glide, FlexX, and Surflex and comparing the obtained docking poses with the poses of their counterparts in the original hDHFR-carboranyl antifolate crystal structures. Under optimized conditions, AutoDock and Glide were equally good in docking of the closo-carboranyl antifolates followed by Surflex and FlexX whereas Autodock, Glide, and Surflex proved to be comparably efficient in the docking of nido-carboranyl antifolates followed by FlexX. Differences in geometries and partial atom charges in the structures of the carboranyl antifolates resulting from different data sources and/or optimization methods did not impact the docking performances of AutoDock or Glide significantly. Scoring functions generated by all four programs were in accordance with experimental data.

Keywords: Carboranes, Carboranyl Antifolates, Docking, Dihydrofolate Reductase, AutoDock, FlexX, Glide, Surflex

Introduction

Chemical, physicochemical, and structural versatility combined with high stability under physiological conditions are distinctive features of carboranes and other boron clusters.1,2 They have been used for decades in the design and synthesis of therapeutics for Boron Neutron Capture Therapy (BNCT) of cancer and, more recently, in other areas of drug design.1,2 Hydrophobic closo-carboranes, comparable in dimensions to adamantane,3 were used as bioisosteric replacements for (hetero)aromatic and (hetero)aliphatic ring systems and other bulky entities in the design and synthesis of carboranyl derivatives of various amino acids and peptides,2 estrogen receptor modulators,4 androgen receptor antagonists,5 retinoids,3 benzolactamic protein kinase C inhibitors,6 thalidomide,7 flufenamic acid,8 diflunisal,8 thrombin inhibitors,9 and trimethoprim.10 Many of these boronated derivatives displayed biological activities comparable or even superior to those of their non-boronated counterparts. In a different type of application, negatively charged boron clusters have been utilized as “prosthetic groups” for radiohalogens in the design and synthesis of radiotherapeutics and imaging agents.11,12 This type of cages is readily halogenated and the boron-halogen bonds within these structures appear to be less susceptible to in vivo cleavage than carbon-halogen bonds.11,12 In addition, metallocarboranes were found to be effective inhibitors of HIV-1 protease,13 and therapeutics containing single boron atoms, such as the proteasome inhibitor bortezomib,14 have attracted considerable attention in recent years.

Drug design involving boron-containing agents has two major disadvantages compared with conventional drug design: (1) There is a lack of compound libraries containing boron agents for virtual screening,15 and (2) many software packages available for structure/ligand-based drug design do not have inbuilt parameters for boron. Several strategies to circumvent the latter problem have been reported in recent years. These include the substitution of boron with carbon16-19,21,22 and the calculation of suitable boron parameter for specific applications.20 Other reports dealing with docking studies of boron compounds do not provide specific information on computational strategies addressing this problem.9,23,24

Crystal structures of proteins complexed with carborane-containing agents, including those of closo- and nido-carboranyl antifolates with human dihydrofolate reductase (hDHFR), have been reported recently.10,13 These provide for the first time the opportunity to verify docking strategies that have been previously developed for carboranyl compounds.17 In the present paper, we evaluate the docking of these carboranyl antifolates into the active site of the hDHFR crystal structures using the docking programs AutoDock, Glide, FlexX, and Surflex. The computationally determined docking poses are compared with the poses of the corresponding carboranyl antifolates in the original hDHFR crystal structures, and differences between the docking programs are discussed. The described computational studies will be of great interest for scientists involved in the design and development of boron-containing therapeutics and diagnostics who have had very limited access to computational tools in the past.

Methods

Docking algorithms

AutoDock 4 is based on a Lamarckian genetic algorithm (LGA) method. Basically, this program determines total interaction energies between random pairs of ligands and various selected portions of protein to determine docking poses.25,26 FlexX is a fragment based docking algorithm, which builds putative poses of the ligands using an incremental construction approach.27-29 The modeling of protein-ligand interactions and their binding energy predictions is based on the de novo design tool LUDI.28 Surflex generates putative poses for molecular fragments using a surface based molecular similarity method. This program employs the Hammerhead docking system for scoring.30,31 Glide docking uses a series of hierarchical filters to find the best possible ligand binding locations in a pre-built receptor grid space. The filters include a systematic search approach, which samples the positional, conformational, and orientational space of the ligand before evaluating the energy interactions of the ligand with the protein.32,33

Crystal Structures

PDB ID # 2C2S: Human dihydrofolate reductase (hDHFR) complexed with [5-(1,2-closo-dicarbadodecarboran-1-yl)methyl]-2,4-diamino-6-methylpyrimidine (Compound 1 in Figure 1). This form of compound 1 will be referred to as 1 C-prot throughout the paper. PDB ID # 2C2T: Human DHFR complexed with a racemic mixture of [5-(7,8-nido-dicarbaundecarboran-7-yl)methyl]-2,4-diamino-6-methylpyrimidine (Compound 2 in Figure 1). This form of compound 2 will be referred to as 2 N-prot throughout the paper.

Figure 1
Trimethoprim and carboranyl antifolates 1 and 2. The positioning of the “extra hydrogen” between B9-B10 corresponds to that in 2 N-crys (see Ligand Construction and Optimization for further information)

Small molecule crystal structure data for 1 and 2 not complexed with hDHFR were generously provided by Dr. Steven Ealick, Cornell University, Ithaca, NY, and Dr. David Borhani, Harvard Medical School, Boston, MA.10 These forms of 1 and 2 will be referred to as 1 C-crys and 2 N-crys, respectively, throughout the paper.

Ligand Construction and Optimization

1 C-prot and both enantiomers of 2 N-prot were extracted from their hDHFR crystal structures. The coordinates of the “extra hydrogen”40 of 2 N-crys were inserted into the structure of 2 N-prot. Other hydrogen atoms were added to both 1 C-prot and 2 N-prot using the “add valence” option of GaussView. Mulliken-34 and APT-35 charges were calculated with Gaussian at AM1-36,37 and HF/6-31+G* levels.38,39 Espfit charges were calculated for 1 C-prot and 2 N-prot at HF/6-31+G* level (Table 1). Compounds 1 and 2 were also constructed de novo without using crystallographic coordinates with HyperChem. These forms of 1 and 2 will be referred to as 1 C-con and 2 N-con throughout the paper. In the case of 2 N-con, the nido-carboranyl moieties were aligned with the closo-carborane moiety of 1 C-prot, the coordinates for the “missing” boron atom were transferred, and then changed to hydrogen. Optimization and Mulliken charge calculations for 1 C-con and 2 N-con were carried at HF/6-31+G* level (Table 1). Following de novo construction and optimization of 2 N-con, the “extra hydrogen” of the nido-cluster positioned itself between B10 and B11, while in the case of and 2 N-crys, the “extra hydrogen” was located between B9 and B10 (see Figure 1 and Supporting Information, Figure 1SI). In solution, the “extra hydrogen” appears to be in a fluxional equilibrium between B9/B10 and B10/B11.40 Mulliken charges for 1 C-crys and 2 N-crys were calculated at HF/6-31+G* level (Table 1).

Table 1
Optimization characteristics and AutoDock docking patterns of carboranyl antifolates 1 and 2

Ligand Preparation

All Gaussian generated mol2 files of the ligands were aligned with the ligand coordinates of 1 C-prot and 2 N-prot using Sybyl before ligand preparation. AutoDock does not provide parameters that recognize boron atoms. Therefore, the boron atoms were changed to “C” by modifying the pdbqt (text) files. In the case of FlexX and Surflex, boron atoms were changed to “C.3” and “C”, respectively. FlexX automatically converts boron atoms into “Du” (dummy) atoms if they are not changed to “C.3” by the user. Surflex does not have “Du” atom parameters and it recognizes both boron and Du atoms as “funky atoms” and automatically changes their atom type to “C”. Glide v 5. (OPLS2001) does not have boron and “Du” atom parameters. Therefore all the borons were replaced with “C.3”.

AutoDock: All files generated for 1 C-prot, 2 N-prot, 1 C-crys, 2 N-crys, 1 C-con, and 2 N-con, as described under ‘Ligand Construction and Optimization’, were used for docking. All non polar hydrogens were merged before saving the file into pdbqt format. Bonds of the carborane cages that were recognized as rotatable by AutoTors in ADT were changed to non-rotatable bonds. The total number of torsions (TORSDOF) for the written output files of all ligands was set to 2. FlexX: Only 1 C-crys, 2 N-crys, 1 C-con, and 2 N-con were used for docking (see Results and Discussion for details). FlexX automatically assigned formal charges on the ligands when they were imported into the FlexX environment. For docking in “user-defined” mode, either the 2,4-diamino-5-methyl pyrimidine portion (referred to as pyrimidyl portion [PP] throughout this manuscript), the 1,2-closo-dicarbadodecarboran-1-yl portion, or 7,8-nido-dicarbaundecarboran-7-yl) portion were selected as a “base fragments”. The latter two are referred to as carboranyl portions (CPs) throughout this manuscript. Surflex: Only 1 C-crys, 2 N-crys, 1 C-con, and 2 N-con were used for docking (see Results and Discussion for details). The input files generated for FlexX docking were also used for Surflex docking. As in the case of FlexX, Surflex operates with formal charges. For docking using the fragment placement method within the Surflex environment, PPs- or CPs were prepared from 2 N-crys and 2 N-con with HyperChem and imported into Surflex. Glide: 1 C-prot, 2 N-prot, 1 C-crys, 2 N-crys, 1 C-con, and 2 N-con, as described under ‘Ligand Construction and Optimization’, were used for docking

Protein Preparation

For docking with AutoDock, FlexX, and Surflex, the protein structures in pdb format were prepared with the structure preparation tool of Sybyl. Monomers were separated from both crystal structures and the blocking groups AMI and CXC were added to the N- and C-terminis, respectively, for neutralization. Water molecules were removed, hydrogen atoms were added, and side chain amides and side chains bumps were fixed. In the case of AutoDock, pdb files along with the NADPH (Nicotinamide adenine dinucleotide phosphate) cofactor were imported into the ADT environment, atom types were assigned, and Gasteiger charges were added. For FlexX, the NADPH cofactor was removed from the pdb files and saved as a separate mol2 file before defining the active site. The receptor description files (rdf) for docking were created with Sybyl using both hDHFR crystal structures. The amino acids within a 6.5 Å radius of 1 C-prot and 2 N-prot were selected to define the active site. The pdb, rdf, and NADPH mol2 files were then imported into the FlexX environment. For Surflex, pdb files containing NADPH were imported into the Surflex environment. Protomol files 30, 31 were generated using 1 C/2 N-crys and 1 C/2 N-con separately.

In the case of Glide, protein coordinates were pre-processed for docking using the Protein Preparation Wizard provided in the Schrodinger Maestro environment. Hydrogen atoms were added and the bond orders were assigned after deleting the monomer B as well as water molecules. Assignment of protonation states was carried out followed by hydrogen bond optimization for hydroxyl groups as well as Asn, Gln, and His residues. The hydrogen atoms were then minimized with the OPLS2001 force field. Grid calculations were performed for the protein residues of the active sites with 1 C-prot and 2 N-prot coordinates as the center with default box size (Grid box center x,y,z coordinates = 2.285766, 29.832586, -3.303069; Grid box cube size in x,y,z direction = 21.816994 Å).

Docking

AutoDock: The Lamarckian genetic algorithm (LGA) was selected for the ligand conformational search. The docking area was defined using AutoGrid 4. A 40×40×40 3-D affinity grid centered around the antifolate binding site with a 0.375 Å grid point spacing was calculated for each of the following atom types: (a) protein: A (aromatic C), C, HD, N, NA, OA, SA; (b) ligand: C, A, HD, N, NA, e (electrostatic), and d (desolvation). Additional docking parameters were: Dockings trials: 100, Population size: 100, random starting position and conformation, translation step ranges: 2.0 Å, rotation step ranges: 50, elitism: 1, mutation rate: 0.02, crossover rate: 0.8, local search rate: 0.06, and energy evaluations: 250,000. Higher energy evaluations (2,500,000) or higher population sizes (250) did not alter docking performances significantly. Final docked conformations were clustered using a tolerance of 2Å root-mean-square deviations (RMSD). FlexX: Docking was performed in command line mode. The cofactor NADPH was read using a separate mol2 file as a part of the active site. Docking in the automated mode proceeded with an automated base fragment selection, which is followed by placement algorithm 3 (PA3). User-defined base fragment selection was followed by PAs 1, 2 or 327,29 and incremental build up of the entire ligand. Surflex: Default-, multistart 5-, and fragment placement (only 2 N-crys and 2 N-con) modes were explored. Glide: Final calculations were performed with the Standard-Precision (SP) rigid docking protocol. Ten thousand poses were kept in the initial phase of the docking keeping the default scoring window cutoff level 100. Ligand van der Waals radii were scaled to a factor of 0.80 for non-polar atoms with a partial charge cut-off level of 0.15 (absolute value). A maximum of 100 conformations with the best binding energies was retained for the final analysis while discarding poses with less than 0.01 Å RMSD and 0.01 Å atomic displacement as duplicates.

Docking comparisons

As stated above, the four docking programs are based on different algorithms and it may be difficult to directly compare the results obtained with each of the programs, as has been discussed previously by Cole et al.41 Apart from AutoDock, other docking programs did not generate large numbers of duplicate poses within a RMSD of 2 Å. Therefore, a comparison of the avg. RMSDs of all docked poses would be meaningless for FlexX, Surflex, and Glide. However, in the case of the latter three programs, approximately the top 25 % of all docked poses (based on the binding energies) were below 4Å of RMSD. In addition, the top ranked poses accurately reproduced the binding mode of the ligands in the crystal structure. Therefore, we used the top 25% of the poses and the top ranked poses to compare the docking results from FlexX, Surflex, and Glide.

Scoring

AutoDock, FlexX, and Glide predict the free energy of binding in terms of kcal/mol.25,28,31 Surflex predicts the binding affinity of the ligand-protein complex in the form of −log (Kd).42 In order to compare scoring features appropriately, the Surflex scoring values were converted into free energy of binding values using the following equation: kcal/mol = 0.59 loge (10 −pKd).42

Results and Discussion

The binding interactions of 1 C-prot and both enantiomers of 2 N-prot with amino acid residues in the original hDHFR crystal structures have been discussed previously by Reynolds et al.10 1 C-prot and 2 N-prot showed almost identical binding poses in the hDHFR crystal structures, as shown in Figure 3 A. The additional BH vertex of the closo-cage of 1 C-prot is positioned slightly above the center of the open faces of the nido-cages of 2 N-prot. This binding pattern is similar to that of trimethoprim in the active site of hDHFR. The carborane cages of 1 C-prot and 2 N-prot bind to the same hydrophobic pocket of hDHFR as the trimethoxyphenyl group of trimethoprim. The 2,4-diamino-5-methylpyrimidine moieties of 1 C-prot and 2 N-prot form hydrogen bonds with various amino acid residues in the active site (Glu 30, Ile 7, Val 115), as shown in Figure 2 using the example of 1 C-prot. The binding between amino acid residues of the active site with both the neutral closo-cage of 1 C-prot and the presumably negatively charged nido-cages of 2 N-prot are typical for hydrophobic interactions. Similar hydrophobic interactions were observed for a metallocarborane, consisting of two negatively charged nido-carboranes, bound to HIV protease.13 Proton–hydride type bonds43,44 may not play major roles in the interactions of the carboranyl moieties of 1 C-prot and 2 N-prot with active site amino acid residues, as is evident from the orientations of the Thr56 hydroxyl groups in the original hDHFR crystal structures, which point away from the cages (see also Figure 2).

Figure 2
Binding of 1C-prot with various active site amino acid residues of hDHFR. The 2N-prot enantiomers have similar binding interactions in the hDHFR active site. 10
Figure 3
(A) Overlap of the crystal structures poses of 1C-prot (red) and 2 N-prot (green).10 (B) Overlap of the docked pose (AutoDock) of 1C-con (magenta) with the corresponding pose of 1C-prot (cyan) in the original crystal structure (avg. RMSD: 0.666 Å) ...

An analysis of the docking studies with 1 C-crys, 2 N-crys, 1 C-prot, 2 N-prot, 1 C-con, and 2 N-con using AutoDock is shown in Table 1 and Figure 3 B-C. The number of docked poses for the closo-carboranyl structures with a RMSD below 2 Å was in all cases high (83 - 100) with minimal differences in avg. RMSDs (0.53 Å– 0.86 Å). With the exception of 1 C-con, the nido-carboranyl structures showed similar docking patterns with 93-100 poses having an RMSD below 2 Å and avg. RMSDs between 0.5 and 0.74 Å. In the case of 2 N-con, 68% of the docked poses had RMSD's lower than 2 Å while 32% had an avg. RMSD of 5.2 Å. A detailed analysis of this minor cluster revealed that the pose of the nido-o-carboranyl portion did not show any significant deviation from its counterpart in the original crystal structure while the pose of the 2,4-diaminopyrimidyl portion was altered forming a strong hydrogen bond with Asp 21 rather than with Glu 30, as in the case of crystallized ligand (Figure 3D). In the case of 2 N-con, the de novo construction and optimization resulted in a different location of the “extra hydrogen” compared with 2 N-crys and 2 N-prot (see Figure 1 and Supporting Information, Figure 1SI), which may have contributed to the occurrence of the minor cluster. In general, both enantiomers of 2 showed similar docking patterns (data not shown).

The boron-bound hydrogen atoms (1-6, 9-11) (see Figure 1 and Supporting Information, Figure 1SI) in carboranes are predicted to be slightly more electronegative than carbon-bound hydrogen42,43 and may have on average negative partial charges. According to ESPfit calculations, the average partial charges for these hydrogens are about -0.24 for two-fold negatively charged nido-carborane [C2B9H11]2-, -0.15 for one-fold negatively charged nido-carborane [C2B9H12]-, and -0.07 for neutral nido-carborane [C2B9H13].43,44 The charges generated in our own ESPfit calculations for 1 C-prot and 2 N-prot (Table 1) are consistent with those reported in the literature.43,44 However, Mulliken charges obtained for 1 C-prot and 2 N-prot were generally more positive.

Overall, the docking results obtained with AutoDock indicate that minor differences in geometries, stemming from different data sources and/or optimization methods/calculations (e.g. “extra hydrogen” position), and partial atomic charges did not have a major impact on docking accuracy although they seemed to effect the avg. binding energies to some extent (Table 1). In general the binding energies were somewhat lower for the nido-carboranyl antifolates compared with their closo-carboranyl counterparts.

In an attempt to further validate our strategy to enable docking of biomolecules containing boron clusters by replacing boron with carbon, we tried to implement known boron parameters into AutoDock before docking. Borons atoms in cage structures, such as closo-o-carborane and nido-carboranes, are hexavalent. To the best of our knowledge, no force field parameters for this boron atom type are available. So far only MM2 force field parameters for sp3 boron, as in borates, have been described.20,45 These parameters include non-bonded interaction of sp3 borons, such the van der Waals energy well depth, and the van der Waals radius. Information related to the atomic solvation parameter for sp3 borons is not available. Therefore, only partial implementation of sp3 boron parameters into AutoDock has been carried out. With respect to values related to atomic solvation and other parameters, the corresponding carbon atom type parameters were used (See Supporting Information). Overall, the docking patterns for 1 C/2 N-crys, 1 C/2 N-prot and 1 C/2 N-con following partial implementation of sp3 boron parameter were found to be comparable to those using exclusively carbon parameters (See Table 1SI in Supporting Information for details).

The Glide docking program uses ConfGen, a systematic sampling method to generate several conformations of the ligand prior to the actual docking calculation. Initially, the switch of boron atom type to “C.3” atom type using the ConfGen failed to generate conformations of the ligands because the atomic bonds in the cage structures of the ligands were considered as torsions. Therefore, we adopted a rigid docking protocol, which did not directly use the conformational sampling capability in Glide. The docking results for 1 C-crys and 2 N-crys using this approach are shown in Table 2. Overall, Glide reproduced the crystal binding modes closely as can be seen by the average RMSDs of the top 25 % poses. The average binding energy of the top 25% poses was also very close to the binding energy of the top ranked pose.

Table 2
Glide docking of 1 C-crys and 2 N-crys

The impact of charge differences on the carborane clusters as a result of different data sources and/or optimization methods/calculations was also explored with Glide (see Supporting Information, Table 2SI). As in the case of AutoDock docking, the binding energies were lower for the nido-carboranyl antifolates compared to the nido-carboranyl counterparts. Overall, differences in avg. binding energies, as a result of different geometries and partial atom charges, were not as pronounced as in the case of AutoDock while RMSDs and pose clustering were comparable.

Since FlexX and Surflex operate with formal charges,27,28,30,31 only docking of 1 C/2 N-crys and 1 C/2 N-con was evaluated with these programs. The latter carboranyl antifolates were chosen to address the influence of geometric differences on the docking performances. Automated docking of 1 C-crys and 2 N-crys by FlexX produced top ranked poses with RMSDs of 3.19 Å and 3.18 Å, respectively. The top 25 % of the docked poses for both 1 C-crys and 2 N-crys had avg. RMSD's of 5.97 Å and 5.89 Å respectively (Tables 3 and and4).4). The automated mode selected the PPs of 1 C-crys and 2 N-crys as the base fragments and PA3 before incremental construction of the ligands. User defined base fragment placement (BFP) of the PP of 1 C-crys followed by PA1 and PA2 improved the RMSDs of the top-ranking poses by factors of 1.6 and 2, respectively, and the avg. RMSDs of top 25% poses by factors of 3.2 and 1.4, respectively. The corresponding improvement factors in the case 2 N-crys were 2.5, 1.6, 1.7, and 1.0. As expected, user defined BFP of PP, followed by PA3, produced the same results as automated docking. 1 C-con and 2 N-con produced similar FlexX docking results indicating that in the absence of crystal structure information, binding interactions of virtual carboranyl compounds with proteins can be predicted accurately (see Supporting Information, Table 3SI).

Table 3
The effects of automated docking mode, selection of base fragment for placement (BFP), and placement algorithm (PA) on FlexX docking of 1C-crys and 2 N-crys
Table 4
Surflex docking of 1C-crys and 2 N-crys

For BFP, both in automated and user defined mode, FlexX favors PP. User defined BFP of CPs, in particular in the cases of 2 N-crys and 1 C/2 N-con, was inferior to BFP of PP indicating that hydrogen bonds and/or ionic interactions may be favored in this initial phase of the FlexX docking.27,28 This may also explain the fact that user defined BFP of CPs, produced poor docking solutions with high RMSDs. PA1 was found to be superior to PA2 and PA3 following initial BFP of PP. PA1 only utilizes interactions that are hydrophobic and unspecific in nature.27 This indicates that the second phase of incremental construction in FlexX docking, involving the PAs, is mainly guided by the hydrophobicity of the CPs.

The results of Surflex docking of 1 C-crys and 2 N-crys are summarized in Table 4. Using default docking parameters without fragment placement, the top 25% poses of 1 C-crys have an avg. RMSD of 0.78 Å and the top ranking pose had an RMSD of 0.89 Å. Docking of 1 C-crys using the “multistart5” option (5 different starting positions) without fragment placement decreased the docking quality of both top 25% poses and top ranked pose by a factor ~ 2. Docking with fragment placement was not explored because the default settings produced satisfactory docking results. In the case 1 C-con, there were no significant differences in docking quality between default settings and the “-multistart5” option (see Supporting Information, Table 4SI).

Docking of 2 N-crys using default parameters without fragment placement was unsatisfactory resulting in an avg. RMSD of the top 25 % poses of 3.25 Å. The “multistart5” option without fragment placement improved the docking by a factor of 1.65 (1.97 Å). The difference in RMSDs of the top ranked poses was insignificant (0.741 vs. 0.755 Å). These values improved even further for docking with fragment placement. When PP was selected as a placed fragment, the corresponding avg. RMSD and the top ranking pose RMSD decreased to 0.731 and 0.671 Å, respectively. The corresponding values for fragment placement of CPs were 1.26 Å and 0.463 Å. This is in contrast to FlexX, where BFP of the nido-CP led to a significant deterioration of docking quality. Overall, the Surflex docking results for 2 N-con showed similar patterns as for those of 2 N-crys (see Supporting Information, Table 4SI).

In preliminary studies, we also explored the Molecular Operating Environment (MOE) program for docking of carboranyl antifolates into hDHFR. Docking was carried out in standard mode. MOE is based on simulated annealing, which calculates the grid based interaction energy between the docked ligand and the receptor/enzyme.46,47 MOE has parameter sets for sp2 and sp3 boron atoms but not for hexavalent boron. If at all, the docking performance of MOE seemed to be comparable with that of FlexX in automated mode. Both replacement of boron with “C.3” and the use of the boron atom types provided by MOE produced similar docking results.

Enzymatic studies indicated that antifolate 1 is approximately ten times more potent than 2 as an inhibitor of hDHFR,10 which is consistent with the binding energies that were predicted by all four docking programs for 1 and 2 (Figure 4). However, these free energy values should be viewed with caution. Scoring functions are in general not very accurate with significant differences between various docking algorithms.48 Additionally, the hydrophobic interactions were calculated without applying appropriate force fields for hexavalent boron.

Figure 4
Average total binding energies of docked 1 C-crys (black bar) and 2 N-crys (grey bar) with hDHFR. For AutoDock and Glide, Mulliken charges at AM1 level were considered for 1C-crys and 2 N-crys. For Surflex, FlexX and Glide, the avg. binding energies for ...

Summary and Conclusions

This is the first report of docking studies with compounds containing the negatively charged nido-carboranyl cluster. Using 1 C-crys and 2 N-crys as examples, the docking performances of AutoDock, Glide, FlexX, and Surflex are summarized in Figure 5. Under optimized conditions (e.g. rigid docking and fragment placements), docking of 1 C-crys with AutoDock and Glide were comparably good followed by Surflex and then FlexX. In the case of 2 N-crys, AutoDock, Glide, and Surflex showed comparably good docking followed by FlexX. Under optimized conditions, docking of 1 C-crys and 2 N-crys is comparable when AutoDock, Glide, and FlexX were used. In the case of Surflex, however, 2 N-crys docked better than 1 C-crys. Overall, a comparative evaluation of 1 C-con and 2 N-con produced similar results (Supporting Information, Figure 2SI).

Figure 5
Comparative docking evaluation of 1 C-crys (A) and 2 N-crys (B). The black bars indicate the percentage of the total poses with RMSD below 2Å and the grey bars indicate the percentage of the total poses with RMSDs below 4 Å. Data used ...

Differences in geometries and partial atom charges resulting from different data sources and/or optimization methods/calculations did not impact the docking performances of AutoDock and Glide significantly. A careful selection of placement algorithm and/or fragment placement improved the docking quality of both FlexX and Surflex. The scoring values generated by all four programs were in accordance with experimental data. The observed docking patterns of closo- and nido-carboranyl antifolates with hDHFR validate previously developed methods for structure and ligand-based design of carborane-containing compounds.16, 49

The force field parameters developed for sp3 borons by Otkidach et al.20,45 were partially implemented in AutoDock for docking of the carboranyl antifolates. The described docking model appears to remain unaffected by this substitution, which further validates our strategy to enable docking of biomolecules containing boron clusters by replacing boron with carbon. It should be noted, however, that apart from AutoDock the implementation of these parameters in the other softwares packages discussed in the present study is not straightforward.

Carborane cages are rigid, hydrophobic entities and, intuitively, the docking of these scaffolds could be visualized as the docking of a single hydrophobic entity into the hydrophobic pocket of hDHFR. From this perspective, it is difficult to envision how the interactions of this single entity with its amino acid environment are changed during the docking process by replacing hexavalent boron either with carbon or sp3 boron. Although the reported docking protocols could be validated successfully, there is certainly further need for the development of accurate forcefield parameters, probably for the entire cage structures rather than for individual atom types.

A drawback of the present study was that the carboranyl antifolates used for docking only contained two rotatable bonds at the methylene bridge between the carboranyl- and the pyrimidyl portions. Consequently, the torsional degrees of freedom were limited in the present study and evaluating more complex systems may be necessary to confirm that the reported docking protocols are also applicable to other biomolecules that contain carborane clusters. Unfortunately, such systems are currently not available.

Supplementary Material

1_si_001

2_si_002

3_si_003

Acknowledgments

The presented work was financially supported by funds from The Ohio State University College of Pharmacy and, in part, by NIH grant 1 R01 CA127935-01A2. The authors would like to thank Drs. Steven Ealick and for generously providing the small molecule crystal structure parameters of 1 and 2. Dr. David Borhani and Dr. Richard Dixon (Abbott Bioresearch Centre, Worcester, MA), Dr. Holger Claussen (BioSolveIT), Dr. Christopher Williams (Chemical Computing Group) and Dr. Jeffrey Saunders (Schrodinger, LLC) are acknowledged for helpful discussions. Authors of this manuscript are grateful to In Hee Park for their help with implementation of boron parameters in AutoDock Tools.

List of abbreviation

1 C-prot
closo-o-carboranyl antifolate (1) co-crystallized with hDHFR (PDB # 2C2S)
2 N-prot
nido-o-carboranyl antifolate (2) co-crystallized with the hDHFR (PDB # 2C2T)
1 C-crys
small molecule crystal structure of the closo-o-carboranyl antifolate (1)
2 N-crys
small molecule crystal structure of the nido-o-carboranyl antifolate (2)
1 C-con
de novo constructed form of closo-o-carboranyl antifolate (1)
2 N-con
de novo constructed form of the nido-o-carboranyl antifolate (2)
ADT
AutoDock Tools
APT
Atomic Polar Tensors
Avg
Average
BFP
Base Fragment Placement
ESPfit
Electrostatic Potential Fit Method
BNCT
Boron Neutron Capture Therapy
Du
Dummy atom type
hDHFR
human Dihydrofolate Reductase
LGA
Lamarckian Genetic Algorithm
MOE
Molecular Operating Environment
NADPH
Dihydronicotinamide Adenine Dinucleotide Phosphate
PA1
Placement Algorithm 1
PA2
Placement Algorithm 2
PA3
Placement Algorithm 3
PP
pyrimidyl portion used as a base fragment
CP
closo-o-carboranyl- or nido-o-carboranyl portion used as a base fragment
rdf
receptor description file
RMSD
Root Mean Square Deviation
SP
Standard Precision

Footnotes

#Presented in part at the European Science Foundation (ESF) Exploratory Workshop “BioBor-Exploring New Opportunities of Boron Chemistry Towards Medicine”, in Lodz, Poland, May 9-12, 2008.

Supporting Information Available. Additional data, as indicated above, is available in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org.

References

1. Armstrong AF, Valliant JF. The bioinorganic and medicinal chemistry of carboranes: from new drug discovery to molecular imaging and therapy. Dalton Trans. 2007;38:4240–4251. [PubMed]
2. Lesnikowski ZJ. Boron units as pharmacophores - new applications and opportunities of boron cluster chemistry. Collect Czech Chem Commun. 2007;72:1646–1658.
3. Endo Y, Iijima T, Yaguchi K, Kawachi E, Inoue N, Kagechika H, Kubo A, Itai A. Structure-Activity study of retinoid agonists bearing substituted dicarba-closo-dodecaborane. Relation between retinoidal activity and conformation of two aromatic nuclei. Bioorg Med Chem Lett. 2001;11:1307–1311. [PubMed]
4. Ogawa T, Ohta K, Yoshimi T, Yamazaki H, Suzuki T, Ohta S, Endo Y. m-Carborane bisphenol structure as a pharmacophore for selective estrogen receptor modulators. Bioorg Med Chem Lett. 2006;16:3943–3946. [PubMed]
5. Goto T, Ohta K, Suzuki T, Ohta S, Endo Y. Design and synthesis of novel androgen receptor antagonists with sterically bulky icosahedral carboranes. Bioorg Med Chem. 2005;13:6414–6424. [PubMed]
6. Endo Y, Yoshimi T, Kimura K, Itai A. Protein kinase C modulators bearing dicarba-closo-dodecaborane as a hydrophobic pharmacophore. Bioorg Med Chem Lett. 1999;9:2561–2564. [PubMed]
7. Tsuji M, Koiso Y, Takahashi H, Hashimoto Y, Endo Y. Modulators of tumor necrosis factor alpha production bearing dicarba-closo-dodecaborane as a hydrophobic pharmacophore. Biol Pharm Bull. 2000;23:513–516. [PubMed]
8. Julius RL, Farha OK, Chiang J, Perry LJ, Hawthorne MF. Synthesis and evaluation of transthyretin amyloidosis inhibitors containing carborane pharmacophores. Proc Natl Acad Sci U S A. 2007;104:4808–4813. [PubMed]
9. Page MFZ, Jalisatgi SS, Maderna A, Hawthorne MF. Design and synthesis of a candidate alpha -human thrombin irreversible inhibitor containing a hydrophobic carborane pharmacophore. Synthesis. 2008;4:555–563.
10. Reynolds RC, Campbell SR, Fairchild RG, Kisliuk RL, Micca PL, Queener SF, Riordan JM, Sedwick WD, Waud WR, Leung AKW, Dixon RW, Suling WJ, Borhani DW. Novel boron-containing, nonclassical antifolates: Synthesis and preliminary biological and structural evaluation. J Med Chem. 2007;50:3283–3289. [PubMed]
11. Tolmachev V, Koziorowski J, Sivaev I, Lundqvist H, Carlsson J, Orlova A, Gedda L, Olsson P, Sjoeberg S, Sundin A. Closo-Dodecaborate(2-) as a Linker for Iodination of Macromolecules. Aspects on Conjugation Chemistry and Biodistribution. Bioconjugate Chem. 1999;10:338–345. [PubMed]
12. Wilbur DS, Chyan MK, Hamlin DK, Vessella RL, Wedge TJ, Hawthorne MF. Reagents for astatination of biomolecules. 2. Conjugation of anionic boron cage pendant groups to a protein provides a method for direct labeling that is stable to in vivo deastatination. Bioconjugate Chem. 2007;18:1226–1240. [PubMed]
13. Kozisek M, Cigler P, Lepsik M, Fanfrlik J, Rezacova P, Brynda J, Pokorna J, Plesek J, Gruner B, Grantz Saskova K, Vaclavikova J, Kral V, Konvalinka J. Inorganic polyhedral metallacarborane inhibitors of HIV protease: A new approach to overcoming antiviral resistance. J Med Chem. 2008;51:4839–4843. [PubMed]
14. Moore BS, Eustaquio AS, McGlinchey RP. Advances in and applications of proteasome inhibitors. Curr Opin Chem Biol. 2008;12:434–440. [PMC free article] [PubMed]
15. Hawthorne MF, Lee Mark W. A critical assessment of boron target compounds for boron neutron capture therapy. J Neurooncol. 2003;62:33–45. [PubMed]
16. Bandyopadhyaya AK, Tiwari R, Tjarks W. Comparative molecular field analysis and comparative molecular similarity indices analysis of boron-containing human thymidine kinase 1 substrates. Bioorg Med Chem. 2006;14:6924–6932. [PubMed]
17. Johnsamuel J, Byun Y, Jones TP, Endo Y, Tjarks W. A new strategy for molecular modeling and receptor-based design of carborane containing compounds. J Organomet Chem. 2003;680:223–231.
18. Martichonok V, Jones JB. Cysteine proteases such as papain are not inhibited by substrate analog peptidyl boronic acids. Bioorg Med Chem. 1997;5:679–684. [PubMed]
19. Minkkila A, Saario SM, Kasnanen H, Leppanen J, Poso A, Nevalainen T. Discovery of boronic acids as novel and potent inhibitors of fatty acid amide hydrolase. J Med Chem. 2008;51:7057–7060. [PubMed]
20. Tafi A, Agamennone M, Tortorella P, Alcaro S, Gallina C, Botta M. AMBER force field implementation of the boronate function to simulate the inhibition of beta -lactamases by alkyl and aryl boronic acids. Eur J Med Chem. 2005;40:1134–1142. [PubMed]
21. Byun Y, Thirumamagal BTS, Yang W, Eriksson S, Barth RF, Tjarks W. Preparation and biological evaluation of 10B-Enriched 3-[5-{2-(2,3-Dihydroxyprop-1-yl)-o-carboran-1-yl}pentan-1-yl]thymidine (N5-2OH), a new boron delivery agent for boron neutron capture therapy of brain tumors. J Med Chem. 2006;49:5513–5523. [PubMed]
22. Narayanasamy S, Thirumamagal BT, Johnsamuel J, Byun Y, Al-Madhoun AS, Usova E, Cosquer GY, Yan J, Bandyopadhyaya AK, Tiwari R, Eriksson S, Tjarks W. Hydrophilically enhanced 3-carboranyl thymidine analogues (3CTAs) for boron neutron capture therapy (BNCT) of cancer. Bioorg Med Chem. 2006;14:6886–6899. [PubMed]
23. Chazalette C, Riviere-Baudet M, Scozzafava A, Abbate F, Maarouf ZB, Supuran CT. Carbonic anhydrase inhibitors, interaction of boron derivatives with isozymes I and II: a new binding site for hydrophobic inhibitors at the entrance of the active site as shown by docking studies. J Enzyme Inhib. 2001;16:125–133. [PubMed]
24. Endo Y, Iijima T, Yamakoshi Y, Fukasawa H, Miyaura C, Inada M, Kubo A, Itai A. Potent estrogen agonists based on carborane as a hydrophobic skeletal structure: a new medicinal application of boron clusters. Chem Biol. 2001;8:341–355. [PubMed]
25. Huey R, Morris GM, Olson AJ, Goodsell DS. A semiempirical free energy force field with charge-based desolvation. J Comput Chem. 2007;28:1145–1152. [PubMed]
26. Morris GM, Goodsell DS, Halliday RS, Huey R, Hart WE, Belew RK, Olson AJ. Automated docking using a lamarckian genetic algorithm and an empirical binding free energy function. J Comput Chem. 1998;19:1639–1662.
27. Rarey M, Kramer B, Lengauer T. Docking of hydrophobic ligands with interaction-based matching algorithms. Bioinformatics. 1999;15:243–250. [PubMed]
28. Rarey M, Kramer B, Lengauer T, Klebe G. A fast flexible docking method using an incremental construction algorithm. J Mol Biol. 1996;261:470–489. [PubMed]
29. Rarey M, Wefing S, Lengauer T. Placement of medium-sized molecular fragments into active sites of proteins. J Comput-Aided Mol Des. 1996;10:41–54. [PubMed]
30. Jain AN. Surflex: fully automatic flexible molecular docking using a molecular similarity-based search engine. J Med Chem. 2003;46:499–511. [PubMed]
31. Jain AN. Surflex-Dock 2.1: Robust performance from ligand energetic modeling, ring flexibility, and knowledge-based search. J Comput-Aided Mol Des. 2007;21:281–306. [PubMed]
32. 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. J Med Chem. 2004;47:1739–1749. [PubMed]
33. Friesner RA, Murphy RB, Repasky MP, Frye LL, Greenwood JR, Halgren TA, Sanschagrin PC, Mainz DT. Extra precision glide: Docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J Med Chem. 2006;49:6177–6196. [PubMed]
34. Mulliken RS. Electronic population analysis on LCAO-MO [linear combination of atomic orbital-molecular orbital] molecular wave functions. I. J Chem Phys. 1955;23:1833–1840.
35. Cioslowski J. A new population analysis based on atomic polar tensors. J Am Chem Soc. 1989;111:8333–8336.
36. Dewar MJS, Jie C, Zoebisch EG. AM1 calculations for compounds containing boron. Organometallics. 1988;7:513–521.
37. Dewar MJS, Zoebisch EG, Healy EF, Stewart JJP. Development and use of quantum mechanical molecular models. 76. AM1: a new general purpose quantum mechanical molecular model. J Am Chem Soc. 1985;107:3902–3909.
38. Chandrasekhar J, Andrade JG, Schleyer PVR. Efficient and accurate calculation of anion proton affinities. J Am Chem Soc. 1981;103:5609–5612.
39. Clark T, Chandrasekhar J, Spitznagel GW, Schleyer PVR. Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21 + G basis set for first-row elements, lithium to fluorine. J Comput Chem. 1983;4:294–301.
40. Fox MA, Goeta AE, Howard JA, Hughes AK, Johnson AL, Keen DA, Wade K, Wilson CC. The molecular structure of (PSH+)(nido-7,8-C2B9H12-) determined by neutron diffraction (PS = proton sponge, 1,8-bis(dimethylamino)naphthalene) Inorg Chem. 2001;40:173–175. [PubMed]
41. Cole JC, Murray CW, Nissink JW, Taylor RD, Taylor R. Comparing protein-ligand docking programs is difficult. Proteins. 2005;60:325–332. [PubMed]
42. Holt PA, Chaires JB, Trent JO. Molecular docking of intercalators and groove-binders to nucleic acids using AutoDock and Surflex. J Chem Info Model. 2008;48:1602–1615. [PMC free article] [PubMed]
43. Fanfrlik J, Hnyk D, Lepsik M, Hobza P. Interaction of heteroboranes with biomolecules. Part 2. The effect of various metal vertices and exo-substitutions. Phys Chem Chem Phys. 2007;9:2085–2093. [PubMed]
44. Fanfrlik J, Lepsik M, Horinek D, Havlas Z, Hobza P. Interaction of carboranes with biomolecules: formation of dihydrogen bonds. Chemphyschem. 2006;7:1100–1105. [PubMed]
45. Otkidach DS, Pletnev IV. Conformational analysis of boron-containing compounds using Gillespie-Kepert version of molecular mechanics. Theochem. 2001;536:65–72.
46. Esposito EX, Baran K, Kelly K, Madura JD. Docking of sulfonamides to carbonic anhydrase II and IV. J Mol Graphics Modell. 2000;18:283–289. 307–288. [PubMed]
47. MOE (The Molecular Operating Environment) version 2007.09. Chemical Computing Group Inc.; Montreal, Canada H3A 2R7: 2007.
48. Warren GL, Andrews CW, Capelli AM, Clarke B, LaLonde J, Lambert MH, Lindvall M, Nevins N, Semus SF, Senger S, Tedesco G, Wall ID, Woolven JM, Peishoff CE, Head MS. A critical assessment of docking programs and scoring functions. J Med Chem. 2006;49:5912–5931. [PubMed]
49. Johnsamuel J, Byun Y, Jones TP, Endo Y, Tjarks W. A convenient method for the computer-aided molecular design of carborane containing compounds. Bioorg Med Chem Lett. 2003;13:3213–3216. [PubMed]