|Home | About | Journals | Submit | Contact Us | Français|
Proper simulation of dynamic properties, including molecular diffusion, is an important goal of empirical force fields. However, the widely used TIP3P water model does not reproduce the experimental viscosity of water. Consequently, scaling of simulated diffusion constants of solutes in aqueous solutions is required to effectively compare them with experiment. It is proposed that scaling by the ratio of viscosities of model and real water is appropriate in the regime where the concentration dependence of simulated and experimental solution viscosities is parallel. With this ansatz, viscosity scaling can be carried out for glucose and trehalose up to 20 wt% for simulations carried out with the CHARMM additive carbohydrate force field C35 and TIP3P water; above this value, the concentration dependence of simulated viscosities lags that of experiment, and scaling is not advised. Scaled translational diffusion constants for glucose, and the disaccharides trehalose, maltose, and melibiose at low concentration agree nearly quantitatively with experiment, as do NMR 13C T1’s for glucose, trehalose, and maltose; these results support the use of C35 for simulations of sugar transport properties at low concentration. At high concentrations the scaled diffusion constants for glucose and trehalose underestimate and overestimate experiment, respectively. Hydrodynamic bead model calculations indicate a hydration level of approximately 1 water/hydroxyl for glucose. Patterns for the disaccharides are more complicated, though trehalose binds 0.5 to 1 more water than does maltose depending on the analysis.
This paper concerns the comparison of diffusion constants obtained from simulation based on empirical force fields and experiment. While many considerations are general, the focus is on validating the recently developed CHARMM additive carbohydrate1,2 force field (FF), C35, using the translational diffusion constant, D, and NMR spin lattice relaxation time, T1, for mono and disaccharides. Both of these properties are sensitive to solute-solvent interactions in that bound water increases the effective hydrodynamic size of the solute, and therefore the translational and rotational friction. T1’s also contain contributions from internal motions such as libration in torsional minima and isomerization, and can thereby provide a sensitive probe of the barriers for the dihedral potentials. In practice, a direct comparison of simulation and experimental D and T1 is confounded by inaccuracies in the viscosity of the water model, and sensitivity to cutoff methods. For example, η20 (the viscosity η at 20 °C) for TIP3P water was shown to equal 0.62 and 0.35 cP when force-switching (12 Å cutoff) and Particle Mesh Ewald (PME), respectively, were used to treat the electrostatic interactions;3 the experimental value is 1.00 cP. More recent evaluations of η20 with PME and 3D-IPS/DFFT to treat the Lennard Jones (LJ) interactions yielded a value of 0.33 cP.4 Hence, a direct comparison of simulated and experimental results is clearly of limited usefulness because the results will be dominated by the errors in the water model.
Rotational and translational diffusion constants for polar molecules in water may be assumed to scale as η−1 under most conditions. Hence, it is reasonable to propose that the viscosity in the simulation system could be scaled to experiment when comparing diffusion constants. This paper considers cases where a scaling approach is appropriate, and when it may lead to fortuitous agreement with experiment. The idea is simple: simulated diffusion constants and rotational relaxation times for solutes in water may be scaled by the ratio of pure experimental and simulated water in the regime where the concentration dependences of viscosities match. The underlying assumption of such a scaling is that differences with experiment arise from deficiencies in the water model, not the solute model.
Following a description of the Methods (Section 2), Section 3.1 reports the simulated viscosities for pure TIP3P water at a range of temperatures, and for glucose and trehalose at a range of concentrations. These provide upper limits for the analysis of simulated and experimental D and T1 presented for glucose, trehalose, maltose and melibiose (Fig 1) in Sections 3.2 and 3.3. Section 3.4 presents the results of hydrodynamic calculations with varying levels of bound water to examine differences in the preceding carbohydrates. The Discussion and Conclusions (Section 4) includes examples where scaling is not advised.
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 ps−2 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 (Table 6).
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.
Viscosities were obtained from
where Pαβ 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αβ = 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:
where L 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 13C nucleus and its N attached protons, the 13C NMR relaxation rates 1/NT1 for each carbon are given by19
where is Plank’s constant divided by 2π, rC-H is the effective C-H bond length, γH, γC, ωH, and ωC are the gyromagnetic ratios and Larmor frequencies, respectively, of the 13C and 1H nuclei; ωC = γCH and ωH = γHH, where H is the field strength. J(ω) is the spectral density of the second rank reorientational correlation function:
where P2 is the second order Legendre polynomial (t) 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) reduces to
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.
Table 1 compares the simulated viscosities for TIP3P water from 20 to 50 °C. The values at 20 °C, 0.329 cP for PME and 0.332 for IPS, are within statistical error to those obtained previously.4 New results here demonstrate that the temperature dependence of TIP3P also underestimates experiment. The activation energies Ea obtained from a plot of −lnη vs. 1/T (Fig 2) are 1.8 and 3.8 kcal/mol for PME and experimental viscosities, respectively. Interactions between real waters are stronger than those in the TIP3P water model. Consequently, different scalings must be applied at different temperatures to compensate for deficiencies in the model.
The viscosity of a solution increases with the concentration of solute. Figure 3 plots experimental and simulated viscosities (scaled by the factors in Table 1) for solutions of glucose at 20 °C (top) and trehalose at 30 °C (bottom). Errors for glucose are less than 17% up to 27 wt% (2 m), but increase to 39% for 48 wt% (5 m). Agreement for scaled simulated values and experiment is comparable for trehalose as a function of wt%. Densities of glucose solutions are within 1% of experiment up to 5 m,1 indicating that the structural properties of the FF are more robust than the dynamical ones.
Based on these results, it is reasonable to propose that simulated diffusion constants and rotational relaxation times for solutes in water may be scaled by the ratio of pure experimental and simulated water in the regime where the concentration dependences of viscosities match; i.e., up to about 20 wt% (2 m for glucose and 1 m for trehalose) for mono and disaccharides. The underlying assumption of such a scaling is that differences with experiment arise from deficiencies in the water model, not the solute model. Scaling at higher concentrations is not advisable because agreement with experiment may be fortuitous. Concentrations used for NMR measurements are typically less than 10 wt%, so the preceding restriction is not severe.
As expected from the preceding discussion of viscosity, simulated diffusion constants for glucose and trehalose evaluated substantially overestimate experiment. However, as evident in Figure 4 and Table 2, the agreement at low concentration is excellent when the simulated results are scaled by the viscosity factors listed in Table 1. In particular, the D0 values (the diffusion constant extrapolated to zero concentration) are within 3% of experiment. Table 2 contains additional results for melibiose and maltose, confirming this good agreement. Simulation and experiment deviate for glucose and trehalose as concentration increases, though in opposite directions (Fig 4). The errors for trehalose appear to be primarily associated with the viscosity; i.e., the (scaled) simulated viscosity at 30 wt% is too low by a factor of 1.5 (Fig 3) and the simulated diffusion constant is too high by a factor of 1.7. In contrast, the simulated viscosities and diffusion constants of glucose are both too low.
Concentration dependent viscosities and diffusion constants (data not shown) for glucose at 20 °C are nearly identical for PME and IPS. This parallels results obtained earlier,4 and and further supports the notion that that long-range Lennard-Jones forces do not play a significant role in the dynamics of these systems. Long range electrostatics, of course, is critical, as evidenced by the factor of two difference in water viscosities with and without PME.3
Eq. (5) has been used directly to evaluate T1’s of alkanes from simulation.21 However, most measurements of carbohydrates are obtained in D2O, so in addition to the scaling applied above for diffusion constants, an additional scaling for the increased viscosity of D2O is required; this is 1.23 near room temperature (Table 1).
Table 3 lists τ for α and β-glucose obtained by integration of the reorientational correlation function, simulated 1/T1 calculated from Eq. (5) with the scaling as noted above, and the experimental 1/T1; the statistical error in the simulated relaxation times, and therefore the 1/T1, is approximately 4%. Results are reasonable for all carbons, with an average difference between simulation and experiment of 12%, not substantially larger than the error of the measurements (usually taken to be approximately 10%). Note that the relaxation times for the two CH vectors of C6 (the exocyclic methyl group) are only slightly less than those of the ring carbons. This implies that on the tumbling time of glucose, there is libration, but no isomerization of the exocyclic torsion. This is consistent with the relatively high barrier leading to isomerizations on the ns, but not ps, time scale.1
The NMR study by Engelsen et al.22 of maltose and trehalose provides both diffusion constants and T1’s, allowing a comparison of both observables obtained under the same conditions. Table 4 lists the T1 results. As for the translational diffusion constants (Table 2), agreement of simulation and experiment is near quantitative. Maltose was also simulated at the experimental conditions of Benesi and Brant23 (the same study reporting the glucose 1/T1). Agreement is comparable to that in Table 1 (data not shown).
The excellent results obtained for translational and rotational diffusion at low concentration for glucose and the disaccharides allow an estimate of the number of bound waters for each carbohydrate. Hydrodynamic calculations of diffusion constants of polar solutes in water based on Stokes’ Law require the addition of bound waters, or a hydration shell, to yield agreement with experiment.24 The same requirement holds for bead model calculations,15,16 as is evident in column 3 of Table 5: the calculated averaged diffusion constants of carbohydrates with no bound water overestimate experiment by approximately 50%. The lower value of <D> for melibiose compared to trehalose and maltose is consistent in that the former has a 1–6 linkage, and is therefore more extended than the others.
Waters may be placed on a solute in a variety of ways, though an energy cutoff is most natural when using conformations from simulations. In previous estimates, a range of energy cutoffs (Ecut) for carbohydrate-water interaction based on integer multipliers of Ecut/kBT was presented,16,17 with values at 3 and 4 usually bracketing experiment. A slightly different approach is adopted here, in that non-integer multipliers between -3 and -4 kBT with an interval of 0.1 were varied to provide the best agreement (a difference of 0.02×10−6 cm2/s or less) to experiment for each molecule; the averages are denoted <D>E and are listed in Table 5. For glucose, <Nb> =5.1, almost exactly one water per free hydroxyl group. <Nb> varies from 9 to 12 for the disaccharides, though each has only 8 free hydroxyl groups. Trends in <SASA>, Ecut, and <Nb> are readily apparent when these results are listed in order of decreasing D. Melibiose, with the more extended and flexible 1–6 linkage, exhibits the largest value for <Nb>, consistent with the slower diffusion and larger surface area. The noticeable difference in Ecut suggests that the additional bound waters have a weaker interaction with the sugar than those directly interacting with the hydroxyl groups. <Nb> is larger for trehalose than maltose, which is consistent with the lower experimental values for D.25–27 While a range of values is reported in the preceding experimental studies, depending on the method employed, in every case trehalose bound more water than maltose. Fig. 5 compares the hydration shells of these two disaccharides. It is possible that the rigidity of trehalose (with its 1–1 α linkage) allows the formation of a more structured water layer than maltose, allowing the inclusion of more water. This observation is consistent with the average number of waters in different size shells for two disaccharides (Table 6), and the average number of hydrogen bonded waters: 16.2 for maltose and 16.7 for trehalose. It is also consistent with the results from MD simulations listed in Tables 2 and and44 for translational and rotational diffusion, respectively. In each case, trehalose has the larger friction, even though the translational diffusion constant with no bound water is nearly the same as that of maltose (Table 5).
The empirical potentials presently used in atomic level simulations are approximate, carbohydrate-water interactions are complex, and transport properties are generally more difficult to model than equilibrium properties. In the present case, the low viscosity of the TIP3P water model (Table 1 and Fig. 2) complicates the validation of the carbohydrate parameter set using experimental results for translational diffusion and rotational relaxation. Simply scaling the viscosity of the solution to experiment at all concentrations is not advisable because discrepancies with experiment may not only arise from the water model: as illustrated in Fig. 3, the viscosities of glucose and trehalose solutions deviate from experiment at concentrations above 20 wt% even after scaling. Given the wide use of TIP3P in the simulation community, it is essential to understand its deficiencies and to develop a procedure for comparing experiment and simulation.
Based on the results in Fig. 3, it is proposed that viscosity scaling is permissible for solutions of mono and disaccharides of less than 20 wt% for simulations using TIP3P water This leads to excellent agreement with experiment for diffusion constants (Table 2) and NMR 13C T1’s (Tables 3 and and4),4), and thereby further validates the CHARMM parameter set C35. NMR experiments are usually carried at low concentration, so the limit of 20 wt% is not a severe restriction, and C35 can be applied with confidence to analyze relaxation measurements. While the approach presented here can be applied to other water models, the upper range of scaling will need to be determined for each.
The preceding results allow a hydrodynamic analysis of bound waters at low solute concentration (Table 5). While glucose binds one water per free hydroxyl, binding varies for disaccharides. It is proposed that the rigidity of trehalose allows the formation of a more structured hydration shell than maltose (Fig 5), allowing for binding of additional water (Table 5 and and6).6). Simulations of trehalose solutions at higher concentrations could in principle yield insight into the cryopreservative behavior of this molecule.28 However, the calculated diffusion constants for glucose and trehalose are incorrect at high concentration even after scaling (Fig 4). This result indicates deficiencies in the parametrization, and implies that simulations at higher concentration should be approached with extra caution.
It is possible to avoid viscosity scaling in some circumstances, and still reasonably compare simulation and experiment. For example, the hydrodynamic radius of nearly spherical particles could be evaluated because the different viscosities can be incorporated into Stokes’ Law. However, this is unsatisfactory for nonspherical particles. Likewise, relative relaxation times for intramolecular vectors can be compared without scaling, though reporting calculated T1’s without scaling is potentially confusing. There are circumstances in which scaling should be avoided. In addition to the high concentration regime already discussed, the possibility that the internal motions do not track viscosity should be considered. While there are interesting examples of translation29,30 and rotation31 in some organic solvents with more complicated viscosity dependencies, it is reasonable to assume that the translation and rotational frictions are proportional to water viscosity near room temperatures. However, the situation is less clear for isomerization, where a rapid change in an internal degree of freedom may take place within a cage of solvent, allowing the mode to partially uncouple from viscosity. For example, folding of some α-helical peptides follows a η−0.6 (or “non-Kramers”) dependence in water/sucrose solutions.32 Examples from simulation include spinning of a model methyl group in water,33 and gauche-trans isomerization in alkanes.34
In conclusion, viscosity scaling is reasonable to apply in a defined set of circumstances. It is helpful for validating parameter sets and for obtaining physical insight. It should not be applied where the simulated concentration dependent viscosity departs significantly from experiment, or where internal motions might not track viscosity.
This research was supported in part by the Intramural Research Program of the NIH, National Heart, Lung and Blood Institute, by NIH grant GM070855 to ADM Jr., and utilized the high-performance computational capabilities at the National Institutes of Health, Bethesda, MD (NHLBI LoBoS cluster).