Biomolecular Simulations of Halogen Bonds with a GROMOS Force

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Biomolecular Simulations of Halogen Bonds with a GROMOS Force Field Rafael Nunes,†,‡,§ Diogo Vila-Viçosa,†,‡ Miguel Machuqueiro,*,†,‡ and Paulo J. Costa*,†,‡ †

Downloaded via UNIV OF SUNDERLAND on September 29, 2018 at 17:26:34 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Centro de Química e Bioquímica, Departamento de Química e Bioquímica, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal ‡ BioISI - Biosystems & Integrative Sciences Institute, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, C8 bdg, 1749-016 Lisboa, Portugal § Centro de Química Estrutural, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal S Supporting Information *

ABSTRACT: Halogen bonds (XBs) are non-covalent interactions in which halogens (X), acting as electrophiles, interact with Lewis bases. XBs are able to mediate protein−ligand recognition and therefore play an important role in rational drug design. In this context, the development of molecular modeling tools that can tackle XBs is paramount. XBs are predominantly explained by the existence of a positive region on the electrostatic potential of X named the σ-hole. Typically, with molecular mechanics force fields, this region is modeled using a charged extra point (EP) linked to X along the R−X covalent bond axis. In this work, we developed the first EP-based strategy for GROMOS force fields (specifically GROMOS 54A7) using bacteriophage T4 lysozyme in complex with both iodobenzene and iodopentafluorobenzene as a prototype system. Several EP parametrization schemes were tested by adding a virtual interaction site to ligand topologies retrieved from the Automated Topology Builder (ATB) and Repository. Contrary to previous approaches using other force fields, our analysis is based on the capability of each parametrization scheme to sample XBs during MD simulations. Our results indicate that the implementation of an EP at a distance from iodine corresponding to Rmin provides a good qualitative description of XBs in MD simulations, supporting the compatibility of our approach with the GROMOS 54A7 force field.

1. INTRODUCTION Traditionally, halogenation of molecules has occupied a prominent role in drug discovery, mostly as a strategy for tuning drug-like properties such as membrane permeability and pharmacokinetic stability.1,2 It has been also recognized that halogenated molecules can directly bind biomolecular targets3,4 via a specific and directional non-covalent interaction, R−X···B (X = Cl, Br, I; B = Lewis base; R = substituent), known as a halogen bond (XB).5 This interaction arises from the anisotropic distribution of electron density in covalently bound halogens (except for fluorine because of its weak polarizability), which leads to the existence of an electrophilic site at their outermost region, known as the σ-hole.6 Their distinctive features have encouraged their widespread application,7,8 among others, in supramolecular chemistry9 and anion recognition,7,10,11 where some of us have been involved.12−15 Since the initial account of the relevance of XBs in biological systems,3 subsequent comprehensive surveys of Protein Data Bank (PDB) structures showed the occurrence of structurally diversified biomolecular XBs.4,16−19 These feature a wide range of acceptors from both the main chain and side chains of proteins, such as sulfur, oxygen, and nitrogen atoms and π systems, although the backbone carbonyl oxygens are typically involved and the side chains are notably under-represented. © XXXX American Chemical Society

Recently, this odd feature was explained by the fact that force fields (FFs) used to refine the structures cannot handle XBs, and therefore, those interactions with the more flexible side chains are lost.20 Curiously, the same type of analysis showed that backbone XBs are also overlooked, with XBs being lost in going from high- to low-resolution structures.21 XBs have been increasingly exploited as tools for rational drug design2,17,19,22,23 and biomolecular engineering,24−26 often modulating drug recognition. Indeed, heavier halogens persist throughout drug discovery and development, while fluorinated leads are more readily eliminated, even though the latter are dominant in earlier stages of the process.4 As for X-ray structure refinement,20,21 XB-based strategies in medicinal chemistry require an appropriate representation of XBs by computational or molecular modeling methodologies.2,27 In this scope, FFbased methods are privileged tools to study large biomolecular systems, but they are typically unable to account for the anisotropic charge distribution displayed by halogens. To circumvent this issue (and excluding polarizable FFs or similar approaches),28−30 a common strategy relies on a positively charged extra point (EP) placed at a defined distance from the Received: March 19, 2018 Published: September 14, 2018 A

DOI: 10.1021/acs.jctc.8b00278 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. (a) X-ray structure of bacteriophage T4 lysozyme L99A in complex with 5fibz (PDB ID 3DN3), featuring a very short I···S halogen bond with the side chain of Met102 within the hydrophobic cavity of the protein. (b, c) The crystal structure of the complex with ibz (PDB ID 3DN4) features two alternative binding modes: one (b) featuring an I···S halogen bond with the side chain of Met102, while in the alternative pose (c) the ligand is rotated by approximately 180°, leading to an I···O halogen bond with the main-chain carbonyl of Leu84 instead.56

In this work, we used a protein−ligand model to evaluate the effect of the EP in the overall network of intermolecular interactions captured throughout MD simulations. Moreover, we intended to evaluate the potential transferability of different EP parametrization schemes to a united-atom GROMOS FF, allowing the scope of XB simulations to be broadened to other FF families. For that purpose, we used bacteriophage T4 lysozyme (T4L) as a prototype system. This system has been extensively characterized and frequently employed to probe fundamental properties related to the structure, dynamics, and stability of proteins52,53 and ligand binding,54 some of them in the context of XBs.24 We focused on the Leu99 → Ala (L99A) mutant of phage T4L, which exhibits an unusually large internal cavity within the hydrophobic core of the enzyme (∼150 Å3; Figure S1)55 and is known to accommodate apolar molecules such as benzene and its halogenated derivatives,56 rendering it a suitable model to study ligand−protein interactions.54 Indeed, X-ray crystallographic data obtained for T4L L99A complexed with molecules of the C6F5X (X = H, F, Cl, Br, I) or C6H5X (X = H, I) family show that iodinated derivatives bind to the protein via I···S XB interactions. The structure containing iodopentafluorobenzene (5fibz), shown in Figure 1a, features a halogen bond between the iodine atom of the ligand and the sulfur atom in the side chain of Met102. In contrast, the weaker XB-donor iodobenzene (ibz) interacts with the protein by two distinct orientations (see Figure 1b,c), thus showing that this system provides an appropriate case study to be tackled by MD simulations.

halogen atom along the R−X covalent bond axis in order to explicitly represent the σ-hole,31−33 thereby establishing the directionality and energetics of the XB. In the first implementation, reported by Ibrahim31 using the General Amber Force Field (GAFF),34 atomic partial charges (including that of the EP) were determined by fitting the molecular electrostatic potential (ESP) according to the restrained ESP (RESP) procedure35 at a given X···EP distance. Optimal X···EP distances, corresponding to the atomic radius of X derived from the Rmin Lennard-Jones (LJ) parameter, were determined by comparing the resulting XB lengths with those from density functional theory (DFT) calculations for a series of model halobenzene···formaldehyde dimers. Almost simultaneously, a very similar scheme was proposed by Sironi and co-workers32 in which the optimal X···EP distances were assigned according to the best fitting to the quantum-mechanical (QM) ESP generated by RESP charges, meaning that each molecule has a corresponding X···EP distance. Later, the group of Hobza33 suggested the use of an universal set of EP charges and X···EP distances to model brominated compounds that was later expanded to include iodine and chlorine in the context of a docking implementation.36 In this model, the charges of both the EP and X are assigned on the basis of previously calculated ones, partitioning the charge of X into those two contributions. Although this strategy offers numerous advantages in terms of assigning charges for large data sets of compounds, the performance of the model was inferior to that obtained by fitting all of the charges as in the above-mentioned models.31,32 Even though related implementations exist for the OPLS37,38 and CHARMM39 FF families, the aforementioned strategies developed in the framework of AMBER-flavored FFs, particularly the one from Ibrahim,31 are more common in molecular dynamics (MD) simulations of XB-mediated enzyme−inhibitor or receptor−(ant)agonist associations.40−45 In spite of these recent advances, to the best of our knowledge a comparative analysis of the performance of distinct EP-based methods in long MD simulations has not been performed to date. Additionally, none of the reported studies used the popular united-atom GROMOS FFs.46−49 The use of these FFs might be hindered by the fact that drug-like molecules possess a wide variety of chemical groups and thus are not trivial to parametrize. The recently developed Automated Topology Builder (ATB) platform50,51 tackles this issue, allowing the user to extract GROMOS-compatible topologies for virtually any drug-like molecule. However, if the molecule is halogenated, no specific term for describing XBs is added.

2. METHODS 2.1. Ligand Parametrization. We selected iodobenzene (ibz) and iodopentafluorobenzene (5f ibz) as probe molecules, representing mild and strong XB donors, respectively. Topologies were taken from the ATB repository50,51 and were subsequently curated by deriving new sets of atomic partial charges as described below. All of the other parameters (bonded and nonbonded) were kept as provided in the original ATB topology files with the exception of the atomic mass of iodine, which was corrected to 126.9044 u. A new set of atomic charges was calculated following a procedure similar to the one described for ATB that uses GAMESS-US57 for all of the QM calculations. Herein we used Gaussian 0958 to optimize the molecular geometries with the B3LYP functional59,60 and the 6-31G(d) basis set61,62 for all atoms, with the exception of iodine, for which the 6-311G(d) basis set63 was used. The calculations were performed in an implicit solvent (water) using the C-PCM64,65 variant of the B


Article

Journal of Chemical Theory and Computation polarizable continuum model (PCM)66 and the set of Bondi atomic radii.67 Frequency calculations were carried out to ensure that the optimized geometries corresponded to minima. ESPs were obtained at the same level of theory. In those calculations, default Merz−Kollman (MK) atomic radii were employed for all elements except iodine, for which an atomic radius of 2.3 Å was used.68 Atomic partial charges were then generated according to the RESP35 procedure using the Antechamber module as implemented in AmberTools14.69 A massless particle (EP) was added during the RESP charge fitting procedure at a variable distance from iodine along the C−I covalent bond axis direction, with the C−I···EP angle being fixed at 180.0°. For each halobenzene ligand, three I···EP distances were tested: a distance obtained using the procedure of Ibrahim,31 a distance found using the procedure of Sironi and co-workers,32 and an intermediate distance. For the Ibrahim model, we used an Rmin value of 2.16 Å, derived from the C6 and C12 LJ parameters for iodine provided by ATB, yielding ligand models ibz2.16 and 5fibz2.16. It should be noted that this Rmin value is very similar to that employed in the latest versions of GAFF (2.15 Å).34 In the case of Sironi’s approach, I···EP distances were determined as those optimizing the quality of the fit to the QM ESP based on the relative root-mean-square (RRMS) error,35 which was achieved at 1.15 and 1.88 Å for ibz and 5f ibz, respectively, yielding ibz1.15 and 5f ibz1.88. The third charge scheme corresponds to an I···EP distance intermediate between those from the other two methods, i.e., 1.65 Å for ibz (ibz1.65) and 2.02 Å for 5f ibz (5f ibz2.02). Models without the addition of an EP (ibznoEP and 5fibznoEP) were also tested. All of the models are summarized in Table 1, and the full sets of charges are provided

Water molecules as well as cocrystallized phosphate ions and 2hydroxyethyl disulfide molecules present in the original structure were removed. All of the side chains in the protein were simulated with the typical protonation states at pH 7. The two alternative binding modes observed in the ibz−L99A complex (see Figure 1b,c) were used as starting conformations for distinct simulation replicates for both ibz and 5fibz. Each ibz−L99A or 5f ibz−L99A complex was solvated with 8540 single point charge (SPC)74 water molecules, and eight chloride anions were added to neutralize the net charge of each system. A rhombic-dodecahedral simulation box was used, applying periodic boundary conditions in all directions with the minimum image convention. All of the simulations were performed in the isothermal−isobaric (NPT) ensemble. The temperature was kept constant at 300 K by separately coupling the solute and solvent (including anions) to external thermostats using the velocity-rescale algorithm75 with a coupling constant of 0.1 ps, while an isotropic Parrinello−Rahman barostat76,77 with a coupling constant of 2.0 ps and an isothermal compressibility of 4.5 × 10−5 bar−1 was used to maintain the pressure at 1 bar. Long-range electrostatic interactions were treated using the smooth particle mesh Ewald (PME) method78,79 with a Fourier grid spacing of 0.12 nm and a cutoff of 1.0 nm for direct contributions. Lennard-Jones interactions were treated using a nonbonded neighbor pair list with a cutoff of 1.0 nm, allowing the use of the more efficient GPU implementation of the buffered Verlet list scheme,80 which requires the same cutoff value for all nonbonded interactions. Although the LJ cutoff is different from the original implementation of the GROMOS 54A7 FF (1.4 nm), this should have only a negligible effect since the LJ interactions of isotropic systems usually converge at a very short cutoff radius. All of the protein and ligand bonds were constrained using the parallel version of the linear constraint solver (P-LINCS)81 algorithm, while SETTLE82 was used to constrain water molecules. Energy minimization was carried out in two steps using the steepest-descent algorithm until machine precision was reached: an initial minimization was performed without constraints, followed by another minimization with all of the bond lengths constrained. Initialization was performed by means of a 100 ps simulation with the positions of all atoms restrained using a force constant of 1000 kJ nm−2 mol−1, followed by a second 100 ps simulation with the protein Cα and ligand atoms restrained (1000 kJ nm−2 mol−1). Initial velocities were randomly taken from a Maxwell−Boltzmann distribution at 300 K in the first initialization step. For each system, and using both starting conformations, a total of five simulation replicates were run for 100 ns. The equations of motion were numerically integrated using the leapfrog algorithm with a time step of 2 fs. For analysis purposes, conformations were saved every 10 ps, and only the last 80 ns of each trajectory was considered. 2.3. Analysis. We computed the distances between the iodine atom of the ligand and all of the potential XB acceptors (A) in the protein (dI···A) along with the respective C−I···A bond angles (θ), as shown in Figure 2a. For each configuration, only one acceptor was selected, corresponding to the shortest of all possible I···A distances. We considered as XB acceptors all oxygen atoms from backbone carbonyl groups83 and side-chain hydroxy (Tyr, Ser, Thr), carboxylate (Asp, Glu), and carboxamide (Asn, Gln)84 groups as well as sulfur atoms from Met residues.85 Nitrogen atoms from both the protein backbone and amino acid side chains (Gln, Asn, His, Trp) have often been addressed as potential acceptors in biomolecular XBs21,84,86,87

Table 1. Summary of the Halobenzene Ligands Parametrized and the Relevant Equilibrium Distances and Charges for Each Parametrization Schemea ligand

rI···EP/Å

qC/e

qI/e

qEP/e

ibznoEP ibz1.15 ibz1.65 ibz2.16 5f ibznoEP 5f ibz1.88 5f ibz2.02 5f ibz2.16

− 1.15 1.65 2.16 − 1.88 2.02 2.16

−0.262 +0.290 +0.207 +0.115 −0.645 −0.184 −0.205 −0.228

−0.072 −0.536 −0.384 −0.288 +0.122 −0.132 −0.110 −0.090

− +0.241 +0.132 +0.076 − +0.093 +0.080 +0.070

a

Partial charges are shown only for iodine (qI), its covalently bonded carbon (qC), and (when applicable) the EP (qEP). For the full set of charges, see Figure S2.

in Figure S2. The EPs and charges were incorporated into the respective ATB topology files by means of a virtual interaction site defined by the C−I bond, as implemented in GROMACS.70−72 It should be noted that the parametrization schemes reported in the literature31,32 are herein reproduced in a modified version since the EP is no longer represented explicitly as a dummy particle but instead is represented by a massless virtual site without LJ parameters. 2.2. MM/MD Settings. Molecular mechanics/molecular dynamics (MM/MD) simulations were carried out using the GROMACS software package, version 2016.1.70−72 The GROMOS 54A7 force field47 was used for the protein, together with the modified ATB ligand topologies as described above. The crystal structure of bacteriophage T4L L99A in complex with ibz was retrieved from the PDB73 (PDB ID 3DN456). C


Article


the corresponding C−I···COM angle (θ). In the case of the indole group from Trp residues, the two rings forming the bicyclic structure were treated separately, i.e., as two distinct acceptors, for simplicity. It should be noted that in the literature,4,16 besides the C−X···COM angle (θ), the angle α formed by X, the COM of the aromatic ring, and the vector normal to the ring plane passing through the centroid (Figure 2b) is also employed as a criterion for the identification of C− X···π XBs. However, careful inspection of the data from our simulations showed that this additional criterion is met in practically all conformations of interest (see Interaction Preferences by Acceptor Type in the Supporting Information) and thus, for the purpose of our analysis, the geometrical description could be simplified to the two-coordinate (dI···COM, θ) space. From a total of 451 acceptors in the protein, only 84 were effectively engaged as the nearest acceptor, encompassing the more frequent oxygens (49) and nitrogens (28) but also sulfurs (3) and aromatic π systems (4). Probability density surfaces and free energy landscapes were used to describe the two-dimensional (2D) space generated by

Figure 2. Geometrical criteria for defining (a) C−I···A (A = O, S, N) and (b) C−I···π XBs. In this work, the geometrical parameter α was not considered, thereby allowing a two-dimensional description of C−I···π interactions (see the text).

and therefore were also screened, with the exception of those positively charged (Arg and Lys). π systems from Phe, Tyr, Trp, and His were also included by computing the distance between iodine and the center of mass (COM) of the aromatic ring and

Figure 3. Free energy landscapes obtained from the simulations of L99A T4 lysozyme in complex with the ibz (left panel) and 5f ibz (right panel) ligands parametrized with different EP models or without EP addition, generated using the shortest I···A distance (d) and the C−I···A angle (θ) for all possible acceptors. The dash-limited regions identify typical halogen bonds. D


Article


charge variation and the existence of a positively charged iodine atom in 5f ibznoEP show the usefulness of these models as probes for MD simulations, as we will discuss in the following section. 3.2. ibz/5f ibz−L99A Intermolecular Interactions. The structure of T4L is mainly helical (>60% α-helix content), comprising two domains connected by a long helix, and the hydrophobic cavity of the L99A mutant is located within the Cterminal domain.52 Our simulations captured a well-characterized91,92 large-amplitude hinge-bending motion between opened and closed states associated with substrate binding. Nevertheless, the hydrophobic pocket does not undergo any major structural variation, with the ibz/5f ibz ligands freely rotating inside this pocket, allowing the use of the two crystallographic ligand binding modes as initial configurations. The network of intermolecular interactions established by the halobenzene ligands was represented by a two-dimensional landscape (Figures 3 and S4) as detailed in Methods. Although it does not provide a full description of the configurational space, this approach was devised with the aim of capturing not only XBs but also other potential interactions involving iodine and the range of acceptors considered, namely, orthogonal hydrogen bonds (HBs) as well as other unspecific intermolecular contacts. The landscapes obtained in the absence of an EP (Figure 3, top panel) are quite homogeneous, exhibiting a single free energy minimum centered at relatively long I···A distances (∼4 Å) and interaction angles of ∼140° and ∼110° for ibznoEP and 5f ibznoEP, respectively. This indicates that XBs are not sampled in these simulations, since these are typically characterized by shorter I··· A distances and larger angles (dI···A < 3.8 Å and θ > 140°; see Methods), with the exception of those involving aromatic acceptors (dI···A < 4.5 Å and θ > 120°), for which the contribution is low (ca. 5%). While inability of the models without an EP to describe XBs is not surprising, the free energy landscape obtained for 5f ibznoEP shows that a positively charged iodine or the existence of stable MM-optimized halogen-bonded dimers is not sufficient for the description of XBs in MD simulations. In turn, when this ligand is modeled with an EP, more scattered surfaces are obtained where several minima can be clearly distinguished, including in the region corresponding to XB interactions. In fact, as the I···EP distance increases, the XB region becomes more favorable while other local minima are also visible, indicating that alternative intermolecular interactions are sampled as well. In regard to EP-based ibz simulations, XBs are mostly not captured when using shorter I···EP distances (ibz1.15, ibz1.65), whereas the use of a larger distance (ibz2.16) again enables sampling at the target XB region. Hence, in 5fibz2.16 and ibz2.16 the resulting configurational sampling is consistent with ibz being a relatively weak halogenbonding ligand compared with 5f ibz and therefore is more physically realistic. In contrast, for the ibz1.15 and 5f ibz1.88 models, a poorer description of XB interactions was obtained, as no XBs were sampled for ibz. This disagrees with the X-ray structure, indicating that minimization of the error associated with the RESP fitting is not an appropriate criterion for assigning optimal I···EP distances and charges in the GROMOS 54A7 force field. In the range of I···EP distances considered in our study (1.15−2.16 Å and 1.88−2.16 Å for ibz and 5f ibz, respectively), the RRMS variation is quite small (values in Figure S2), i.e., even though the difference in the errors is negligible, the resulting sets of charges may vary considerably (Table 1), ultimately leading to considerably distinct scenarios regarding the accessible configurational space.

the calculated pairs of geometrical coordinates. The probability density functions were estimated from the simulation data using a Gaussian kernel,88 and the resulting probability density surfaces were converted (when necessary) into free energy surfaces according to E(r) = −RT ln

P(r) Pmax

(1)

where r is the coordinate along the 2D space and Pmax is the maximum value of the probability density function P(r).89 The presence of XB interactions was evaluated according to geometrical criteria commonly employed for XB assignment in the context of crystallographic studies,4,16,20,21 namely, a C−I··· A angle (θ) larger than 140° and an I···A distance (dI···A) shorter than the sum of the respective van der Waals radii (dI···O < 3.50 Å, dI···S < 3.78 Å, and dI···N < 3.53 Å). For aromatic acceptors, an I···COM distance (dI···COM) less than 4.5 Å along with a C−I··· COM angle (θ) larger than 120° was employed. In order to exclude potential false positives pertaining to hydrogen bonds formed with side-chain hydroxy groups (Tyr, Ser, Thr), conformations with dI···H < dI···O were discarded. All of the reported error bars in the figures correspond to the standard error of the mean of the replicates.

3. RESULTS AND DISCUSSION 3.1. Molecular Mechanical Description of Model Systems. Before carrying out biomolecular simulations, we checked the ability of the GROMOS 54A7 force field and the derived charge sets to predict halogen-bonded geometries with and without the addition of an EP. In line with previous reports,31,33,90 we performed gas-phase MM calculations for a set of dimers formed between ibz/5f ibz and several Lewis bases that are representative mimics of XB-acceptor motifs present in proteins (Figure S3). These include acetone and acetamide (backbone, Gln, Asn), methanol (Ser, Thr, Tyr), dimethyl sulfide (Met), imidazole (His), and toluene (Phe, Tyr, Trp, His). As expected, the incorporation of an EP improves the MM description of halogen bonds since for all of the models a halogen-bonded local minimum is obtained (Table S1). Interestingly, halogen-bonded dimers can also be obtained in the absence of an EP for the specific case of 5f ibznoEP. This arises from the positive charge derived for iodine (see Table 1) due to the electron-withdrawing effect of the fluorine substituents, contrasting with the negatively charged iodine in ibznoEP. Additionally, as the I···EP distance increases for each molecule, the charge of the EP becomes less positive (Table 1). This invariably leads to shorter I···A distances (and typically larger interaction energies) closer to those obtained from QM calculations (data not shown), which might appear counterintuitive. However, it should be noted that the decrease in positive charge of the EP with increasing I···EP distance is concomitant with the iodine becoming less negatively charged. The magnitude of this charge appears to have a predominant role in the XB interaction distance and strength. Systems containing toluene are an exception since all of the calculated dimers display rather short I···A distances compared with the defined criteria (see Methods). The interaction energies are significantly larger for 5f ibz compared with ibz, as expected. In summary, the inclusion of an EP leads to a more correct description of XBs in MM-minimized dimers, in agreement with previous studies for other force fields. The mentioned EP/I E


Article


nitrogen in the 5fibz simulations, which have frequencies above and below their abundances, respectively. These frequencies should not be confused with the frequency associated with each acceptor type being effectively engaged in halogen bonding (see below). The observed variability of the nearest acceptor type cannot be dissociated from the spatial distribution of each acceptor around the ligand, which eventually determines its distinct propensity for engaging in halogen bonding. In order to gain qualitative insight into such preferences, we generated a spatial distribution of the nearest acceptor around iodine obtained from our MD simulations for each model (Figure 5). Visual inspection of that spatial distribution, particularly the top-view representations (right panels), shows that the region located head-on to iodine along the C−I bond axis is more efficiently sampled when a larger I···EP distance is employed; this effect is more clear in the sequence 5f ibznoEP → 5f ibz2.16. Indeed, while for 5f ibznoEP only a negligible number of acceptors are observed in this region, as an EP is added at incremental I···EP distances, the same area becomes increasingly populated, mostly with oxygen atoms. On the other hand, sulfur acceptors apparently prefer to increasingly occupy the vicinity of the aromatic ring in the 5fibz models with an EP. 3.4. Halogen Bond Frequency Evaluation. The number of I···A contacts represented in the free energy landscapes that effectively correspond to XB interactions can be quantified by applying the dI···A/θ geometrical criteria of each acceptor type (see Methods). These criteria are compatible with the interaction distances and angles that characterize XB sampling in our GROMOS simulations. Figure 6 shows the frequencies of XB interactions in the ibz and 5f ibz simulations, decomposed by backbone and side-chain acceptors (a further decomposition by acceptor subtype is shown in Figure S6). Very small numbers of XBs are captured in the simulations without addition of an EP (∼0.7% and ∼2.7% for ibznoEP and 5f ibznoEP, respectively). In contrast, when the ligand is modeled with an EP, increased numbers of interactions are obtained, with the results depending on the specific ligand/model system. 5f ibz simulations invariably yield more XB occurrences (9.9−37.2%) compared with ibz (1.7−4.5%), in line with their distinct XB donor propensities, as previously discussed. The use of the larger X··· EP distances results in a more representative and balanced sampling of the target interactions. Indeed, XB interactions account for ∼4.5% and ∼37.2% of the I···A contacts in ibz2.16 and 5f ibz2.16 simulations, respectively. A pronounced prevalence is observed for XB interactions with backbone acceptors (∼85.2%), in agreement with recently reported database surveys.4,20 The contribution of each acceptor type to the total number of XBs captured in the EP-based simulations is shown in Figure S7. XB preferences vary depending on the model but the most important feature is that the previously observed nearest acceptor preferences (see Figure 4) are not necessarily associated with their respective XB acceptor propensities. For instance, nitrogen acceptors, which were previously shown to frequently act as the nearest acceptor, are only marginally engaged in XB interactions. We further highlight that in the case of 5f ibz2.16, where oxygen acceptors account for ∼95% of the XBs, a significant number of residues (28) are involved, showing that the observed preferences are not driven by sampling limitations.

The presented free energy landscapes (Figure 3) can be deconvoluted into contributions pertaining to each acceptor type, thus allowing their individual geometrical preferences to be discriminated. This analysis together with a thorough discussion is presented in Interaction Preferences by Acceptor Type in the Supporting Information. Overall, the deconvolution by acceptor supports the better performance of the 5f ibz2.16 and ibz2.16 models in describing not only XBs but also the geometrically orthogonal HBs.93,94 Altogether, these results pinpoint the use of longer X···EP distances, in particular those corresponding to the Rmin value of the halogen,31 as a suitable strategy for performing biomolecular simulations of XBs using GROMOS 54A7. As we have already discussed, increasing the I···EP distance leads to a less positively charged EP but also to a less negatively charged iodine (Table 1). This effect seems to modulate the depiction of XBs in the simulations, since for each ligand only the charges of the EP, the iodine atom, and its covalently bonded carbon vary significantly among models (see Figure S2). We did not observe any instabilities during the MD simulations, even when larger I···EP distances were employed (ibz2.16, 5f ibz2.16). Numerical instabilities occurred only for I···EP distances > Rmin as a result of clashes between the EP and the acceptor (results not shown). Therefore, the Rmin value, calculated from the interaction between two iodine atoms, provides a good estimation of the maximum acceptable I···EP distance. It should be noted that Kolár ̌ and Hobza have proposed that the inclusion of a repulsive LJ parameter on the EP is most probably required to guarantee the stability of the simulations when the EP is introduced at longer X···EP distances, although the authors also argued that the relevance of this effect should be further investigated.27,33 3.3. Nearest Acceptor Distribution. We quantified the frequency with which each acceptor type effectively contributes to the surface representations, i.e., the nearest acceptor at each time frame (Figure 4). The frequencies are slightly different

Figure 4. Frequencies (in %) of the nearest acceptor types averaged over the simulation trajectories for each system/parametrization scheme.

among different systems and/or parametrization models. Oxygen atoms are invariably the most prevalent acceptors, in agreement with their superior occurrence in the protein, except for 5f ibznoEP, though within the error bars. Figure S5 shows the frequency of each acceptor weighted by the respective abundance in the protein (relative frequency). The relative frequency is ∼1 for all of the acceptors, except for sulfur and F


Article


Figure 5. Representations of the spatial distributions of the nearest acceptor around iodine obtained from the simulations of L99A T4 lysozyme in complex with the (left) ibz or (right) 5f ibz ligands parametrized with different models. The acceptor types are colored as follows: oxygen in red, sulfur in yellow, nitrogen in blue, and π systems (represented by the corresponding center of mass) in black. For each ligand system, a side view is shown at the left and a top view at the right. All of the rendered structures were generated using a total of 800 equally spaced conformations from each simulation trajectory, which were superimposed by RMS fitting of all ligand ring carbon atoms. For visualization purposes, only one reference structure of each ligand is shown for clarity.

4. CONCLUSIONS A proper description of halogen bonds is mandatory to correctly model halogenated drugs in MD simulations. We have implemented several models to describe this intermolecular interaction via the addition of an extra point (EP) to an iodine atom. The analysis of their respective free energy surfaces has clarified relevant requirements for the implementation of this approach toward XB modeling in MD simulations using GROMOS force fields. Although we specifically used GROMOS 54A7,47 we expect this strategy to hold for the related 53A646 and 54A848,49 force fields. The presence of a positively charged halogen atom is not sufficient for an accurate description of XBs in MD simulations, notwithstanding the stability of the corresponding MM-optimized structure. In contrast, when an EP is added, XBs are more sampled, but different X···EP distances/charge sets yield significantly different free energy surfaces. Only when the EP is introduced at Rmin does the relative XB strength of each halobenzene ligand seem to be qualitatively reproduced. With the inclusion of an EP, particularly using this setup, XBs and other relevant interactions

Figure 6. Frequencies (in %) of XB interactions averaged over the total number of I···A contacts in the simulations for each system/ parametrization scheme. XB assignments were based on different geometrical criteria according to each acceptor type (see Methods).

G


Article


Regional de Lisboa (Lisboa 2020), Portugal 2020, FEDER/FN, and the European Union under Project 28455 (LISBOA-010145-FEDER-028455, PTDC/QUI-QFI/28455/2017).

(in particular HBs) are sampled and can be clearly distinguished by our representation of the configurational space. The total number of I···A XBs was quantified by applying typical crystallographic criteria. The number of XB interactions captured in the simulations varies among different ligand models, but the use of larger I···EP distances affords a more balanced sampling between the two XB donors. The distinct propensity of each acceptor type to engage in XB interactions was also characterized, and the backbone carbonyl oxygen was found to be the preferred interaction partner in the case of 5fibz2.16 simulations. Interactions involving aromatic π systems or sulfur acceptorsincluding those present in the crystallographic binding posesare also sampled, whereas nitrogen species do not contribute. Although the introduction of an EP at Rmin seems to provide a correct description of halogen bonding, we note that fine-tuning of the LJ parameters of the halogen atom and, eventually, of the EP might lead to improved results, as noted by Kolár ̌ and Hobza.33 Such an analysis is beyond the scope of this work, but nonetheless, herein we have provided the first comprehensive study of the behavior of EP-based methods in describing halogen bonds using MD simulation data, which can have significant implications for computer-aided drug discovery.

■

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Ricardo J. M. Rosa for fruitful discussions on the initial steps of the project.

■

(1) Gerebtzoff, G.; Li-Blatter, X.; Fischer, H.; Frentzel, A.; Seelig, A. Halogenation of drugs enhances membrane binding and permeation. ChemBioChem 2004, 5, 676−684. (2) Ford, M. C.; Ho, P. S. Computational tools to model halogen bonds in medicinal chemistry. J. Med. Chem. 2016, 59, 1655−1670. (3) 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. (4) Xu, Z.; Yang, Z.; Liu, Y.; Lu, Y.; Chen, K.; Zhu, W. Halogen bond: its role beyond drug-target binding affinity for drug discovery and development. J. Chem. Inf. Model. 2014, 54, 69−78. (5) 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−1713. (6) Clark, T.; Hennemann, M.; Murray, J. S.; Politzer, P. Halogen bonding: the σ-hole. J. Mol. Model. 2007, 13, 291−296. (7) Cavallo, G.; Metrangolo, P.; Milani, R.; Pilati, T.; Priimagi, A.; Resnati, G.; Terraneo, G. The halogen bond. Chem. Rev. 2016, 116, 2478−2601. (8) Costa, P. J. The halogen bond: nature and applications. Phys. Sci. Rev. 2017, 2, 20170136. (9) Gilday, L. C.; Robinson, S. W.; Barendt, T. A.; Langton, M. J.; Mullaney, B. R.; Beer, P. D. Halogen bonding in supramolecular chemistry. Chem. Rev. 2015, 115, 7118−7195. (10) Brown, A.; Beer, P. D. Halogen bonding anion recognition. Chem. Commun. 2016, 52, 8645−8658. (11) Beale, T. M.; Chudzinski, M. G.; Sarwar, M. G.; Taylor, M. S. Halogen bonding in solution: thermodynamics and applications. Chem. Soc. Rev. 2013, 42, 1667−1680. (12) Caballero, A.; Zapata, F.; White, N. G.; Costa, P. J.; Félix, V.; Beer, P. D. A halogen-bonding catenane for anion recognition and sensing. Angew. Chem., Int. Ed. 2012, 51, 1876−1880. (13) Zapata, F.; Caballero, A.; White, N. G.; Claridge, T. D. W.; Costa, P. J.; Félix, V.; Beer, P. D. Fluorescent charge-assisted halogen-bonding macrocyclic halo-imidazolium receptors for anion recognition and sensing in aqueous media. J. Am. Chem. Soc. 2012, 134, 11533−11541. (14) Gilday, L. C.; Lang, T.; Caballero, A.; Costa, P. J.; Félix, V.; Beer, P. D. A catenane assembled through a single charge-assisted halogen bond. Angew. Chem., Int. Ed. 2013, 52, 4356−4360. (15) Nunes, R.; Costa, P. J. Ion-pair halogen bonds in 2-halofunctionalized imidazolium chloride receptors: substituent and solvent effects. Chem. - Asian J. 2017, 12, 586−594. (16) Lu, Y.; Wang, Y.; Zhu, W. Nonbonding interactions of organic halogens in biological systems: implications for drug discovery and biomolecular design. Phys. Chem. Chem. Phys. 2010, 12, 4543−4551. (17) 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, 314−318. (18) Parisini, E.; Metrangolo, P.; Pilati, T.; Resnati, G.; Terraneo, G. Halogen bonding in halocarbon-protein complexes: a structural survey. Chem. Soc. Rev. 2011, 40, 2267−2278.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00278.

■

REFERENCES

Representation of the X-ray structure of the ligand-free cavity-containing L99A mutant of phage T4 lysozyme, atomic partial charges for all of the parametrized halobenzene ligands, representative geometries of the model complexes studied by gas-phase MM calculations, interaction energies and representative geometrical parameters for the calculated complexes, P(r) surfaces associated with all of the free energy landscapes, relative frequency of each nearest acceptor type, frequency of XB interactions split for all of the different acceptor subtypes, frequency of XB interactions by acceptor type, and results and discussion for the Interaction Preferences by Acceptor Type section (PDF) Topology files in ITP format for all of the parametrized halobenzene ligands (ZIP)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. Phone: +351-217500112. *E-mail: [email protected]. Phone: +351-21-7500845. ORCID

Rafael Nunes: 0000-0002-9014-0570 Miguel Machuqueiro: 0000-0001-6923-8744 Paulo J. Costa: 0000-0002-0492-6666 Funding

We thank Fundaçaõ para a Ciência e a Tecnologia (FCT), Portugal, for the Investigador FCT Program IF/00069/2014 and Exploratory Project IF/00069/2014/CP1216/CT0006 (P.J.C.), Project PTDC/QEQ-COM/5904/2014 (M.M.), Doctoral Grant SFRH/BD/116614/2016 (R.N.), and Strategic Projects UID/MULTI/00612/2013 and UID/MULTI/04046/ 2013. This work was also cofinanced by Programa Operacional H


Article

Journal of Chemical Theory and Computation (19) Scholfield, M. R.; Zanden, C. M. V.; Carter, M.; Ho, P. S. Halogen bonding (X-bonding): a biological perspective. Protein Sci. 2013, 22, 139−152. (20) Zhang, Q.; Xu, Z.; Zhu, W. The underestimated halogen bonds forming with protein side chains in drug discovery and design. J. Chem. Inf. Model. 2017, 57, 22−26. (21) Zhang, Q.; Xu, Z.; Shi, J.; Zhu, W. Underestimated halogen bonds forming with protein backbone in protein data bank. J. Chem. Inf. Model. 2017, 57, 1529−1534. (22) 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. (23) Lu, Y.; Shi, T.; Wang, Y.; Yang, H.; Yan, X.; Luo, X.; Jiang, H.; Zhu, W. Halogen bonding-a novel interaction for rational drug design? J. Med. Chem. 2009, 52, 2854−2862. (24) Scholfield, M. R.; Ford, M. C.; Carlsson, A.-C. C.; Butta, H.; Mehl, R. A.; Ho, P. S. Structure-energy relationships of halogen bonds in proteins. Biochemistry 2017, 56, 2794−2802. (25) El Hage, K.; Pandyarajan, V.; Phillips, N. B.; Smith, B. J.; Menting, J. G.; Whittaker, J.; Lawrence, M. C.; Meuwly, M.; Weiss, M. A. Extending halogen-based medicinal chemistry to proteins: iodoinsulin as a case study. J. Biol. Chem. 2016, 291, 27023−27041. (26) Carter, M.; Voth, A. R.; Scholfield, M. R.; Rummel, B.; Sowers, L. C.; Ho, P. S. Enthalpy-entropy compensation in biomolecular halogen bonds measured in DNA junctions. Biochemistry 2013, 52, 4891−4903. (27) Kolár,̌ M. H.; Hobza, P. Computer modeling of halogen bonds and other σ-hole interactions. Chem. Rev. 2016, 116, 5155−5187. (28) Mu, X.; Wang, Q.; Wang, L.-P.; Fried, S. D.; Piquemal, J.-P.; Dalby, K. N.; Ren, P. Modeling organochlorine compounds and the σhole effect using a polarizable multipole force field. J. Phys. Chem. B 2014, 118, 6456−6465. (29) Adluri, A. N. S.; Murphy, J. N.; Tozer, T.; Rowley, C. N. Polarizable force field with a σ-hole for liquid and aqueous bromomethane. J. Phys. Chem. B 2015, 119, 13422−13432. (30) Lin, F.-Y.; MacKerell, A. D., Jr. Polarizable empirical force field for halogen-containing compounds based on the classical drude oscillator. J. Chem. Theory Comput. 2018, 14, 1083−1098. (31) Ibrahim, M. A. A. Molecular mechanical study of halogen bonding in drug discovery. J. Comput. Chem. 2011, 32, 2564−2574. (32) Rendine, S.; Pieraccini, S.; Forni, A.; Sironi, M. Halogen bonding in ligand-receptor systems in the framework of classical force fields. Phys. Chem. Chem. Phys. 2011, 13, 19508−19516. (33) Kolár,̌ M.; Hobza, P. On extension of the current biomolecular empirical force field for the description of halogen bonds. J. Chem. Theory Comput. 2012, 8, 1325−1333. (34) 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. (35) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269−10280. (36) Kolár,̌ M.; Hobza, P.; Bronowska, A. K. Plugging the explicit σholes in molecular docking. Chem. Commun. 2013, 49, 981−983. (37) Jorgensen, W. L.; Schyman, P. Treatment of halogen bonding in the OPLS-AA force field: application to potent anti-HIV agents. J. Chem. Theory Comput. 2012, 8, 3895−3901. (38) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L.; Kaus, J. W.; Cerutti, D. S.; Krilov, G.; Jorgensen, W. L.; Abel, R.; Friesner, R. A. OPLS3: a force field providing broad coverage of drug-like small molecules and proteins. J. Chem. Theory Comput. 2016, 12, 281−296. (39) Soteras Gutiérrez, I.; Lin, F.-Y.; Vanommeslaeghe, K.; Lemkul, J. A.; Armacost, K. A.; Brooks, C. L., III; MacKerell, A. D., Jr. Parametrization of halogen bonds in the CHARMM general force field: improved treatment of ligand−protein interactions. Bioorg. Med. Chem. 2016, 24, 4812−4825.

(40) Celis-Barros, C.; Saavedra-Rivas, L.; Salgado, J. C.; Cassels, B. K.; Zapata-Torres, G. Molecular dynamics simulation of halogen bonding mimics experimental data for cathepsin L inhibition. J. Comput.-Aided Mol. Des. 2015, 29, 37−46. (41) Rosa, M.; Caltabiano, G.; Barreto-Valer, K.; Gonzalez-Nunez, V.; Gómez-Tamayo, J. C.; Ardá, A.; Jiménez-Barbero, J.; Pardo, L.; Rodríguez, R. E.; Arsequell, G.; Valencia, G. Modulation of the interaction between a peptide ligand and a G protein-coupled receptor by halogen atoms. ACS Med. Chem. Lett. 2015, 6, 872−876. (42) Luchi, A. M.; Angelina, E. L.; Andujar, S. A.; Enriz, R. D.; Peruchena, N. M. Halogen bonding in biological context: a computational study of D2 dopamine receptor. J. Phys. Org. Chem. 2016, 29, 645−655. (43) Rana, N.; Conley, J. M.; Soto-Velasquez, M.; Léon, F.; Cutler, S. J.; Watts, V. J.; Lill, M. A. Molecular modeling evaluation of the enantiomers of a novel adenylyl cyclase 2 inhibitor. J. Chem. Inf. Model. 2017, 57, 322−334. (44) González-Vera, J. A.; Medina, R. A.; Martín-Fontecha, M.; Gonzalez, A.; de la Fuente, T.; Vázquez-Villa, H.; García-Cárceles, J.; Botta, J.; McCormick, P. J.; Benhamú, B.; Pardo, L.; López-Rodríguez, M. L. A new serotonin 5-HT6 receptor antagonist with procognitive activity-importance of a halogen bond interaction to stabilize the binding. Sci. Rep. 2017, 7, 41293. (45) Jiang, S.; Zhang, L.; Cui, D.; Yao, Z.; Gao, B.; Lin, J.; Wei, D. The important role of halogen bond in substrate selectivity of enzymatic catalysis. Sci. Rep. 2016, 6, 34750. (46) 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. (47) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; van Gunsteren, W. F. Definition and testing of the GROMOS force-field versions 54A7 and 54B7. Eur. Biophys. J. 2011, 40, 843−856. (48) Reif, M. M.; Hünenberger, P. H.; Oostenbrink, C. New interaction parameters for charged amino acid side chains in the GROMOS force field. J. Chem. Theory Comput. 2012, 8, 3705−3723. (49) Reif, M. M.; Winger, M.; Oostenbrink, C. Testing of the GROMOS force-field parameter set 54A8: structural properties of electrolyte solutions, lipid bilayers, and proteins. J. Chem. Theory Comput. 2013, 9, 1247−1264. (50) Malde, A. K.; Zuo, L.; Breeze, M.; Stroet, M.; Poger, D.; Nair, P. C.; Oostenbrink, C.; Mark, A. E. An automated force field topology builder (ATB) and repository: version 1.0. J. Chem. Theory Comput. 2011, 7, 4026−4037. (51) Koziara, K. B.; Stroet, M.; Malde, A. K.; Mark, A. E. Testing and validation of the Automated Topology Builder (ATB) version 2.0: prediction of hydration free enthalpies. J. Comput.-Aided Mol. Des. 2014, 28, 221−233. (52) Baase, W. A.; Liu, L.; Tronrud, D. E.; Matthews, B. W. Lessons from the lysozyme of phage T4. Protein Sci. 2010, 19, 631−641. (53) Vallurupalli, P.; Chakrabarti, N.; Pomès, R.; Kay, L. E. Atomistic picture of conformational exchange in a T4 lysozyme cavity mutant: an experiment-guided molecular dynamics study. Chem. Sci. 2016, 7, 3602−3613. (54) Xie, B.; Nguyen, T. H.; Minh, D. D. L. Absolute binding free energies between T4 lysozyme and 141 small molecules: calculations based on multiple rigid receptor configurations. J. Chem. Theory Comput. 2017, 13, 2930−2944. (55) Eriksson, A. E.; Baase, W. A.; Zhang, X.-J.; Heinz, D. W.; Blaber, M.; Baldwin, E. P.; Matthews, B. W. Response of a protein structure to cavity-creating mutations and its relation to the hydrophobic effect. Science 1992, 255, 178−183. (56) Liu, L.; Baase, W. A.; Matthews, B. W. Halogenated benzenes bound within a non-polar cavity in T4 lysozyme provide examples of I··· S and I···Se halogen-bonding. J. Mol. Biol. 2009, 385, 595−605. (57) 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., Jr. General atomic I


Article

Journal of Chemical Theory and Computation and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347−1363. (58) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision A.2; Gaussian, Inc.: Wallingford, CT, 2009. (59) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (60) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (61) Hariharan, P. C.; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. Theor. Chim. Acta 1973, 28, 213−222. (62) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. J. Chem. Phys. 1982, 77, 3654−3665. (63) Glukhovtsev, M. N.; Pross, A.; McGrath, M. P.; Radom, L. Extension of Gaussian-2 (G2) theory to bromine-and iodinecontaining molecules: use of effective core potentials. J. Chem. Phys. 1995, 103, 1878−1885. (64) Barone, V.; Cossi, M. Quantum calculation of molecular energies and energy gradients in solution by a conductor solvent model. J. Phys. Chem. A 1998, 102, 1995−2001. (65) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model. J. Comput. Chem. 2003, 24, 669−681. (66) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilizaion of AB initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117−129. (67) Bondi, A. van der Waals volumes and radii. J. Phys. Chem. 1964, 68, 441−451. (68) The MK radius for iodine in GAMESS-US is unknown, and therefore, a value of 1.8 Å is attributed by default. This radius is clearly too small, and therefore, a value of 2.3 Å (used in Gaussian 09 for bromine) was assigned. (69) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (70) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: a message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43−56. (71) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718. (72) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1−2, 19−25. (73) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. Nucleic Acids Res. 2000, 28, 235−242.

(74) Hermans, J.; Berendsen, H. J. C.; van Gunsteren, W. F.; Postma, J. P. M. A consistent empirical potential for water-protein interactions. Biopolymers 1984, 23, 1513−1518. (75) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (76) Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (77) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: a new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (78) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: an N· log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (79) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577−8593. (80) Pall, S.; Hess, B. A flexible algorithm for calculating pair interactions on SIMD architectures. Comput. Phys. Commun. 2013, 184, 2641−2650. (81) Hess, B. P-LINCS: a parallel linear constraint solver for molecular simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (82) Miyamoto, S.; Kollman, P. A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952−962. (83) Wilcken, R.; Zimmermann, M. O.; Lange, A.; Zahn, S.; Boeckler, F. M. Using halogen bonds to address the protein backbone: a systematic evaluation. J. Comput.-Aided Mol. Des. 2012, 26, 935−945. (84) 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. (85) 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. (86) 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. (87) 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. (88) Silverman, B. W. Density Estimation for Statistics and Data Analysis; Chapman and Hall: London, 1986. (89) Campos, S. R. R.; Baptista, A. M. Conformational analysis in a multidimensional energy landscape: study of an arginylglutamate repeat. J. Phys. Chem. B 2009, 113, 15989−16001. (90) Ibrahim, M. A. A. Molecular mechanical perspective on halogen bonding. J. Mol. Model. 2012, 18, 4625−4638. (91) de Groot, B. L.; Hayward, S.; van Aalten, D. M. F.; Amadei, A.; Berendsen, H. J. C. Domain motions in bacteriophage T4 lysozyme: a comparison between molecular dynamics and crystallographic data. Proteins: Struct., Funct., Genet. 1998, 31, 116−127. (92) Chen, J.-L.; Yang, Y.; Zhang, L.-L.; Liang, H.; Huber, T.; Su, X.C.; Otting, G. Analysis of the solution conformations of T4 lysozyme by paramagnetic NMR spectroscopy. Phys. Chem. Chem. Phys. 2016, 18, 5850−5859. (93) 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. (94) Rowe, R. K.; Ho, P. S. Relationships between hydrogen bonds and halogen bonds in biological systems. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2017, 73, 255−264.

J


Biomolecular Simulations of Halogen Bonds with a GROMOS Force

Recommend Documents