Predicting the Chemical Potential and Osmotic Pressure of

Article pubs.acs.org/JCTC

Predicting the Chemical Potential and Osmotic Pressure of Polysaccharide Solutions by Molecular Simulations Jörg Sauter and Andrea Grafmüller* Theory and Bio-Systems, Max Planck Institute of Colloids and Interfaces, 14424 Potsdam, Germany S Supporting Information *

ABSTRACT: Differences in the chemical potential of water and the resulting osmotic pressure across semipermeable membranes are of fundamental importance for many biological systems. Here, we calculate the osmotic pressure and the chemical potential of water for polysaccharide solutions by molecular simulations. We set up a method to measure the osmotic pressure in polysaccharide systems at different concentrations and found that for monomers the experimental trend with respect to the solute concentration is reproduced correctly. However, the calculated osmotic pressure values are systematically too low, and two common carbohydrate force fields (FFs) cannot correctly describe the relationship between the osmotic pressure and the degree of polymerization. Therefore, we reparametrized parts of the GLYCAM06 TIP5P FF based on osmotic pressure data. The predictive power of the resulting GLYCAM06TIP5P OSMOr14 FF is demonstrated for two different sugar molecules over a wide range of concentrations, and additional evaluations for other solution properties show improved agreement with experimental data. Finally, we discuss different methods to obtain the chemical potential of water in solutions. the GLYCAM06 FF6 with TIP3P7 water and the CHARMM36 FF8 have trouble reproducing different solution properties. The best agreement with experimental data was found for the GLYCAM06 FF with TIP5P9 water. With increasing computational resources, and, more importantly, more elaborate simulation methods becoming available, MD simulations are successfully applied to address a wider range of problems. Several methods have been devised that allow the computation of the chemical potential of water or the osmotic pressure of a solution by MD simulations. Common methods to calculate the chemical potential are Widom Test Particle insertion10 (TPI) and enhancements based on the same principle that are designed to correct the shortcomings of the TPI method or to increase the sampling efficiency.11,12 Other methods include Inverse Widom13 as well as free energy methods such as Thermodynamic Integration14 (TI) and Bennets Acceptance Ratio15 (BAR). In addition, a “noble water method” combines TPI and TI16 in order to predict the water uptake in polymer networks. In addition, a class of methods that allow the determination of the osmotic pressure directly from molecular simulations has been developed. These methods mimic a semipermeable membrane and measure the interactions of the solutes with the virtual membrane. Methods of this type were introduced by Murad and Powles17 and Luo and Roux.18 Previous studies have shown that measurements obtained with this type of method can reproduce experimental data18 and offer a promising way to link the molecular interactions to

1. INTRODUCTION Polysaccharides in animals and plant cells predominantly exist in aqueous environments, and their interactions with water can profoundly influence the organisms at larger scales. Differences in the chemical potential of water and the resulting osmotic pressure across semipermeable membranes play a central role in many living organisms, and elaborate osmo-sensing and osmo-signaling mechanism are found in the cell membrane. The impact of glucose on the osmotic balance is so pronounced that special mechanisms to store glucose, in the form of starch in plants and glycogen in animals, have developed to minimize its osmotic impact. An unbalanced glucose concentration can lead to dehydration, as happens for example in diabetes, and can contribute to hyperosmotic stress which can be linked to many diseases.1 Another example of an elaborate osmosensitive mechanism is the interaction of hemicellulose polysaccharides in the plant cell wall, which play an important role in the directed swelling of plant cell walls driven by water uptake.2,3,4 A molecular level understanding of the interactions of polysaccharides with water may help to gain a better understanding of diseases associated with osmosis or may be useful in the development of actuators or adaptive materials and new methods to overcome the plant cell wall recalcitrance. Molecular dynamics (MD) simulations have been successfully used to gain insights into the molecular details of biomolecules for many years. As the role of polysaccharides is recognized in an increasing number of biological systems, MD simulations of polysaccharides are gaining both importance and accuracy in reproducing molecular conformations. However, a comparison of several carohydrate force fields (FF) in simulations of simple saccharide solutions has shown5 that © XXXX American Chemical Society

Received: March 22, 2016

A

DOI: 10.1021/acs.jctc.6b00295 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

where G is the Gibb’s free energy. The difference of μw from the ideal gas contribution is the excess chemical potential μex w , which is equivalent to the hydration free energy of water ΔG obtained from vapor pressure experiments or simulations.24,25 Knowledge of the chemical potential of water allows the computation of the activity of water aw by using

experimentally accessible observables. Osmotic data have been used to improve FF parameters for solutions of small ions18−20 and amine−phosphate interactions.21 The trend found in measurements of the osmotic pressure of GLY1, GLY2, and GLY322 was that the agreement with experimental data becomes worse with increasing molecule size, which would suggest poor performance for polysaccharide solutions. Nonetheless, very recently, it has been applied successfully to improve the performance of two carbohydrate FFs for reproducing solution properties.23 Elcock et al.23 addressed the aggregation problem reported5 for the GLYCAM TIP3P FF and show that for both the GLYCAM TIP3P FF and the CHARMM36 FF the correct osmotic coefficients for 1 M glucose solutions can be matched by systematically decreasing the ε Lennard-Jones parameters for the intersolute C−C and C−O interactions. The parameters obtained by matching the osmotic coefficient at 1 M were found to also improve the agreement with experiments at other concentrations, as well as for sucrose, and lead to good agreement with experimental data for the size of a dextran polymer; however, the shape of the concentration dependence still deviates from experimental data, and the osmotic coefficient for sucrose, i.e., a molecule with a degree of polymerization (DP) of two is significantly overestimated. To be able to predict the effects of small structural or chemical differences on the osmotic pressure of the system, a FF that can reproduce these dependencies closely is required. Here, we first determine the appropriate method parameters to apply the method of Luo and Roux to polysaccharide solutions to determine suitable method parameters for the application to larger solutes. Then, we test the GLYCAM06 TIP5P FF, which was best able to reproduce solution properties of small molecules, with respect to its ability to reproduce experimental osmotic pressure data. Similar limitations as discussed by Elcock et al.23 become apparent also for this FF. Therefore, we introduce a partial reparametrization for β-Dglucose and maltose for concentrations ranging from 1.0 to 6.17 m. The predictive power of the resulting GLYCAM06TIP5P OSMOr14 FF is demonstrated over a wide range of concentrations for β-Dxylose and maltotriose. To further assess the reliability of the new FF, the molecular conformations are compared to those found with the original FF. Solution properties such as mass density and diffusion coefficients are compared to experimental data. In addition, the aggregation behavior in terms of the saccharide−saccharide radial distribution functions (RDFs) in the original and the new FF are discussed. Finally, to relate the osmotic pressure in the system to the chemical potential of water μw, knowledge of the chemical potential of water in the pure water reference state is calculated. In that context, we compare two methods to calculate the chemical potential and briefly discuss the application of these methods to directly calculate μw in polysaccharide solutions.

μw − μw⊖ = RT ln aw

where μw is the chemical potential of water in the solute− solvent system, μ⊖ w the chemical potential of water in the pure water system in the standard state, and R denotes the gas constant. To compare to experimental data, the activity of water aw can be obtained from experimental osmotic coefficients Φ by Φ=−

p

ln(aw) M wm

(3)

where Mw is the molar mass of water, and m the molality of the solute. The osmotic pressure Π is related to the water activity by Π=−

RT ln(aw) Vw

(4)

where Vw is the molar volume of water. This can be directly compared to the osmotic pressure data obtained from the method of Luo and Roux. Given the chemical potential of water in the pure water reference state μ⊖ w , the chemical potential of water in the solute−solvent system μw can be calculated using osmotic pressure obtained from the method of Luo and Roux.

3. METHODS 3.1. Molecular Dynamics. All computations were performed in a modified version of GROMACS 5.1.,26 which changes the wall force behavior such that all solute particles experience the same wall force irrespective of the particle type. This also corrects a bug (GROMACS Bug No. 1912) which gave incorrect results when using tabulated potentials. In the original implementation, the wall force is scaled according to the atom Lennard-Jones parameters so that effective force on the solutes is not obvious. The group cutoff-scheme had to be used, as otherwise a floating point exception is encountered (GROMACS Bug No. 1843). Note that here TPI calculations are not affected by GROMACS Bug No. 1895, ignoring PME Mesh contributions in TPI, since the TPI systems do not contain charges. All covalent bonds involving hydrogen atoms were constrained with LINCS,27 and all water models use Settle28 constraints. Electrostatic interactions were treated with the Particle Mesh Ewald method29 and using a 1.4 nm cutoff. For the osmotic pressure setup, the EW3DC30 method was used for long-range electrostatic interactions. This adds a correction term to the standard Ewald summation for systems with a slab geometry that have a finite length in one dimension. A 1.4 nm cutoff was also used for the Lennard-Jones interactions. Except for the free energy calculations, we used the LeapFrog integrator31 with a 2 fs time-step size in all simulations, and the temperature of 298.15 K was controlled using the Nosé−Hoover thermostat.32,33 For the free energy calculations, we used Langevin dynamics with the stochastic dynamics integrator34 and a thermostat with 1 fs time step size. The GLYCAM06h topology and coordinate files were generated in Amber 1235 using pdb files from the GLYCAM builder36 and were converted to the GROMACS format using

2. THEORY First, we review the definition of the chemical potential and its relationship to activity- and osmotic coefficients as well as the osmotic pressure. In a system with Nw water and Np solute molecules, the chemical potential of water at constant temperature T and pressure P is ⎛ ∂G ⎞ μw = ⎜ ⎟ ⎝ ∂Nw ⎠T , P , N

(2)

(1) B


Article


Figure 1. (a) Osmotic pressure as a function of the wall force for a simulation time of 200 ns for 6.17 m β-D-glucose solutions. (b) Mass density profile dependence on the wall force in kJmol−1nm−1. In these simulations, the GLYCAM06TIP5P OSMOr14 FF (see section 4.2) was used since the largest errors are expected with that FF because of the larger osmotic pressure values.

Figure 2. Measured osmotic pressure as a function of simulation time with 32 kJmol−1nm−1 wall force. Note that in these simulations the GLYCAM06TIP5P OSMOr14 FF (see section 4.2) was used.

the glycam2gmx script, 37,38 followed by hydration in GROMACS with TIP5P water. For the GLYCAM06TIP5P OSMOr14 simulations, the parameters were altered in the GROMACS topology. Energy minimization and equilibration were performed using standard protocols. For visualization, we used Visual Molecular Dynamics (VMD).39 The pressure in the direct BAR simulations was set to 1 bar with the Parinello-Rahmann barostat.40,41 For the noble water method, we used constant NVT simulation for TPI and BAR using the average volume obtained from a previous NPT simulation using the Parinello-Rahmann barostat. For the osmotic pressure method, we used weak-pressure coupling.42 Dispersion correction was used for energy and pressure, where applicable. 3.2. Osmotic Pressure Method. To measure the osmotic pressure, we use a method similar to that developed by Luo and Roux. However, in contrast to the original method18 and the setup used in several previous studies,19−22 here periodic boundary conditions are only used in the xy-direction, as the GROMACS implementation of the wall force does not allow periodic boundary conditions. In the z-direction, the water atoms within 1 nm of the nearest wall experience a force of 10(1 − d)12 kJmol−1nm−1 to keep them contained in the box. Here, d is the distance to the nearest wall.

Virtual walls exert a constant wall force on the atoms of the solute molecules if they are outside the central region of the box. The potential shape of the wall is thus different from the previously used harmonic potential.19−23 We have tested several potential shapes for the wall forces, including the halfharmonic potential used in previous studies; however, significantly higher force constants were required to retain the solute molecules within the central region so that here the constant force version was chosen. The force exerted by each wall on all solute atoms is summed up and divided by the wallarea. The osmotic pressure is obtained as an average over both walls and all time frames, which are taken every 2 ps. Half of the difference of the osmotic pressure measured separately for both virtual walls gives an error estimate. A water bath of 2 nm outside both virtual walls was used. Test simulations using a 6 nm water bath on each side yielded identical results within the error estimate, as shown in the Supporting Information (SI) so that the smaller water bath may be used. To apply the method, which was originally devised for ionic solutions, to polysaccharides, we ran an analysis of 1 and 6.17 m β-D-glucose solutions and investigated wall forces between 1 kJmol−1nm−1 and 64 kJmol−1nm−1 in terms of the measured osmotic pressure and the mass density. For 6.17 m β-D-glucose solutions, Figure 1a shows that, starting from 32 kJmol−1nm−1, convergence of the osmotic pressure was observed. Figure 1b C


Article


Figure 3. Results from the osmotic pressure method using the unmodified GLYCAM06 TIP5P FF in comparison to experimental data.46−48 The results are shown in (a) for β-D-glucose solutions for various concentrations and in (b) for β-D-glucose (DP1), maltose (DP2), and maltotriose (DP3) solutions for a concentration of 1 m. In addition to the result using the standard protocol, the results using a reduced cutoff of 0.9 nm for all interactions is shown. Note that here we used a simulation time of 1 μs for the simulation with a 1.4 nm cutoff and only 300 ns for the simulations using the 0.9 nm cutoff, hence the larger error estimate.

shows convergence of the mass density at about 16 kJmol−1nm−1 and thus 32 kJmol−1nm−1 was used in all subsequent simulations as the wall force. The bulk density in the pure water region of ≈980 kg/m3 is consistent with the bulk density of the TIP5P water model found in a pure water simulation. In addition, the convergence analysis of the osmotic pressure over time shown in Figure 2 reveals that convergence was achieved within 400 ns for 6.17 m β-D-glucose solutions. The result does not significantly change up to a simulation time of 1 μs. Even for only 100 ns, the result is already very close to the final value, though with a significantly larger error estimate. It should be noted that the error estimate is expected to converge more slowly compared to the osmotic pressure. The slow convergence is related to the slow dynamics of the carbohydrate solutions. For a concentration of 1 m, convergence is reached significantly faster. All other solute molecules are studied at a lower concentration (mass and number based), and thus, converged results can be assumed. Nonetheless, an analysis of the osmotic pressure as a function of simulation time was routinely performed in the same way as in Figure 2, and all figures indicated convergence. In addition, for the two highest concentrations of each molecule type, longer simulation times of up to 1.5 μs were used to achieve a smaller error estimate. 3.3. Chemical Potential Calculations. The excess chemical potential μex w is calculated with the “noble water” method introduced by Knopp et. al.16 as well as by directly applying BAR to insert a molecule into a given system. In the first case, noble water, for each water model, is defined by deactivation of the electrostatic interactions for hydrogen and oxygen and, where applicable, the deactivation of the Lennard-Jones interactions for hydrogen. In addition, the oxygen σ parameter is scaled by 0.8. The free energy difference ΔAex BAR for the mutation from regular water to noble water is calculated using the BAR method (see SI for details). The phase space overlap, measured in terms of relative entropy,43 for this mutation is better than that for the inverse direction. Then, TPI is performed on the noble water system, using a 50 ns trajectory with a frame interval of 2 ps and performed 2 million insertions per frame to obtain μex TPI. The total excess chemical potential μex of water is then obtained from w μwex = −

ex ΔABAR ex − μTPI Nw

An error estimate for the BAR contribution is obtained by block averaging. For the TPI method, using twice the amount of insertions and a doubled simulation time gave virtually indistinguishable results, i.e., a deviation smaller than 0.001 kJ mol−1. Instead of this combination of methods, a single water molecule can be directly inserted into the system and the corresponding free energy calculated using the BAR method. For the BAR simulations, the system was sampled for 5 ns, and the first 500 ps were neglected as an additional equilibration. The free energy output was written every 0.1 ps. To ensure overlap in phase space, the transition from the starting state to the end state was done by introducing 25 intermediate steps defined by the coupling parameter λ ∈ [0,1]. We insert the molecules into the system as the phase space overlap, measured in terms of relative entropy,43 is slightly better. Error estimates are obtained from block averaging using five blocks.

4. RESULTS AND DISCUSSION 4.1. Osmotic Pressure Results. The results for the osmotic pressure of β-D-glucose solutions as a function of concentration in the range of 1 to 6.17 m are shown in Figure 3a for the GLYCAM06 TIP5P FF in comparison to experimental data. Despite the high complexity of the molecule, the correct trend for the osmotic pressure versus the concentration is found. However, the osmotic coefficient at 1 m is 0.67, which is significantly lower than the experimental finding, and the simulated osmotic pressures are by approximately a factor of 2 lower than the experimental values, suggesting that as discussed for the GLYCAM06 TIP3P FF and the CHARMM36 FF,23 the solute solute interactions are overly attractive. Furthermore, the deviation from linearity at higher concentration is significantly more distinct than in the experimental data. We now evaluate the osmotic pressure for maltose dimers and maltotriose trimers. In those simulations, the osmotic pressure of maltose at the same molality, was lower than that of glucose monomers over a range of concentrations. The same relationship was found between maltose and maltotriose. This trend contradicts experimental data where the opposite is observed, as summarized in Figure 3b for a concentration of 1 m. To investigate the FF dependence, the osmotic pressure for the GROMOS56A6CARBO44 SPC and the CHARMM368,45 TIP3P FF was calculated for a concentration of 4 m. The

(5) D


Article


lowered, which can be achieved by reducing the van der Waals parameters. The attraction between the saccarides can be controlled by the ε parameters, and we limit our modifications to these parameters since the σ parameters are expected to be less molecule specific. As shown before,5 the saccharide−water interactions are well described; therefore, we only modify the parameters for the nonbonded saccharide−saccharide interactions, i.e., for the nonbonded saccharide−water interactions all parameters remain unchanged. Elcock et al.23 have shown, that similar modifications can be used to tune and obtain the correct osmotic coefficient for a solution. Here, we employ an empirical parametrization process to reproduce the osmotic pressure of β-D-glucose and maltose solutions, aiming to reproduce the entire concentration range of the experimental data. To be readily usable with standard MD packages, parameter modifications are applied to all solute−solute interactions beyond the 1−4 intramolecular nonbonded interactions. 4.2.2. Parametrization. A set of new nonbonded saccharide−saccharide parameters was obtained by parametrizing for β-D-glucose and maltose solutions over a range of different concentrations. The 1−4 nonbonded intramolecular interactions were not altered. A uniform scaling of all ε parameters leads to a reduced aggregation of β-D-glucose as well as for maltose and can achieve agreement with experimental data for low concentrations. However, the uniform scaling leads to a significantly overpronounced deviation from linearity at high concentrations. In addition, although the differences between different DPs showed the correct trend, they were significantly overestimated. Therefore, the Os and Ho parameters were scaled separately. The location of Os further on the outside of the molecule makes it more relevant for direct solute−solute interactions at intermediate concentrations as compared to Cg. By varying the relative magnitude of the Os and Cg ε parameters, the concentration dependence of the osmotic pressure can be precisely tuned. To address the over-representation of the differences for different DP, consider that when comparing two β-D-glucose monomers to a maltose dimer, in terms of the number of atoms, the dimer misses one hydrogen and a hydroxyl group. In addition, intramolecular hydrogen bonds are formed in the dimer which limit the accessibility of the atoms involved for intermolecular interactions. We found that the DP dependence could be adjusted by a separate scaling of the Ho ε parameter. The new parameters are shown in Table 1. Most ε parameters of the original GLYCAM06 FF needed only a minor scaling by 0.94 to match the osmotic pressure at low concentrations. The Os ε parameter was scaled by 0.71 and the Ho ε parameter by 4.0 to obtain the correct concentration and

calculated osmotic pressure for these FFs differs even more from experimental data; it is lower by about 15% for the GROMOS56A6CARBO SPC FF than with the GLYCAM06 TIP5P FF and lower by about 25% for the CHARMM36 TIP3P FF. GLYCAM06 TIP3P was found by Elcock et al.23 to be even further from the experimental osmotic coefficient, in agreement with the nonphysical aggregation of the molecules reported previously.5 In addition, a similar decrease of the osmotic pressure with DP is also found for the CHARMM36 TIP3P FF. Since different cut-offs are frequently used, the results for the GLYCAM06 TIP5P FF with a reduced cutoff of 0.9 nm for all interactions are also shown in Figure 3. It can be seen that the cutoff has a profound impact on the osmotic pressure. However, clearly, it does not resolve the incorrect trend of the osmotic pressure as a function of the DP. Furthermore, the modification of the cutoff is likely to also affect the good performance of the GLYCAM06 TIP5P FF, when it comes to the prediction of the solution properties,5 and it does not resolve the reported aggregation problems with the GLYCAM06 TIP3P FF. Therefore, we continue to use the larger cutoff value of 1.4 nm in the following. Nonetheless, its strong effect on the osmotic pressure underlines the large effect of the choice of cutoff parameter. 4.2. Force Field Parametrization. In order to gain insights into the subtle effects of the molecular architecture on the osmotic pressure, a FF which can correctly reproduce the available experimental data is required. The problems described in ref 23 and in the last section indicate that a reparametrization of the FF parameters will be required for a predictive model. Because of its good overall performance and the closest agreement with experimental data for the osmotic data, here we focus on the GLYCAM06 TIP5P FF for this optimization. 4.2.1. Considerations. Previously,5 we found that the GLYCAM06 TIP5P FF is best suited for simulating the solution properties of polysaccharides, while with the GLYCAM06 TIP3P FF nonphysical aggregation was observed. However, although no nonphysical aggregation at low concentration with the GLYCAM06 TIP5P FF is observed, the osmotic data suggest that the FF does not describe the quantitative aggregation behavior correctly. For the GLYCAM06 FF, charges were obtained specifically for saccharides, however, based only on intramolecular structural considerations. The van der Waals parameters used in the GLYCAM06 FF, however, are partly taken from the PARM9449 FF which used densities and enthalpies of vaporization of small molecules for the parametrization of the C and H atoms. O and Oh were taken from the OPLS FF,50 which was also parametrized for densities and solution properties. Thus, the nonbonded saccharide−saccharide interactions have not been explicitly included in the parametrization process of the GLYCAM06 TIP5P FF. The same holds for the CHARMM36 FF. However, the results for the osmotic pressure illustrate that it is important to include molecule specific data into the FF parametrization, to properly balance the solute−solute and solute−solvent interactions. In addition, the aggregation is primarily driven by the van der Waals interactions and counteracted by electrostatic interactions.5 Therefore, to parametrize for the correct osmotic trends, we focus on the van der Waals parameters. Since the original GLYCAM06 TIP5P FF leads to an osmotic pressure that is systematically too low, the aggregation needs to be

Table 1. New Nonbonded ε Parameters for the Saccharide− a Saccharide Interactions of the GLYCAM06TIP5P OSMOr14 FF

a

E

atom type

ε (kJ mol−1)

H1 Ho Os H2 Oh Cg

0.061747472 0.50208 0.5050088 0.061747472 0.62502294 0.4302662

The 1−4 nonbonded intra-molecular interactions were not altered. DOI: 10.1021/acs.jctc.6b00295 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article


Figure 4. Calculated osmotic pressure of the model (top) and evaluation molecules (bottom) using the GLYCAM06TIP5P OSMOr14 FF in comparison to experimental data.46−48

Figure 5. Free energy maps for β-D-glc(1−4)β-D-glc with both the original GLYCAM06 TIP5P FF and the GLYCAM06TIP5P OSMOr14 FF.

because for each parameter set, the osmotic pressure of the whole concentration range has to be calculated. Another possibility to explore, to further reduce the changes of individual parameters, is to apply the introduced parametrization procedure homogeneously to both solute−solute and solute−solvent interactions. A strategy that might reduce the computational cost for such a balancing act is presented by the Kirkwood-Buff integrals, which relate the thermodynamic properties of a solution to the RDFs.22,51 4.3. Osmotic Pressure Results for GLYCAM06TIP5P OSMOr14. The osmotic pressure results of β-D-glucose and maltose solutions for the new GLYCAM06TIP5P OSMOr14 FF at different concentrations are shown in Figure 4. In addition, the osmotic pressure results for β-D-xylose and maltotriose solutions are shown to evaluate the performance of the new parameters with

DP dependence, respectively. As before, all mixed parameters are obtained by the standard combination rule. As discussed in ref 23, changing the Os ε parameters can affect the strength of hydrogen bonds. The same can be expected for the changed Ho parameters. As differences in the molecular architecture or the type of linkage may change the probabilities for solute−solute hydrogen bond formation, more in-depth testing is required to evaluate the transferability of the new FF. The much more highly occupied solute−solvent hydrogen bonds, however, should not be affected since only solute−solute interactions were modified. It should be noted that the effect of these parameter modification are interdependent, making the balancing procedure highly challenging. Furthermore, this parametrization procedure is demanding in terms of computational cost F


Article


Figure 6. (a) Mass density of β-D-glucose solutions at different concentrations. (b) Diffusion coefficients of β-D-glucose solutions at different concentrations.

Figure 7. RDFs between β-D-glucose monomers (COM) in water at different concentrations with both the original GLYCAM06 TIP5P FF and the GLYCAM06TIP5P OSMOr14 FF.

properties. As solute−solvent interactions are not altered by our reparametrization, only observables that involve nonbonded saccharide−saccharide interactions were tested. We compare the mass density and diffusion coefficients and dielectric constants to those found with the unmodified FF and experimental data. In addition, we compare the aggregation behavior in terms of the saccharide−saccharide RDFs for β-Dglucose solutions between the original and the new FF. The data for mass density and diffusion coefficient are shown in Figure 6. Especially for the mass density, the GLYCAM06TIP5P OSMOr14 FF gives almost identical results to those of the GLYCAM06 TIP5P FF. Only at high saccharide concentrations does a small difference appear where the density is shifted toward lower values as the saccharide aggregation is reduced. The diffusion coefficients shown in Figure 6b reflect the reduced aggregation which results in faster diffusion. Although the absolute increase is larger at lower concentrations, where the diffusion is much faster, the relative increase is most significant for higher saccharide concentrations. For the dielectric constants, like for the optimized GLYCAM06 TIP3P and CHARMM36 parameters,23 no significant change was found with the new parameters. To get a more direct idea of the effect of the parameters on the aggregation behavior, the RDFs obtained with the TIP5P GLYCAM06OSMOr14 FF and GLYCAM06 TIP5P FF are shown in Figure 7a and b for various concentrations. The strong concentration dependence of the RDFs, observed with the GLYCAM06 TIP5P FF, is reduced with the new

molecules that were not used in the parametrization process. The results are in good agreement with experimental data in all cases. For β-D-glucose as well as maltotriose, the osmotic pressure is slightly overestimated at high concentrations, while the osmotic pressure is slightly too low for β-D-xylose because the effect of the missing atoms in β-D-xylose compared to β-Dglucose is still somewhat over-represented. Nonetheless, these effects are small, and in all cases, the results are very close to the experimental data, and all trends with respect to concentration as well as DP are correctly reproduced. 4.4. Further Force Field Evaluation. First, to evaluate the effect of the parameter modification on the internal conformations of the molecules, we compared the conformations of the most flexible degrees of freedom, the dihedrals in the glycosidic link, for a β-D-glc(1−4)β-D-glc (cellobiose) molecule. The free energy maps for the glycosidic link were calculated with metadynamics52 using the Plumed 2.2.253 implementation with GROMACS 5.1.2. (See SI for details of the simulation setup.) These maps are shown in Figure 5 for both the original GLYCAM06 TIP5P FF and the GLYCAM06TIP5P OSMOr14 FF. Figure 5 shows that the conformations of the most flexible degrees of freedom are essentially preserved. Therefore, the separation of intra- and intermolecular Lennard-Jones parameters beyond 1− 4 is not necessary. Then, to ensure that the new ε parameters do not affect the FFs general applicability we evaluate the performance of the new GLYCAM06TIP5P OSMOr14 FF for the prediction of other solution G


Article

Journal of Chemical Theory and Computation Table 2. Hydration Free Energies of Water for Different Water Models at 298.15 Ka water model 54

SPC SPC/E55 OPC56 TIP3P TIP4P757 TIP4P-EW58 TIP5P TIP5P-EW59

ex ΔABAR Nw

(kJ mol−1) 24.77 28.38 31.66 24.24 24.2 27.89 22.26 22.66

−1 μex TPI (kJ mol )

−1 μex w (kJ mol )

ΔG (kJ mol−1)

|26.46 − ΔG| (kJ mol−1)

−1 |μex w − ΔG| (kJ mol )

3.36 3.36 1.96 2.72 2.79 2.72 2.70 2.18

−28.13 −31.74 −33.62 −26.96 −26.99 −30.61 −24.97 −24.84

−26.13 ± 0.05 −29.7 ± 0.05 −33.62 ± 0.07 −25.62 ± 0.05 −25.59 ± 0.04 −29.22 ± 0.04 −23.75 ± 0.05 −24.47 ± 0.05

0.33 3.24 7.16 0.84 0.87 2.76 2.71 1.99

2.0 2.04 0.003 1.34 1.40 1.39 1.21 0.37

The experimental value24 is −26.46 kJ mol−1. The values μex w are found using eq 5 and the values listed in the second and third columns. The error estimate for the noble water method (BAR < 0.01 kJ mol−1 and TPI < 0.001 kJ mol−1) is negligble. ΔG lists the values obtained by the direct BAR method. The last two columns compare ΔG to the experimental values and the noble water results, respectively. a

difference between a 4 m β-D-glucose solution and pure water is found as −0.22 kJ mol−1. In principle, one could think about measuring μex w of the saccharide solutions directly by a similar method. In fact, the noble water method was originally devised for the calculation of the excess chemical potential of water in dense polymer structures. However, the experimental excess chemical potential difference between pure water and a 4 m β-D-glucose solution is only approximately −0.2 kJ mol−1 at 298.15 K. Thus, the magnitude of the TPI error found for pure water exceeds the precision required to observe the solute effects on the excess chemical potential. In addition, the mutation to noble water perturbs the fragile balance between solute−solute and solute− solvent interactions and causes the saccharides to aggregate. This leads to massive sampling and convergence problems, as shown in the SI. Although the direct BAR method does not include the TPI associated error, sufficient sampling of the complex system is clearly infeasible, leaving the osmotic pressure method as the only feasible approach.

GLYCAM06TIP5P OSMOr14 FF. A comparison of the RDFs, especially at 4 m, reveals that the structural patterns, best recognizable in the two main coordination peaks, are essentially preserved. However, especially the first peak is significantly less pronounced, reflecting the lower solute−solute aggregation. The lower tendency to aggregate means that the reduction of the solvent entropy will be larger, which is reflected in the higher osmotic pressure measured for these parameters. In addition, for the new FF the peaks are slightly shifted to longer distances. In contrast to the mass density and diffusion coefficients, the RDF changes most significantly for low concentrations because at high concentration, the system is more homogeneously filled with saccharides. Therefore, correlations have a much higher impact at low concentrations. 4.5. Chemical Potential. As discussed in section 2, the osmotic pressure can be related to the chemical potential and the activity of water in the solution. To make that connection, the chemical potential of water in the pure water reference ex state, μ⊖ w , is required. Here, we calculate μw in a pure water system and compare the values obtained using the noble water method or the direct BAR method. The results are summarized in Table 2 for different water models. The experimental value24 is −26.46 kJ mol−1. The results from the direct BAR method reproduce the experimental value reasonably well, and the results for SPC, TIP3P, and TIP4P water are in excellent agreement with the value obtained previously.25,60 Good agreement is also found for TIP3P, SPC, TIP4P, and TIP5P with results from a completely different method,61 whereas the values found here for SPCE and TIP4P-EW are significantly larger. Comparing the two methods, it is noticeable that the absolute values found with the noble water method are systematically larger than those from the direct BAR method. Only for the OPC model do the results from both methods agree. The magnitude of the difference depends on the water model. In general, it is largest where the TPI contribution is more significant, and thus, it appears that the discrepancy stems from the TPI method. The uniform insertions result in a bias toward closer contacts and thus more negative values of μex w for dense systems.62 It is surprising, however, that a pronounced bias is already apparent for such a low density system. The value for TIP5P water with the BAR method is μex w = −23.75 ± 0.05 kJ mol−1. Using eqs 4 and (2), μex w can be −1 for a 5 calculated, yielding, for instance, μex w = −24.03 kJ mol −1 ex m β-D-glucose solution and μw = −23.99 kJ mol for β-Dxylose at the same concentration. The chemical potential

5. CONCLUSIONS The osmotic pressure of saccharide solutions can be calculated in MD simulations using the GLYCAM06 TIP5P FF, and the results can reproduce the correct trends for the concentration dependence with monomers. However, osmotic pressure values are too low, and the experimental trend of the osmotic pressure as a function of the DP is not reproduced with either the CHARM36 TIP3P FF or the GLYCAM06 TIP5P FF. Therefore, we introduce a new set of ε van der Waals parameters for nonbonded saccharide−saccharide interactions, which were not parametrized with saccharide−saccharide specific data for either of the two FFs in the GLYCAM06 TIP5P FF. The predictive power of the resulting GLYCAM06TIP5P OSMOr14 FF was demonstrated for the osmotic pressure trend of different saccharide solutions. In addition, mass density and diffusion coefficients were slightly improved in comparison to the original GLYCAM06 TIP5P FF. Nonetheless, more indepth testing should be performed on saccharides with different linkages and more complex structures to asses the transferability of the new FF. Another path that may be worthwhile exploring is to apply the introduced parametrization procedure homogeneously to both solute−solute and solute−solvent interactions. The excess chemical potential of water for saccharide solutions can be calculated from the osmotic pressure data, using the excess chemical potential of water for pure water H


Article


(25) Jorgensen, W. L.; Blake, J. F.; Buckner, J. K. Chem. Phys. 1989, 129, 193−200. (26) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D. Bioinformatics 2013, 29, 845−854. (27) Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. J. Comput. Chem. 1997, 18, 1463−1472. (28) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952− 962. (29) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089−10092. (30) Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155− 3162. (31) Birdsall, C. K.; Langdon, A. B. Plasma Physics via Computer Simulation; McGraw-Hill: New York, 1985. (32) Nosé, S. J. Chem. Phys. 1984, 81, 511−519. (33) Hoover, W. G. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695. (34) Van Gunsteren, W.; Berendsen, H. Mol. Simul. 1988, 1, 173− 185. (35) Case, D.; Darden, T.; Cheatham, T., III; Simmerling, C.; Wang, J.; Duke, R.; Luo, R.; Walker, R.; Zhang, W.; Merz, K. Amber 12; University of California, San Francisco, CA, 2012. (36) Polysaccharides: Structural Diversity and Functional Versatility (2005−2016); Woods Group GLYCAM Web. Complex Carbohydr. Res. Center, University of Georgia: Athens, GA. http://www.glycam. com. (37) Sorin, E. J.; Pande, V. S. Biophys. J. 2005, 88, 2472−2493. (38) Wehle, M.; Vilotijevic, I.; Lipowsky, R.; Seeberger, P. H.; Varon Silva, D.; Santer, M. J. Am. Chem. Soc. 2012, 134, 18964−18972. (39) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38. (40) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182−7190. (41) Nose, S.; Klein, M. Mol. Phys. 1983, 50, 1055−1076. (42) Berendsen, H. J.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. J. Chem. Phys. 1984, 81, 3684−3690. (43) Wu, D.; Kofke, D. A. J. Chem. Phys. 2005, 123, 054103. (44) Hansen, H. S.; Hünenberger, P. H. J. Comput. Chem. 2011, 32, 998−1032. (45) Guvench, O.; Greene, S. N.; Kamath, G.; Brady, J. W.; Venable, R. M.; Pastor, R. W.; Mackerell, A. D. J. Comput. Chem. 2008, 29, 2543−2564. (46) Uedaira, H.; Uedaira, H. Bull. Chem. Soc. Jpn. 1969, 42, 2137− 2140. (47) Miyajima, K.; Sawada, M.; Nakagaki, M. Bull. Chem. Soc. Jpn. 1983, 56, 1620−1623. (48) Miyajima, K.; Sawada, M.; Nakagaki, M. Bull. Chem. Soc. Jpn. 1983, 56, 1954−1957. (49) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179−5197. (50) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657−1666. (51) Ploetz, E. A.; Bentenitis, N.; Smith, P. E. Fluid Phase Equilib. 2010, 290, 43−47. Proceedings of the 17th Symposium on Thermophysical Properties. (52) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (53) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. Comput. Phys. Commun. 2014, 185, 604−613. (54) Berendsen, H. J.; Postma, J.; Van Gunsteren, W.; Hermans, J. Intermolecular Forces; Springer: Dordrecht, Netherlands, 1981; pp 331−342. (55) Berendsen, H.; Grigera, J.; Straatsma, T. J. Phys. Chem. 1987, 91, 6269−6271. (56) Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. J. Phys. Chem. Lett. 2014, 5, 3863−3871. (57) Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985, 56, 1381− 1392.

systems, which can be calculated by alchemical free energy methods.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00295. Method details (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Funding

The project was funded by the Deutsche Forschungsgemeinschaft (GR 3661/2-1). Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors thank Reinhard Lipowsky and Ana Celia Vila Verde for helpful discussions. REFERENCES

(1) Brocker, C.; Thompson, D. C.; Vasiliou, V. Biomol. Concepts 2012, 3, 345−364. (2) Burgert, I.; Eder, M.; Gierlinger, N.; Fratzl, P. Planta 2007, 226, 981−987. (3) Harrington, M. J.; Razghandi, K.; Ditsch, F.; Guiducci, L.; Rueggeberg, M.; Dunlop, J. W. C.; Fratzl, P.; Neinhuis, C.; Burgert, I. Nat. Commun. 2011, 2, 337. (4) Dawson, C.; Vincent, J. F. V.; Rocca, A.-M. Nature 1997, 390, 668−668. (5) Sauter, J.; Grafmüller, A. J. Chem. Theory Comput. 2015, 11, 1765−1774. (6) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; GonzálezOuteiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. J. Comput. Chem. 2008, 29, 622−655. (7) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (8) Guvench, O.; Hatcher, E.; Venable, R. M.; Pastor, R. W.; MacKerell, A. D., Jr J. Chem. Theory Comput. 2009, 5, 2353−2370. (9) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910−8922. (10) Widom, B. J. Chem. Phys. 1963, 39, 2808−2812. (11) Jedlovszky, P.; Mezei, M. J. Am. Chem. Soc. 2000, 122, 5125− 5131. (12) Hörstermann, H.; Hentschke, R.; Amkreutz, M.; Hoffmann, M.; Wirts-Rütters, M. J. Phys. Chem. B 2010, 114, 17013−17024. (13) Shing, K.; Gubbins, K. Mol. Phys. 1982, 46, 1109−1128. (14) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300−313. (15) Bennett, C. H. J. Comput. Phys. 1976, 22, 245−268. (16) Knopp, B.; Suter, U. W.; Gusev, A. A. Macromolecules 1997, 30, 6107−6113. (17) Murad, S.; Powles, J. J. Chem. Phys. 1993, 99, 7271−7272. (18) Luo, Y.; Roux, B. J. Phys. Chem. Lett. 2010, 1, 183−189. (19) Yoo, J.; Aksimentiev, A. J. Phys. Chem. Lett. 2012, 3, 45−50. (20) Yoo, J.; Wilson, J.; Aksimentiev, A. Biopolymers 2016, 105, 752− 763. (21) Yoo, J.; Aksimentiev, A. J. Chem. Theory Comput. 2016, 12, 430−443. (22) Karunaweera, S.; Gee, M. B.; Weerasinghe, S.; Smith, P. E. J. Chem. Theory Comput. 2012, 8, 3493−3503. (23) Lay, W. K.; Miller, M. S.; Elcock, A. H. J. Chem. Theory Comput. 2016, 12, 1401−1407. (24) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016−2027. I


Article

Journal of Chemical Theory and Computation (58) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665− 9678. (59) Rick, S. W. J. Chem. Phys. 2004, 120, 6085−6093. (60) Hermans, J.; Pathiaseril, A.; Anderson, A. J. Am. Chem. Soc. 1988, 110, 5982−5986. (61) Henchman, R. H. J. Chem. Phys. 2007, 126, 064504. (62) Berendsen, H. J. Simulating the Physical World: Hierarchical Modeling from Quantum Mechanics to Fluid Dynamics; Cambridge University Press: Cambridge, UK, 2007.

■

NOTE ADDED AFTER ASAP PUBLICATION This paper was published ASAP on August 29, 2016 with errors in Table 2. The corrected version was reposted on August 29, 2016.

J


Predicting the Chemical Potential and Osmotic Pressure of

Recommend Documents