Hybrid Molecular Dynamics for Elucidating Cooperativity Between Halogen Bond and Water Molecules During the Interaction of p53-Y220C and the PhiKan5196 Complex

Introduction

Halogen bonds have become a trending subject for studies aiming at contributing to rational drug design (Lu et al., 2012). These bonds are an important non-covalent interaction in a protein–ligand binding process (Lu et al., 2009) and are similar to hydrogen bonds (Metrangolo et al., 2005). Although the contributions of different parameters vary, the electrostatic property, dispersion, polarization, and charge transfer effects are all known to contribute to the halogen bond interaction (Hanna, 1968; Riley et al., 2013; Esrafili et al., 2017). Previously, a σ-hole model was developed to characterize halogen bonds, wherein a crown of positive charge on the outer side of a halogen atom could form a halogen bond and a belt of negative charge surrounding the halogen atom could form a hydrogen bond (Clark, 2013, 2017; Politzer et al., 2013, 2017). Therefore, the co-existence of both halogen as well as hydrogen bonds due to the presence of a σ-hole in a halogen atom is one of the important characteristics of halogen bonds (Li et al., 2008; Alkorta et al., 2010; Zhao et al., 2011; Yan et al., 2014; Zhou et al., 2016; Molina et al., 2017; Ciancaleoni, 2018).

In recent years, numerous studies have reported the positive cooperativity between halogen bonds and hydrogen bonds (Li et al., 2008; Grabowski, 2013; Wu et al., 2013; Esrafili and Mousavian, 2017; Esrafili and Vakili, 2017; Carlsson et al., 2018) as well as the effects of the substituents on the cooperativity of halogen bonds (Solimannejad et al., 2013). In view of the potential application prospects of halogen bonds in the field of drug design, Adasme-Carreno et al. (2016) performed calculations on 126 complexes of drug-like molecules; consequently, a positive cooperativity effect was observed in N-methylacetamide complexes with di-, tri-, and tetrafluoro-iodobenzenes, which led to the strengthening of halogen bonds. The water molecules are ubiquitous in protein complexes; particularly, those in close proximity to halogen bonds exhibited reciprocal interactions (De Santis et al., 2018). For example, this was observed during the enhancement of anion recognition (Langton et al., 2014) and strengthening binding affinity for CK2 ligands containing halogen atoms (Deepa, 2017; Deepa et al., 2018).

Robinson et al. (2016) employed explicit solvent molecular dynamics and solvent thermodynamic analysis in solvent-exposed ligand modifications to improve ligand selectivity, which was determined by water thermodynamics. Breiten et al. (2013) reported that the water molecules fill the active site of the protein and surround the ligand and are as important as the ligands that possess biomolecular recognition and determine the thermodynamics of bindings. Therefore, we are interested in validating the existence of the hydrogen bonds that are formed by the water molecules with nearby halogen atoms, elucidating the cooperative effects between the hydrogen and the halogen bonds, and investigating the effects of the thermodynamic behavior of water molecules on the stability of the halogen bonds.

In order to be able to reproduce the σ-hole characteristics in molecular modeling, several force fields for halogen bonds have been developed, such as extra-point charge (Ibrahim, 2011), force field for biological halogen bonds (Carter et al., 2012), OPLS-AAx and OPLS/CM1Ax (Jorgensen and Schyman, 2012), and multipole- (El Hage et al., 2016) and quantum-based generalized empirical potential (Santos et al., 2017). The polarizable ellipsoidal force field (Du et al., 2013; Liu et al., 2015) represented the anisotropic charge distribution with the combination of a negatively charged sphere and a positively charged ellipsoid. This model could correctly reproduce the potential energy surface of halogen bonds at the MP2 level. However, the performances of these force fields have not been completely tested with respect to the cooperative systems of hydrogen and halogen bonds; this is partially because the dynamics governing the cooperativity of the hydrogen and the halogen bonds are still unknown (Evangelisti et al., 2011; Chen et al., 2018; Lu et al., 2019).

The utilization of the Y220C mutant of the p53 tumor suppressor with a ligand complex is a good approach to study the dynamics of the cooperativity of the hydrogen and the halogen bonds (Yang et al., 2019). In several cancers, p53 is inactivated by the direct mutation or perturbation of its associated pathways. The Y220C mutant destabilized p53 and generated a drugable surface crevice, which has led to approximately 75,000 new cancer cases per year. Wilcken et al. (2012) developed a strategy that utilized halogen bonds to stabilize the mutant with small molecules, which use halogen-enriched fragment libraries as the starting points. After a careful screening of their solved structures, we selected the Y220C–PhiKan5196 complex (PDB ID: 4AGQ) as the model. As shown in the crystal structure (Figure 1), a halogen bond was formed between the iodine atom of PhiKan5196 (P96) and the carbonyl oxygen of Leu50. Three water molecules were observed to be around the iodine atom, which has not been addressed previously. Therefore, it is necessary to explore the existence of hydrogen bonds that have been formed by these water molecules by incorporating iodine atoms and observing their cooperative effects with the existing halogen bonds.

FIGURE 1

Figure 1. Crystal structure of the Y220C mutant of the p53 core domain bound to a stabilized small molecule PhiKan5196 (PDB ID: 4AGQ), in which the protein is shown as its secondary structure and the ligand is illustrated by the ball-and-stick model. The major interactions of P96 with Leu50 and three water molecules are indicated.

In this work, we employed a method governed by hybrid quantum mechanical molecular mechanical molecular dynamics (QM/MM MD) to simulate the Y220C–PhiKan5196 complex; this method considered the polarization effects of the halogen and the hydrogen bonds and was successful in striking a balance between computational demand and accuracy (Ke et al., 2009; Ahmadi et al., 2018; Awoonor-Williams and Rowley, 2018; Steinmann et al., 2018; Visscher et al., 2018). For comparison purposes, we designed two models, wherein one comprised three water molecules (model 2) and the other was rendered devoid of any water molecule by the removal of three water molecules in the QM region (model 1) (see “MATERIALS AND METHODS”). We compared the dynamic performance of halogen bonds between the two models and concluded that the hydrogen bonds formed by the three spatial correlated water molecules and the iodine atom were able to not only synergize but also strengthen the halogen bond. Furthermore, the contribution of each hydrogen bond to the halogen bond was observed to be different.

Materials and Methods

Computational Methods and Model Setup

The workflow of this study is shown in Figure 2. Starting from the crystal structure, we first carried out the classical molecular dynamic simulation to render the system equalized; two models were built for the QM/MM MD simulation. The trajectories were then classified for further analysis. The specific details of the method are described in the following discussion.

FIGURE 2

Figure 2. Schematic illustration of the workflow of the system setup and the QM/MM MD simulation.

The preliminary structure of the classical MD simulation was selected from the crystal structure (PDB ID: 4AGQ); this structure was of the p53 core domain mutant, Y220C, which was bound to the stabilizing small molecule of PhiKan5196 (P96). It comprised two biological units: chain A and chain B; chain A was selected and chain B was manually deleted. The AMBER ff14SB force field (A. Maier et al., 2015) was used and the parameters of the ligands were derived from the generalized AMBER force field (Wang et al., 2004). The unbonded model (Pang, 1999) was applied to zinc and the surrounding amino acids. The complex was solvated in a cubical TIP3P water box with periodic boundary conditions, with a margin of at least 10Å from any edge of the box to any atom of the complex. The solvated systems were neutralized by adding Na⁺ and Cl⁻ ions and the final concentration for the solvated system was observed to be 0.15 M.

The system energy was minimized through conjugate gradient minimization. After that, the system was heated to 300 K at a constant volume over 2 ns. A 5-ns MD simulation was extended as production. An integral time step of 2 fs was used, and electrostatics were treated by the particle-mesh Ewald method (Essmann et al., 1995) with a 12-Å Coulomb interaction distance cutoff; further, the cutoff for the van der Waals interactions was set to 12 Å. The simulation was performed using the NAMD package (Phillips et al., 2005).

QM/MM MD Simulation

The frame of the production trajectory that had the smallest root mean square deviation value with respect to the crystal structure was selected for the QM/MM MD simulation. We set up two models in this work. For model 2, the QM region was composed of the following atoms: all the atoms of Cys134, all side chain atoms and some of the backbone atoms of Val52, Leu50, Thr135, and Pro128, all the backbone atoms of Trp51, some of the backbone atoms of Asp133 and Pro127, W1 to W3, and a part of PhiKan5196 (P96) (Figure 3). Some of the backbone atoms were included to obtain valid QM/MM cuts (Murphy et al., 2000). The QM region was saturated with capping hydrogen atoms. The QM region of model 1 was the same as that of model 2, except for the absence of W1 to W3 (see Figure S1). In other words, there are no water molecules in the QM region of model 1. The QM region of model 2 comprised 98 atoms, while the QM region of model 1 comprised 89 atoms; the remaining MM parts were composed of 32,375 atoms in both models.

FIGURE 3

Figure 3. Structure of the active site in p53 Y220C that is used as the QM region in QM/MM MD simulations. The protein is shown as a gray cartoon representation and selected residues are highlighted as stick models. The halogen bond between the iodine and the carbonyl oxygen of Leu50 is indicated by a yellow dash line. Additional polar interactions between iodine and water are shown as orange dash lines. Water 1 to water 3 were named as W1 to W3 in this work.

The MM region of both systems was described using the Amber ff14SB force field in MM calculations. In the QM/MM implementation, E_MM is calculated classically from the MM atom positions, whereas H_QM is evaluated using the B3LYP and LANL2DZ bases that are set for all the atom types. The AMBER16 (Case et al., 2016) integrated Gaussian 09 (Frisch et al., 2009) package was used to carry out all the QM/MM MD simulations.

The QM/MM minimization process included 500 steps of steepest descent, which was followed by 200 steps of conjugate gradient minimization. In this process, any atom within 10Å of the ligand was allowed to move freely and the remaining part of the system was kept restrained. Subsequently, the entire system was simulated for 2.5 ps by using the NPT ensemble with a time step of 0.5 fs, and the temperature was regulated using a Langevin thermostat at 300 K (Loncharich et al., 1992). The QM part of the system was not constrained by the SHAKE algorithm (Ryckaert et al., 1977). The non-bonded interactions were calculated using a cutoff of 10Å. We carried out six identical simulations for each model; therefore, a total of 3,000 snapshots (500 for each trajectory) for each model were used for subsequent analysis.

Clustering Analysis

The snapshots of each model were classified into several clusters using the K-means clustering method (Arthur and Vassilvitskii, 2007). As a simple and robust algorithm, the K-means method is suitable for identifying significant structural features. An arbitrarily low number of clusters may lead to important information being missed, whereas too many clusters would be unnecessarily redundant. Thus, it is necessary to set an appropriate number of premeditated clusters. The number of clusters was set to three after testing the cluster numbers of two to 10 (i.e., k = 2~10); this was aimed at achieving a better structural similarity within each cluster. Five variables were used to serve as clustering descriptors, namely: length and angle of the halogen bond [symbolized by d(I⋯O) and ∠(C-I⋯O)], length of the three hydrogen bonds formed by three water molecules [symbolized by d(I⋯O_w1), d(I⋯O_w2), and d(I⋯O_w3)]. To ensure that the descriptors were comparable to each other, we normalized the five kinds of descriptors. This protocol has been implemented with the scikit-learn (Pedregosa et al., 2011) module of Python.

Interaction Energy Calculation

We used the Moller—Plesset second order (MP2) (Headgordon et al., 1988) method with an aug-cc-pVTZ-PP pseudo-potential basis set for iodine atoms and the aug-cc-pVTZ basis set for hydrogen, carbon, nitrogen, and oxygen atoms; this was aimed at calculating the interaction energy. The calculated interaction energies (ΔE_Corr) of the complexes were corrected for the basis set as superposition error that employed the counterpoise method of Boys and Bernardi (1970). All ab initio calculations were carried out using the Gaussian 09 package (Frisch et al., 2009).

For convenience, the hydrogen-bonded complexes are denoted by A–B, and the halogen-bonded complexes are denoted by B–C; A–B–C represents the entire complexes. Here A represents the three water molecules around the iodine atom, B stands for the ligand fragment, and C denotes the leucine residue fragment. The halogen bonding interaction energies are calculated here according to the method employed in earlier studies pertaining to cooperativity effects (Zhao et al., 2011).

$\begin{array}{ll}{\Delta E}_{BC}=E_{\text{ABC}}-({\text{E}}_{\text{C}}+E_{AB})& (1)\end{array}$

where E_ABC is the total energy of the entire system and E_A, E_C and _{EBC, EAB} are the energies of the monomeric and the binary systems. It should be noted that the interaction of B and C (ΔE_BC) also includes the interaction energy between fragments A and C. In this work, we consider ΔE_BC as a small constant since the distance of A and C is relatively long.

Correlation Coefficient Matrix Analysis

To evaluate the effect of the three water molecules on the halogen bond, correlation coefficients matrices (Dawson and Trapp, 2004) were calculated between two pairs of the three hydrogen bond lengths, based on Eq. (2):

$\begin{array}{ll}{R}_{ij}=\frac{C_{ij}}{\sqrt{C_{ii}*C_{jj}}}& (2)\end{array}$

$\begin{array}{ll}{C}_{ij}=\frac{1}{F}{\displaystyle \sum _{i=1}^F}(ri-{\overline{r}}_i)(r_j-{\overline{r}}_j)& (3)\end{array}$

In Eq. (2), C is the covariance matrix calculated between two variables, as shown in Eq. 3. F is the number of the total frames (or structures) and r is the hydrogen bond length. The values of R are between−1 and 1, both inclusive. Udovicic et al. (2007) stated that the absolute values of correlation coefficients from 0 to 0.25 indicate an absence of correlation, those from 0.25 to 0.50 indicate poor correlation, those from 0.50 to 0.75 indicate moderate to good correlation, and those above 0.75 indicate good to excellent correlation.

Results

Comparison of QM/MM MD Trajectories of Two Models

In order to compare the difference in the halogen bonds between the two models, we calculated the distribution of 3,000 frames of the I⋯O distance and the C-I⋯O angle for each model, as shown in Figures 4A,B. Evidently, the mean values of the I⋯O distance and the C-I⋯O angle (bond length and bond angle of halogen bond) are different. The average bond length was shortened from 3.70 to 3.45Å (Figure 4A); meanwhile, the bond angle was widened from 156 to 165° (Figure 4B); both of these values suggested that the three water molecules tend to enhance the halogen bond between the P96 ligand and the p53 protein (Cavallo et al., 2016). In addition to the overall distribution, the average value of the halogen bond length and the bond angle of each trajectory (Tables S1, S2) can also support the result.

FIGURE 4

Figure 4. Normalized distributions of halogen bond length and angle. (A) Distribution of halogen bond length that is labeled by I⋯O distances. (B) Distribution of halogen bond angle that is labeled by C-I⋯O angles. The color codes are red for model 1 and black for model 2.

Clustering of the Trajectories and the Roles of Water Molecules

In this section, we focused on analyzing the function of the three water molecules. Therefore, the majority of the analyses were based on the trajectories of model 2. The K-mean algorithm was applied with the aim of partitioning six trajectories of model 2 into distinct clusters by their geometries. The final clustering results of the 3,000 frames of these dynamic trajectories are listed in Table 1, and the representative structures of the three clusters are shown in Figure 5. We classified the frames of the trajectories into three distinct classes, after comparing its results with trials in different numbers of classes.

TABLE 1

Table 1. Average values and standard deviations of key parameters and number of structures of three clusters.

FIGURE 5

Figure 5. Representation of the clustering results. (A) Superposition of the three representative structures of each cluster in a ball-and-stick model. It is worth noting that the iodine atom, the oxygen atoms in the three water molecules, and the oxygen atom in the protein that formed a halogen bond are seen in the opaque material. The oxygen atoms in the water molecules in each cluster are shown as points. Color code: red, cluster A; gray, cluster B; blue, cluster C. The simplified representative structure of cluster A (B), cluster B (C), and cluster C (D), saturating the broken bond with hydrogen atoms. The values of the five clustering descriptors are indicated.

Several geometric characteristics of these three similarly sized clusters can be seen in Table 1. Firstly, in terms of the length and the angle of the halogen bond, the I⋯O distance and the C-I⋯O angle are similar for cluster A and B; this implies that clusters A and B have a similar halogen bond strength. The I⋯O distance is increased by 0.4Å for both average values (Table 1). The representative structure (Figures 5B,D) and the C-I⋯O angle are decreased by 7.7 and 9°, respectively, on average (Figures 5B,D) in cluster C, as compared to those in cluster A, which indicated that the average bond strength of cluster C is weaker than that of both A and B. Secondly, in terms of the length of the hydrogen bond, in cluster A W1 is significantly closer to the iodine atom than that in W2 and in W3. The main difference between clusters A and B is the I⋯O_w3 distance (hydrogen bond length), which was reduced by 0.6Å on average (Table 1 and Figures 6B,C) and by 0.4° in the representative structure (Figures 5B,C) from clusters A to B. The main difference between clusters A and C is reflected on the distances of I⋯O_w3 and I⋯O_w2, which were reduced by 0.7 and 0.5Å on average (Table 1 and Figures 6B,D) and by 0.5Å in the representative structure (Figures 5B,D) from clusters A to C; therefore, it is reasonable to consider cluster A as a representation of W1 near the halogen atom, cluster B as a representation of W3 near the halogen atom, and cluster C as the representation of W2 and W3 near the iodine atom.

FIGURE 6

Figure 6. Normalized distribution of the locations of the three water molecules near the iodine atom for all the six trajectories of model 2. (A) Overall distributions without clustering, (B) cluster A, (C) cluster B, and (D) cluster C. Color code: red, W1; black, W2; blue, W3.

Interestingly, the cluster with all the three water molecules observed to be simultaneously close to the halogen atom was not observed. In principle, the polarization of multiple hydrogen bonds should facilitate the enhancement of the halogen bond, which is determined by the nature of the σ-hole. Therefore, we further analyzed the function of the three water molecules from the perspective of energy.

The halogen bond interaction energy contributions of different amounts of water molecules in each cluster are shown in Table 2. In order to calculate the interaction energy, the protein environment was omitted and only the simplified QM region (Figure 5) was considered. The contribution of each water molecule was calculated by removing other water molecules in the representative structures. When adding W1, W2, and W3 consecutively back to their corresponding locations, the contribution of each addition to the interaction energy was −1.89, −0.01, and −0.75 kcal/mol in cluster A, −0.97, −0.34, and −1.24 kcal/mol in cluster B, and −1.64, −0.13, and −0.74 kcal/mol in cluster C, respectively (Table 2). Hence, the addition of water molecules can indeed enhance the strength of the halogen bond. This observation is in accordance with those obtained in the previous studies by Carlsson et al. (2018), which stated that a hydrogen bond formed with a halogen atom can enhance the halogen bond by increasing the positive electrostatic potential of σ-hole. Moreover, it can be seen in Table 2 that the contribution of each water molecule to the halogen bond in different clusters is not consistent: the largest one in clusters A and C is from W1, whereas that in cluster B is from W3.

TABLE 2

Table 2. Halogen bond interaction energy ΔE (kcal/mol) of the representative structure for each cluster in the presence of different water molecules and the change in the halogen bond interaction energy, ΔE (kcal/mol), due to the presence of additional water molecules (at the MP2/aug-cc-pVTZ-PP level of theory for I and at the MP2/aug-cc-pVTZ level of theory for H, C, N, and O).

The observations performed with regard to energy are correlated with the previous distance analyses; however, these are also different in some aspects. In cluster C, the distance between the three water molecules and the halogen atom are similar, while their energy contributions vary from −0.13 to −1.64 kcal/mol (Table 2). Therefore, we concluded that the functions of the three water molecules in enhancing the halogen bond are not equivalent, and the order of their general contribution is given by W1 > W3 > W2; this can be reflected in the overall locations of the three water molecules in Figure 6A. This is because the position of W1 depicts that it can form a double hydrogen bond not only with the hydroxyl oxygen on the benzene ring but also with iodine atom, which in turn is stabilized by hydrogen bonds with the main-chain atoms of Val52 and Asp133 (Wilcken et al., 2012). Nevertheless, this order is not absolute. In nearly one-third of all the structures from cluster B, the contribution of W3 is observed to be greater than that of W1, which is also reflected in Figure 6C.

Relationship of the Water Molecules

In this section, we performed various analyses to elucidate the dynamic polarization effect exerted by the three water molecules; this is achieved by analyzing the representative trajectory. The hydrogen bond lengths formed by the iodine atom and the three water molecules for each trajectory are shown in Table S1. The existence of a correlation between the movements of the three water molecules was first validated. The correlation coefficients (Eq. 3) between the pairwise hydrogen bond length of the three water molecules were calculated for 500 frames in each trajectory (Table 3).

TABLE 3

Table 3. Cross-correlation coefficients between two pairs of three hydrogen bond lengths in six trajectories.

The absolute values of the pairwise correlation coefficients of the three water molecules are not >0.6, which implied that their motions are weakly correlated. In addition, the correlation between the hydrogen bonds varies with the trajectory. (Parra and Ohlssen, 2008) confirmed that two bifurcated H bonds observed to be simultaneously present in the same molecule were found to cooperatively reinforce each other from geometrical, energetic, and topological analyses, although according to the results obtained in the study by Hogan, there is no such situation wherein two water molecules led to the simultaneous formation of hydrogen bonds with halogen atoms such as W2 and W3 (Hogan and van Mourik, 2019). Structurally, in our system, W2 and W3 are on the same side of the benzene ring, and a triangular interaction network is formed between the two water molecules and the iodine atom; this facilitates the formation of this structure. Therefore, the relationship between multiple water molecules around the halogen atom is complicated and is also related to the surrounding chemical environment that needs to be studied further.

Effect of Protein-Rich Environment on the Location of Water Molecules

Based on the previous analysis, it can be established that there is no such situation wherein three water molecules are simultaneously near the iodine atom in a protein-rich environment. For comparison, we optimized the QM structure in a protein-rich environment (QM/MM-minimized) as well as under vacuum (Figure 7). It was discovered that the three water molecules can be near the iodine atom at the same time under vacuum. In the protein environment, except for W1, the distances between the remaining two water molecules and the iodine atom are observed to be larger than those in vacuum. Therefore, the protein environment could weaken the hydrogen bonds formed by three water molecules with the incorporation of an iodine atom. Additionally, compared to the structure obtained in the protein-rich environment (Figure 7A), we found that the position of W2 in vacuum remains unchanged; however, W3 is observed to rearrange to the back of the benzene ring that is occupied by Thr135 and vice versa (Figure 7A). Thus, the illustrated steric constraint effects and the weakening of hydrogen bonds caused by the protein environment cause a reduction in the strength of the halogen bond in the protein-rich environment (the interaction energy of the halogen bond in vacuum is −11.91 kcal/mol and that in protein is −6.78 kcal/mol).

FIGURE 7

Figure 7. Changes in the hydrogen bonds formed by the three water molecules and the iodine atom in a protein-rich environment and under vacuum. (A) The QM/MM-minimized structure of model 2 in the protein environment. (B) The optimized simplified QM structure in vacuum (at the MP2/aug-cc-pVTZ-PP level of theory for I and at the MP2/aug-cc-pVTZ level of theory for H, C, N, and O).

Consequently, we wonder whether the presence of water molecules also strengthen the halogen bonds in other halides. Therefore, we carried out a simple treatment to replace I in the optimized and simplified QM structure for both environments (protein-rich environment and under vacuum) with F, Cl, and Br. For the structure obtained after substitution, the interaction energies of halogen bonds with and without the three water molecules in both environments were calculated. It can be seen in Table 4 that other halogenated compounds exhibit the same behavior as that observed in this work, with respect to the halogen bonds; the protein-rich environment will weaken the halogen bonds, and the presence of water molecules under the same environmental condition will strengthen the halogen bonds.

TABLE 4

Table 4. Halogen bond interaction energy ΔE (kcal/mol) of halides in the protein-rich environment and in vacuum (at the MP2/aug-cc-pVTZ-PP level of theory for I and at the MP2/aug-cc-pVTZ level of theory for H, C, N, O, F, Cl, and Br).

Discussion

Herein, the role of the three water molecules around the halogen bond in the Y220C–PhiKan5196 complex of p53 protein was analyzed. In order to elucidate the polarization effect of the hydrogen bond on the halogen atom, we adopted the QM/MM MD simulation method. By comparing this method to a model that was free of water, it was revealed that the strength of the halogen bond in the model with three water molecules was stronger. This confirmed that the water molecules formed weak hydrogen bonds with halogen atoms and augmented the halogen bonds by cooperation.

The clustering and energy analyses indicated that the roles of the three water molecules were not similar, and the contribution to the cooperativity of hydrogen and halogen bonds was primarily in the following order: W1 > W3 > W2. At the same time, the protein-rich environment played an important role in the location of the water molecules.

We believe that these results provide insights into the dynamic cooperation of the hydrogen and the halogen bonds, which can facilitate the development of the force field for halogen bonds. Meanwhile, these results are also useful for rational drug design.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

JG conceived and designed the study. TD, JG, and XW wrote the manuscript. TD, HP, and XH performed QM/MM MD simulation and data analysis.

Funding

This work was supported by the National Key R&D Program of China (Grant No. 2017YFB0203405), the National Natural Science Foundation of China (No. 21873034), the Fundamental Research Funds for the Central Universities (Project 2662018JC027), and the Huazhong Agricultural University Scientific and Technological Self-innovation Foundation (Program No. 2015RC008).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2020.00344/full#supplementary-material

References

Adasme-Carreno, F., Munoz-Gutierrez, C., and Alzate-Morales, J. H. (2016). Halogen bonding in drug-like molecules: a computational and systematic study of the substituent effect. Rsc Adv. 6, 61837–61847. doi: 10.1039/c6ra14837a

CrossRef Full Text | Google Scholar

Ahmadi, S., Herrera, L. B., Chehelamirani, M., Hostas, J., Jalife, S., and Salahub, D. R. (2018). Multiscale modeling of enzymes: QM-cluster, QM/MM, and QM/MM/MD: a tutorial review. Int. J. Quant. Chem. 118:e25558. doi: 10.1002/qua.25558

CrossRef Full Text | Google Scholar

Alkorta, I., Blanco, F., Deyà, P. M., Elguero, J., Estarellas, C., Frontera, A., et al. (2010). Cooperativity in multiple unusual weak bonds. Theor. Chem. Acc. 126, 1–14. doi: 10.1007/s00214-009-0690-1

CrossRef Full Text | Google Scholar

Arthur, D., and Vassilvitskii, S. (2007). “k-means++: the advantages of careful seeding,” in: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms (New Orleans, LA: Society for Industrial and Applied Mathematics), 1027–1035.

Google Scholar

Awoonor-Williams, E., and Rowley, C. N. (2018). The hydration structure of methylthiolate from QM/MM molecular dynamics. J. Chem. Phys. 149:045103. doi: 10.1063/1.5038010

PubMed Abstract | CrossRef Full Text | Google Scholar

Boys, S. F., and Bernardi, F. D. (1970). The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. Mol. Phys. 19, 553–566. doi: 10.1080/00268977000101561

CrossRef Full Text | Google Scholar

Breiten, B., Lockett, M. R., Sherman, W., Fujita, S., Al-Sayah, M., Lange, H., et al. (2013). Water networks contribute to enthalpy/entropy compensation in protein–ligand binding. J. Am. Chem. Soc. 135, 15579–15584. doi: 10.1021/ja4075776

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlsson, A. C. C., Scholfield, M. R., Rowe, R. K., Ford, M. C., Alexander, A. T., Mehl, R. A., et al. (2018). Increasing enzyme stability and activity through hydrogen bond-enhanced halogen bonds. Biochemistry 57, 4135–4147. doi: 10.1021/acs.biochem.8b00603

PubMed Abstract | CrossRef Full Text | Google Scholar

Carter, M., Rappé, A. K., and Ho, P. S. (2012). Scalable anisotropic shape and electrostatic models for biological bromine halogen bonds. J. Chem. Theory Comput. 8, 2461–2473. doi: 10.1021/ct3001969

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Betz, R. M., Cerutti, D. S., Cheatham, T. E. III., Darden, T. A., Duke, R. E., et al. (2016). AMBER 2016. (San Francisco CA: University of California).

Google Scholar

Cavallo, G., Metrangolo, P., Milani, R., Pilati, T., Priimagi, A., Resnati, G., et al. (2016). The halogen bond. Chem. Rev. 116, 2478–2601. doi: 10.1021/acs.chemrev.5b00484

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J. H., Wang, J., Zheng, Y., Feng, G., and Gou, Q. (2018). Halogen bond in the water adduct of chloropentafluoroethane revealed by rotational spectroscopy. J. Chem. Phys. 149:154307. doi: 10.1063/1.5052132

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciancaleoni, G. (2018). Cooperativity between hydrogen- and halogen bonds: the case of selenourea. Phys. Chem. Chem. Phys. 20, 8506–8514. doi: 10.1039/c8cp00353j

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, T. (2013). σ-Holes. wiley interdisciplinary reviews. Comput. Mol. Sci. 3, 13–20. doi: 10.1002/wcms.1113

CrossRef Full Text | Google Scholar

Clark, T. (2017). Halogen bonds and sigma-holes. Faraday Discuss. 203, 9–27. doi: 10.1039/c7fd00058h

PubMed Abstract | CrossRef Full Text | Google Scholar

Dawson, B., and Trapp, R. (2004). Basic and Clinical Biostatistics. New York, NY: Lange MedicalBooks/McGraw-Hill

Google Scholar

De Santis, M., Nunzi, F., Cesario, D., Belpassi, L., Tarantelli, F., Cappelletti, D., et al. (2018). Cooperative role of halogen and hydrogen bonding in the stabilization of water adducts with apolar molecules. New J. Chem. 42, 10603–10614. doi: 10.1039/c8nj00422f

CrossRef Full Text | Google Scholar

Deepa, P. (2017). An Overview of the factors influencing CK2 ligands and the impact of crystal waters: a theoretical study. Crystal Growth Design 17, 1299–1315. doi: 10.1021/acs.cgd.6b01711

CrossRef Full Text | Google Scholar

Deepa, P., Pandiyan, B. V., and Kolandaivel, P. (2018). Does the presence of water clusters induce the binding affinity of CK2 halogen ligands?: A quantum chemical perspective study. Int. J. Quant. Chem. 118:25644. doi: 10.1002/qua.25644

CrossRef Full Text | Google Scholar

Du, L. K., Gao, J., Bi, F. Z., Wang, L. L., and Liu, C. B. (2013). A polarizable ellipsoidal force field for halogen bonds. J. Comput. Chem. 34, 2032–2040. doi: 10.1002/jcc.23362

PubMed Abstract | CrossRef Full Text | Google Scholar

El Hage, K., Bereau, T., Jakobsen, S., and Meuwly, M. (2016). Impact of quadrupolar electrostatics on atoms adjacent to the sigma-hole in condensed-phase simulations. J. Chem. Theory. Comput. 12, 3008–3019. doi: 10.1021/acs.jctc.6b00202

PubMed Abstract | CrossRef Full Text | Google Scholar

Esrafili, M. D., Moharnrnadian-Sabet, F., and Vessally, E. (2017). An ab initio study on substituent and cooperative effects in bifurcated fluorine bonds. Mol. Phys. 115, 278–287. doi: 10.1080/00268976.2016.1257828

CrossRef Full Text | Google Scholar

Esrafili, M. D., and Mousavian, P. (2017). Unusual cooperativity effects between halogen bond and donor-acceptor interactions: the role of orbital interaction. Chem. Phys. Lett. 678, 275–282. doi: 10.1016/j.cplett.2017.04.052

CrossRef Full Text | Google Scholar

Esrafili, M. D., and Vakili, M. (2017). The effect of hydrogen-bonding cooperativity on the strength and properties of sigma-hole interactions: an ab initio study. Mol. Phys. 115, 913–924. doi: 10.1080/00268976.2017.1292013

CrossRef Full Text | Google Scholar

Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A smooth particle mesh ewald method. J. Chem. Phys. 103, 8577–8593. doi: 10.1063/1.470117

CrossRef Full Text | Google Scholar

Evangelisti, L., Feng, G., Ecija, P., Cocinero, E. J., Castano, F., and Caminati, W. (2011). The halogen bond and internal dynamics in the molecular complex of CF3Cl and H2O. Angew. Chem. Int. Ed. Engl. 50, 7807–7810. doi: 10.1002/anie.201102544

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2009). Gaussian 09. Wallingford, CT: Gaussian, Inc.

Google Scholar

Grabowski, J. S. (2013). Cooperativity of hydrogen and halogen bond interactions. Theor. Chem. Acc. 132:1347. doi: 10.1007/s00214-013-1347-7

CrossRef Full Text | Google Scholar

Hanna, M. W. (1968). Bonding in donor-acceptor complexes. electrostatic contributions to ground-state properties of benzene-halogen complexes. J. Am. Chem. Soc. 90, 285–291. doi: 10.1021/ja01004a014

CrossRef Full Text | Google Scholar

Headgordon, M., Pople, J. A., and Frisch, M. J. (1988). Mp2 energy evaluation by direct methods. Chem. Phys. Lett. 153, 503–506. doi: 10.1016/0009-2614(88)85250-3

CrossRef Full Text | Google Scholar

Hogan, S. W. L., and van Mourik, T. (2019). Halogen bonding in mono- and dihydrated halobenzene. J. Comput. Chem. 40, 554–561. doi: 10.1002/jcc.25733

PubMed Abstract | CrossRef Full Text | Google Scholar

Ibrahim, M. A. A. (2011). Molecular mechanical study of halogen bonding in drug discovery. J. Comput. Chem. 32, 2564–2574. doi: 10.1002/jcc.21836

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., and Schyman, P. (2012). Treatment of halogen bonding in the OPLS-AA force field: application to potent anti-HIV agents. J. Chem. Theory Comput. 8, 3895–3901. doi: 10.1021/ct300180w

PubMed Abstract | CrossRef Full Text | Google Scholar

Ke, Z. H., Wang, S. L., Xie, D. Q., and Zhang, Y. K. (2009). Born-oppenheimer ab initio QM/MM molecular dynamics simulations of the hydrolysis reaction catalyzed by protein arginine deiminase 4. J. Phys. Chem. B 113, 16705–16710. doi: 10.1021/jp9080614

PubMed Abstract | CrossRef Full Text | Google Scholar

Langton, M. J., Robinson, S. W., Marques, I., Felix, V., and Beer, P. D. (2014). Halogen bonding in water results in enhanced anion recognition in acyclic and rotaxane hosts. Nat. Chem. 6, 1039–1043. doi: 10.1038/nchem.2111

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Lin, Q., Li, W., Cheng, J., Gong, B., and Sun, J. (2008). Cooperativity between the halogen bond and the hydrogen bond in H₃N···XY···HF complexes (X, Y=F, Cl, Br). Eur. Chem. Soc. 9, 2265–2269. doi: 10.1002/cphc.200800467

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Du, L. K., Gao, J., Wang, L. L., Song, B., and Liu, C. B. (2015). Application of polarizable ellipsoidal force field model to pnicogen bonds. J. Comput. Chem. 36, 441–448. doi: 10.1002/jcc.23819

PubMed Abstract | CrossRef Full Text | Google Scholar

Loncharich, R. J., Brooks, B. R., and Pastor, R. W. (1992). Langevin dynamics of peptides: the frictional dependence of isomerization rates of N-acetylalanyl-N′-methylamide. Biopolymers 32, 523–535. doi: 10.1002/bip.360320508

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, T., Zhang, J. Q., Chen, J. H., Gou, Q., Xia, Z. N., and Feng, G. (2019). Structure and non-covalent interactions of 1,3-difluoropropane and its complex with water explored by rotational spectroscopy and quantum chemical calculations. J. Chem. Phys. 150:064305. doi: 10.1063/1.5079564

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y., Liu, Y., Xu, Z., Li, H., Liu, H., and Zhu, W. (2012). Halogen bonding for rational drug design and new drug discovery. Expert Opin. Drug Disc. 7, 375–383. doi: 10.1517/17460441.2012.678829

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y., Shi, T., Wang, Y., Yang, H., Yan, X., Luo, X., et al. (2009). Halogen bonding-a novel interaction for rational drug design? J. Med. Chem. 52, 2854–2862. doi: 10.1021/jm9000133

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier, A., Martinez, J., Kasavajala, C., Wickstrom, K., Hauser, K., et al. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11:3696–3713. doi: 10.1021/acs.jctc.5b00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Metrangolo, P., Neukirch, H., Pilati, T., and Resnati, G. (2005). Halogen bonding based recognition processes: a world parallel to hydrogen bonding. ChemInform 36, 386–395. doi: 10.1002/chin.200529294

CrossRef Full Text | Google Scholar

Molina, P., Zapata, F., and Caballero, A. (2017). Anion recognition strategies based on combined noncovalent interactions. Chem. Rev. 117, 9907–9972. doi: 10.1021/acs.chemrev.6b00814

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, R. B., Philipp, D. M., and Friesner, R. A. (2000). A mixed quantum mechanics/molecular mechanics (QM/MM) method for large-scale modeling of chemistry in protein environments. J. Comput. Chem. 21, 1442–1457. doi: 10.1002/1096-987X(200012)21:16<1442::AID-JCC3>3.0.CO;2-O

CrossRef Full Text | Google Scholar

Pang, Y.-P. (1999). Novel zinc protein molecular dynamics simulations: steps toward antiangiogenesis for cancer treatment. Mol. Model. Ann. 5, 196–202. doi: 10.1007/s008940050119

CrossRef Full Text | Google Scholar

Parra, R. D., and Ohlssen, J. (2008). Cooperativity in intramolecular bifurcated hydrogen bonds: an ab initio study. J. Phys. Chem. A 112, 3492–3498. doi: 10.1021/jp711956u

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., et al. (2011). Scikit-learn: machine learning in python. J. Machine Learn. Res. 12, 2825–2830.

Google Scholar

Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., et al. (2005). Scalable molecular dynamics with NAMD. J. Comput. Chem. 26, 1781–1802. doi: 10.1002/jcc.20289

PubMed Abstract | CrossRef Full Text | Google Scholar

Politzer, P., Murray, J. S., and Clark, T. (2013). Halogen bonding and other [sigma]-hole interactions: a perspective. Phys. Chem. Chem. Phys. 15, 11178–11189. doi: 10.1039/c3cp00054k

PubMed Abstract | CrossRef Full Text | Google Scholar

Politzer, P., Murray, J. S., Clark, T., and Resnati, G. (2017). The sigma-hole revisited. Phys. Chem. Chem. Phys. 19, 32166–32178. doi: 10.1039/c7cp06793c

PubMed Abstract | CrossRef Full Text | Google Scholar

Riley, K. E., Murray, J. S., Fanfrlik, J., Rezac, J., Sola, R. J., Concha, M. C., et al. (2013). Halogen bond tunability II: the varying roles of electrostatic and dispersion contributions to attraction in halogen bonds. J. Mol. Model. 19, 4651–4659. doi: 10.1007/s00894-012-1428-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, D., Bertrand, T., Carry, J.-C., Halley, F., Karlsson, A., Mathieu, M., et al. (2016). Differential water thermodynamics determine PI3K-Beta/delta selectivity for solvent-exposed ligand modifications. J. Chem. Inform. Model. 56, 886–894. doi: 10.1021/acs.jcim.5b00641

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryckaert, J.-P., Ciccotti, G., and Berendsen, H. J. C. (1977). Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. doi: 10.1016/0021-9991(77)90098-5

CrossRef Full Text | Google Scholar

Santos, L. A., da Cunha, E. F. F., and Ramalho, T. C. (2017). Toward the classical description of halogen bonds: a quantum based generalized empirical potential for fluorine, chlorine, and bromine. J. Phys. Chem. A 121, 2442–2451. doi: 10.1021/acs.jpca.6b13112

PubMed Abstract | CrossRef Full Text | Google Scholar

Solimannejad, M., Malekani, M., and Alkorta, I. (2013). Substituent effects on the cooperativity of halogen bonding. J. Phys. Chem. A 117, 5551–5557. doi: 10.1021/jp405211p

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinmann, C., Olsson, M. A., and Ryde, U. (2018). Relative ligand-binding free energies calculated from multiple short QM/MM MD simulations. J. Chem. Theory Comput. 14, 3228–3237. doi: 10.1021/acs.jctc.8b00081

PubMed Abstract | CrossRef Full Text | Google Scholar

Udovicic, M., Bazdaric, K., Bilic-Zulle, L., and Petrovecki, M. (2007). What we need to know when calculating the coefficient of correlation? Biochem. Med. 17, 10–15. doi: 10.11613/BM.2017.047

CrossRef Full Text | Google Scholar

Visscher, K. M., Swope, W. C., and Geerke, D. P. (2018). A QM/MM derived polarizable water model for molecular simulation. Molecules 23:3131. doi: 10.3390/molecules23123131

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J. M., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi: 10.1002/jcc.20035

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilcken, R., Liu, X. R., Zimmermann, M. O., Rutherford, T. J., Fersht, A. R., Joerger, A. C., et al. (2012). Halogen-enriched fragment libraries as leads for drug rescue of mutant p53. J. Am. Chem. Soc. 134, 6810–6818. doi: 10.1021/ja301056a

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, W. J., Zeng, Y. L., Li, X. Y., Zhang, X. Y., Zheng, S. J., and Meng, L. P. (2013). Interplay between halogen bonds and hydrogen bonds in OH/SH center dot center dot center dot HOX center dot center dot center dot HY (X = Cl, Br; Y = F, Cl, Br) complexes. J. Mol. Model. 19, 1069–1077. doi: 10.1007/s00894-012-1657-z

CrossRef Full Text | Google Scholar

Yan, X. C., Schyman, P., and Jorgensen, W. L. (2014). Cooperative effects and optimal halogen bonding motifs for self-assembling systems. J. Phys. Chem. A 118, 2820–2826. doi: 10.1021/jp501553j

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Mahmoud, A. H., and Lill, M. A. (2019). Modeling of halogen-protein interactions in co-solvent molecular dynamics simulations. J. Chem. Inform. Model. 59, 38–42. doi: 10.1021/acs.jcim.8b00806

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Q., Feng, D., and Hao, J. (2011). The cooperativity between hydrogen and halogen bond in the XY···HNC···XY (X, Y = F, Cl, Br) complexes. J. Mol. Model. 7, 2817–2823. doi: 10.1007/s00894-011-0974-y

CrossRef Full Text | Google Scholar

Zhou, P. P., Yang, X., Ye, W. C., Zhang, L. W., Yang, F., Zhou, D. G., et al. (2016). Competition and cooperativity of sigma-hole and pi-hole intermolecular interactions between carbon monoxide and bromopentafluorobenzene. New J. Chem. 40, 9139–9147. doi: 10.1039/c6nj01904h

CrossRef Full Text | Google Scholar