Validation and Comparison of Force Fields for Native Cyclodextrins in

Dec 29, 2017 - ... of Force Fields for Native Cyclodextrins in Aqueous Solution. Julia Gebhardt†, Catharina Kleist‡, Sven Jakobtorweihen‡, and N...
0 downloads 0 Views 4MB Size
Download PDF
Article Cite This: J. Phys. Chem. B 2018, 122, 1608−1626

pubs.acs.org/JPCB

Validation and Comparison of Force Fields for Native Cyclodextrins in Aqueous Solution Julia Gebhardt,† Catharina Kleist,‡ Sven Jakobtorweihen,*,‡ and Niels Hansen*,† †

Institute of Thermodynamics and Thermal Process Engineering, University of Stuttgart, D-70569 Stuttgart, Germany Institute of Thermal Separation Processes, Hamburg University of Technology, D-21073 Hamburg, Germany

‡

Downloaded via UNIV OF SOUTH DAKOTA on June 30, 2018 at 16:25:43 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Molecular dynamics simulations of native α-, β-, and γ-cyclodextrin in aqueous solution have been conducted with the goal to investigate the performance of the CHARMM36 force field, the AMBER-compatible q4md-CD force field, and five variants of the GROMOS force field. The properties analyzed are structural parameters derived from Xray diffraction and NMR experiments as well as hydrogen bonds and hydration patterns, including hydration free enthalpies. Recent revisions of the torsional-angle parameters for carbohydrate systems within the GROMOS family of force fields lead to a significant improvement of the agreement between simulated and experimental NMR data. Therefore, we recommend using the variant 53A6GLYC instead of 53A6 and 56A6CARBO_R or 2016H66 instead of 56A6CARBO to simulate cyclodextrins in solution. The CHARMM36 and q4mdCD force fields show a similar performance as the three recommended GROMOS parameter sets. A significant difference is the more flexible nature of the cyclodextrins modeled with the CHARMM36 and q4md-CD force fields compared to the three recommended GROMOS parameter sets.

■

developed for noncyclic (oligo)saccharides.12 The latter study resulted in an AMBER-compatible force field, termed q4mdCD, that reproduces successfully structural, dynamical, and hydrogen bonding aspects of cyclodextrins. Therefore, this force field was used in different cyclodextrin studies,13−16 but also other AMBER versions were used.17,18 For other families of biomolecular force fields such as CHARMM,9,19−23 GROMOS,24−29 or OPLS,30−32 such a validation is lacking, although they are used frequently in the context of cyclodextrin simulations. The parameter set CHARMM36 was applied in some studies,33−36 but also older CHARMM versions were used.37−43 For GROMOS, the parameter set 53A6 is rather often applied.44−51 The OPLS-AA force field was used to model β-cyclodextrin with implicit solvent.52 Studies performed with other force fields are discussed in ref 53. Besides the problem of transferring force field parameters from small (oligo)saccharides to cycloamyloses within a given force field family, it is interesting to compare the molecular properties resulting from different force fields. Historically, force fields such as AMBER,54,55 CHARMM,9,19−23 and GROMOS24−29 have been primarily developed in the context of MD software packages56−58 having the same name as the force field. Today, these force fields are accessible through a variety of other MD packages allowing various force fields to be compared for one system using a single MD software.59 Such a comparison has been conducted for liquids,60 short peptides,61 proteins,62 lipids,63−65 nucleic acids,66,67 and various carbohy-

INTRODUCTION Cyclodextrins (CDs) are cyclic oligosaccharides with a hydrophilic outer surface and a hydrophobic inner cavity. This makes them soluble in water and facilitates the formation of host−guest complexes with hydrophobic organic compounds. Since cyclodextrins are furthermore nontoxic, they are used in diverse applications including the pharmaceutical, food, and textile industries. They can be used as carriers, for example, for drug molecules, vitamins, and fragrances or for masking unwanted odors and flavors. Knowledge of cyclodextrin properties (e.g., complexation processes with other molecules) is therefore of great interest for both research and industry.1−4 On the basis of their numerous practical applications, cyclodextrins have always been the subject of computer simulation studies.5 With today’s computational power, classical molecular dynamics (MD) simulations can routinely generate trajectories on the microsecond time scale for cyclodextrins in explicit solvent.6 These time scales allow one to uncover differences in the force fields used to model cyclodextrins that may be hidden when using much shorter simulations. In addition, it allows one to monitor events that occur relatively rarely, such as the rotation of a pyranose residue around the glycosidic linkage, but may have an impact on the computed values that are compared to experimental data. Force field development for carbohydrates is usually carried out at the mono- or disaccharide level,7−11 leaving open the question of how the parameter sets perform for cyclic oligosaccharides such as cyclodextrins. Indeed, a comparison of several force fields available within the AMBER10 package based on structural parameters of native and permethylated CDs solvated in water revealed issues with the transferability of parameter sets © 2017 American Chemical Society

Received: November 30, 2017 Published: December 29, 2017 1608

DOI: 10.1021/acs.jpcb.7b11808 J. Phys. Chem. B 2018, 122, 1608−1626

Article

The Journal of Physical Chemistry B

Figure 1. Nomenclature and definitions used in this Article for structural parameters and atom labels. The glycosyl dihedral angles ϕ and ψ are defined as ϕ = O5nC1nO1nC4n−1 and ψ = C1nO1nC4n−1C3n−1. Alternate dihedral angles ϕ′ and ψ′ (not shown) are defined as ϕ′ = H1nC1nO1nC4n−1 and ψ′ = C1nO1nC4n−1H4n−1. The transformation between the two conventions is ϕ ≈ ϕ′ + 120° and ψ ≈ ψ′ + 120°.94 The exocyclic dihedral angle ω (O5nC5nC6nO6n) describes the orientation of the primary hydroxyl group O6−H6 relative to the pyranoid ring, whereby the three staggered conformations are referred to as gauche−gauche (gg, ω = −60°), gauche−trans (gt, ω = +60°), and trans−gauche (tg, ω = ± 180°). The dihedral angle ζ defined as O2nC2nC3nO3n describes the orientation of the secondary hydroxyl groups within one glycosidic unit relative to one another. The dihedral angle ξ defined as O2nC1nC4n−1O3n−1 measures the relative orientation of a pair of successive glucose units. The angle δ is defined through three adjacent glucosidic oxygen atoms. Interatomic distances were measured between atoms O1nO1n−1, O2nO3n−1, and O6nO6n−1 of adjacent glucopyranose units.

■

drates.68−73 However, since biomolecular force fields have different functional forms, different combination rules for nonbonded interaction parameters, different ways to treat nonbonded interaction exclusions, different ways to handle long-range electrostatic forces, and different ways of defining torsional-angle interactions, to mention just a few nontrivial technical aspects of the definition of a force field, it cannot be taken for granted that a particular well-defined force-field version is actually implemented as such in software developed by another research group than the one that defined the force field.74 It is therefore of interest to compare not only different force fields using one software package but also different software packages implementing one particular force field.75−78 For the AMBER-compatible parameter set q4md-CD,12 a comparison between its implementation in the GROMACS simulation package and the AMBER10 package, respectively, showed very little differences.16 In light of the ongoing interest in binding thermodynamics of cyclodextrin-based host−guest systems both on the experimental side79−81 as well as on the simulation side,13−18,38,82−85 it is desirable that each of the main biomolecular force fields offers validated parameter sets for cyclodextrins. Recently, we reported differences between binding free energies of primary alcohols to α-cyclodextrin calculated using various flavors of the GROMOS force field.86 Additionally, we found discrepancies in binding free energies calculated with the CHARMM36 and GROMOS 53A6GLYC force fields using the same system.87 These variations may originate from differences in the host− guest and guest−solvent interactions but may also be encoded in the host molecules themselves. Therefore, we present an extensive comparison between properties of native cyclodextrins in aqueous solution modeled using five variants of the GROMOS force field, the CHARMM36 force field, as well as the q4md-CD force field.

COMPUTATIONAL DETAILS Force Fields. The GROMOS force fields used were 53A6,26 53A6GLYC,88 56A6CARBO,10 56A6CARBO_R,11 and 2016H66.89 In the context of pure carbohydrate systems, the force field 53A6 is equivalent to 53A526 as well as recent modifications referred to as 54A728 and 54A829 and nearly identical to the set 45A4.27 The only minor difference between 45A4 and 53A6 force fields was a change in the repulsive coefficient of the Lennard-Jones interactions between atom type OA (used for all sugar oxygen atoms in these two force fields) and nonpolar carbon atoms (intermolecular or intramolecular beyond third covalent neighbors).10 The set labeled 53A6GLYC is based on 53A6 and includes modified torsional-angle parameters for hexopyranoses. The set labeled 56A6CARBO represents a full reparametrization of hexopyranose-based carbohydrates. It is also based on 53A6 and includes three additional atom types for ring atoms. A revised version 56A6CARBO_R was also proposed recently to improve the representation of ring conformational properties in carbohydrate chains. Finally, the most recent GROMOS-compatible parameter set 2016H66 was tested which was parametrized following slightly different principles and procedures than those used for the other sets. However, for the systems studied here, the only differences between 2016H66 and 56A6CARBO_R are the Lennard-Jones parameters for the atom type OA (hydroxyl oxygen) and OE (ether oxygen) that are adopted from the set labeled 53A6OXY.90 Note that for the systems studied in the present work the set 2016H66 is identical to the set termed 56A6CARBO_R+ in our recent study on cyclodextrin−alcohol binding.86 Differences between the force fields in terms of dihedral-angle parameters are reported in Table S13 of the Supporting Information. Topological details of the employed GROMOS cyclodextrin models are provided by Figures S8 and S9 of the Supporting Information. In the official CHARMM36 force field, cyclodextrins are not present. However, all parameters needed are available. Guvench et al.9 introduced the CHARMM parame1609

DOI: 10.1021/acs.jpcb.7b11808 J. Phys. Chem. B 2018, 122, 1608−1626

Article

The Journal of Physical Chemistry B Table 1. Overview of the Simulations Performed in This Worka code A01 A02 A03 A04 A05 A06 A07 A08 A09 A10 A11 A12 A13 A14 A15 A16 A17 A18 B01 B02 B03 B04 B05 B06 B07 B08 G01 G02 G03 G04 G05 G06 G07 G08

starting structure 1CXF_1c 1CXF_2c 2ZYMc 4FEMc C1 d

4

5351 5365 5350 5358 5358 5366 5352 5356 4230 4230 4230 4230 4230 5353 5351 5351 5351 5350

4RERc C1 f

6320 6320 5151 5151 5151 5151 5151 6320

1P2Gc C1 g

7185 7185 5151 5151 5151 5151 5151 7185

4

4

force field

Nwatb

α-Cyclodextrin CHARMM36e 53A6GLYC CHARMM36e 53A6GLYC CHARMM36e 53A6GLYC CHARMM36e 53A6GLYC 53A6 53A6GLYC 56A6CARBO 56A6CARBO_R 2016H66 53A6GLYC CHARMM36e CHARMM36 CHARMM36 q4md-CD β-Cyclodextrin CHARMM36e CHARMM36 53A6 53A6GLYC 56A6CARBO 56A6CARBO_R 2016H66 q4md-CD γ-Cyclodextrin CHARMM36e CHARMM36 53A6 53A6GLYC 56A6CARBO 56A6CARBO_R 2016H66 q4md-CD

water model

software

TIP3P SPC TIP3P SPC TIP3P SPC TIP3P SPC SPC SPC SPC SPC SPC SPC TIP3P TIP3P TIP3P TIP3P

GROMACS 4.6.5 GROMOS11 GROMACS 4.6.5 GROMOS11 GROMACS 4.6.5 GROMOS11 GROMACS 4.6.5 GROMOS11 GROMOS11 GROMOS11 GROMOS11 GROMOS11 GROMOS11 GROMACS 4.6.2 GROMACS 4.6.5 GROMACS 4.6.5 GROMACS 5.1.1 GROMACS 4.6.7

TIP3P TIP3P SPC SPC SPC SPC SPC TIP3P

GROMACS 4.6.5 GROMACS 5.1.1 GROMOS11 GROMOS11 GROMOS11 GROMOS11 GROMOS11 GROMACS 4.6.7

TIP3P TIP3P SPC SPC SPC SPC SPC TIP3P

GROMACS 4.6.5 GROMACS 5.1.1 GROMOS11 GROMOS11 GROMOS11 GROMOS11 GROMOS11 GROMACS 4.6.7

For the simulations of α-cyclodextrin using the 53A6GLYC force field (A10 and A14), different schemes to calculate the long-range electrostatics interactions were investigated, as detailed in the text. bNumber of water molecules. c1CXF,95 2ZYM,96 4FEM,97 4RER,98 and 1P2G99 are Protein Data Bank IDs, where 1CXF contains two α-cyclodextrin molecules in the asymmetric unit. dInitial structure taken from previous work.100 All residues in the 4C1 conformation. eAlternative cutoff set compared to other CHARMM36 simulations with GROMACS; see main text. fInitial structure based on coordinates obtained from the Automated Topology Builder (ATB) and Repository101 after energy minimization in vacuum using the 53A6GLYC force field. All residues in the 4C1 conformation. gInitial structure based on X-ray data102 after energy minimization in vacuum using the 53A6GLYC force field. All residues in the 4C1 conformation. a

Simulation Parameters. MD simulations were performed for systems consisting of one cyclodextrin molecule in a box of water. Starting coordinates were taken from crystal structures or previous work, respectively, as indicated in Table 1. All MD simulations were performed under minimum image periodic boundary conditions based on cubic boxes. The equations of motion were integrated using the leapfrog algorithm103 with a time step of 2 fs. The simulations employing the GROMOS11 program package58,104,105 were carried out with the latest release version.106 Bond lengths and the bond angle of water molecules were constrained by applying the SHAKE algorithm107 with a relative geometric tolerance of 10−4. The center of mass motion of the computational box was removed every 20 ps. Equilibration of the system was started by generating velocities

ters for glycosidic linkages which were used in the present work. The oxygen atom O1 (see Figure 1) at the glycosidic linkage was modeled with the CHARMM atom type OC301 and a partial charge of −0.36; further details are shown in Figure S10 of the Supporting Information. On the basis of the GLYCAM0491−93 force field, Cézard et al.12 introduced the q4md-CD force field for native and modified cyclodextrins. Among the investigated force fields in this study, it is the only one specially designed for cyclodextrins. A comparison of q4md-CD to other GLYCAM and AMBER versions is provided by Cézard et al.12 Therefore, in the present work, we focus on q4md-CD as an AMBER-compatible force field. An overview over the simulations conducted in the present work is given in Table 1. 1610

DOI: 10.1021/acs.jpcb.7b11808 J. Phys. Chem. B 2018, 122, 1608−1626

Article

The Journal of Physical Chemistry B

barostat132,133 by isotropic coupling using a coupling constant of τp = 5 ps. If not otherwise mentioned, for all simulations with the CHARMM force field, the same cutoff and electrostatic settings were used. However, to rule out an influence of these settings, one simulation was carried out with alternative settings. In the following, the alternative settings are given in brackets. The neighbor list cutoff was 1.2 nm (1 nm), and the list was updated every 10th step. The Lennard-Jones potential was switched off between 1.0 (0.8) and 1.2 nm; long-range corrections were not used. Hence, only for the alternative settings, a twin-range cutoff with 1/1.2 nm was used for the Lennard-Jones potential, where the long-range forces were updated during neighbor searching. Electrostatic interactions were truncated at 1.2 nm (1.0 nm) with interactions beyond the cutoff taken into account with the PME method.127,128 For this method, a grid spacing of 0.18 nm (0.15 nm) with fourth (sixth) order B-spline interpolation was employed. While parametrizing the CHARMM force field, a force-based switching function for the Lennard-Jones interactions was applied,76 which is now available in recent GROMACS versions. This method is actually recommended for lipid bilayer simulations with the CHARMM force field.76 Therefore, we also carried out simulations with the CHARMM36 force field using GROMACS version 5.1.1134 employing the forcebased switching function and the Verlet cutoff scheme. Prior to the production simulations, the structures were equilibrated by energy minimization in a vacuum where the heavy atoms were restrained with 1000 kJ mol−1 nm−2. Thereafter, water was added followed by a NPT equilibration simulation of 200 ps with restraints on the heavy atoms of cyclodextrin (1000 kJ mol−1 nm−2). These equilibration simulations used a weak pressure coupling;108 all other parameters are the same as those for the production runs. For the simulations with the q4md-CD force field,12 the same equilibration steps as those for the CHARMM36 force field simulations were carried out. Zhang et al.16 have validated the performance of the q4md-CD force field with GROMACS, and we used the same simulation parameters, with the exception of the temperature coupling method. That is, the Lennard-Jones potential was cut off at 1.4 nm with no dispersion corrections, whereas the electrostatic interactions were cut off at 1.0 nm in combination with the PME method. The neighbor list cutoff was 1.0 nm with a list update every 10th step. All bonds were constrained, and a time step of 2 fs was used. Temperature and pressure coupling schemes were the same as those used for the simulations with the CHARMM force field. The analyzed simulations were carried out for 500 ns while saving configurations every 10 ps. Our results correspond well to the results of Zhang et al.16 and therefore also to the results obtained with AMBER10.12 Hydration Free Enthalpy. The calculation of the hydration free enthalpy using the GROMOS software relied on thermodynamic integration135 over the average of the derivative of the Hamiltonian with respect to a scaling parameter λ applied to the interactions of CD with the solvent, from full interaction at λ = 0 to no interaction at λ = 1. The λdecoupling was chosen to be linear and involves a soft-core scaling136 with αLJ = 0.5 and αCRF = 0.5 nm2 for the LennardJones and electrostatic interactions, respectively. The solute− solvent interactions were gradually deactivated using 21 equispaced λ-points with a sampling period of 3.6 ns after 0.4 ns equilibration at each λ-point. Free enthalpies calculated in this way are based on a unit molarity (1 mol dm−3) standard

from a random Maxwell−Boltzmann distribution at 60 K. Position restraints were introduced on the solute with a force constant of 2.5 × 104 kJ mol−1 nm−2. In five 20 ps simulations, the temperature was raised in steps of 60 K, until the final temperature of 298 K was reached. With every increase in temperature, the force constant of the position restraints was reduced by a factor of 10. In the final equilibration step at 298 K, no position restraints were used. All production simulations were performed at constant pressure and temperature for 100− 200 ns. The temperature was maintained at 298 K by weak coupling to an external bath108 with a relaxation time of 0.1 ps. Solute and solvent were coupled to separate heat baths. The pressure was calculated using a group-based virial and held constant at 1 atm using weak coupling of the atomic coordinates and box dimensions to a pressure bath with a relaxation time of 0.5 ps108 and an isothermal compressibility of 7.513 × 10−4 (kJ mol−1 nm−3)−1 representing the value for water.109 van der Waals and electrostatic interactions were handled using the Lennard-Jones potential110−112 and the Barker−Watts reaction field scheme113 within a triple-range cutoff approach114 applied on the basis of distances between charge group centers,115 with short- and long-range cutoff radii of 0.8 and 1.4 nm, respectively, and an update frequency of five time steps for the short-range pair list and intermediate-range interactions. The reaction field scheme was applied with εRF of 61 for the relative permittivity of the dielectric continuum surrounding the cutoff sphere, as appropriate116 for the SPC water model.117 The reaction-field self-term and excludedatom-term contributions118 to the energy, forces, and virial were included as described in ref 119. For α-cyclodextrin modeled with the 53A6GLYC force field,88 two additional simulations were carried out with the GROMOS program using either the Barker−Watts reaction field scheme113 in combination with an atom-based cutoff or the particle−particle−particle−mesh (P3M) method120,121 for the treatment of electrostatic interactions. The latter was applied with120,122,123 a spherical hat charge-shaping function of width 1.0 nm, a triangular shaped cloud assignment function, a finitedifference (FD) scheme of order 2, and a grid spacing of about 0.08 nm. Three additional simulations for α-cyclodextrin modeled with the 53A6GLYC force field88 were carried out using the GROMACS program package124 (version 4.6.2). All bond lengths were kept fixed at their equilibrium values using either SETTLE125 (for water) or LINCS126 (for CD). The bond lengths were constrained by application of the LINCS procedure using a lincs order of 4. The number of iterations to correct for rotational lengthening in LINCS was set to 4. The long-range electrostatics was treated either by the smooth particle-mesh Ewald (PME) summation,127,128 using different grid spacings of 0.15 and 0.09 nm, or the Barker−Watts reaction field scheme,113 respectively. The simulations using the CHARMM force field were conducted with the GROMACS program package 75,124 (version 4.6.5). Water was modeled with the CHARMM TIP3P force field, which is a slight modification19,129 of the original TIP3P model,130 and was kept rigid with the SETTLE algorithm.125 For cyclodextrin, the LINCS algorithm126 was used to constrain all bonds to hydrogen atoms. Each simulation was carried out for 500 ns saving configurations every 10 ps. The temperature was kept constant at 298 K with the Nosé− Hoover thermostat131 using a coupling constant of τT = 1 ps. The pressure was set to 1 bar with the Parrinello−Rahman 1611

DOI: 10.1021/acs.jpcb.7b11808 J. Phys. Chem. B 2018, 122, 1608−1626

Article

The Journal of Physical Chemistry B state for the gas and solution phases.137 To improve configurational sampling of the molecule as it is decoupled from its surroundings, a stochastic thermostat138 was employed with a friction coefficient of 91 ps−1 mimicking water.139 Note that a choice for a specific value of the friction coefficient affects the dynamics of the system but not the averages of the calculated thermodynamic properties. The Hamiltonian derivatives were stored every 0.5 ps. Alchemical simulations with the CHARMM36 force field were carried out with GROMACS version 5.1.1134 and with GROMACS version 4.6.5124 for the q4md-CD force field. In both cases, the decoupling parameters were the same. Twentysix λ-points were used (9 for switching off electrostatic interactions followed by 17 states for switching off LennardJones interactions) with a soft-core potential with αsc = 0.5, σsc = 0.3, and a soft core power of 1. Each λ-point was simulated for 7 ns, where only the last 6 ns was used for analysis. The Hamiltonian derivatives were stored every 0.4 ps. The first λpoint (interacting cyclodextrin) was started from the final snapshot of the corresponding 500 ns simulation (see above), and the other λ-points were started from its preceding simulations after it had reached 1 ns. For the leapfrog stochastic dynamics integrator, a friction coefficient of 1 ps−1 was used. For all force fields, the integration was performed using the trapezoidal rule. Errors in the individual ensemble averages for the Hamiltonian derivative were calculated by block-averaging,140 and propagated into an error estimate for the free enthalpy according to the trapezoidal rule. Trajectory Analysis. The trajectories were analyzed in terms of geometrical parameters illustrated in Figure 1, structural patterns such as ring puckering and rotamer populations, hydrogen bonding patterns, and partial molar volumes. Hydrogen bonds were identified according to geometric critera: A two-center hydrogen bond was assumed to exist if the hydrogen−acceptor distance is smaller than 0.25 nm and the donor−hydrogen−acceptor angle is larger than 135°. Occurrences of three-center hydrogen bonds are defined for a donor atom D, hydrogen atom H, and two acceptor atoms A1 and A2 if (i) the distances H−A1 and H−A2 are smaller than 0.27 nm; (ii) the angles D−H−A1 and D−H−A2 are larger than 90°; (iii) the sum of the angles D−H−A1, D−H−A2, and A1−H− A2 is larger than 340°; and (iv) the dihedral angle defined by the planes through the atoms D−A1−A2 and H−A1−A2 is smaller than 15°.141 3 J coupling constants were calculated using the Karplus relation142 3

J = a cos2 β + b cos β + c

Figure 2. Karplus curves describing the relation between 3J coupling and the torsional angles C−O−C−H, that is, ϕ′ and ψ′ (left panel) and the torsional angles Hi−Ci−Oi−Hi with i = 2, 3, 6 (right panel). The curves shown for 3JCOCH couplings are based on the calibration constants reported by Mulloy et al.144 (red), Tvaroška et al.143 (black), Tvaroška145 (blue), Tafazzoli et al.146 (green), Säwen et al.147 (purple), and Yang et al.148 (orange). The curves shown for 3JOHCH couplings are based on the calibration constants reported by Fraser et al.149 and Zhao et al.150 (red, eq 2 in ref 150; blue, eq 4 in ref 150).

The assignment of configurations to the three distinct classes of conformers relied on intervals of ±60° centered at −60, +60, or 180° for the attribution of the ω-dihedral angle to the gg, gt, or tg wells, respectively. The puckering parameters for pyranoid rings were calculated according to Cremer and Pople.151 With this convention, a perfect 4C1 chair conformation of the six-membered ring has an azimuthal angle θ of 0°.152 Following previous work,153−155 conformations of the pyranoid rings were classified as chair (C), boat (B), or skew-boat (S) along with indices n and m representing the atom numbers of two atoms pointing upward (superscript) or downward (subscript), where the “upper” face of the pyranoid ring is defined by the “screwdriver rule” for a counterclockwise rotation (i.e., following ring atoms in order of decreasing atom numbers). The hydration pattern was investigated on the basis of the number of direct hydrogen bonds between cyclodextrin and water as well as in terms of the radial and cylindrical pair correlation function with respect to the center of the cyclodextrin cavity. The radial distribution function was calculated in the usual way as140 g (r ) =

1 ρW

NW

∑ δ(r − |rCi|) (2)

i=1

where ρW is the number density of water and rCi the vector connecting the center of the cyclodextrin cavity to water molecule i. The cylindrical distribution function was calculated as156,157 g (r , z ) =

(1)

where a, b, and c depend on the type of coupling. To calculate the H1′−C4 and H4−C1′ constants from the ϕ′ and ψ′ dihedral angles, respectively, we have used a = 5.7 Hz, b = −0.6 Hz, and c = 0.5 Hz according to Tvaroška et al.143 To calculate the OH−CH coupling constants of hydroxy protons, we used eq 1 with a = 10.4 Hz, b = −1.5 Hz, and c = 0.2 Hz according to Fraser et al.149 Alternative parametrizations of the Karplus curve are displayed in Figure 2. The differences in the curves suggest an uncertainty of about 1 Hz for the relation between the ϕ′ and ψ′ angles, respectively, and the corresponding coupling constant. For the OH−CH couplings, the uncertainty may even exceed 1 Hz. For the analysis of the rotamer distribution, three staggered conformations of the hydroxymethyl group were considered.

1 ρW δ(r −

NW

∑ δ(z − nzrCi) i=1

|rCi|2 − (nzrCi)2 ) (3)

where the z-axis is aligned along the symmetry axis of the cyclodextrin, the unit vector nz points toward the primary rim, and the origin (r = z = 0) is defined as the center of geometry of the glycosidic oxygen atoms. Because there is no unambiguous discrimination between water molecules that are inside and those that are outside the cavity and different criteria are used in the literature leading to different conclusions,12,158,159 such an analysis was not attempted in the present work. Partial molar volumes were obtained from the difference in average volume between two aqueous systems involving the 1612

DOI: 10.1021/acs.jpcb.7b11808 J. Phys. Chem. B 2018, 122, 1608−1626

Article

The Journal of Physical Chemistry B

Figure 3. Occurrence of glucopyranose ring conformation types154,155 in α-cyclodextrin simulations.

Figure 4. Time series of the flip dihedral angle ξ for α-cyclodextrin. The different colors denote the different glucopyranose units.

force field did not affect the results of the simulations. Therefore, if not otherwise mentioned, we show and discuss in the following the results of simulations A09-A13, A17, and A18 for α-cyclodextrin, simulations B02−B08 for β-cyclodextrin, and simulations G02−G08 for γ-cyclodextrin (see Table 1). Glucopyranose Ring Conformations. Figure 3 shows the occurrence of glucopyranose ring conformation types in αcyclodextrin. As described before, possible conformations are divided into 14 types.153−155 Conformations which cannot be assigned to any of these types are referred to as unassigned. Three of the GROMOS force fields, 53A6GLYC, 56A6CARBO_R, and 2016H66, the CHARMM36 force field, as well as q4mdCD show a very similar occurrence of conformation types. Throughout a vast majority of the simulation time, the glucopyranose units adopt the chair conformation 4C1. In the simulation with GROMOS force field 53A6, however, the

same number of water molecules, either in the absence or in the presence of one cyclodextrin molecule.

■

RESULTS AND DISCUSSION

The five different starting structures (see Table 1) used for the CHARMM36 simulations lead to very similar results. For the GROMOS force fields 53A6GLYC, 56A6CARBO_R, and 2016H66, the starting structure also did not affect the outcome of the simulations because the 4C1 conformation was readily adopted for all residues during energy minimization. These results are in agreement with the picture derived from NMR data.160,161 Therefore, in the following sections, only those results are reported which are based on starting coordinates in which all glucopyranose units are in the 4C1 conformation. Also, the alternative cutoff set in the case of CHARMM36 and using GROMACS instead of GROMOS in the case of the 53A6GLYC 1613

DOI: 10.1021/acs.jpcb.7b11808 J. Phys. Chem. B 2018, 122, 1608−1626

Article

The Journal of Physical Chemistry B

Table 2. Structural Properties of α-Cyclodextrin Obtained from MD Simulations Using Different Force Fields and from Experimenta property

53A6

56A6CARBO

53A6GLYC

56A6CARBO_R

2016H66

CHARMM36

q4md-CD

experiment

ϕ dihedral angle (deg) ψ dihedral angle (deg) ζ dihedral angle (deg) δ angle (deg) dO1 (nm) dO23 (nm) dO6 (nm) circularityb puckering amplitude Q (nm) θ azimuthal angle (deg)

74.2 160.7 148.8 119.4 0.482 0.620 0.507 0.875 0.052 164.0

80.0 167.4 116.6 118.7 0.451 0.729 0.576 0.766 0.055 110.0

115.3 118.6 79.7 119.7 0.434 0.287 0.566 0.925 0.056 8.9

117.2 116.0 80.5 119.6 0.421 0.301 0.551 0.894 0.068 14.2

115.6 118.2 52.4 119.6 0.424 0.304 0.551 0.916 0.058 12.6

107.5 129.0 56.2 119.4 0.426 0.321 0.531 0.901 0.057 10.3

111.7 126.3 69.6 119.4 0.417 0.303 0.527 0.911 0.055 12.1

107.4,c 109.2d 130.7,c 128.8d 68.1d 119.9d 0.4235d 0.306,c 0.2981d 0.596d 0.88−0.93e 0.0577c 5.7c

a Experimental data are for α-cyclodextrin in crystalline form. The simulations using the various GROMOS force fields were conducted with the standard settings; see main text. bRatio of the smallest to the largest distance between any pair of glucose O1 atoms that lie across the ring from each other. cReference 94. dReference 161. eReference 168.

Figure 5. Distribution of the flip dihedral angle ξ for α-cyclodextrin. The different colors denote the different glucopyranose units.

for the whole trajectory are 110.35 and 126.22° (see also Table S1). If the frames showing a rotated unit are excluded, these values change to 112.54 and 125.09°. For γ-cyclodextrin, more frames have rotated units (68.50%) than frames without any rotated unit. For this simulation (G02), ϕ calculated over the whole trajectory is 100.93° (see also Table S2) and excluding the frames with rotated units 114.87°. The experimental value is 108.9°. ψ over the whole trajectory is 131.15°, with excluding frames 122.82° and the experimental value is 127.1°. Hence, a comparison to these experimental values alone cannot explain if the frequent glucose unit rotations in the CHARMM γcyclodextrin represents a realistic behavior. Simulations with the GROMOS force fields show less conformational variability. With the force field 53A6GLYC, almost no rotations occur, while for 56A6CARBO_R and 2016H66 occasional tumbling events take place for α-cyclodextrin. Consequently, for force field 53A6GLYC, the distribution of flip dihedral angle values is more narrow than that for 56A6CARBO_R and 2016H66; see Figure 5. For β-cyclodextrin, 40% of the frames with force field 56A6CARBO_R have at least one glucose unit rotated (see Figure S3a). For γ-cyclodextrin, more frames have rotated units and force field 2016H66 also shows some rotations (see Figure S3b). With force field 53A6GLYC, no rotations occur for β- and γ-cyclodextrin. With force fields 53A6 and 56A6CARBO, all frames show rotated glucose units. For 53A6, the peak of the flip-dihedral angle distribution is shifted to −90°, while for 56A6CARBO every glucose unit exhibits a different distribution. However, as shown below, these two force fields lead to considerable disagreement with experimentally derived structural data. Tumbling is influenced by two main contributions, the dihedral angle potentials of the glycosidic linkages and the strength of the intramolecular hydrogen bond network within the cyclodextrins. Both will also be discussed below. For future force-field refinement, experimental evidence of tumbling

dominant conformation type is the inverted chair conformation 1 C 4. With the force field 56A6 CARBO, the 4C 1 and 1C 4 conformations occur with alternating population between neighboring units. Some other conformation types also occur but to a significantly smaller extent. As shown in the Supporting Information, the observations made for β- and γ-cyclodextrin are very similar (see Figures S1 and S2). Tumbling of Glucopyranose Units. Especially in the simulations with the CHARMM36 and q4md-CD force fields, a glucopyranose unit may turn upside down. For α- and β-CD, these configurations are not persistent (see Figures 4 and S3a, respectively), while, for γ-CD, the flipped state remains stable over several nanoseconds (see Figure S3b). It is known that for the hexahydrate form of α-cyclodextrin one glucose unit may rotate out of register with the other five.161,162 Unfortunately, we do not have any experimental information on how frequently these rotations occur in native cyclodextrins in solution. In cyclodextrin derivatives carrying altropyranose units or permethylated ones, tumbling has been observed in experiment163−165 as well as simulation using the CHARMM36 force field.35 Of course, if one glucose unit is turned, some properties at these instances change (e.g., the RMSD of α-cyclodextrin and the angle ξ). We define a rotated glucose unit by ξ < −70° or ξ > 70°. In the 500 ns CHARMM36 simulation of α-cyclodextrin started from 4C1 (simulation A17), 6.17% of the frames show a glucose rotation; see Figure 4. The averaged properties shown in this Article would change very little if we exclude these frames from the analysis. For example, ϕ and ψ would be 108.59 and 128.46° without rotated units compared to 107.53 and 129.00° (see also Table 2). Since the simulations show more glucose unit rotations for β-cyclodextrin and γ-cyclodextrin, the influence on the properties is larger. For βcyclodextrin (simulation B02), 12.97% of the frames have at least one unit rotated. For this simulation, the ϕ and ψ values 1614

DOI: 10.1021/acs.jpcb.7b11808 J. Phys. Chem. B 2018, 122, 1608−1626

Article

The Journal of Physical Chemistry B

Table 3. NMR Coupling Constants (3J), Proton−Proton Distances (rNOE), and Relative Populations of the Three Staggered Conformers of the Hydroxymethyl Group (ω) for α-Cyclodextrin Obtained from MD Simulations with Different Force Fields and from Experimenta property 3

c

JC,H (Hz) JOH,CHc (Hz)

3

rNOEg (nm) ωtor. pop. (%)

H1′−C4 H4−C1′ O(2)H O(3)H O(6)H H1−H2 H1′−H4 gg gt tg

53A6

56A6CARBO

53A6GLYC

56A6CARBO_R

2016H66

CHARMM36

q4md-CD

experiment

3.3 3.3 4.4 3.9 4.8 0.24 0.22 14 85 1

3.4 4.6 4.5 5.0 5.4 0.23 0.21 27 41 32

5.4 5.4 5.2 2.6 4.9 0.24 0.21 62 38 0

5.3 5.4 4.6 3.2 5.6 0.24 0.19 53 45 2

5.3 5.4 3.7 3.1 5.6 0.24 0.19 57 42 1

5.0 4.9 8.0 5.4 6.6 0.25 0.22 55 44 1

5.2 5.0 6.4 3.3 5.9 0.24 0.21 47 52 1

5.2,d 6.4e 5.0d 6.6f