Improved Description of Sulfur Charge Anisotropy in OPLS Force

Subscriber access provided by University of Massachusetts Amherst Libraries

Article

Improved Description of Sulfur Charge Anisotropy in OPLS Force Fields. Model Development and Parameterization Xin Cindy Yan, Michael J Robertson, Julian Tirado-Rives, and William L Jorgensen J. Phys. Chem. B, Just Accepted Manuscript • Publication Date (Web): 19 Jun 2017 Downloaded from http://pubs.acs.org on June 20, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 39

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

The Journal of Physical Chemistry

Improved Description of Sulfur Charge Anisotropy in OPLS Force Fields. Model Development and Parameterization Xin Cindy Yan, Michael J. Robertson, Julian Tirado-Rives, and William L. Jorgensen* Department of Chemistry, Yale University, New Haven, Connecticut 06520-8107, United States Received May 00, 2017 ABSTRACT: The atomic point-charge model used in most molecular mechanics force fields does not represent well the electronic anisotropy that is featured in many directional noncovalent interactions. Sulfur participates in several types of such interactions with its lone pairs and σ-holes. The current study develops a new model, via the addition of off-atom charged sites, for a variety of sulfur compounds in the OPLS-AA and OPLS/CM5 force fields to address the lack of charge anisotropy. Parameter optimization is carried out to reproduce liquid-state properties, torsional and non-covalent energetics from reliable quantum mechanical calculations, and electrostatic potentials. Significant improvements are obtained for computed free energies of hydration, reducing the mean unsigned errors from ca. 1.4 to 0.4-0.7 kcal/mol. Enhanced accuracy in directionality and energetics is also obtained for molecular complexes with sulfurcontaining hydrogen and halogen bonds. Moreover, the new model reproduce the unusual conformational preferences of sulfur-containing compounds with 1,4-intramolecular chalcogen bonds. Transferability of the new force field parameters to cysteine and methionine is verified via molecular dynamic simulations of blocked dipeptides. The study demonstrates the effectiveness of using off-atom charge sites to address electronic anisotropy, and provides a parameterization methodology that can be applied to other systems.

ACS Paragon Plus Environment

1


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 2 of 39

INTRODUCTION Anisotropy of molecular electrostatic potentials is widely observed in chemical compounds and is responsible for directional non-covalent interactions. Quantum mechanics (QM) calculations and statistical analyses of protein and small-molecule crystal structures indicate strong orientation dependence of hydrogen bonding.1-4 Halogen bonding formed between an electropositive σ-hole on a halogen atom X and a Lewis base L also exhibits a clear C-X···L linear preference.5,6 Similar to halogen bonding, chalcogen bonding is another type of directional non-covalent interactions that involves σ-holes on atoms in the chalcogen group such as sulfur or selenium.7-9 The electronic anisotropy of sulfur, whose role in non-covalent interactions was largely underestimated in the past, has been highlighted in a number of recent studies.7,10-15 In many compounds, the electrostatic potential around sulfur has both positive and negative regions, which allows it to act as hydrogen and halogen bond acceptors in lone-pair directions, and chalcogen bond donors along the extension of covalent bonds with sulfur. In a series of spectroscopic studies from Biswal and coworkers,11,16,17,18 large red-shifts of the N-H and O-H stretching frequencies were observed for N-H···S and O-H···S hydrogen bonds in analogues of tryptophan-methionine and tyrosine-methionine complexes. The magnitudes of the red-shifts were comparable to the N/O-H···O counterparts, indicating similar hydrogen-bonding strength. In benchmark studies using CCSD(T) calculations for biologically relevant complexes involving sulfur,10,13 strong hydrogen bonds were identified in complexes where dimethylsulfide acts as the hydrogen-bond acceptor. In the formation of halogen bonds, sulfur can function as a Lewis base. Protein crystal structures and QM calculations have revealed that halogen bonds formed between halobenzenes and methionine can be exploited to fine-tune binding affinity, geometry, and


2

Page 3 of 39

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


selectivity in structure-based drug design.19-21 Combining computational, FTIR and Raman spectroscopy results, Hauchecorne et al. estimated that the strengths of halogen bonding interactions between CF3Br, CF3I and dimethylsulfide are -2.3 and -4.3 kcal/mol in liquid krypton solution.22 In contrast to halogen bonding, which is mostly associated with intermolecular interactions, interesting chalcogen bonding involving sulfur has been found for 1,4- and 1,5-intramolecular interactions.23,24 Such interactions have served as useful directional forces for conformational control in crystal engineering25 and drug design.15 Despite the important role that electronic anisotropy plays in non-covalent interactions involving sulfur, most biomolecular force fields represent sulfur electrostatics with a single isotropic point-charge, which fails to account for electropositive σ-holes, and underestimates electronegativity in lone-pair directions.4,14,26-28 For example, the anisotropy is not accounted for in the OPLS-AA force field,29 whose charges were obtained largely from fitting to reproduce properties of organic liquids; nor in the OPLS/CM5 force field,30 whose charges were based on Charge Model 5 that utilizes gas-phase atomic charges form a Hirshfeld population analysis.31 Due to this inadequacy, the representations of sulfur compounds in the OPLS force fields lead to substantial errors in free energies of hydration, despite the fact that heats of vaporization and molecular volumes are well reproduced.30 Past approaches to improve the treatment of electronic anisotropy include replacing point-charges with multipoles, which can be derived from distributed multipole analysis32,33 or molecular electrostatic potential fitting.34,35 Nevertheless, multipole-based force fields are accompanied with considerable computational cost and lack of transferability, which make their use uncommon for biomolecular simulations. Alternatively, charge anisotropy can be approximated by placing extra off-atom charge sites on the targeted heavy atom. Previous studies have added negative charge sites to mimic lone-pairs in anions and


3


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 4 of 39

hydrogen bond acceptors (N, O, S etc.),36-39 and positively charged sites on halogen atoms to represent σ-holes.40-42 With adequate optimization, the strength and directionality of noncovalent interactions in these models can be properly represented with little additional computational cost. As reported here, a new model for the representation of sulfur in the OPLS-AA and OPLS/CM5 force fields have been developed to yield an improved description of electronic anisotropy. Parameter optimization was carried out for 14 sulfur-containing organic molecules. In most cases, off-atom point charges have been added to reproduce molecular electrostatic potentials (ESP) and multipole moments determined from QM calculations. Force field parameters were optimized with an emphasis on reproducing thermodynamic properties of pure liquids, minimizing mean-unsigned errors (MUEs) for free energies of hydration, and maximizing transferability among similar compounds. The study also includes a detailed comparison of the energetics and geometries of molecular complexes with sulfur participating in hydrogen or halogen bonding interactions. Consideration of 1,4-intramolecular S···N and S···O chalcogen bonding included fitting torsional parameters to high-level QM torsional energy scans, and the results are compared to a statistical survey of small molecule crystal structures in the Cambridge Structure Database (CSD). Furthermore, the transferability of the new sulfur parameters from small organic compounds to methionine and cysteine was validated by fitting new χ1 and χ2 torsion parameters and comparing rotamer distributions to QM results and experimental NMR measurements. Consistent with the treatment of halogen bonding,42 the force fields with the added interaction sites are referred to as OPLS-AAx and OPLS/CM5x.

METHODS


4

Page 5 of 39

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


The parameterization process entailed a mostly automated workflow facilitated by the Nelder– Mead simplex algorithm.43,44 First, each extra site (XS) is modeled as a charge site attached to the sulfur atom by a stiff bond. The initial guesses for the geometry, including relative distances, angles, and dihedral angles of extra sites with respect to the sulfur atom, were made based on past studies.37,39,42,45 Then point-charges on the extra sites, sulfur, and the atoms covalently bonded to sulfur were fit by an in-house variation of the restrained electrostatic potential (RESP) methodology.46 The details of the RESP fitting are provided in the Supporting Information. The resulting charges were further adjusted by the program to reproduce the dipole and quadrupole moments scaled up by 10-20% from gas-phase values to account for solvent polarization effects. The geometry of extra sites and partial charges were optimized in an iterative fashion, until good agreement was obtained for both the ESP and multipole moments. Next, sulfur Lennard-Jones (LJ) parameters were optimized to reproduce experimental thermodynamic properties of pure liquids, a central concept in the development of OPLS force fields. The final model was further tested through Monte Carlo/Free Energy Perturbation (MC/FEP)47,48 calculations to determine free energies of hydration, as well as for transferability among related compounds. After all of the results were satisfactory, the needed torsional parameters were fit to QM torsional energy profiles. The parameterization process was repeated iteratively to ensure all parameters were self-consistent and the overall model was optimal. Metropolis Monte Carlo (MC) simulations of pure liquids were carried out using the BOSS program.49 The simulations were conducted in the isothermal-isobaric ensemble (NPT) at 1 atm pressure and at 5.96 ºC for methanethiol and 25 ºC for all other molecules. Simulations were performed using periodic cubic cells consisting of ca. 400 molecules. 8 million MC configurations were used for both the equilibration and averaging stages (8M/8M) in each gas-


5


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 6 of 39

phase run, and 30M/50M in each liquid-phase run. Volume changes were attempted every 2500 configurations. Intermolecular interactions were truncated at 12 Å, with quadratic feathering to zero over the last 0.5 Å. The ranges for attempted moves were adjusted to yield ca. 40% acceptance rate. The molecules were treated as completely flexible. To compute absolute free energies of hydration, MC/FEP calculations were carried out with the BOSS program following standard procedures.49,50 Simulations were conducted at 25 ºC and 1 atm in the NPT ensemble, using periodic cubic cells consisting of a single solute molecule surrounded by 512 TIP4P waters.45 Solute moves and volume changes were attempted every 100 and 3125 configurations, respectively. Ranges for translations and rotations of ±0.10 and ±6.0º for the solutes and ±0.22 and ±15.0º for water were chosen to produce MC acceptance rates of 30-50%. Full annihilations of the solute molecules were achieved in two successive steps, first neutralizing the atomic charges, then decreasing the LJ parameters to zero while simultaneously shrinking the molecules.31 The FEP calculation of each phase of the annihilation was carried out in 25 double-wide-sampling windows, with 8M/8M configurations in the gas phase and 30M/70M configurations in the aqueous phase for each window. The uncertainties of the calculated absolute free energies are below 0.15 kcal/mol, which are smaller than the typical experimental uncertainties of 0.3 kcal/mol.51-53 All hydration free energies reported in this study have been corrected for the neglected LJ interactions beyond the cutoff distance.54 In order to develop the required new torsional parameters, QM torsional energy profiles were calculated using the Gaussian 09 software revision C.0155 at the MP2/aug-ccpVTZ//MP2/6-311G(d,p) level in intervals of 10º. Calculations at similar levels have been found to yield accurate conformational energetics.56,57 The molecular mechanics (MM) torsional energy profiles were computed with the “dihedral drive” script included in the BOSS software. The


6

Page 7 of 39

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


torsional energy term in the force field consists of a Fourier series summed over all dihedral angles in the molecule (eq 1). The torsional parameters were determined using the Etorsion = ∑ {V1,i (1 + cos ϕi ) / 2 + V2,i (1 − cos 2ϕi ) / 2 + V3,i (1 + cos3ϕi ) / 2 + V4,i (1 − cos 4ϕi ) / 2}

(1)

i

multivariate least-squares-fitting method, aiming at minimizing residual errors using the smallest number of fitting coefficients. For the dihedral angles considered in this study, all of the V4 terms and some of the V1-V3 terms were too small to affect the quality of fitting, so they were excluded. Interaction energies and geometries were computed for various molecular complexes with methanethiol, dimethyl sulfide, or dimethyl disulfide. For each complex, full geometry optimizations were carried out with the BOSS program using the original and the new force field parameters. The previous implementation of the OPLS-AAx force field with an extra pointcharge representing the σ-hole on halogen was used for the calculations of halogen-containing compounds.42 For complexes involving ions, the effects of intermolecular polarization with inducible dipoles were considered for all four variations of the force fields following the protocols documented in the OPLS-AAP force field.58,59 The polarizabilities, α, for carbon and sulfur were assigned as 1.0 and 2.3 Å3, respectively. In addition, intermolecular interaction scans of sulfur-containing compounds with water were conducted by QM and MM calculations. As illustrated in Figure S1, the in-plane (IP) scans were performed by changing the orientation of the water oxygen atom with respect to sulfur in the plane consisting of sulfur and its two covalently bonded atoms in intervals of 10º. In out-of-plane (OOP) scans, water monomers were placed in the perpendicular plane bisecting the covalent bonds. The corresponding QM scans were carried out with MP2/aug-cc-pVQZ//MP2/cc-pVTZ calculations. For all molecular


7


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 8 of 39

complexes, basis set superposition errors (BSSE) were removed via the Boys-Bernardi counterpoise method.60 New torsion parameters for cysteine and methionine were also developed with the addition of the extra sites. The fitting was performed with a steepest descent algorithm to minimize the Boltzmann-weighted energy differences between QM scans at the B2PLYPD3BJ/aug-cc-pVTZ//ωB97X-D/6-311++G(d,p) level and the corresponding MM scans, as described previously.61 The sole difference is that in the case of cysteine, additional χ1 QM scans were performed with the backbone held in a polyproline II conformation and used in the parameterization. This was found to improve the agreement for the computed rotamer populations and experimental data. Aqueous-phase simulations for blocked cysteine and methionine dipeptides were performed, again using the previously described protocols.61

RESULTS New Models for Sulfur-Containing Compounds. The current work initially focused on 14 sulfur-containing molecules including thiols, sulfides, disulfides, a sulfoxide, and heterocycles (Figure S2). They serve as prototypes for sulfur-containing components of proteins and drug-like molecules. For each type of compound, a full set of parameters was developed for the OPLSAAx force field. Then the torsional and LJ parameters were directly transferred to the OPLS/1.20*CM5x force field,30,54 which is appropriate when the OPLS-AAx charges are unavailable. Three different arrangements of extra sites were tested (Figure 1). Model A places a single extra site on sulfur along the bisector of its covalent bonds. This model is reminiscent of the TIP4P water,45 which improved the representation of the quadrupole moment and liquid


8

Page 9 of 39

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


structure. In model B, two extra sites were appended in the manner of lone pairs. Alternatively, model C features a geometry with two positively charged extra sites along the extension of the covalent bonds with sulfur in a way analogous to the representation of halogen σ-holes in the implementation of the OPLS-AAx force field.42 Optimizations were carried out to refine the geometries of the extra sites, and to obtain the best partial charges on extra sites, sulfur,

Figure 1. Illustration of the three different arrangements of extra sites (in purple) in three models of methanethiol, dimethyl sulfide, and dimethyl disulfide.

and all atoms covalently bonded to sulfur atoms. Sulfur LJ parameters were also included in the optimization process. Other force-field parameters were left unchanged to preserve maximally the original OPLS-AA force field. For each model, the top two sets of parameters were compared for computing thermodynamic properties, electrostatic potentials, and interaction energies of complexes (Table S1).


9


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 10 of 39

The best results were obtained when the extra sites were placed in the lone-pair manner, though the reference XS-S-XS angle optimized to 160° (Table S2), and each site is positioned 0.70 Å from the sulfur atom. A stiff bond-stretching force constant of 600 kcal/(mol·Å2) is used to restrain the extra sites from deviating significantly from the equilibrium distance. The extra sites are placed symmetrically with respect to the plane spanned by sulfur covalent bonds. All angle bending force constants involving the extra site are 200 kcal/(mol·deg2). For dihedral angles involving the extra site, the Fourier coefficients V1-V4 in the torsional potential (eq 1) are set to zero. The old and new non-bonded parameters are listed in Table 1. The partial charge for each extra site in both OPLS-AAx and OPLS/CM5x models is -0.20 e, comparable to other studies that take a similar approach.37-39 In the new OPLS-AAx force field, partial charges on sulfur become positive, though the sum of charges on sulfur and the extra sites remains similar to the original OPLS-AA force field. In addition, charges on carbon and hydrogen atoms covalently bonded to sulfur are made less positive in most of the new model to accommodate for the increase in dipole moment due to the addition of extra sites. In the OPLS/CM5x model, the CM5 charges scaled by 1.20 are used for all atoms except sulfur, whose charge is obtained from the 1.20*CM5 charge + 0.40 e to balance the extra sites. The new model for dimethyl sulfoxide (DMSO) is an exception that does not contain extra sites. QM calculations demonstrate that the electron-withdrawing effect of oxygen significantly reduces the negative electrostatic potential on the sulfur lone pair. However, charges on carbon and sulfur are made more positive and oxygen more negative in the OPLSAAx model to increase the dipole moment.


10

Page 11 of 39

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


To accommodate the new charge distributions, sulfur LJ parameters were adjusted to reproduce the properties of pure liquids. The new LJ parameters result in reduced van der Waals (vdW) attractions for thiols and DMSO, and increased ones for sulfides. For disulfides, the modified parameters σ and ε impact the vdW interactions in opposite directions. Overall, the new LJ parameters exhibit good transferability among sulfur compounds, and can be used with both the OPLS-AAx and CM5x charge models.

Table 1. Non-Bonded Parameters for Sulfur Compounds in the Old and New Models OPLS-AA

OPLS-AAx

type

q

σ

ε

q

σ

ε

SH

-0.335

3.600

0.425

0.085

3.639

0.381

S in alkyl thiols

SH

-0.220

3.600

0.425

0.190

3.639

0.381

S in thiophenol

HS

0.155

0

0

0.147

0

0

H on S in thiols

CT

0.060

3.500

0.066

0.048

3.500

0.066

C(S): RCH2SH

CT

0

3.500

0.066

-0.012

3.500

0.066

C(S): CH3SH

HC

0.060

2.500

0.030

0.060

2.500

0.060

H on alkyl CT

CA

0.065

3.550

0.070

0.063

3.550

0.070

C(S), thiophenol

S

-0.335

3.600

0.355

0.132

3.343

0.392

S in sulfides

CT

0.0475

3.500

0.066

0.014

3.500

0.066

C(S): CH2, sulfides

CT

-0.0125

3.500

0.066

-0.046

3.500

0.066

C(S): CH3, sulfides

S

-0.220

3.600

0.355

0.243

3.343

0.392

S in thioanisole

CA

0.0525

3.550

0.070

0.023

3.550

0.070

C(SMe),thioanisole

S

-0.2175

3.600

0.355

0.236

3.523

0.313

S in disulfides

CT

0.0975

3.500

0.066

0.044

3.500

0.066

C(S): CH2, disulfides

CT

0.0375

3.500

0.066

-0.016

3.500

0.066

C(S): CH3, disulfides

atom name


11


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 39

SZ

0.130

3.560

0.395

0.203

3.682

0.237

S in sulfoxide

OY

-0.420

2.930

0.280

-0.510

2.861

0.161

O in sulfoxide

CT

0.025

3.500

0.066

0.0335

3.500

0.066

C(S): CH2, sulfoxide

CT

-0.035

3.500

0.066

0.0265

3.500

0.066

C(S): CH3, sulfoxide

SA

-0.030

3.600

0.355

0.454

3.343

0.392

S in thiazole

CR

0.242

3.550

0.070

0.177

3.550

0.070

C2 in thiazole

CW

-0.299

3.550

0.070

-0.318

3.550

0.070

C5 in thiazole

XS

˗

˗

˗

-0.200

0

0

extra site, all

As illustrated in Figure S3, new torsional parameters were developed for 12 types of sulfur-containing dihedral angles using QM torsion scans at the MP2/aug-cc-pVTZ level as fitting references. In the torsional energy profiles of thiol compounds, the original OPLS-AA torsional parameters from 1996, which were developed based on HF/6-31G* calculations,29 tend to overestimate the heights of rotational barriers. Better agreement is obtained between the new MM and QM results. Improvements in torsional energetics are also obtained for sulfides. In particular, the new parameterization lowers the torsional barrier for thioanisole by ca. 1.5 kcal/mol, leading to close alignment with the QM profile. Similar refinements were made for disulfides. All the torsional parameters demonstrate good transferability to the OPLS/CM5x force field.

Thermodynamic properties of pure liquids. Thermodynamic properties determined from the pure liquid MC simulations of the 14 sulfur compounds are compared to the experimental data in Figure 2 and Tables S3 and S4. High accuracy for the liquid densities is preserved in the new model with MUE of 0.011 (OPLS-AAx) and 0.018 (OPLS/CM5x) g/cm3, which are comparable to 0.012 and 0.015 g/cm3 obtained from the original force fields. The most significant


12

Page 13 of 39

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


improvement is observed for disulfides, where the original OPLS-AA force field overestimates the densities by ca. 0.03 g/cm3. Notably, improved accuracy is found for heats of vaporization with MUEs decreasing from 0.61 (OPLS-AA) to 0.36 kcal/mol (OPLS-AAx) and from 0.68 (OPLS/CM5) to 0.49 kcal/mol (OPLS/CM5x). The original OPLS-AA force field tends to yield larger heats of vaporization for the tested compounds except thiols; in particular, 1-2 kcal/mol discrepancies are found for the disulfides and DMSO. The new OPLS-AAx force field is able to reduce the error to less than 0.5 kcal/mol for a majority of the compounds. For DMSO and thiazole, the error remains around 1 kcal/mol due to trade-offs for other properties. With CM5 charges, the largest errors with original model are overestimates of 1.2 kcal/mol for diethyl disulfide and 4.2 kcal/mol for DMSO, while the new model reduces these errors by half. Heat capacities were also computed for six liquids with available experimental data (Table S5). There are no large discrepancies. The new OPLS-AAx model gives a MUE of 2.2 cal/(mol·K), which is 1 cal/(mol·K) lower than with OPLS-AA, and the MUEs are 1.9 and 2.9 cal/(mol·K) using OPLS/CM5 and OPLS/CM5x.


13


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 14 of 39

Figure 2. Densities (left, g/cm3) and heats of vaporization (right, kcal/mol) for 14 sulfurcontaining compounds from experiment and MC simulations. The blue and orange lines are obtained from least square fitting of the data with the OPLS-AAx and OPLS/CM5x force fields, and the red dashed line denotes an ideal correlation.

Free energies of hydration. Faithful reproduction of free energies of hydration is of great importance for simulations of biomolecular systems including protein-ligand complexes. With MUEs of 1.32 and 1.54 kcal/mol with the original OPLS-AA and OPLS/CM5 force fields, the accuracy of hydration free energies for sulfur compounds has been a nagging problem.31 The computed values are systematically too positive, indicating insufficiently favorable interactions


14

Page 15 of 39

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


Table 2. Free Energies of Hydration ∆Ghyd (kcal/mol) of Sulfur-Containing Compounds molecule

OPLS-AA

OPLS-AAx

OPLS/CM5

OPLS/CM5x

exptl.a

methanethiol

-0.40

-1.04

0.19b

-0.82

-1.22

ethanethiol

-0.42

-0.82

-0.19

-1.01

-1.14

propanethiol

-0.11

-0.56

-0.06

-0.64

-1.10

thiophenol

-1.46

-3.11

-1.24b

-2.46

-2.55

dimethyl sulfide

-0.12

-1.47

0.75b

-1.10

-1.61

diethyl sulfide

0.05

-1.44

0.30

-2.30

-1.43

ethyl methyl sulfide

0.16

-1.64

0.49

-1.85

-1.50c

thiane

-0.44

-2.12

-0.09

-3.06

-1.73d

thioanisole

-0.92

-1.79

-1.16b

-3.37

-2.73

dimethyl disulfide

-1.91

-1.09

-0.25

-0.96

-1.54c

diethyl disulfide

-1.13

-1.57

-0.39

-1.87

-1.43c

DMSO

-6.91

-10.91

-12.42b

-12.90

-10.11

thiophene

N/A

N/A

-0.24

-1.67

-1.40

thiazole

-2.34

-4.05

-2.55

-3.41

(-4.25)d

MUE

1.32

0.37

1.54

0.69

a

Ref 53. bRef 54. cRef 62. dComputed in implicit water using the Polarizable Continuum Model (PCM) at the MP2/aug-cc-pVQZ level.

with water, though this discrepancy is at odds with the overly large heats of vaporization. Further analyses reveal that the seemingly contradictory results stem from the lack of charge anisotropy. Indeed, the results in Table 2 confirm that the new OPLS-AAx and OPLS/CM5x force fields substantially reduce the MUEs for free energies of hydration to 0.37 and 0.69 kcal/mol. The most notable improvements are found for sulfides for which the previous results were too positive by


15


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 16 of 39

ca. 2 kcal/mol. The discrepancy is effectively resolved by the implementation of the extra sites, which made sulfides better hydrogen bond acceptors. Another key repair is for DMSO, where the OPLS-AAx model reduces the error from 3.74 to 0.80 kcal/mol. Nevertheless, due to an overestimation of the dipole moment for DMSO, similar improvement is not seen for the CM5 charge models. QM Self-Consistent Reaction Field (SCRF) calculations were carried out using the Polarizable Continuum Model (PCM) to provide an alternative estimate of the free energy of hydration for thiazole, whose experimental value is unavailable.63 The PCM calculations predict a ∆Ghyd of -4.25 kcal/mol, which is in reasonable agreement with the FEP results from the OPLS-AAx and OPLS/CM5x models, but it is ca. 2 kcal/mol more negative than from the original force fields. For thiophene, the free energy of hydration determined using the new OPLS/CM5x model is in close agreement with the experimental value. In this case, OPLS-AA parameters have not been reported, so the corresponding entry is missing in Table 2. Thus, it is recommended to use OPLS/CM5x for modeling thiophenes and other sulfur-containing heterocycles.

Electrostatics. Dipole moments of both the new and the old models are compared to the QM gas-phase values in Table S6. Although a 15-20% enhancement of dipole moments is normal to implicitly account for solvent polarization effects, disulfides are outliers in the OPLS-AA force with ca. 60% increases. The larger dipole moments are in line with the overestimation of heats of vaporization for these compounds, while reduction of the dipole moments results in larger deviations for the free energies of hydration. Incorporation of extra sites in the new model solves the dilemma and provides improvements in both thermodynamic and electrostatic properties. Dipole moments obtained from unscaled CM5 charges are in excellent accord with experimental


16

Page 17 of 39

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


values, since they were developed to reproduce this quantity.64 Thus, the OPLS/CM5 charges with the 1.2 scaling yield dipole moments 20% greater than the observed gas-phase values (Table S6). However, with the addition of the extra charged sites, the dipole moments in the CM5x model are further increased by ca. 0.2 D. The relative root mean square (RRMS) errors of molecular electrostatic potential fitting are also tabulated in Table S6. Almost all OPLS-AAx and OPLS/CM5x results demonstrate better fitting to the QM molecular electrostatic potentials in comparison to the old model. The biggest improvements are found for thiophenol and disulfides in the OPLS-AAx force field and thiol and thiophene in the OPLS/CM5x force field, for which the RRMS errors are reduced by 20-35%. Overall, the results confirm that the addition of off-atom charge sites provides a more accurate description of the anisotropic nature of sulfur electrostatics.

Binding profiles of sulfur compound – water complexes. Interaction energies of water with methanethiol, dimethyl sulfide and dimethyl disulfide were determined in the geometries of Figure S1. Figures 3A-D show the in-plane (Figure S1) interaction energy profiles of methanethiol (MeSH), dimethyl sulfide (Me2S), dimethyl disulfide (DMDS), and thiophene with a water molecule. Although the binding profiles determined with all four force fields correctly


17


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 18 of 39

Figure 3. Binding energy (kcal/mol) profiles of molecular complexes determined by the MP2/aug-cc-pVQZ calculations, and MM calculations with the OPLS-AA, OPLS-AAx, OPLS/CM5, and OPLS/CM5x force field. In H and I, explicit polarization with induced dipoles was included in the MM calculations.

capture the shape of the MP2 scans, better quantitative agreement is obtained using OPLS-AAx and OPLS/CM5x. In addition, the interaction energies determined by the OPLS-AAx force field are uniformly more attractive than with OPLS/CM5x. From Figure 3A, the strongest hydrogen bonding interaction is found for the methanethiol – water complex at ca. 270°. The hydrogen


18

Page 19 of 39

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


bonding strengths predicted by all MM calculations are within 0.5 kcal/mol of the QM results. The computed binding-energy profiles using the OPLS-AA and OPLS/CM5 force fields display larger errors for the dimethyl disulfide – water complex (Figure 3C). The attractive interactions between the water molecule and sulfur are underestimated in the 200-230° range, and the interactions between water and the sulfur σ-hole region in the range of 80-180° are too strong with OPLS-AA. The out-of-plane scans (Figures 3E-G) with the OPLS-AA and OPLS/CM5 force fields demonstrate significant divergence from the QM results. The QM profiles indicate substantially stronger binding as the water molecules moves out-of-plane, reaching maximal binding near the perpendicular geometries with θ close to 90º. The old results are particularly poor for dimethyl sulfide, for which the interactions weaken rather than strengthen as θ increases. The results in Figure 3 with OPLS-AAx and OPLS/CM5x are significantly better, which can be attributed to the improved description of the electronic anisotropy around sulfur.

Binding profiles of sulfur compound – sodium ion complexes. Out-of-plane interaction energy profiles were determined for complexes of dimethyl sulfide and dimethyl disulfide with Na+. Explicit polarization with induced dipoles was included in the MM calculations to better account for the nature of the electrostatic interactions. In their absence, the interaction energies are too weak by ca. 10 kcal/mol. As shown in Figure 3H-I, the binding profiles obtained from the OPLS-AA and OPLS/CM5 force fields differ considerably from the QM results. In the QM results, there exists one local minimum around 60° in the dimethyl sulfide – Na+ interaction profile, and there are two minima at -40 and 60° for dimethyl disulfide. However, the profiles determined by the old force field model are monotonic and favor the co-planar geometry (θ = 0).


19


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 20 of 39

With the extra sites installed, the local minimum for the dimethyl sulfide – Na+ interaction surface is properly reproduced, though the strength of interactions is still systematically underestimated by ca. 3.5 kcal/mol. Likewise, the scans for dimethyl disulfide – Na+ with the new force fields yield correctly two energy-minima, though the minima near 90° are now too deep by 3 kcal/mol. Simply put, reasonable representation of the interaction of sulfur with cations requires explicit treatment of polarization effects for strength and the inclusion of extra interaction sites for geometry. Even then, challenges remain.

Hydrogen-bonded complexes. The interaction energies and geometries were determined for six molecular complexes, in which sulfur acts either as a hydrogen-bond donor or acceptor (Figure 4). The computed results obtained from the force fields are contrasted with those from QM calculations at the CCSD(T)/CBS level in Table 3.13

Figure 4. Hydrogen-bonded complexes with methanethiol and dimethyl sulfide.


20

Page 21 of 39

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


When methanethiol serves as the hydrogen-bond donor, the interaction energies from the QM calculations are attractive by 2-3 kcal/mol. The force fields yield similar results with the interactions stronger by ca. 0.5 kcal/mol using OPLS-AAx and OPLS/CM5x. The force-field model also gives good geometrical results, though hydrogen-bond lengths are generally 0.3 Å shorter than the QM results. This pattern is considered to be necessary to provide accurate densities for liquids.45 For the complexes in which sulfur acts as the hydrogen-bond acceptor, it is interesting to notice that the interactions can be as attractive as 5-6 kcal/mol according to the CCSD(T) calculations. In the past, divalent sulfur was not generally considered to be a good hydrogenbond acceptor in view of its larger atomic radius than for oxygen and nitrogen. Consistently, the hydrogen-bond strengths in these cases from the OPLS-AA and OPLS/CM5 force fields are too weak by an average of 1.6 kcal/mol compared to the QM references. The results are improved with the OPLS-AAx and OPLS/CM5x models, reducing the difference to 0.5 kcal/mol.

Halogen-bonding complexes. Halogen bonding has emerged as an alternative to hydrogen bonding in the design of materials and drugs.40-42,64-67 The current study examined halogen bonding in complexes between halobenzenes (C6H5X, X = Cl, Br, I) and methanethiol and dimethyl sulfide, which can serve as cysteine and methionine analogues. For MM geometry optimizations, the OPLS-AA force field with an additional point-charge representing the σ-hole on the halogen was adopted for the halobenzenes,42 and the original and new force field models were used to describe the sulfur containing compounds. The results are tabulated in Table 4 together with the available QM data at the CCSD(T)/CBS level.20,64 OPLS-AAx and OPLS/CM5x here means that the extra sites are included on both halogens and sulfur.


21


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 22 of 39

Table 3. Binding Energies (kcal/mol) and Geometries (distance in Å, angle in deg.) of Hydrogen Bonding Complexes Determined by MM and QM Calculations. method

Ebind

Rxy

θABC

Ebind

Rxy

θABC

MeSH-H2O(xy=SO, ABC=CSO)

formamide-MeSH(xy=NS,ABC=NSC)

QMa

-2.38

3.63

89

-6.37

3.52

85

OPLS-AA

-2.58

3.25

95

-5.00

3.35

85

OPLS-AAx

-2.94

3.24

93

-5.97

3.22

91

OPLS/CM5

-2.32

3.27

97

-4.67

3.38

79

OPLS/CM5x

-2.82

3.24

96

-5.77

3.22

90

MeSH-NH3(xy=SN, ABC=CSN)

H2O-Me2S(xy=OS, ABC=OSC)

QMa

-2.98

3.68

94

-5.39

3.30

79

OPLS-AA

-2.49

3.40

94

-3.90

3.20

86

OPLS-AAx

-2.86

3.38

94

-6.15

2.91

93

OPLS/CM5

-2.26

3.42

97

-3.61

3.26

82

OPLS/CM5x

-2.77

3.38

96

-6.17

2.93

91

MeSH-CH2O(xy=SO, ABC=CSO)

MeOH-Me2S(xy=OS, ABC=OSC)

QMa

-3.13

3.56

83

-5.72

3.28

80

OPLS-AA

-2.64

3.24

94

-4.08

3.21

85

OPLS-AAx

-2.91

3.22

91

-5.96

2.94

92

OPLS/CM5

-2.41

3.25

98

-3.86

3.25

81

OPLS/CM5x

-2.79

3.22

95

-6.03

2.96

90

a

CCSD(T)/CBS results from ref. 13.

The QM calculations for all complexes show that the optimized structures have S···X distances shorter than the sum of the van der Waals radii, indicating the formation of halogen bonds. The binding energies and halogen bonding distances determined by the OPLS-AAx and


22

Page 23 of 39

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


OPLS/CM5x force fields are in close agreement with the QM results. In addition, the doubly extra site model is able to correctly account for the relative strength of halogen bonding (Cl < Br < I) indicated by the QM calculations. With respect to the geometrical variables, θ(C-S···X) and θ(S···C-X), the results from all force fields are in good accord with the QM findings. However, the singly extra site force fields yield significantly too weak halogen bonding by 0.7-1.5 kcal/mol, and the S···X distances are 0.2-0.4 Å longer than in the QM geometries. In all, force fields with extra sites on both halogen and sulfur are preferred to represent sulfur-containing halogen bonding interactions with good accuracy.

Table 4. Binding Energies (kcal/mol) and Geometries (distance in Å, angle in deg.) of Halogen Bonding Complexes Determined by MM and QM Calculations Eint

R(SX)a

θ(SXC) θ(CSX)

QMb

-2.32

3.54

167

81

OPLS-AA

-1.53

3.66

179

79

OPLS-AAx

-2.11

3.61

173

80

OPLS/CM5

-1.73

3.80

172

70

OPLS/CM5x

-2.11

3.64

173

78

QMb

-3.08

3.56

171

85

OPLS-AA

-2.11

3.86

170

74

OPLS-AAx

-2.70

3.67

173

83

OPLS/CM5

-2.06

3.89

174

72

OPLS/CM5x

-2.66

3.70

173

82

method PhBr-MeSH

PhI-MeSH


23


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 24 of 39

PhCl-Me2S QMb

-2.41

3.48

163

81

OPLS-AA

-1.81

3.69

167

79

OPLS-AAx

-2.27

3.36

179

86

OPLS/CM5

-1.76

3.73

172

78

OPLS/CM5x

-2.28

3.39

177

85

QMb

-3.11

3.43

170

86

OPLS-AA

-2.24

3.69

169

80

OPLS-AAx

-2.96

3.38

164

88

OPLS/CM5

-2.18

3.72

174

79

OPLS/CM5x

-2.96

3.39

175

87

QMb

-4.13

3.39

168

89

OPLS-AA

-2.67

3.78

171

83

OPLS-AAx

-3.93

3.40

176

92

OPLS/CM5

-2.56

3.82

176

81

OPLS/CM5x

-3.91

3.42

176

90

PhBr-Me2S

PhI-Me2S

a

X denotes chlorine, bromine, or iodine. bCCSD(T)/CBS results from refs. 20 and 64.

1,4-Intramolecular Chalcogen Bonds with Sulfur. Sulfur is associated with some remarkable conformational effects. 1,4-O···S and 1,4-N···S intramolecular chalcogen bonds are fascinating conformational regulators that are seeing use in medicinal chemistry.15 In this study, statistical analyses of crystals of sulfur-containing heterocyclic compounds in the CSD were conducted using ConQuest.68 1869 and 242 crystal structures were identified to contain 1,4 O···S and 1,4


24

Page 25 of 39

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


N···S chalcogen bonds with interaction distance shorter than the sum of the van der Walls radii, and with the 1,4-dihedral angle within the range of ±20º. The distributions of distances and dihedral angles are illustrated in Figure S4, where the most probable values are found around 3.00 Å/0.0º for 1,4 O···S, and 2.90 Å/0.0º for 1,4 N···S. Two specific systems, 2-acetylthiophene and 2-(2’-thienyl)pyridine, were chosen to investigate 1,4-intramolecular chalcogen bonding using the CM5 charge model due to their ubiquity in the crystal database and drug-like molecules. Torsional energy profiles for the S-C-CO or S-C-C-N dihedrals were obtained with MP2/aug-cc-pVTZ//MP2/6-311G(d,p) calculations and MM scans (Figure 5). The QM results for 2-acetylthiophene clearly show a preference for the syn conformer by 0.90 kcal/mol ( or 0.82 kcal/mol at the MP2(full)/aug-ccpVTZ//MP2(full)/aug-cc-pVTZ level) due to chalcogen bonding, which corroborates the preference in the crystal database. However, dihedral scans with the original OPLS/CM5 force field do not differentiate the stability of the syn and anti conformers. Refitting the S-C-C-O torsional parameters using the new thiophene model in the OPLS/CM5x force field is able to reproduce the QM results, and the accompanying 3.03 Å S···O distance agrees well with the QM result of 2.98 Å. Similarly, the S-C-C-N torsional profile was examined for 2-(2’-thienyl)pyridine. Consistent with the crystal analysis, the QM results indicate that the syn conformer is 0.68 kcal/mol (or 0.65 kcal/mol at the MP2(full)/aug-cc-pVTZ//MP2(full)/aug-cc-pVTZ level) more stable than the anti conformer. In contrast, the profile obtained with the OPLS/CM5 force field gives the wrong preference favoring the anti conformer by 1.48 kcal/mol. After refitting, the new OPLS/CM5x torsional profile matches the QM torsional energetics. The S···N distance obtained with OPLS/CM5x is 2.94 Å, in close agreement with the QM value of 2.91 Å. Furthermore, the


25


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 26 of 39

S-C-C-N torsional profile for 2-(thiophen-2-yl)pyrimidine was examined as an additional check. The OPLS/CM5x implementation reduces the error in the barrier height by about one-half compared to OPLS/CM5, and there is good agreement for the optimal S···N distance between OPLS/CM5x (2.97 Å) and the QM geometry (2.96 Å). It is clear that in the absence of specific consideration of such systems and chalcogen bonding with sulfur, force fields can not correctly represent such novel conformational effects.

Figure 5. Torsional energy profiles for S-C-C-O and S-C-C-N dihedral angles in three compounds determined by MP2/aug-cc-pVTZ calculations, and MM calculations with OPLS/CM5 and OPLS/CM5x force fields.

Fitting Side-Chain Torsional Parameters for Cysteine and Methionine. In order for the updated sulfur model to be adopted in simulations of biomolecular systems, new side-chain torsional parameters were developed for cysteine and methionine for the OPLS-AA/M protein force field.61 Specifically, torsional parameters were determined for χ1 (N-Cα-Cβ-Sγ) and χ1’ (CCα-Cβ-Sγ) in cysteine and χ1 (N-Cα-Cβ-Cγ), χ1’ (C-Cα-Cβ-Cγ), and χ2 (Cα-Cβ-Cγ-Sδ) in methionine.


26

Page 27 of 39

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


Molecular dynamics (MD) simulations of blocked dipeptides for cysteine and methionine (Acetyl-X-NMe) were carried out in aqueous-phase simulations to validate the parameterization. Simulated and experimental 3J couplings values together with the χ1 rotamer populations are reported in Table S7. The simulated 3J coupling values obtained with and without the extra sites are similar and compare favorably with the NMR measurements. For the χ1 angle in cysteine, the MUE in rotamer populations is 10.8% with the new sulfur model, only marginally higher than with the original force field. Similar performance is also observed for the χ1 angle in methionine, where the MUE in rotamer distributions are 9.4 and 13.7% before and after the addition of the extra sites. Given that the error bars for individual χ1 rotamer populations are up to 2.6% in some cases, these results can be considered to be roughly equivalent. The close agreement between the MD results and experimental measurements confirm the transferability of the new sulfur model from small organic molecules to protein residues.

CONCLUSION The current study has addressed sulfur electronic anisotropy via the addition of off-atom charged sites to in the OPLS-AA and OPLS-AA/CM5 force fields. Testing included computation of properties of pure liquids and aqueous solutions. In particular, significant improvement is obtained for free energies of hydration of 14 diverse compounds with average errors reduced from 1.32 and 1.54 kcal/mol with the original OPLS-AA and OPLS/CM5 models to 0.37 and 0.69 kcal/mol with OPLS-AAx and OPLS/CM5x. The addition of extra sites on sulfur clearly results in much improved representation of charge anisotropy, as also reflected in better fitting to molecular electrostatic potentials and close agreement for the geometries and energetics of various hydrogen- and halogen-bonded complexes in comparison to benchmark QM


27


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 28 of 39

calculations. The study also examined 1,4-S···O and 1,4-S···N chalcogen bonds in three sulfurcontaining compounds. Although the original force fields do not reproduce the novel conformational preferences, new torsional parameters were derived in the context of OPLS/CM5x that resolve the problems. Lastly, new side-chain torsional parameters were derived for cysteine and methionine to be used with incorporation of the new sulfur model. Testing on the dynamics for dipeptides in aqueous solution revealed similarly satisfactory results with OPLS/AA-M and the augmented force field. In summary, the improvements seen in the current study indicate that off-atom charged sites, which lead to better representation of electronic anisotropy and directional non-covalent interactions, are important for the development of high quality force fields. Further validation and applications of the new sulfur model in simulations of biological systems will be reported.

ASSOCIATED CONTENT Supporting Information Mathematical derivation of a modified restrained electrostatic potential (RESP) fitting method is given. Computed densities, heats of vaporization, heat capacities, and electrostatics of sulfurcontaining molecules are tabulated. Torsional energy profiles of dihedral angles involving sulfur are plotted. The small molecule crystal survey results of the Cambridge Structural Database featuring heterocycles containing 1,4-intramolecular chalcogen bonds are shown. Results of MD simulations for the cysteine and methionine dipeptides are reported. This material is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION


28

Page 29 of 39

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


Corresponding Author *E-mail: [email protected] Notes The authors declare no competing financial interests. Acknowledgments Gratitude is expressed to the National Institute of Health (GM32136) for support of this work, and the Yale Center of Research Computing for computational resource. The authors thank Drs. Daniel Cole, Jonah Vilseck, and Israel Cabeza de Vaca for helpful discussions.

REFERENCES (1)

Taylor, R.; Kennard, O.; Versichel, W. Geometry of the N-H···O=C Hydrogen Bond. 1.

Lone-Pair Directionality. J. Am. Chem. Soc. 1983, 105, 5761-5766. (2)

Taylor, R.; Kennard, O. Hydrogen-Bond Geometry in Organic Crystals. Acc. Chem. Res.

1984, 17, 320-326. (3)

Platts, J. A.; Howard, S. T.; Bracke, B. R. F. Directionality of Hydrogen Bonds to Sulfur

and Oxygen. J. Am. Chem. Soc. 1996, 118, 2726-2733. (4)

Morozov, A. V.; Kortemme, T.; Tsemekhman, K.; Baker, D. Close Agreement between

the Orientation Dependence of Hydrogen Bonds Observed in Protein Structures and Quantum Mechanical Calculations. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 6946-6951.


29


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

(5)

Page 30 of 39

Ramasubbu, N.; Parthasarathy, R.; Murray-Rust, P. Angular Preferences of

Intermolecular Forces around Halogen Centers: Preferred Directions of Approach of Electrophiles and Nucleophiles around the Carbon ˗ Halogen Bond. J. Am. Chem. Soc. 1986, 108, 4308-4314. (6)

Metrangolo, P.; Meyer, F.; Pilati, T.; Resnati, G.; Terraneo, G. Halogen Bonding in

Supramolecular Chemistry. Angew. Chem., Int. Ed. 2008, 47, 6114-6127. (7)

Murray, J. S.; Lane, P.; Clark, T.; Politzer, P. σ-Hole Bonding: Molecules Containing

Group VI Atoms. J. Mol. Model. 2007, 13, 1033-1038. (8)

Wang, W.; Ji, B.; Zhang, Y. Chalcogen Bond: A Sister Noncovalent Bond to Halogen

Bond. J. Phys. Chem. A 2009, 113, 8132-8135. (9)

Lu, J.; Lu, Y.; Yang, S.; Zhu, W. Theoretical and Crystallographic Data Investigations of

Noncovalent S···O Interactions. Struct. Chem. 2011, 22, 757-763. (10) Wennmohs, F.; Staemmler, V.; Schindler, M. Theoretical Investigation of Weak Hydrogen Bonds to Sulfur. J. Chem. Phys. 2003, 119, 3208. (11) Zhou, P.; Tian, F.; Lv, F.; Shang, Z. Geometric Characteristics of Hydrogen Bonds Involving Sulfur Atoms in Proteins. Proteins: Struct., Funct., Bioinf. 2009, 76, 151-163. (12) Biswal, H. S.; Wategaonkar, S. Nature of the N-H···S Hydrogen Bond. J. Phys. Chem. A 2009, 113, 12763-12773. (13) Mintz, B. J.; Parks, J. M. Benchmark Interaction Energies for Biologically Relevant Noncovalent Complexes Containing Divalent Sulfur. J. Phys. Chem. A 2012, 116, 1086-1092.


30

Page 31 of 39

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


(14) Kramer, C.; Spinn, A.; Liedl, K. R. Charge Anisotropy: Where Atomic Multipoles Matter Most. J. Chem. Theory Comput. 2014, 10, 4488-4496. (15) Beno, B. R.; Yeung, K.-S.; Bartberger, M. D.; Pennington, L. D.; Meanwell, N. A. A Survey of the Role of Noncovalent Sulfur Interactions in Drug Design. J. Med. Chem. 2015, 58, 4383-4438. (16) Biswal, H. S.; Chakraborty, S.; Wategaonkar, S. Experimental Evidence of O-H···S Hydrogen Bonding in Supersonic Jet. J. Chem. Phys. 2008, 129, 184311. (17) Biswal, H. S.; Gloaguen, E.; Loquais, Y.; Tardivel, B.; Mons, M. Strength of NH···S Hydrogen Bonds in Methionine Residues Revealed by Gas-Phase IR/UV Spectroscopy. J. Phys. Chem. Lett. 2012, 3, 755-759. (18) Biswal, H. S.; Bhattacharyya, S.; Bhattacherjee, A.; Wategaonkar, S. Nature and Strength of Sulfur-Centred Hydrogen Bonds: Laser Spectroscopic Investigations in the Gas Phase and Quantum-Chemical Calculations. Int. Rev. Phys. Chem. 2015, 34, 99-160. (19) Alam, M.; Beevers, R. E.; Ceska, T.; Davenport, R. J.; Dickson, K. M.; Fortunato, M.; Gowers, L.; Haughan, A. F.; James, L. A.; Jones, M. W.; et al. Synthesis and SAR of Aminopyrimidines as Novel C-Jun N-Terminal Kinase (JNK) Inhibitors. Bioorg. Med. Chem. Lett. 2007, 17, 3463-3467. (20) 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.


31


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 32 of 39

(21) Lange, A.; Günther, M.; Büttner, F. M.; Zimmermann, M. O.; Heidrich, J.; Hennig, S.; Zahn, S.; Schall, C.; Sievers-Engler, A.; Ansideri, F.; et al. Targeting the Gatekeeper MET146 of C-Jun N-Terminal Kinase 3 Induces a Bivalent Halogen/Chalcogen Bond. J. Am. Chem. Soc. 2015, 137, 14640-14652. (22) Hauchecorne, D.; Moiana, A.; van der Veken, B. J.; Herrebout, W. A. Halogen Bonding to a Divalent Sulfur Atom: an Experimental Study of the Interactions of CF3X (X = Cl, Br, I) with Dimethyl Sulfide. Phys. Chem. Chem. Phys. 2011, 13, 10204-10213. (23) Nagao, Y.; Hirata, T.; Goto, S.; Sano, S.; Kakehi, A.; Iizuka, K.; Shiro, M. Intramolecular Nonbonded S···O Interaction Recognized in (Acylimino)thiadiazoline Derivatives as Angiotensin II Receptor Antagonists and Related Compounds. J. Am. Chem. Soc. 1998, 120, 3104-3110. (24) Iwaoka, M.; Takemoto, S.; Tomoda, S. Statistical and Theoretical Investigations on the Directionality of Nonbonded S···O Interactions. Implications for Molecular Design and Protein Engineering. J. Am. Chem. Soc. 2002, 124, 10613-10620. (25) Kusamoto, T.; Yamamoto, H. M.; Kato, R. Utilization of σ-Holes on Sulfur and Halogen Atoms for Supramolecular Cation···Anion Interactions in Bilayer Ni(dmit)2 Anion Radical Salts. Cryst. Growth Des. 2013, 13, 4533-4541. (26) Freddolino, P. L.; Harrison, C. B.; Liu, Y.; Schulten, K. Challenges in Protein-Folding Simulations. Nat. Phys. 2010, 6, 751-758.


32

Page 33 of 39

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


(27) Lu, Z.; Zhou, N.; Wu, Q.; Zhang, Y. Directional Dependence of Hydrogen Bonds: A Density-Based Energy Decomposition Analysis and Its Implications on Force Field Development. J. Chem. Theory Comput. 2011, 7, 4038-4049. (28) Lupyan, D.; Abramov, Y. A.; Sherman, W. Close Intramolecular Sulfur–Oxygen Contacts: Modified Force Field Parameters for Improved Conformation Generation. J. Comput.Aided Mol. Des. 2012, 26, 1195-1205. (29) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. (30) Vilseck, J. Z.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation of CM5 Charges for Condensed-Phase Modeling. J. Chem. Theory Comput. 2014, 10, 2802-2812. (31) Marenich, A. V.; Jerome, S. V.; Cramer, C. J.; Truhlar, D. G. Charge Model 5: An Extension of Hirshfeld Analysis for the Accurate Description of Molecular Interactions in Gaseous and Condensed Phases. J. Chem. Theory Comput. 2012, 8, 527-541. (32) Stone, A. J.; Alderton, M. Distributed Multipole Analysis Methods and Applications. Mol. Phys. 2002, 100, 221-233. (33) Ren, P.; Wu, C.; Ponder, J. W. Polarizable Atomic Multipole-Based Molecular Mechanics for Organic Molecules. J. Chem. Theory Comput. 2011, 7, 3143-3161. (34) Williams, D. E. Representation of the Molecular Electrostatic Potential by Atomic Multipole and Bond Dipole Models. J. Comput. Chem. 1988, 9, 745-763.


33


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 34 of 39

(35) Kramer, C.; Bereau, T.; Spinn, A.; Liedl, K. R.; Gedeck, P.; Meuwly, M. Deriving Static Atomic Multipoles from the Electrostatic Potential. J. Chem. Inf. Model. 2013, 53, 3410-3417. (36) Jorgensen, W. L.; Briggs, J. M.; Gao, J. A Priori Calculation of pKa’s for Organic Compounds in Water. The pKa of Ethane. J. Am. Chem. Soc. 1987, 109, 6857-6858. (37) Dixon, R. W.; Kollman, P. A. Advancing Beyond the Atom-Centered Model in Additive and Nonadditive Molecular Mechanics. J. Comput. Chem. 1997, 18, 1632-1646. (38) Wennmohs, F.; Schindler, M. Development of a Multipoint Model for Sulfur in Proteins: A New Parametrization Scheme to Reproduce High-Level Ab Initio Interaction Energies. J. Comput. Chem. 2005, 26, 283-293. (39) Zhu, X.; Mackerell, A. D. Polarizable Empirical Force Field for Sulfur-Containing Compounds Based on the Classical Drude Oscillator Model. J. Comput. Chem. 2010, 31, 23302341. (40) Ibrahim, M. A. A. Molecular Mechanical Study of Halogen Bonding in Drug Discovery. J. Comput. Chem. 2011, 32, 2564-2574. (41) Kolář, 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. (42) 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. (43) Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Comput. J. 1965, 7, 308-313.


34

Page 35 of 39

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


(44) Faller, R.; Schmitz, H.; Biermann, O.; Müller-Plathe, F. Automatic Parameterization of Force Fields for Liquids by Simplex Optimization. J. Comput. Chem. 1999, 20, 1009-1017. (45) Jorgensen, W. L.; Chandrasekhar, J.; Madura J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. (46) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269-10280. (47) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087-1092. (48) Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method. 1. Nonpolar Gases. J. Chem. Phys. 1954, 22, 1420-1426. (49) Jorgensen, W. L.; Tirado-Rives, J. Molecular Modeling of Organic and Biomolecular Systems Using BOSS and MCPRO. J. Comput. Chem. 2005, 26, 1689-1700. (50) Jorgensen, W. L.; Thomas, L. L. Perspective on Free Energy Perturbation Calculations for Chemical Equilibria. J. Chem. Theory Comput. 2008, 4, 869-876. (51) Wolfenden, R. Interaction of Peptide Bond with Solvent Water: A Vapor Phase Analysis. Biochemistry 1978, 17, 201-204. (52) Radzicka, A.; Pedersen, L.; Wolfenden, R. Influences of Solvent Water on Protein Folding: Free energies of Solvation of Cis and Trans Peptides are Nearly Identical. Biochemistry 1988, 27, 4538-4541.


35


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 36 of 39

(53) Abraham, M. H.; Andonian-Haftvan, J.; Whiting, G. S.; Leo, A.; Taft, R. S. HydrogenBonding. 34. The Factors that Influence the Solubility of Gases and Vapours in Water at 298 K, and a New Method for Its Determination. J. Chem. Soc., Perkin Trans. 2 1994, 1777-1791. (54) Dodda, L. S.; Vilseck, J. Z.; Cutrona, K. J.; Jorgensen, W. L. Evaluation of CM5 Charges for Nonaqueous Condensed-Phase Modeling. J. Chem. Theory Comput. 2015, 11, 4273-4282. (55) 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.; et al. Gaussian 09, revision C.01. Gaussian, Inc., Wallingford, CT, 2010. (56) Murphy, R. B.; Pollard, W. T.; Friesner, R. A. Pseudospectral Localized Generalized Moller-Plesset Methods with a Generalized Valence Bond Reference Wave Function: Theory and Calculation of Conformational Energies. J. Chem. Phys. 1997, 106, 5073-5084. (57) Dahlgren, M. K.; Schyman, P.; Tirado-Rives, J.; Jorgensen, W. L. Characterization of Biaryl Torsional Energetics and its Treatment in OPLS All-Atom Force Fields. J. Chem. Inf. Model. 2013, 53, 1191-1199. (58) Jorgensen, W. L.; Jensen, K. P.; Alexandrova, A. N. Polarization Effects for HydrogenBonded Complexes of Substituted Phenols with Water and Chloride Ion. J. Chem. Theory Comput. 2007, 3, 1987-1992. (59) Schyman, P.; Jorgensen, W. L. Exploring Adsorption of Water and Ions on Carbon Surfaces Using a Polarizable Force Field. J. Phys. Chem. Lett. 2013, 4, 468-474.


36

Page 37 of 39

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


(60) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553566. (61) Robertson, M. J.; Tirado-Rives, J.; Jorgensen, W. L. Improved Peptide and Protein Torsional Energetics with the OPLS-AA Force Field. J. Chem. Theory Comput. 2015, 11, 34993509. (62) Mobley, D. L. Experimental and Calculated Small Molecule Hydration Free energies. UC Irvine: Department of Pharmaceutical Sciences, UCI, 2013. (63) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378-6396. (64) Řezáč, J.; Riley, K. E.; Hobza, P. Benchmark Calculations of Noncovalent Interactions of Halogenated Molecules. J. Chem. Theory Comput. 2012, 8, 4285-4292. (65) Voth, A. R.; Hays, F. A.; Ho, P. S. Directing macromolecular conformation through halogen bonds. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 6188-6193. (66) 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. (67) Meyer, F.; Dubois, P. Halogen Bonding at Work: Recent Applications in Synthetic Chemistry and Materials Science. CrystEngComm 2013, 15, 3058-3071. (68) Bruno, I. J.; Cole, J. C.; Edgington, P. R.; Kessler, M.; Macrae, C. F.; McCabe, P.; Pearson, J.; Taylor, R. New Software for Searching the Cambridge Structural Database and


37


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 38 of 39

Visualizing Crystal Structures. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2002, B58, 389-397.

TOC Graphic


38

Page 39 of 39

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


TOC Graphic


Improved Description of Sulfur Charge Anisotropy in OPLS Force

Recommend Documents