Physical Mechanisms Governing Substituent Effects on Arene-Arene Interactions in a Protein Milieu.

Arene-arene interactions play important roles in protein-ligand complex formation. Here, we investigate the characteristics of arene-arene interactions between small organic molecules and aromatic amino acids in protein interiors. The study is based on X-ray crystallographic data and quantum mechanical calculations using the enzyme acetylcholinesterase and selected inhibitory ligands as a model system. It is shown that the arene substituents of the inhibitors dictate the strength of the interaction and the geometry of the resulting complexes. Importantly, the calculated interaction energies correlate well with the measured inhibitor potency. Non-hydrogen substituents strengthened all interaction types in the protein milieu, in keeping with results for benzene dimer model systems. The interaction energies were dispersion-dominated, but substituents that induced local dipole moments increased the electrostatic contribution and thus yielded more strongly bound complexes. These findings provide fundamental insights into the physical mechanisms governing arene-arene interactions in the protein milieu and thus into molecular recognition between proteins and small molecules.

Introduction

Noncovalent molecular interactions involving two or more aromatic hydrocarbon moieties, so-called arene–arene interactions, are important in diverse chemical and biological processes: they help stabilize the 3D structure of macromolecules and play vital roles in molecular recognition, supramolecular chemistry, organic solar cells, and catalyst design.[1−9] The arene–arene interactions in protein–inhibitor systems where aromatic amino acid side chains (Phe, Tyr, Trp, and His) interact with aromatic moieties of an inhibitor are of particular interest to medicinal chemists working in drug discovery. A deeper understanding of such interactions could provide valuable guidance when optimizing inhibitors into drug candidates. Here, we investigate the strengths and geometries of arene–arene interactions between inhibitors and aromatic residues within a protein interior milieu using the enzyme acetylcholinesterase (AChE). Conformations observed in crystal structures were studied by performing coupled-cluster singles, doubles, and noniterative triples correction (CCSD(T)), dispersion-corrected density functional theory (DFT), and symmetry-adapted perturbation theory (SAPT) calculations whose results were then compared to measurements of inhibitory activity.

Most studies on arene–arene interactions have focused on small model systems consisting of benzene dimers with geometries including parallel stacked “sandwich” dimers (aligned “face-to-face”), parallel-displaced dimers (offset “face-to-face”), and T-shaped “edge-to-face” dimers. Highly accurate CCSD(T)[10−12] calculations have shown that the T-shaped and parallel-displaced dimers are more energetically favorable than the sandwich dimer,[13−16] with stabilization energies of 2.5–2.8 kcal mol–1,[13−16] compared to 1.5–1.7 kcal mol–1 for the latter.[13−15]

Substituents on the aromatic rings influence both the geometries and the interaction energies of interacting dimers and heterodimers.[15,17−21] Two main physical mechanisms have been proposed to explain substituents’ effects on the interaction energies of the sandwich and parallel-displaced benzene dimers. One suggests that the substituents affect the π-electron density of the aromatic ring and thus influence the arene–arene electrostatic interactions.[22,23] The other states that the effects result from local dipole moments induced by the substituents, which form local direct interactions with the other arene.[19,24] This model is supported by results showing that monosubstitution of sandwiched and parallel displaced benzene dimers increased the interaction strengths regardless of the electron-donating or -withdrawing character of the substituent.[15,19,25] For example, Kim et al.[15] showed that the interaction strengths of parallel-displaced monosubstituted heterodimers decreased in the order NO2 > NH2 > CN > CH3 > OH > F > H, showing that the nonsubstituted system had weaker interactions than systems substituted with either electron-donating or electron-withdrawing groups. This was recently confirmed by experiments using a guest–host system designed to study parallel displaced arene–arene interactions, although in this system the interaction strength for the CH3-substituted complex was lower than that for the unsubstituted case.[26] Furthermore, Parrish and Sherrill dissected the details of substituent effects in sandwiched benzenes, revealing a more complex picture in which both hypotheses had some validity but local direct interactions generally dominated.[27]

The strengths of T-shaped benzene interactions depend on the electron-withdrawing or -donating nature of the substituent and on whether the substitution is on the arene C–H donor or on the acceptor benzene.[18,21,28] Electron-withdrawing substituents on the C–H donor reduce the electron density around the hydrogens at the interface between the two rings, resulting in a stronger interaction.[21,28,29] The interaction is also strengthened if the substituent on the H-donating ring is in close proximity to the accepting ring.[28]

Stacked benzenes are primarily stabilized by dispersion,[14,30] but electrostatics and induction must be considered to understand the effects of substituents.[20,28,31] A study of the contributions of the electrostatic energies in sandwiched arene systems with varying distances between the aryl rings revealed that at large intermolecular separations (5 Å < r < 7 Å) electron-donating substituents reduced the electrostatic term while electron-withdrawing substituents enhanced it. On the other hand, at smaller intermolecular separations (3 Å < r < 4 Å), all substituents increased the electrostatic interaction contributions.[32] The Hammett constant quantifies the electron-withdrawing or -donating effects of para/meta-substituents in terms of their effects on the acidity of benzoic acid and thus indirectly measures their inductive and mesomeric effects. A few studies have shown that the sum of Hammett constants (∑σ) can be used as a predictor for the interaction energy of benzene sandwich dimers with electron-withdrawing substituents (∑σ > 0),[31,33] although others have questioned its usefulness.[19,32,34] Recently, Wheeler et al.[35] developed models that could predict the strength of stacking interactions between heterocycles and amino acids based on the size and polarizability of the heterocycles together with electrostatics descriptors. This shows that these interactions can be estimated by using quite simple metrics.

Most computational studies on substituents’ effects on arene–arene interactions in small model systems have focused on systems where the rings are symmetrically aligned or have undergone full geometry optimization. However, the conformations observed in crystal structures represent low-energy conformational states relevant to biological systems. In biological systems, the geometries of arene interactions are controlled by factors that depend on surrounding structures and solvent molecules. For example, most stacking interactions in protein crystal structures have “suboptimal” geometries, with tilt angles above 30°, and 50% of complexes with tilt angles of 10°–50° have large (>3.5 Å) offsets.[36] The high abundance of such arene–arene geometries in proteins justifies a detailed study of the interactions in these complexes to support efforts in fields such as molecular recognition and drug discovery.[37−39] Here, we investigate arene–arene interactions in systems with geometries observed in the interior of the protein AChE by X-ray crystallography. AChE is an essential enzyme in the nervous system that terminates cell signaling by catalyzing acetylcholine hydrolysis. Its active site (Figure ) is a deep gorge lined with aromatic amino acid residues that can host aromatic and cationic inhibitors, providing an excellent system for studying interactions of these types. The gorge can be divided into the peripheral aromatic site (PAS) at the entrance and the catalytic site (CAS) in the interior. AChE inhibitors are used in the symptomatic treatment of conditions including myasthenia gravis and Alzheimer’s disease.[40,41] AChE is also the target of nerve agents; antidotes to these toxic substances reactivate AChE via noncovalent binding followed by a reaction to release the functional enzyme. We have previously demonstrated the formation of nonclassical hydrogen bonds between cationic inhibitors and aromatic amino acid residues using crystal structures and calculations.[42] Here, we focus on arene–arene interactions to clarify the physical mechanisms underpinning substituents’ effects on their geometries and interaction energies in protein milieus using AChE inhibitors (Table and Figure ) as a study case.

Figure 1

(a) AChE (PDB code 6TD2) shown in ribbon form with the inhibitor AL037 (green) and the active site gorge with the PAS in blue and the CAS in red. (b) Close-up of the aromatic moiety of AL037 (green) and the interacting amino acid side chains in the PAS. (c) Noncovalent interactions investigated for the benzylic inhibitors. (d) AChE (PDB code 4B84) in complex with AL040 (purple) and the interacting amino acid side chains in the PAS. (e) Noncovalent interactions investigated for the non-benzylic compounds.

Table 1

Inhibitors Included in the Study with X-ray Crystal Structure PDB Codes for Complexes with AChE and Inhibition Dataa

Inhibitor synthesis and IC50 data were published previously by Andersson et al.[45]

Determined for Mus musculus AChE; 95% confidence interval within parentheses based on at least four replicates.

Modeled from 4B7Z.

Methods

X-ray Crystallography

Wild-type AChE from Mus musculus was expressed in HEK293F cells and purified and crystallized as previously described.[43] Briefly, HEK293F cells expressing secreted AChE were grown in suspension using Freestyle 293 and Glutamax (Gibco) media containing 20 μg mL–1 Gentamicin (Gibco). The AChE-containing supernatant was centrifuged, and the protein was purified by using affinity and size exclusion chromatography. Protein crystallization was done by the hanging drop vapor diffusion method at a protein concentration of 10 mg mL–1 with a well solution containing 27–30% (w/v) PEG750MME and 0.1 M HEPES, pH 6.9–7.1. The AChE–AL037 complex was generated by soaking the compound into AChE crystals before flash freezing in liquid nitrogen as previously described.[44] X-ray data collection was performed at the MAX-lab synchrotron (Lund, Sweden) on beamline I911-5 using a MAR Research CCD detector. The data processing and refinement of AChE–AL037 was performed by using the same software, rejection criteria, and strategy as the previously reported crystal structures included in this study.[45] The intensity data were indexed and integrated by XDS[46] and scaled by using Scala in the CCP4 suite.[47] The structure was determined by difference Fourier methods with a modified apo structure of AChE (PDB code 1J06) as a starting model using Refmac.[48] Further crystallographic refinement and manual rebuilding were performed by using Phenix[49] and COOT.[50] The quality of the model was validated by using MolProbity (in Phenix[49]) and the wwPDB Validation Service. Data collection and refinement statistics are listed in Table S1 of the Supporting Information. The coordinates and structure factors have been deposited in the RCSB Protein Data Bank with accession code 6TD2.

Preparation of AChE–Inhibitor Complexes

X-ray crystal structures of AChE in complex with inhibitors (PDB codes 4B81, 4B80, 4B7Z, 4B84, and 4B83)[45] and the new structure obtained in this work (6TD2) with resolutions between 2.3 and 2.8 Å were prepared in Maestro[51] by adding hydrogen atoms, assigning bond orders, and defining disulfide bonds. The protonation states of amino acids and ligands were set as they would be at pH 7. The AChE–inhibitor complexes were energy minimized by using the MMFF94s[52,53] force field as implemented in MacroModel.[54] A constrained energy minimization was performed where hydrogens were allowed to move freely while all heavy atoms were restrained to a maximum deviation of 0.2 Å from their initial positions by applying a force of 100 kcal mol–1 Å–2. Wherever structures were truncated for force field or quantum mechanical calculations, hydrogens were added and minimized with the same method.

Quantum Mechanical Geometry Optimizations and Interaction Energy Calculations

Interaction energies for the compounds’ interactions with AChE residues were estimated by using a chemical cluster approach[55,56] in which only the enzyme’s core amino acid residues flanking the inhibitor were included. Generally, the core encompassed 400–423 atoms, and amino acid residue atoms were included if they belonged to a residue in direct contact with the inhibitor or were needed to control the movements of residues in contact with the inhibitor. The inhibitor and residues in direct contact with it were allowed to move during the geometry optimization (121–138 atoms), while the remaining atoms in the core were frozen at the coordinates obtained from the constrained energy minimization. The molecular coordinates are given in the Supporting Information. DFT geometry optimizations were performed by using BLYP-D3[57−59] with the 6-31G** basis set as implemented in Jaguar.[60] We have previously shown that this method offers a good trade-off between speed and accuracy when performing geometry optimizations of large systems.[42] The self-consistent field calculations were run at the ultrafine level using the ultrafine pseudospectral grid type. All optimizations were run until convergence using the direct inversion of the iterative subspace method,[61] with energy and root-mean-square density matrix change cutoffs set to 5 × 10–5 and 5 × 10–6 hartree, respectively.

Single-point energies were calculated by using three methods. Given the system size involved in this study, standard CCSD(T) calculations were not viable. We therefore instead applied domain-based local pair natural orbital coupled cluster theory with single, double, and perturbative triple excitations (DLPNO-CCSD(T)) as implemented in ORCA.[62,63] DLPNO-CCSD(T) is a linear scaling method and can thus be applied efficiently to large systems, and its accuracy is comparable to that of the “gold standard” method CCSD(T). PBE[64]-D3[59] was chosen based on a benchmarking study performed by us,[42] which showed that this functional yielded an error of only 0.58 kcal mol–1 (6%) when compared to benchmark values obtained at the CCSD(T)/CBS level. Additionally, since dispersion contributes significantly to the interaction energies of interest here, the use of dispersion correction was expected to increase accuracy.[65] We also assessed the accuracy of the implicitly dispersion-corrected MN15 method recently developed by Truhlar and co-workers,[66] which yielded a mean unsigned error of only 0.25 kcal mol–1 for 87 noncovalent data sets. Gas phase electronic interaction energies (ΔE) were calculated by using the def2-TZVP basis set according towhere E(complex), E(protein), and E(inhibitor) correspond to the electronic energy of the two-body complex (selected protein atoms plus selected inhibitor atoms), the selected protein atoms, and the selected inhibitor atoms, respectively, as shown in Figures and 4–6. Electrostatic potential maps were generated by using B3LYP-D3[57−59] with the 6-31G** basis set using Jaguar[60] and were analyzed at isovalues between 0.001 and 0.01 electrons (bohr3)−1.

Figure 2

Parallel-displaced stacking interactions between benzylic inhibitors and Tyr341 and non-benzylic inhibitors and Phe338. DLPNO–CCSD(T) interaction energies are shown in kcal mol–1 (PBE-D3 values are given in parentheses). The given distances are ring center-to-center, and the given angles are between rings. The complexes’ electrostatic potentials are shown at an isovalue of 0.005 electrons (bohr3)−1; red/yellow and blue/purple indicate electron-rich and -poor areas, respectively, on a scale from −50 to 90 kcal mol–1.

Figure 4

(a) T-shaped interactions between non-benzylic inhibitors and Tyr341 with DLPNO-CCSD(T) interaction energies in kcal mol–1 (PBE-D3 energies in parentheses), ring center distances (in Å and also shown using dashed lines), and the angles between ring planes (in degrees). The electrostatic potential surface of each complex is shown at an isovalue of 0.005 electrons (bohr3)−1; red/yellow and blue/purple indicate electron-rich and -poor areas, respectively, on a scale ranging from −50 to 90 kcal mol–1 (b) Interaction energy components calculated by using DFT-SAPT. (c) Contributions (in percent) of the attractive energy components to the interaction energy.

Figure 6

Hydrogen bonds between inhibitors and tyrosine with DLPNO-CSD(T) interaction energies in kcal mol–1 (PBE-D3 values in parentheses). Classical hydrogen bonds (H-bond) are indicated by black dashed lines; their lengths are given in Å, and their N–H···O bond angles are given in deg. Hydroxyl contacts are indicated by dashed blue lines, and the distance between the hydroxyl hydrogen and the ring center (O–H···π) is given in Å. The complexes’ electrostatic potentials are shown at an isovalue of 0.005 electrons (bohr3)−1; red/yellow and blue/purple indicate electron-rich and -poor areas, respectively, on a scale ranging from −50 to 90 kcal mol–1.

Interaction Energy Decomposition

To understand the fundamental nature of the arene–arene interactions, the SAPT method[67] was used to decompose the total interaction energy (Etot) into electrostatic (Eelec), dispersion (Edisp), exchange–repulsion (Eexch–repul), induction (Eind), and higher-order energy terms (δHF) according to the expressionGiven the size of the fragments involved in this study, we used density-fitting approximations to DFT-SAPT (DF-DFT-SAPT)[68] as implemented in the MOLPRO package.[69] The LPBE0AC exchange-correlation potential (which uses the localized Hartree–Fock scheme in the exchange part of the PBE0AC functional) was used in conjunction with the aug-cc-pVDZ basis set.[68] The asymptotic correction was incorporated by using the ionization potential of the respective monomers and the energy of the corresponding highest occupied molecular orbital.

Results and Discussion

Selected AChE Inhibitors

A congeneric set of seven AChE inhibitors differing only in their aromatic moieties was selected to focus on the arene–arene interactions between the inhibitors and the enzyme (Table ). The binding poses of the inhibitors in these complexes exhibited only minor differences,[45] enabling direct detailed comparisons of interactions with amino acid residues. The chosen inhibitors had half-maximum inhibitory concentration (IC50) values ranging from 7 to 139 μM. The p-Cl, p/m-CF3, and m-OCH3 compounds were the most potent inhibitors; the p-F and p-CH3 compounds were less potent, and the unsubstituted inhibitor (p-H) was the weakest. The crystal structures showed that the sulfonamide diethylamine moieties of all the inhibitors had overlapping positions and interactions in the lower region of the active site gorge (CAS), but there were differences in the upper region (PAS), mainly between benzylic and non-benzylic compounds but also within each class (Figure ). We therefore hypothesized that the differences in the inhibitors’ potency could be explained by considering their interactions in the PAS region, assuming similar entropic and desolvation effects within each inhibitor class.

Arene–Arene Interactions in the Protein Milieu

DFT optimizations of the crystal structure geometries were performed to determine positions of hydrogens and to allow relaxation of the amino acid side chains and inhibitor. To this end, a chemical cluster approach[55,56] was applied to the crystal structure conformations, with the whole inhibitor and the residues flanking the binding site gorge being included in the DFT geometry optimization. Minor movements of inhibitors and amino acid residues occurred during this step; the ring center distances of the arenes in the optimized inhibitors moved by 0.52 ± 0.15 Å on average compared to the crystal structures.

Four different types of arene–arene interactions between inhibitors and amino acid residues in AChE were identified (Figure ). The first type was a parallel-displaced interaction (i) observed in heterodimers consisting of a para-substituted benzylic arene (p-H, p-F, p-Cl, p-CH3, and p-CF3) and Tyr341 in which the planes of the two aryl rings were “tilted” (i.e., nonparallel). A similar tilted parallel-displaced heterodimer was also formed between non-benzylic meta-substituted arenes (m-CF3 and m-OCH3) and Phe338. The second type was a T-shaped interaction (ii) formed between the hydrogen(s) of the meta-substituted aromatic moieties and the aromatic ring of Tyr341 (Figure d,e). The third (iii) and fourth types (iv) were only formed between the benzylic compounds and AChE; type iii interactions involved contact between the benzenes’ para-substituents and the indole of Trp286, while type iv interactions involved an arene–hydroxyl contact with Tyr124 (Figure b,c). Tyr124 also formed a classical hydrogen bond with the sulfonamide moieties of all of the inhibitors. Notably, the angles and distances of all these interactions (other than the classical hydrogen bonds) differ from those in the corresponding free heterodimers at energy minima but are consistent with the specified interaction types.[36]

Interaction energies for the identified amino acid residue/inhibitor pairs were estimated by using coupled cluster theory (CCSD(T)), which is often termed the gold standard of quantum chemistry,[62,63] and by using two DFT methods (PBE[64]-D3[59] and MN15[66]). The interaction energy trends estimated by the DFT calculations were consistent with those of the coupled cluster calculations, although there were differences in absolute energy values (Figure S1 and Table S2). For example, the PBE-D3 method underestimated the energies of stacking and T-shaped interactions by ∼9% and ∼15%, respectively, while the MN15 method underestimated both by ∼23% compared to CCSD(T). The differences in the inhibitors’ interaction energies were analyzed based on three factors: the distances and angles of the arene interactions, substituent-related effects on the electron density distribution of the aromatic rings, and direct local interactions of the substituents with the interacting arenes.

Parallel-Displaced Stacking Interactions

The interaction energies of the parallel displaced stacking interactions calculated by using CCSD(T) were significantly attractive, ranging from −5.7 to −3.8 kcal mol–1 (Figure ), and correlated to the measured pIC50 values of the inhibitors toward AChE (Figure S2). The interaction energies were also in the same range as those reported previously for similar benzene heterodimers.[15,16,70] All inhibitors in the same class exhibited similar ring-stacking geometries (Figure ), and the interaction energies did not correlate with the center-to-center distances. Geometric differences therefore do not contribute greatly to the observed differences in interaction energies. However, the arenes of the non-benzylic compounds were more parallel to the amino acid residues’ arenes than the benzylic ones, leading to slightly stronger interactions for inhibitors with comparable substituents. The interaction strengths of the protein-milieu geometries thus clearly depended on the substituents: Cl/CF3/OCH3 yielded stronger interactions than F/CH3, which in turn yielded stronger interactions than H (Figure ). The finding that any substituent increased the interaction strength is consistent with results obtained for optimized sandwiched or parallel-displaced benzenes.[15,19,25] The substituents’ effect on the π-electron density distribution of the aromatic rings also did not seem to affect the strengths of these arene–arene interactions: in non-benzylic inhibitors, the electron-withdrawing group CF3 and the electron-donating group OCH3 had similar effects on interaction strength. Likewise, CF3 and Cl substituents on benzylic inhibitors yielded similar interaction strengths despite their different electronic effects. This is also supported by the weak correlation (R2 = 0.74) between the sum of the substituents’ Hammett σ constants (∑σ)[71,72] and the calculated interaction energies (Figure S3).

Inspection of the electrostatic potential maps revealed that all the benzylic inhibitors formed complexes featuring an attractive interaction between a hydrogen atom on the inhibitor’s arene and the proximal vertex of the Tyr341 next to the hydroxyl group (Figure ). This interaction was stronger in inhibitors with more electronegative substituents (i.e., F, Cl, and CF3) than in those substituted with H or CH3, suggesting that it was enhanced by local dipole moments that increased the partial positive charge of the interacting hydrogens. Furthermore, the geometries of arenes with the electron-rich substituents Cl or CF3 enabled the formation of an additional attractive interaction between these substituents and the proximal H-vertex of the Tyr341, but F appears too small to form this interaction (Figure ).

The two non-benzylic inhibitors differed slightly in binding mode (Figure ): the m-OCH3- substituted arene formed a better stacking overlap, while the m-CF3-substituted arene was positioned more out of plane, enabling the formation of a localized interaction between a sulfonamide oxygen and a hydrogen of Phe338.

To elucidate the different theoretical energy contributions to the interaction strengths, we investigated the van der Waals (vdW) area of the inhibitors’ substituents and found a weak correlation (R2 = 0.72) between interaction energies and this simple estimate of dispersion. However, as expected,[20,32] the vdW area alone could not fully explain the observed differences: the substituents CH3 and F yielded similar interaction energies despite having very different vdW areas. We therefore performed a decomposition of the interaction energies using DFT-SAPT (Figure a, Figures S1 and S4, and Table S3). The interactions were mainly stabilized by dispersion, which accounted for 62–72% of the total attractive energy. The electrostatic term was the second biggest contributor, accounting for 20–30% of the total (Figure b). Induction (2–4%) and higher order energy terms (5–6%) made minor contributions.

Figure 3

(a) Interaction energy components for parallel-displaced stacking interactions between benzylic inhibitors and Tyr341 and between non-benzylic inhibitors and Phe338 calculated by using DFT-SAPT. (b) Contributions (in percent) of the different attractive energy components to the total attractive interaction energies for each heterodimer.

The DFT-SAPT results for the benzylic inhibitors showed that the electrostatic contribution to the attractive energy was slightly greater for interactions involving halogens than for those involving methyl groups or hydrogen, supporting the conclusion that local dipole moments strengthened attractive interactions. The different electronic nature of Cl and F also appeared to influence the theoretical energy contributions: the “harder” F in CF3 resulted in a stronger electrostatic interaction but higher exchange repulsion, while the “softer” Cl yielded weaker electrostatics but also reduced exchange repulsion. For the non-benzylic inhibitors, m-OCH3 had a larger dispersion contribution, in keeping with its better arene–arene overlap, while the larger electrostatic contribution of m-CF3 is indicative of an additional direct, local substitution effect that is also visible in the electrostatic potential maps.

There is still controversy regarding the usefulness of Hammett constants as predictors of the interaction energy of sandwiched configurations of benzene dimers in optimal systems,[19,32,34] but recent studies suggest that they are reasonably good predictors for electron-withdrawing substituents (∑σ > 0).[31,33] However, the results obtained here suggest that Hammett constants are not useful for this purpose in the case of suboptimally stacked heterodimers in a protein milieu because they only correlated to the electrostatic component. It remains to be seen whether this is generally true for these kinds of tilted offset parallel stacking interactions.

T-Shaped Interactions

T-shaped interactions were formed between the two non-benzylic inhibitors’ arene hydrogens and the π-electrons of Tyr341. The interaction energies for the m-CF3 and m-OCH3 inhibitors were −4.7 and −3.6 kcal mol–1, respectively, by using the coupled cluster method. This is ∼1–2 kcal mol–1 stronger than the values reported for T-shaped benzene–phenol complexes.[73] The center-to-center distances and angles of the two dimers differed only slightly, but a shift in ring position allowed two H atoms to interact with Tyr341 in the m-CF3 complex compared to one in the m-OCH3 complex (Figure a). The strong electron-withdrawing effect of CF3 compared to OCH3 reduced the partial charges of the donating hydrogens and thus enhanced the attractive interaction and contributed to the 1.1 kcal mol–1 stronger interaction of the former. Electron-withdrawing substituents on the hydrogen-donating partner were previously shown to stabilize complexes of this type.[74] Direct interactions with the substituents on the other ring had only minor effects on T-shaped interactions because of the long distances involved. The DFT-SAPT energy profiles showed that the T-shaped interactions were dispersion-dominated (61% and 67%), with electrostatics contributing 30% and 24%, respectively, while induction (∼4%) and higher order terms (∼6%) again made minor contributions (Figure b, Figure S1, and Table S4). The electrostatic component was more pronounced in the case of m-CF3, further supporting the conclusion that electron-withdrawing substituents strengthen T-shaped interactions (Figure c). Liu et al.[73] reported dispersion and electrostatic contributions of 71% and 23%, respectively, to the total attractive energy for the T-shaped benzene–phenol complex. These values are more similar to those we obtained for the m-OCH3 complex than for the m-CF3 complex. The contributions of dispersion and electrostatics to the T-shaped interactions were comparable to those for the parallel-displaced case and also showed that strengthening the electrostatic component increased the interaction strength.

Substituent–Arene Interactions

The substituents on the benzylic inhibitors formed an interaction with the indole of Trp286 in AChE. This is the only interaction between the inhibitors and AChE in the PAS that depends solely on a direct contact with the substituent. The interaction strengths were between −2.7 and −1.4 kcal mol–1, with Cl substitution yielding a significantly stronger interaction than H and F (Figure a). The arene rings’ center-to-center distances were similar for all complexes (6.15 ± 0.05 Å). Consequently, the substituents’ distances to the tryptophan varied with their vdW volumes, with CH3 and CF3 and Cl being closest while F and H were further away (Figure a and Figure S5). The electrostatic potential maps (Figure S5) showed that both Cl and CH3 interacted directly with the proximal vertex of the indole moiety, accompanied by an additional contact between the hydrogen next to the substituent and the proximal edge of the indole moiety. CF3 had only one close contact with the closest H-vertex of the indole, while the H and F substituents were more distant, allowing only for hydrogen–hydrogen contacts (Figure S5).

Figure 5

(a) Substituent interactions between benzylic inhibitors and Trp286 with DLPNO-CCSD(T) interaction energies in kcal mol–1 (and PBE-D3 values in parentheses) and the distances between the indole ring center and the closest atom in the substituent (given in Å and shown using dashed lines). (b) Interaction energy components calculated using by DFT-SAPT. (c) Contributions (in percent) of the different attractive energy components to the total attractive interaction energies.

The DFT-SAPT analysis showed that these interactions were also dispersion-dominated (Figure b,c, Figure S1, and Table S5). Accordingly, the larger substituents (Cl, CH3, and CF3) all had similar contributions to the attractive energies, although CF3 had a higher exchange repulsion suggesting that it was too close to Trp286 to maximize the interaction strength. It appears that Cl represents a beneficial trade-off between size and electronic properties, allowing it to come close enough to the indole moiety to form a strong contact without too much repulsion. Our findings are coherent with a previous study by Imai et al.,[75] where they reported the presence of Cl–arene interactions in a selected set of protein–ligand crystal structures and concluded that the interactions were clearly attractive and dispersion dominated.

Hydroxyl–Arene Interactions and Hydrogen Bonding

The arene–arene interactions between the inhibitors and the Tyr124 of AChE were the strongest of those studied here, with energies between −9.7 and −7.2 kcal mol–1 (Figure ). The classical hydrogen bond formed between the inhibitors’ sulfonamides (donor) and the tyrosine hydroxyl (acceptor) had similar geometries in terms of distances and angles, which varied by at most 0.6 Å and 6.7°, respectively. The estimated hydrogen-bonding strengths of the isolated sulfonamide–phenol complexes were also similar for the inhibitors; the average interaction energy was −6.9 ± 0.19 kcal mol–1 (PBE-D3). The additional O–H···arene contacts of the benzylic inhibitors thus strengthened their interactions with AChE by 0.5–4.0 kcal mol–1, depending on the substituents (Figure ). The distances between the O–H and the inhibitors’ ring center ranged from 2.48 to 2.90 Å and correlated to the interaction energies; shorter distances accompanied stronger interactions (R2 = 0.85, see Figure S6). Liu et al.[73] obtained a hydroxyl–arene interaction energy of −4.8 kcal mol–1 for a similar benzene–phenol complex, confirming that interactions of this type can be quite strong.

The strengthening of the interaction due to substituents decreased in the order Cl > F/CF3 > H/CH3. Inspection of geometries and electrostatic potential maps revealed that for the p-Cl system the O–H···arene contact was hydrogen-bond-like, with the hydroxyl hydrogen projecting toward the π-electrons of the inhibitor’s arene moiety, unlike in the p-H and p-CH3 cases (Figure ). A similar but less pronounced tendency was observed in the p-F case, while the p-CF3 case was more complicated because CF3 formed additional direct, local interactions. The electron-rich Cl increased the attraction between the π-electrons and the hydroxyl hydrogen, resulting in a shorter distance and a stronger interaction. The differences in the interaction energies of the inhibitors in complex with Tyr124 were thus attributable to the substituents’ differing effects on the nature of the hydroxyl–arene contacts.

The DFT-SAPT energy decomposition analysis revealed only minor differences between the inhibitors (Figure a, Figures S1 and S7, and Table S6). The interactions were driven slightly more by electrostatics (44–46%) than dispersion (36–43%), unlike the other interactions studied here. Additionally, the induction and higher order terms were more significant for this type of interaction, contributing 8–11% and 6–9% of the total energy, respectively (Figure b). Domination by the electrostatic and induction components would be expected for a classical hydrogen bond,[76] but dispersion also contributed significantly to this mixed-type interaction. In the equivalent interaction between phenol and benzene, the contributions of dispersion, electrostatics, and induction were 55–62%, 29–32%, and 9–13%, respectively, indicating a more dispersion-dominated interaction.[73]

Figure 7

(a) Interaction energy components of the hydrogen bond/hydroxyl–arene contact calculated by using DFT-SAPT. (b) Contributions (in percent) of the attractive energy components to the interaction energy.

Collective Interactions, Summed Energies, and Predictions

The strong interaction energies of the classical hydrogen bond formed by all inhibitors indicate that this interaction is vital for the positioning of inhibitors in the PAS. However, the strengths and geometries of this interaction were similar for all inhibitors, so it cannot be responsible for the differences in potency. Instead, these differences were attributed to the different electronic properties of the substituents, which generated local dipole moments that enabled beneficial arene–arene interactions. These interactions were strongest in the p-Cl case, where the interacting amino acid residues and the arene of the inhibitor itself had different conformations compared to the p-H case. Geometrical differences compared to p-H were also seen for p-CF3, although here the effect was mainly due to the need to optimize the contacts of the three fluorines with the arene hydrogens of the amino acid residues. This shows that there is a delicate balance between maximizing electrostatic interactions and minimizing exchange–repulsion. In contrast to the direct fluorine contacts of the benzylic inhibitors, the CF3 group had a strong electron-withdrawing effect in the non-benzylic case, enabling the vertex hydrogens to participate in a stronger T-shaped interaction.

The summed DLPNO-CCSD(T) interaction energies for the studied complexes ranged from −17.7 to −13.2 kcal mol–1, going from strongest to weakest in the order p-Cl/m-CF3 > m-OCH3/p-CF3 > p-CH3/p-F > p-H. The order of interaction strengths correlated to the order of inhibitor potencies and plotting the summed interaction energies against pIC50 revealed a correlation (R2 = 0.71), although the non-benzylic inhibitors deviated from the trend (Figure S8). Plotting the results for the benzylic inhibitors alone yielded a strong correlation (R2 = 0.89, Figure S9), indicating that the summed energy is a good predictor for ranking this compound class in terms of binding strength. It also confirms our initial hypothesis that differences in inhibitor potencies are due to interactions in the PAS region.

Conclusions

Arene–arene interactions between inhibitors and AChE have been studied by modeling complexes observed in the protein interior in crystal structures. The geometries of these heterodimers differed substantially from those of small optimized model systems but were consistent with the established definitions of parallel stacking, T-shaped, and hydroxyl–arene interactions. Despite these geometric differences, the interaction strengths of the studied interactions in the protein milieu were in the same range as those reported for optimal benzene (hetero)dimers. Additionally, the energy decompositions of the arene–arene interactions were similar to those for geometry-optimized model systems. However, the studied complexes differed from previously reported model systems in some important ways. Notably, it appears that substituents have an even larger impact on arene interactions in complexes with tilted and offset geometries such as those in proteins. These substituent effects were electrostatics-driven, although the dominating component in the interactions was dispersion. The results also showed that arenes substituted with halogens interacted beneficially with aromatic amino acid residues. The halogens’ high electronegativity induced local dipole moments, which enhanced the electrostatic interactions. Additionally, their electron distributions enabled direct interactions with the electron-poor arene hydrogens of the amino acid residues. These findings provide a physical explanation for the impact of halogens such as Cl and F as aromatic substituents (or parts thereof, as in CF3) when optimizing protein–ligand interactions in drug discovery. Finally, the calculated summed interaction energies of the benzylic inhibitors correlated well with experimental potency data, highlighting the potential of using quantum mechanical calculations to support the design of new analogues.