Validation and Comparison of Force Fields for Native Cyclodextrins in

Subscriber access provided by AUSTRALIAN NATIONAL UNIV

Article

Validation and Comparison of Force Fields for Native Cyclodextrins in Aqueous Solution Julia Gebhardt, Catharina Kleist, Sven Jakobtorweihen, and Niels Hansen J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b11808 • Publication Date (Web): 29 Dec 2017 Downloaded from http://pubs.acs.org on December 31, 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 47 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

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 E-mail: [email protected]; [email protected] Phone: +49 (0)40 42878 2411; +49 (0)711 685 66112. Fax: +49 (0)40 42878 4072; +49 (0)711 685 66140

1

ACS Paragon Plus Environment

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

Abstract Molecular dynamics simulations of native α-, β- and γ-cyclodextrin in aqueous solution have been conducted with the goal to investigate the performance of the CHARMM36, the AMBER-compatible q4md-CD, and five variants of the GROMOS force field. The properties analyzed are structural parameters derived from X-ray 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 to use the variant 53A6GLYC instead of 53A6 and 56A6CARBO R or 2016H66 instead of 56A6CARBO to simulate cyclodextrins in solution. The CHARMM36 and q4md-CD 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.

Introduction Cyclodextrins (CDs) are cyclic oligosaccharids 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 non-toxic 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 Based on their numerous practical applications, cyclodextrins have always been the subject of computer simulation studies. 5 With todays’ computational power classical molecular 2


Page 2 of 47

Page 3 of 47 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


dynamics (MD) simulations can routinely generate trajectories on the microsecond time scale for cyclodextrins in explicit solvent. 6 These time scales allow to uncover differences in the force fields used to model cyclodextrins that may be hidden when using much shorter simulations. In addition, it allows 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 developed for non-cyclic (oligo)saccharides. 12 The latter study resulted in an AMBER-compatible force field, termed q4md-CD, 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 GROMOS 24–29 have been primarily developed in the context of MD software packages 56–58 having the same name as the force field. Today, these force fields are accessible through a variety of other MD packages allowing to compare various force

3



fields 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 carbohydrates. 68–73 However, since biomolecular force fields have different functional forms, different combination rules for non-bonded interaction parameters, different ways to treat non-bonded interaction exclusions, different ways to handle long-range electrostatic forces and different ways of defining torsional-angle interactions, to mention just a few non-trivial 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 hostguest systems both on the experimental side 79–81 as well as on the simulation side, 13–18,82–85 it is desirable that each of the main biomolecular force fields offer 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 the GROMOS 53A6GLYC force field 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.

4


Page 4 of 47

Page 5 of 47 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


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 53A5 26 as well as recent modifications referred to as 54A7 28 and 54A8 29 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 neighbours). 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 used for the other sets. However, for the systems studied here the only difference between 2016H66 and 56A6CARBO R are the LennardJones 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 Figs 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 parameters for glycosidic linkages which were used in the present work. The oxygen atom O1

5



Page 6 of 47

OH O

HO

OH

O

O

δ

HO OH O

HO O OH

HO6

HO

O6

O OH

O6

ω

C6

O HO

HO6

dO6

OH

C6

C5

O5

C3

C2

φ

C5

O5

C3

C2

OH

C4

O

dO1

O OH OH

O

O3 O

O1

C1

ψ

OH

OH

ζ

flip ξ O2

OH O

HO3

HO2

n

HO

C4

dO23

C1

O3

O2

HO3

HO2

O1

n-1

Figure 1: Nomenclature and definitions used in this article for structural parameters and atom labels. The glycosyl dihedral angles φ and ψ are defined as φ = O5n C1n O1n C4n−1 and ψ = C1n O1n C4n−1 C3n−1 . Alternate dihedral angles φ′ and ψ ′ (not shown) are defined as φ′ = H1n C1n O1n C4n−1 and ψ ′ = C1n O1n C4n−1 H4n−1 . The transformation between the two conventions is φ ≈ φ′ + 120◦ and ψ ≈ ψ ′ + 120◦ . 94 The exocyclic dihedral angle ω (O5n C5n C6n O6n ) describes the orientation of the primary hydroxyl group O6-H6 relative to the pyranoid ring, whereby the three staggered conformations are referred to as gauchegauche (gg, ω = −60◦ ), gauche-trans (gt, ω = +60◦ ), and trans-gauche (tg, ω = ±180◦ ). The dihedral angle ζ defined as O2n C2n C3n O3n describes the orientation of the secondary hydroxyl groups within one glycosidic unit relative to one another. The dihedral angle ξ defined as O2n C1n C4n−1 O3n−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 O1n O1n−1 , O2n O3n−1 and O6n O6n−1 of adjacent glucopyranose units. (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 Fig. S10 of the Supporting Information. Based on the GLYCAM04 91–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 AMBER-compatible force field. An overview over the simulations conducted in the present work is given in Table 1.

6


Page 7 of 47 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 1: Overview of the simulations performed in this worka code

Starting structure

Nwat b

force field water model software α-cyclodextrin A01 1CXF 1c 5351 CHARMM36e TIP3P GROMACS 4.6.5 A02 5365 53A6GLYC SPC GROMOS11 A03 1CXF 2c 5350 CHARMM36e TIP3P GROMACS 4.6.5 A04 5358 53A6GLYC SPC GROMOS11 A05 2ZYMc 5358 CHARMM36e TIP3P GROMACS 4.6.5 A06 5366 53A6GLYC SPC GROMOS11 A07 4FEMc 5352 CHARMM36e TIP3P GROMACS 4.6.5 A08 5356 53A6GLYC SPC GROMOS11 A09 4 C1 d 4230 53A6 SPC GROMOS11 A10 4230 53A6GLYC SPC GROMOS11 A11 4230 56A6CARBO SPC GROMOS11 GROMOS11 A12 4230 56A6CARBO R SPC A13 4230 2016H66 SPC GROMOS11 A14 5353 53A6GLYC SPC GROMACS 4.6.2 A15 5351 CHARMM36e TIP3P GROMACS 4.6.5 A16 5351 CHARMM36 TIP3P GROMACS 4.6.5 A17 5351 CHARMM36 TIP3P GROMACS 5.1.1 A18 5350 q4md-CD TIP3P GROMACS 4.6.7 β-cyclodextrin B01 4RERc 6320 CHARMM36e TIP3P GROMACS 4.6.5 B02 4 C1 f 6320 CHARMM36 TIP3P GROMACS 5.1.1 B03 5151 53A6 SPC GROMOS11 B04 5151 53A6GLYC SPC GROMOS11 B05 5151 56A6CARBO SPC GROMOS11 B06 5151 56A6CARBO R SPC GROMOS11 B07 5151 2016H66 SPC GROMOS11 B08 6320 q4md-CD TIP3P GROMACS 4.6.7 γ-cyclodextrin G01 1P2Gc 7185 CHARMM36e TIP3P GROMACS 4.6.5 4 g G02 C1 7185 CHARMM36 TIP3P GROMACS 5.1.1 G03 5151 53A6 SPC GROMOS11 G04 5151 53A6GLYC SPC GROMOS11 G05 5151 56A6CARBO SPC GROMOS11 GROMOS11 G06 5151 56A6CARBO R SPC G07 5151 2016H66 SPC GROMOS11 G08 7185 q4md-CD TIP3P GROMACS 4.6.7 a 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. b Number of water molecules. c 1CXF, 95 2ZYM, 96 4FEM, 97 4RER, 98 and 1P2G 99 are protein databank IDs, where 1CXF contains two α-cyclodextrin molecules in the asymmetric unit. d Initial structure taken from previous work. 100 All residues in 4 C1 conformation. e Alternative cut-off set compared to other CHARMM36 simulations with GROMACS, see main text. f Initial structure based on coordinates obtained from the Automated Topology Builder (ATB) and Repository 101 after energy minimization in vacuum using the 53A6GLYC force field. All residues in 4 C1 conformation. g Initial structure based on X-ray data 102 after energy minimization in vacuum using the 53A6GLYC force field. All residues in 4 C1 conformation.

7



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 algorithm 103 with a timestep of 2 fs. The simulations employing the GROMOS11 program package 58,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 algorithm 107 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 from a random MaxwellBoltzmann 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 to 200 ns. The temperature was maintained at 298 K by weak coupling to an external bath 108 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 ps 108 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 potential 110–112 and the Barker-Watts reaction field scheme 113 within a triple-range cutoff approach 114 applied on the basis of distances between charge group centers, 115 with shortand 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 8


Page 8 of 47

Page 9 of 47 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


scheme was applied with εRF of 61 for the relative permittivity of the dielectric continuum surrounding the cutoff sphere, as appropriate 116 for the SPC water model. 117 The reactionfield self-term and excluded-atom-term contributions 118 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 scheme 113 in combination with an atom-based cutoff or the particle-particle-particle-mesh (P3M) method 120,121 for the treatment of electrostatic interactions. The latter was applied with 120,122,123 a spherical hat charge-shaping function of width 1.0 nm, a triangular shaped cloud assignment function, a finite-difference (FD) scheme of order two and a grid spacing of about 0.08 nm. Three additional simulations for α-cyclodextrin modeled with the 53A6GLYC force field 88 were carried out using the GROMACS program package 124 (version 4.6.2). All bond lengths were kept fixed at their equilibrium values using either SETTLE 125 (for water) or LINCS 126 (for CD). The bond lengths constrained by application of the LINCS procedure are 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 modification 19,129 of the original TIP3P model, 130 and was kept rigid with the SETTLE algorithm. 125 For cyclodextrin the LINCS algorithm 126 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 thermostat 131 using a coupling constant of τT =1 ps. The pressure was set to 1 bar with the Parrinello-Rahman barostat 132,133 by isotropic coupling using a coupling constant

9



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 4th (6th) order B-spline interpolation was employed. While parametrizing the CHARMM force field a force-based switching function for the Lennard-Jones interactions were 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.1 134 employing the force-based switching function and the Verlet cutoff-scheme. Prior to the production simulations the structures were equilibrated by energy minimization in vacuum where the heavy atoms were restrained with 1000 kJ mol−1 nnm−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 for the production runs. For the simulations with the q4md-CD force field 12 the same equilibration steps as 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-

10


Page 10 of 47

Page 11 of 47 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


Jones potential was cutoff at 1.4 nm with no dispersion corrections, whereas the electrostatic interactions were cutoff 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 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 integration 135 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 scaling 136 with αLJ = 0.5 and αCRF = 0.5 nm2 for the Lennard-Jones 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 state for the gas and solution phases. 137 To improve configurational sampling of the molecule as it is decoupled from its surroundings a stochastic thermostat 138 was employed with a friction coefficient of 91ps−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.1 134 and with GROMACS version 4.6.5 124 for the q4md-CD force field. In both cases the decoupling parameters were the same. 26 λ-points were used (nine for switching off electrostatic interactions followed by 17 states for switching off Lennard-Jones 11



Page 12 of 47

interactions) with a soft-core potential with αsc =0.5, σsc =0.3, and a soft core power of one. Each λ-point were simulated for 7 ns, where only the last 6 ns were 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 has reached 1 ns. For the leap-frog stochastic dynamics integrator a friction coefficient of 1ps−1 were 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 blockaveraging, 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 Fig. 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 relation, 142 3

J = a cos2 β + b cos β + c

12


(1)


−180.0

JOH,CH / Hz

10.0 8.0 6.0 4.0 2.0 0.0

3

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

JC,H / Hz

Page 13 of 47

−90.0 0.0 90.0 C−O−C−H / deg

180.0

16.0 14.0 12.0 10.0 8.0 6.0 4.0 2.0 0.0 −2.0 −180.0

−90.0 0.0 90.0 H−C−O−H / deg

180.0

Figure 2: The Karplus curves describing the relation between 3 J-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 3 J COCH -couplings are based on the calibration constants reported by Mulloy et al. 144 (red), Tvaroˇska et al. 143 (black), Tvaroˇska 145 (blue), Tafazzoli et al. 146 (green), Säwen et al. 147 (purple) and Yang et al. 148 (orange). The curves shown for 3 J OHCH -couplings are based on the calibration constants reported by Fraser et al. 149 and Zhao et al. 150 (red, eq. 2; blue, eq. 4 in ref. 150). 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ˇska 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. 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 4 C1 chair conformation of the six-membered ring has an azimuthal angle θ of 0◦ . 152 Following previous work 153–155 conformations of the pyranoid

13



Page 14 of 47

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 upwards (superscript) or downwards (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 based on the number of direct hydrogen bonds between cyclodextrin and water as well as in terms of the radial and a cylindrical pair correlation function with respect to the center of the cyclodextrin cavity. The radial distribution function was calculated in the usual way as 140 W X 1 N δ(r − |rCi |) g(r) = ρW i=1

*

+

(2)

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 as 156,157 q W X 1 N 2 2 δ(z − nz rCi )δ r − |rCi | − (nz rCi ) g(r, z) = ρW i=1 *

+

(3)

where the z-axis is aligned along the symmetry axis of the cyclodextrin, the unit vector nz points towards 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 same number of water molecules, either in the absence or in the presence of one cyclodextrin molecule.

14


Page 15 of 47 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


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 did also not affect the outcome of the simulations because the 4 C1 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 4 C1 conformation. Also the alternative cut-off set in case of CHARMM36 and using GROMACS instead of GROMOS in case of the 53A6GLYC 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 , 2016H66, the CHARMM36 force field as well as q4md-CD show a very similar occurence of conformation types. Throughout a vast majority of the simulation time the glucopyranose units adopt the chair conformation 4 C1 . In the simulation with GROMOS force field 53A6, however, the dominant conformation type is the inverted chair conformation 1 C4 . With the force field 56A6CARBO , 4 C1 as well as 1 C4 conformation occur with alternating population between neighbouring 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

15



53A6

56A6CARBO

53A6GLYC

56A6CARBO_R

2016H66

CHARMM36

q4md−CD

1.0 0.8 glucosidic unit 1

0.6 0.4 0.2 1.0 0.8

glucosidic unit 2

0.6 0.4 0.2 1.0 0.8

glucosidic unit 3

0.6 0.4 fraction

0.2 1.0 0.8 glucosidic unit 4

0.6 0.4 0.2 1.0 0.8

glucosidic unit 5

0.6 0.4 0.2 1.0 0.8

glucosidic unit 6

0.6 0.4

4C 1 1C B3 4 ,O B2 ,5 B1 ,4 3,O B 2,5 B 1,4 B 1S 3 1S 5 3S 5 3S 1 5S 1 un as 5S3 sig ne d






0.2 4C 1 1C B3 4 ,O B2 ,5 B1 ,4 3,O B 2,5 B 1,4 B 1S 3 1S 5 3S 5 3S 1 5S 1 un as 5S3 sig ne d

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 47

conformation type

Figure 3: Occurence of glucopyranose ring conformation types 154,155 in α-cyclodextrin simulations. (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 Figs. 4 and S3a, respectively) while for γ-CD the flipped state remains stable over several nanoseconds (see Fig. 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 Unfortunely, we do

16


Page 17 of 47 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


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 experiment 163–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 4 C1 (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). As for β-cyclodextrin and γ-cyclodextrin the simulations show more glucose unit rotations 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 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 a 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 for 56A6CARBO R and 2016H66, see Figure 5. For β-cyclodextrin 40% of the frames with force field 56A6CARBO R

17



53A6

56A6CARBO

53A6GLYC

56A6CARBO_R

2016H66

CHARMM36

q4md−CD

180

ξ [°]

90 0 −90 −180 0

20 40 60 80 100 0 time [ns]

20 40 60 80 100 0 time [ns]

20 40 60 80 100 0 time [ns]

20 40 60 80 100 0 time [ns]

20 40 60 80 100 0 100 200 300 400 500 0 100 200 300 400 500 time [ns]

time [ns]

time [ns]

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

56A6CARBO

53A6GLYC

56A6CARBO_R

2016H66

CHARMM36

q4md−CD

0.05 0.04 Probability

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 47

0.03 0.02 0.01 0 −90 0 90 180 flip dihedral angle / °

−90 0 90 180 flip dihedral angle / °






Figure 5: Distribution of the flip dihedral angle ξ for α-cyclodextrin. The different colors denote the different glucopyranose units. have at least one glucose unit rotated (see Fig. S3a). For γ-cyclodextrin more frames have rotated units and force field 2016H66 also shows some rotations (see Fig. 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 an experimental evidence of tumbling glucopyranose units would be desirable due to the impact on inclusion mechanisms. 35,163–165 However, below we will discuss that some NMR derived data leads to the conclusion that at least for α-cyclodextrin the frequency of glucopyranose unit tumblings seems to be more realistic for CHARMM36 and q4md-CD compared to the 18


Page 19 of 47 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


GROMOS force fields.

Structural Properties Table 2 shows averaged structural properties obtained with the different force fields compared to experimental data for crystalline α-cyclodextrin reported in the literature. For the dihedral angles along the glycosidic linkage φ (O5 C1 O1 C’4 ) and ψ (C1 O1 C’4 C’3 ) the results with CHARMM36 are very close to the experimental data. The q4md-CD results show also good agreement. With force fields 53A6GLYC , 56A6CARBO R and 2016H66 the values for φ and ψ differ less from each other and lie above the experimental data for φ and below for ψ. For force fields 53A6 and 56A6CARBO the results are far from the experimental values. For the dihedral angle ζ (O3 C3 C2 O2 ), which describes the rotation of the secondary hydroxyl groups of one glucose unit against each other, the values obtained with force fields 53A6 and 56A6CARBO are significantly larger than in experiment. Likewise the distances between the secondary hydroxyl groups of neighbouring glucose units dO23 obtained with these force field versions is significantly larger than the experimental value. A comparison between 56A6CARBO R and 2016H66 reveals that the different Lennard-Jones parameter used for OA and OE have an impact on the dihedral angle ζ, which is significantly smaller in case of 2016H66. The puckering parameters show which conformation type prevails for the glucopyranoside units and how far the ring atoms are deflected from the ring plane. For force fields 53A6GLYC , 56A6CARBO R and 2016H66 as well as CHARMM36 and q4md-CD the small azimuthal angle θ indicates the chair conformation 4 C1 . We note that 56A6CARBO R and CHARMM36 predict similar ring inversion properties for unfunctionalized monomers while 53A6GLYC shows significantly higher ring-inversion free energies for α-d-glucose in water. 166,167 When inserting a hexopyranose residue into a chain the force fields differ qualitatively and quantitatively in their predictions of the effect on ring flexibility, 166,167 which is also reflected in the results obtained in the present work, in particular for β-cyclodextrin, which exhibits a larger azimuthal angle (Table S1) and ring flexibility (Figure S1), respec19



Page 20 of 47

Table 2: Structural Properties of α-Cyclodextrin obtained from MD Simulations using Different Force Fields and from Experimenta . Property φ dihedral angle [◦ ] ψ dihedral angle [◦ ] ζ dihedral angle [◦ ] δ angle [◦ ] dO1 [nm] dO23 [nm] dO6 [nm] circularityb Puckering amplitude Q [nm] θ azimuthal angle [◦ ] a

53A6 74.21 160.68 148.78 119.41 0.482 0.620 0.507 0.875 0.052 164.03

56A6CARBO 80.00 167.35 116.63 118.68 0.451 0.729 0.576 0.766 0.055 110.01

53A6GLYC 115.27 118.55 79.74 119.68 0.434 0.287 0.566 0.925 0.056 8.86

56A6CARBO R 117.17 115.99 80.46 119.60 0.421 0.301 0.551 0.894 0.068 14.21

2016H66 115.57 118.17 52.39 119.64 0.424 0.304 0.551 0.916 0.058 12.60

CHARMM36 107.53 129.00 56.18 119.42 0.426 0.321 0.531 0.901 0.057 10.27

q4md-CD 111.73 126.27 69.62 119.38 0.417 0.303 0.527 0.911 0.055 12.08

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

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. b ratio of the smallest to the largest distance between any pair of glucose O1 atoms that lie across the ring from each other. c Ref. 94. d Ref. 161. e Ref. 168. tively. The average θ angles for force fields 53A6 and 56A6CARBO are very high due to the prevalent occurrence of the inverted chair conformation 1 C4 . The experimental value

for θ is very small which indicates a dominant chair conformation. The corresponding tables for β- and γ-cyclodextrin are presented in the Supporting Information as Tables S1 and S2, respectively. For β-cyclodextrin all force fields except 53A6 and 56A6CARBO show very similar average φ-angles while the average ψ-angle obtained with 2016H66 is smaller than for 53A6GLYC and 56A6CARBO R . For γ-cyclodextrin the three GROMOS force fields 53A6GLYC , 56A6CARBO R and 2016H66 show very similar φ and ψ angles that are close to the experimental values. Compared to these GROMOS force fields, CHARMM36 shows better agreement for ψ and worse agreement for φ. Note that the CHARMM36 values include the effects of frequently occurring tumbling events (see above). However, q4md-CD shows the best agreement to the experimental values for these angles in case of γ-cyclodextrin. Table S3 shows averages of the same structural properties as seen in Table 2 obtained from MD simulations using the 53A6GLYC force field with different schemes to treat the electrostatic interactions in GROMOS as well as in GROMACS. For the structural properties considered here there is no significant difference observed showing that these properties do not depend on subtle details of how the nonbonded interactions are evaluated. 169 In summary, the comparison to structural properties showed that 53A6 and 56A6CARBO 20


Page 21 of 47 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


are not suitable to model the structural features of native cyclodextrins. The glycosyl dihedral angles φ and ψ are best described by CHARMM36. Among the tested force fields, q4md-CD is the only one especially optimized for cyclodextrins. However, the only structural property were it outperforms the other tested force fields is the ζ dihedral angle, which describes the orientation of the secondary hydroxyl groups within one glycosidic unit (see Figure 1). It has to be considered, however, that the experimental data from X-ray analysis represents the crystalline structure of α-cyclodextrin which may differ significantly from the solution structure. 170–172 Cyclodextrin in aqueous solution is expected to exhibit a more flexible character such that comparison to NMR data will be more informative and thus will be discussed in the next section.

Comparison to NMR Experiments Table 3 shows NMR coupling constants, proton-proton distances and relative populations of the three staggered conformers of the hydroxymethyl group obtained with different force fields in comparison to experiment. The NMR coupling constants along the glycosidic linkage, 3 JC4,H1′ and 3 JC1′ ,H4 , show good agreement with the experimental values for force fields 53A6GLYC , 56A6CARBO R and 2016H66 as well as CHARMM36 and q4md-CD. Figure 6 shows a contour map of the respective dihedral angles φ’ and ψ’. The φ’-ψ’ distribution is broader and shifted to lower φ’-values for force fields 53A6 and 56A6CARBO . The CHARMM36 force field exhibits a broader distribution compared to the force fields 53A6GLYC , 56A6CARBO R and 2016H66, indicating more structural flexibility. The distribution obtained with q4md-CD indicates more flexibility compared to the GROMOS force fields but less flexibility compared to CHARMM36. The overall shape of the contour map for the CHARMM36 force field is in agreement with the one reported in ref. 173 derived from NMR data, including the two peaks in the distribution of ψ’-angles. In contrast, the force fields 53A6GLYC , 56A6CARBO R and 2016H66 show only one peak in the φ’-ψ’ distribution while q4md-CD exhibits a small second peak. Note that the latter was not reported in the original work, possibly due to 21



Page 22 of 47

Table 3: NMR Coupling Constants (3 J), 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 JC,H c [Hz] 3

JOH,CH c [Hz]

rNOE [nm]g ωtor.pop. [%]

H1’-C4 H4-C1’ O(2)H O(3)H O(6)H H1-H2 H1’-H4 gg gt tg

53A6 3.28 3.25 4.38 3.92 4.84 0.24 0.22 14 85 1

56A6CARBO 3.35 4.58 4.52 5.04 5.42 0.23 0.21 27 41 32

53A6GLYC 5.38 5.41 5.17 2.64 4.92 0.24 0.21 62 38 0

56A6CARBO R 5.31 5.39 4.59 3.23 5.57 0.24 0.19 53 45 2

2016H66 5.34 5.39 3.67 3.12 5.57 0.24 0.19 57 42 1

CHARMM36 4.99 4.88 7.96 5.40 6.55 0.25 0.22 55 44 1

q4md-CD Experiment 5.16 5.2d ,6.4e 5.02 5.0d 6.44 6.6f 3.32 < 3f 5.91 5.6f 0.24 0.2512d 0.21 0.2245d 47 65b 52 34b 1 1b

a

The simulations using the various GROMOS force fields were conducted with the standard settings, see main text. b Ref. 12. c In the GROMOS force field non-polar hydrogens are not represented explicitly. Therefore a phase shift of −120◦ has been used to calculate φ′ from φ and ψ ′ from ψ, respectively. In case of O(2)H,CH and O(3)H,CH-couplings phase shifts of +120◦ have been used to calculate the required dihedral angles from the C3n C2n O2n HO2n and C2n C3n O3n HO3n dihedral angles, respectively. In case of O(6)H,CH coupling two calculations with phase shifts of −120◦ and +120◦ have been performed and subsequently averaged because the C6-atom carries two virtual hydrogen atoms. d Ref. 173. e Ref. 174. f Ref. 175. g For simulations with the GROMOS force field the positions of the hydrogen atoms were calculated based on a standard geometry. 115 shorter trajectories (see Figure S2 in ref. 12). A recent simulation study 6 of native cyclodextrins employing the GLYCAM-06j force field 55 also showed a multimodal character of the ψ-angle distribution, with two well separated peaks for β- and γ-cyclodextrin and thus confirm the interpretation of the NMR experiments conducted by Thaning et al., 173 who concluded that a bimodal distribution fits best to their measurements. If the φ’-ψ’ distributions for CHARMM36 and q4md-CD were recalculated by excluding frames that contain a rotated glucose unit, compared to the corresponding distributions in Figure 6, both the “island” at ca. -60◦ /150◦ and the “peninsula” at ca. -60◦ /-50◦ vanish. That is, the distributions would show only one peak. Considering that two peaks were derived from NMR data 173 and here a second peak occurs due to tumbling of glucose units, we could conclude that at least for α-cyclodextrin the frequency of unit tumblings (see above) are better described by CHARMM36 and q4md-CD than the GROMOS force fields. However, here we find that the NMR parameters calculated with force fields 53A6GLYC , 56A6CARBO R 22


Page 23 of 47 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


56A6CARBO

53A6

53A6GLYC

56A6CARBO_R

2016H66

CHARMM36

q4md−CD

240

0.01

180

0.008

120

0.006

ψ´ / ° 60

0.004

0

0.002

−60 −120

0 −60

0

60

120 −60

0

60

120 −60

0

60

120 −60

0

60

120 −60

0

60

120 −60

0

60

120 −60

0

60

120

φ´ / °

Figure 6: Contour maps of ψ ′ and φ′ -angles for α-cyclodextrin. The dihedral angles of all six glucose units contributed to the distributions. and 2016H66 agree well with experiment despite being based on unimodal φ-ψ-distributions. We note, however, that Thaning et al. 173 also considered residual dipolar couplings in their analysis that were not calculated in the present work. The NMR coupling constants for the hydroxyl groups, 3 JO(2)H,CH , 3 JO(3)H,CH and 3 JO(6)H,CH exhibit significant differences between the GROMOS force fields and CHARMM36, the latter showing higher values. Indeed the distribution of the corresponding dihedral angles show substantial differences between CHARMM36 and the GROMOS force fields 53A6GLYC , 56A6CARBO R and 2016H66 which show similar distributions with peaks at the same positions, see Fig. 7. CHARMM36 shows significantly different distributions for each of the dihedral angles, while those obtained with q4md-CD are more similar to the CHARMM36 distributions compared to the GROMOS ones. The proton-proton distances rNOE show good agreement with the experimental values for all force fields applied. Comparing the GROMOS force fields, the experimental values for the relative populations of the three staggered conformers of the hydroxymethyl group are reproduced significantly better with force fields 53A6GLYC , 56A6CARBO R and 2016H66 than with force fields 53A6 and 56A6CARBO . Among the tested force fields only 53A6GLYC results in populations in good agreement with the experimental data. For 56A6CARBO R , 2016H66, CHARMM36, and q4md-CD the gt population is too high. For β- and γ-CD the three GROMOS force fields 53A6GLYC , 56A6CARBO R and 2016H66 slightly underestimate 3

JO(2)H,CH , while CHARMM36 overestimates it. For 3 JO(3)H,CH and 3 JO(6)H,CH CHARMM36

predicts too large values while the GROMOS values are in good agreement (see Tables S4

23



56A6CARBO

53A6

53A6GLYC

56A6CARBO_R

2016H66

CHARMM36

q4md−CD

HO2CH

0.02 0.01 0.00

HO3CH

0.02 0.01 0.00

HO6CH (1)

Probability

0.02 0.01 0.00

HO6CH (2)

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 47

0.02 0.01 0.00 −90

0

90

180

−90

0

90

180

−90

0

90

180

−90 0 90 dihedral angle / °

180

−90

0

90

180

−90

0

90

180

−90

0

90

180

Figure 7: Distributions of the dihedral angles Hi − Ci − Oi − Hi with i = 2,3,6 for αcyclodextrin. The different colors denote the different glucopyranose units. and S5). Table S6 shows averages of the same parameters as shown in Table 3 obtained with the 53A6GLYC force field with different schemes to treat electrostatic interactions in the GROMOS code and the GROMACS code, respectively. Again the properties are insensitive to the treatment of electrostatic interactions although some minor effect on the rotamer population is visible. The comparison to NMR experiments shows again that among the GROMOS force fields 53A6GLYC , 56A6CARBO R and 2016H66 are to be preferred compared to 53A6 and 56A6CARBO . CHARMM36, q4md-CD and the three recommended GROMOS versions show similar performance for the NMR data, except that CHARMM36 slightly overestimates the 3 J-coupling constants related to the primary and secondary hydroxyl groups.

Intramolecular Hydrogen Bonds Intramolecular hydrogen bonds are being formed between all O2 and O3 hydroxy groups of adjacent glucose units. 176 Table 4 shows the intramolecular O-H...O hydrogen bonds oc24


Page 25 of 47 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


curring in the simulation trajectories. For force fields 53A6 and 56A6CARBO hardly any intramolecular hydrogen bonds are formed. Force fields 53A6GLYC , 56A6CARBO R and 2016H66 show high population of about 60% of hydrogen bonds for donors of type O3 and less populated bonds for donors of type O2 with decreasing population in the order 53A6GLYC , 56A6CARBO R and 2016H66. The simulation time averages of the H...A distances are distinctly larger for donors of type O2 and the simulation time averages of the D-H...A angles are smaller. With force fields CHARMM36 and q4md-CD both hydrogen bond types are equally populated. Tables S7 and S8 show the intramolecular O-H...O hydrogen bonds occurring in β- and γ-cyclodextrin. Again, the force fields 53A6 and 56CARBO exhibit hardly any hydrogen bonds. As for α-cyclodextrin, force fields 53A6GLYC , 56A6CARBO R and 2016H66 show higher population for hydrogen bonds involving the O3-donor compared to those involving the O2-donor, in contrast to CHARMM36 and q4md-CD which show slightly higher population for the latter. Moreover, the decreasing population of O2-donor hydrogen bonds in the order 53A6GLYC , 56A6CARBO R and 2016H66 is also present in β- and γ-cyclodextrin. For β-cyclodextrin a comparison between 56A6CARBO R and 2016H66 shows that the O3donor hydrogen bonds are much higher populated for the latter force field (approx. 60% vs. 42%). For γ-cyclodextrin this difference is not as pronounced. For the three-centered hydrogen bonds, shown in Table 5, the 53A6GLYC force field shows the highest overall populations with the hydrogen bonds involving a O3 donor having a more than three times higher population than those involving a O2 donor. The higher population of O3-type donor hydrogen bonds is also seen for the force fields 56A6CARBO R and 2016H66 while CHARMM36 predicts the same occurrence for both types three-centered hydrogen bonds. In contrast, q4md-CD predicts a higher occurence of O2-type donor hydrogen bonds. The simulations using the force fields 53A6 and 56A6CARBO show hardly any three-centered hydrogen bonds. Tables S9 and S10 show the three-centered hydrogen bonds observed in βand γ-cyclodextrin, respectively. For β- and γ-cyclodextrin the 56A6CARBO R and 2016H66 force fields show slightly less populated O3-donor hydrogen bonds compared to 53A6GLYC .

25



Page 26 of 47

Table 4: Intramolecular Two-centered Hydrogen Bonds in α-Cyclodextrina . Donor-H-Acceptor

53A6

56A6CARBO

O21 − HO21 − O32 O22 − HO22 − O33 O23 − HO23 − O34 O24 − HO24 − O35 O25 − HO25 − O36 O26 − HO26 − O31 O32 − HO32 − O21 O33 − HO33 − O22 O34 − HO34 − O23 O35 − HO35 − O24 O36 − HO36 − O25 O31 − HO31 − O26

0.10 0.15 0.28 0.07 0.16 0.01 1.18 1.18 1.16 0.35 0.59 1.10

0.55 1.21 0.36 0.79 0.02 2.50 3.23 7.30 0.77 2.21 0.32 3.19


0.5879 0.5842 0.5765 0.5898 0.5888 0.5861 0.6569 0.6540 0.6415 0.6536 0.6617 0.6566

0.5435 0.4848 0.5949 0.5493 0.5959 0.5484 0.5494 0.4774 0.5660 0.5913 0.6232 0.5871


107.85 108.70 106.33 104.19 109.47 108.38 65.84 66.11 66.62 65.42 64.92 65.41

85.25 78.09 59.71 99.25 89.14 97.18 84.16 84.93 80.65 74.68 73.95 74.12

53A6GLYC 56A6CARBO R Population in % 25.88 14.59 25.22 16.05 25.52 15.71 26.38 15.55 24.28 15.95 24.92 16.04 62.85 60.35 63.38 60.18 63.54 61.37 62.16 56.30 64.29 58.88 63.90 60.67 Distance H...A in nm 0.3044 0.3465 0.3070 0.3331 0.3063 0.3321 0.3049 0.3513 0.3075 0.3336 0.3064 0.3334 0.2405 0.2581 0.2387 0.2480 0.2383 0.2433 0.2407 0.2674 0.2376 0.2519 0.2379 0.2483 Angle D-H-A in degree 77.15 61.64 75.44 64.40 75.60 62.55 77.04 62.49 74.92 65.31 75.37 63.35 123.60 124.78 125.27 125.02 125.11 126.06 123.72 122.14 125.83 123.60 125.27 124.17

2016H66

CHARMM36

q4md-CD

9.50 9.13 9.63 9.34 9.25 9.44 65.50 66.16 64.90 66.66 65.58 64.66

29.69 28.34 29.80 27.95 28.96 28.69 31.01 32.69 30.21 32.24 30.90 31.76

41.22 39.96 40.13 40.10 39.77 38.69 37.68 36.50 37.29 37.62 36.42 37.21

0.3523 0.3519 0.3533 0.3492 0.3522 0.3550 0.2423 0.2399 0.2444 0.2372 0.2412 0.2459

0.3114 0.3129 0.3128 0.3153 0.3152 0.3151 0.2989 0.2964 0.3018 0.2981 0.3019 0.2994

0.2786 0.2801 0.2808 0.2821 0.2825 0.2825 0.2952 0.2990 0.2975 0.2960 0.2975 0.2963

54.63 54.07 55.00 54.45 54.14 54.50 133.44 134.31 132.93 134.61 133.57 132.58

93.75 91.95 94.30 91.35 93.01 92.12 102.15 103.46 101.50 103.33 101.95 102.89

104.71 103.31 103.11 103.64 102.83 102.27 95.14 93.86 94.52 95.13 93.37 94.24

a

A hydrogen bond is present when the distance H–O is smaller than 0.25 nm and the angle Donor–H–Acceptor is larger than 135◦ . Like for α-cyclodextrin, CHARMM36 exhibits little difference in population between O2and O3-donor type hydrogen bonds, while q4md-CD shows significant higher population of three-centered hydrogen bonds involving an O2-donor. Tables S11 and S12 show the influence of the treatment of electrostatic interactions on the population of two-centered and three-

26


Page 27 of 47 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 5: Intramolecular Three-centered Hydrogen Bonds in α-Cyclodextrina . Donor-H-A1/A2

53A6

56A6CARBO

53A6GLYC

O21 − HO21 − O32 /O11 O22 − HO22 − O33 /O12 O23 − HO23 − O34 /O13 O24 − HO24 − O35 /O14 O25 − HO25 − O36 /O15 O26 − HO26 − O31 /O16 O32 − HO32 − O21 /O11 O33 − HO33 − O22 /O12 O34 − HO34 − O23 /O13 O35 − HO35 − O24 /O14 O36 − HO36 − O25 /O15 O31 − HO31 − O26 /O16

0.07 0.06 0.17 0.04 0.07 0.05 0.83 0.91 0.89 0.27 0.46 0.89

0.28 0.60 0.21 0.25 0.01 1.50 2.20 4.85 0.59 1.05 0.24 2.39

7.98 7.38 7.76 8.05 8.03 7.64 28.04 27.70 28.21 26.92 28.33 28.00

56A6CARBO R 2016H66 Population [%] 3.10 1.72 3.86 1.89 3.54 1.89 3.30 1.65 3.89 1.89 3.45 1.88 22.69 23.78 23.20 24.17 22.74 24.45 20.65 24.94 22.93 24.74 23.01 24.68

CHARMM36

q4md-CD

17.47 16.97 17.63 16.51 17.21 17.04 17.80 18.49 17.43 18.12 17.57 18.00

34.03 32.69 32.77 32.69 32.41 31.27 24.26 23.22 23.90 23.99 23.62 24.24

a

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 centered hydrogen bonds calculated for α-cyclodextrin with the 53A6GLYC force field. No significant differences are observed between different methods to handle these interactions. Compared to the other properties discussed so far, the tested force fields show significant differences in the intramolecular hydrogen bonds. Unfortunately, due to a lack of experimental data, it is not possible to judge which observed hydrogen bonds are the most realistic ones. However, we note that differences between the GROMOS force fields on the one hand and CHARMM36/q4md-CD on the other may not only be caused by the intramolecular interactions but also by the different water models. 177

Hydration Pattern Figure 8 shows cylindrical distribution functions describing the water density around αcyclodextrin calculated with all investigated force fields. The normalized density is sampled with respect to their cylinder coordinates r and z. The origin (r = z = 0) is chosen as the center of mass of the O1 atoms. The z-coordinate points along the central cyclodextrin axis with positive values in the direction of the primary rim. The dark blue regions in

27



Figure 8 reveal an enhanced water molecule density inside the cyclodextrin cavity as well as in a layer around the cyclodextrin molecule. The dark-red and red colored areas indicate the zones where the probability encounter a water molecule is lowest. This is the case where the water is replaced by the atoms of the cyclodextrin molecule and also at the openings of the cyclodextrin cavity. For force fields 53A6GLYC , 56A6CARBO R and 2016H66 the cyclindrical distribution functions of water are very similar. The hydrophobic areas are slightly more pronounced for 53A6GLYC . For 53A6 and 56A6CARBO the distorted structure of the cyclodextrin molecule is reflected in the distribution functions. The structures modeled with CHARMM36 and q4md-CD appear slightly less bulky which is also reflected in their partial molar volume (see below). Figure 9 shows the numbers of cumulated water molecules as well as the radial distribution functions (RDFs) of water molecules around the center of mass of the O1-atoms of α-cyclodextrin for all investigated force fields. The force fields 53A6GLYC , 56A6CARBO R , 2016H66, CHARMM36 and q4md-CD show three peaks, the first sharp one corresponds to water molecules located inside the cavity while the two other peaks at r ≈ 0.85 and r ≈ 1.1 nm indicate the first and second hydration shell around the α-cyclodextrin molecule. For the force fields 53A6 and 56A6CARBO for which the cyclodextrin has a more distorted shape the first peak is displaced from the origin and the second peak is broader and shifted to the right. For β- and γ-CD the position of the first peak is displaced from the origin by about 2 ˚ A for the force fields 53A6GLYC , 56A6CARBO R , 2016H66, CHARMM36 and q4md-CD (see Figures S7 and S8) as reported before. 12 The number of cumulated water molecules are very similar among the force fields, except for 53A6, which shows a slightly larger number of water molecules for r

Validation and Comparison of Force Fields for Native Cyclodextrins in

Recommend Documents