Article pubs.acs.org/jcim
Statistical Analysis, Investigation, and Prediction of the Water Positions in the Binding Sites of Proteins Wei Xiao,†,‡ Zenghui He,‡ Meijian Sun,‡ Shiliang Li,‡ and Honglin Li*,†,‡ †
School of Information Science and Engineering, East China University of Science and Technology, Shanghai 200237, China Shanghai Key Laboratory of New Drug Design, School of Pharmacy, East China University of Science and Technology, 130 Meilong Road, Shanghai 200237, China
‡
S Supporting Information *
ABSTRACT: Water molecules play a crucial role in biomolecular associations by mediating a hydrogen bond network or filling spaces with van der Waals interactions. Although current drug design technologies have taken water molecule interactions into account, their applications are still limited to their reliance on either excessive computer resources or a particular potential energy model. Here, we introduce a statistical method that is based on experimentally determined water molecules in the binding sites of high-resolution X-ray crystal structures to predict the potential hydration sites in the binding sites of crystal structures of interest. By clustering and analyzing the various interaction patterns of water molecules in the training data set, we derived a tetrahedronwater-cluster model based on a series of residue group triplets that form feature triangles of different shapes. In the tetrahedralwater-cluster model, a triplet of three polar atoms in the residue group triplet acts as the vertices of the bottom triangle of the tetrahedron, and the water molecule that interacts with these three polar atoms is set as the top vertex of the tetrahedron. By comparing the shapes of the bottom triangles in the training data set with the shape of the triangle in the residue group triplets in the crystal structure of interest, we can identify the bottom triangle that is most similar to the one in the residue group triplet of the crystal structure of interest. According to the tetrahedron-water-cluster model, the hydration site for the residue group triplet in the crystal structure of interest can be predicted based on the height of the tetrahedron that has the most similar bottom triangle in the training data set. A test set containing 193 crystal structures was used to evaluate model performance, and extensive comparison with the recently published program Dowser++ revealed that our model is at least as good at providing an accurate set of the potential hydration sites in crystal structures of interest. improving docking performance8−13 or facilitating receptor− ligand recognition.4 In general, valuable information about the location of water molecules can be obtained not only from primary X-ray crystallography but also from neutron diffraction and nuclear magnetic resonance. In many cases, however, the location of water molecules is uncertain14 or the information is missing.15 Hence, computational approaches allowing efficient prediction of potential hydration sites are of great value.16 Generally speaking, three categories of applications have been developed to date to deal with the problem of hydration sites in the crystal structures of interest, particularly in protein− ligand docking procedures. First, empirically based methods,
1. INTRODUCTION Biomolecules typically exist in an aqueous environment, and interactions with water molecules are of key importance in their structures and functions.1 Aside from mediating the formation of hydrogen bonds between proteins and their ligands,2 water molecules can also modify the shape and flexibility of the ligand binding site in crystal structures.3−6 A study of 392 highresolution complexes retrieved from the Protein Data Bank7 (PDB) found that over 85% of the protein−ligand complexes included at least one water molecule bridging the protein and the ligand. Specifically, the average number of ligand-bound water molecules was found to be 4.6, and 76% of these water molecules were identified as participating in polar interactions with both the protein and the ligand.4 In addition, many studies have shown that considering water molecules may result in © 2017 American Chemical Society
Received: October 12, 2016 Published: June 12, 2017 1517
DOI: 10.1021/acs.jcim.6b00620 J. Chem. Inf. Model. 2017, 57, 1517−1528
Article
Journal of Chemical Information and Modeling such as Consolv,17 WaterScore,18 and WATGEN,19 focus on crystallographically resolving hydration sites and evaluating them with a set of descriptors. These methods can provide quick solutions to distinguish important and potentially conserved water molecules from those that are easily displaced. Of note, these empirically based methods rely heavily on the training data set used. Methods in the second category, such as GRID,20,21 AcquaAlta,22 AQUARIUS,23 SuperStar,24 the force field Fold-X,25 and wPMF,26 are knowledge-based methods that aim to predict hydration sites ab initio and rely on the training data set to determine favorable hydrated sites. The drawback of knowledge-based applications is that they do not directly indicate which predicted water molecules are likely to be displaced by the ligand or remain bound to the protein. Different from above, the third category of physics-based methods, including Dowser++,27 JAWS,28 double decoupling,29 WaterMap,30 WATsite,31 WaterDock,32 and 3D-RISM,33 model explicit water molecules and extract hydration sites from peaks in water density constructed in the protein, or averaged over the location of the water molecule,34 by running molecular dynamics or Monte Carlo simulations. The advantage of this category of methods is that they include entropic effects in the prediction; however, they are timeconsuming to run, especially in scenarios with buried cavities. In addition to the diverse methods described above, many other well-known approaches, such as AutoDock35 and 3D-QSAR,36 full atom models GOLD,12 GLIDE,37,38 and SLIDE,39 and spherical particle models FlexX13,40,41 and FITTED,42 attempt to explore different aspects of hydration properties. While the above three categories of applications confirm the importance of considering water molecules, capturing water molecules in the crystal structure is still challenging.16 In this study, we propose a statistical predictive method that relies on experimentally identified water molecules in the binding sites of high-resolution crystal structures. The interaction patterns of water molecules in the binding sites of crystal structures in the training data set were analyzed in order to develop a more effective molecular docking method that considers water molecules. A modified four-body statistical pseudopotential was also included in the method to analyze the hydrophilic characteristics of residue groups. We propose a tetrahedron-water-cluster model based on a series of residue group triplets that form feature triangles of different shapes. In the tetrahedral-water-cluster model, a triplet of three polar atoms acts as the vertices of the bottom triangle of the tetrahedron, and the water molecule that interacts with these three polar atoms is set as the top vertex of the tetrahedron. By classifying the bottom triangles and the tetrahedrons based on similarity of shape and establishing the relationship between them, the experimentally determined water sites in each type of tetrahedron can be superposed. Furthermore, by comparing the shapes of the bottom triangles in the training data set with the shapes of the triangles in the residue group triplets in the crystal structure of interest, we can identify the bottom triangle that is most similar to the one in the residue group triplet of crystal structure of interest. According to our tetrahedron-water-cluster model, the hydration sites for the residue group triplet in the crystal structure of interest can be predicted based on the height of the tetrahedron that has the most similar bottom triangle in the training data set. The performance of the model in predicting potential hydration sites was evaluated using a test set and compared directly with the program Dowser++. Results showed that the tetrahedron-water-cluster model is at least as
good as Dowser++ at providing an accurate set of the potential hydration sites in crystal structures of interest.
2. EXPERIMENTAL METHODS 2.1. Training Data set. All available protein−ligand complexes in PDB (accessed April 10, 2013) were collected. Selection criteria were as follows: neutron diffraction or nuclear magnetic resonance structures and X-ray crystal structures with a resolution larger than 2.0 Å4,12,17,19,26,43 were removed. Any complexes lacking water molecules in their binding sites were also not retained. The resulting training data set consisted of 2003 complexes (see Table S1 for a comprehensive overview of the complexes included). 2.2. Training Process. Key features of our training process are described below: (I) For each complex, the size of the binding site was defined by the protein atoms and water molecules within 7 Å of any ligand atoms.12,18 (II) On average, each interfacial water molecule had three polar interactions with the surrounding atoms and was less mobile than protein atoms when they had more than three polar interactions.4 Hence, for each water molecule in a binding site of the crystal structure of interest, we only considered at most three surrounding polar atoms within a hydrogen bonding range of 2.5 Å−3.5 Å. Then, according to their sources, we divided the selected polar atoms into three types, that is, ligands (including halogen atoms and metal ions), water molecules, and amino acid residues. (III) To determine if the selected polar atoms from the amino acid residues are exposed to the surfaces of the binding sites, we introduced a program named NACCESS44 to calculate the atomic accessible surface by rolling a probe of a default size around the van der Waals surface and following the coordinates of the center of the probe.45 Only atoms with an atomic accessible surface greater than zero Å2 were retained as interface atoms. The selected polar atoms from the amino acid residues were then combined, forming sets of different residue groups: triplets, two-tuples, and unaries, for further analysis. 2.3. Design of the Prediction Model. As demonstrated in many biophysical investigations, water molecules can make hydrogen bonds with the surrounding atoms and form tetrahedron water clusters.46,47 In order to detect potential hydration sites in the binding sites of the crystal structures of interest, a tetrahedron-water-cluster model (see Figure 1) was proposed based on a series of residue group triplets. In Figure 1, a triplet of three polar atoms serves as the vertices of a bottom triangle of a tetrahedron, and the water molecule that interacts with these atoms is set as the top vertex of a tetrahedron. The shapes of tetrahedrons formed can be classified according to the shape of the bottom triangle and the height of the tetrahedron. First, two ratios are adopted to characterize the shape of the bottom triangle and are used in the similarity calculation. The first ratio is calculated by dividing the longest edge of the bottom triangle a by the shortest edge c in each triplet. The first ratio is then classified into six sectors, i.e., 1.00−1.20, 1.20−1.40, 1.40−1.60, 1.60−1.80, 1.80−2.00, and 2.00−3.00 (Figure 2A), referred to here as fractions. The bottom triangles in each sector are then further divided based on the second ratio, calculated by dividing the middle edge of the same triangle b by the shortest edge c. Subsequently, the bottom triangles are subdivided into more specific subfractions (Figure 2B-2G). 1518
DOI: 10.1021/acs.jcim.6b00620 J. Chem. Inf. Model. 2017, 57, 1517−1528
Article
Journal of Chemical Information and Modeling
retained residue group triplets was calculated. According to the relationship in Table 1, the average height of the tetrahedron in the training data set that has the most similar bottom triangle was set as the distance between the potential hydration site and the center of the triangle in the retained residue group triplet. Two sites were then predicted based on symmetry. Predicted sites whose orientation to the center of the triangle in the retained residue group triplet was the same as that of the experimentally determined water molecule in the crystal structure to the center of the triangle were kept as potential hydration sites. (V) Note that not every experimentally determined water site matched consistently with the predicted potential hydration site in a crystal structure of interest. In practical applications, the predicted sites may collide with the surrounding atoms. In this regard, steric clashes between predicted sites were defined as occurring when the oxygen atom of the water molecule and the surrounding atoms were separated by a distance of less than 75% of the sum of their van der Waals radii,49,50 and the standard deviations of the heights of the tetrahedrons in Table 1 were set as the maximum allowable deviations for predicted sites. In addition, considering the effects of hydrogen bonds and steric hindrance, any potential hydration sites that were predicted within the minimum value of the range of the hydrogen bond interaction, 2.5 Å, were clustered together. 2.5. Analysis of Water-Bound Residue Groups. To investigate the interaction patterns of residue groups, PF values were employed to analyze their preferences using eq 1 (taking the residue group triplets as examples)
Figure 1. Architecture of the tetrahedron-water-cluster model, where a, b, and c stand for the length of the longest, middle, and shortest edges of the bottom triangle, respectively, O represents the center of the bottom triangle, and h represents the height of the tetrahedron.
Having categorized the shapes of the bottom triangles, the shapes of the tetrahedrons are then determined according to the height, h, from the water molecule to the bottom triangle. By observing the shapes of the tetrahedrons in each sector, the sum of the lengths of the three edges (a + b + c) of the bottom triangle in each case is used to establish the relationship between the bottom triangle and the height h of the tetrahedron. The first two columns of Table 1 give the two ratios of the bottom triangles; the sums of the lengths of the three edges of the bottom triangle are divided into different ranges according to the average height of the tetrahedrons. 2.4. Prediction Process. Determining the distribution of water molecules in the binding sites of the crystal structures is a challenging yet important task and is of great practical value for modeling biomolecular structures and their interactions.16 Based on the aforementioned tetrahedron-water-cluster model, we applied the following prediction procedures: (I) The crystal structures of interest for predicting the potential hydration sites in their binding sites were selected using the same selection criteria as for the training data set. A cutoff of 7 Å from any ligand atoms was defined as the binding site. All polar protein atoms in the binding sites of the crystal structures of interest were extracted, and the atomic accessible surface area was calculated to rule out noninterface atoms. (II) Three polar protein atoms from three different amino acid residues in the binding sites were selected to constitute the triplets. Water molecules that interacted with two rather than three amino acid residues were ignored; these data points will be addressed in future work. Residue group triplets for which each edge of the triangles formed was greater than the maximum range of the edge between corresponding residue pairs in the training data set were also filtered. The threshold for corresponding residue pairs that were not in the training data set was set as 7 Å (2 times the maximum range of hydrogen bonds). (III) A preference factor48 (PF) value was introduced to determine whether or not to retain the residue group triplets in crystal structures of interest. (IV) The shapes of triangles formed by the retained residue group triplets were matched to the shapes of triangles formed by residue group triplets in the training data set. First, the two ratios of the triangles formed by the retained residue group triplets were searched to match the shapes of the triangles. Then the sum of the lengths of the three edges of the triangle in
Npattern_i
PF =
∑i Npattern_i N Nres1 N * res2 * res3 ∑j Nresj ∑j Nresj ∑j Nresj
(1)
where Npattern_i designates the occurrence of interaction pattern i in a data set; ∑iNpattern_i designates the total occurrence of all types of interaction patterns in this data set; Nres1, Nres2, and Nres3 are the occurrences of the amino acid residues of the interaction pattern i; and ∑jNresj is the total occurrence of all amino acid residues in this data set. A comprehensive overview of PF values obtained is presented in Tables S2 and S3. Ratios for residue group unaries were calculated by counting the occurrence of the amino acid residue against the total number of the corresponding amino acid residues in the data set (Table S4). To gain a better understanding of the preference of residue group triplets, we introduced four-body statistical pseudopotentials51−55 based on the potential mean force (PMF) in statistical mechanics. Four-body statistical pseudopotentials in the four-body residue neighbors were calculated by means of a statistical geometrical method.53 Statistical potentials are based on the assumption that “contact” energies between amino acid residues in native proteins are related to their observed frequencies in a representative structural database. The detailed calculation process for calculating four-body statistical pseudopotentials is shown in the Supporting Information Appendix 1. The success of the four-body statistical pseudopotentials in predicting the binding affinity of four-body residue neighbors gave us confidence to apply it to study the tetrahedron-watercluster model in order to provide a quantitative measurement of the interaction potentials between the residue group triplets and water molecules. By considering the actual tetrahedron1519
DOI: 10.1021/acs.jcim.6b00620 J. Chem. Inf. Model. 2017, 57, 1517−1528
Article
Journal of Chemical Information and Modeling
Figure 2. Distribution of the ratios used to characterize the shape of the bottom triangle: (A) The bottom triangles are divided into six sectors according to the first ratio. (B−G) Each of the six sectors is then subdivided by the second ratio. The two ratios are then aligned to match the shape of the bottom triangle.
water-clusters formed, we modified the calculation for the fourbody statistical pseudopotentials with the previously calculated
PF values. The water molecule in our tetrahedron-water-cluster model can be treated as one residue in the calculation. Hence, 1520
DOI: 10.1021/acs.jcim.6b00620 J. Chem. Inf. Model. 2017, 57, 1517−1528
1.00−1.05 1.05−1.10 1.10−1.15 1.15−1.20 1.00−1.10 1.10−1.20 1.20−1.30 1.30−1.40 1.00−1.10 1.10−1.20 1.20−1.30 1.30−1.40 1.40−1.50 1.50−1.60 1.00−1.20 1.20−1.40 1.40−1.60 1.60−1.80 1.00−1.20 1.20−1.40 1.40−1.60 1.60−1.80 1.80−2.00 1.00−1.50 1.50−2.00 2.00−2.50 2.50−3.00
1.00−1.20
1521
2.00−3.00
1.80−2.00
1.60−1.80
1.40−1.60
1.20−1.40
ratio2
ratio1
≥13.90 ≥14.00 ≥14.00 ≥13.70 ≥13.54 ≥13.87 ≥13.90 ≥14.00 ≥13.46 ≥13.76 ≥13.70 ≥13.90 ≥13.84 ≥13.74 ≥12.10 ≥13.50 ≥13.54 ≥13.98 ≥11.89 ≥12.60 ≥12.98 ≥13.43 ≥14.00 ≥12.30 ≥13.50 ≥14.10 / 13.90−12.39 14.00−13.00 14.00−13.10 13.70−13.17 13.54−11.66 13.87−12.40 13.90−12.67 14.00−13.17 13.46−11.87 13.76−12.23 13.70−12.20 13.90−12.50 13.84−12.67 13.74−12.48 12.10−10.87 13.50−12.30 13.54−12.60 13.98−13.37
