2.1 Simulation details
All simulations and primary data analysis or data extractions were performed with the CHARMM5
program (version c35b3), using the recent C35 carbohydrate FF.1,2
Starting with the c35b1 general release, the program includes a feature to save the pressure tensor at every integration step, essential for the accurate calculation of viscosity. All systems were simulated with a 1 fs time step and the Particle Mesh Ewald (PME) method6
for electrostatics, and a Lennard-Jones (LJ) switching function for the van der Waals interactions;7
the PME real space and LJ cutoff was 12 Å, the parameter set default. The PME calculations used a 5th
order polynomial spline lookup, κ = 0.32 Å−1
, and grid spacing ≈ 1 Å. Additional water and glucose simulations employed the 3D-IPS/DFFT (Three-Dimensional Isotropic Periodic Sum with the Discrete Fast Fourier Transform extension) non-bond method8
for both electrostatics and LJ interactions, with a real space cutoff of 10 Å. This method, henceforth denoted IPS, obtains both long-range electrostatic and LJ interactions; it is an alternative to PME with a long-range LJ correction. Cubic periodic boundary conditions were used for all simulations, in conjunction with the extended system methods9,10
for the NVT and NPT ensembles (N is particle number, V is volume, T is temperature, and P is pressure).
A box of 1340 TIP3P11,12
water molecules, with a cube edge length of ca. 34 Å, was simulated in triplicate for 12.5 ns at T = 20, 25, 30, 40 and 50 °C, using both PME and IPS non-bond methods (30 simulations total) in the NVT ensemble. By changing the random seed, each of the 3 replicate simulations was assigned different initial velocities. The edge length was set to match the experimental density for water at each T. The Hoover thermostat with a coupling constant of 5000 kcal mol−1
was used to regulate the temperature. The last 12 ns were used for data analysis, with standard errors computed over the replicate simulations.
Glucose (D-glucopyranose) solutions were prepared by positioning molecules on a grid based on a perfect cube (27, 64, or 125), with grid locations skipped at random, based on the number of molecules. One third of the molecules were randomly assigned an α anomeric conformation, and the rest were β, resulting in a 2:1 β:α anomeric ratio. In addition, each molecule was given a random rotation, and a small random off-lattice translation. Solutions of ca. 0.3, 0.5, 1, 2, 3, 4 and 5 molal (m) were prepared by first placing 7, 12, 24, 48, 72, 96 and 120 glucose molecules, then combining with a cube of TIP3P water large enough to encompass all of the molecules. Waters overlapping the glucose molecules were then deleted, the water box trimmed if necessary, and finally random water molecules were deleted to obtain a final count of 1320 water molecules. The coordinates and the cubic periodic boundary lattice parameters were then each briefly optimized via minimization. Each solution was then equilibrated via 500 ps of NPT dynamics at 20 °C (30 °C for 0.3 m) using the IPS method, with a pressure piston mass of 500 amu. The average volume was calculated over the last 100 ps, and the stored trajectory frame with a unit cell edge that best matched the average volume was selected for the NVT simulations. Two separate simulations, one using the PME method and one using IPS, were performed at each concentration for at least 16 ns; in order to improve convergence, those at 0.3 m and the PME simulations at 3 m and above were extended to 20 ns. The simulations at 20 °C were used to compute viscosity at the experimental temperature; for translational diffusion, PME simulations were performed via the same protocol at 25 °C for the 0.5 m and higher concentrations. Coordinate sets were saved every ps for all but the 0.3 m simulations; these were used for the rotational relaxation calculations, and coordinate sets were saved every 0.2 ps. Data analysis intervals were 15 and 18 ns, with blocks of 5 or 6 ns used to estimate standard errors.
Three melibiose (D-galactopyranosyl-α-1,6-D-glucopyranoside) solutions of 0.15, 0.20, and 0.25 m were prepared by combining 3, 4, or 5 disaccharides with 1100 waters. The solute molecules were placed on a cubic grid with grid points centered approximately 1/3 the size of the side of the cubic box then TIP3P waters were added to the box. The melibiose molecules were generated with a 2:1 β:α anomeric ratio. Following extensive NPT equilibration, the average volume at each concentration was used to run 37 ns NVT simulations with the PME method, saving coordinate sets every ps. The final 36 ns were used for analysis, with standard errors estimated from four 9 ns blocks.
Five trehalose (D-glucopyranosyl-α,α-1,1-D-glucopyranoside) solutions of ca. 0.10, 0.30, 0.56, 1.11, and 1.26 m were prepared by combining 2, 6, 11, 22, or 25 disaccharides with 1100 waters, analogous to the procedure used for melibiose. A total of seven simulations were run, with three at 25 °C at 0.10, 0.30, and 1.26 m, and four at 30 °C at 0.10, 0.30, 0.56, and 1.11 m. As for melibiose, the average NPT volume at each concentration was used for subsequent NVT simulations, and coordinate sets were saved every ps. Concentrations below 1.0 m were simulated for 25 ns each, while the two high concentrations were simulated for 37 ns each. The first ns of each simulation was used as equilibration, and four blocks of either 6 or 9 ns were used to estimate standard errors.
Solutions of ca. 0.12 m (4 wt%) maltose (D-glucopyranosyl-α-1,4-D-glucopyranoside) and solutions of trehalose at the same concentration were prepared by placing 8 molecules uniformly spaced about the origin, and combining with 3648 TIP3P waters. For each disaccharide, five replicate solutions were made, with random rotations and small random displacements used to generate five unique starting conditions. Each was simulated for 1 ns under NPT conditions; the volume for the NVT simulations was computed as the average over the last 500 ps from all five replicates for each disaccharide. A stored coordinate set which best matched the chosen volume was chosen from the last 500 ps of each replicate, and each was then simulated for 12.5 ns under NVT conditions. For the maltose replicates, 3 of the 8 molecules were assigned α anomeric conformations, and the other 5 were assigned β conformations. Coordinate sets were saved every 0.25 ps for these simulations to improve the accuracy of the short time behavior of the correlation functions used to compute rotational relaxation. The last 12 ns of each simulation were used for analysis, and standard error estimates were calculated over the five replicates. These simulations were also used for hydration shell analysis for these two disaccharides, which counted water O atoms within a shell thickness from any sugar atom ().
Number of watersa as a function of shell thickness for maltose and trehalose.
The lowest concentration simulations for glucose, maltose, melibiose and trehalose were used to extract heavy atom coordinate sets of a sugar with nearby waters, for a series of calculations of the diffusion constant using a hydrodynamic bead model (described in the next subsection). The nearby waters were defined on the basis of an energy cutoff for the interaction of a water molecule with the sugar molecule; the cutoff was specified as a multiple of kBT, where kB is the Boltzmann constant. For each extraction, the sugar molecule was centered at the origin, and the sugar-water interaction energy computed. The sugar molecule, along with water molecules with an interaction energy more favorable than the cutoff, were written to a formatted coordinate file. Energy cutoffs of -2, -3, -4, and -5 kBT were initially evaluated for each of the four sugars, with a second series of four sugar dependent cutoffs at intervals of 0.1 kBT in the range between -3 and -4 kBT to obtain more exact matches to experimental values. Different sugar molecules were sampled from each simulation trajectory, at intervals of 250 or 500 ps between time points; ca. 200 coordinates sets were extracted for each disaccharide at each multiple of kBT, and ca. 380 coordinates sets were extracted for glucose. The same sugar coordinates were used for the calculation of D and the solvent accessible surface area (SASA) with no hydration.
2.2. Calculation of viscosities, diffusion constants and NMR T1
Viscosities were obtained from
are the off-diagonal elements of the instantaneous pressure tensor. The integrands typically converged between 2 and 5 ps at low concentrations (< 10 wt%), and 10–20 ps at higher concentrations (10–50 wt%). Visual inspection of the integrals was required to confirm the location of the stable plateau region; the viscosity estimate is an average over a window of width 0.5 to 1.5 ps, depending on the width of the stable plateau. The pressure tensor was assumed to be symmetric (Pαβ
); the CHARMM program was used to compute the correlation function of the 3 unique off-diagonal elements, average them, and then compute the integral.
Diffusion constants evaluated from a simulation in periodic boundary conditions, DPBC
, must be corrected for periodic boundary condition artifacts (a correction for viscosity is not required). The formula of Yeh and Hummer 13
is highly accurate:
is the cubic box length, and ξ = 2.837297, and η is the viscosity of the solution calculated from the simulation. The continuous (unfolded) center of mass time series for each sugar was used to compute the mean squared displacement <r2
>, and DPBC
was obtained from a linear fit of <r2
> vs. t
based on the Einstein relation for diffusion. A FORTRAN program was used to calculate D
from the <r2
> vs. t
data, computing the slope and finite size correction.
Diffusion constants were also computed using a hydrodynamic bead model methodology developed by Garcia de la Torre and Bloomfield,14
and later tailored to included bound water.15–18
For the present application, single sugar molecules with varying amounts of bound water (including none) were prepared as described in the previous subsection. <D
> and <Nb
> were computed over the results at a specific energy cutoff for bound water, where Nb
is the number of bound waters in each coordinate set. Translational and rotational friction tensors were evaluated with the FORTRAN program denoted here HI6 (for Hydrodynamic Interaction, version 6). The original implementation is described in ref. 15
. Subsequent improvements include automated recognition of coordinate file type (original free format, PDB, or CHARMM), calculation of vector relaxation in the diffusive frame, and automatic array allocation based on input file size. HI6 is available for download at URL http://www.lobos.nih.gov/mbs/scriptprog.shtml
Assuming that relaxation is due to dipolar interactions between the 13
C nucleus and its N
attached protons, the 13
C NMR relaxation rates 1/NT1
for each carbon are given by19
is Plank’s constant divided by 2π, rC-H
is the effective C-H bond length, γH
, and ωC
are the gyromagnetic ratios and Larmor frequencies, respectively, of the 13
C and 1
H nuclei; ωC
, where H
is the field strength. J(ω)
is the spectral density of the second rank reorientational correlation function:
is the second order Legendre polynomial
) and is the unit vector along the C-H bond direction at time t
. If tumbling is rapid (all decays of C(t)
are much less than 1/ωH
), then Eq. (3)
where τ is the integral of C(t)
. Setting rC-H
=1.117 Å 20
yields the numerical factor in Eq. (5)
. Time series of individual C-H bond vectors were extracted from the CHARMM trajectories, normalized to the unit vector, and the time correlation function computed; the C(t)
were averaged over molecules (and replicate simulations when available), and then integrated numerically to obtain the relaxation time.