Structure of Aqueous Solutions of Monosodium Glutamate - American

Apr 27, 2009 - Christopher D. Daub,† Kevin Leung,‡ and Alenka Luzar*,†. Department of Chemistry, Virginia Commonwealth UniVersity, Richmond, Vir...
0 downloads 0 Views 5MB Size
Download PDF
J. Phys. Chem. B 2009, 113, 7687–7700

7687

Structure of Aqueous Solutions of Monosodium Glutamate Christopher D. Daub,† Kevin Leung,‡ and Alenka Luzar*,† Department of Chemistry, Virginia Commonwealth UniVersity, Richmond, Virginia 23284, and Sandia National Laboratories, MS 1415, Albuquerque, New Mexico 87185 ReceiVed: NoVember 26, 2008; ReVised Manuscript ReceiVed: February 13, 2009

We studied monosodium glutamate (MSG) in aqueous solution using molecular dynamics simulations and compared the results with recent neutron diffraction with isotope contrast variation/empirical potential structure refinement (EPSR) data obtained on the same system (McLain et al. J. Phys. Chem. B 2006, 110, 21251-21258). We used classical simulations with empirical force fields to study both dilute and more concentrated (1.40 M) solutions. To gauge the importance of polarization and other quantum effects, we carried out first-principles molecular dynamics in the dilute case. The glutamate structure was well reproduced by the OPLS/AA and SPC/E force fields: we found a reasonable agreement between the simulations and the experimental data with respect to the hydration numbers for glutamate carboxylate and amine groups and the observation of significant sodium ion-carboxylate binding. However, none of our simulations could reproduce the dramatic reduction in water-water correlations observed experimentally. Simulations showed a large amount of carboxylate-amine binding, as well as segregation of water and glutamate, at moderately high concentrations of MSG. We attribute this result to the breakdown of currently available classical force fields when applied to concentrated ionic solutions, especially large polyatomic ions. We also did not observe the sharing of a water proton by two carboxylate oxygens simultaneously, and we argue against this interpretation of EPSR data on a variety of physical grounds. We offer several suggestions to resolve these discrepancies between simulation and the current interpretation of neutron diffraction data, which should advance the understanding of aqueous ionic solutions in general. I. Introduction Understanding the structure of aqueous solutions and hydrogenbonded mixtures has long been an active area of research. Among the most useful experimental techniques for studying solution structure in such systems is neutron diffraction with isotope substitution (NDIS). Through substitution of isotopes in the solute and/or the solvent, usually replacement of hydrogen with deuterium, data for different isotopic makeups can be analyzed to extract correlations between different individual pairs of atoms. Traditionally, the isotope method has been used to examine water structure around small and spherical ions in solution,1,2 in which case the symmetry of the problem allows the local water structure to be derived directly from the radial distribution functions. However, for molecules in solution, the site-site radial distribution functions are averaged over molecular orientations and so, by themselves, cannot be interpreted in the same way as for ions. Furthermore, because of limitations on the available isotopic contrast, it is normally not possible to obtain a complete separation of all N(N + 1)/2 site-site correlations that are needed to define the structure of a system containing N distinct atomic components. To obtain a full set of correlations for a system in question, Soper developed a modeling technique known as empirical potential structure refinement (EPSR).3 It combines accurate structural data from neutron diffraction with the use of Monte Carlo simulation to generate atomic configurations that are consistent with the experimental data. In this way, pair correlations between all of the atom types in the system can be computed using the same * To whom correspondence should be addressed. E-mail: [email protected]. † Virginia Commonwealth University. ‡ Sandia National Laboratories.

statistical techniques as are used routinely in other computer simulations. Because the generated configurations are compared with several isotopically substituted neutron diffraction results, a high degree of confidence can be attached to the structural data so generated. Neutron diffraction in tandem with EPSR has been used with great success to describe simple fluids such as neat water,3 amphiphilic molecules such as acetone4 and dimethyl sulfoxide,4,5 and liquid mixtures.6,7 More recently, practitioners of the method have extended their reach to aqueous ionic solutions, mainly of hydroxides8-10 and alkali halides.11,12 These studies have led to new insights into the impact of ions on the network of hydrogen-bonded solvating water. For example, studies of aqueous NaOH8,10 have revealed that, even at a modest concentration of 2 M, the water-water radial distribution functions (RDFs) are quite strongly affected by the presence of the ions, with the reduction in the first peak height in gOWOW(r) corresponding to a significant reduction in the water hydration number (i.e., the number of water molecules in the solvation shell of a particular water molecule). Studying the hydration of amino acids is an important step toward understanding the hydration of proteins. It has long been known that the amount of water hydrating a protein can be determined by simply adding up the water hydrating its constituent amino acids.13 Amino acids show a wide variety of solvation behaviors in water, from near-insolubility for Ltyrosine to solubility in excess of 200 g/L for L-serine and glycine.14 The relative solubility of amino acids depends on the degree of solute-solute association; when 1:1 amino acid/water mixtures are frozen and analyzed, the water in most of the samples has an icelike structure, whereas in highly soluble

10.1021/jp810379m CCC: $40.75  2009 American Chemical Society Published on Web 04/27/2009

7688

J. Phys. Chem. B, Vol. 113, No. 21, 2009

L-serine, the water structure does not resemble ice, showing that the two species are well-mixed.15 To understand the structure of solutions, such as amino acids in aqueous solutions, ideally, one needs to consider three aspects: first, how the solute molecules, or important chemical groups on those molecules, hydrate (solute-solvent correlations); second, how the structure of the hydration water is perturbed from that in the bulk solvent (solvent-solvent correlations); and third, what the structure is from the solute molecule’s point of view (solute-solute correlations). Neutron diffraction studies have now been performed on aqueous solutions of L-proline,16 zwitterionic glutamate ions and sodium counterions [monosodium glutamate (MSG)],17,18 and zwitterionic dipeptides.19,20 In both the proline study and the dipeptide study, the overall neutral solutes do not require the addition of counterions to neutralize the system. This difference has an impact on all aspects of the solution structure. In proline, as the concentration is increased, the degree of water-water coordination is slightly reduced but remains very close to that of bulk water.16 Neutron diffraction/ EPSR studies of zwitterionic dipeptides in aqueous solution19 have shown some reduction in the water-water coordination numbers, not to the extent seen in ionic solutions, but to a greater extent than in the proline study. Studies have also detected the formation of clusters of dipeptides,20 in contrast to the proline study, where only slight clustering was observed.16 In the aqueous MSG study, the solutes were found to have a much greater effect on the solvent structure, even at a lower solute concentration, showing a reduction by nearly one-half of a water-water hydrogen bond per water molecule on average.17 The presence of the charge on the glutamate also increases the number of solvating water molecules compared to proline. The number of water molecules solvating each amine nitrogen rises from 1.7 in aqueous proline16 to 4.3 in aqueous MSG.17 The corresponding number for each carboxylate oxygen rises from 2.8 in aqueous proline16 to 4.5 in aqueous MSG.17 This large number of solvating water molecules in MSG is interpreted as indicating the sharing of a solvating water molecule by both carboxylate oxygens at once.17 The clustering of dipeptides also suggests a lower number of solvating water molecules in these cases. No significant clustering was observed in aqueous MSG.17 In addition to experimental data from neutron diffraction studies, molecular simulations are a common technique for studying aqueous solutions, providing direct insight into the microscopic physics underlying the experimental measurements. Comparison of computer simulations with neutron diffraction data on amino acids in aqueous solutions can also provide a benchmark with which to improve force fields that are widely used in molecular dynamics software packages by academics and the pharmaceutical industry. A recent study21 investigated a concentrated aqueous solution of proline via classical molecular dynamics (MD) with the CHARMM2222 force field in conjunction with several different empirical water models23,24 and with the first principle MD. Both simulation data were compared with the neutron diffraction experiment on the same system16 that we mentioned above. The match between the simulation and the experiment was generally good, and the agreement improved when approximate quantum corrections were applied. Given the difficulty in accurately simulating polyatomic ions in solution with classical methods,25 the natural question to ask is whether classical force fields can be reliable for modeling a charged amino acid in aqueous solution. Glutamic acid (C5H9NO4) is a particularily interesting amino acid; it has an unusually low solubility,26 which would not be expected a priori

Daub et al. because of its negatively charged side chain. Glutamic acid is the most abundant amino acid in the human diet,27 making up 20-40% by weight of typical plant proteins, and mammalian taste receptors associated with the “umami” taste respond strongly to MSG,28,29 demonstrating the importance of this molecule in human biochemistry. If concentrated amino acid solutions are to be simulated accurately, it is crucial that the degree of association between solutes be correctly reproduced. The charge on glutamate and the presence of ions provides a stringent test of classical force fields for simulation of aqueous amino acids and proteins. In this article, we report an extensive MD simulation study of MSG solutions and a comparison of the results with the aforementioned recent high-resolution neutron diffraction/EPSR study of the same system.17,18 We examined structural data including RDFs and their logical extension, three-dimensional spatial distribution functions (SDFs), which yield details of the spatial relationships among water molecules around the glutamate ions, as well as several distributions of dihedral angles that turned out to be particularly force-field-dependent. We considered the possible role of electronic polarizability effects by comparing simulations of a dilute (one Glu- ion) MSG solution made using both classical and ab initio molecular dynamics (AIMD), and we compared our results with previous studies of carboxylate and/or amine hydration of amino acids and other molecules, both theoretical21,25,30-40 and experimental.16,19,41-43 We also increased the concentration of the solution to ∼1.4 M to investigate the effect of increased concentration on the hydration structure by classical simulations only. Although recent simulations of concentrated NaOH44 and LiCl45 solutions point toward the possibility, concentration effects that require large simulation cells containing multiple polyatomic solutes cannot yet be directly studied with ab initio methods. AIMD data can also be used to accurately describe the intramolecular correlations in individual water molecules, which allowed us to directly compare our results with the experimental structure factors F(Q). Changing the solute force field greatly affects the molecular conformation of glutamate, with the OPLS/AA46-48 force field being closer to the neutron diffraction results.18 Changing the concentration has appreciable effects on the entire solution structure. We found classical-force-field simulations, AIMD, and force-field-assisted empirical potential structure refinement (EPSR)3 of neutron diffraction data17 to be in qualitative agreement in many respects but differing in others. Our simulation results highlight possible shortcomings of both EPSRpredicted solution structures and current classical force fields for amino acids and ions. In particular, our results disagree on the observation of water sharing by carboxylate oxygens and on the weakening of water-water correlations reported in ref 17. The present work emphasizes the importance of all three types of studies and suggests that detailed comparison among them might pave the way toward a full, quantitative understanding of the hydration structures of aqueous amino acid solutions in particular and ionic solutions in general. This study of solution structure should also be relevant to an improved understanding of protein hydration, mirroring recent work on the range of solute-driven perturbations of water dynamics, a problem that is still hotly debated.49,50 II. Theoretical Methods Potential Models. Water. We used the SPC/E model23 to describe the interactions between water molecules. We chose this model among the many available force fields for water

Structure of Aqueous Solutions of MSG

J. Phys. Chem. B, Vol. 113, No. 21, 2009 7689

TABLE 1: Summary of Simulation Conditions in This Study and in Ref 17a system f

1 2g 3g EPSRh

nH2Ob

nGlu-b

nNa+b

Lc (Å)

concentrationd (M)

atomic densitye (Å-3)

95 248 216 580

1 1 6 20

0 1 6 20

14.13 19.486 19.215 27.68

0.59 0.22 1.40 1.57

0.107 0.103 0.107 0.100

a See text for other details. b na is the number of molecules of species a in the simulation box. c L is the linear size of the cubic simulation box. d Concentration is given in moles of glutamate per liter of solution. e For reference, the experimental value of the atomic density is 0.1043 atoms/Å3 at a concentration of 1.37 M.55 f Conditions for system 1 are the same as for the AIMD simulations. g Conditions for systems 2 and 3 are classical MD with the OPLS/AA force field46-48 for the glutamate(s) and the Smith and Dang force field53 for the counterion(s). h Reference 17.

because it reproduced the experimental diffusion constant of water. Although an accurate diffusion constant is not crucial for structural studies, we are also exploring dynamic properties in aqueous amino acid solutions,51 so for consistency, we used SPC/E for the current study as well. Another consideration is that EPSR simulations generally use SPC/E as their reference potential for water.17 Glutamate. For interactions within and between glutamate ions, we applied the OPLS/AA all-atom force field developed by Jorgensen and co-workers.46-48 Even though OPLS/AA is parametrized for use with the TIP3P water model, it has generally been found to work well with the SPC/E model as well. In a recent study of aqueous proline, it was found that SPC/E gave good results when used in combination with another biomolecular force field (CHARM22) that was parametrized with TIP3P water.21 For comparison, we also performed some runs using the Amber ff94 force field.52 We studied the zwitterionic glutamate ion, with terminal NH3+ and COO- groups, as it is the equilibrium structure in aqueous solution. Neither the OPLS/AA force field nor the Amber ff94 force field directly simulates a zwitterionic amino acid, so some modifications were necessary. Both force fields provide parameter sets for amino acids present as C- and N-termini in a polypeptide, as well as for an amino acid in the middle of a peptide chain. In effect, a zwitterionic amino acid is both a Cand N-terminus at the same time, and this can guide the necessary force field modifications. The most recent OPLS/AA force field46-48 was used with little modification. The charge on CR was set to the average of the values prescribed for the carbons in the R2CHCOO- and R2CHNH3+ groups (-0.16e and 0.25e, respectively), plus an extra 0.1e to maintain the overall charge on the glutamate at -1e. The resultant charge on CR is virtually the same as in the force field for the capped amino acid (0.15e versus 0.14e). In the Amber ff94 force field, the partial charges on each atom differ between the C- and N-termini. The charges were therefore set to the average of the C- and N-terminal charges in all cases where the atoms were common to both, and then all charges were scaled to maintain the total system charge at -1e. The total charge before scaling was within 5% of -1e, so this scaling should not have caused a noticeable effect. There are no specific dihedrals prescribed for individual amino acids in the Amber force field, making the construction of the dihedral potentials for the zwitterionic glutamate ion straightforward. Sodium Ions. For the sodium counterions, we used potential parameters due to Smith and Dang.53 These parameters were developed to be used with the RPOL polarizable water model, but they can also be used with little modification with SPC/E water.53 There are, however, a variety of ionic force field parameters available. We have found (results not shown) that the alternate parametrization due to Aqvist54 yields qualitatively similar results.

It is crucial to apply the appropriate mixing rules for unlike Lennard-Jones centers when simulating with different force fields, as these are often parametrized with differing mixing rules. Water-glutamate, glutamate-sodium, and water-sodium Lennard-Jones interaction parameters were obtained with the geometric mixing rules used throughout the OPLS/AA force field [σij ) (σiσj)1/2, εij ) (εiεj)1/2]. For the Amber force field, we used the Lorenz-Berthelot combining rules. Simulation Methodology. Classical Simulations. Initial configurations for the low-concentration classical simulations were prepared by inserting a glutamate zwitterion molecular configuration prepared by using LeAP in Amber 7 into an equilibrated configuration of 256 water molecules, along with a Na+ counterion. After the water molecules overlapping with the glutamic acid and counterion had been removed, 248 water molecules remained. We refer to this system as system 2 (see Table 1 for a summary of the systems under study). We prepared the initial configuration for the high concentration by inserting five more copies of the equilibrated glutamate configuration (and five more counterions) into the simulation box and again deleting overlapping waters. This gave us a system of 6 Glu- ions, 6 Na+ ions, and 216 water molecules, to which we refer as system 3. The concentration of this Glu- solution (∼1.4 M) is comparable to the experimental concentration of 2 M and the concentration of the EPSR simulation of ∼1.6 M.17 The linear size of the cubic simulation box varied with the force field(s), but was between 19 and 20 Å in all cases. The integration step size used was 0.5 fs. Ewald summation was used to compute the long-range electrostatic interactions. The SHAKE algorithm was used to maintain the rigid geometry of the water molecules. All simulations were executed at a target temperature of 300 K. The box size was adjusted by running in the NPT ensemble with a Nose-Hoover barostat at 1 atm for 100 ps. The resultant system densities for each simulation are given in Table 1, expressed in atoms/Å3 to aid comparison with neutron diffraction experiments. We also give the experimental density of MSG solutions55 at a comparable concentration. The simulation densities are somewhat higher than the experimental data, most likely because of slight inaccuracies in the force fields that we used. Constant pressure-simulations were followed by an equilibration run of 4 ns in the NVT ensemble before data collection was begun. We completed at least two 4-ns classical NVT simulations of each of the systems that we studied. When the force field was changed, we re-equilibrated the system by following the above steps starting from the final configuration of a simulation using the original force field. The very longlived nature of the ion-glutamate and glutamate-glutamate coordinations make equilibrations of at least 1 ns necessary, to ensure that all simulations can be considered independent and uncorrelated.

7690

J. Phys. Chem. B, Vol. 113, No. 21, 2009

Daub et al.

AIMD Simulations. We conducted two 20-ps AIMD trajectories of an isolated glutamate anion solvated in 95 water molecules starting from different initial configurations generated with OPLS/AA and SPC/E force fields. They were performed using the revised Perdew-Burke-Ernzerhof (RPBE) exchangecorrelation functional,56 the Vienna Ab-initio Simulation Package (VASP),57,58 the projector augmented wave method,59 and associated pseudopotentials.60 PBE pseudopotentials were applied in the RPBE calculations. The deuterium mass was substituted for all protons, and a 0.25-fs time step was used in the Born-Oppenheimer dynamics AIMD simulations. A 400 eV energy cutoff and a 10-6 eV accuracy convergence criterion were enforced at each time step. These settings led to energy drifts of less than 1 K/ps. A thermostat kept the systems at T ) 300 K; under these conditions, the RPBE approach predicts a reasonable liquid water structure.61 The initial configurations were equilibrated for ∼1 ps before AIMD statistics were collected. We refer to this system as system 1. Radial and Spatial Distribution Functions. Many useful functions for understanding the structure of amorphous systems such as ours can be computed from simulation data.62 One of the most well-known is the radial distribution function (RDF), also denoted gRβ(r). The RDF can also be integrated to find the coordination number nRβ(r), which defines the average number of atoms of type β found out to a distance r around each atom of type R. If type β is the water oxygen atom, this is also referred to as the hydration number. The radial distribution function, however, has some disadvantages. It is dependent only on the distance r, so that all orientational information is averaged out. To include orientational information, several researchers63-66 have proposed the extension of the concept of the radial distribution function to three dimensions by defining the spatial distribution function (SDF) as

〈 g(x, y, z) )

∑∑ i

Nbin

molecule, i, forming the apex of the angle must be in the solvation shell of (a) given solute atom(s), and the other two water molecules, j and k, must both be near neighbors of the central water. Mathematically ni

δ(x - xij) δ(y - yij) δ(z - zij)〉

j

Figure 1. Distributions of side-chain dihedral angles at both low MSG concentrations [top, systems 1 (AIMD, dotted lines) and 2 (OPLS/ AA, solid lines)] and high MSG concentrations [bottom, system 3 (OPLS/AA, solid lines) and experiment/EPSR (ref 18, dotted lines)]. (See Table 1 for details on each system.) The dihedrals shown are identical to those in Figure 9 of ref 18; they consist of the dihedral angles between CR and the side-chain carboxylate carbon (black), between the amine nitrogen and the second side-chain carbon Cγ (red), and between the backbone carboxylate carbon and Cγ (green). We computed the OPLS/AA distributions over 8 ns of total simulation time and the AIMD distributions over 40 ps.

P{cos[θ(rij, rik, rjk)]} )

∑∑∑ i)1 j)1 k 5 Å-1, demonstrating that AIMD gives a good description of the intramolecular contributions to the structure factor. However, overstructuring of AIMD data beyond the first solvation shell leads to a more pronounced peak

Structure of Aqueous Solutions of MSG

Figure 13. Real-space structure factor F(r) calculated from the weighted sum (eq 3) of all atom-atom gRβ(r)’s for the fully deuterated MSG solution (top) and for the fully protiated solution (bottom). EPSR curve (dotted line) from the weighted sum of atom-atom gRβ(r)’s obtained from EPSR analysis of the neutron diffraction data in ref 17. Other line colors are identical to those in Figure 12. Note that only the AIMD simulation and experimental raw data include the intramolecular contributions to F(r).

in the deuterated F(Q) distribution near Q ) 2 Å-1. When the longer-range correlations are instead derived from the classical simulations, a very good match to the raw experimental data is obtained for all Q > 1.5 Å-1 or so. Except for the peak near Q ) 3.5 Å-1 in the protiated case, which is higher in the simulations, our fits to the raw data do not appear any worse than the EPSR fits shown in ref 17. The peaks appear at roughly the same values of Q, and the differences in amplitude could be a result of the broader range of glutamate conformations seen in the neutron diffraction/EPSR data (Figure 1). Features in F(Q) for aqueous solutions at values of Q ≈ 2-3.5 Å are associated with extended hydrogen-bonded water networks,44 so even minor deviations in these peaks indicate that there are differences in the water network seen in the neutron diffraction experiment and in our simulations. The effect of increased concentration on the simulated F(Q) distribution is very small, in agreement with our observation of significant ion-solute and solute-solute correlations, suggesting that some clustering of solutes and ions is occurring in the MD simulations. Interestingly, a Fourier transform of the “raw” experimental F(r) distribution regenerates the original F(Q) distribution very closely (dashed orange line in Figure 12), right down to the statistical noise seen in the protiated case. This points out the difficulty of direct comparison to experimental data expressed in real space, as some of the features in F(r) arise as spurious Fourier components of the noisy data in F(Q). This is another advantage of the EPSR technique, as it generates solution structures that might be analyzed to obtain F(Q) [and associated atom-atom g(r) distributions] using the same statistical methods that are used in other simulations. To understand solution structure it is still important to view data in real space, since the peaks seen in F(Q) at particular values of Q cannot be interpreted as intuitively as a similar peak in F(r), although experimentalists trained to interpret features in reciprocal space may disagree. On the other hand, as soon as EPSR is brought to bear on the problem, it can no longer be claimed that the results are

J. Phys. Chem. B, Vol. 113, No. 21, 2009 7697 from a direct experimental probe of the solution. It is possible that quite different solution structures can reproduce the same F(Q), as is the case here. Despite the large differences in water-water correlations especially, EPSR and MD lead to very similar expressions for the total structure factor. The differences are mainly confined to the two peaks near Q ) 3 Å-1 and Q ) 4.5 Å-1 in the protiated cases, and these are likely a signature of differences in the extended hydration network.44 All of our simulated F(r) distributions bear a closer resemblance to the neutron diffraction/EPSR data for pure water80 than to the experimental results for aqueous monosodium glutamate. In agreement with the water-water RDFs shown in Figure 7, the apparent huge disruption of the water network observed in the neutron diffraction experiments is not reflected by any of our simulation methods. The fact that we can fit the experimental structure factor, F(Q), quite well with our simulation data, except for the range of Q associated with the network of solvent water, suggests that most aspects of the solution structure are being accurately modeled. The most likely explanation of the differing results is that the degree of Na+-COO- and NH3+-COO- association is somewhat overestimated in the classical MD simulations, so that, although the first solvation shell of the glutamate is accurately modeled, at higher concentrations, the water remains somewhat partitioned from the solute-ion pairs and, hence, the water-water correlations change only slightly compared with the experimental results. Another possibility we considered is that the underestimated system density in the EPSR simulations might, in turn, cause the degree of water-water correlation to be underestimated. Although certainly true, the degree to which this would change the results is likely less than the systematic error in the experiments,81 and hence, it is probably not to blame for the difference in the water-water correlations. Another issue preventing unambiguous interpretation of the experimental results is the noise in the F(Q) distribution measured by neutron diffraction for the protiated system, which causes spurious peaks to appear in F(r), especially seen in the raw F(r) data between 2 and 3 Å. In this region, the background and inelasticity correlations to the data are most difficult to remove when light hydrogens are present in the sample. AIMD simulations have allowed us to include accurate intramolecular water-water correlations in our structural determinations, and we have demonstrated the utility of this approach. We might also hope that an ab initio method such as AIMD would be a better model in general and that, if it were possible with current computational resources to extend these calculations to larger concentrated systems, the results would more closely match the experiments. However, the fact that the AIMD calculation of F(r) appears farther from the experimental results than the classical results does not give us confidence that this would be the case. One possibility that we reiterate here is that the atomic density of the AIMD system (0.107/Å3) is likely to be too large. This would lead to an overemphasis on the longer-range correlations observed in the system, leading to the more structured appearance of the AIMD-derived F(r) in Figure 13. Our neglect of the counterion in the AIMD simulations will also cause the results to be somewhat different. Further development of more accurate DFT functionals, as well as techniques to do quantummechanics-level simulations on larger systems more efficiently, should lead to better agreement between simulations and experiments. One promising idea is the LocalSCF semiempirical quantum mechanical method,82 which should make long (>100ps) quantum mechanical simulations of concentrated aqueous

7698

J. Phys. Chem. B, Vol. 113, No. 21, 2009

solutions tractable. Also, the inclusion of dispersive interactions in the LocalSCF method, and in other DFT-based approaches,83 should make these quantum mechanical methods more accurate at resolving long-range structures. Even without theoretical advances, recent success in simulating highly concentrated LiCl45 and NaOH44 aqueous solutions suggests that an ab initio simulation of a concentrated MSG (or another amino acid) solution cannot be too long in coming. IV. Conclusions We have completed a thorough study of the structure of hydration and association in monosodium glutamate solutions, with both classical and ab initio molecular dynamics simulations in a dilute solution and with classical methods in a concentrated solution. We found that the classical simulations using the OPLS/AA46-48 force field predict backbone and side-chain dihedral angles that are in close agreement with neutron diffraction experiments in combination with empirical potential structure refinement (EPSR) analysis.3 Simulations with the Amber ff94 force field52 were not in as close agreement, and we attribute this difference to the lack of carefully parametrized side-chain dihedral potentials in the ff94 force field, although it is also possible that our approximate method to arrive at the partial charges of zwitterionic glutamate based on the average of C- and N-terminal glutamate could lead to poor performance, given that the conformational energies are sensitive to small changes in the charges. More accurate parametrization of the partial charges for the zwitterion using the RESP method,84 for example, could answer this question; however, it was not the purpose of this work to argue for or against different biomolecular force fields, so we did not perform such an analysis. Our use of SPC/E water instead of TIP3P could also have an effect on the results, although we reiterate that generally good results were obtained with SPC/E water in combination with both the Amber and OPLS/AA force fields. The improved treatment of water-water interactions in SPC/E is more important for comparison with the neutron diffraction data than small discrepancies in water-peptide interactions; it is unlikely that we would obtain such a good comparison with the experimental F(Q) distribution if we had used TIP3P water. We compared our structural results with available neutron diffraction data.17,18 We executed simulations at conditions similar to the experimental conditions. However, there were some differences between the systems that should be kept in mind when interpreting the results. Most importantly, the densities in both our MD simulations and the EPSR simulations of McLain et al.17 could be in closer agreement with the experimental data.55 In addition, the computationally challenging AIMD simulations are limited to a single glutamate ion in solution, without a counterion, forcing us to neglect glutamateglutamate and glutamate-ion correlations. These differences notwithstanding, we are convinced that our classical MD results at two different concentrations and including AIMD data as well provided an adequate quantitative and qualitative picture of the MD results at both high and low concentration that can be compared with available neutron diffraction/EPSR results.17 For the most part, there is general agreement between the spatial distribution of glutamate hydrating water detected in our simulations and in the neutron diffraction/EPSR experiments.17 Spatial distribution functions (SDFs) and hydration numbers are similar, showing significant hydration of both the cationic and anionic residues. Both neutron diffraction data and simulation results also show that ion-solute interactions are present and perturb the hydration structure of the system. However, these

Daub et al. interactions appear to cause different perturbations depending on the analysis method. In the experiments, the high solute and ion concentrations mainly perturb the water-water interactions, leading to a dramatic lowering of gOWOW(r). By contrast, in our simulations, raising the concentration leads to a drastic decrease in all hydration numbers as the presence of the ions excludes hydration water, but little difference in the water-water RDFs is seen. Although the density difference between our two studies might be blamed for some of these contrasting results, it is unlikely to be the major culprit. The EPSR results17 suggest a marked disruption of the water-water correlations, including a lowering of the first peak in gOWOW(r) by almost a full unit and an inward shift of the second peak. The simulations show only a modest change in gOWOW(r) from results obtained for pure water, and most of the changes appear outside the first coordination shell. Other neutron diffraction experiments have also reported large disruptions in water structure in concentrated ionic solutions,9,10 but not in a concentrated solution of neutral proline molecules,16 suggesting that the presence of counterions is an important factor. Unfortunately, there is a dearth of simulation studies on concentrated ionic solutions. One simulation study74 found that the first peak in gOWOW(r) changed only slightly until the ionic concentration became quite large (∼5-10 M), and this is consistent with the only two studies of which we are aware that simulate highly concentrated ionic solutions with ab initio methods.44,45 Another simulation study85 did show some significant decrease in the first peak, but not to nearly the same degree as the experimental data suggest, both in the current system17 and in concentrated hydroxide solutions.10 The main reasons for the discrepancies between our MD results and the neutron diffraction/EPSR data appear to be related to the inadequacy of using classical potentials to describe (a) zwitterionic amino acids, which differ in important ways from amino acids in a protein or polypeptide,25 and (b) higher solute concentrations, where intersolute effects arise and ion-solute interactions become more prevalent. The EPSR simulations were run using very similar force fields, but in EPSR, the additional input of the experimental structure factor(s) and the imposition of other constraints80 makes exact parametrization of the force field less crucial. Clearly, more simulation work and potential model development for concentrated ionic solutions is necessary. The classical models that we used lack explicit polarization effects, and this omission could be the cause for the overestimation of binding between COO- and NH3+ and between COO- and Na+. Interestingly, the calculated structure of a concentrated ionic solution using a polarizable water model88 did not demonstrate large changes in the water-water correlations, suggesting that it is more important to include polarizable ions87,88 than polarizable water. Alejandre and Hansen have suggested that a nonadditive modification of the Lennard-Jones cross terms is necessary to correctly model the transition from saturated solution to the crystallization of the ionic solid in a model of NaCl.89 The observation of solute clustering through NH3+-COO- binding in neutron diffraction and MD studies of aqueous dipeptide solutions20 demonstrates that reliable empirical force fields for amino acid solutions must accommodate very different hydration behaviors in quite similar systems. Another area of discrepancy between EPSR analysis of the neutron diffraction experiments and our simulations resides in the hydration of the carboxylate oxygens. Previous34,35 and present AIMD simulations of the carboxylate anion group yield hydration numbers in good agreement with experiments,17,41,42 and it is consistently found that, on average, less than 0.1 water

Structure of Aqueous Solutions of MSG proton simultaneously resides in the hydration shell of the two oxygen atoms in each carboxylate group. This finding is consistent with our classical-force-field results. It is also consistent with earlier quantum chemistry calculations, which show that sharing a water proton among two carboxylate oxygen atoms is not among the most energetically favorable configurations.38 A shared-water configuration also restricts water entropy. Thus, both energetic and entropic considerations are consistent with the AIMD-predicted scarcity of hydration water sharing. Our simulation results thus appear to be at odds with EPSR analysis of recent neutron diffraction data.17 Hence the improvement of at least one aspect of EPSR analysis might be advantageous. We suggest that experimental investigation of simpler examples of carboxylate-containing ions, such as formate or acetate, might be a useful step forward in this regard. This would reduce the number of pair correlations that are not directly measured by the neutron diffraction experiments, thus removing possible artifacts introduced by the EPSR analysis. We have demonstrated that all of our simulations (systems 1-3 in Table 1) lead to very similar total structure factors that match the neutron diffraction data quite well, except for the small differences associated with discrepancies in the extended hydrogen-bond network. It is possible that the EPSR-derived structure of the carboxylate hydration shell represents only one of several possible structures that are consistent with the total structure factor. Soper examined this issue in other neutron diffraction/EPSR work.90 Stillinger and Torquato also addressed it by analytical work on lattice models, showing that even in very simple model systems, multiple particle configurations can lead to identical pair correlation functions.91 The experimental picture of liquid structure should benefit greatly from the addition of X-ray diffraction data from these systems.44 By combining the neutron and X-ray data together in the EPSR analysis, statistical noise should be lessened, and a detailed picture of the microscopic structure of aqueous solutions free of ambiguities should emerge.81 Acknowledgment. We are grateful to Gren Patey, Peter Kusalik, Alan Soper, and Silvia Imberti for helpful discussions over the course of writing this article. We are especially thankful to Sylvia McLain for sending us original data files of raw detector data, experimental RDFs, and dihedral distributions. This work was supported by a U.S. National Science Foundation grant (CHE-0718724) to A.L. K.L. acknowledges support from the U.S. Department of Energy. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the U.S. Department of Energy’s National Nuclear Security Administration under Contract DE-AC04-94AL8500. References and Notes (1) Howe, R. A.; Howells, W. S.; Enderby, J. E. J. Phys. C 1974, 7, L111. (2) Texler, N. R.; Holdway, S.; Neilson, G. W.; Rode, B. M. J. Chem. Soc., Faraday Trans. 1998, 94, 59. (3) Soper, A. K. Mol. Phys. 2001, 99, 1503. (4) McLain, S. E.; Soper, A. K.; Luzar, A. J. Chem. Phys. 2006, 124, 074502. (5) Luzar, A.; Soper, A. K.; Chandler, D. J. Chem. Phys. 1993, 99, 6836. (6) McLain, S. E.; Soper, A. K.; Luzar, A. J. Chem. Phys. 2007, 127, 174515. (7) Soper, A. K.; Luzar, A. J. Phys. Chem. 1996, 100, 1357. (8) Botti, A.; Bruni, F.; Imberti, S.; Ricci, M. A.; Soper, A. K. J. Chem. Phys. 2004, 120, 10154. (9) Imberti, S.; Botti, A.; Bruni, F.; Cappa, G.; Ricci, M. A.; Soper, A. K. J. Chem. Phys. 2005, 122, 194509. (10) McLain, S. E.; Imberti, S.; Soper, A. K.; Botti, A.; Bruni, F.; Ricci, M. A. Phys. ReV. B 2006, 74, 094201.

J. Phys. Chem. B, Vol. 113, No. 21, 2009 7699 (11) Mancinelli, R.; Botti, A.; Bruni, F.; Ricci, M. A.; Soper, A. K. J. Phys. Chem. B 2007, 111, 13570. (12) Mancinelli, R.; Botti, A.; Bruni, F.; Ricci, M. A.; Soper, A. K. Phys. Chem. Chem. Phys. 2007, 9, 2959. (13) Kuntz, I. D. J. Am. Chem. Soc. 1971, 93, 514. (14) Ji, P. J.; Feng, W. Ind. Eng. Chem. Res. 2008, 47, 6275. (15) Zhang, Y.; Zhang, P.; Ford, R. C.; Han, S. H.; Li, J. C. J. Phys. Chem. B 2005, 109, 17784. (16) McLain, S. E.; Soper, A. K.; Terry, A. E.; Watts, A. J. Phys. Chem. B 2007, 111, 4568. (17) McLain, S. E.; Soper, A. K.; Watts, A. J. Phys. Chem. B 2006, 110, 21251. (18) McLain, S. E.; Soper, A. K.; Watts, A. J. Phys. Chem. B 2006, 110, 24275. (19) McLain, S. E.; Soper, A. K.; Watts, A. Eur. Biophys. J. 2008, 37, 647. (20) McLain, S. E.; Soper, A. K.; Daidone, I.; Smith, J. C.; Watts, A. Angew. Chem., Int. Ed. 2008, 47, 9059. (21) Troitzsch, R. Z.; Martyna, G. J.; McLain, S. E.; Soper, A. K.; Crain, J. J. Phys. Chem. B 2007, 111, 8210. Troitzsch, R. Z.; Tulip, P. R.; Crain, J.; Martyna, G. J. Biophys. J. 2008, 95, 5014. (22) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (23) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (25) Hugosson, H. W.; Laio, A.; Maurer, P.; Rothlisberger, U. J. Comput. Chem. 2006, 27, 672. (26) Jin, X. Z.; Chao, K. C. J. Chem. Eng. Data 1992, 37, 199. (27) Christensen, H. N. Physiol. ReV. 1990, 70, 43. (28) Bellisle, F. Neurosci. BiobehaV. ReV. 1999, 23, 423. (29) Li, X. D.; Staszewski, L.; Xu, H.; Durick, K.; Zoller, M.; Adler, E. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 4692. (30) Alagona, G.; Ghio, C.; Kollman, P. J. Am. Chem. Soc. 1986, 108, 185. (31) Alagona, G.; Ghio, C.; Kollman, P. A. J. Mol. Struct. (THEOCHEM) 1988, 43, 385. (32) Bruge, F.; Bernasconi, M.; Parrinello, M. J. Am. Chem. Soc. 1999, 121, 10883. (33) Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986, 90, 2174. (34) Leung, K.; Rempe, S. B. J. Am. Chem. Soc. 2004, 126, 344. (35) Leung, K.; Rempe, S. B. J. Chem. Phys. 2005, 122, 184506. (36) Liang, T.; Walsh, T. R. Phys. Chem. Chem. Phys. 2006, 8, 4410. (37) Liang, T.; Walsh, T. R. Mol. Simul. 2007, 33, 337. (38) Sagarik, K.; Dokmaisrijan, S. J. Mol. Struct. (THEOCHEM) 2005, 718, 31. (39) Campo, M. G. J. Chem. Phys. 2006, 125, 114511. (40) Leenders, E. J. M.; Bolhuis, P. G.; Meijer, E. J. J. Chem. Theory Comput. 2008, 4, 898. (41) Kameda, Y.; Fukuhara, K.; Mochiduki, K.; Naganuma, H.; Usuki, T.; Uemura, O. J. Non-Cryst. Solids 2002, 312, 433. (42) Kameda, Y.; Mori, T.; Nishiyama, T.; Usuki, T.; Uemura, O. Bull. Chem. Soc. Jpn. 1996, 69, 1495. (43) Kameda, Y.; Sasaki, M.; Amo, Y.; Usuki, T. Bull. Chem. Soc. Jpn. 2007, 80, 1746. (44) Megyes, T.; Balint, S.; Grosz, T.; Radnai, T.; Bako, I.; Sipos, P. J. Chem. Phys. 2008, 128, 044501. (45) Petit, L.; Vuilleumier, R.; Maldivi, P.; Adamo, C. J. Chem. Theory Comput. 2008, 4, 1040. (46) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (47) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (48) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. J. Phys. Chem. B 2001, 105, 6474. (49) Born, B.; Kim, S. J.; Ebbinghaus, S.; Gruebele, M.; Havenith, M. Faraday Discuss. 2009, 141, 161. (50) Qvist, J.; Persson, E.; Mattea, C.; Halle, B. Faraday Discuss. 2009, 141, 131. (51) Daub, C. D.; Luzar, A., to be submitted. (52) 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. (53) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (54) Aqvist, J. J. Phys. Chem. 1990, 94, 8021. (55) Suzuki, Y.; Matsuo, H.; Koga, Y.; Mukae, K.; Kawakita, T.; Sawamura, S. High Pressure Res. 2001, 21, 93.

7700

J. Phys. Chem. B, Vol. 113, No. 21, 2009

(56) Hammer, B.; Hansen, L. B.; Norskov, J. K. Phys. ReV. B 1999, 59, 7413. (57) Kresse, G.; Furthmuller, J. Phys. ReV. B 1996, 54, 11169. (58) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15. (59) Blochl, P. E. Phys. ReV. B 1994, 50, 17953. (60) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758. (61) Asthagiri, D.; Pratt, L. R.; Kress, J. D. Phys. ReV. E 2003, 68, 041505. (62) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (63) Bergman, D. L.; Laaksonen, L.; Laaksonen, A. J. Mol. Graphics Modell. 1997, 15, 301. (64) Andreani, C.; Menzinger, F.; Ricci, M. A.; Soper, A. K.; Dreyer, J. Phys. ReV. B 1994, 49, 3811. (65) Svishchev, I. M.; Kusalik, P. G. J. Chem. Phys. 1993, 99, 3049. (66) Shelley, J. C.; Sprik, M.; Klein, M. L. Langmuir 1993, 9, 916. (67) Vaisman, I. I.; Berkowitz, M. L. J. Am. Chem. Soc. 1992, 114, 7889. (68) Mountain, R. D.; Thirumalai, D. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 8436. (69) Soper, A. K.; Castner, E. W.; Luzar, A. Biophys. Chem. 2003, 105, 649. (70) Aziz, E. F.; Ottosson, N.; Eisebitt, S.; Eberhardt, W.; JagodaCwiklik, B.; Vacha, R.; Jungwirth, P.; Winter, B. J. Phys. Chem. B 2008, 112, 12567. (71) Lukovits, I.; Karpfen, A.; Lischka, H.; Schuster, P. Chem. Phys. Lett. 1979, 63, 151. (72) Luzar, A.; Chandler, D. J. Chem. Phys. 1993, 98, 8160. (73) Ferrario, M.; Haughney, M.; McDonald, I. R.; Klein, M. L. J. Chem. Phys. 1990, 93, 5156. (74) Chowdhuri, S.; Chandra, A. J. Chem. Phys. 2001, 115, 3732.

Daub et al. (75) Schwegler, E.; Grossman, J. C.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 121, 5400. (76) Petersen, P. B.; Saykally, R. J. Annu. ReV. Phys. Chem. 2006, 57, 333. (77) Soto, A.; Arce, A.; Khoshkbarchi, M. K. Biophys. Chem. 1998, 74, 165. (78) Xenides, D.; Randolf, B. R.; Rode, B. M. J. Mol. Liq. 2006, 123, 61. (79) Sears, V. F. Neutron News 1992, 3, 29. (80) Soper, A. K. Chem. Phys. 2000, 258, 121. (81) Soper, A. K. , Personal communication, ISIS Facility, Rurtheford Appleton Laboratories, Chilton, Didcot, Oxon OX11 0QX, UK, 2008. (82) Anikin, N. A.; Anisimov, V. M.; Bugaenko, V. L.; Bobrikov, V. V.; Andreyev, A. M. J. Chem. Phys. 2004, 121, 1266. (83) Lin, I. C.; Rothlisberger, U. Phys. Chem. Chem. Phys. 2008, 10, 2730. (84) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (85) Mountain, R. D.; Thirumalai, D. J. Phys. Chem. B 2004, 108, 19711. (86) Zhu, S. B.; Robinson, G. W. J. Chem. Phys. 1992, 97, 4336. (87) Harder, E.; Kim, B. C.; Friesner, R. A.; Berne, B. J. J. Chem. Theory Comput. 2005, 1, 169. (88) Kaminski, G. A.; Stern, H. A.; Berne, B. J.; Friesner, R. A.; Cao, Y. X. X.; Murphy, R. B.; Zhou, R. H.; Halgren, T. A. J. Comput. Chem. 2002, 23, 1515. (89) Alejandre, J.; Hansen, J. P. Phys. ReV. E 2007, 76, 061505. (90) Soper, A. K. J. Phys.-Condens. Matter 2007, 19, 15. (91) Stillinger, F. H.; Torquato, S. J. Phys. Chem. B 2004, 108, 19589.

JP810379M