Proposed Hydrogen-Bonding Index of Donor or Acceptor Reflecting Its

Article pubs.acs.org/jcim

Proposed Hydrogen-Bonding Index of Donor or Acceptor Reflecting Its Intrinsic Contribution to Hydrogen-Bonding Strength Suqing Zheng,† Shaofang Xu,† Guitao Wang,‡ Qing Tang,† Xiaonan Jiang,‡ Zhanting Li,‡ Yong Xu,§ Renxiao Wang,*,‡ and Fu Lin*,† †

School of Pharmaceutical Sciences, Wenzhou Medical University, Wenzhou 325035, P. R. China State Key Laboratory of Bioorganic Chemistry, Shanghai Institute of Organic Chemistry, Chinese Academy of Sciences, 345 Lingling Road, Shanghai 200032, P. R. China § Institute of Chemical Biology, Guangzhou Institutes of Biomedicine and Health, Chinese Academy of Sciences, 190 Kaiyuan Avenue, Guangzhou Science Park, Guangzhou, Guangdong 510530, P. R. China ‡

S Supporting Information *

ABSTRACT: In this work, we tentatively propose that the hydrogenbonding strength EHB (referring to the minimal hydrogen-bonding energy) and its corresponding hydrogen-bond (HB) distance (referring to the optimal HB distance dHB) for simple mono-HB systems have an exponential relationship on the basis of MP2 and DFT computational results. We take a step further and propose that the hydrogen-bonding indices of the donor (Idonor) and acceptor (Iacceptor), reflecting their intrinsic contributions to hydrogen-bonding strength, also have an exponential relation with the hypothetical effective hydrogen-bond radii of the donor (rdonor) and acceptor (racceptor), respectively. On the basis of extensive quantum-mechanical calculations, relevant assumptions about the hydrogen-bonding index are rationalized. Moreover, the hydrogen-bonding index is also suggested as an additional prefiltering criterion for virtual screening besides the widely accepted Lipinski’s rule of five. Finally, a “Hydrogen-Bond Index Estimator (HBIE)” module has been implemented in our Visual Force Field Derivation Toolkit (VFFDT) program to approximately and rapidly estimate the hydrogen-bonding indices of any small molecules in batch and screen possible stronger donors or acceptors from the small-molecule database. To the best of our knowledge, the concept of the hydrogen-bonding index and its potential application are proposed here for the first time. pKBHX.11 In the past, many theoretically calculated properties have been proposed to predict the hydrogen-bonding free energy by quantitative structure−activity or structure−property relationships,9,18,20,22−31 including the molecular electrostatic potential,27,31,32 local molecular parameters,22 and ISIDA fragment descriptors.28 Hydrogen-bonding energy generally refers to the energy of the interaction between the donor and acceptor for any specific hydrogen-bond geometry from the computational perspective.7,14,15 In most cases, the hydrogen-bonding energy can be practically computed by quantum mechanics (QM) with the supermolecule approach,33 where the hydrogen-bonding energy often refers to the potential energy difference between the hydrogen-bond complex and its corresponding donor and acceptor in vacuum.7 Occasionally, however, the term hydrogen-bonding energy also refers to the enthalpy or Gibbs free energy difference in vacuum or solvent, depending on the userspecified keyword in the QM computation.34

1. INTRODUCTION A hydrogen bond (HB) is an attractive interaction between a hydrogen atom attached to an electronegative atom (donor) and another electronegative atom (acceptor).1−6 Hydrogenbond is an important type of inter- and intramolecular interaction that often plays a crucial role in biological events. Hence the energetics of hydrogen-bond between the donor and acceptor is of great interest to both experimental and computational chemists. Basically three concepts are commonly used to describe the energetics of hydrogen-bond,7 namely, the hydrogen-bonding free energy,8−13 the hydrogen-bonding energy,7,14,15 and the hydrogen-bonding strength.7,14−22 The hydrogen-bonding free energy is a rigorous physical concept intimately connected with the equilibrium constant. It usually refers to Gibbs free energy including entropy and solvation/desolvation effects, which can be experimentally measured by the equilibrium constant for HB complex formation.8−13 This experimental quantity can be represented on the pKHB scale against 4-fluorophenol,12 the log Kβ scale against 4-nitrophenol,8 the log KBH scale against a number of reference HB donors,13 or the hydrogen-bond basicity scale, © 2017 American Chemical Society

Received: January 13, 2017 Published: June 1, 2017 1535

DOI: 10.1021/acs.jcim.7b00022 J. Chem. Inf. Model. 2017, 57, 1535−1547

Article

Journal of Chemical Information and Modeling

In this work, we test whether such a correlation between QM-computed hydrogen-bonding strengths and the corresponding optimal hydrogen-bond distances still holds. More importantly, the ultimate goal of this project is to try to develop a simple and intuitive quantity to approximately separate the intrinsic contributions of the HB donor and acceptor to the hydrogen-bonding strength. Thus, we tentatively adopt the similarly shaped alternative exponential function −exp(k1dHB + k2) to describe such correlation, since this function can be conveniently manipulated with logarithm operation for the purpose of easier separation. On the basis of this exponential function, we propose and rationalize the new concept of the hydrogen-bonding index of the donor or acceptor, which can be derived from the hydrogen-bonding strength. Finally, in order to present the preliminary applicability of this concept, the hydrogen-bonding index is suggested as an additional prefiltering criterion for virtual screening commonly used in the modern drug discovery process, and we have developed a module called the “Hydrogen-Bonding Index Estimator” (HBIE) and included it in our Visual Force Field Derivation Toolkit (VFFDT) software42 to help medicinal chemists screen or select compounds with strong donors/acceptors from a small-molecule database.

Hydrogen-bonding strength frequently means the minimal hydrogen-bonding energy (usually a negative value) between the donor and acceptor when the hydrogen-bond geometry is located at the local/global minimum.7,14−22 However, in some literature the notions of hydrogen-bonding strength and hydrogen-bonding energy seem interchangeable.14,15,26,35 For clarity, the term hydrogen-bonding strength used in this work specifically refers to the minimal hydrogen-bonding energy based on the potential energy scan (PES) method, which is a definition similar to the one used by Hao.7 Thus, in this work, hydrogen-bonding strength and hydrogen-bonding energy are assumed to be two different concepts. In the past two decades, the relationship between the experimental hydrogen-bonding free energy and hydrogenbonding strength has been well-studied.7,8,12,13 Abraham and co-workers experimentally demonstrated that there is a linear relationship between the hydrogen-bonding free energy, as measured by the equilibrium constant for hydrogen-bond complex formation, and the hydrogen-bonding strength derived from the spectroscopic wavelength shift for the hydrogen-bond system.8,13 On the basis of this observation, it is postulated that the entropy term of the hydrogen-bonding free energy either changes linearly with the enthalpy or is trivial relative to the enthalpy.7 The linear relationship has also been addressed in three typical computational works.7,16,36 Hao7 calculated the hydrogen-bonding strength using density functional theory (DFT)37 at the B3LYP/6-31++G** level38,39 corrected for basis set superposition error (BSSE).33 They found that the hydrogenbonding strength is linearly correlated with the logarithm of the experimental hydrogen-bonding equilibrium constant obtained by isothermal titration calorimetry with an R2 value of 0.94, while the logarithm of the hydrogen-bonding equilibrium constant is proportional to the hydrogen-bonding free energy. Moreover, Nocker et al.36 reported that hydrogen-bonding strength computed using the B3LYP functional combined with Dunning’s augmented double-ζ basis set (aug-cc-pVDZ)40 correlates well with the logarithm of the experimental hydrogen-bonding equilibrium constant. In addition, El Kerdawy et al.16 recommended the use of DFT calculations at the X3LYP/aug-cc-pVDZ level41 with one angle constraint to derive the hydrogen-bonding strengths of diverse molecules in complex with a standard donor or acceptor in vacuum, and their work showed a good linear correlation between the hydrogen-bonding strength and the logarithm of the experimental hydrogen-bonding equilibrium constant. In summary, all of these works illustrate that the experimental hydrogenbonding free energy is linearly correlated with the theoretically computed hydrogen-bonding strength. In addition to the linear relationship between the experimental hydrogen-bonding free energy and the theoretically computed hydrogen-bonding strength,7,8,12,16,36 Lamarche and Platts34 found the existence of a relationship between the experimental hydrogen-bonding free energy and the theoretically computed optimal HB distance, dHB . According to their study of about 40 hydrogen-bond systems, the relationship can be described by the formulation k1 + k2/ (dHB + k3), where k1, k2, and k3 are constant parameters, with an R2 value of 0.85.34 Thus, their work leads to our hypothesis that the relationship between the theoretically computed hydrogenbonding strength (EHB) and dHB can also be described by the function k1 + k2/(dHB + k3).

2. METHODS 2.1. Selection of Hydrogen-Bond Donor and Acceptor Model Systems. We selected 10 and 11 small model compounds to represent typical common HB donors and acceptors, respectively. The compact size of the model compounds aims to alleviate undesirable secondary interactions between the HB donor and acceptor. All these structures were first optimized at the B3LYP/6-31++G(d,p) level with Grimme’s D3 correction43,44 in Gaussian 0945 and subsequently used to construct 110 HB complexes, which considers all of the possible donor−acceptor pairs between the model molecules in Table 1. For each HB complex, only a distance-dependence scan of the hydrogen-bonding energy was conducted to explore the minimal hydrogen-bonding energy. The two hydrogen-bond angles (the X−H···A angle, θ, and the H···A−Y angle, ϕ, in Figure 1) were fixed at 180° and were not scanned systematically since a simultaneous scan along the hydrogenbond distance and two hydrogen-bond angles would have consumed large amounts of computational resources and such scans often encounter failure of convergence during QM optimization. Hence the constraint of the hydrogen-bond angles to 180° is usually adopted in theoretical calculations on hydrogen bonding7,16 to maximally reduce the secondary interactions between the donor and acceptor. For some special cases in which the acceptor atom is attached to two neighboring atoms (e.g., the O atom in water (wat), the O atom in dimethyl ether (met), the O atom in furan (frn), the N atom in pyridine (py), etc.), a dummy atom was created midway between the two neighboring atoms. This virtual dummy atom was treated as atom Y (Figure 1) to define the hydrogen-bond angle ϕ = 180° (Figure 1) so as to keep the symmetry during the geometry optimization. Additionally, Hao7 concluded that the hydrogen-bonding energy is not sensitive to the directionality of the hydrogen bond. Hence, it was reasonable to perform the distance-dependence scan of the hydrogen-bonding energy with the two hydrogen-bond angles constrained. More specifically, the hydrogen-bond distance (d) in Figure 1 was changed from 1.6 to 4.5 Å with an interval of 0.1 Å, 1536


Article


energy was computed for each configuration with the supermolecule approach33 as described by eq 1:

Table 1. Model Molecules for Hydrogen-Bond Donor and Acceptor

E = E DA − E D − EA + E BSSE

(1)

where E is the hydrogen-bonding energy between the donor (D) and acceptor (A), EDA is the total energy of the donor− acceptor (DA) complex, ED and EA are the energies of D and A, respectively, with the configurations of D and A taken directly from the structure of the HB complex, and EBSSE is the optional BSSE correction energy33 depending on the different computational methods as follows. In order to yield more robust results, two computational levels (MP247 and DFT37), three functionals (M06-2X,48 B3LYP,38,39 and X3LYP41), and two basis sets (6-31++G(d,p) and aug-cc-pVDZ40) along with Grimme’s D3 correction43,44 and BSSE correction33 were considered in the single-point energy calculations. Hence the computational methods we tried in this work include MP2/6-31++G(d,p)-BSSE, X3LYP/augcc-pVDZ, B3LYP/aug-cc-pVDZ, B3LYP/aug-cc-pVDZ-D3, M06-2X/aug-cc-pVDZ, M06-2X/aug-cc-pVDZ-D3, B3LYP/631++G(d,p)-BSSE, B3LYP/6-31++G(d,p)-BSSE-D3, M06-2X/ 6-31++G(d,p)-BSSE, and M06-2X/6-31++G(d,p)-BSSE-D3. In this step, the BSSE correction was applied to 6-31++G(d,p) basis set but not the aug-cc-pVDZ basis set since the latter is more complete. Grimme’s D3 correction was not applied in our DFT computations with X3LYP functional because X3LYP is intentionally parametrized to significantly improve the accuracy for van der Waals and hydrogen-bond complexes44 and has no compatible option to support Grimme’s D3 correction in Gaussian 09. It should be noted that all of the QM calculations in this work were carried out using the Gaussian 09 software package.45 Upon completion of the single-point energy calculations, the hydrogen-bonding energies were computed for each configuration using eq 1. The hydrogen-bonding energies (E) as functions of the hydrogen-bond distance (d) are plotted in Figures S1−S110 in the Supporting Information. Values of the hydrogen-bonding strength (EHB), which refers to the lowest hydrogen-bonding energy, and the corresponding hydrogenbond distance, termed the optimal hydrogen-bond distance (dHB), were obtained from the curves in Figures S1−S110. The obtained values for the 110 HB complexes are summarized in Table 2 for the MP2 method and in Tables S1−S9 for various DFT methods. To evaluate the effect of different computational approaches on EHB and dHB, relative standard deviations (RSDs) were calculated and are summarized in Table 3. The RSDs were obtained using the expression RSD = |100% × SD/average value|, where SD is the standard deviation and the absolute value sign is necessary because of the negative value of the hydrogen-bonding strength. 2.3. Estimation of the Hydrogen-Bonding Index (I) from the QM-Computed Hydrogen-Bonding Strength. This study hypothesizes that the relationship between the hydrogen-bonding strength EHB and its optimal hydrogen-bond distance dHB can be formulated by the exponential function shown in eq 2:

Figure 1. Geometrical parameters for describing a hydrogen bond: d is the hydrogen-bond distance; and θ and ϕ are the hydrogen-bond angles; X is the donor atom; H is the donor hydrogen directly involved in the hydrogen bond; A is the acceptor atom; Y is the neighboring atom of A or a dummy atom defined to keep the symmetry during the geometry optimization if there is more than one neighboring atom.

resulting in 30 configurations. In each configuration, the two hydrogen-bond angles θ and ϕ in Figure 1 were fixed. Structural optimization of each configuration was conducted at the B3LYP/6-31++G(d,p)-D3 level with the geometric constraints of d, θ, and ϕ, which generated a total of 3300 (or 110 × 30) optimized configurations. In this optimization stage, this computationally affordable yet relatively accurate QM calculation level was recommended.7,46 However, different higher levels of QM computation were utilized in the subsequent single-point energy calculations. 2.2. Computation of the Hydrogen-Bonding Strength and Its Corresponding Optimal Hydrogen-Bond Distance. After structural optimization, the hydrogen-bonding

E HB = −E θ ekdHB+ b

(2) θ

where k and b are constants, the quantity E = 1 kcal/mol has been included to give EHB the desired energy units (since the exponential function is dimensionless), and the negative sign ensures that EHB is negative (since the exponential factor is 1537


Article

dHB

2.3 2.1 2.4 2.3 2.2 2.3 2.4 2.2 2.9 3.0 3.1 −3.36 −4.30 −2.67 −3.54 −4.40 −2.71 −3.41 −5.15 −2.26 −1.50 −1.25 2.3 2.2 2.4 2.3 2.3 2.4 2.5 2.3 2.9 3.1 3.1 −2.16 −2.77 −1.64 −2.16 −2.71 −1.64 −2.12 −3.30 −1.25 −0.77 −0.68

ln( −E HB /E θ ) = kdHB + b

(3)

For ease of description, the concept of “hydrogen-bond magnitude scale” (Zdonor−acceptor) is proposed and defined as the natural logarithm of the additive inverse of the hydrogenbonding strength (eq 4): Zdonor−acceptor ≡ ln( −E HB /E θ )

2.2 2.1 2.2 2.2 2.1 2.2 2.3 2.1 2.8 2.9 2.9

dHB ett EHB dHB

positive). It should be noted that Eθ will be omitted in eqs 10−12 below for the sake of conciseness. The function shown in eq 2 can be transformed into eq 3 by applying the operation of natural logarithm:

(4)

2.4 2.2 2.4 2.3 2.3 2.4 2.5 2.3 3.0 3.1 3.2 −2.10 −2.98 −1.88 −2.43 −3.19 −1.88 −2.42 −3.42 −1.20 −0.89 −0.66

−3.91 −4.82 −3.05 −3.94 −4.79 −3.04 −3.94 −5.63 −2.29 −1.31 −1.17

dHB

2.3 2.2 2.3 2.3 2.2 2.3 2.4 2.3 2.9 3.1 3.1 −2.66 −3.11 −2.13 −2.77 −3.54 −1.74 −2.83 −3.67 −1.27 −0.78 −0.53

EHB

−2.32 −2.92 −1.78 −2.32 −2.72 −1.82 −2.26 −3.32 −1.30 −0.77 −0.71

dHB

2.0 1.8 2.1 2.0 1.9 2.0 2.1 1.8 2.6 2.9 2.9

eta

Zdonor−acceptor = kdHB + b

(5)

The optimal hydrogen-bond distance (dHB) in eq 5 is assumed to be the sum of the hypothetical effective hydrogenbond radii of the donor (rdonor) and acceptor (racceptor) in the hydrogen-bond complex X−H···A (eq 6). dHB = rdonor + racceptor

(6)

The magnitude scales for the donor (Zdonor) and acceptor (Zacceptor) are also defined as the natural logarithms of the additive inverses of the hydrogen-bonding indices Idonor and Iacceptor, respectively (eqs 7 and 8):

2.3 2.2 2.4 2.3 2.3 2.3 2.4 2.3 2.9 3.1 3.1

dHB mea EHB

pri EHB

EHB

aln

From eqs 2−4, eq 5 can be simply obtained as follows:

dHB

Zdonor ≡ ln( −Idonor) = k1rdonor + b1

(7)

Zacceptor ≡ ln( −Iacceptor) = k 2racceptor + b2

(8)

Similar to eq 5, Zdonor and Zacceptor are also hypothesized to depend linearly on rdonor and racceptor, respectively (eqs 7 and 8). Moreover, Zdonor−acceptor is assumed to consist of Zdonor and Zacceptor (eq 9):

−5.67 −7.64 −4.68 −6.07 −7.63 −4.64 −6.39 −9.62 −3.21 −1.69 −1.21 2.0 1.9 2.1 2.1 2.0 2.1 2.1 2.0 2.7 3.0 2.9 −4.44 −5.35 −3.24 −4.18 −5.24 −3.13 −4.43 −6.56 −2.07 −0.91 −0.87 2.0 1.9 2.1 2.1 2.0 2.1 2.1 2.0 2.7 2.9 2.9 2.1 2.0 2.2 2.1 2.0 2.2 2.2 2.0 2.8 3.1 3.0 −4.26 −4.83 −3.27 −4.16 −5.20 −2.65 −4.40 −5.91 −1.69 −0.84 −0.56

(9)

In the preceding equations, k1, k2, b1, and b2 are constants and Zdonor−acceptor, Zdonor, Zacceptor, Idonor, and Iacceptor are dimensionless quantities. In light of eq 9, the hydrogen-bond magnitude scale, Zdonor−acceptor, is additive and is simply the sum of the magnitude scales for the donor (Zdonor) and acceptor (Zacceptor). Thus, Zdonor and Zacceptor can be derived through standard multilinear regression (MLR) if Zdonor−acceptor for a number of donor− acceptor pairs is available. In this work, 110 hydrogen-bond complexes formed from the different 21 model compounds in Table 1 were employed in the MLR, while the Zdonor−acceptor values used in the MLR were calculated from the hydrogenbonding strengths (Tables 2 and S1−S9) using eq 4. From the MLR method, the magnitude scales for the donor (Zdonor) and acceptor (Zacceptor) were obtained and are listed in Table S10. From them, the corresponding hydrogen-bonding indices (Idonor and Iacceptor) were calculated using their definitions (eqs 7 and 8, respectively) and are provided in Table S11. 2.4. Estimation of the Hypothetical Effective Hydrogen-Bond Radii (rdonor and racceptor) from QM-Computed Optimal Hydrogen-Bond Distances. To demonstrate the validity of eqs 7 and 8, the hypothetical effective hydrogenbond radii of the donor (rdonor) and acceptor (racceptor) are required. Eq 6 suggests that the optimal hydrogen-bond distance (dHB) is also additive and simply the sum of rdonor and racceptor. Thus, the values of rdonor and racceptor were derived

−4.51 −5.38 −3.28 −4.21 −5.29 −3.14 −4.48 −6.58 −2.08 −0.90 −0.85

dHB

acceptors

wat EHB

EHB

meo

dHB

EHB

eto

dHB

EHB

phy

Zdonor−acceptor = Zdonor + Zacceptor

wat met fad act ace frn cn py dms prt tph

donors

Table 2. Hydrogen-Bonding Strengths (EHB, in kcal/mol) and the Corresponding Optimal Hydrogen-Bond Distances (dHB, in Å) Computed by the MP2 Method

bnt EHB


1538


Article


Table 3. Relative Standard Deviations (RSDs, in %) of EHB and dHB among 10 Different QM (MP2 and DFT) Computational Methods donors wat

acceptors

wat met fad act ace frn cn py dms prt tph

meo

eto

phy

eta

mea

pri

aln

ett

bnt

EHB

dHB

EHB

dHB

EHB

dHB

EHB

dHB

EHB

dHB

EHB

dHB

EHB

dHB

EHB

dHB

EHB

dHB

EHB

dHB

8 8 6 6 5 11 4 7 11 18 28

1 1 3 1 1 2 2 2 2 3 2

7 10 8 8 7 14 16 9 15 25 31

2 1 2 2 2 2 2 2 2 3 2

8 10 8 8 7 15 6 9 16 27 33

2 1 2 2 2 2 2 2 2 3 2

9 12 10 9 9 16 8 10 18 26 44

2 2 3 1 2 2 2 3 1 4 4

14 20 17 16 11 27 13 17 30 42 50

2 2 1 2 2 2 0 2 3 4 4

11 16 11 10 9 23 8 13 23 30 48

1 2 2 1 1 2 1 2 2 3 3

16 18 15 19 16 24 16 16 42 32 49

2 3 2 2 3 2 2 2 5 4 4

10 15 12 12 10 21 9 13 22 31 39

1 1 2 2 2 2 1 2 2 3 3

17 20 17 15 13 29 13 18 30 40 51

2 2 2 1 2 3 1 2 2 5 4

15 20 16 16 14 28 13 18 27 34 45

1 3 3 2 2 3 2 3 3 4 4

from the QM-computed dHB values (Tables 2 and S1−S9) by the MLR method and are listed in Table S12.

evidently manifested. For further clarity, the hydrogen-bonding strength and hydrogen-bonding energy are intentionally denoted by the symbols EHB and E, respectively, and the optimal hydrogen-bond distance and hydrogen-bond distance are purposely symbolized by dHB and d, respectively. In this work, the hydrogen-bonding strength and optimal hydrogenbond distance, rather than the hydrogen-bonding energy and hydrogen-bond distance, are probed because the relation between EHB and dHB is not well-addressed, whereas the relation between E and d is quite straightforward from Figures 2 and S1−S110. Before the relation between EHB and dHB was explored, the relative standard deviation (RSD) was used to evaluate the influence of MP2 and different DFT computational methods on the values of EHB and dHB. On the basis of the RSDs in Table 3, the differences among the dHB values derived with various computational approaches are unanimously small, which can be observed from the low RSD values (0−5%). Especially for the eta···cn hydrogen-bond system, the RSD is 0%, which means that all of the computational methods gave the same optimal hydrogen-bond distance (2.4 Å). However, the hydrogen-bonding strengths exhibit large discrepancies among different QM methods, since the RSD values range from 4% to 51%. Therefore, the relatively accurate MP2 computational results will be the focus of the discussion, while various DFT computational results will be also employed to test the robustness of the proposed exponential relationship between hydrogen-bonding strength and optimal hydrogen-bond distance as described below. 3.2. Relationship between the Hydrogen-Bonding Strength and Its Corresponding Optimal HydrogenBond Distance. The hydrogen-bonding strength has been theoretically demonstrated to be linearly correlated with the experimental hydrogen-bonding free energy,7,8,12,16,36 whereas the hydrogen-bonding free energy itself was shown to be related to the computed optimal hydrogen-bond distance (dHB) by the formulation k1 + k2/(dHB + k3) with R2 = 0.85 in the work of Lamarche and Platts.34 Thus, we hypothesize that a similar relationship should also hold between the theoretical hydrogen-bonding strength (EHB) and the computed optimal hydrogen-bond distance (dHB). The scatter plot of EHB and its corresponding dHB in Figure 3 demonstrates that our hypothesized relationship indeed holds and can be welldescribed (R2 = 0.85) by eq 10:

3. RESULTS AND DISCUSSION 3.1. Derivation of Hydrogen-Bonding Strengths from the Distance-Dependent Hydrogen-Bonding Energies. The graphical illustrations of all of the hydrogen-bonding energies as functions of the different hydrogen-bond distances, with 1100 curves in total, are given in Figures S1−S110. According to Figures S1−S110, the profiles of the hydrogenbonding energy and hydrogen-bond distance are Morse-like curves (one typical exponential curve), which have a very similar shape as the curves obtained by Hao.7 Also, these profiles show that there is an energy minimum in each curve, which corresponds to the hydrogen-bonding strength and optimal hydrogen-bond distance (Tables 2 and S1−S9). In order to make a distinction between the two pairs of concepts hydrogen-bond distance/optimal hydrogen-bond distance and hydrogen-bonding energy/hydrogen-bonding strength, aln···act, one of hydrogen-bond systems, is taken as a typical example. In Figure 2, the difference between the hydrogen-bond distance and optimal hydrogen-bond distance is vividly shown, and the distinction between the hydrogenbonding strength and hydrogen-bonding energy is also

Figure 2. Profile of the hydrogen-bonding energy and hydrogen-bond distance for the hydrogen-bond complex aln···act. 1539


Article


To better reproduce the constants k1 and k2 from Wendler et al., eq 12 was used to fit our MP2 computational results (green curve in Figure 3): E HB = k1/dHB3.78

(12)

The fit offered a k1 value of −71.05 and an R value of 0.85. Obviously, this k1 value (−71.05) is close to the value of −78.99 in the work of Wendler et al., while the R2 value (0.85) is also same as its counterpart in the work of Lamarche and Platts. Thus, the work of Wendler et al. is also in accord with our hypothesis. In short, eqs 10−12 can well describe the relationship between hydrogen-bonding strength and optimal hydrogenbond distance with comparable accuracies. Nevertheless, these equations cannot be simply manipulated to separate the contributions of the donor and acceptor to the hydrogenbonding strength. Therefore, we explored an alternative equation to meet the following two criteria: (1) the equation should share a similar shape with eqs 10−12, and (2) the equation should be easily transformed to separate the intrinsic capabilities of the donor and acceptor, which can be defined by a simple and intuitive quantity. Fortunately, the scatter plot (Figure 3) also indicates the exponential relation (eq 2) between hydrogen-bonding strength (EHB) and its corresponding optimal hydrogen-bond distance (dHB) with the same R2 value of 0.85. By the use of this exponential function (eq 2), two parameters (k1 = −1.91 and k2 = 5.48) are sufficient to achieve the same R2 value (0.85) as for eq 10 based on our MP2 computational results (red curve in Figure 3). Furthermore, the natural logarithm of the additive inverse of the hydrogen-bonding strength, ln(−EHB/Eθ), denoted by the symbol Zdonor−acceptor, as a function of the optimal hydrogen-bond distance dHB is plotted in Figure 4, 2

Figure 3. Scatter plot of hydrogen-bonding strength (EHB) and optimal hydrogen-bond distance (dHB) from our MP2 computation. The red, black, pink, and green curves are fits to eqs 2, 10, 11, and 12 respectively. The orange curve is plotted by eq 11 with the two parameters k1 (−78.99) and k2 (3.78) originating from Wendler et al.49 The square of the correlation coefficient (R2 value) is 0.85 for all five curves.

E HB = k1 + k 2/(dHB + k 3)

(10)

where k1, k2, and k3 are constants, which have the values 1.34, −3.77, and −1.42, respectively, based on our MP2 computational results (black curve in Figure 3). However, the specific values of k1, k2, and k3 were not reported by Lamarche and Platts,34 and thus, it is impossible for us to make a direct comparison. Nevertheless, the R2 value reported by Lamarche and Platts (0.85) is the same as ours, which can partially verify our hypothesis from this perspective. In addition, Wendler et al.49 performed DFT optimization at the BLYP-D/def2-TZVPP level without any hydrogen-bond angle constraints for 256 hydrogen-bond dimers, including peptides and water. Interestingly, they proposed eq 11 to describe the relation between hydrogen-bonding strength and optimal hydrogen-bond distance: E HB = k1/dHBk 2

(11)

where k1 and k2 are constants, which have the values −78.99 and 3.78, respectively, based on the work of Wendler et al. (orange curve in Figure 3). It should be mentioned that k1 = −1.2 × 1010 and k2 = 3.78 in the original work,49 since in their work kilojoules per mole and picometers were adopted as the units of EHB and dHB, respectively, rather than kilocalories per mole and angstroms as used in our work. To verify whether eq 11 also works in our model systems, it was used to fit the hydrogen-bonding strength and optimal hydrogen-bond distance from our MP2 computation (pink curve in Figure 3), which gave the constants k1 and k2 with values of −96.69 and 4.19, respectively, and an R2 value of 0.85. The value of k2 (4.19) is close to the value of 3.78 from the work of Wendler et al., while the value of k1 (−96.69 vs −78.99) seemingly differs much. However, as can be seen from Figure 3, our fitting curve (pink curve) is also close to the curves from Wendler et al. (orange curve) and Lamarche and Platts (black curve), and the R2 value (0.85) is also the same as its counterpart in the work of Lamarche and Platts.

Figure 4. Correlation between magnitude scale of hydrogen-bond (Zdonor−acceptor) and optimal hydrogen-bond distance (dHB). The square of correlation coefficient (R2-value) is 0.86.

which demonstrates that there is a robust linear relationship between Zdonor−acceptor and dHB with R2 = 0.86. Thus, it makes sense to assume that the hydrogen-bonding strength (EHB) is exponentially related to its corresponding optimal hydrogenbond distance (dHB) on the basis of the MP2 computational results. More interestingly, Zdonor−acceptor is additive and can be decomposed into the contributions of the donor (Zdonor) and acceptor (Zacceptor) on the basis of eq 9. Hence, our proposed exponential relation satisfies the two aforementioned criteria. 1540


Article


Table 4. Correlations (R2 Values) between Hydrogen-Bond Magnitude Scale (Zdonor−acceptor) and Its Corresponding Optimal Hydrogen-Bond Distance (dHB)a method

QM ZQM donor−acceptor vs dHB

MLR ZMLR donor−acceptor vs dHB

MLR ZQM donor−acceptor vs dHB

MLR Zdonor−acceptor vs dQM HB

001_MP2_6-31++G(d p)-BSSE 002_X3LYP_aug-cc-pVDZ 003_B3LYP_aug-cc-pVDZ 004_B3LYP_aug-cc-pVDZ-D3 005_M062X_aug-cc-pVDZ 006_M062X_aug-cc-pVDZ-D3 007_B3LYP_6-31++G(d,p)-BSSE 008_B3LYP_6-31++G(d,p)-BSSE-D3 009_M062X_6-31++G(d,p)-BSSE 010_M062X_6-31++G(d,p)-BSSE-D3

0.862 0.880 0.893 0.821 0.832 0.812 0.886 0.849 0.872 0.866

0.876 0.890 0.903 0.831 0.852 0.829 0.899 0.858 0.886 0.883

0.851 0.881 0.893 0.800 0.831 0.804 0.854 0.835 0.872 0.865

0.864 0.879 0.891 0.811 0.832 0.812 0.875 0.843 0.868 0.865

a QM Definitions: ZQM donor−acceptor, hydrogen-bond magnitude scale given by QM computation; dHB , optimal hydrogen-bond distance given by QM MLR , hydrogen-bond magnitude scale predicted by the MLR method; d computation; ZMLR donor−acceptor HB , optimal hydrogen-bond distance predicted by the MLR method.

not very promising, with an R2 value of 0.574 (purple curve in Figure 5), which may be ascribed to our relatively diverse

To further test the reliability of this relationship, different computational methods including MP247 and DFT37 were assessed for the calculation of hydrogen-bonding strength (EHB), hydrogen-bond magnitude scale (Zdonor−acceptor), and optimal hydrogen-bond distance (dHB). All of the methods achieved consistent results with R2 values ranging from 0.80 to 0.90 (Table 4). In Table 4, there are four sets of correlations QM QM between Z don or− acceptor and d HB : Z d ono r− acceptor /d HB , MLR MLR QM MLR MLR Zdonor−acceptor/dHB , Zdonor−acceptor/dHB , and Zdonor−acceptor/dQM HB . MLR The R2 values between ZMLR donor−acceptor and dHB are consistently larger than the others. This is probably due to the removal of systematic error from the QM computation by the MLR method. Consequently, the exponential relation between hydrogen-bonding strength and optimal hydrogen-bond distance holds irrespective of the computational method. Furthermore, in the work of Grabowski,26 both hydrogenbonding strength and optimal hydrogen-bond distance were computed at the MP2/6-311++G** level for nine simple donor···acceptor pairs, including HF···H2O, HF···NH3, HF··· LiH, (H2O)2, (HCOOH)2, C2H2···OH2, HOH···C2H2, HOH··· NH3, and (C2H2)2. Grabowski26 indicated that there is no nice linear relationship between the hydrogen-bond strength and optimal hydrogen-bond distance. However, he found that a quadratic function with the modified distance parameter ΔrH···A as independent variable (eq 13) can excellently describe the relation between hydrogen-bonding strength and ΔrH···A with an R2 value of 0.945 for several hydrogen-bond systems: −E HB = k1(ΔrH···A + k 2)2 + k 3

Figure 5. Scatter plot of the additive inverse of the hydrogen-bonding strength (−EHB) vs the modified distance parameter (ΔrH···A). Black spheres refer to hydrogen-bond (HB) systems without sulfur as the acceptor atom, and red spheres stand for HB systems with sulfur as the acceptor atom. Both sets of spheres are from our MP2 computational results in this work. The black quadratic curve was fitted to the black spheres, the red quadratic curve to the red spheres, and the purple quadratic curve to all of the data including both the black and red spheres. The R2 values for the black, red, and purple curves are 0.814, 0.684, and 0.574 respectively. The quadratic curve is described by the equation −EHB = k1 (ΔrH···A + k2)2 + k3. The respective values of k1, k2, and k3 are 2.16, 0.65, and 0.75 for the purple curve; 6.97, 0.18, and 0.83 for the black curve; and 4.07, −0.01, and 1.15 for the red curve.

(13)

where k1, k2, and k3 are constants and ΔrH···A is the sum of the van der Waals radii of the hydrogen atom (1.2 Å) and the acceptor atom A (O, 1.4 Å; N, 1.5 Å; S, 1.8 Å) minus the optimal hydrogen-bond distance (dHB). Additionally, Zou et al.50 adopted the MP2/aug-cc-PVTZ method to calculate hydrogen-bonding strengths and optimal hydrogen-bond distances for 34 hydrogen-bond systems with NH3 as the acceptor (X−H···NH3). They also display that the quadratic form shown in eq 13 can also fit nicely with an R2 value of 0.943. Nevertheless, unlike the exponential function, this quadratic function cannot be transformed handily to separate the intrinsic contributions of donor and acceptor to the hydrogen-bonding strength. Furthermore, a quadratic fit of eq 13 to our MP2 computational results in this work was also attempted, but the overall performance of the quadratic fit for all of the hydrogen-bond systems considered in our work was

model compounds, including ones with sulfur as the acceptor atom. Figure 5 indicates that the R2 value for hydrogen-bond complexes with sulfur as the acceptor atom (red spheres in Figure 5) is 0.684, while the R2 value for hydrogen-bond complexes without sulfur as the acceptor atom (black spheres in Figure 5) is 0.814. Therefore, the quadratic function may not work very well for various hydrogen-bond systems used in this work. However, in another work by Grabowski,51 the hydrogenbonding strength was reported to depend linearly on the optimal hydrogen-bond distance for the hydrogen-bond complexes R−CN···H−X (R = H, F, Li, CH3, CH2F, CHF2, 1541


Article


was demonstrated to hold an exponential relation with optimal hydrogen-bond distance. Similarly, we also suggest that the hydrogen-bonding index (I) is exponentially related to the hypothetical effective hydrogen-bond radius (r) as described by eqs 7 and 8. The hypothetical effective hydrogen-bond radii of the donor (rdonor) and acceptor (racceptor), derived by MLR method on the basis of eq 6, are tabulated in Table S12. The optimal hydrogen-bond distance for a donor−acceptor complex can be simply predicted from the sum of rdonor and racceptor (Table S12), which is linearly related to the QM-computed optimal hydrogen-bond distance dHB (Tables 2 and S1−S9) with R2 values ranging from 0.97 to 0.99 (Table 5). Hence eq 6 is reasonable from this aspect.

CF3, NO2, NH2, BH2, OH; X = F, Cl). It seems that the relationship between hydrogen-bonding strength and optimal hydrogen-bond distance varies among different hydrogen-bond systems. Therefore, the data from all of the studies (Grabowski26,51 and Zou et al.50) were combined with our MP2 computational results and plotted in Figure 6 so as to probe whether their data support our proposed exponential relationship between EHB and dHB.

Table 5. Correlations (R2 Values) between QM-Computed and MLR-Predicted Results (Zdonor−acceptor and dHB)a method

ZMLR donor−acceptor vs ZQM donor−acceptor

dMLR HB vs dQM HB


0.971 0.989 0.990 0.963 0.976 0.979 0.950 0.973 0.984 0.980

0.986 0.987 0.987 0.976 0.976 0.970 0.973 0.983 0.979 0.979

Figure 6. Scatter plot of hydrogen-bonding strength (EHB) and optimal hydrogen-bond distance (dHB). Blue, green, and orange spheres are from our MP2 computational results, Grabowski,26,51 and Zou et al.,50 respectively. The blue spheres were used in the fitting to generate the red curve, while the blue, green, and orange spheres were combined together to obtain the black curve. Both the red and black curves were fitted using eq 2. The R2 values are 0.85 and 0.89 for the red and black curves, respectively.

a Definitions: ZMLR donor−acceptor, hydrogen-bond magnitude scale predicted by the MLR method; ZQM donor−acceptor, hydrogen-bond magnitude scale given by QM computation; dMLR HB , optimal hydrogen-bond distance predicted by the MLR method; dQM HB , optimal hydrogen-bond distance given by QM computation.

As illustrated in Figure 6, the data of Grabowski26,51 (green spheres) and the data of Zou et al.50 (orange spheres) are located at the left half of our data (blue spheres). When all of the data sets (green, orange, and blue spheres) were combined for curve fitting, the R2 value increased from 0.85 to 0.89. The newly fitted curve (black curve in Figure 6) is still highly similar to the original fitted curve (red curve) generated by our data (blue spheres). In this sense, all of the results from Grabowski and Zou et al. corroborate our hypothesis that the hydrogenbonding strength has an exponential relationship with the optimal hydrogen-bond distance (eq 2). Therefore, the works of Lamarche and Platts, Wendler et al., Grabowski, and Zou et al. all agree with our proposed exponential relationship, despite the different equations used in those works. It is worth mentioning that in their works no hydrogen-bond angles were constrained during QM optimization, while in our computational protocol the hydrogen-bond angles were constrained to minimize the undesirable secondary interactions. Nevertheless, from the discussion in this section, it seems that the absence of hydrogen-bond angle constraints does not obviously affect our proposed exponential relationship. This may be attributed to the small model systems used in all of these works, which probably do not exhibit severe secondary interactions. 3.3. Relationship between the Hydrogen-Bonding Index and the Hypothetical Effective Hydrogen-Bond Radius. In the previous section, hydrogen-bonding strength

At the same time, both Zdonor and Zacceptor (Table S10) can be obtained from QM-computed hydrogen-bonding strengths (Tables 2 and S1−S9) via the MLR approach based on eq 9. The MLR-predicted values ZMLR donor−acceptor obtained from Zdonor and Zacceptor are nicely linear with the QM-computed values 2 ZQM donor−acceptor with R values ranging from 0.95 to 0.99 (Table 5). Similarly, the MLR-predicted hydrogen-bonding strength (EHB) is also well-correlated to the QM-computed EHB with R2 values varying from 0.98 to 0.99 (Table S13). For example, Figure 7 clearly depicts an excellent linear relationship between the MLR-predicted and MP2-computed values of EHB with an R2 value of 0.98. Thus, the good R2 value supports the additivity property of Z described by eq 9. With the data in Tables S10−S12, the relationships between the hydrogen-bonding indices I (or donor/acceptor magnitude scales Z) and hypothetical effective hydrogen-bond radii r are ready to be explored. According to Table 6, the HB donor magnitude scale Zdonor is linearly correlated to its corresponding hypothetical effective hydrogen-bond radius rdonor with R2 values ranging from 0.80 to 0.93, while the HB acceptor magnitude scale Zacceptor is linearly associated with racceptor with R2 values varying from 0.88 to 0.92. In other words, Idonor is exponentially related to rdonor, whereas Iacceptor is exponentially related to racceptor. From this perspective, the hypothesized eqs 7 and 8 are rational. Thus, the exponential relation also holds between the hydrogen-bonding index (Idonor or Iacceptor) and hypothetical effective hydrogen-bond radius (rdonor or racceptor). 1542


Article


Actually this intuitive meaning of the hydrogen-bonding index is one of the major reasons to adopt eq 2. In addition, relative to the functions k1 + k2/(dHB + k3) and k1/dHBk2, the exponential function is easier for the logarithm transformation, since we aim to separate the individual contributions of the donor and acceptor to the hydrogen-bonding strength. More importantly, the hydrogen-bonding index can be potentially used as a prefiltering criterion for virtual screening in modern drug discovery projects, which will be discussed below. Built upon these considerations, the exponential function has been used in this work to describe the relation between hydrogenbonding strength and optimal hydrogen-bond distance. 3.5. Interpretation of the Proposed Concepts and Assumptions. For easier description, several concepts have been proposed, including Zdonor−acceptor, Zdonor, Zacceptor, Idonor, Iacceptor, rdonor, and racceptor. Among them, Zdonor−acceptor, Zdonor, and Zacceptor are the magnitude scales for the hydrogen-bond complex, donor, and acceptor, respectively. Idonor and Iacceptor are the hydrogen-bonding indices of the donor and acceptor and can be interpreted as the hypothetical effective charges of the donor and acceptor respectively, which reasonably explains the multiplicativity of Idonor and Iacceptor. The quantities rdonor and racceptor are thought to be the hypothetical effective hydrogenbond radii of the donor and acceptor respectively, while the sum of rdonor and racceptor equals the optimal hydrogen-bond distance (dHB). Furthermore, two assumptions have been rationalized. First, the hydrogen-bonding strength is correlated with the optimal hydrogen-bond distance by the exponential function. Second, the hydrogen-bonding index I is also exponentially related to its corresponding hypothetical effective hydrogen-bond radius r. In other words, the magnitude scale for the donor, Zdonor (or the acceptor, Zacceptor) is linearly correlated with its corresponding hypothetical effective hydrogen-bond radius, rdonor (or acceptor, racceptor). Both assumptions have been verified in this work on the basis of the hydrogen-bonding strengths and optimal hydrogen-bond distances computed using various QM methods. 3.6. Qualitative Application of Hydrogen-Bonding Index. In this study, we have derived the hydrogen-bonding indices of the donor and acceptor for measuring their intrinsic contributions to the hydrogen-bonding strength. The value of the hydrogen-bonding index (I) is negative and reflects the maximum contribution of the donor or acceptor to the hydrogen-bonding strength. As the hydrogen-bonding index of the donor (Idonor) or acceptor (Iacceptor) decreases, the contribution of the donor or acceptor to the hydrogen-bonding strength increases. Figure 8 illustrates that the order of the donor hydrogen-bonding indices is Idonor(phy) < Idonor(meo) ≈ Idonor(eto) ≈ Idonor(aln) < Idonor(bnt) < Idonor(wat) < Idonor(mea) ≈ Idonor(pri) ≈ Idonor(eta) ≈ Idonor(ett), while the rank of the acceptor hydrogen-bonding indices is Iacceptor(py) < Iacceptor(ace) ≈ Iacceptor(met) < Iacceptor(cn) ≈ Iacceptor(act) ≈ Iacceptor(wat) < Iacceptor(fad) ≈ Iacceptor(frn) < Iacceptor(dms) < Iacceptor(prt) ≈ Iacceptor(tph). On the basis of these indices, sulfur-containing model compounds (dms, prt, and tph) are consistently and expectedly poor hydrogen-bond acceptors relative to oxygen/nitrogen-containing compounds, while sulfur-containing donors (ett and bnt) have donor capacities comparable to those of nitrogen-containing compounds, which seems different from our common sense that the sulfur atom does not possess strong electronegativity. More specific trends

Figure 7. Correlation between multilinear regression (MLR)predicted and MP2-computed hydrogen-bonding strengths (EHB). The R2 value is 0.98.

Table 6. Correlations (R2 Values) between Hydrogen-Bond Magnitude Scales (Z) and Hypothetical Effective Radii (r) Zdonor vs

Zacceptor vs

method

rdonor

racceptor


0.822 0.910 0.926 0.821 0.842 0.803 0.913 0.885 0.878 0.864

0.908 0.906 0.918 0.875 0.886 0.875 0.910 0.894 0.901 0.903

3.4. Interpretation of the Hydrogen-Bonding Index. To interpret the concept of hydrogen-bonding indices for the donor (Idonor) and acceptor (Iacceptor), we eliminate the logarithm operation for eq 9, which leads to eq 14: −E HB = ( −Idonor) × ( −Iacceptor) = Idonor × Iacceptor

(14)

The donor and acceptor in the hydrogen-bond complex are generally thought to interact mainly via electrostatic interaction,52 which can be roughly described by eq 15: E HB = −Q donor × Q acceptor /dHB

(15)

where Qdonor and Qacceptor are the hypothetical effective charges of the donor and acceptor, respectively, and dHB is the optimal hydrogen-bond distance for the given hydrogen-bond complex in equilibrium, which is constant and can be included implicitly, as shown in eq 16: ′ ′ −E HB = Q donor × Q acceptor

(16)

where both Qdonor ′ and Qacceptor ′ are assumed to bear negative values for the intuitive description that follows. A direct comparison between eqs 14 and 16 indicates that the hydrogenbonding indices can be interpreted as the hypothetical effective charges of the donor and acceptor and also suggests that our assumption about eq 9 is reasonable. 1543


Article


The order of hydrogen-bonding indices are also consistent with the findings of Hao,7 who computed the hydrogenbonding strength of three acceptors and showed that ketone (act) is stronger than furan (frn) but weaker than amide (ace). This rank agrees with the order of hydrogen-bonding indices ace (Iacceptor = −2.52) < act (Iacceptor = −2.02) < frn (Iacceptor = −1.48) in our study. Additionally, they concluded that an sp2hybridized nitrogen atom is a stronger donor than an sp2hybridized oxygen atom in the oxazole group. The hydrogenbonding indices in our study indeed illustrate that Iacceptor(py) is lower than Iacceptor(frn). Thus, hydrogen-bonding indices can be used as straightforward and qualitative indicators of the capabilities of various common donors and acceptors in hydrogen-bond systems. 3.7. Implementation of the Hydrogen-Bonding Index as a Prefiltering Criterion for Virtual Screening. In the drug discovery process, medicinal chemists often use the strategy of maximizing the hydrogen-bonding interactions between the protein and ligand to improve the binding affinity or specificity of the ligand.36 Besides the widely used Lipinski’s rule of five,53 the hydrogen-bonding index, reflecting the individual contribution of the donor or acceptor to the hydrogen-bonding strength, is suggested to be used for additional prefiltering of the ligand database in virtual screening studies. For instance, in the case that the key residue in the protein structure is identified as a hydrogen-bond donor, ligands with lower Iacceptor, possibly conferring stronger hydrogen-bonding strength, will be preferred in the smallmolecule database for the subsequent virtual screening. Thus, it will help us save computational resources in advance if the hydrogen-bonding index is adopted as a prefiltering criterion to keep the small molecules with lower hydrogen-bond indices, which will possibly maximize the hydrogen-bonding strength between the protein and ligand. For the convenience of the user, we have developed a module called the “Hydrogen-Bonding Index Estimator (HBIE)” to rapidly and approximately estimate the hydrogenbonding indices of small molecules. The HBIE module has been integrated into our previous software, the Visual Force Field Derivation Toolkit (VFFDT), which was originally designed to derive metal-related force field parameters for AMBER in the 3D molecular viewer.42 In this work, we implemented the following procedure in the HBIE module: (1) The donor or acceptor types in Table S14 were defined on the basis of the model compounds in this work; for types of donors or acceptors that were not considered in our current model systems, the hydrogen-bonding indices of those missing types were estimated from those of similar types. (2) The donor or acceptor types in the ligand in which the user is interested are automatically checked. (3) If one or several types of donors or acceptors exist in the single ligand, the maximum hydrogenbonding index is returned, and the matching donor hydrogen atom or acceptor atom is highlighted in the 3D viewer; otherwise, the hydrogen-bonding index of this ligand is assigned as 0.0 (Figure 9). (4) A batch function is provided to automatically screen the ligand database to obtain the hydrogen-bonding index for each ligand (Figure 9). The program, PDF manual, and all of the related tutorials can be freely downloaded from the Dropbox shared folder at the following Web site: https://www.dropbox.com/sh/ hlefz8xy8o23q44/AADZGTBN-aJ07flZ-kM0s4U6a?dl=0. The hydrogen-bonding index offers a complementary prefiltering criterion to the number of hydrogen bonds, which

Figure 8. Hydrogen-bonding indices for donor (Idonor) and acceptor (Iacceptor). These indices were derived from MP2-computed hydrogenbonding strengths by the MLR method.

observed from the hydrogen-bonding index will be given as follows. The hydrogen-bonding indices in Figure 8 can explain the major tendencies observed in the MP2-computed hydrogenbonding strengths (Table 2) for the hydrogen-bond model systems (Tables 1). As can be seen from Table 2, hydrogen bonds formed by sp2-hybridized nitrogen atoms as acceptors are stronger than the corresponding ones formed by sp2- or sp3hybridized oxygen atoms as acceptors. For example, the acceptor py (Iacceptor = −2.96) is stronger than the acceptors met (Iacceptor = −2.48), act (Iacceptor = −2.02), and wat (Iacceptor = −1.98). An opposite trend can be also found in hydrogen-bond donors. Hydrogen bonds formed by sp3-hybridized nitrogen atoms as donors are weaker than those formed by sp3hybridized oxygen atoms as donors. For instance, the donors eto (Idonor = −2.03) and meo (Idonor = −2.04) are stronger than the donor eta (Idonor = −1.20). In addition, as indicated by the Idonor values, a nearby phenyl group can enhance the capability of donors. Figure 8 shows that donors containing aryl groups, phy (Idonor = −2.99), aln (Idonor = −2.04), and bnt (Idonor = −1.90), are stronger than the corresponding ones with aliphatic groups, meo (Idonor = −2.04), eta (Idonor = −1.20), and ett (Idonor = −1.14). 1544


Article


ln(−Idonor)) and Zacceptor (or ln(−Iacceptor)), which can be derived by MLR. These concepts such as Zdonor−acceptor, Zdonor, Zacceptor, Idonor, and Iacceptor have been interpreted and rationalized. Furthermore, the hydrogen-bonding indices (Idonor and Iacceptor) can be used as qualitative indicators. Moreover, the hydrogen-bonding index is also suggested as an additional prefiltering criterion for virtual screening besides Lipinski’s rule of five. Finally, we have implemented the “Hydrogen-Bonding Index Estimator (HBIE)” module in our VFFDT software to quickly estimate the hydrogen-bonding indices of ligands for database screening.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00022. Figures S1−S110 and Tables S1−S14 (PDF)

Figure 9. Hydrogen-Bonding Index Estimator (HBIE) module to rapidly estimate the hydrogen-bonding index (I) with two approaches. (1) When the user clicks on the menu “Protocol → HydrogenBonding Index Estimator (HBIE) → HBI(donor)”, Idonor is automatically printed in the window titled “Console Log” and the donor hydrogen is highlighted in yellow. (2) When the user inputs the command “cal HBI donor” in the console window, the same function as in the first method is performed. The Mol List module is implemented to obtain the hydrogen-bonding indices for the ligand database in batch.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. ORCID

Renxiao Wang: 0000-0003-0485-0259 Fu Lin: 0000-0002-9680-1536

is one part of Lipinski’s rule of five.53 However, two points are worth mentioning about the hydrogen-bonding index in prefiltering for virtual screening. First, the geometrical average of the hydrogen-bonding indices of some similar types of donors or acceptors is used instead for the missing types of donors or acceptors, since the types of donor and acceptor are so diverse in large small-molecule databases such as SPECS and ZINC54 that they may be not fully covered by our limited hydrogen-bond model systems in this work. In the future, we will adopt this protocol to derive more hydrogen-bonding indices to consider versatile donor or acceptor types, which requires much more effort and computational resources. (2) Second, the hydrogen-bonding index reflects the intrinsic contribution of the donor or acceptor to the hydrogen-bonding strength in the ideal geometry of the hydrogen bond with the optimal hydrogen-bond distance and the two hydrogen-bond angles both fixed at 180°. In reality, the ligand may not form the ideal hydrogen bond with the protein in the crystal structure. Thus, it should be noted that the hydrogen-bonding index reflects the maximum contribution of the donor or acceptor to the hydrogen-bonding strength; the more negative hydrogen-bonding index of a ligand means that the ligand may possibly form stronger hydrogen-bonds with a given donor or acceptor from the protein side.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the support of startup funding of Wenzhou Medical University, National Natural Science Foundation of China (Grant No. 21502144), and Zhejiang Provincial Natural Science Foundation of China (Grant No. LY17B030007), Ministry of Science and Technology of China (Grant No. 2016YFA0502302 to Prof. Wang), National Natural Science Foundation of China (Grant No. 21673276 to Prof. Wang), Chinese Academy of Sciences (Grant No. XDB20000000 to Prof. Wang). The authors thank Prof. Arthur Grollman at Stony Brook University and Dr. Yibing Wang at UT Southwestern Medical Center for their kind suggestions. The authors also thank the reviewers for their valuable comments and suggestions.

■

REFERENCES

(1) Huggins, M. L. 50 Years of Hydrogen Bond Theory. Angew. Chem., Int. Ed. Engl. 1971, 10, 147−152. (2) Pimentel, G. C.; McClellan, A. L. Hydrogen Bonding. Annu. Rev. Phys. Chem. 1971, 22, 347−385. (3) Kollman, P. A.; Allen, L. C. Theory of The Hydrogen Bond. Chem. Rev. 1972, 72, 283−303. (4) Šponer, J.; Leszczynski, J.; Hobza, P. Hydrogen Bonding And Stacking of DNA Bases: A Review of Quantum-chemical ab initio Studies. J. Biomol. Struct. Dyn. 1996, 14, 117−135. (5) Perrin, C. L.; Nielson, J. B. Strong” Hydrogen Bonds in Chemistry and Biology. Annu. Rev. Phys. Chem. 1997, 48, 511−544. (6) Martin, T. W.; Derewenda, Z. S. The Name is Bond-H Bond. Nat. Struct. Biol. 1999, 6, 403−406. (7) Hao, M. H. Theoretical Calculation of Hydrogen-Bonding Strength for Drug Molecules. J. Chem. Theory Comput. 2006, 2, 863− 872. (8) Abraham, M. H.; Duce, P. P.; Prior, D. V.; Barratt, D. G.; Morris, J. J.; Taylor, P. J. Hydrogen Bonding. Part 9. Solute Proton Donor And

4. CONCLUSIONS In this study, we have conducted intensive QM computations with MP2 and DFT methods and have shown that an exponential relationship is held between the hydrogen-bonding strength (EHB) and the optimal hydrogen-bond distance (dHB) and also exists between the hydrogen-bonding index (I) and the hypothetical effective hydrogen-bond radius (r). For more intuitive demonstration, the magnitude scale for the hydrogenbond complex, Zdonor−acceptor, holds the linear relationship with the optimal hydrogen-bond distance dHB with R2 values ranging from 0.80 to 0.90. Our work has indicated that Zdonor−acceptor is an additive quantity and is simply the sum of Zdonor (or 1545


Article

Journal of Chemical Information and Modeling Proton Acceptor Scales For Use In Drug Design. J. Chem. Soc., Perkin Trans. 2 1989, 1355−1375. (9) Arnaud, V.; Berthelot, M.; Le Questel, J. Y. Hydrogen-bond Accepting Strength of Protonated Nicotine. J. Phys. Chem. A 2005, 109, 3767−3770. (10) Laurence, C.; Berthelot, M. Observations on the Strength of Hydrogen Bonding. Perspect. Drug Discovery Des. 2000, 18, 39−60. (11) Laurence, C.; Brameld, K. A.; Graton, J.; Le Questel, J.-Y.; Renault, E. The pKBHX Database: Toward a Better Understanding of Hydrogen-Bond Basicity for Medicinal Chemists. J. Med. Chem. 2009, 52, 4073−4086. (12) Taft, R. W.; Gurka, D.; Joris, L.; Schleyer, P. v. R.; Rakshys, J. W. Studies of Hydrogen-Bonded Complex Formation With p-fluorophenol. V. Linear Free Energy Relationships With OH Reference Acids. J. Am. Chem. Soc. 1969, 91, 4801−4808. (13) Abraham, M. H.; Grellier, P. L.; Prior, D. V.; Morris, J. J.; Taylor, P. J. Hydrogen Bonding. Part 10. A Scale of Solute HydrogenBond Basicity Using log K Values For Complexation In Tetrachloromethane. J. Chem. Soc., Perkin Trans. 2 1990, 521−529. (14) Oliveira, B. G.; Pereira, F. S.; de Araújo, R. C. M. U.; Ramos, M. N. The Hydrogen Bond Strength: New Proposals to Evaluate the Intermolecular Interaction Using DFT Calculations and the AIM Theory. Chem. Phys. Lett. 2006, 427, 181−184. (15) Reynisson, J.; McDonald, E. Tuning of Hydrogen Bond Strength Using Substituents on Phenol and Aniline: a Possible Ligand Design Strategy. J. Comput.-Aided Mol. Des. 2004, 18, 421−431. (16) El Kerdawy, A.; Tautermann, C. S.; Clark, T.; Fox, T. Economical And Accurate Protocol For Calculating Hydrogen-BondAcceptor Strengths. J. Chem. Inf. Model. 2013, 53, 3262−3272. (17) Grabowski, S. J. An Estimation of Strength of Intramolecular Hydrogen Bonds - ab initio and AIM Studies. J. Mol. Struct. 2001, 562, 137−143. (18) Grabowski, S. J. Hydrogen Bonding Strength-Measures Based on Geometric and Topological Parameters. J. Phys. Org. Chem. 2004, 17, 18−31. (19) Ireta, J.; Neugebauer, J.; Scheffler, M. On the Accuracy of DFT for Describing Hydrogen Bonds: Dependence on the Bond Directionality. J. Phys. Chem. A 2004, 108, 5692−5698. (20) Leo, A. J. Evaluating Hydrogen-bond Donor Strength. J. Pharm. Sci. 2000, 89, 1567−1578. (21) Mohajeri, A.; Nobandegani, F. F. Detection and Evaluation of Hydrogen Bond Strength in Nucleic Acid Base Pairs. J. Phys. Chem. A 2008, 112, 281−295. (22) Schwöbel, J.; Ebert, R.-U.; Kühne, R.; Schüürmann, G. Prediction of the Intrinsic Hydrogen Bond Acceptor Strength of Organic Compounds by Local Molecular Parameters. J. Chem. Inf. Model. 2009, 49, 956−962. (23) Abraham, M. H.; Ibrahim, A.; Zissimos, A. M.; Zhao, Y. H.; Comer, J.; Reynolds, D. P. Application of Hydrogen Bonding Calculations in Property Based Drug Design. Drug Discovery Today 2002, 7, 1056−1063. (24) Gancia, E.; Montana, J. G.; Manallack, D. T. Theoretical Hydrogen Bonding Parameters for Drug Design. J. Mol. Graphics Modell. 2001, 19, 349−362. (25) Glavatskikh, M.; Madzhidov, T.; Solov’ev, V.; Marcou, G.; Horvath, D.; Graton, J.; Le Questel, J. Y.; Varnek, A. Predictive Models for Halogen-bond Basicity of Binding Sites of Polyfunctional Molecules. Mol. Inf. 2016, 35, 70−80. (26) Grabowski, S. J. A New Measure of Hydrogen Bonding StrengthAb Initio and Atoms in Molecules Studies. Chem. Phys. Lett. 2001, 338, 361−366. (27) Kenny, P. W.; Montanari, C. A.; Prokopczyk, I. M.; Ribeiro, J. F. R.; Sartori, G. R. Hydrogen Bond Basicity Prediction for Medicinal Chemistry Design. J. Med. Chem. 2016, 59, 4278−4288. (28) Ruggiu, F.; Solov’ev, V.; Marcou, G.; Horvath, D.; Graton, J.; Le Questel, J.-Y.; Varnek, A. Individual Hydrogen-Bond Strength QSPR Modelling with ISIDA Local Descriptors: a Step Towards Polyfunctional Molecules. Mol. Inf. 2014, 33, 477−487.

(29) Schwobel, J.; Ebert, R. U.; Kuhne, R.; Schuurmann, G. Prediction of the Intrinsic Hydrogen Bond Acceptor Strength of Chemical Substances From Molecular Structure. J. Phys. Chem. A 2009, 113, 10104−10112. (30) Green, A. J.; Popelier, P. L. A. Theoretical Prediction of Hydrogen-Bond Basicity pKBHX Using Quantum Chemical Topology Descriptors. J. Chem. Inf. Model. 2014, 54, 553−561. (31) Graton, J.; Le Questel, J. Y.; Maxwell, P.; Popelier, P. HydrogenBond Accepting Properties of New Heteroaromatic Ring Chemical Motifs: A Theoretical Study. J. Chem. Inf. Model. 2016, 56, 322−334. (32) Murray, J. S.; Politzer, P. Correlations Between the Solvent Hydrogen-Bond-Donating Parameter.Alpha. And The Calculated Molecular Surface Electrostatic Potential. J. Org. Chem. 1991, 56, 6715−6717. (33) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by The Differences of Separate Total Energies. Some Procedures With Reduced Errors. Mol. Phys. 1970, 19, 553−566. (34) Lamarche, O.; Platts, J. A. Theoretical Prediction of the Hydrogen-Bond Basicity pKHB. Chem. - Eur. J. 2002, 8, 457−466. (35) Grabowski, S. J. Ab Initio Calculations on Conventional And Unconventional Hydrogen BondsStudy of The Hydrogen Bond Strength. J. Phys. Chem. A 2001, 105, 10739−10746. (36) Nocker, M.; Handschuh, S.; Tautermann, C.; Liedl, K. R. Theoretical Prediction of Hydrogen Bond Strength for Use in Molecular Modeling. J. Chem. Inf. Model. 2009, 49, 2067−2076. (37) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (38) Becke, A. D. Density-functional thermochemistry. I. The Effect of The Exchange-Only Gradient Correction. J. Chem. Phys. 1992, 96, 2155−2160. (39) Lee, C.; Yang, W.; Parr, R. G. Development of The ColleSalvetti Correlation-Energy Formula Into A Functional of The Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (40) Kendall, R. A.; Dunning, J. T. H.; Harrison, R. J. Electron Affinities of The First-row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (41) Xu, X.; Goddard, W. A. The X3LYP Extended Density Functional for Accurate Descriptions of Nonbond Interactions, Spin States, and Thermochemical Properties. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 2673−2677. (42) Zheng, S.; Tang, Q.; He, J.; Du, S.; Xu, S.; Wang, C.; Xu, Y.; Lin, F. VFFDT: A New Software for Preparing AMBER Force Field Parameters for Metal-Containing Molecular Systems. J. Chem. Inf. Model. 2016, 56, 811−818. (43) Grimme, S. Density Functional Theory With London Dispersion Corrections. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 211−228. (44) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate ab initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (45) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. 1546


Article

Journal of Chemical Information and Modeling (46) Kruse, H.; Goerigk, L.; Grimme, S. Why The Standard B3LYP/ 6-31G* Model Chemistry Should Not Be Used In DFT Calculations of Molecular Thermochemistry: Understanding And Correcting The Problem. J. Org. Chem. 2012, 77, 10824−10834. (47) Møller, C.; Plesset, M. S. Note On An Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618−622. (48) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, And Transition Elements: Two New Functionals And Systematic Testing of Four M06-class Functionals And 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (49) Wendler, K.; Thar, J.; Zahn, S.; Kirchner, B. Estimating The Hydrogen Bond Energy. J. Phys. Chem. A 2010, 114, 9529−9536. (50) Zou, J. W.; Huang, M.; Hu, G. X.; Jiang, Y. J. Toward a Uniform Description of Hydrogen Bonds and Halogen Bonds: Correlations of Interaction Energies with Various Geometric, Electronic and Topological Parameters. RSC Adv. 2017, 7, 10295−10305. (51) Grabowski, S. J. Ab initio and AIM Studies on Measures of Hydrogen Bonding Strength R-CN···HF and R-CN···HCl Complexes. J. Mol. Struct. 2002, 615, 239−245. (52) Pauling, L. The Nature of the Chemical Bond and the Structure of Molecules and Crystals: An Introduction to Modern Structural Chemistry; Cornell University Press: Ithaca, NY, 1960. (53) Lipinski, C. A. Lead- and Drug-like Compounds: the Rule-ofFive Revolution. Drug Discovery Today: Technol. 2004, 1, 337−341. (54) Irwin, J. J.; Sterling, T.; Mysinger, M. M.; Bolstad, E. S.; Coleman, R. G. ZINC: A Free Tool to Discover Chemistry for Biology. J. Chem. Inf. Model. 2012, 52, 1757−1768.

1547


Proposed Hydrogen-Bonding Index of Donor or Acceptor Reflecting Its

Recommend Documents