Optimizing Nonbonded Interactions of the OPLS Force Field for Aqueous Solutions of Carbohydrates: How to Capture Both Thermodynamics and Dynamics.

Knowledge on thermodynamic and transport properties of aqueous solutions of carbohydrates is of great interest for process and product design in the food, pharmaceutical, and biotechnological industries. Molecular simulation is a powerful tool to calculate these properties, but current classical force fields cannot provide accurate estimates for all properties of interest. The poor performance of the force fields is mainly observed for concentrated solutions, where solute-solute interactions are overestimated. In this study, we propose a method to refine force fields, such that solute-solute interactions are more accurately described. The OPLS force field combined with the SPC/Fw water model is used as a basis. We scale the nonbonded interaction parameters of sucrose, a disaccharide. The scaling factors are chosen in such a way that experimental thermodynamic and transport properties of aqueous solutions of sucrose are accurately reproduced. Using a scaling factor of 0.8 for Lennard-Jones energy parameters (ϵ) and a scaling factor of 0.95 for partial atomic charges ( q), we find excellent agreement between experiments and computed liquid densities, thermodynamic factors, shear viscosities, self-diffusion coefficients, and Fick (mutual) diffusion coefficients. The transferability of these optimum scaling factors to other carbohydrates is verified by computing thermodynamic and transport properties of aqueous solutions of d-glucose, a monosaccharide. The good agreement between computed properties and experiments suggests that the scaled interaction parameters are transferable to other carbohydrates, especially for concentrated solutions.

Introduction

Saccharides define a class of carbohydrates considered vital for a wide range of biological and industrial processes.[1,2] Examples include the role of (poly-)sacharides or glyco-proteins in cryo- and lyopreservation of biomaterials and foods,[3−8] the regulation of the osmotic pressure in living cells,[9] the interplay between saccharides and proteins in the cell membrane,[10,11] and use as a feedstock for the production of biofuels and renewable chemicals in biotechnology.[1] Either in biology or in industrial applications, saccharides predominantly occur in aqueous solutions. To obtain a better understanding of the underlying mechanisms defining the biological function of a saccharide, or to obtain a better description of the thermophysical properties of saccharide solutions (as needed in technological applications), requires a molecular-level description of the interactions between sugars, water, and water–sugar. Classical molecular simulations provide a valuable tool for this; however, it is essential that the used force fields reproduce relevant volumetric, structural, and dynamic properties of the system of interest. For aqueous solutions of saccharides, this has proven a difficult task.

As several studies[12−17] have shown, a common problem when applying carbohydrate force fields such as CHARMM,[18−21] GLYCAM06,[22] GROMOS,[23,24] or OPLS-AA[28] to saccharide solutions is the overestimation of the sugar–sugar interactions, leading to strong sugar aggregation at elevated sugar concentrations. As a result, both thermodynamic properties (densities, activity coefficients, thermodynamic factors, second virial coefficients, etc.) and transport properties (shear viscosities, diffusion coefficients, etc.) of concentrated solutions are generally not well reproduced.[12−15,17] Recent work showed that rescaling the sugar–sugar nonbonded interactions to lower values improves the description.[13,15,29] Basically, two strategies have been proposed. In the first approach (e.g., the work of Sauter and Grafmüller[14] and Lay et al.[15]), only the dispersion interactions—more precisely, the energy parameters for some of the Lennard-Jones (LJ) interactions between sugars—are modified. In the second approach (e.g., Batista et al.[13]), only the partial atomic charges of the sugars are scaled. The first approach improves the description of thermodynamic properties, such as the osmotic pressure and second virial coefficient.[14,15] However, transport properties appear to be less affected by solely modifying the dispersion interactions (see the Supporting Information of Lay et al.[15] for results on the viscosity of glucose solutions). Results obtained using the second approach indicate that rescaling partial atomic charges leads to considerable improvements in the description of transport properties and suggest that thermodynamic properties (e.g., density) are also better described. The results of Batista et al.[13] are not conclusive, however, as mixture properties known to be more sensitive to changes in the partial charges of the atoms (e.g., thermodynamic factors and osmotic coefficients) were not taken into account. Also, the only transport property studied was the viscosity, and diffusion coefficients were not considered.

In this work, we first present a more complete analysis on how the properties of aqueous carbohydrate solutions change when the LJ or electrostatic interactions are scaled. The force fields of Batista et al.[13] (OPLS with scaled partial atomic charges) and Lay et al.[15] (GLYCAM06 with scaled LJ interactions) are analyzed in more detail. We focus on the disaccharide sucrose (Figure ) as a relevant test case, as this molecule has many applications, ranging from biology to food science.[3,4,6,7] Our results show that only scaling partial charges is insufficient to reproduce both thermodynamic and transport properties of concentrated aqueous saccharide solutions. Scaling only the dispersion interactions leads to significantly better agreement, with qualitative improvements in transport properties, as well as thermodynamic properties. Although scaling the LJ interactions improves the force field, our analysis indicates that obtaining quantitative agreement between experiments and computed properties of saccharide solutions may require a global optimization procedure of the force field, with different scaling factors for LJ interactions of different atom pairs. Since the choice of atom pairs is rather arbitrary, and an optimization of many force field parameters is difficult, we here propose a different approach, in which all LJ energy parameters and all partial atomic charges of the sugars are scaled simultaneously.[29] This procedure only requires optimization of two parameters: one scaling factor for LJ interactions and one scaling factor for partial atomic charges. For the optimization of sucrose nonbonded interaction parameters, the OPLS force field[25,26] is used as a basis. The transferability of the refined OPLS force field is analyzed for d-glucose as a prototypical example of monosaccharides.

Figure 1

Schematic and atomistic representations of sucrose. The atomistic representation was created by using iRASPA.[79]

This manuscript is organized as follows. In section 2, the details of simulations and equations used for computing transport and thermodynamic properties are provided. In section 3, a set of optimum scaling factors for nonbonded interaction parameters of the OPLS force field is obtained and the accuracy of the refined force field is investigated. Our conclusions are summarized in the last section.

Simulation Details

Force-field-based equilibrium MD simulations were carried out to compute thermodynamic properties (i.e., liquid densities and thermodynamic factors) and transport properties (i.e., shear viscosities, self-, and Fick diffusion coefficients) of water–sucrose mixtures at 298 K and 1 atm. The software package LAMMPS[30] (version of February 16, 2016) was used to perform the MD simulations. The molecular configurations of aqueous sucrose solutions were made in Packmol.[31] We considered five systems, with sucrose mass fractions (wsucrose) of 20%, 30%, 40%, 50%, and 60%. These systems contained 6, 9, 12, 16, and 20 sucrose molecules combined with 456, 399, 342, 304, and 253 water molecules, respectively. These systems yield homogeneous liquid solutions, as sucrose crystallizes in mixtures with a mass fraction above 67% at 298 K.[32] The input files for these configurations are produced using VMD.[33] The same procedure was followed to construct LAMMPS input files for four aqueous glucose solutions. To create mixtures with glucose mass fractions (wglucose) of 20%, 30%, 40%, and 50%, the systems contained 12, 18, 24, and 30 glucose molecules combined with 480, 420, 360, and 300 water molecules, respectively. At 298 K, α-d-glucose monohydrate forms in aqueous glucose solutions with mass fractions above 51%.[34]

Two groups of force fields are considered for sucrose: GLYCAM06[22] and OPLS-AA.[25,26] The three-site SPC/Fw water model is used.[35] Nonbonded LJ interactions are truncated at a cutoff radius of 9 Å, and analytic tail corrections are included for the computation of energies and pressures.[36] The standard mixing rules of the OPLS force field (see the Supporting Information) are applied for nonbonded LJ interactions between different atoms in different molecules, or atoms of the same molecule separated by at least three bonds.[25,36] The Lorentz–Berthelot mixing rules are considered for the GLYCAM06 force field.[36] The 1–4 nonbonded intramolecular interactions are scaled by a factor of 0.5 for the OPLS force field[25] and a factor of 0 for the GLYCAM06 force field.[22] Long-range electrostatic interactions are taken into account using the PPPM method with a relative precision of 10–6.[36] The velocity-Verlet algorithm with a time step of 1 fs is used to integrate the equations of motion.[36]

Liquid densities are computed from an NPT ensemble, in which the temperature and pressure of the system are controlled by using the Nosé–Hoover thermostat and barostat.[36] All transport properties and thermodynamic factors are computed in an NVT ensemble, where the temperature is kept fixed by using the Nosé–Hoover thermostat.[36] The lengths of each simulation for computing transport and thermodynamic properties are 200 ns and (at least) 10 ns, respectively. To assert the statistical uncertainties, five independent simulations were performed for each condition.

The Einstein approach to calculating transport properties in equilibrium MD simulations is used.[36−38] The shear viscosity is obtained from the off-diagonal components of the stress tensor (P), according to[36,38,39]where kB and t are the Boltzmann constant and time and V and T are the volume and temperature of the system. The angle brackets ⟨...⟩ denote an ensemble average. The self-diffusion coefficient of species i is calculated based on the mean-square displacement, as[36,37]in which N is the number of molecules of species i and r is the position of jth molecule of species i. The finite-size effects of self-diffusion coefficients are taken into account by including the analytic correction proposed by Yeh and Hummer.[40,41]

There are two approaches to defining the mutual diffusion coefficient of a binary mixture:[42−44] the well-known Fick diffusion coefficient (DFick) is mainly used in industrial applications, since Fick’s law relates the mass flux to gradients in concentrations, which can be measured experimentally. An equivalent formulation is based on the Maxwell–Stefan (MS) diffusion coefficient,[45] which relates the mass flux to gradients in chemical potentials. For binary mixtures, Fick and MS diffusion coefficients are related by the thermodynamic factor, Γ, according to[42,44,46−49]with Γ defined as[42,46,49−51]where x1 and γ1 are the mole fraction and activity coefficient of species 1, respectively. For most binary mixtures (including those considered in this work), the deviations of the thermodynamic factor and solute activity coefficient from 1 are of opposite sign.[50,52] For such cases, Γ > 1 determines the case that the interactions between different species are favored, whereas 0 < Γ < 1 determines the case that interactions between similar species dominate. The thermodynamic factor can thus be used as a measure for the degree of self-association of solutes. In some studies on water–carbohydrate solutions,[15,16] the osmotic coefficient (Φ = 1 + ln(γsolvent)/ln(xsolvent)) is used instead. For the solutions studied, the osmotic coefficient contains the same information as the thermodynamic factor, with Φ > 1 determining the case that the interactions between different species are favored, while Φ < 1 determines the case that the interactions between similar species are strongest.

Fick diffusion coefficients can be obtained both from equilibrium or from nonequilibrium MD simulations.[36,44,53−55] Here, we calculate Fick diffusivities based on eq , with MS diffusion coefficients and thermodynamic factors obtained from equilibrium MD simulations.[44,56,57] Further details can be found in the review article of Liu et al.[44]

The MS diffusivity can be computed by calculating the Onsager coefficients (Λ):[44,56]The MS diffusivity of a binary system then follows from[44,51]The finite-size effects of Fick and MS diffusion coefficients are corrected using the analytic relation proposed by Jamali et al.[58] The finite-size effects of MS diffusion coefficients can be as large as the magnitude of the computed MS diffusivities, depending on the nonideality of the mixture.[58] As the modified force field is verified by comparing the mutual diffusivities with experiments, it is crucial to consider these finite-size effects and compute the mutual diffusion coefficients in the thermodynamic limit. To improve the sampling of the correlation functions in eqs , 2, and 5, an order-n algorithm is used.[37,59]

The thermodynamic factor is calculated from the so-called Kirkwood–Buff coefficients (G), which are obtained by integrating radial distribution functions over volume[60,61] and applying appropriate finite-size corrections.[62−65] For a binary mixture, the thermodynamic factor can be written as[44]where c1 is the number density of species 1. The finite-size effects of radial distribution functions are corrected for according to the work of van der Vegt and co-workers.[62,63] Kirkwood–Buff coefficients (G) are computed in the thermodynamic limit using the method of Krüger et al.[64,65]

Results and Discussion

In this section, we first analyze several procedures for optimizing force fields of carbohydrate solutions. The quality of a procedure is measured by the accuracy with which thermodynamic factors and shear viscosities are estimated. Once the optimal procedure is established, we continue by developing and verifying a refined OPLS force field for aqueous solutions of sucrose. The transferability of the refined force field to aqueous solutions of other carbohydrates is examined by considering glucose. All computed properties presented in this section, and the force field parameters of both the original and refined OPLS force fields, are listed in the Supporting Information.

Analysis of Different Optimization Procedures

Batista et al.[13] studied the effects of nonbonded electrostatic interactions in the GROMOS 56ACARBO and OPLS force fields for aqueous solutions of d-glucose. This study concluded that the properties of diluted solutions can be predicted accurately based on the 56ACARBO force field. At high glucose mass fractions (wglucose > 0.60), the overestimated self-association of glucose molecules was shown to result in overestimation of shear viscosities by up to 650%.[13] To decrease the solute–solute interactions, all partial atomic charges of the glucose molecule were scaled by a factor 0.8. The estimation of shear viscosities of concentrated water–glucose mixtures was thereby improved significantly, with a maximum deviation from experiments of 27%.[13] Besides an analysis on the density of the solution, these authors did not consider thermodynamic properties such as osmotic coefficients and thermodynamic factors to verify the performance of the optimized force field.

Inspired by the promising results of Batista et al.,[13] we study a modified OPLS force field for which the partial atomic charges of sucrose are scaled by a factor 0.8. Similar to the work of Batista et al.,[13] the LJ interactions parameters remain unchanged. In Figure , thermodynamic factors and shear viscosities computed based on the original OPLS force field and the OPLS force field with scaled partial atomic charges are compared to experimental data at 298 K and 1 atm. In agreement with the results of Batista et al.,[13] the overestimation of shear viscosities is decreased by scaling the electrostatic interactions. This is especially the case for mixtures of high concentrations of sucrose, where sucrose–sucrose interactions mainly determine the properties of the mixture. Thermodynamic factors are less accurately reproduced. Our results show qualitatively incorrect changes as compared to the original OPLS force field, with larger deviations of thermodynamic factors from experiments. Estimated thermodynamic factors are closer to zero, meaning the self-association of sucrose has increased—not decreased. This suggests that merely scaling the partial atomic charges is insufficient to reproduce both thermodynamic and transport properties of carbohydrate solutions.

Figure 2

Computed (a) shear viscosities and (b) thermodynamic factors for water–sucrose mixtures as a function of the mass fraction of sucrose. The properties are computed based on the OPLS force field[25] (blue circles) and the OPLS force field with scaled partial atomic charges[13] (red squares). The SPC/Fw water model was used.[35] Lines represent (a) experimental shear viscosities (solid,[80] dashed[81]) and (b) thermodynamic factors (solid,[82] dashed[83]). Error bars indicate 95% confidence intervals. No error bars are reported for the thermodynamic factors.

Another approach to optimizing a force field is to modify the LJ parameters of solute molecules. As a good rule of thumb, the size (σ) and energy (ϵ) parameters can be modified to improve the predictions of volumetric (e.g., liquid densities) and thermal (e.g., heat of vaporization) properties, respectively.[2,23,27,66] In recent modifications of the MARTINI, CHARMM, and GLYCAM06 force fields for carbohydrate solutions,[15,17] only the energy parameters were changed. An example is the work of Lay et al.,[15] in which the CHARMM and GLYCAM06 force fields were modified by altering some of the LJ energy parameters (ϵ) defining the solute–solute interactions in aqueous solutions of glucose. On the basis of introducing modifications to carbon–carbon and carbon–oxygen interactions, accurate predictions of osmotic coefficients were obtained. However, shear viscosities were underestimated by about 60%.[15]

The modified GLYCAM06 force field of Lay et al. is analyzed in more detail in Figure . In this figure, computed thermodynamic factors and shear viscosities of water–sucrose mixtures are compared to results obtained based on the original GLYCAM06 force field. The modified energy parameters considerably improve the description of both properties, for a wide range of sucrose concentrations. Shear viscosities are especially well reproduced, showing a good agreement with experiments. We note that the improved results for viscosity as compared to those reported by Lay et al.[15] may be due to using the SPC/Fw water model instead of TIP3P. The TIP3P water model is known to provide poor estimates for the shear viscosity of water.[67−69] The SPC/Fw water model is one of the most accurate three-site models for estimating transport properties of water under ambient conditions.[35,69,70] We chose a three-site model due to the smaller computational requirements compared to more complex models such as four-site (e.g., TIP4P/2005[71]), five-site (e.g., TIP5P-E[72]), and polarizable (e.g., HBP[73]) water models. An extensive comparison of the performance of different water models can be found in the literature.[66,68,69,74] Although thermodynamic factors are not as well described as shear viscosities, the observed increase in the thermodynamic factor is at least qualitatively correct. These results suggest that modifying the LJ energy parameters of solute molecules could be a sufficient procedure to develop a force field that accurately reproduces both thermodynamic and transport properties of carbohydrate solutions.

Figure 3

Computed (a) shear viscosities and (b) thermodynamic factors for water–sucrose mixtures as a function of the mass fraction of sucrose at 298 K and 1 atm. The properties are computed based on the GLYCAM06 force field[22] (blue circles) and a modified GLYCAM force field as proposed by Lay et al.[15] (red squares), in which some of the LJ sugar–sugar interactions are modified. The SPC/Fw water model was used.[35] Lines represent (a) experimental shear viscosities (solid,[80] dashed[81]) and (b) thermodynamic factors (solid,[82] dashed[83]). Error bars indicate 95% confidence intervals. No error bars are reported for the thermodynamic factors.

While the approach of Lay et al. seems effective, this approach may require the modification of the LJ energy parameters of many different atom pairs, each with its own scaling factor. The choice of atoms whose ϵ is changed requires a profound understanding of the interactions between the atoms of the carbohydrate, which is not always known a priori. This, combined with the fact that the optimization of many force field parameters usually proceeds by trial and error (and is thus computationally expensive), leads us to suggest a different approach. On the basis of the results shown in Figures and 3, we propose scaling the partial atomic charges and LJ energy parameters of all atoms of the carbohydrate. This can simultaneously improve the description of thermodynamic and transport properties, without requiring different scaling factors for different atom pairs. We propose to scale the nonbonded interaction parameters in such a way that experimental data for at least one thermodynamic property (e.g., thermodynamic factors) and one transport property (e.g., shear viscosities) can be reproduced by the force field. The advantage of this approach over the method of Lay et al.[15] is a considerable reduction in the number of modification parameters to two: a single scaling factor for all LJ energy parameters of the carbohydrate and a single scaling factor for all partial atomic charges of the carbohydrate.

Refined OPLS Force Field

The optimization procedure proposed in the previous section is here applied to the OPLS force field, for aqueous solutions of sucrose. As in the previous section, thermodynamic factors and shear viscosities are considered as target properties. The following scaled LJ energy parameters (first number) and scaled partial atomic charges (second number) are considered: 1.0 and 1.0 (corresponding to the original OPLS force field), 1.0 and 0.8 (corresponding to the work of Batista et al.[13]), 0.8 and 1.0, and 0.8 and 0.8.

In Figure , computed thermodynamic factors and shear viscosities are compared to experimental data. In agreement with what is shown in Figure , scaling all partial atomic charges by a factor less than 1 decreases both the thermodynamic factors and shear viscosities. Scaling all LJ interactions by a factor less than 1 decreases the shear viscosities while increasing the thermodynamic factors. Adjusting LJ energy parameters can thus correct for a strong decrease in thermodynamic factors caused by scaling down partial atomic charges (as in Figure ). Of these four combinations, the best combination is for the scaling factors of 0.8 and 1.0 for the LJ energy parameters and partial atomic charges, respectively. By interpolating between the computed thermodynamic factors and shear viscosities, corresponding to the four combinations of scaling factors, the optimum combination is a scaling factor of 0.8 for the LJ energy parameters and 0.95 for the partial atomic charges. Hereafter, the OPLS force field with this combination of scaling factors will be referred to as the refined OPLS force field. The interaction parameters of this force field (as well as the original OPLS force field) are listed in the Supporting Information.

Figure 4

Computed (a) shear viscosities and (b) thermodynamic factors for water–sucrose mixtures as a function of the mass fraction of sucrose. All LJ energy (ϵ) parameters and partial atomic charges (q) of the OPLS force field[25] are scaled by respectively 0.8 and 0.8 (green diamonds), 0.8 and 1.0 (magenta pentagons), 1.0 and 0.8 (red squares), and 1.0 and 1.0 (blue circles). The SPC/Fw water model was used.[35] Lines represent (a) experimental shear viscosities (solid,[80] dashed[81]) and (b) thermodynamic factors (solid,[82] dashed[83]). Error bars indicate 95% confidence intervals. No error bars are reported for thermodynamic factors.

To verify that the refined OPLS force field outperforms the original OPLS force field, several thermodynamic and transport properties of the water–sucrose mixtures are computed at 298 K and 1 atm. In Figure , computed (a) liquid densities, (b) thermodynamic factors, (c) shear viscosities, (d) self-diffusivities of sucrose, (e) water, and (f) Fick diffusion coefficients are compared to experimental data for a wide range of concentrations. As compared to the original OPLS force field, thermodynamic properties (i.e., liquid densities and thermodynamic factors) computed using the refined OPLS force field show a considerably better agreement with experiments. In Figure a, excellent agreement between computed liquid densities and experimental data is observed. Since the size (σ) parameters of the LJ potentials are unaltered, the differences between liquid densities computed using the original and refined OPLS force fields are small. These differences reach a maximum value of 1%, at a sucrose mass fraction of 60%. Since the original OPLS force field overestimates solute–solute interactions at high concentrations of sucrose, the sucrose molecules tend to aggregate and form more packed liquid structures. This increases the density of the liquid, leading to overestimation of computed densities. Due to this overestimation of solute–solute interactions by the original OPLS force field, the thermodynamic factors are underestimated (see Figure b). The good agreement between computed thermodynamic factors and experimental data confirms an accurate description of the solute–solute and solute–water interactions by the refined OPLS force field.

Figure 5

Computed (a) densities, (b) thermodynamic factors, (c) shear viscosities, (d) self-diffusivities of sucrose, (e) self-diffusivities of water, and (f) Fick (mutual) diffusion coefficients for water–sucrose mixtures as a function of the mass fraction of sucrose. Temperature and pressure are 298 K and 1 atm, respectively. The results of the original OPLS force field[25] are shown by blue circles. The results of refined OPLS force field (scaled LJ energy parameters and partial atomic charges of all atoms of sucrose by factors of 0.8 and 0.95, respectively) are shown by red squares. Lines represent experimental data for (a) densities (solid,[84] dashed[85]), (b) thermodynamic factors (solid,[82] dashed[83]), (c) shear viscosities (solid,[80] dashed[81]), (d) self-diffusivities of sucrose (solid,[81] dashed[86]), (e) self-diffusivities of water (solid[81]), and (f) Fick diffusion coefficients (solid,[87] dashed[88]). Error bars indicate 95% confidence intervals.

As shown in Figure c–f, transport properties computed using the refined OPLS force field are in excellent agreement with experimental data. At low sucrose concentrations, both the original and refined OPLS force fields yield similar values for the transport properties. Deviations become significant as the concentration of sucrose increases. In agreement with the work of Batista et al.,[13] the shear viscosities (Figure c) computed based on the original OPLS force field are overestimated by up to 300% at a sucrose mass fraction of 0.6. For the refined OPLS force field, this deviation reaches a maximum of 30%. Figure d shows that the refined force field provides better estimates for the self-diffusion coefficients of sucrose at high sucrose mass fractions. This is due to less self-association, leading to higher mobility of sucrose molecules. In Figure e, no significant differences between self-diffusion coefficients of water computed using the original and refined OPLS force fields are observed. Both force fields were combined with the SPC/Fw water model, which mainly determines the self-diffusivity of water. As the concentration of sucrose increases, the sucrose–water interactions gradually become more important in determining the properties of the medium in which water molecules diffuse. This leads to the small differences observed between computed self-diffusivities at high sucrose mass fractions. In Figure f, Fick diffusion coefficients are compared with experiments. Similar to what is observed for the shear viscosity, Fick diffusivities computed using the refined OPLS force field are significantly closer to experimental values than those computed based on the original OPLS force field.

To verify whether the optimum scaling factors for LJ energy parameters (0.8) and partial atomic charges (0.95) are transferable to other carbohydrates, MD simulations were performed for aqueous solutions of the monosaccharide d-glucose at 298 K and 1 atm. Since mixtures of glucose with mass fractions above 51% produce α-d-glucose monohydrate, the following mass fractions of water–glucose mixtures are considered: 20%, 30%, 40%, and 50%. d-Glucose forms two anomers in an aqueous solution.[75] Simulation input files are constructed for a ratio 1:2 of α-d-glucose to β-d-glucose.[15,27,76,77] Lay et al.[15] observed that the exact value of this ratio does not affect the outcome of MD simulations.

In Figure , computed thermodynamic and transport properties of water–glucose mixtures are compared with experiments. The order of the subfigures is the same as in Figure . The results shown in Figure a and b clearly show that, compared to the original OPLS force field, the refined OPLS force field leads to a significantly better description of thermodynamic properties. Computed thermodynamic factors are larger than 1 and close to experimental values; the unphysical self-association of glucose molecules at elevated concentrations is thus resolved by using the refined OPLS force field. Figure c–f shows the transport properties for glucose–water mixtures (i.e., the shear viscosity, self-diffusivity of glucose and water, and Fick diffusion coefficient). A good agreement between experimental data and computed transport properties is observed. At low mass fractions of glucose (wglucose ≤ 0.4), the original OPLS force field shows slightly better agreement with experimental data. However, it is important to note that, for these diluted solutions, the difference between the results obtained based on these two force fields is at most 40% (for self-diffusivity of glucose at a mass fraction of 40%). This is mainly due to the low concentration of glucose, which decreases the effect of the carbohydrate force field on mixture properties. As the mass fraction of glucose rises above 40%, the refined OPLS force field provides better estimates for transport properties. At these concentrations, the results computed based on the original OPLS force field deviate significantly from experiments.[13]

Figure 6

Computed (a) densities, (b) thermodynamic factors, (c) shear viscosities, (d) self-diffusivities of glucose, (e) self-diffusivities of water, and (f) Fick (mutual) diffusion coefficients for water–glucose mixtures as a function of the mass fraction of glucose. Temperature and pressure are 298 K and 1 atm, respectively. The results of the original OPLS force field[25] are shown by blue circles. The results of the refined OPLS force field (scaled LJ energy parameters and partial atomic charges of all atoms of glucose by factors of 0.8 and 0.95, respectively) are shown by red squares. Solid lines represent experimental data for (a) densities (solid,[89] dashed[90]), (b) thermodynamic factors (solid,[91] dashed[92]), (c) shear viscosities (solid,[78] dahsed[89]), (d) self-diffusivities of glucose (solid[90]), (e) self-diffusivities of water (solid,[90] dashed[93]), and (f) Fick diffusion coefficients (solid,[94] dahsed[95]). Error bars indicate 95% confidence intervals.

By increasing the temperature of the water–glucose mixture, the solubility of glucose in water increases, and mixtures with higher mass fractions of glucose can be realized. For these concentrated solutions, the refined OPLS force field is expected to lead to significant improvements over the original OPLS force field. To verify this, we computed the shear viscosity of an aqueous solution with a glucose mass fraction of 60%, at 313.15 K and 1 atm. The shear viscosities computed based on the original and refined OPLS force fields are 32 ± 4 cP and 15.0 ± 1.5 cP, respectively. Compared to an experimental value of 13.08 cP for this mixture,[78] the original OPLS force field leads to a considerable overestimation of the shear viscosity (of about 150%). The refined OPLS force field significantly improves the description at such a high glucose mass fraction.

Conclusions

For aqueous solutions of carbohydrates, most force fields predict self-aggregation of solute molecules due to the overestimation of solute–solute interactions. This overestimation results in large deviations of computed properties from experimental data, especially for concentrated solutions, which are relevant to food and biotechnological industries. To decrease the tendency of solute molecules to self-aggregate, we proposed scaling the Lennard-Jones energy parameters (ϵ) and partial atomic charges (q) of the OPLS force field. In this way, accurate estimates were obtained for thermodynamic and transport properties of aqueous solutions of sucrose with mass fractions in the range 20%–60%. For this optimization, the three-site SPC/Fw water model was used. No modification in the bonded and nonbonded interaction parameters of the water model was applied. The optimum scaling factors for the Lennard-Jones energy parameters and partial atomic charges were obtained by reproducing experimental thermodynamic factors and shear viscosities of aqueous sucrose solutions. These factors are 0.8 and 0.95, respectively. For both the original and refined OPLS force fields, MD simulations were performed to calculate thermodynamic properties (i.e., liquid densities and thermodynamic factors) and transport properties (i.e., shear viscosities, self-diffusivities of water and sucrose, and Fick diffusion coefficients). Excellent agreement is observed between the properties computed based on the refined OPLS force field and experiments, for a wide range of sucrose concentrations. The transferability of the optimum scaling factors was verified by performing MD simulations of aqueous solutions of d-glucose. The computed thermodynamic and transport properties agree well with available experimental data, especially at high concentrations of glucose. This suggests that the scaling factors are transferable to other carbohydrates. By using the refined OPLS force field, accurate estimates for thermodynamic and transport properties can be obtained. While the proposed method for optimizing the nonbonded interactions of a force field was verified for the OPLS force field, this method may also be used for other force fields combined with other water models.