Scaffold Effects on Halogen Bonding Strength - Journal of Chemical

Subscriber access provided by United Arab Emirates University | Libraries Deanship

Pharmaceutical Modeling

Scaffold Effects on Halogen Bonding Strength Andreas Lange, Johannes Heidrich, Markus O. Zimmermann, Thomas Exner, and Frank M. Boeckler J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00621 • Publication Date (Web): 10 Jan 2019 Downloaded from http://pubs.acs.org on January 11, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Scaffold Effects on Halogen Bonding Strength Andreas Lange†,§, Johannes Heidrich†, Markus O. Zimmermann†, Thomas E. Exner†,‡,#, and Frank M. Boeckler†,‡,* † Molecular Design and Pharmaceutical Biophysics, Institute of Pharmaceutical Sciences, Eberhard Karls University Tuebingen, Auf der Morgenstelle 8, 72076 Tuebingen, Germany ‡ Center for Bioinformatics Tuebingen (ZBIT), Eberhard Karls University Tuebingen, Sand 1, 72076 Tuebingen, Germany HALOGEN BONDING, VMAX, ADDUCT FORMATION, CORRELATION, PREDICTION

ABSTRACT: Halogen bonds have become increasingly popular interactions in molecular design and drug discovery. One of its key features is the strong dependence of the size and magnitude of the halogen's σ-hole on the chemical environment of the ligand. The term σ-hole refers to a region of lower electronic density opposite to a covalent bond, e.g., the C-X bond. It is typically (but not always) associated with a positive electrostatic potential in close proximity to the extension of the covalent bond. Herein, we use a variety of 30 nitrogen-bearing heterocycles, halogenated systematically by chlorine, bromine, or iodine, yielding 468 different ligands that are used to exemplify scaffold effects on halogen bonding strength. As a template interaction partner, we have chosen N-methylacetamide representing the ubiquitously present protein backbone. Adduct formation energies were obtained at a MP2/TZVPP level of theory. We used the local maximum of the electrostatic potential on the molecular surface in close proximity to the -hole, VS,max, as a descriptor for the magnitude of the positive electrostatic potential characterizing the tuning of the -hole. Free optimization of the complexes gave reasonable correlations with VS,max, but was found to be of limited use, because considerable numbers of chlorinated and brominated ligands lost their halogen bond or showed significant secondary interactions. Thus, starting from a close to optimal geometry of the halogen bond, we used distance scans to obtain the best adduct formation energy for each complex. This approach provided superior results for all complexes exhibiting correlations with R2 > 0.96 for each individual halogen. We evaluated the dependence of VS,max from the molecular surface onto which the positive electrostatic potential is projected, altering the isodensity values from 0.001 au to 0.050 au. Interestingly, the best overall fit using a third-order polynomial function (R2 = 0.99, RMSE = 0.562 kJ/mol) with rather smooth transitions between all halogens was obtained for VS,max calculated from an isodensity surface at 0.014 au.

Introduction Halogen bonding, an attractive interaction between chlorine, bromine, or iodine and any kind of nucleophile (bearing n- or -electrons),1 has evolved in the past 10 years as a useful tool in life sciences and drug discovery.2-11 To date, we have studied halogen bonding contacts with different interaction partners, such as the backbone carbonyl6, 10, 12, the sulfur of methionine5, 9, the nitrogen of histidine8, or carboxylates and carboxamides (aspartate, glutamate, asparagine, and glutamine)13 in protein binding sites to facilitate their use in molecular design. For the successful application of halogen bonding in drug discovery, not only an improved knowledge about favorable geometries for targeting the various interaction partners is key, but also a better understanding of the impact of various decorated scaffolds on halogen bonding strength is likewise required.6, 7, 14

The aromatic or heteroaromatic scaffold and the attachment position of a halogen within the scaffold will certainly influence the strength of a halogen bond and, therefore, the size and shape of the positive electrostatic potential in proximity to the σ-hole.15 Polarization by electron-withdrawing atoms or groups in the scaffold and shifts in the electrostatic potential by hydrogen bond donors or acceptors in close proximity to the halogen can tune the σ-hole, impose changes in the geometry and help to form “binding motifs” (a cooperative network of interactions).7, 14 Of course, similar effects are observed when different substituents (e.g., F, NO2, OH, etc.) are introduced into the same scaffold. Riley et al. pioneered to study tuning effects on halogen bonding in fluorinesubstituted halobenzenes.16, 17 They have shown that in extremely tuned molecules, like in a pentafluorohalobenzene, the interaction energy can be approximately doubled.

1 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Thus, the interesting question arises, what degree of tuning of a certain halogen is favorable, while still retaining drug-like molecular properties? Introducing fluoro substituents does not only affect the σ-hole, but also can generate exciting new binding motifs based on the synergistic modulation of both individual interactions.7, 18, 19 In a recent study combining computation and crystallography, Fanfrlík et al. demonstrated the potential of halogen bond tuning in a designed series of aldose reductase inhibitors based on a systematic variation of the fluorine substitution pattern.20 Kolar et al. analyzed a subset of the ZINC database containing fragment-like compounds to provide a simplified description of the hole dependent on its magnitude, size, angular deviation, and range. They conclude that there is a low effect of intramolecular polarization on the spatial position of VS,max, whereas the other characteristics are more strongly influenced.21 In this systematic study, we strive to gain further insights into tuning effects and geometry changes by investigating the electrostatics and the adduct formation energies of more than 450 fragments of 30 different nitrogen-bearing heterocycles substituted by chlorine, bromine or iodine at the MP2 level of theory using a TZVPP basis set. Such systems are quite common ligand systems used in drug discovery and medicinal chemistry.

Results and Discussion Strategy overview. Starting from a 1H-indole system, we transformed the nitrogen to every possible and reasonable position, gaining four more ligand systems: indolizine, 2Hcyclopenta[c]pyridine, 1H-cyclopenta[b]pyridine and 2Hisoindole. To further expand our field of heterocycles, we introduced a second nitrogen atom into every reasonable position yielding another 16 heterocycles bearing two nitrogen atoms (Figure 1).

Figure 1. Overview of the model systems of nitrogen-bearing heterocycles used in this study. For the structures marked with an asterisk (*) tautomeric forms are differentiated.

Page 2 of 14

For 1H-indazole, 1H-pyrrolo[3,2-b]pyridine, 1H-pyrrolo[2,3b]pyridine, 1H-pyrrolo[3,2-c]pyridine, 2H-pyrrolo[3,4c]pyridine, and 3H-cyclopenta[d]pyrimidine (marked with an asterisk in Figure 1), both tautomeric forms were separately evaluated. As a comparison, we added the three 5-ring N-heterocycles 1H-pyrrole, 1H-pyrazole and 1Himidazole to our series of model systems. A far more extensive, but less systematic diversity of heterocycles is used in a parallel study aiming at building machine learned models for the fast and efficient prediction of the strength of σ-hole tuning.22 The isolated geometry of every hetrocycle was energyminimized using MP2/TZVPP (as implemented in Turbomole 6.4). After optimization, all carbon-bonded hydrogen atoms were iteratively substituted by chlorine, bromine, or iodine (Figure 2a). In this study, we only considered mono-halogenated ligand systems. We received 30 different nitrogen-bearing heterocycles (including tautomers) as basis structures (Figure 1) and a total of 468 halogenated scaffold model systems. All halogenated structures were again optimized at the MP2/TZVPP level. Subsequently, the electrostatic potential was calculated using a 2D and 3D grid for each of the 468 fragments to visualize the difference of the σ-hole in-plane and in three-dimensional space (Figure 2b+c). Using the 3D grid results, we automatically assessed the local maximum of the positive electrostatic potential on the isodensity surface in close proximity to the -hole, the so-called Vmax value (or better VS,max). From the VS,max value, we wanted to gain insights how the scaffold and the attachment position of the halogen influences the tuning of the -hole (Figure 2d). To reduce the immediate effects of the negative electrostatic potentials of neighboring nitrogen atoms, distorting the electrostatics around the halogen, we initially decided to use an isodensity surface more proximal to the nucleus (at 0.02 au) than usually applied (at 0.001 au) for the calculation of VS,max. It should be noted that we eventually made a systematic comparison, investigating, how the choice of isodensity surface (between 0.001 au and 0.05 au) and the respective VS,max values influence the finally obtained correlations. In addition, we calculated the interaction energy between each halogenated scaffold and N-methylacetamide (BB), serving as a small model system for the protein backbone.6 To generate starting structures for further optimization, we matched each (hetero)arylhalide onto the halobenzene of the optimized halobenzene···BB complex (dO···X: 302 pm for iodine, 304 pm for bromine, and 312 pm for chlorine, with a σ-hole αO···X-C angle of 175.6° for iodine, 177.4° for bromine, and 171.2° for chlorine), as published in Wilcken et al.6 Each oriented and superposed complex was optimized without constraints at the MP2/TZVPP level of theory (Figure 2e). Finally, we analyzed the data and obtained correlations of the calculated adduct formation energies (ΔE = (Eadduct – (EBB + Escaffold)) and the VS,max value of each halogenated scaffold (Figure 2f). To evaluate positive or negative tuning effects, we compared the resulting energies and VS,max values to the simple halobenzene systems, as a reference.


Page 3 of 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling generated in this way, were subjected to a free geometry optimization.

Figure 2. Overview of the workflow applied for the systematic evaluation of σ-hole tuning effects. (a) represents the systematic generation of halogenated scaffolds, illustrated using pyrazolo[1,5-a]pyridine. The resulting electrostatic features around the halogen are shown for 2-iodopyrazolo[1,5-a]pyridine using a 2D (b) and 3D grid (c). The most positive area on the isodensity surface of the σ-hole, VS,max, is automatically selected and highlighted in (d). Adduct formation energies are calculated on a MP2/TZVPP level of theory (e) and correlations between VS,max and these interaction energies are derived (f).

Calculating the Adduct Formation Energy. Model systems for quantum mechanical calculations are usually as small as possible and as large as necessary to characterize the essential chemistry of the interaction. Based on Riley et al., who used acetone and bromobenzene23 as model systems, we employed Nmethylacetamide to represent the protein backbone6. Despite known tendencies that MP2 can significantly overestimate the halogen bond strength24, 25, we found in our initial study of this system6 that the combination of MP2 with a triple-ζ basis set (TZVPP) without Counterpoise correction yielded comparable results to CCSD(T) reference calculations using a complete basis set (CBS) extrapolation scheme. As a consequence, we have employed MP2/TZVPP for all further calculations in this study. Based on a representative selection of examples, we suggest that the BSSE obtained by using Counterpoise correction is about 2,5 to 3 kJ/mol for chlorine and bromine, while for iodine it is slightly higher (4 to 5 kJ/mol). It should be further noted that this model system represents the most ubiquitous interaction partner for forming halogen bonds in protein binding sites and, thus, is of most relevance to medicinal chemistry and molecular design. Small model systems tend to have a lower risk of forming unwanted secondary interactions, which can significantly limit the interpretability of the relative comparison intended in this study. In addition, the subsequent protocol aims at further minimizing this risk. After full optimization of each halogenated scaffold, starting geometries for the optimization of the complexes were generated by reconstruction of the halobenzene’s orientation in the halobenzene···N-methylacetamide complex6, i.e. matching of the rings and the C-X vector (see Supporting Information). Consecutively, all complexes

After optimization, the vast majority of the iodine heterocyclic complexes remain virtually unchanged. For bromine, and even more pronounced for chlorine, a change in the interaction geometries occurs more frequently due to secondary interactions. In some cases, the chlorine halogen bonds were completely lost due to newly formed hydrogen bonds between the heteroatoms of the scaffold and the backbone model system. This trend was somewhat expected due to the decreasing magnitude of the σ-hole for the lighter halogens. Consequently, reasonable interaction geometries could not be obtained for all 468 complexes after optimization. To automatically assess the loss of the halogen bond, we employed XBScoreQM, a scoring function for halogen bonds addressing the carbonyl oxygen of the N-methylacetamide model system10. We checked the persistency of the -hole interaction by analyzing every step of each free optimization with this scoring function. The heterocyclic complexes started with XBScoreQM values of 0.96, 0.87, and 0.67 for iodine, bromine, and chlorine, respectively. Small deviations in the XBScoreQM do not necessarily indicate a weakened -hole interaction, as the scoring function was developed using untuned, symmetric halobenzenes. In contrast, the herein investigated ligands are expected to show significant tuning effects and can exhibit asymmetric VS,max positions based on electrostatic influences of their scaffold nitrogens. However, complexes with an XBScoreQM ≤ 0.3 typically indicate insufficient σ-hole contacts with the BB model system and, thus, have been removed from the data set. Based on this crude exclusion criterion, in 11 iodine complexes (7%), 35 bromine complexes (22%), and 80 chlorine complexes (51%) the halogen bond was lost. This is the best result obtained, even after careful reconsideration of the starting geometries. The finally accepted optimized complexes had average XBScoreQM values (± standard deviations) of 0.95 ± 0.03 for iodine (n = 145), 0.83 ± 0.11 for bromine (n = 121), and 0.66 ± 0.07 for chlorine (n = 76). For most of the retained systems, the score is hardly impaired during optimization. Still, in some of the complexes containing mainly bromine and chlorine, a perceivable decrease of the score indicates weakening of the -hole contact and the onset of secondary interactions. Correlation between VS,max (at 0.02 au electron isodensity) and the Adduct Formation Energy. We perceive VS,max as a simple, but efficient descriptor encoding the characteristic positive electrostatic potential representing the -hole (and neighboring influences) of the ligand. Thus, it appears to be a valid hypothesis that it can indicate the strength of the resulting halogen bond involving the -hole of the ligand and a nucleophile, such as the carbonyl n-electrons in N-methylacetamide. Politzer et al. have originally reported, that the more positive VS,max becomes, the more attractive the resulting interaction energy (ΔE) will be11.



Page 4 of 14

Figure 3. Linear correlation plots of adduct formation energies vs. VS,max and depiction of selected outliers. All complexes showing halogen bonds after optimization are plotted for (a) iodine (n = 145), (b) bromine (n = 121) and (c) chlorine (n = 76). Geometries of labelled outliers (A – I) are depicted underneath the individual plots. Red dotted lines indicate the halogen bond interaction, whereas yellow dotted lines highlight secondary interactions. Some complexes, where the ligands visibly engage in secondary interactions, were omitted from the plots (d - f), as discussed. 135 complexes remain for iodine (d), 100 for bromine (e), and 73 for chlorine (f). Plots were prepared with KaleidaGraph26.

Here, we first focused on linear correlations between VS,max and the adduct formation energies. Figure 3 shows the results for all complexes featuring halogen bonds after the optimization. All three halogenated systems show a reasonable linear correlation with coefficients of determination R2 of 0.88 for iodine, 0.66 for bromine, and 0.60 for chlorine, respectively. Figure 3d-f shows the modified linear correlations after removing outliers that were caused by visible secondary interactions.

Unsurprisingly, the overall R2 is improved for all three halogen systems. To illustrate the challenge of secondary interactions in the correlations, we highlight and depict some of these outliers. Most of these outliers were found to be involved in secondary interactions between the halogen and the methyl groups of the BB model system such as shown for the iodine outliers A, B, and C in Figure 3a.


Page 5 of 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


While the VS,max values of these structures are rather weak for iodine (0.166 au for A and B, 0.169 au for C), close contacts of each iodine to the CH3-groups of Nmethylacetamide and, possibly, attractive electrostatic interactions of the nitrogen in ortho position of each iodine towards N-methylacetamide, are explanations for the significantly improved total interaction energies. These additional interactions are enabled by an orthogonal binding mode towards the plane of the backbone amide model system. It should be noted that this binding mode is still one of the hotspots for the -hole interaction, reflecting that there are also attractive interactions possible with the -electron density of the peptide backbone.12, 27 For bromine, similar secondary interactions are found. In the complexes close to outlier D, deviations are all caused by the same CH3 contacts as exemplified for iodine. In contrast, the binding mode of the ligands in complex E and F is not orthogonal to the plane of the backbone model. In these examples, the bromine atoms still engage in close contacts to the CH3 group next to the carbonyl oxygen. For chlorine, only three significant outliers remain. H shares the same binding mode and secondary interactions as discussed for E and F. For G and I, a σ-hole angle C-ClO of about 147° indicates that the hole interaction may be very weak at best. Still, there are significant interactions between hydrogen atoms of the heteroaromatic scaffold and the -electrons of the backbone amide model (particularly close to the nitrogen). Overall, we observed 10 outliers for iodine and 3 for chlorine. In the case of chlorine, this is likely the result of the previously described frequent loss of the -hole interaction, leading to an exclusion of already more than 50% of all the investigated systems.

Figure 4. Polynomial curve fit of all three halogens including (a) all complexes showing halogen bonds after free optimization (n = 343) and (b) all complexes without visible secondary interactions (n = 308). The larger, differently colored dots represent the halobenzene references (yellow: chlorobenzene [VS,max = 0.117 au, ∆E = -8.2 kJ/mol], orange: bromobenzene [VS,max = 0.150 au, ∆E = -11.8 kJ/mol], cyan: iodobenzene [VS,max = 0.182 au, ∆E = -18.8 kJ/mol]).

For bromine, where fewer complexes had to be discarded due to the loss of their halogen bond, 21 outliers were identified showing secondary interactions. It is self-evident that the coefficient of determination increases considerably for all halogens: R2 of 0.95 for iodine and bromine; R2 of 0.81 for chlorine. However, the total number of omitted structures clearly indicates that the free optimization is not the ideal method. Next, we chose to investigate the hypothesis that there is one general equation comprising the relationship of VS,max to adduct formation energies for all three halogens (Figure 4). A simple third order polynomial function (y = ax3 + b) provides a rather reasonable curve fit for the entire dataset. Judged by the Akaike information criterion (AICc) the quality of this simple model is preferable over more complex third order polynomials or second order polynomials. The coefficient of determination R2 is 0.918 including secondary interactions (Figure 4a, degrees of freedom (df) = 341, RSME = 1.45 kJ/mol) and 0.971 without secondary interactions (Figure 4b, df = 306, RSME = 0.881 kJ/mol). The overall curve fit equations are y = -5.66 – 2274x3 (Figure 4a) and y = -4.95 – 2346x3 (Figure 4b), respectively. As references, the values of the halobenzene complexes are emphasized as differently colored squares in Figure 4 (iodobenzene in cyan, bromobenzene in orange and chlorobenzene in yellow). Almost all scaffolds show interaction energies stronger than their respective halobenzenes. Focusing on the systems without visible secondary interactions (Figure 4b, n = 308), the maximum increase in complex formation energy with respect to the halobenzene complexes is +54% for the iodine ligands, +73% for the bromine ligands, and + 83% for the chlorine ligands. Thus, the tuning effects of neutral, nitrogenbearing heterocycles as commonly used in medicinal chemistry and drug discovery on the halogen bonding strength can be quite significant. It is also clearly visible from these plots that a transition between different halogens is possible based on the scaffold. Chlorine can be tuned to show stronger halogen bonds than some of the bromine systems. Bromine can be tuned to show halogen bonds stronger than some of the iodine systems. Correlation between VS,max and Adduct Formation Energies, Calculated by a Distance Scan Approach. Based on the issues we experienced with some of the free optimizations of the complexes, we additionally considered another approach to derive the complex formation energies trying to reduce the influence of secondary interactions even more. The halogenated ligand models were superposed. The angle between the carbonyl group of N-methylacetamide and the halogen was set to 120° and the dihedral angle between the aromatic ring system of the σ-hole donor and the amide plane was set to 90°.28



Page 6 of 14

Figure 5. Linear and polynomial correlations for adduct formation energies calculated by a distance scan approach and 𝐦𝐚𝐮 VS,max values obtained from two different isodensity surfaces. Linear correlations using 𝐕𝟐𝟎 are depicted for (a) iodine 𝑺,𝐦𝐚𝐱 ligands, (b) bromine ligands, and (c) chlorine ligands. Corresponding linear correlations for the individual halogens using 𝐦𝐚𝐮 𝟐𝟎 𝐦𝐚𝐮 𝐕𝟏𝟒 𝑺,𝐦𝐚𝐱 are presented below (e-g). Best global fits of all halogens using a polynomial function of third order at (d) 𝐕𝑺,𝐦𝐚𝐱 𝟏𝟒 𝐦𝐚𝐮 and (h) 𝐕𝑺,𝐦𝐚𝐱 This standardized interaction geometry is very close to the global optimum of the halobenzene complexes, except for the orthogonal orientation of the two planes, which was chosen to minimize the risk of clashes and secondary interactions. We then performed distance scans, changing the distance between Car of the ligand system (to which the halogen is attached to) and the oxygen of Nmethylacetamide from 400 pm to 680 pm in steps of 10 pm. Each geometry obtained in this way, was subjected to a single point energy calculation at the MP2/TZVPP level of theory. The minimal energy was determined by regional regression analysis of the global minimum of energies of the distance scan, approximated by a quadratic equation in R.29 Without the free optimization, the risk of additional secondary interactions is much reduced. Figure 5 summarizes the results. In contrast to the free optimizations, where only 343 complexes (73.3%) showed -hole contacts and only 308 complexes were devoid of significant secondary interactions (65.8%), reasonable results were obtained for all ligands (n=468, 100%) as intended. The quality of the linear correlations (Figure 5ac) is much improved, yielding R2 values of 0.969, 0.979, and 0.988 for the chlorine, bromine, and iodine ligands, respectively. The best second or third order polynomial

curve fit for all ligand systems together (Figure 5d) was selected based on the AICc weights of the statistical model. A general polynomial function of third order (y = ax3 + bx2 + cx + d) with a = -7140.6 ± 1171.3, b = 2592.1 ± 553.5, c = 459.17 ± 85.43, and d = 22.521 ± 4.306 provides the best results for the entire data set (df = 464, R2 = 0.985, RMSE = 0.687 kJ/mol). Systematic Assessment of VS,max-dependence on the Electron Isodensity Surface and its Impact on Correlation Quality. Because of the significantly reduced scattering in the different halogens, it becomes apparent in the overall fit that the transition of chlorine to bromine ligands is much more heterogeneous than between bromine and iodine ligands. Thus, we reconsidered to investigate the dependence of VS,max from the electron isodensity surface onto which the electrostatic potential is mau projected. As stated before, we chose to use the V20 𝑆,max (VS,max value detected at an isodensity surface of 0.02 au) mau instead of the common V1𝑆,max , trying to limit the immediate effects of the negative electrostatic potentials of neighboring nitrogen atoms. Smaller values of the electron isodensity surface lead to smaller VS,max values. However, the three halogens are differently affected.


Page 7 of 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Dependence of linear, quadratic, and cubic curve fits on the definition of VS,max. Isodensity surfaces between 0.001 au and 0.05 au were used to calculate VS,max. The minimal and maximal isodensity surface, as well as the isodensity surface at 0.014 au, yielding the best overall results, are depicted in the two-dimensional electrostatic plots of the three halobenzenes. Correlations of determination and RMSE values are plotted for all curve fits using VS,max values obtained for the individual isodensity surfaces. Results for linear correlations of individual halogens are shown as crossed squares for chlorine (green), bromine (brown), and iodine (purple). Results for linear correlations (y = ax + b) of all halogen complexes are indicated by cyan upward pointing triangles. Results for quadratic correlations (y = ax2 + bx + c) of all halogen complexes are represented as blue diamonds. Results for cubic correlations (y = ax3 + bx2 + cx + d) of all halogen complexes are shown as dark-blue downward pointing triangles.

In consequence, we hypothesized that there might be an optimal isodensity surface, yielding a more homogeneous transition between the different halogens. We performed an analysis of all VS,max values for isodensities from 0.001 au mau mau (V1𝑆,max ) to 0.05 au (V50 𝑆,max ) and repeated all correlation analyses for first, second, and third order polynomials, including different numbers of fitted variables. The full results (plots and AICc analyses) can be found in the Supporting Information. In summary, the linear correlation for each individual halogen is favored for VS,max obtained from smaller isodensities (Figure 6). For the mau iodine ligands, the best model is found for V7𝑆,max (R2 = 0.990, RMSE = 0.344 kJ/mol). For the bromine ligands, mau V3𝑆,max yields the optimal linear model (R2 = 0.984, RMSE = mau 0.368 kJ/mol) and for the chlorine ligands, V2𝑆,max is the best 2 basis for a linear model (R = 0.982, RMSE = 0.347 kJ/mol). In contrast, the best overall model providing rather smooth transitions between the three halogens is observed mau for intermediate isodensities (V14 𝑆,max ).

Judged by Akaike weights, the best polynomial function to fit this data is y = ax3 + cx + d with a = -2717.2 ± 111.6, c = 104.58 ± 4.96, and d = 3.566 ± 0.386 (df = 465, R2 = 0.990, RMSE = 0.562 kJ/mol, Figure 5h), which is slightly preferred over other quadratic or cubic functions. The RMSE of the individual models should be considered with respect to the data range. This normalized RMSE is 2.38% for the final model comprising all halogens. In comparison, the normalized RMSE for the best individual models is 2.49% for iodine, 3.52% for bromine, and 2.96% for chlorine. Of course, the reported results are specifically reflecting the individual characteristics of the analyzed data set. Still, we found similar dependencies based on the definition of VS,max in a parallel study for ligands obtained from the Protein Data Bank, where the best overall model mau 22 was obtained for V18 𝑆,max .



Page 8 of 14

Figure 7. Key electrostatic features of the -hole and neighboring scaffold effects based on different attachment positions of iodine to pyrazolo[1,5-a]pyridine. 2D-plots for (a) 2-iodo-pyrazolo[1,5-a]pyridine, (b) 3-iodo-pyrazolo[1,5-a]pyridine, (c) 4-iodopyrazolo[1,5-a]pyridine, (d) 5-iodo-pyrazolo[1,5-a]pyridine, (e) 6-iodo-pyrazolo[1,5-a]pyridine, and (f) 7-iodo-pyrazolo[1,5a]pyridine. Positive to negative electrostatic potentials are color-coded by a red–white–blue gradient. Tuning effects and shape changes of the -hole and other positive electrostatic potentials (red) as interaction hotspots are indicated by differently sized orange arrows. Variations of the negative electrostatic potentials (blue) as interaction hotspots, based on the equatorial belt surrounding the halogen and the n-electrons of the heterocyclic nitrogen atoms are depicted by differently sized blue arrows. Plots were made using MOLCAD30, 31.

CONCLUSION In this systematic evaluation of scaffold effects of different nitrogen-bearing heterocycles on halogen bonding strength, we have focused on VS,max as a descriptor for halogen bonding strength. While free optimizations of the complexes are problematic in many cases of tuned systems, we suggest that the calculation of the adduct formation energies via distance scans avoiding secondary interactions is a very well-suited method. Most of the studied heterocycles had the potential for stronger binding to Nmethylacetamide than their respective halobenzenes. This tuning effect can be as high as +54% for iodine, +73% for bromine, and +83% for chlorine. The increase in halogen bonding strength can be easily correlated to VS,max for each halogen separately. Finding reasonable correlations for all ligands containing chlorine, bromine, and iodine together is much more challenging. Eventually, we have scrutinized the definition basis of VS,max and found that the best results were obtained mau for an isodensity surface of 0.014 au and the derived V14 𝑆,max . Variations from the classical 0.001 au isocontour level have been recently discussed.32 Polynomial functions of third order could best be fitted to the complete data set. In addition to characterizing and predicting22 halogen bonding strength, it is important to understand, how dramatically the electrostatic features surrounding the

halogen can deviate, even based on the attachment position to a particular scaffold. This is exemplified using all possible iodo-pyrazolo[1,5a]pyridines in Figure 7. In 4 of the 6 ligands, the tuning effect is strong enough that the positive potential predominates the entire 2D-plot. The neutral white line is not anymore enclosing the positive potential like in 2-iodopyrazolo[1,5-a]pyridine (Figure 7a) or 7-iodo-pyrazolo[1,5a]pyridine (Figure 7f), but the remainder of the bluecolored negative potential (Figures 7b-e). Even in the tuned state, the -hole and the halogen-surrounding negative belt can be rather symmetrical (Figure 7d). Particularly neighboring nitrogen atoms of the heteroaromatic scaffold can lead to fairly asymmetric systems (e.g. Figure 7f). These geometric features of the hole and their integration into more general binding motifs of the ligand are crucial to consider when aiming for a strong interaction network with the binding site. Hence, a better understanding of halogen bond tuning essentially includes VS,max as a descriptor for the magnitude of the hole and the resulting strength of the halogen bond, but also has to take into consideration the geometric complexity exemplified in Figure 7. When aspiring to implement halogen bonding in molecular design and drug discovery strategies, this diversity of aspects should be incorporated. All these efforts will lead to a broader recognition of halogen bonding as a useful tool for


Page 9 of 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

medicinal chemists in pharmaceutical research.

Journal of Chemical Information and Modeling industrial

and

academic

function fitted locally to the smallest energy values using R.29

ASSOCIATED CONTENT Methods MP2 Structure Optimization Quantum mechanical calculations were performed using TURBOMOLE Version 6.4, 7.1 and 7.233-36. A triple-ζ basis set (def2-TZVPP)37, 38 was used throughout the study. MP2 calculations were done in combination with the resolution of identity (RI) technique28, 38-40 and the frozen core approximation. The frozen core orbitals were defined using default settings by which all orbitals possessing energies below -3.0 au are considered core orbitals. The SCF convergence criterion was increased to 10-8 hartree for all calculations. 2D, 3D-Plots and VS,max Calculations All ligands were oriented after their optimization employing in-house scripts by placing the ring systems in the same plane and by aligning the C-X bond. 2D grid calculations were conducted for visualization purposes only. 3D grid calculations were conducted by expanding the surrounding box of a given molecule by 5 au in each Cartesian direction to ensure that the grid is big enough for the subsequent analysis. A grid point distance of 0.2 au was chosen for all directions. MOLCAD II v1.430, 31 was used for visualization purposes of electrostatic potentials and electron densities. Adduct formation energy calculation For the adduct formation energy calculations, we first superposed each ligand system to its corresponding, optimized halobenzene…N-methylacetamide structure, e.g. iodine ligands superposed on iodobenzene, bromine ligands superposed on bromobenzene and chlorine ligands superposed on chlorobenzene. Details are illustrated in Figure S1 (see Supporting Information). Afterwards each structure was optimized without constraints using TURBOMOLE 6.4.34 Distance scans A distance scan was applied to calculate a close-to-globalminimum adduct formation energy, while avoiding secondary interactions. The X-C axis of (hetero)aryl halides was oriented 120° to the carbonyl function of Nmethylacetamide. The dihedral angle defined by the plane of the amide of the backbone model system and the plane of the aromatic ring of the -hole donor was set to 90°. The distance between the carbon of C-X and the oxygen of the backbone model system was increased from 4.0 Å to 6.8 Å in increments of 0.1 Å, while maintaining the other geometrical parameters. Single point energies of the resulting geometries were obtained. Global minima were approximated by calculating the minimum of a quadratic

Supporting Information. Details on generation of starting structures for free QM calculations (Figure S1), goodness of fit of electron isodensity-dependent correlation of adduct formation energy and VS,max (Figures S2-S51) and their AICc values (Table S1-S3). This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION Corresponding Author * [email protected]

Present Addresses § Institute for Evolution and Biodiversity, University of Münster, Hüfferstrasse 1, D-48149 Münster, Germany # Douglas Connect GmbH, Technology Park Basel, Hochbergerstrasse 60C, CH-4057 Basel, Switzerland

Author Contributions A.L., J.H., and F.M.B. designed research and conducted data analysis; A.L. performed free QM optimization of the ligands and complexes, as well as 2D and 3D grid calculations; A.L. analyzed free optimizations and prepared correlations of VS,max and adduct formation energies; M.O.Z. provided tools for setting up scaffold orientation; J.H. performed QM distance scans, provided isodensity-dependent VS,max assessments and conducted full statistical analysis of all respective results; T.E.E. provided tools for calculating the 3D grid and helped with data refinement using MOLCAD; A.L., M.O.Z., and F.M.B. wrote the manuscript. All authors have given approval to the final version of the manuscript.

ACKNOWLEDGMENT The authors acknowledge support by the state of BadenWürttemberg through bwHPC.

REFERENCES 1. Desiraju, G. R.; Ho, P. S.; Kloo, L.; Legon, A. C.; Marquardt, R.; Metrangolo, P.; Politzer, P.; Resnati, G.; Rissanen, K. Definition of the halogen bond (IUPAC Recommendations 2013). Pure Appl. Chem. 2013, 85, 1711-1733. 2. Auffinger, P.; Hays, F. A.; Westhof, E.; Ho, P. S. Halogen bonds in biological molecules. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 16789-16794. 3. Politzer, P.; Lane, P.; Concha, M. C.; Ma, Y. G.; Murray, J. S. An overview of halogen bonding. J. Mol. Model. 2007, 13, 305311. 4. Hardegger, L. A.; Kuhn, B.; Spinnler, B.; Anselm, L.; Ecabert, R.; Stihle, M.; Gsell, B.; Thoma, R.; Diez, J.; Benz, J.; Plancher, J.-M.; Hartmann, G.; Banner, D. W.; Haap, W.; Diederich, F. Systematic Investigation of Halogen Bonding in Protein–Ligand Interactions. Angew. Chem. Int. Ed. 2011, 50, 314318. 5. Wilcken, R.; Zimmermann, M. O.; Lange, A.; Zahn, S.; Kirchner, B.; Boeckler, F. M. Addressing Methionine in Molecular



Design through Directed Sulfur–Halogen Bonds. J. Chem. Theory Comput. 2011, 7, 2307-2315. 6. Wilcken, R.; Zimmermann, M.; Lange, A.; Zahn, S.; Boeckler, F. Using halogen bonds to address the protein backbone: a systematic evaluation. J. Comput. Aided Mol. Des. 2012, 26, 935-945. 7. Wilcken, R.; Zimmermann, M. O.; Lange, A.; Joerger, A. C.; Boeckler, F. M. Principles and Applications of Halogen Bonding in Medicinal Chemistry and Chemical Biology. J. Med. Chem. 2013, 56, 1363–1388. 8. Lange, A.; Zimmermann, M. O.; Wilcken, R.; Zahn, S.; Boeckler, F. M. Targeting Histidine Side Chains in Molecular Design through Nitrogen–Halogen Bonds. J. Chem. Inf. Model. 2013, 53, 3178-3189. 9. Lange, A.; Gunther, M.; Buttner, F. M.; Zimmermann, M. O.; Heidrich, J.; Hennig, S.; Zahn, S.; Schall, C.; Sievers-Engler, A.; Ansideri, F.; Koch, P.; Laemmerhofer, M.; Stehle, T.; Laufer, S. A.; Boeckler, F. M. Targeting the Gatekeeper MET146 of C-Jun NTerminal Kinase 3 Induces a Bivalent Halogen/Chalcogen Bond. J. Am. Chem. Soc. 2015, 137, 14640-14652. 10. Zimmermann, M. O.; Lange, A.; Boeckler, F. M. Evaluating the Potential of Halogen Bonding in Molecular Design: Automated Scaffold Decoration Using the New Scoring Function XBScore. J. Chem. Inf. Model. 2015, 55, 687-699. 11. Politzer, P.; Murray, J. S.; Clark, T. Halogen bonding and other sigma-hole interactions: a perspective. PCCP 2013, 15, 1117811189. 12. Zimmermann, M. O.; Boeckler, F. M. Targeting the protein backbone with aryl halides: systematic comparison of halogen bonding and π⋯π interactions using N-methylacetamide. MedChemComm 2016, 7, 500-505. 13. Zimmermann, M. O.; Lange, A.; Zahn, S.; Exner, T. E.; Boeckler, F. M. Using Surface Scans for the Evaluation of Halogen Bonds toward the Side Chains of Aspartate, Asparagine, Glutamate, and Glutamine. J. Chem. Inf. Model. 2016, 56, 1373-1383. 14. Zimmermann, M. O.; Lange, A.; Wilcken, R.; Cieslik, M. B.; Exner, T. E.; Joerger, A. C.; Koch, P.; Boeckler, F. M. Halogenenriched fragment libraries as chemical probes for harnessing halogen bonding in fragment-based lead discovery. Future Med. Chem. 2014, 6, 617-639. 15. Kolar, M.; Hostas, J.; Hobza, P. The strength and directionality of a halogen bond are co-determined by the magnitude and size of the [sigma]-hole. PCCP 2014, 16, 99879996. 16. Riley, K. E.; Murray, J. S.; Fanfrlik, J.; Rezac, J.; Sola, R. J.; Concha, M. C.; Ramos, F. M.; Politzer, P. Halogen bond tunability I: the effects of aromatic fluorine substitution on the strengths of halogen-bonding interactions involving chlorine, bromine, and iodine. J. Mol. Model. 2011, 17, 3309-3318. 17. Riley, K. E.; Murray, J. S.; Fanfrlik, J.; Rezac, J.; Sola, R. J.; Concha, M. C.; Ramos, F. M.; Politzer, P. Halogen bond tunability II: the varying roles of electrostatic and dispersion contributions to attraction in halogen bonds. J. Mol. Model. 2013, 19, 4651-4659. 18. Bissantz, C.; Kuhn, B.; Stahl, M. A medicinal chemist's guide to molecular interactions. J. Med. Chem. 2010, 53, 5061-84. 19. Varadwaj, A.; Varadwaj, P. R.; Jin, B.-Y. Fluorines in tetrafluoromethane as halogen bond donors: Revisiting address the nature of the fluorine's σhole. Int. J. Quantum Chem 2015, 115, 453-470. 20. Fanfrlík, J.; Kolář, M.; Kamlar, M.; Hurný, D.; Ruiz, F. X.; Cousido-Siah, A.; Mitschler, A.; Řezáč, J.; Munusamy, E.; Lepšík, M.; Matějíček, P.; Veselý, J.; Podjarny, A.; Hobza, P. Modulation of Aldose Reductase Inhibition by Halogen Bond Tuning. ACS Chem. Biol. 2013, 8, 2484-2492. 21. Kolar, M. H.; Carloni, P.; Hobza, P. Statistical analysis of [sigma]-holes: a novel complementary view on halogen bonding. PCCP 2014, 16, 19111-19114.

Page 10 of 14

22. Heidrich, J.; Exner, T. E.; Boeckler, F. M. Predicting the Magnitude of σ-holes Using VmaxPred, a Fast and Efficient Tool Supporting the Application of Halogen Bonds in Drug Discovery. Unpublished. 23. Riley, K. E.; Murray, J. S.; Politzer, P.; Concha, M. C.; Hobza, P. Br···O Complexes as Probes of Factors Affecting Halogen Bonding: Interactions of Bromobenzenes and Bromopyrimidines with Acetone. J. Chem. Theory Comput. 2009, 5, 155-163. 24. Kozuch, S.; Martin, J. M. L. Halogen Bonds: Benchmarks and Theoretical Analysis. J. Chem. Theory Comput. 2013, 9, 19181931. 25. Anderson, L. N.; Aquino, F. W.; Raeber, A. E.; Chen, X.; Wong, B. M. Halogen Bonding Interactions: Revised Benchmarks and a New Assessment of Exchange vs Dispersion. J. Chem. Theory Comput. 2018, 14, 180-190. 26. Boys, S. F.; Bernardi, F. The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. Mol. Phys. 1970, 19, 553-566. 27. Voth, A. R.; Khuu, P.; Oishi, K.; Ho, P. S. Halogen bonds as orthogonal molecular interactions to hydrogen bonds. Nat. Chem. 2009, 1, 74-79. 28. Hattig, C. Optimization of auxiliary basis sets for RIMP2 and RI-CC2 calculations: Core-valence and quintuple-[small zeta] basis sets for H to Ar and QZVPP basis sets for Li to Kr. PCCP 2005, 7, 59-66. 29. Team, R. C. R: A Language and Enviroment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria: 2014. 30. Brickmann, J., Exner, T. E., Gimmler, J., Lautenschläger, P., Heiden, W., Moeckel, G., and Zahn, D. MOLCAD II, V1.4; MOLCAD GmbH: Darmstadt, Germany. 31. Brickmann, J.; Exner, T. E.; Keil, M.; Marhofer, R. J. Molecular graphics - Trends and perspectives. J. Mol. Model. 2000, 6, 328-340. 32. Politzer, P.; Murray, J. S.; Clark, T.; Resnati, G. The sigma-hole revisited. PCCP 2017, 19, 32166-32178. 33. Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic structure calculations on workstation computers: The program system turbomole. Chem. Phys. Lett. 1989, 162, 165-169. 34. TURBOMOLE V6.4 2012, a Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 19892007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. 35. TURBOMOLE V7.1 2016, a Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. 36. TURBOMOLE V7.2 2017, a Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 19892007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. 37. Hobza, P.; Šponer, J. Toward True DNA Base-Stacking Energies: MP2, CCSD(T), and Complete Basis Set Calculations. J. Am. Chem. Soc. 2002, 124, 11802-11808. 38. Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. PCCP 2005, 7, 32973305. 39. Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of approximate integrals in ab initio theory. An application in MP2 energy calculations. Chem. Phys. Lett. 1993, 208, 359-363. 40. Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: optimized auxiliary basis sets and demonstration of efficiency. Chem. Phys. Lett. 1998, 294, 143-152.


Page 11 of 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




Page 12 of 14

For Table of Contents Only:

ACS Paragon Plus Environment

12

Page 13 of 14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Overview of the workflow applied for the systematic evaluation of σ-hole tuning effects. (a) represents the systematic generation of halogenated scaffolds, illustrated using pyrazolo[1,5-a]pyridine. The resulting electrostatic features around the halogen are shown for 2-iodo-pyrazolo[1,5-a]pyridine using a 2D (b) and 3D grid (c). The most positive area on the isodensity surface of the σ-hole, VS,max, is automatically selected and highlighted in (d). Adduct formation energies are calculated on a MP2/TZVPP level of theory (e) and correlations between VS,max and these interaction energies are derived (f). 84x56mm (600 x 600 DPI)



88x34mm (600 x 600 DPI)


Page 14 of 14

Scaffold Effects on Halogen Bonding Strength - Journal of Chemical

Recommend Documents