AutoDock4(Zn): an improved AutoDock force field for small-molecule docking to zinc metalloproteins.

Zinc is present in a wide variety of proteins and is important in the metabolism of most organisms. Zinc metalloenzymes are therapeutically relevant targets in diseases such as cancer, heart disease, bacterial infection, and Alzheimer's disease. In most cases a drug molecule targeting such enzymes establishes an interaction that coordinates with the zinc ion. Thus, accurate prediction of the interaction of ligands with zinc is an important aspect of computational docking and virtual screening against zinc containing proteins. We have extended the AutoDock force field to include a specialized potential describing the interactions of zinc-coordinating ligands. This potential describes both the energetic and geometric components of the interaction. The new force field, named AutoDock4Zn, was calibrated on a data set of 292 crystal complexes containing zinc. Redocking experiments show that the force field provides significant improvement in performance in both free energy of binding estimation as well as in root-mean-square deviation from the crystal structure pose. The new force field has been implemented in AutoDock without modification to the source code.

Introduction

Zinc is present in numerous biological structures and is found in virtually all aspects of metabolism across multiple species.[1] It can play a structural role as in zinc finger proteins, the most prevalent proteins in eukaryotic genomes,[2] and is present in all enzyme classes,[3] usually in the form of coordinated Zinc(II) or Zn2+ ion. Zinc metalloenzymes are therapeutically relevant targets in many diseases, like heart disease,[4] cancer,[5,6] bacterial infections,[7] and Alzheimer’s.[8,9] In most cases, a drug molecule establishes coordination bonds with the zinc ion[10] present in the protein; thus, an accurate description of this interaction is crucial for drug design.

To properly model the zinc coordination interactions, two issues should be addressed: the coordination geometry and the interaction strength.

Most force fields describe metal coordination using descriptions derived from the original Stote and Karplus nonbonded model,[11] where the interaction is described using Lennard-Jones and Coulomb potentials. This description relies on assignment of partial charges,[11−13] and thus accuracy becomes strongly dependent on the choice of charge model. Also, an electrostatic model based on the filled valence orbital of the Zn2+ ion fails to explain the prevalence of histidine and cysteine over the more electronegative carboxylate groups of glutamate and aspartate as the most frequent zinc coordinating residues.[14−17]

Moreover, some high potency inhibitors coordinate Zn2+ via uncharged nitrogens with electron lone pairs, such as those found in sulfonamides[18] and imidazoles[6] (see Figure 13), that seem to interact more strongly than negatively charged nitro groups.[7] Recently, DFT calculations were used to calibrate a nondirectional zinc coordination force field independent of atomic partial charges.[19]

Figure 13

Comparison of redocking accuracy of 1s63 using (a) standard AutoDock4, (b) AutoDock Vina, and (c) AutoDock4Zn force fields (RMSD are shown in parentheses). Zinc-coordinating residues and experimental ligand pose are shown as thin gray sticks; zinc is cyan; docked poses are shown as green thick sticks; the location and the optimal radius of the TZ pseudoatom potential is shown as semitransparent sphere (red). Hydrogens are not shown for sake of clarity.

Polarization and charge transfer models[20−22] could provide a more accurate description, although their computational complexity makes them unsuitable for dockings, which typically involve a large number of energy estimations over the course of the calculation.

The coordination geometry issue is addressed differently by bonded and nonbonded models. It has been recently demonstrated that zinc exhibits a strong preference for the tetrahedral geometry, with some of the previously observed variability in coordination spheres being artifactual.[17] Bonded models, such as the Zinc AMBER Force Field,[23] describe the tetrahedral coordination with harmonic potentials and angle terms for an explicit bond that provides directionality. Due to the requirement of the explicit bonds where ligands coordinate with zinc, bonded models are not suitable for docking calculations.

In nonbonded models, few force fields provide directional potentials. Two examples that do are the cationic dummy atom model[5] and the scoring function implemented in FlexX.[24] However, in this latter case, while improving coordination geometry accuracy, no improvements in binding energy prediction were reported.

To be suitable for docking, and virtual screens in particular, modeling the interaction with zinc must provide a description of the geometry that is computationally efficient, and good accuracy in the estimation of the interaction strength.

In this paper, we report the development of a directional, charge-independent model for zinc-coordination force field for AutoDock4 which provides higher accuracy than the standard force field. Interactions are modeled independently for different atom types, providing specific potentials for each one.

Over the years, AutoDock and its force field were modified by us and others to improve scoring of low affinity ligands[25] or to obtain a scoring function tailored to specific targets, like kinases.[26] Other methods added new features, like receptor flexibility models,[27,28] flexible macrocycle docking,[29] and docking with waters.[30] The highly customizable architecture of the program allows implementation of substantial modifications relatively easily and without requiring source code changes.

Methods

To implement the new zinc-coordination model, first we identified a data set of high-quality complexes for which experimental affinity values had been determined. The data set analysis enabled determination of the parameters for geometrical terms that were then calibrated to fit within the AutoDock force field. The new force field, named AutoDock4Zn, was then cross-validated on the data set.

Data Set Creation

To design the new force field, a suitable set of zinc metalloprotein–ligand complexes was defined. The ligands cover a wide range of structure diversity and binding affinity, thus providing an optimal calibration set for generic applicability of the force field for drug design.

In order to build the data set, the Binding MOAD[31] was filtered using the following criteria: (a) presence of at least one zinc ion; (b) experimentally determined inhibition (Ki) or dissociation (Kd) constants; (c) no alternate conformation or missing atoms for the ligand, and (d) no alternate side chain conformations in receptor residues within 5 Å from any ligand atom. This filtering led to a set of 510 complexes, which were downloaded from the Protein Data Bank.[32]

These complexes were then analyzed to isolate and characterize the zinc coordination geometry within the receptor and its interaction with ligands. Each complex was classified accordingly to the number of receptor (r) and ligand (l) atoms within coordination distance (≤2.8 Å for sulfur atoms, ≤2.5 Å for all others) from the zinc ion,[24] and denoted as Zn.

A specific treatment was used to analyze the coordination geometry of carboxylic acids. Carboxylic acids from aspartate or glutamate side chains have been described to coordinate zinc mainly with bidentate, monodentate, syn or anti modes. However, it has been demonstrated that carboxylate groups can adopt any coordination geometry ranging between mono- and bidentate,[33,34] which is poorly described by a discrete classification scheme. To address this issue, carboxylic acid groups on receptors were always considered as monovalent and represented by a weighted average of the position of the two coordinating oxygens (see Figure 6); the method used to calculate the weighted average for carboxylic acids is described in the Supporting Information. The distribution of different coordination geometries is summarized in Table 1. The most represented coordination geometry in our data set is the tetrahedral one (Zn3,1, Zn4,0), that was indeed found to be the most common in biological systems.[17] Other geometries, like five- (Zn4,1, Zn3,2) and six-coordinated (Zn3,3, Zn4,2), were also found, but were much rarer.

Figure 6

Definition of carboxyl group average atom. Details about the methods are reported in the Supporting Information.

Table 1

Number of Zn Classes for Each Zinc Ion Found in the Initial Data Seta

	l
r	0	1	2	3
0	14	3	1	0
1	56	1	0	0
2	72	2	5	0
3	57	244	43	3
4	214	12	15	0
5	8	0	0	0

Complexes with at least one the classes in bold were selected for the final data set.

Unoccupied coordination geometry locations are usually engaged by a of the resides coordinating the metal modulates association energies of the ligands.[35] To our knowledge, no experimental structures have been reported where ligands interacted with zinc through the mediation of a water molecule. In fact, dockings performed with the hydrated ligand protocol[30] consistently predicted the displacement of the waters. Therefore, zinc was considered always desolvated, while desolvation energy was implicitly accounted for during the calibration of the force field.

Complexes where ligands were not directly involved in zinc coordination (i.e., l = 0) were discarded. This included also Zn4,0 cases, where zinc plays a structural role helping protein folding,[14,15] coordinating four cysteine side chains. Some of the Zn4,0 cases were misclassified as Zn2,0 because the zinc ion bridges two monomeric units that were split during the analysis process.

Geometries where receptor atoms were not involved, or only partially involved in zinc interactions (0 ≤ r ≤ 2) were also discarded upon visual inspection. In particular, the Zn0,0 class contains complexes where Zn is used as an aid in crystallization and has no biological significance, surrounding the protein structures often in large number and at toxic concentrations.[17] Finally, 5 complexes involving serine protease inhibitors from Zn2,2 class, were discarded because zinc is known to be recruited transiently as coinhibitor only, and it is not consistently present in the binding site.[36]

This left four coordination classes, Zn3,1, Zn3,2, Zn3,3, Zn4,1, and Zn4,2, resulting in a calibration set of 292 unique complexes. A summary of ligand properties in the set is shown in Figure 1.

Figure 1

Summary of the distributions of ligand properties in the final data set: molecular weight (a), LogP (b), number of heavy atoms (c), torsional degrees of freedom (d), experimental free energy of binding (e).

For complexes where the tetrahedral coordination geometry is possible (Zn3,), we analyzed the distribution of the ligand atoms coordinating zinc, using the AutoDock atom types: NA (nitrogen HB acceptor), N (nitrogen non-HB acceptor), OA (oxygen HB acceptor), and SA (sulfur HB acceptor). The ideal zinc tetrahedral geometry was calculated with respect to the averaged position of receptor atoms. The tetrahedral plane was defined as the plane calculated between average coordinating receptor atoms and the zinc atom (Figure 7). QM optimizatons performed on few representative cases confirmed high quality of experimental tetrahedral geometries (see the Supporting Information).

Figure 7

Tetrahedral zinc geometry. (a) Ligand and receptor atoms are shown as sticks colored by atom type. The tetrahedral plane defined by three receptor atoms (black spheres) is determined. The TZ pseudoatom is located at unoccupied corner of the ideal tetrahedral geometry. Coordination geometry is calculated on weighted average oxygen positions from carboxylic side chains. (b) The potential for atom type NA (nitrogen acceptor) is shown as iso-contour surfaces (cyan).

Then, we measured the deviation of ligand atoms from the ideal position in the tetrahedral geometry, defined as the angle between the vector Zn-ligand atom and the tetrahedral plane.

In Figures 2–5 are shown the tridimensional scattering coordinating ligand atoms with respect to the tetrahedral zinc (a and b) and their angle deviations (c). More details about the alignment method and analysis are reported in the Supporting Information.

Figure 2

Distribution of 137 NA atom types coordinating zinc: (a) perspective projection; (b) top view; (c) angle histogram. Atoms are shown as spheres: receptor atoms (black), zinc (green), NA atoms (blue). Tetrahedral geometries are colored in gray; tetrahedral plane is shown as semitransparent polygon; pseudoatom location is shown as wireframe sphere.

Figure 5

Distribution of 27 SA atom types coordinating zinc: (a) perspective projection; (b) top view; (c) angle histogram. Atoms are shown as spheres: receptor atoms (black), zinc (green), SA atoms (yellow). Tetrahedral geometries are colored in gray; tetrahedral plane is shown as semitransparent polygon; pseudoatom location is shown as wireframe sphere.

Distribution of 151 OA atom types coordinating zinc: (a) perspective projection; (b) top view; (c) angle histogram. Atoms are shown as spheres: receptor atoms (black), zinc (green), OA atoms (red). Tetrahedral geometries are colored in gray; tetrahedral plane is shown as semitransparent polygon; pseudoatom location is shown as wireframe sphere.

The deviation analysis showed that nitrogen HB acceptor (NA) is consistently found very close to the ideal position (>80% within ≤10°, Figure 2). On the other hand, the placement of nitrogen non-HB acceptor (N) and oxygen (OA) and sulfur (SA) is less well-defined, appearing to be dependent solely on the accessibility of the zinc atom in the receptor (Figures 3–5).

Figure 3

Distribution of 15 N atom types coordinating zinc: (a) perspective projection; (b) top view; (c) angle histogram. Atoms are shown as spheres: receptor atoms (black), zinc (green), N atoms (blue). Tetrahedral geometries are colored in gray; tetrahedral plane is shown as semitransparent polygon; pseudoatom location is shown as wireframe sphere.

New Force Field

The standard AutoDock force field supports several ligand–metal interactions.[37] Similar to all other pairwise interactions in the force field, the interaction between ligand atoms and metals contained in the receptor is described mainly by van der Waals (ΔHvdW) and Coulomb electrostatic (ΔHelec) terms and, to a smaller degree, by the desolvation term (ΔGdesolv). This approach has several limitations. First, the van der Waals equilibrium distances for the atoms involved in zinc coordination are significantly larger than the coordination distances[13,19] (i.e., for nitrogen, the vdW equilibrium distance is 2.49 Å, compared to coordination distance of 2.0 Å). Second, due to the lack of a specialized terms for the metal coordination, directionality is not accounted for. Finally, while the electrostatic term is very effective in describing interactions involving partial charges, it makes the energy function highly sensitive to strongly charged groups, such as metals with formal charges. Also, in the Gasteiger[38] charge model used in AutoDock,[37] oxygen atoms are systematically assigned a more negative charge than nitrogen and sulfur, thus resulting in the preferred candidates for chelating positively charged metal. While this approach is accurate enough for magnesium ion interactions, it is not sufficient to properly describe zinc coordination preferences.

From our data set analysis, we found that the coordination of zinc requires a specialized treatment, so we modified the standard AutoDock force field. The standard force field includes the following terms (eq 1):[37]where the free energy of binding (FEB) is calculated as a sum of van der Waals (vdW), hydrogen bond (Hbond), Coulomb electrostatic (Elec), desolvation (desolv) and ligand torsional entropy (torsDoF); each term is weighted by a specific value (Wterm) estimated using a linear regression model.[37] To extend the force field, we first disabled the electrostatic potential for zinc by setting its partial charge to zero. Then, the pairwise interactions of each atom types involved in zinc coordination was defined as a new potential energy term. For N, OA, and SA atom types, spherical potentials VZn,N, VZn,OA, and VZn,SA were defined to reflect the known coordination distances, by adapting the van der Waals potential in the AutoDock force field (eq 2):The pairwise equilibrium distance r between zinc and N, OA, and SA atom types was set to 2.0, 2.1, and 2.25 Å, respectively, and independent ε well-depth values were estimated. Spherical potentials are particularly suitable for accurately reproducing hydroxamate coordination geometries.[39,40]

For the NA type a new directional tetrahedral potential VTZ,NA was defined and the interaction with zinc was split in two separate components. The repulsive component is mediated by the zinc atom, while the attractive component is mediated by a new pseudoatom TZ that has been added to the standard force field table.[41] The pseudoatom interacts only with NA, therefore no interaction is defined with any other atom type. The pseudoatom is added in the receptor structure for all complexes where the tetrahedral coordination geometry is present, i.e., all Zn3, classes, where only three receptor atoms are coordinating zinc. The pseudoatom is placed at the unoccupied vertex of the tetrahedral geometry, located at the optimal coordination distance for nitrogen (r = 2.0 Å) (Figure 7a), and an attractive 12–6 potential with a corresponding ε is defined (Figure 7b). Finally, the zinc–hydrogen pairwise interaction was eliminated to prevent clashes that would interfere with the proper interaction between groups like sulfonamide −NH2, or hydroxyl, with zinc. This allows ligands to establish the proper coordination interaction independent of the orientation of the hydrogen with respect to the heavy atom.

Therefore, the following potential was added to eq 1:and the FEB becomes the linear combination of the five standard AutoDock terms plus the new zinc coordination pairwise potential.

All modifications to the AutoDock force field were made by adapting the force field table and parameter files, without source code modifications. The details of the implementation are described in the Supporting Information.

The ε values for eq 3 were then calibrated independently from each other and from the other terms in eq 1.

Calibration Protocol

The new force field was calibrated with an iterative least-squares scheme. Initial attempts to calibrate combined terms from eqs 1 and 3 led to performance degradation in nonzinc complexes. Optimization of different term combination were tried, and best results were obtained by optimizing only terms in eq 3, while keeping the standard terms (eq 1) unmodified.

The calibration protocol consisted of the following steps: (a) crystallographic structures of the ligands were minimized with the current version of the force field, using Solis-Wet local search implemented in AutoDock;[37] (b) unweighted terms were calculated from minimized structures; (c) a regression model was built; (d) weights from the new regression model were used in the next minimization step. The protocol iterated through steps c and d five times to achieve convergence; stable weight values were achieved after the first two iterations.

Initial calibration results and cross-validation tests showed that no statistical significance could be achieved for the VZn,N term. This is likely due to insufficient experimental data for the N atom type. Therefore, standard force field term for this interaction (i.e., van der Waals) was restored, while keeping the correct equilibrium radius (2.0 Å) identified in the analysis. Then the calibration was repeated omitting the VZn,N term from eq 3. Final force field weights were selected from the last iteration, with a residual standard error was 2.804 kcal/mol. Final coefficients and extended analysis of the iterative calibration are described in the Supporting Information.

Results and Discussion

Predictive capabilities of the regression model were assessed with 5-fold cross-validation. The data set was divided in five bins containing an approximately uniform distribution of ligand atom types coordinating zinc, then redocking calculations were performed. Cross-validation docking results are summarized in Table 2. Details on docking preparation and RMSD calculations are available in the Supporting Information. Reproducing proper metal-coordination geometries and accurate energy estimations are notoriously difficult, especially for zinc.[42] The performance of AutoDock4Zn was evaluated accordingly to three different criteria: FEB estimation error, ligand pose RMSD calculated on all heavy atoms, and deviation from ideal zinc coordination geometry. Overall, the new AutoDock4Zn force field performed significantly better than both standard AutoDock and Vina force fields. These force fields provide roughly the same prediction errors in FEB estimations, while AutoDock4Zn consistently improved success rate (+50% with <1.0 kcal/mol) (Figure 8). The new force field also improves RMSD accuracy over the standard AutoDock force field in pose prediction accuracy (RMSD), producing results comparable with Vina (Figure 9). Not surprisingly, a remarkable improvement was achieved in reproducing the proper zinc-coordination geometry (RMSDZn), where AutoDock4Zn outperforms the two other force fields by a large amount (+127% success rate <1 Å, Figure 10).

Table 2

Cross-validation of Docking Performances and FEB Estimation Accuracy

	FEB error (kcal/mol)			RMSD (Å)		RMSDZn (Å)
	<1.0	<2.0	<3.0	<2.0	<2.5	<1.0	<1.5
AutoDock4Zn	32%	64%	81%	45%	51%	75%	80%
AutoDock4	18%	34%	53%	36%	42%	33%	46%
Vina	20%	38%	64%	45%	52%	37%	52%

Figure 8

Comparison of FEB prediction errors of the new force field with (a) standard AutoDock4 force field and (b) AutoDock Vina.

Figure 9

Comparison of RMSD error of the new force field with (a) standard AutoDock4 force field and (b) AutoDock Vina.

Figure 10

Comparison of RMSD error on zinc coordination geometry of the new force field with (a) standard AutoDock4 force field and (b) AutoDock Vina.

The use of specific potentials for describing the interaction is the main factor responsible for such an improvement, as shown by redocking experiments.

Examples

In some cases, where sulfur is directly involved in coordinating zinc, neither AutoDock4 nor Vina force fields were able to establish the proper interactions between the ligand and the zinc ion. A key example of this is provided by reredocking results of a potent inhibitor of neutral endopeptidase (NEP)[43] in the crystallographic complex with the PDB ID 1r1j (Figure 11). Results are summarized in Table 3. Both AutoDock4 and Vina predicted zinc to be coordinated by the carboxylate group, resulting in a misalignment of the ligand with respect to the receptor. The AutoDock4Zn, on the other hand, predicted the proper coordination by sulfanyl and carbonyl groups and provided more accurate FEB estimation. Similar results were found when redocking a potent aryl sulfonamide TACE inhibitor in the crystallographic complex with PDB ID 1oi0 (Figure 12). The increase in docking pose prediction and accuracy in the coordination geometry indentification resulted in a more precise FEB estimation (Table 4). It must also be noted also that no performance degredation was found when docking compounds without ligand-zinc interactions using the new force field. Improvements provided by AutoDock4Zn makes it suitable for virtual screening campaigns involving zinc.

Figure 11

Comparinson of redocking accuracy with 1r1j using (a) standard AutoDock4, (b) AutoDock Vina, and (c) AutoDock4Zn force fields (RMSD are shown in parentheses). Zinc-coordinating residues and experimental ligand pose are shown as thin gray sticks; docked poses are shown as green thick sticks. Hydrogens are not shown for sake of clarity.

Table 3

Docking Performances and FEB Estimation Accuracy on NEP (1r1j)

	FEB error (kcal/mol)	RMSD (Å)	RMSDZn (Å)
AutoDock4Zn	+0.74	1.21	1.04
AutoDock4	+1.76	4.85	6.38
Vina	+4.68	9.19	4.78

Figure 12

Comparison of redocking accuracy of 2oi0 using (a) standard AutoDock4, (b) AutoDock Vina, and (c) AutoDock4Zn force fields (RMSD are shown in parentheses). Zinc-coordinating residues and experimental ligand pose are shown as thin gray sticks; zinc is cyan; docked poses are shown as green thick sticks. Hydrogens are not shown for sake of clarity.

Table 4

Docking Performances and FEB Estimation Accuracy on TACE (2oi0)

	FEB error (kcal/mol)	RMSD (Å)	RMSDZn (Å)
AutoDock4Zn	+0.33	2.10	0.65
AutoDock4	+1.56	7.98	11.43
Vina	+2.50	2.34	5.03

Conclusions

We extended the standard AutoDock force field to include a specialized potential describing interactions of zinc-coordinating ligands. The potential has a description for both energetic and geometric components of the interaction

The new force field, named AutoDock4Zn, was calibrated on 292 complexes containing zinc using an iterative linear regression model. Redocking experiments showed that the force field provides a considerable improvement in performance when compared to both standard AutoDock4 and Vina force fields. Improvements are particularly relevant in the accuracy of FEB estimations, as well as in reproducing the proper coordination geometry. In fact, AutoDock4Zn provides a significant advantage when docking zinc-coordinating ligands, as shown in the examples described.

Due to the fully configurable nature of AutoDock4, the new potential was implemented by modifying the standard force field tables and few Python helper scripts available at http://autodock.scripps.edu/.

Moreover, the potential itself does not add any overhead to the docking calculation. There is no increase in the search complexity nor in the computational power requirement, therefore docking speed is completely unaffected. Also, performance on dockings not involving zinc-coordination are unchanged. Therefore, accuracy increase and lack of computational overhead make the AutoDock4Zn suitable for virtual screening campaigns, particularly when coordination of zinc is important.

Table 5

Docking Performances and FEB Estimation Accuracy on Farnesyltransferase (1s63)

	FEB error (kcal/mol)	RMSD (Å)	RMSDZn (Å)
AutoDock4Zn	+0.07	0.68	0.34
AutoDock4	+3.52	7.62	7.78
Vina	+2.93	0.50	0.54