Elucidating Hyperconjugation from Electronegativity to Predict Drug

Article pubs.acs.org/jcim

Elucidating Hyperconjugation from Electronegativity to Predict Drug Conformational Energy in a High Throughput Manner Zhaomin Liu,† Joshua Pottel,† Moeed Shahamat,† Anna Tomberg,† Paul Labute,‡ and Nicolas Moitessier*,† †

Department of Chemistry, McGill University, 801 Sherbrooke St. W., Montréal, QC, Canada H3A 0B8 Chemical Computing Group Inc., 1010 Sherbrooke St. W., Montréal, QC, Canada H3A 2R7

‡

S Supporting Information *

ABSTRACT: Computational chemists use structure-based drug design and molecular dynamics of drug/protein complexes which require an accurate description of the conformational space of drugs. Organic chemists use qualitative chemical principles such as the effect of electronegativity on hyperconjugation, the impact of steric clashes on stereochemical outcome of reactions, and the consequence of resonance on the shape of molecules to rationalize experimental observations. While computational chemists speak about electron densities and molecular orbitals, organic chemists speak about partial charges and localized molecular orbitals. Attempts to reconcile these two parallel approaches such as programs for natural bond orbitals and intrinsic atomic orbitals computing Lewis structures-like orbitals and reaction mechanism have appeared. In the past, we have shown that encoding and quantifying chemistry knowledge and qualitative principles can lead to predictive methods. In the same vein, we thought to understand the conformational behaviors of molecules and to encode this knowledge back into a molecular mechanics tool computing conformational potential energy and to develop an alternative to atom types and training of force fields on large sets of molecules. Herein, we describe a conceptually new approach to model torsion energies based on fundamental chemistry principles. To demonstrate our approach, torsional energy parameters were derived on-the-fly from atomic properties. When the torsional energy terms implemented in GAFF, Parm@Frosst, and MMFF94 were substituted by our method, the accuracy of these force fields to reproduce MP2-derived torsional energy profiles and their transferability to a variety of functional groups and drug fragments were overall improved. In addition, our method did not rely on atom types and consequently did not suffer from poor automated atom type assignments.

■

INTRODUCTION Drug Conformational Space. Drugs are built around a variety of heterocycles (e.g., imidazole, oxazole, carbohydrates) and functional groups (e.g., hydroxamic acids, polyols). However, this diversity still covers only a small fraction of the drug structure space with privileged structures being overrepresented. As a result, the diversity of these privileged structures and libraries has been looked at computationally and experimentally.1−9 Although increasing the diversity of the chemical databases is expected to increase the hit discovery rate, this ever-increasing diversity is a major concern in computational chemistry and more specifically molecular mechanics as large diversity correlates with a need for highly transferable methods. Computationally Rationalizing Observations. On the computational chemistry side, it is critical to accurately and quickly predict the preferred conformations of a given molecule and to compute their relative potential energy. Currently, many medicinal chemistry projects involve modeling tools such as docking of small molecules to proteins10 and molecular dynamics simulations.11 To carry out these calculations in a high throughput manner, molecular mechanics (MM) has been © XXXX American Chemical Society

employed. In MM, potential energy is computed through a set of mathematical equations accounting for bond stretching, angle bending, torsional rotation (e.g., eq 1) and nonbonded interactions. These equations rely on precomputed parameters describing the atoms and their interactions. The combination of equations and parameters are referred to as force fields (FFs). E = E bonds + Eangles + Etorsions + EvanderWaals + Eelectrostatics (1)

There are a number of widely used FFs including the AMBER12−15 series, CHARMM,16,17 MMFF94,18−22 GROMOS,23−28 and OPLS.29−31 These are built around a limited number of atom types (e.g., hydroxyl group oxygen, sp3 carbon) which are assigned to each atom of the molecules under investigation. However, several functional groups within uncommon atomic environments are not specifically parametrized and the accuracy may drop when such functional groups are modeled by FFs. To address this issue, research groups categorize and specifically parametrize classes of Received: January 11, 2016

A

DOI: 10.1021/acs.jcim.6b00012 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

rationalize the outcome of chemical transformations or the properties of some molecules. Interestingly, most are implying hyperconjugation. Hybrid Method. While principles used by organic chemists are well transferable, they are qualitative. While force fields are quantitative, they are poorly transferable. An alternative approach to MM is to quantify these chemical principles. Thus, we thought to revisit the use of FFs in MM and to develop a method that personalizes parameters for molecules by quantifying the principles discussed above. This will ultimately lead to a method that would not rely on very large set of parameters or atom types. Methods such as natural bond orbitals52 and intrinsic bond orbitals53 have emerged, which connect quantum mechanical techniques to more intuitive aspects of reaction mechanisms and Lewis structures. While bond stretching and angle bending cannot be ignored, the primary factor controlling conformational preferences are torsions and accurately computing torsional energy profiles is critical. As of today, torsional energy is computed as a combination of van der Waals (vdW) and electrostatic interactions corrected with an empirical torsion energy term (Figure 2). The aim of this study is not to develop another FF through fitting parameters to experimental or high-level calculations data. Rather, we thought to investigate the physical origins of the empirical term used in FF, and to quantify it to serve as a parameter-free method for MM. The transferability of this method is, in theory, much higher than the conventional methodology. Herein, we wish to report our efforts toward an accurate non-atom type-based MM method named H-TEQ. We first identified hyperconjugation to be the major contributor to the empirical torsion terms found in FFs. Then we developed a method that quantifies hyperconjugation contribution to the torsional energy and a mechanistic approach to compute it.

molecules: for example, carbohydrates feature several polar, hydrogen bond donor, and acceptor hydroxyl groups that are close to each other and may form a complex hydrogen bond network as observed in cyclodextrin.32 These strong electronic effects together with the anomeric effect are not properly modeled with common parameters, thus, specialized force fields (e.g., GLYCAM33,34) were developed or carbohydrate parameters were added to the existing FFs.35−37 Similarly, remote electron-withdrawing or donating substituents can significantly modulate the potential energy surface, and, consequently, the MM descriptors (parameters) trained on a limited number of existing molecules would not be suitable.38 One way of rectifying the transferability issue is to constantly complement the existing FFs as shown by the development of Parm@Frostt spanning nearly two decades.39 To cover the enormous set of functional groups, these atom type-based methods have to rely on an ever-growing number of parameters. The total number of possible small organic molecules is likely to exceed 1060 compounds with a large diversity.40 Training FFs to cover the entire chemical space is expectedly impossible. To avoid this uncertainty, others have attempted to develop FFs covering more of the chemical space, by developing parameters on-thefly: DREIDING (1990),41 Universal Force Field (UFF, 1992),42 MAB (1998),43,44 and Extensible and Systematic Force Field (ESFF, 2003),45 for example. On one hand, these FFs extend to a wider range of systems: UFF covers most of the periodic table, and ESFF was designed to model organic, inorganic, and organometallic systems; on the other hand, the accuracy of modeling organic molecules, most common in drug discovery, is not satisfactory.46−48 Notable Chemical Principles. On the experimental side, a number of molecular conformational preferences have been rationalized by dipole−dipole interaction, hyperconjugation, and steric clashes. For example, the preference of glucose for the more sterically encumbered α configuration has been explained by the anomeric effect, which results from a favorable overlap between the endocyclic oxygen lone pair and the antibonding orbital of the exocyclic C−O bond. This hyperconjugation stabilization together with the more favorable dipole−dipole interaction overcomes the steric preference for the β configuration (Figure 1a). A similar hyperconjugation

■

UNDERSTANDING THE CHEMISTRY ORIGINS Fundamental Rationales. In the past, we have shown that encoding chemical principles such as the Hammond−Leffler postulate and Curtin−Hammett principle to predict the outcome of asymmetric reactions was possible and effective.54,55 We have also shown that implementing chemistry knowledge such as proton transfers56 and covalent bond formation 57 in the process of drug−enzymes binding significantly improves the prediction of small ligand−protein recognition. Following the same strategy, we thought that refining MM guided by chemistry knowledge would fulfill our goal of both high accuracy of torsional profile prediction and transferability to a wide range of drug molecules. More specifically, chemical and physical properties, such as hyperconjugation (explaining the gauche and anomeric effect), conjugation (explaining the conformation preference of butadiene for the s-trans conformation), steric hindrance, polarizability (such as used in hard and soft acids and bases (HSAB) theory58), inductive effects (explaining the regioselectivity of aromatic substitutions together with resonance effects), and electrostatics, have been qualitatively used by chemists for decades to explain conformational behaviors of small molecules or complexes. Some of them, such as vdW, electrostatics, and polarizability are already included in FFs in a quantitative way. Consequently, we hypothesized that quantifying other interactions and properties could be used to more accurately model torsional energy profiles. As the first part of

Figure 1. Conformational preferences rationalized by chemical principles: (a) anomeric effect; (b) gauche effect.

effect has been proposed to explain the gauche effect (e.g., preference for the gauche conformation of 1,2-difluoroethane, Figure 1b). Predictive models have been developed which are based on conformational preferences such as the wellestablished Felkin−Anh model predicting the stereochemical outcome of additions to carbonyl derivatives49,50 or the Cieplak model rationalizing the stereoselective reduction of cyclohexanone with “small” hydrides.51 These principles and models have been taught and used for decades by organic chemists to B


Article


Figure 2. H-TEQ, an alternative to MM force fields using chemical principles.

compared to profiles computed in MM using a widely used force field, GAFF. The torsion term within GAFF was then replaced by a hyperconjugation term while keeping the electrostatic and vdW unchanged. In this first part, the hyperconjugation energy was computed for each conformation using the Natural Bond Orbital (NBO) program. As shown in Figure 3, while both approaches predicted a gauche effect, the

our investigations, we describe herein our proof-of-principle studies on non-conjugated systems. Hyperconjugation as a Major Torsional Energy Term. Since the conformations of molecules are primarily described by dihedral angles, an accurate computation of the torsion energy profiles is critical and will be the focus of this study. Ethane derivatives were a good starting case to study in detail. The major causes of the ethane rotational barrier have been previously investigated: Pophristic and Goodman claimed that hyperconjugation is, for the most part, responsible for the preferred staggered conformation by applying natural bond orbital (NBO) analysis;59 whereas another study from Mo and Gao stated the conventional steric repulsions dominate this behavior based on the block-localized wave function (BLW) method.60 We believe that the torsion energy is actually controlled by a combination of hyperconjugation, electrostatic, and steric effects and that the impact of each component would vary depending on the molecule. The latter two have been implemented in most current FFs (fourth and fifth terms in eq 1) and complemented by a specific torsion term (third term in eq 1). Parameters for this additional term are traditionally trained on sets of molecules. We hypothesize that incorporating a hyperconjugation term, as a substitution to the current implemented empirical torsion term in most FFs, would improve torsion modeling. It should be noted that this manuscript is intended to focus on hyperconjugation, one component of the torsional energy. In order to test our idea, we use the current FF platform/scheme and take the other terms as available. It may be suggested that atom type-based FFs implicitly include some hyperconjugation effects through their training/parametrization. Although MM4 encoded hyperconjugation as an extra correction for bond stretching,61 to the best of our knowledge, hyperconjugation, as traditionally applied to evaluate conformational preferences, has not been explicitly considered for describing torsional energies in FFs, and could lead to improved results. As a first test of our hypothesis, 1,2-difluoroethane, exhibiting a well-known gauche effect, was taken as an example. The rotational energy profile was computed in quantum mechanics (QM) at the MP2/6-311+G** level as the reference and

Figure 3. Torsional energy for 1,2-difluoroethane computed with a FF (GAFF), QM (MP2/6-311+G**), and our MM method using our hyperconjugation term.

hyperconjugation method predicted the high rotation barrier at F−C−C−F angle of 0° while GAFF poorly reproduced the MP2 torsional energy profile. The hyperconjugation contribution has been tested on a set of molecules and will be discussed below. As shown in Figure 3, hyperconjugation appeared as a major player in the torsional energy. Using a schematic representation of the sp3−sp3 rotatable bond (Figure 4), hyperconjugation takes place between the bonding orbital σ of D−X and the antibonding orbital σ* of Y−A and is denoted σ → σ*. Thus, (1) the role of donor (D) and acceptor (A) must be fully understood; (2) the chemical nature of the central bond X−Y can also impact the hyperconjugation and should also be C


Article


Figure 4. Hyperconjugation: interaction between orbitals between a donor (D) and an acceptor (A).

considered; (3) last, the influence of the chemical environment (functional groups R1−4), should be accounted for as well. HyperconjugationQualitative Evaluation. Organic chemists often relate hyperconjugation to electronegativity (χ): more electronegative atoms are better hyperconjugation acceptors and worse hyperconjugation donors. A previously reported computational study has confirmed this correlation.62 Moreover, the size, the softness, and the polarizability of the element are likely other factors impacting hyperconjugation (the acceptor ability increases for the halogens when going down the periodic table). This constituted the basis of our approach: the hyperconjugation donating/accepting ability can be represented by simple chemical propertyelectronegativity of the atoms constituting the torsion. However, previous studies and current chemistry knowledge are not enough to help us to build it into MM: while the hyperconjugation in trans conformation was scrutinized, the impact of the conformation (rotation) on the hyperconjugation energy has never been carefully evaluated. The next task was therefore to understand the relationship between hyperconjugation and conformation. The hyperconjugation acceptor and donor abilities were primarily studied on substituted C−C bonds (ethane derivatives62 and cyclohexyl cations)63 while, the central bonds made of atoms other than carbons have been seldom investigated.62 A previous study showed that the H−C−X−C hyperconjugation energy (Ehyp) in the trans conformation was similar whether X was a carbon or an oxygen (Ehyp = −3.20 to −3.36 kcal/mol) but dropped when the C−X bond length increased (X = Se, Ehyp = −1.21 kcal/mol; X = S, Ehyp = −1.60 kcal/mol).62 To further understand the central bond impact on hyperconjugation, we proposed to compare H−C−C−X systems to the highly polarized H−C−N(+)−X and H−N(+)− C−X counterparts as test cases. These molecules will be used first to evaluate frontier molecular orbital (FMO)-derived principles and to propose trends. Origins of Hyperconjugation: The Orbital Overlap and Energy Gap. The strength of the hyperconjugation energy was dictated by the orbital overlap and the energy gap between the σ (donor) bond and the σ* (acceptor) bonds. When considering the FMO theory, as more electronegative atoms have atomic orbitals lower in energy, the presence of an electronegative A within D−C1−C2−A lowers C2−A σ* making it a strong hyperconjugation acceptor. On the donor side, an electronegative D, for example fluorine, lower the energy of σ (D−C1), making a larger energy gap between σ and σ*, ultimately leading to a weaker interaction. However, the impact of substitution on the central bond is not straightforward: the substitution of a central bond carbon at the acceptor side by a more electronegative N(+) atom is expected to lower the energy of σ*, making N(+)−X a better acceptor (Figure 5). In the meantime, the presence of an electron-withdrawing atom (N(+)) bound to the carbon in H−C−N(+)−X should also make it electron poor (through inductive effects), hence more

Figure 5. Expected impact (qualitative) of the central bond polarization on the orbital energies.

electronegative. This should in turn result in a σ lower in energy making H−C a poorer donor. While the former effect is expected to increase the hyperconjugation stabilization, the latter is expected to reduce it and the net impact of electronegativity is difficult to predict from these qualitative observations (Figure 5). In order to better comprehend the interplay between these effects, we selected ethane, monofluoroethane, and monofluoropropane, substituted the central carbon to nitrogen, and computed the necessary data. In the past, it has been proposed that the hyperconjugation can be predicted by looking at the overlap between FMOs hence quantifying a well-established principle. Thus, the orbital overlap, bonding/antibonding energy, and hyperconjugation energy were computed using the NBO package,64 which provides comprehensive tools for the analysis of molecular orbitals (Table 1). The collected data for terminal atom substitutions was consistent with FMO predictions that more electronegative atoms lower the orbital energies: H → F substitution on the acceptor site decreased εσ* from 0.62 to 0.49; H → CH3 on the donor side lowers εσ (−0.81 vs −0.69). The presence of a more electronegative carbon and a more electronegative N(+) lowered both the σ and σ* relatively to the C−C analog, as expected. It also appeared that the presence of a more electronegative carbon in H−C−N−F and H−N−C−F and compared to H−C−C−F (Δεσ H−C = −0.27; Δεσ* C−F = −0.18) was less impacting than the presence of a highly electronegative nitrogen (Δεσ H− N = −0.47; Δεσ* N−F = −0.40). This data also revealed that the donor ability (Δεσ = −0.47, −0.26) was more affected than the acceptor ability (Δεσ* = −0.18, −0.40). If σ and σ* involved in this process are the HOMO and LUMO, this electron gap has been referred to as the absolute hardness and can be directly related to the HSAB theory.65 Overall, this data confirmed that the hyperconjugation is inversely proportional to the energy gap between the σ and σ* involved and could be used to quantify the hyperconjugation difference between different molecules. The orbital overlap is also expected to vary as a function of the bond polarization. According to the FMO theory, the bonding orbital resides more on the most electronegative atom. Thus, it is expected to be more localized (dense) near N(+) of H−N(+)−C−X and near the carbon of H−C−N(+)−X. This increased compactness of the bonding orbital would affect the overlap between the σ donor and σ* acceptor. Our computations confirmed that introducing N(+), on either the donor (H−N(+)−C−F) or acceptor site (H−C−N(+)−F), reduced the orbital overlap (Sij values in Table 1). One may notice that H−C−N(+)−F has lower orbital overlap than H− N(+)−C−F (Sij = 0.122 vs 0.168), rationalized by the polarization effect on σ*. In H−C−C/N(+)−F, the C−F D


Article

Journal of Chemical Information and Modeling Table 1. Energy Gap and Overlap for Hyperconjugation around C−C vs C−N(+) Bonds

a Bonding (εσ) or antibonding (εσ*) orbital energies computed using NBO. bChange in (εσ) or antibonding (εσ*) orbital energies relative to H−C− C−F or CH3−C−C−F. cEnergy gap between εσ and εσ* computed for a given molecule and εσ(H−C or H3C−X) or εσ*(C−F) computed for H/ CH3−C−C−F. dChange in energy gap relative to H−C−C−F or CH3−C−C−F. eSij term of the overlap matrix as computed with NBO. f Hyperconjugation energy computed through the second-order perturbation of the bonding−antibonding interactions in NBO.66

observed, resulting in an energy minimum. As the angle reached 0°, the conformation would be a local energy minimum, due to the partial orbital overlap (Sij = 0.089, Table 1). This directionality of hyperconjugation makes it crucial describing the torsion rotation. Hyperconjugation in Different Chemical Environments. With different functional groups around, hyperconjugation might behave differently. We used pentafluoroethane and monofluoroethane and evaluate the same CH → CF* with different Rs at the two ends of the central bond. It can be seen that the amplitude of the hyperconjugation was smaller, in the presence of the other four fluorine atoms in pentafluoroethane (Figure 7).

bond is highly polarized and the antibonding orbital resides predominantly on the least electronegative atom (C) and hyperconjugation with the H−C bonding σ bond is strong. The substitution of the C−F bond by N(+)−F and the higher electronegativity of N(+) leads to a less polarized σ* with a higher contribution from F. This modulation of the relative atomic contributions to the antibonding orbital is expected to reduce the hyperconjugation stabilization. This polarization effect on the energy gap and orbital overlap also appeared when the donor was a methyl group. This first study of the interplay of σ−σ* energy gap and overlap on the hyperconjugation stabilization overall confirmed that our FMO-derived qualitative principles were predictive and should now be quantified. Hyperconjugation upon Rotation. The hyperconjugation energy was computed upon torsion rotation (Figure 6), when the bonding and antibonding orbitals were perpendicular (±90° rotation) and the overlap was minimal (Sij = 0.028, Table 1), the donation was at a minimum, giving the highest energy barrier, whereas when they were parallel, more specifically antiperiplanar, the overlap was maximal (Sij = 0.182, Table 1) and the strongest hyperconjugation was

■

IMPLEMENTATION Torsional Energy As Combination of Hyperconjugation, Sterics, and Electrostatics. To further assess our hypothesis that hyperconjugation can substitute the torsional energy terms commonly used in MM, we first examined whether a simple combination of hyperconjugation energy with vdW and electrostatics can model torsional energy profiles using 261 ethane derivatives including common functional E


Article


by N(+) to generate 76 ammonium derivatives which were added to the set for a total of 337 small molecules. The selection was driven by the following objective: decoupling of the σ → σ* hyperconjugation effect from more complex phenomena such as conjugation, intramolecular interactions, and n → σ* hyperconjugation, would make our analysis simpler. Hyperconjugation profiles for each torsion in this set were computed. We first generated a set of conformations for each molecule by rotating the central bond in 10° increments and optimized them using a high level of quantum mechanics (MP2/6-311+G**). This MP2 calculation also provided torsional energy profiles for comparison, which were used as reference for evaluating the quality of our method. We then computed their corresponding potential energies using GAFF. In order to test the use of hyperconjugation energies in place of traditional torsion terms, we replaced the GAFF-computed torsional energies of all the torsions (i.e., 9 torsions per C−C/N(+) bond) involved in the central bond rotation by NBO-computed hyperconjugation energies. All the other energy contributions, which include electrostatic, vdW, bond stretching, angle bending, out-of-plane, and the torsions not involved in this torsional rotation computed by the GAFF were retained. The average root-mean-square deviation between GAFF and MP2 torsional energy profiles (RMSDGAFF) was computed as well as the RMSD between our hyperconjugation corrected FFs and MP2 profiles (RMSDNBO corrected FF) for the 337 small molecules. The measured deviations (RMSDGAFF = 1.62 kcal/mol and RMSDNBOcorrectedFF = 1.49 kcal/mol) demonstrated that substituting the torsion term by NBO-derived hyperconjugation improved the accuracy of GAFF when MP2 energies were used as reference. Other energy contributions are involved in the GAFF derived torsion energy profiles (vdW and electrostatics). In order to evaluate how the specific torsion parameters improved the torsional energy, we recomputed the deviation considering the vdW and electrostatic terms only, while excluding the torsion terms from the calculation (RMSDotherterms = 2.29 kcal/mol). This data demonstrated that the torsion term implemented in GAFF improved the RMSD by 0.67 kcal/mol, while the incorporation of NBO derived hyperconjugation improved the accuracy by 0.80 kcal/mol (19% more accurate energetic description for torsion parameters). Hyperconjugation Energy from NBO. It is noteworthy that a weighting scheme of 0.6 has been applied to the hyperconjugation term. Initially, no scaling was used when we assembled NBO data with vdW and electrostatics terms. However, we noticed that the gauche conformation was often favored. We proposed two explanations for this observation. As mentioned above, there is still a debate on how much of the torsional energy should be attributed to hyperconjugation and sterics when applying different theoretical methods for example NBO vs BLW. Unfortunately, while computational methods still disagree, no experimental methods enabled us to measure the hyperconjugation as a function of the torsion angles.67 Second, we added hyperconjugation to the vdW and electrostatics as implemented in current available MM platforms (for instance the comparison shown above is based on GAFF). FFs are developed in a self-consistent manner: different terms such as vdW, electrostatics as well as torsional parameters are balanced with respect to each other in their own system. An overpreferred gauche conformation can also be caused by the

Figure 6. Hyperconjugation energy profile. (top) From NBO analysis (CH → CH* in ethane). (bottom) Fourier transform into eq 2 with V1 = 2.17, V2 = −2.43, and V3 = 0.37; see text below for definitions of V1−3.

Figure 7. Hyperconjugation energy with same donor−acceptor pair can vary when placed in different molecules (hyperconjugation parameters for pentafluoroethane V1 = 2.45, V2 = −2.64, V3 = 0.21; for monofluoroethane: V1 = 2.27, V2 = −4.99, V3 = 0.32). See text for definitions of V1−3.

groups (Figure 8). A subset of these ethane derivatives was selected and one of the two central carbon atoms was replaced

Figure 8. Ethane derivatives investigated. F


Article

Journal of Chemical Information and Modeling Table 2. V1−3 Terms Computed for Representative Structures

underestimation of the vdW and electrostatics contributions when they are trained to optimize the overall performance on the molecules in their training set. As the objective of the current study is not developing a new FF from scratch, parameters of terms other than torsion are taken as they are. Different weights applied to the NBO-derived hyperconjugation energy (from 0.4 to 1 with 0.1 increments) have been examined (see the SI). An optimal value of 0.6 was selected as it leads to accurate results across three widely used FFs. Hyperconjugation as a FF Term. As the connection between hyperconjugation and torsion energy was elucidated, we sought to implement this term into existing FFs with the hope to remove the need for atom types and precomputed parameters. NBO analysis required the computation of molecular orbitals and significant computational resources. One cannot efficiently derive the parameters, using NBO, for each torsion of a given molecule, if screening more than one compound. Thus, deriving hyperconjugation without the need for NBO computations is the key to encode it into a useful method and a simple mathematical format was required. Hyperconjugation profiles were first transformed through Fourier expansion into an equation corresponding to torsion angles (eq 2) using the program R,68 followed by a leastsquares fitting to obtain V1−3. The form of this torsion term has been widely used in FF development making our approach applicable to all of these MM programs. The transformation did not lead to a significant accuracy change (ΔRMSD = 0.04 kcal/ mol), indicating that the torsional energy function commonly used in FFs can be used to model hyperconjugation. E hyp =

V1 V (1 + cos θ ) + 2 (1 + cos(2θ )) 2 2 V3 + (1 + cos(3θ )) 2

a

(2)

V1−3 values derived from NBO data.

other molecules and will certainly lead to improved accuracy for a variety of molecules. Generating the V1−3 TermsDonor/Acceptor Electronegativities. The simple mathematical form (eq 2) for hyperconjugation has been purposely developed in the same format of the torsion parameters in many FFs and the meanings of the V values have been elucidated. For increased transferability, V values should not be static and precomputed for any specific molecule, thus a method to quickly derive these V values considering the different chemical environments should be developed. Our hypothesis that qualitative chemical principles can be used to solve this problem was evaluated. As described in a previous section, hyperconjugation donating/ accepting abilities are often related to electronegativities, thus offering a fast and efficient way to compute the corresponding V1−3 values. The averaged V1−3 values for each hyperconjugation donor− acceptor pair, such as the average value for all the H−C → C− F* pairs appearing in 49 ethane derivatives in our set, were computed and compared to the electronegativity differences (Δχ) between donors and acceptors. The use of values averaged over donor−acceptor pairs (e.g., H−C−C−F) in different molecules reduced the dependence on the selected set of molecules and gave us a more general approximation to compute V1−3 from atomic properties. A previously reported study showed that the period of the atoms was critical for proper correlation between hyperconjugation and electronegativity.62 In agreement with this work, we observed

While V1−3 terms were traditionally perceived as tools to model the symmetry of the energy associated with torsions, we were also interested in what these V1−3 values represented in the context of hyperconjugation (Figure 6). Most of the torsional parameters in GAFF only have one such term (either V1, V2, or V3)12 and, based on our analysis, less than 10% of Parm@Frosst torsional parameters have more than one nonzero V term. A nonzero value provides the amplitude of the energy barriers with the corresponding periodicity. In contrast, the V1/V2/V3 hyperconjugation parameters derived herein were all nonzero (i.e., V1 = 2.17, V2 = −2.43, and V3 = 0.37 for the CH → CH* in ethane, Table 2). Hyperconjugation is defined by the electron donation occurring from bonding σ or nonbonding p orbitals to neighboring antibonding σ* orbitals. It is modulated by the degree of orbital overlap upon bond rotation. While the value of V1 controls the preference for trans or cis conformation, the V2 term ensures that hyperconjugation is minimal when the donor and acceptor bonds are perpendicular and maximal when parallel. We believe that V3 is a correcting factor to model the overall shape of the torsional energy profile more accurately. Thus, investigating hyperconjugation explicitly unveils the principles guiding the torsional behaviors rather than just providing the “apparent” symmetry. While the V1 and V2 energies would cancel out if all the nine torsions are identical and separated by 120° (ethane has nine equally contributing H−C−C−H torsions), their sum would not be zero for any G


Article


the chemical environment: we included the convergence point into the data set as an extreme condition and recomputed the linear regression for this correlation (Figure 10). The advantage

significant correlations between these averaged electronegativity differences and V1−3 values when considering element periodicity (Figure 9). This shift from one period to the next

Figure 10. Development of dynamic V2 parameters based on the sum of electronegativity of neighboring atoms including the convergence point. PC → CN* hyperconjugation, with only four data points, is also shown as an example. (See the SI for detailed equations of the regression.)

of this refinement was the reduction of the data set dependency when defining the V2neighboring ∑χ correlation: if only a few molecules contain the desired type of hyperconjugation, then the correlation observed for it might be misleading. The strong influence from neighboring atoms had only been observed on V2 and not on V1 and V3; besides, the V2 term primarily controls the strength (amplitude in the profile) of hyperconjugation, thus we applied this refinement method only to the V2 term while the above-mentioned relationship with electronegativity differences with the V1 and V3 terms were deemed acceptable. This advanced method considering the V2 variation is referred to as H-TEQ1.1. Influence of the Central Bond Atoms. While the C−C central bond was not polarized, the molecules featuring a C− N(+) bond were highly polarized and this polarization was expected to modulate the hyperconjugation profile. As shown in Table 2, with N(+) at the acceptor side, the V1−3 values increased, whereas with N(+) at the donor site, they dropped, indicating a stronger or weaker hyperconjugation contribution to the torsional energy, respectively. More interestingly, we observed a puzzling small negative V1 value for CH3/H−N(+)− C−F molecules. As V1 provides the cis−trans preference: this small negative value indicated a slight preference for the cis conformation previously computed (Table 1). Our first thought was to investigate the orbital overlap for these conformers as we did before. During the course of our investigations, we found that the computed NBO overlap correlated with the hyperconjugation energy with C−C-containing derivatives as previously reported and that this approach also predicted no hyperconjugation when the two orbitals were perpendicular (Table 1, Figure 6). However, it was not predicting the slight preference for the cis over the trans conformation in H−N(+)− C−F, the computed orbital overlaps being inconsistent with the hyperconjugation stabilizing energy of the trans- and cisconformers for the C−N(+) molecules. In order to explain this phenomenon, we reverted back to chemical principles. In the HSAB theory developed by Pearson, electronegativity and hardness are used to explain the preference of hard and soft acids for hard and soft bases,

Figure 9. Hyperconjugation V1 and V2 as a function of electronegativity difference with each data point representing the average V1−2 value. See complete data for V1, V2, and V3 in the Supporting Information (SI).

implicitly includes the roles of the size and softness of the elements on hyperconjugation. Thus, these trend lines can be used to approximate V1−3 values from simple electronegativity differences. We refer to the method which approximates hyperconjugation V1−3 terms only based on the Δχ(acceptor/donor) as H-TEQ1.0 (Hyperconjugation for Torsional Energy Quantification). Equations were developed based on these correlations; the V3 terms did not vary much within the same donor−acceptor combination, thus a constant value was applied. Developing Dynamic V2 TermsIncluding Neighboring Effects. While V1 and V3 values were not affected much by the chemical environment, V2 showed large variances among molecules. For example, CH → CF* is much weaker in pentafluoroethane than in monofluoroethane, with a smaller V2 but fairly similar V1 and V3 values (Figure 7). We found that in a highly electronegative environment, the V2 values were smaller and showed good correlation when plotted against the electronegativity sum (∑χ) of the neighboring atoms. Interestingly, all the trends, regardless of donor and acceptor period, converged at approximately the same point with high ∑χ (20.2) and low V2 value (−1.1). This data indicated that, once the neighboring atoms’ electronegativities were extremely high (hypothetically), electron donation from these specific bonding to antibonding orbitals was diminished and the hyperconjugation was minimized. This provided a way to develop more precise parameters for the V2 term in response to H


Article


hyperconjugation amplitude is controlled by the σ−σ* energy gap, which, in turn, is controlled by the electronegativity of the atoms making up the bond (i.e., correlates with the sum of their electronegativities) (Table 1). However, the polarization of the central bond may have the same (or opposite) direction as the σ → σ* hyperconjugation and may match or compete with the hyperconjugation donation (Figure 11c). Thus, we proposed to compute V2 with consideration that hyperconjugation is a combination of donation from bonding to antibonding orbitals as well as through central atoms (eq 4).

respectively.65 We thought to evaluate the impact of hardness on hyperconjugation energy. A method that considers both these properties is the electronegativity equalization method used to assign partial atomic charges.69 This method was applied to the molecules in Table 1, and while all of the trends can be somewhat predicted from this data, the difference in hyperconjugation between the cis and trans conformations of the fourth and seventh molecule in Table 1 could not be explained. We then switched to other electronic effects which could affect the cis/trans preference. We hypothesized that the central polarization strongly affects this preference through dipole−dipole interactions; while this effect is minimal in functionalized C−C bonds, it is maximized in H−N(+)−C−F (Figure 11a and b). This effect can be directly incorporated into V1.

Δχ = abs((χcenter1 + χterminal1 ) − (χcenter2 + χterminal2 )) − abs(χcenter1 − χcenter2 )

With this refinement, the V2 values of the nonpolarized ethane derivatives and polarized ammonium derivatives were aligned (H-TEQ1.3, Figure S4). The method H-TEQ1.4 is a combination of both refinement given to V1 (as in H-TEQ1.2) and V2 (as in H-TEQ1.3). Developing V1−3 on the Fly. With all the developed methods discussed above, unique parameters for each hyperconjugation pair in any molecule could be derived. Figure 12 offers a summarized sample of the full derivation of torsion parameters on the fly. The concept of group electronegativity was applied to consider functional groups70 (see the SI, section 5). Regardless of the molecule, using H-TEQ1.0, the V1−3 terms of each hyperconjugation pair (i.e., H−C → C−F*) were calculated using the linear fit shown in Figure 9 corresponding to the Δχ (donor−acceptor) and period. In H-TEQ1.1, V1 and V3 were computed in the same way as in H-TEQ1.0, but the V2 term was based on the correlation with neighboring Σχ in response to the element type of donor and acceptor as described in Figure 10. The central bond effects were added to V1 (H-TEQ1.2), V2 (H-TEQ1.3), and both V1 and V2 (HTEQ1.4), within the same frame of hyperconjugation− electronegativity principles. This strategy is quite simple, with a few rules (correlation line equations), and does not require

Figure 11. Polarization: (a) weak dipole−dipole interaction with ethane derivatives; (b) strong dipole−dipole interaction with ammonium derivatives; (c) central polarized bond.

Thus, rather than computing the Δχ between the terminal atoms of the hyperconjugation donor and acceptor, Δχ can be obtained from eq 3 for V1 terms and would consider the central bond as well. This method is referred to as H-TEQ1.2. With this method, the computed values of V1 were closer to those derived from NBO data (Table S6). Δχ = (χcenter1 − χterminal1 ) − (χcenter2 − χterminal2 )

(4)

(3)

The values of V2 could also been refined considering the central atoms’ electronegativities. As discussed above, the

Figure 12. Example of hyperconjugation parameters derived on the fly with H-TEQ1.0 and H-TEQ1.1 derivation based on the rules described in Figures 9 and 10. I


Article

Journal of Chemical Information and Modeling and is not limited toany atom types. Thus, it is expected to be transferable to a large set of molecules. With this method, we computed simple properties (i.e., electronegativity, atomic charges) for each molecule and were then able to derive the hyperconjugation, steric and electrostatic parameters on the fly as opposed to training on large sets in advance.

Table 3. Average RMSD of Energy over the Torsional Profile of FF vs MP2/6-311+G** (kcal/mol) for 337 Small Molecule Sets

■

method

RESULTS, DISCUSSION, AND APPLICATION Results and Discussion. Our approach, based on chemical principles, is compatible with all FFs which evaluate torsional energies as the sum of vdW, electrostatics, and additional torsion parameters, the latter being replaced by our hyperconjugation energy term. The performance of the different versions of our method was examined over three widely used FFs: GAFF as implemented in AMBER 11,71 Merck Molecular Force Field (MMFF94s72,73 as implemented in MOE74), and Parm@Frosst as implemented in MOE. Torsional parameters of given torsions were computed using H-TEQ and the rest of the MM terms and parameters were kept unchanged. At this stage, we investigated whether our hypothesis (“can hyperconjugation energies derived from electronegativities replace trained torsion parameters in force fields?”) was valid. We are aware that our V1−3 parameters were derived from NBO data itself derived from MP2/6-31+G** calculations while GAFF, MMFF94s, and Parm@Frosst parameters were derived using other methods (GAFF12 and MMFF94s75 were trained using MP2 but different basis sets). This may represent a bias as our accuracy is measured using MP2-derived profiles. However, at this proof-of-principle stage, optimal accuracy was not our primary objective and the reported accuracies were not corrected for this potential bias. H-TEQ1.4 improved the accuracy (when compared to MP2 energy profiles) of GAFF and Parm@Frosst FFs and retained the accuracy of MMFF94s (Table 3). More specifically, we evaluated the contribution of the torsional energy to the entire torsion energy profiles (which also includes vdW and electrostatic contributions) and found that H-TEQ improved it by 13% for GAFF and 46% for Parm@Frosst over the 337 small molecules. It should be noted that C−C and C−N rotatable bonds are very common in organic chemistry and are extensively parametrized in FFs. Despite the large training of the original FFs, many profiles were significantly improved by our approach and most notably, hexafluoroethane (reduction of the deviation between FF and MP2 by −1.25, −4.19, −1.25 kcal/mol for GAFF, MMFF94s, and Parm@Frosst, respectively). For several molecules, H-TEQ achieved better results across the board: for fluoroethanol, HTEQ improves GAFF and Parm@Frosst (RMSD lower by 0.93 kcal/mol in both cases) while deteriorating the accuracy of MMFF94s (RMSD increased by 0.46 kcal/mol). These changes in accuracy were likely due to the different training sets used in FFs development, hence covering different chemical spaces. Three out of 337 small molecules demonstrated a drop in accuracy greater than 0.5 kcal/mol and 16 in 337, a drop between 0.2 to 0.5 kcal/mol for all three FFs (H-TEQ1.4). Considering the self-consistency of FFs, where parameters are balanced with respect to each other, replacing one term (torsion) and keeping the others (sterics and electrostatics) unchanged as practiced here may lead to loss in accuracy. Nevertheless, H-TEQ showed overall good accuracy and demonstrated that the use of explicit hyperconjugation was able to lead to an accurate atom type-free torsion energy term. When applying MM modeling methods, it is critical to identify global and local minima. In order to focus on the most

GAFF referencea NBO corrected GAFFb H-TEQ1.0b H-TEQ1.1b H-TEQ1.2b H-TEQ1.3b H-TEQ1.4b MMFF94s reference H-TEQ1.0b H-TEQ1.1b H-TEQ1.2b H-TEQ1.3b H-TEQ1.4b Parm@ Frosst reference H-TEQ1.0b H-TEQ1.1b H-TEQ1.2b H-TEQ1.3b H-TEQ1.4b

small molecule set (337 molecules)

ethane derivatives (261 molecules)

ammonium derivatives (76 molecules)

1.82 2.53

0.90 1.44

1.62 2.29 1.49

1.36 1.77

1.62 1.63 1.62 1.62 1.62 1.31 2.00 1.35 1.35 1.35 1.35 1.35 1.62

1.25 1.26 1.26 1.22 1.23 1.05 1.52 0.95 0.94 0.98 0.85 0.89 1.66

1.53 1.54 1.53 1.53 1.53 1.25 1.90 1.26 1.26 1.26 1.24 1.24 1.63

1.37 1.38 1.37 1.37 1.37 1.76 1.68 1.47 1.48 1.47 1.46 1.46 1.52

2.20 1.39 1.40 1.39 1.39 1.39

1.90 1.43 1.40 1.45 1.33 1.36

2.13 1.40 1.40 1.41 1.38 1.39

1.84 1.46 1.46 1.46 1.46 1.45

drug fragments (59 cases)

a

The deviation of FFs from MP2 without the torsional term is set as a reference. bHyperconjugation term scaled by 0.6.

prevalent conformations, a cutoff of 4.0 kcal/mol from the global minimum in both MP2 and FFs was selected as a second evaluation (Table S7). With this second criterion, H-TEQ holds and/or improves the accuracy over the three FFs and the trends were somewhat similar to what is computed when considering the complete profiles. As summarized in Table 3, we have compared the results obtained with all five versions of H-TEQ to GAFF, MMFF94s, and Parm@Frosst. Overall, all the models displayed an improvement when implemented in GAFF and Parm@Frosst with similar accuracies across the versions. As most of the set did not require the modifications implemented in version 1.1 and above, the overall accuracy was not significantly different across the versions. However, the effect of the neighboring atoms and bond polarization should be taken into account in cases where the electronegativity differences are high. This was demonstrated by 1,2-difluoro-1,2dibromoethane, which features atoms with high electronegativity and the significantly better accuracy of H-TEQ1.4 for C−N(+)-containing molecules. Validation. In order to apply this method explicitly describing the role of hyperconjugation in torsion energy profiles, 59 existing, structurally diverse, drug molecules/ fragments were selected, including a variety of functional groups that were not found in our previous small molecule set. Among those, 11 molecules are built around a C−N(+) rotatable bond. When implementing our method into MMFF94s, we observed an overall slight improvement (RMSD lower by 0.3 kcal/mol) (Table 3). More specifically, 73% (43/59) of the molecules demonstrated an improved accuracy when applying J


Article

Journal of Chemical Information and Modeling H-TEQ1.4 to MMFF94s. Interestingly, the overall accuracy was even better when the torsional parameters of MMFF94s were not applied (as the reference in Table 3). A close look at this unexpected result revealed that this drop of accuracy was due to the use of generic (poorly transferable) torsion parameters for unparameterized torsions. For example, some of the torsion types in N-ethylphenylalanine and 1-(5-hydroxyhexyl)-3,7dimethyl-3,7-dihydro-1H-purine-2,6-dione (Figure 13) were

Figure 13. Example of missing parameters leading to poor torsion energy modeling. MMFF94s atom types are labeled in blue.

not defined in MMFF94s and lead to a low accuracy than not using the torsional parameter. As expected, H-TEQ did not suffer from this lack of transferability and demonstrated robustness on a diverse set of molecules. A slight accuracy improvement was obtained with Parm@ Frosst and a small accuracy drop with GAFF (Table 3). More than half of the molecule torsion profiles have improved when H-TEQ was implemented into GAFF and Parm@Frosst (Figure S6). In addition, we observed an improvement of the prediction of the minima location: Ifosfamide, presented in Figure 14a, showed a better fit to the MP2 energy profile with accurate location of the minima as well as the relative energy if applying H-TEQ to MMFF94s, while MMFF94s alone was not able to predict the gauche effect. Similarly, although the overall RMSD did not change much when applying H-TEQ to Parm@ Frosst for a muraglitazar fragment (a peroxisome proliferatoractivated receptor agonist, Figure 14b), and to GAFF for 2butylimidazole (Figure 14c), minima locations were correctly predicted and the relative energies between minima were wellaligned with MP2 energy. This further demonstrated that our method was transferable to real drug molecules whether or not the fragments were involved in the data set used for the method development while no further correction of the weights for vdW and electrostatic contribution was carried out.

Figure 14. Examples of improving the minima prediction: (a) for MMFF94s, RMSDH‑TEQ1.4 = 1.11 kcal/mol vs RMSDMMFF94s = 1.86 kcal/mol; (b) for Parm@Frosst, RMSDH‑TEQ1.4 = 0.64 kcal/mol vs RMSDParm@Frosst = 0.90 kcal/mol; (c) for GAFF, RMSDH‑TEQ1.4 = 0.39 kcal/mol vs RMSDParm@Frosst = 0.63 kcal/mol.

■

CONCLUSION With these proof-of-principle investigations, we have shown that hyperconjugation is a major contributor to the torsional energy, can be quantified efficiently and, along with vdW and electrostatic parameters, can accurately yield torsional energy profiles. Correlating hyperconjugation parameters to electronegativity makes it computationally accessible and enables MM application. Conformational accuracy has been demonstrated not only on a set of small molecules, where the relations between hyperconjugation and electronegativity were built, but also on real drug molecules/fragments. More importantly, our approach did not require precomputation of torsional parameters for any new torsion (i.e., not covered by the available parameter sets) but generated them on the fly based

on the electronegativity of atoms and functional groups. Furthermore, this atom-type free method also precludes any potential atom type assignment difficulties, another challenge for atom type-based FFs. We expect that this approach, will be a more transferable and reliable alternative to FFs once applied to more diverse molecules, including potential drug candidates. Evaluation over three commonly used FFs indicated that our method was robust and independent of any specific FF, although improvements will be required. It should be noted that the H-TEQ method has been implemented into three selfconsistent FFs without retraining of the other terms and the scaling factors used for the 1,4-vdW and 1,4-electrostatic terms. This work is a first step toward a fully atom-type free MM method that could be used to compute potential energy K



■

surfaces (i.e., conformational preferences) of organic drug molecules.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.6b00012. Computed correlations for the different versions of HTEQ, group electronegativity rules, accuracy of the different versions using two criteria, and evaluation of a scaling factor for NBO-derived hyperconjugation energy (PDF) Small molecule set (TXT) Drug molecule/fragment set (TXT)

■

REFERENCES

(1) Fitzgerald, S. H.; Sabat, M.; Geysen, H. M. Diversity Space and Its Application to Library Selection and Design. J. Chem. Inf. Model. 2006, 46, 1588−1597. (2) Krier, M.; Bret, G.; Rognan, D. Assessing the Scaffold Diversity of Screening Libraries. J. Chem. Inf. Model. 2006, 46, 512−524. (3) MacLellan, P.; Nelson, A. A Conceptual Framework for Analysing and Planning Synthetic Approaches to Diverse Lead-Like Scaffolds. Chem. Commun. 2013, 49, 2383−2393. (4) Burke, M. D.; Schreiber, S. L. A Planning Strategy for DiversityOriented Synthesis. Angew. Chem., Int. Ed. 2004, 43, 46−58. (5) Kim, J.; Kim, H.; Park, S. B. Privileged Structures: Efficient Chemical ″Navigators″ toward Unexplored Biologically Relevant Chemical Spaces. J. Am. Chem. Soc. 2014, 136, 14629−14638. (6) Sen, S.; Prabhu, G.; Bathula, C.; Hati, S. Diversity-Oriented Asymmetric Synthesis. Synthesis 2014, 46, 2099−2121. (7) Dow, M.; Fisher, M.; James, T.; Marchetti, F.; Nelson, A. Towards the Systematic Exploration of Chemical Space. Org. Biomol. Chem. 2012, 10, 17−28. (8) Langdon, S. R.; Brown, N.; Blagg, J. Scaffold Diversity of Exemplified Medicinal Chemistry Space. J. Chem. Inf. Model. 2011, 51, 2174−2185. (9) Polanski, J.; Kurczyk, A.; Bak, A.; Musiol, R. Privileged Structures - Dream or Reality: Preferential Organization of Azanaphthalene Scaffold. Curr. Med. Chem. 2012, 19, 1921−1945. (10) Moitessier, N.; Englebienne, P.; Lee, D.; Lawandi, J.; Corbeil, C. R. Towards the Development of Universal, Fast and Highly Accurate Docking/Scoring Methods: A Long Way to Go. Br. J. Pharmacol. 2008, 153, S7−S26. (11) Durrant, J. D.; McCammon, J. A. Molecular Dynamics Simulations and Drug Discovery. BMC Biol. 2011, 9, 71−71. (12) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (13) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (14) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. An All Atom Force-Field for Simulations of Proteins and Nucleic-Acids. J. Comput. Chem. 1986, 7, 230−252. (15) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Alagona, G.; Profeta, S., Jr; Weiner, P.; Ghio, C. New Force Field for Molecular Mechanical Simulation of Nucleic Acids and Proteins. J. Am. Chem. Soc. 1984, 106, 765−784. (16) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. (17) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, A. D., Jr Charmm General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671−690. (18) Halgren, T. A.; Nachbar, R. B. Merck Molecular Force Field. IV. Conformational Energies and Geometries for Mmff94. J. Comput. Chem. 1996, 17, 587−615. (19) Halgren, T. A. Merck Molecular Force Field. V. Extension of MMFF94 Using Experimental Data, Additional Computational Data, and Empirical Rules. J. Comput. Chem. 1996, 17, 616−641. (20) Halgren, T. A. Merck Molecular Force Field. I. Basis, Form, Scope, Parameterization, and Performance of MMFF94. J. Comput. Chem. 1996, 17, 490−519.

EXPERIMENTAL SECTION Hyperconjugation Energy Calculation from NBO. Molecules initially underwent a full optimization followed by freezing the desired torsion into certain degrees (from −180 to 180 with 10° increment) and reoptimizing in MP2/6-311+G** level using software GAMESS-US.76,77 Basis sets not available in the GAMESS-US package were downloaded from basis set exchange (https://bse.pnl.gov/bse/portal) based on previous studies.78 Reoptimization provided the right geometry and orbital alignment in respond to the defined torsion. NBO calculation was then performed in NBO 6.064 using these conformations with the second order perturbation analysis applied to obtain the hyperconjugation energy in the same level of theory and basis set. Donations from the bonding orbital as well as the available lone pairs of the terminal atoms, for example three lone pairs of fluorine in F−C → C−F* were combined together for this hyperconjugation. MM Calculation. The AMBER11 package was used to perform the calculation with GAFF while the MOE platform was used to compute the energy profiles with MMFF94s and Parm@Frosst. GAFF atom types were assigned using the tLeap routine and the partial charges were assigned using the AM1BCC method on the global minimum structures. These same charges were used for the other conformations of the same molecules. The GAFF-derived potential energy is computed using the sander routine. MMFF94s and Parm@Frosst atom types and atom partial charges were assigned. Partial charges were added prior to energy calculation. A program was written to compute the torsional energy contributions for only the nine torsions related to the central bond as measured by the FFs. These energies were replaced by H-TEQ (1.0−1.4), and all the other energy contributions were retained when evaluate the FF accuracy.

■

Article

AUTHOR INFORMATION

Corresponding Author

*Tel.: [email protected]. Tel.: 514-398-8543. Fax: 514-398-3797. Funding

We thank NSERC (through a CRD grant, CRDPJ 469677) for funding. J.P. thanks FQRNT for a scholarship, and A.T. thanks the CIHR-funded Drug Discovery and Development Training program. Calcul Québec and Compute Canada are acknowledged for generous CPU allocations. Notes

The authors declare no competing financial interest. L


Article


(40) Kirkpatrick, P.; Ellis, C. Chemical Space. Nature 2004, 432, 823−823. (41) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. Dreiding: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897−8909. (42) Rappé, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. Uff, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (43) Gerber, P. R.; Müller, K. Mab, a Generally Applicable Molecular Force Field for Structure Modelling in Medicinal Chemistry. J. Comput.-Aided Mol. Des. 1995, 9, 251−268. (44) Gerber, P. R. Charge Distribution from a Simple Molecular Orbital Type Calculation and Non-Bonding Interaction Terms in the Force Field Mab. J. Comput.-Aided Mol. Des. 1998, 12, 37−51. (45) Shi, S.; Yan, L.; Yang, Y.; Fisher-Shaulsky, J.; Thacher, T. An Extensible and Systematic Force Field, Esff, for Molecular Modeling of Organic, Inorganic, and Organometallic Systems. J. Comput. Chem. 2003, 24, 1059−1076. (46) Gundertofte, K.; Liljefors, T.; Norrby, P. O.; Pettersson, I. A Comparison of Conformational Energies Calculated by Several Molecular Mechanics Methods. J. Comput. Chem. 1996, 17, 429−449. (47) Kuhn, B.; Gerber, P.; Schulz-Gasch, T.; Stahl, M. Validation and Use of the MM-PBSA Approach for Drug Discovery. J. Med. Chem. 2005, 48, 4040−4048. (48) Cramer, C. J. Essentials of Computational Chemistry: Theories and Models, 2nd ed.; John Wiley & Sons: New York, NY, 2004. (49) Anh, N. T.; Eisenstein, O.; Lefour, J. M.; Trân Huu Dâu, M. E. Orbital Factors and Asymmetric Induction. J. Am. Chem. Soc. 1973, 95, 6146−6147. (50) Anh, N. Regio- and Stereo-Selectivities in Some Nucleophilic Reactions. In Organic Chemistry Syntheses and Reactivity; Topics in Current Chemistry; Springer: Berlin Heidelberg, 1980; Vol. 88; pp 145−162. (51) Cieplak, A. S. Stereochemistry of Nucleophilic Addition to Cyclohexanone. The Importance of Two-Electron Stabilizing Interactions. J. Am. Chem. Soc. 1981, 103, 4540−4552. (52) Glendening, E. D.; Landis, C. R.; Weinhold, F. NBO 6.0: Natural Bond Orbital Analysis Program. J. Comput. Chem. 2013, 34, 1429−1437. (53) Knizia, G.; Klein, J. E. M. N. Electron Flow in Reaction MechanismsRevealed from First Principles. Angew. Chem., Int. Ed. 2015, 54, 5518−5522. (54) Corbeil, C. R.; Thielges, S.; Schwartzentruber, J. A.; Moitessier, N. Toward a Computational Tool Predicting the Stereochemical Outcome of Asymmetric Reactions: Development and Application of a Rapid and Accurate Program Based on Organic Principles. Angew. Chem., Int. Ed. 2008, 47, 2635−2638. (55) Weill, N.; Corbeil, C. R.; De Schutter, J. W.; Moitessier, N. Toward a Computational Tool Predicting the Stereochemical Outcome of Asymmetric Reactions: Development of the Molecular Mechanics-Based Program Ace and Application to Asymmetric Epoxidation Reactions. J. Comput. Chem. 2011, 32, 2878−2889. (56) Pottel, J.; Therrien, E.; Gleason, J. L.; Moitessier, N. Docking Ligands into Flexible and Solvated Macromolecules. 6. Development and Application to the Docking of Hdacs and Other Zinc Metalloenzymes Inhibitors. J. Chem. Inf. Model. 2014, 54, 254−265. (57) De Cesco, S.; Deslandes, S.; Therrien, E.; Levan, D.; Cueto, M.; Schmidt, R.; Cantin, L. D.; Mittermaier, A.; Juillerat-Jeanneret, L.; Moitessier, N. Virtual Screening and Computational Optimization for the Discovery of Covalent Prolyl Oligopeptidase Inhibitors with Activity in Human Cells. J. Med. Chem. 2012, 55, 6306−6315. (58) Pearson, R. G. Hard and Soft Acids and Bases. J. Am. Chem. Soc. 1963, 85, 3533−3539. (59) Pophristic, V.; Goodman, L. Hyperconjugation Not Steric Repulsion Leads to the Staggered Structure of Ethane. Nature 2001, 411, 565−568. (60) Mo, Y.; Gao, J. Theoretical Analysis of the Rotational Barrier of Ethane. Acc. Chem. Res. 2007, 40, 113−119.

(21) Halgren, T. A. Merck Molecular Force Field. II. MMFF94 van der Waals and Electrostatic Parameters for Intermolecular Interactions. J. Comput. Chem. 1996, 17, 520−552. (22) Halgren, T. A. Merck Molecular Force Field. III. Molecular Geometries and Vibrational Frequencies for Mmff94. J. Comput. Chem. 1996, 17, 553−586. (23) Daura, X.; Mark, A. E.; Van Gunsteren, W. F. Parametrization of Aliphatic Chn United Atoms of Gromos96 Force Field. J. Comput. Chem. 1998, 19, 535−547. (24) Oostenbrink, C.; Soares, T. A.; Van Der Vegt, N. F. A.; Van Gunsteren, W. F. Validation of the 53a6 Gromos Force Field. Eur. Biophys. J. 2005, 34, 273−284. (25) Soares, T. A.; Daura, X.; Oostenbrink, C.; Smith, L. J.; van Gunsteren, W. F. Validation of the Gromos Force-Field Parameter Set 45a3 against Nuclear Magnetic Resonance Data of Hen Egg Lysozyme. J. Biomol. NMR 2004, 30, 407−422. (26) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The Gromos Force-Field Parameter Sets 53a5 and 53a6. J. Comput. Chem. 2004, 25, 1656−1676. (27) Chandrasekhar, I.; Kastenholz, M.; Lins, R. D.; Oostenbrink, C.; Schuler, L. D.; Tieleman, D. P.; Van Gunsteren, W. F. A Consistent Potential Energy Parameter Set for Lipids: Dipalmitoylphosphatidylcholine as a Benchmark of the Gromos96 45a3 Force Field. Eur. Biophys. J. 2003, 32, 67−77. (28) Eichenberger, A. P.; Allison, J. R.; Dolenc, J.; Geerke, D. P.; Horta, B. A. C.; Meier, K.; Oostenbrink, C.; Schmid, N.; Steiner, D.; Wang, D.; Van Gunsteren, W. F. Gromos++ Software for the Analysis of Biomolecular Simulation Trajectories. J. Chem. Theory Comput. 2011, 7, 3379−3390. (29) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the Opls All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (30) Kaminski, G.; Jorgensen, W. L. Performance of the AMBER94, MMFF94, and OPLS-aa Force Fields for Modeling Organic Liquids. J. Phys. Chem. 1996, 100, 18010−18013. (31) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and Reparametrization of the Opls-Aa Force Field for Proteins Via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105, 6474−6487. (32) Avakyan, V. G.; Nazarov, V. B.; Alfimov, M. V.; Bagaturyants, A. A.; Voronezheva, N. I. The Role of Intra- and Intermolecular Hydrogen Bonds in the Formation of B-Cyclodextrin Head-to-Head and Head-to-Tail Dimers. The Results of Ab Initio and Semiempirical Quantum-Chemical Calculations. Russ. Chem. Bull. 2001, 50, 206− 216. (33) Woods, R. J.; Dwek, R. A.; Edge, C. J.; Fraser-Reid, B. Molecular Mechanical and Molecular Dynamical Simulations of Glycoproteins and Oligosaccharides. 1. Glycam_93 Parameter Development. J. Phys. Chem. 1995, 99, 3832−3846. (34) Kirschner, K. N.; Woods, R. J. Solvent Interactions Determine Carbohydrate Conformation. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10541−10545. (35) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. L. OPLS All-Atom Force Field for Carbohydrates. J. Comput. Chem. 1997, 18, 1955−1970. (36) Momany, F. A.; Willett, J. L. Computational Studies on Carbohydrates: In Vacuo Studies Using a Revised AMBER Force Field, Amb99c, Designed for A-(1→4) Linkages. Carbohydr. Res. 2000, 326, 194−209. (37) Ott, K.-H.; Meyer, B. Parametrization of Gromos Force Field for Oligosaccharides and Assessment of Efficiency of Molecular Dynamics Simulations. J. Comput. Chem. 1996, 17, 1068−1084. (38) Mayne, C. G.; Saam, J.; Schulten, K.; Tajkhorshid, E.; Gumbart, J. C. Rapid Parameterization of Small Molecules Using the Force Field Toolkit. J. Comput. Chem. 2013, 34, 2757−2770. (39) http://www.ccl.net/cca/data/parm_at_Frosst/, (accessed March 3, 2016). M


Article

Journal of Chemical Information and Modeling (61) Allinger, N. L.; Chen, K.; Lii, J.-H. An Improved Force Field (MM4) for Saturated Hydrocarbons. J. Comput. Chem. 1996, 17, 642− 668. (62) Alabugin, I. V.; Zeidan, T. A. Stereoelectronic Effects and General Trends in Hyperconjugative Acceptor Ability of Σ Bonds. J. Am. Chem. Soc. 2002, 124, 3175−3185. (63) Alabugin, I. V.; Manoharan, M. Effect of Double-Hyperconjugation on the Apparent Donor Ability of Σ-Bonds: Insights from the Relative Stability of Δ-Substituted Cyclohexyl Cations. J. Org. Chem. 2004, 69, 9011−9024. (64) Glendening, E. D.; Badenhoop, J. K.; Reed, A. E.; Carpenter, J. E.; Bohmann, J. A.; Morales, C. M.; Landis, C. R.; Weinhold, F. NBO, 6.0 ed.; 2013. (65) Pearson, R. G. Absolute Electronegativity and Hardness Correlated with Molecular Orbital Theory. Proc. Natl. Acad. Sci. U. S. A. 1986, 83, 8440−8441. (66) Weinhold, F.; Landis, C. R. Valency and Bonding: A Natural Bond Orbital Donor−Acceptor Perspective; Cambridge University Press: Cambridge, 2005. (67) Monticelli, L.; Tieleman, D. P. Force Fields for Classical Molecular Dynamics. In Biomolecular Simulations; Monticelli, L., Salonen, E., Eds.; Methods in Molecular Biology; Humana Press: Clifton, NJ, 2013; Vol. 924; pp 197−213. (68) R: A Language and Environment for Statistical Computing; Vienna, Austria, 2008. (69) Mortier, W. J.; Ghosh, S. K.; Shankar, S. ElectronegativityEqualization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315−4320. (70) Campanelli, A. R.; Domenicano, A.; Ramondo, F.; Hargittai, I. Group Electronegativities from Benzene Ring Deformations: A Quantum Chemical Study. J. Phys. Chem. A 2004, 108, 4940−4948. (71) Case, D. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Zhang, W.; Merz, K. M.; Roberts, B.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Wong, K. F.; Paesani, F.; Vanicek, J.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 11; San Francisco, 2010. (72) Halgren, T. A. Mmff Vi. Mmff94s Option for Energy Minimization Studies. J. Comput. Chem. 1999, 20, 720−729. (73) Halgren, T. A. Mmff Vii. Characterization of MMFF94, MMFF94s, and Other Widely Available Force Fields for Conformational Energies and for Intermolecular-Interaction Energies and Geometries. J. Comput. Chem. 1999, 20, 730−748. (74) Molecular Operating Environment (MOE), 2013.08 ed.; Montreal, 2015. (75) Halgren, T. A.; Nachbar, R. B. Merck Molecular Force Field. IV. Conformational Energies and Geometries for Mmff94. J. Comput. Chem. 1996, 17, 587−615. (76) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (77) Gordon, M. S.; Schmidt, M. W. Chapter 41−Advances in Electronic Structure Theory: Gamess a Decade Later. In Theory and Applications of Computational Chemistry; Scuseria, C. E. D. F. S. K. E., Ed.; Elsevier: Amsterdam, 2005; pp 1167−1189. (78) Curtiss, L. A.; McGrath, M. P.; Blaudeau, J. P.; Davis, N. E.; Binning, R. C.; Radom, L. Extension of Gaussian-2 Theory to Molecules Containing Third-Row Atoms Ga-Kr. J. Chem. Phys. 1995, 103, 6104−6113.

N


Elucidating Hyperconjugation from Electronegativity to Predict Drug

Recommend Documents