AutoDock2
AutoDock21 from 1996 is designed to predict the bound conformations of small, flexible ligands to macromolecular targets of known structure. This version combines simulated annealing for conformational searching with a rapid, grid-based method of energy evaluation using the AMBER force field.
Components
The scoring function in AutoDock2 estimates the binding free energy change (\(\Delta G_{binding}\)) using several energy terms that account for different types of interactions between the ligand and the protein. The total binding free energy change is the sum of all individual energy terms:
This empirical scoring function enables AutoDock2 to predict the binding affinity of a ligand to a protein by evaluating the sum of these interactions.
Van der Waals
The van der Waals interactions are modeled using a Lennard-Jones 12-6 potential, which captures both attractive (dispersion) and repulsive interactions between non-bonded atoms:
where \(A_{ij}\) and \(B_{ij}\) are parameters for atom pair \(i, j\), and \(r_{ij}\) is the distance between atoms \(i\) and \(j\). These parameters are essential for accurately modeling the attractive and repulsive forces between non-bonded atoms.
The parameters \(A_{ij}\) and \(B_{ij}\) are calculated using the following relationships, which are based on the well depth (\(\epsilon_{ij}\)) and the equilibrium contact distance (\(r_{eq,ij}\)) for each pair of atoms:
- Equilibrium Contact Distance (\(r_{eq}\)): The equilibrium contact distance \(r_{eq}\) is the distance at which the potential energy is at its minimum, indicating the optimal distance for van der Waals interactions.
- Well Depth (\(\epsilon\)): The well depth \(\epsilon\) represents the depth of the potential well, indicating the strength of the attractive interaction at the equilibrium distance.
Fitted parameters
| Atom Type | Van der Waals Radius (Å) | Van der Waals Well Depth (kcal/mol) |
|---|---|---|
| C | 2.00 | 0.15 |
| P | 2.10 | 0.20 |
| N | 1.75 | 0.16 |
| O | 1.60 | 0.20 |
| S | 2.00 | 0.20 |
| H | 1.00 | 0.02 |
The relationships for calculating \(A\) and \(B\) are:
These formulas ensure that the potential has a minimum value of \(-\epsilon_{ij}\) at the equilibrium distance \(r_{eq,ij}\).
Mixing rules
For heterogeneous pairs of atoms (i.e., atoms of different types), the \(\epsilon\) and \(r_{eq}\) values are calculated using combining rules:
-
Geometric Mean for Well Depth:
\[\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}\] -
Arithmetic Mean for Equilibrium Distance:
\[r_{eq,ij} = \frac{r_{eq,ii} + r_{eq,jj}}{2}\]
These combining rules provide a way to estimate the interaction parameters for pairs of different atom types based on the parameters for homogeneous pairs.
Example
Suppose we have two atom types with the following parameters:
- Atom type 1: \(\epsilon_{11} = 0.15\) kcal/mol, \(r_{eq,11} = 2.00\) Å
- Atom type 2: \(\epsilon_{22} = 0.20\) kcal/mol, \(r_{eq,22} = 2.10\) Å
First, compute the combined parameters for the heterogeneous pair (1-2):
-
Well depth (\(\epsilon_{12}\)):
\[ \epsilon_{12} = \sqrt{0.15 \times 0.20} = \sqrt{0.03} \approx 0.173 \text{ kcal/mol} \] -
Equilibrium contact distance (\(r_{eq,12}\)):
\[ r_{eq,12} = \frac{2.00 + 2.10}{2} = 2.05 \text{ Å} \]
Then, calculate the \(A_{12}\) and \(B_{12}\) parameters:
-
\(A_{12}\):
\[ A_{12} = 4 \times 0.173 \times (2.05)^{12} \approx 4 \times 0.173 \times 458.485 \approx 317.43 \text{ kcal·Å}^{12}/\text{mol} \] -
\(B_{12}\):
\[ B_{12} = 4 \times 0.173 \times (2.05)^{6} \approx 4 \times 0.173 \times 114.89 \approx 79.51 \text{ kcal·Å}^{6}/\text{mol} \]
These values of \(A\) and \(B\) can now be used in the Lennard-Jones 12-6 potential to model the van der Waals interactions between atom type 1 and atom type 2.
Hydrogen Bonding
Hydrogen bonds are represented using a 12-10 potential, which emphasizes the directionality of hydrogen bonds:
Here, \(C_{ij}\) and \(D_{ij}\) are hydrogen bond parameters, and \(r_{ij}\) is the distance between the donor and acceptor atoms. These parameters are also derived from equilibrium distances and well depths.
Fitted parameters
| Atom Type | Hydrogen Bond Radius (Å) | Hydrogen Bond Well Depth (kcal/mol) |
|---|---|---|
| N | 1.9 | 5.0 |
| O | 1.9 | 5.0 |
| S | 2.5 | 1.0 |
The parameters \(C_{ij}\) and \(D_{ij}\) are calculated using the following relationships, which are based on the well depth (\(\epsilon_{ij}\)) and the equilibrium contact distance (\(r_{eq,ij}\)) for each pair of hydrogen-bonded atoms. The relationships for calculating \(C\) and \(D\) are:
These formulas ensure that the potential has a minimum value of \(-\epsilon_{ij}\) at the equilibrium distance \(r_{eq,ij}\).
Electrostatics
Electrostatic interactions are calculated using a Coulomb potential with a distance-dependent dielectric function:
where:
- \(q_i\) and \(q_j\) are the partial charges on atoms \(i\) and \(j\), respectively.
- \(r_{ij}\) is the distance between these atoms.
- \(\epsilon(r_{ij})\) is the dielectric constant, which varies as a function of the distance between the interacting charges.
Dielectric function
The electrostatic interactions between the ligand and the protein are a critical component of the scoring function. To accurately simulate these interactions in a solvated environment, AutoDock employs a sigmoidal distance-dependent dielectric function. This approach allows for a more realistic representation of solvent effects, which vary with distance, compared to a constant dielectric model.
In a biological context, the dielectric constant is not uniform; it changes with the distance due to the presence of water and other solvents. To account for this, AutoDock 2.4 uses a sigmoidal distance-dependent dielectric function, which provides a more accurate simulation of solvent effects. The function is defined as follows:
where:
- \(A\) is the dielectric constant close to the protein surface. This value typically represents the lower dielectric environment near the protein surface, reflecting the lower dielectric constant of the protein interior.
- \(B\) is the dielectric constant in the bulk solvent (e.g., water). This is the dielectric constant of the solvent, usually around 78.4 for water at 25°C.
- \(k\) is a parameter that determines the steepness of the sigmoid curve. This parameter controls how quickly the dielectric constant transitions from \(A\) to \(B\) as the distance increases.
- \(r\) is the distance between the interacting atoms.
- \(r_0\) is the midpoint distance where the dielectric constant is halfway between \(A\) and \(B\). This distance defines the point at which the dielectric constant is the average of \(A\) and \(B\).
The sigmoidal function offers several advantages over a constant dielectric model:
- Realistic Solvent Effects: It provides a smooth transition from the low dielectric environment near the protein surface to the high dielectric environment of the bulk solvent, accurately reflecting the physical environment.
- Improved Accuracy: By varying the dielectric constant with distance, it captures the screening effect of the solvent more precisely, leading to better predictions of electrostatic interactions and binding affinities.
- Flexibility: The parameters \(A\), \(B\), \(k\), and \(r_0\) can be adjusted to model different solvent environments and conditions, enhancing the versatility of the docking simulations.
Internal Ligand Energy
The internal energy of the ligand includes contributions from van der Waals interactions, hydrogen bonds, and electrostatic interactions within the ligand itself. This term ensures that the ligand maintains a realistic conformation during docking:
Evaluation
AutoDock2 employs a grid-based method to precalculate and store interaction energies for rapid evaluation during docking. This method involves creating three-dimensional grid maps of pairwise atomic interaction energies around the protein target, allowing for efficient and accurate energy calculations during the docking simulation.
Precalculation of Grid Maps
The process of generating interaction maps in AutoDock2 involves several key steps:
- Protein Target Placement: The protein target is placed in a three-dimensional grid. Each grid point represents a position in space around the protein.
- Probe Atom Evaluation: A probe atom systematically visits each grid point. For each grid point, the pairwise interaction energy of the probe atom with all protein atoms within a cutoff radius of 8 Å is calculated and summed. These interactions include van der Waals forces, hydrogen bonding, and electrostatic interactions.
- Storage of Interaction Energies: The calculated interaction energies for each grid point are stored in separate grid maps, one for each type of atom in the ligand. This results in multiple grid maps, each corresponding to different interaction types such as dispersion/repulsion, hydrogen bonding, and electrostatic potential.
Grid Sampling
Once the grid maps are precalculated, they are loaded into memory at the beginning of the docking process. These grid maps contain precomputed interaction energies at discrete points in a three-dimensional space surrounding the protein target. During the docking simulation, as the ligand explores different conformations and positions, the interaction energy of each ligand atom with the protein is determined by referencing these grid maps.
To evaluate the interaction energy at any given point in space, the ligand's atoms sample the surrounding grid points. This is achieved through trilinear interpolation, a method that interpolates values within a three-dimensional grid. Trilinear interpolation involves considering the eight grid points that form a cube around the point of interest. The energy values at these eight points are combined based on their distances to the actual position of the ligand atom, providing a smooth and continuous estimate of the interaction energy.
This interpolation process allows for rapid and efficient calculation of the ligand's binding energy, as it avoids the need to compute the interaction energy from scratch at every step. Instead, it leverages the precomputed grid values, significantly speeding up the docking simulation while maintaining accuracy. By using trilinear interpolation, AutoDock2 can swiftly evaluate numerous conformations and positions of the ligand, facilitating an efficient exploration of the binding landscape.
Ligand intramolecular contributions
In addition to evaluating the intermolecular interactions between the ligand and the protein, AutoDock2 also considers the intramolecular interactions within the ligand itself. This is crucial for maintaining realistic ligand conformations during the docking process. The intramolecular evaluation ensures that the ligand does not adopt physically improbable conformations that would lead to inaccurate binding predictions.
Intramolecular energies are calculated at each step of the docking simulation to ensure that the ligand maintains a realistic conformation. The steps involved in calculating intramolecular energies are:
- Torsion Angle Adjustments: Changes in the ligand's torsional angles are applied, allowing the ligand to explore different conformations.
- Energy Calculations: The intramolecular van der Waals, hydrogen bonding, and electrostatic energies are calculated based on the current conformation of the ligand.
- Conformation Evaluation: The total intramolecular energy is evaluated to determine if the current conformation is physically plausible. High intramolecular energy indicates unfavorable conformations, which are less likely to be accepted during the simulated annealing process.
By incorporating intramolecular evaluation, AutoDock2 ensures that the ligand adopts conformations that are energetically favorable and realistic. This leads to more accurate predictions of binding modes and affinities, as the ligand's internal strain is minimized. The consideration of intramolecular energies also helps in distinguishing between different binding poses, favoring those that are not only favorable in terms of ligand-protein interactions but also in terms of internal ligand stability.
Calibration and Validation
The empirical coefficients for each term in the scoring function were calibrated using a set of known protein-ligand complexes. Linear regression analysis was performed to determine the coefficients, ensuring that the scoring function accurately reflects observed binding affinities. Validation was achieved by comparing predicted binding energies with experimental data, showing good correlation.
Fitting Details
The fitting of the scoring function involves determining the optimal parameters for each energy term to match experimental binding free energies. This process includes:
- Data Collection: Gathering a dataset of protein-ligand complexes with known binding constants.
- Parameter Estimation: Using regression techniques to estimate parameters such as \(A_{ij}\), \(B_{ij}\), \(C_{ij}\), \(D_{ij}\), and partial charges.
- Optimization: Iteratively adjusting parameters to minimize the difference between predicted and experimental binding free energies.
- Validation: Comparing the fitted scoring function against a separate validation set to ensure its accuracy and generalizability.
-
Morris, G. M., Goodsell, D. S., Huey, R., & Olson, A. J. (1996). Distributed automated docking of flexible ligands to proteins: parallel applications of AutoDock2. Journal of computer-aided molecular design, 10, 293-304. DOI: 10.1007/BF00124499 ↩