Structure of Aqueous Proline via Parallel Tempering Molecular

8210

J. Phys. Chem. B 2007, 111, 8210-8222

Structure of Aqueous Proline via Parallel Tempering Molecular Dynamics and Neutron Diffraction R. Z. Troitzsch,† G. J. Martyna,*,‡ S. E. McLain,§,| A. K. Soper,§ and J. Crain†,⊥ School of Physics, The UniVersity of Edinburgh, Mayfield Road, Edinburgh EH9 3JZ, United Kingdom, IBM T.J. Watson Research Center, Yorktown Heights, New York 10598, U.S.A., ISIS Facility, Rutherford Appleton Laboratory, Chilton, Didcot, Oxfordshire OX11 0QX, United Kingdom, Biochemistry Department, UniVersity of Oxford, South Parks Road, Oxford, Oxfordshire OX1 3QU, United Kingdom, and National Physical Laboratory, Hampton Road, Teddington TW11 OLW, United Kingdom ReceiVed: February 22, 2007; In Final Form: May 7, 2007

The structure of aqueous L-proline amino acid has been the subject of much debate centering on the validity of various proposed models, differing widely in the extent to which local and long-range correlations are present. Here, aqueous proline is investigated by atomistic, replica exchange molecular dynamics simulations, and the results are compared to neutron diffraction and small angle neutron scattering (SANS) data, which have been reported recently (McLain, S.; Soper, A.; Terry, A.; Watts, A. J. Phys. Chem. B 2007, 111, 4568). Comparisons between neutron experiments and simulation are made via the static structure factor S(Q) which is measured and computed from several systems with different H/D isotopic compositions at a concentration of 1:20 molar ratio. Several different empirical water models (TIP3P, TIP4P, and SPC/E) in conjunction with the CHARMM22 force field are investigated. Agreement between experiment and simulation is reasonably good across the entire Q range although there are significant model-dependent variations in some cases. In general, agreement is improved slightly upon application of approximate quantum corrections obtained from gas-phase path integral simulations. Dimers and short oligomeric chains formed by hydrogen bonds (frequently bifurcated) coexist with apolar (hydrophobic) contacts. These emerge as the dominant local motifs in the mixture. Evidence for long-range association is more equivocal: No long-range structures form spontaneously in the MD simulations, and no obvious low-Q signature is seen in the SANS data. Moreover, associations introduced artificially to replicate a long-standing proposed mesoscale structure for proline correlations as an initial condition are annealed out by parallel tempering MD simulations. However, some small residual aggregates do remain, implying a greater degree of long-range order than is apparent in the SANS data.

I. Introduction and Motivation Among the 20 amino acids, proline is unique in several respects. Structurally, it is the only cyclic amino acid; the side chain is covalently linked to the backbone amine, forming a pyrrolidine ring. As a result, the conformational freedom of the molecule is highly restricted, leading to bends in polypeptide chains. Proline also has the highest solubility of any amino acid in water (up to 6.5 M), and its solutions exhibit unusual colligative properties.2,3 Unlike most low-molecular-weight compounds, the relative viscosity of aqueous proline is anomalously high particularly above concentrations of 3.5 M.4 Proline also appears to have a special functional role as a bioprotectant (osmolyte) in the sense that it inhibits denaturation of proteins under adverse conditions including dehydration and lowtemperature stress. At relatively high concentrations, proline has also been shown to enhance the solubility of hydrophobic compounds (hydrotropism) including certain proteins.5 Finally, there has been some speculation that proline behaves as a socalled chemical chaperone which promotes correct protein refolding.6 †

The University of Edinburgh. IBM T.J. Watson Research Center. § Rutherford Appleton Laboratory. | University of Oxford. ⊥ National Physical Laboratory. ‡

Proline’s unusual properties have led to the suggestion that the mixtures may possess a complex mesoscale structure in which proline molecules form semi-ordered aggregates in aqueous solution. These associations, it is proposed, resemble polymer-like “hydrophilic colloids”.3 The basic structural elements of these assemblies are assumed to comprise a hydrophobic backbone of stacked pyrrolidine rings with solventexposed hydrophilic groups on the surface.3 Some indirect support for this idea has come also from independent measurements including calorimetric data5 which revealed evidence for low-temperature eutectic phase separation in aqueous proline at a range of concentrations. A similar inference is drawn from Fourier-transform infrared (FTIR) spectra which show a splitting of the COO- asymmetric stretch band on increasing proline concentration, implying that the carboxylate group exists in distinct local environments, one of which is possibly the proposed pyrrolidine stack of ref 5 which is also present in proline crystals.7 There are also reports of a significant increase in light scattering in proline solutions above 1.5 M concentration8 which could be explained by proline clustering. The same authors also report an increase in emission intensity of a hydrophobic dye as a function of proline concentration in aqueous solution. This was interpreted as additional evidence for the formation of supramolecular assemblies of proline having distinct hydrophobic regions.8

10.1021/jp0714973 CCC: $37.00 © 2007 American Chemical Society Published on Web 06/26/2007

Structure of Aqueous Proline In contrast to the above evidence, there are also recent reports which directly argue against the formation of ordered supramolecular proline aggregates. In particular, ref 9 suggests that the unusual properties of proline solutions can be well explained by excluded volume considerations alone with no need to invoke solute association as an important or relevant mechanism. They conclude that the protective effects of proline on enzyme stability arise only from the molecule acting as an inert spacefiller. New light-scattering data suggest that proline, similar, to the series of 1-prime alcohols, exhibits enhanced scattering due to concentration fluctuations and not long-range ordering/ assembly formation.10 Also, recent molecular dynamics simulations have suggested that the system is well mixed11 and that proline association in solution even up to concentrations near the solubility limit is unlikely. However, these results were not compared to or tested against any experimental data, sensitivity to empirical parameters was not explored, nor were limitations of the sampling technique evaluated. At present, there is certainly no conclusive nor direct evidence for room-temperature aggregation or microscale segregation in proline mixtures of the type that has been now firmly established in aqueous methanol by neutron diffraction and supported by computer simulation.12,13 Studies of the hydrophobic/hydrophilic character of proline have also failed to offer a consistent picture. For example, Meirovitch et al.14 observed that, on the basis of free energy of transfer from water to organic solvents, proline emerges as a strongly hydrophobic molecule meaning, that one would expect it to show hydrophobic association and to be found typically in protein interiors. This is not supported by evidence from protein structural studies which suggest that proline residues frequently occur on the exterior of proteins in contact with solvent.14 On the basis of the current literature described above, a thorough investigation of aqueous proline is of interest. Therefore, in the present paper, the solution-state structure and properties of aqueous L-proline are examined using replica exchange molecular dynamics simulations to ensure good sampling while comparisons to solution-phase neutron diffraction experiments and X-ray crystal data are employed to probe model accuracy. The X-ray crystal data are compared to Parrinello-Rahman-type constant pressure MD simulations of two molecular crystals, L-proline and L-proline monohydrate. The comparison between the neutron diffraction data and simulation is made directly through the static structure factor. Several empirical water potentials are used, and the effects of nuclear delocalization on S(Q) are approximately included by judiciously employing the results of path integral molecular dynamics simulations on isolated (gas phase) molecules. The objective of the present study is to draw firm conclusions about the local and mesoscale structure of aqueous proline mixtures by combining simulation and diffraction data, to explore the sensitivity of those conclusions to the empirical water potential and the force field, and to assess the effect of quantum corrections on the liquid structure. II. Computational Models and Methods II.A. Classical Empirical Force Field Model. Classical molecular dynamics (MD) simulations were performed to examine the properties of 2.75 M aqueous proline solution, a 20:1 mixture at T ) 300 K, using the CHARM22 (C22) force field15 and three well-known empirical potentials for water (TIP3P, SPC/E, and TIP4P). These water potentials are all rigid (fixed geometry) and nonpolarizable. TIP3P and SPC/E are both three-site potentials with charges and Lennard-Jones parameters

J. Phys. Chem. B, Vol. 111, No. 28, 2007 8211

Figure 1. Naming convention of the atomic sites in proline and water adopted for the purpose of this paper.

assigned to each of the atom positions. TIP4P is a four-site model in which an additional site associated with the oxygen charge is displaced along the bisector of the HOH angle. It should be noted that the C22 force field is designed for use with TIP3P water model, but it is nonetheless informative to explore the effect of changing the water model. A simulation cell, containing 6860 and 343 water and proline molecules, respectively, was prepared by placing the molecules on a cubic lattice with randomized orientations. The system was then equilibrated in the canonical ensemble for 300 ps at a temperature of T ) 300 K in order to anneal out any unphysical high-energy contacts. At this stage, the configuration was replicated three-fold, in order to create initial conditions for independent simulations of the mixture involving each of the three water models, TIP3P, TIP4P, and SPC/E. After another 60 ps relaxation of each system, production runs of 2 ns were performed in the canonical ensemble (TIP4P, TIP3P, SPC/E) and the isothermal-isobaric ensemble at P ) 1 bar (SPC/E). Extended system methods16-19 were used to sample the canonical and isothermal-isobaric ensembles via MD simulation. Multiple time step integrators16,17,20 incorporating constraints on the high-frequency O-H and N-H bond vibrations (SHAKE18 and RATTLE19), together with a 15.999 au hydrogen mass allowed a 6.0 fs “outer” time step to be employed. Periodic boundary conditions were assumed, and Ewald summation was used to compute long-range interactions enabled by the smooth particle mesh21 approximation to the reciprocal space part. In order to present clearly the results of the studies, it is useful to adopt a naming convention for the atoms in the proline molecule. Our convention is given in Figure 1. II.B. Evaluation of the Static Structure Factor from Simulation. Neutron diffraction (particularly on a set of isotopically labeled samples) is perhaps the most powerful technique for the study of liquid structure in aqueous solutions. The experiment yields the total static structure factor, S(Q), directly from the scattered neutron intensity once various corrections for container, multiple, and inelastic scattering are applied to the raw data. Structural information is obtained in the form of the total real-space pair correlation function, G(r), related to S(Q) via Fourier inversion in an infinite system.22,23 The total pair correlation function is a linear combination of partial atomic pair correlation functions defined by

G(r) )

bibjcicjgij(r) ∑ i,j

(1)

where the bi and ci represent the coherent neutron scattering lengths and mole fractions of atomic species i, respectively. The site-site distributions gij(r) are therefore not directly accessible from a single measurement, and angular correlations, which cannot be reconstructed from molecular connectivity, are not included at all.

8212 J. Phys. Chem. B, Vol. 111, No. 28, 2007

Troitzsch et al.

The total distribution is related to S(Q) via a continuous sine transform of

H(r) )

bibjcicj[gij(r) - 1] ∑ i,j

(2)

only in the thermodynamic limit, N f ∞ where N is the total number of particles in the system. Comparison between experiment and simulation at low Q is therefore most properly explored by direct evaluation of S(Q) from the molecular dynamics trajectories and comparison of the result to the diffraction data, as discussed in the Appendix. At high Q, the continuous transform relation represents a good approximation in the finite system (see Appendix). II.C. Nuclear Quantum Effects on the Static Structure Factor: Gas-Phase Path Integral Simulations. Compounds containing light atoms are known to exhibit quantum effects that impact their structural properties,24 and hence the structure factor. Therefore, a simple, approximate method of incorporating intramolecular quantum effects into S(Q) obtained from classical MD has been developed. First, classical gasphase simulations of isolated proline and water molecules are performed and the classical molecular structure factors S(cl,proline)(Q) and S(cl,water)(Q) generated. Assuming the gas-phase structure is not significantly perturbed in the liquid, the classical [intra]molecular structure factors can be subtracted from the total S(Q), and molecular structure factors [obtained from gas-phase simulations in which nuclear delocalization is included using imaginary-time path integration,25 S(qm,proline)(Q) and S(qm,water)(Q)] can be added. Specifically, path integral molecular dynamics (PIMD)26 computations with 64 pseudoparticles have been performed. The one-molecule systems were equilibrated for 20 ps using PIMD. The time step was set to 0.25 fs, and data were collected during 200 ps of production for each isotopic composition. The quantum-corrected structure factor is defined by

S(qm,total)(Q) ) S(cl,total)(Q) + cproline[S˜ (qm,proline)(Q) S˜ (cl,proline)(Q)] + cwater[S˜ (qm,water)(Q) - S˜ (cl,water)(Q)] (3)

TABLE 1: Unit Cell Parameters of the L-Proline Monohydrate Crystal28 experiment SPC/E TIP3P TIP4P

a/Å

b/Å

c/Å

R/deg

β/deg

γ/deg

20.4 20.8 20.8 20.8

6.2 6.4 6.4 6.4

5.1 5.2 5.2 5.2

90.0 90.0 89.7 89.8

95.8 91.0 89.6 89.6

90.0 89.7 90.3 90.3

neutron-scattering measurements. Having established the credibility of the simulations, the nature of the proline-water and proline-proline hydrogen-bond contacts is unravelled. Last, we comment on the long-range order in the solutions. III.A. Molecular Crystals: Force Field Validation and Structural Motifs. Two molecular crystals, the L-proline monohydrate of ref 28 and the L-proline structure of ref 29, were examined at the statepoint T ) 300 K, P ) 1 atm via Parrinello-Rahman MD simulations.17 These test systems with their restricted conformational freedom allow the interaction potentials to be tested and basic structural motifs to be unambiguously discerned. All the water models used in conjunction with the C22 force field predict reasonable cell lengths and angles in comparison to experiment or the monohydrate (see Table 1). This indicates that at least for the limited class of contacts allowed in the crystal, the C22 force field combined with TIP3P, TIP4P, and SPC/E is valid. Let us consider the structural motifs present in the monohydrate crystal (Figure 2). Here, a network of water-proline contacts formed by carboxyl O and water H hydrogen bonds is prominent. The amine H appear to make no water O contacts. Furthermore, proline-proline contacts are primarily formed by HC‚‚‚OC bonds as opposed to ring-ring contacts, or hydrophobic associations; we will discuss the prominence of each of these in the context of aqueous proline simulations later. In contrast, the C22 model does not yield L-proline crystallographic parameters in satisfactory agreement with experiment (see Table 2). The force field model predicts a structure that is much too expanded, particularly along the a and b crystal axes. This indicates that the non-hydrogen-bonded structural motifs involved in crystal packing given in Figure 3 are not properly

where cR is the concentration of molecule R. Discussion is provided in the Appendix, along with a definition of S˜ (Q). II.D. Parallel Tempering Method. The study of the assembly or stability of complex supramolecular aggregates in these solutions by constant-temperature molecular dynamics (MD) simulations is limited by the efficacy with which the method samples the configurational space. The generation of equilibrium configurations over the simulation time is not guaranteed. Therefore, canonical ensemble parallel tempering (or replica exchange) molecular dynamics27 was employed to improve configurational sampling and ensure valid results. We used 36 equispaced replica at temperatures spanning the range 292 < T < 375 K. III. Simulation Results and Comparison to Experiment In order to develop a complete and convincing picture of proline-water solutions via experiment-simulation comparisons, a variety of results are presented. First, the force field models are validated by comparison to structural data on molecular crystals. Second, in selected cases, parallel tempering MD runs were employed to validate the sampling of the force field at room temperature in the more challenging liquid state. Next, the validity of the force field was stringently probed by detailed comparison of aqueous simulation data to different

Figure 2. Slice through the monohydrate crystal of ref 28 in the cacrystal plane. The dashed lines indicate H-bonds. An elaborate network of proline-water associations, as well as some bifurcated H-bonds are evident.

Structure of Aqueous Proline


TABLE 2: Unit Cell Parameters of the L-Proline Crystal29 experiment7 simulation

a/Å

b/Å

c/Å

R/deg

β/deg

γ/deg

11.6 12.2

9.0 9.5

5.2 5.1

90.0 90.1

90.0 89.8

90.0 90.0

described. It is therefore likely that the C22 model will at least slightly underestimate proline-proline ring association. III.B. Effect of Quantum Corrections on S(Q). The total S(Q) for a fully hydrogenous system is shown in Figure 4, using the TIP4P model as an example and with quantum corrections applied as described in section II.C. The effect of the corrections are visible across the entire Q range but are most prominent in the high-Q range, where the intramolecular contributions dominate. In general, the peaks in S(Q) are lowered. As will be shown later, this behavior tends to bring the quantum-corrected structure factor into closer agreement with the experimental diffraction data. While only the TIP4P data is shown here, the results are representative of all the models, and the effects are largest in the fully hydrogenous samples, as expected. The effect of the quantum correction on the intramolecular structure factors is shown in Figure 4, insets a and b. The assumption intrinsic to the quantum correction is that intramolecular liquid- and gas-phase structures are similar. This is tested by extracting an average intramolecular liquid-state structure factor from the MD simulations of the aqueous solution and comparing it to the corresponding gas-phase structure factor. The results of the comparison are shown for proline in Figure 4, inset c, where it is evident that there is negligible perturbation of the gas-phase structure in solution. Therefore, the subtraction and substitution procedure outlined in section 2.3 above is justified. III.C. Comparison of S(Q) to Experiment and Variation between Empirical Water Models. Structure factors obtained from recent experiments1 and from the molecular dynamics simulations (subject to the quantum correction procedure outlined above) for three isotopically distinct mixtures are shown in Figure 5a-c (see Table 3). These correspond to (a) fully hydrogenated samples, (b) fully deuterated samples, and (c) deuteration only on proline amine hydrogen and water sites. Only the amine hydrogens can undergo H/D exchange with the water proton under ambient conditions, and these substitutions eliminate the possibility that S(Q) is influenced by exchange. In general, it is found that results generated using all three water models (TIP3P, SPC/E, and TIP4P) can account for the main features in S(Q) across the entire Q range and for all isotopic mixtures and provide a fair description of the liquidstate features that form the dominant contributions to the total structure factor. The neat liquid water data is given in Figure 6 for comparison. First, based on the discussion in the preceding section, in all cases, the agreement between simulated and experimental S(Q) is improved upon application of the quantum correction procedure. There are, however, a number of notable differences that remain: First, in the fully hydrogenated sample (Figure 5a) the position of the first peak is well described by all the models, but the peak height is best accounted for by the TIP4P solvated system. Hydration by TIP3P and SPC/E water leads to slight over- and under-structuring, respectively. The height of the second peak at Q ≈ 5.5 Å-1 is slightly overestimated, albeit by the same amount for all the models, relative to the diffraction data, but the peak position is in uniformly good agreement with experiment. Second, in the fully deuterated sample Figure 5b the first peak in S(Q) is observed to be distinctly higher in all the

Figure 3. Backbone contacts in the L-proline crystal of ref 29. Note, the stacking of the pyrrolidine rings. The indicated C-C distances (green) are ∼3.8 Å.

Figure 4. Quantum corrections applied to the total S(Q) from a SPC/E-C22 model simulation of aqueous proline. (Top figure): S(qm,total)(Q) and S(cl,total)(Q). (Bottom figures): (a) S˜ (qm,water,g)(Q) and S˜ (cl,water,g)(Q); (b) S˜ (qm,proline,g)(Q) and S˜ (cl,proline,g)(Q); (c) S˜ (cl,proline,aq)(Q) and S˜ (cl,proline,g)(Q). Here, S˜ (cl,proline,aq)(Q) is extracted from a classical solutionphase computation, “g” denotes gas phase and “aq” denotes solution or aqueous phase.

simulated systems than is observed in the experimental data. The other features in S(Q) appear to be reasonably well described by all the models. Third, for the S(Q) data shown in Figure 5c, only proline amine hydrogen and water sites are deuterated. Here, the features in S(Q) appear to be marginally better described by the SPC/E potential. In particular, the first and second peak positions are obtained correctly. All the models underestimate S(Q) at Q ≈ 6.0 Å-1, and TIP4P appears to perform slightly worse than the others for this isotopic composition. Further, an investigation of the neat liquid water structure factors from experiment and from MD (subject to the same quantum correction as the solution results) is depicted in Figure 6a, H2O, and (b), D2O. While (a) displays significant differences in the two main peaks, agreement in (b) between experiment and simulation is uniformly good. This result is analogous to that observed in the proline solution. Comparing Figure 6a to Figure 5a, we conclude that the overestimation in the second


Troitzsch et al.

Figure 6. Comparison between total structure factors for neat liquid H2O (a) and D2O (b), obtained from simulation and experiment (T ) 300 K, F ) 1 g/mL). At low Q agreement is very good for the deuterated case and is significantly lacking for the hydrogenous case.

to the experimental data than in the neat liquid without a straightforward mapping. Subtraction of total S(Q)’s with different isotopic substitutions from each other can reveal contributions to the scattering data arising from particular pairs of atom types (e.g., particular partial structure factors). Let us define

∆Sred(Q) )

[

]

2 [bhA2S(A) hB2S(B) red (Q) - b red (Q)] 2 bhA + bhB 2

(4)

where the weights Figure 5. Total structure factors for three isotopic substitutions of the proline-water mixture; (a) fully hydrogenated samples, (b) fully deuterated samples, (c) deuteration only on proline amine hydrogen and water sites. The inserts in each figure indicate the substitutions pictorially, where purple and white spheres correspond to deuterium and hydrogen, respectively.

TABLE 3: Isotopic Substitutions Studied by Experiment1 and Simulation Figure 5 panel

backbone

amine

water

(a) (b) (c)

hydrogen deuterium hydrogen

hydrogen deuterium deuterium

hydrogen deuterium deuterium

peak in both can be ascribed to the water component of the system. From a similar comparison between the Figures 5b,c and 6b, the solution simulation seems to be less well matched

N

2 bhA/B

)

∑i

2 bA/B,i

N

are required because unitless S(Q)’s are plotted as discussed in the Appendix. Figure 7a emphasizes the polar contacts, obtained by subtracting Figure 5c from 5a and correcting for water contributions by subtracting the difference of Figure 6a and 6b with appropriate weight. Figure 7b shows the backbone interactions, obtained by subtraction of Figure 5c from 5b. The discrepancies between experiment and simulation are evident in the region Q < 2.5 Å-1 and are absent outside this region. This indicates that any misrepresentations originate from intermediary scale structures, for example, rearrangements of solvation shell structures that would arise from a correct treatment of polarization.30



Figure 7. Comparison of appropriate linear combinations of the total S(Q) given in Figure 5 obtained from experiment and simulation. The top plot emphasizes proline-proline and proline-water H-bonding contacts while the bottom plot emphasizes backbone-backbone and backbone-water contacts.

Overall, and despite the differences discussed above, none of the three water models investigated emerges as substantially better than the others on the basis of direct comparison with the experimentally measured S(Q) data. Thus, all appear to generate at least reasonable proline-water mixtures. The general conclusion is that the main features captured in the diffraction data are also present in the simulated systems. The differences between the simulated data might reflect intermediate scale structural motifs in the liquid that are simply not uniformly better accounted for by any one empirical model. It is therefore useful to present, next, a detailed analysis and comparison of the simulated structures with a view to elucidating motifs that are common to all the models and those which differ. Appropriate comparisons are also made to the solids. III.D. Short-Ranged Proline Structure. A visual inspection of the molecular dynamics trajectories of the simulated systems yields a broadly consistent, structural picture characterized by three distinct types of local proline-proline contacts. These are described below. The first contact is driven by conventional hydrogen bonds (H-bonds) between carboxyl oxygen and amino nitrogen atoms (C ) O‚‚‚H-N) and is shown in Figure 8a. This contact leads to the formation of dimers and longer oligomeric chains. Configurations in which a dimer is tightly bound by two H-bonds are also common. Here each proline molecule acts as both a donor and an acceptor as depicted in Figure 8b.

Figure 8. Short-range proline-proline structural motifs : (a) oligomeric chain segment, (b) (doubly bonded) dimer, (c) bifurcated H-bond, (d) (hydrophobic) backbone contact.

A second type of contact evident in the trajectory snapshots is characterized by a proline-proline bifurcated H-bonds (c.f., Figure 8c). These three-center H-bonds31 have been proposed or observed directly in various systems including molecular crystals and biopolymers32 as well as some molecular liquids,33 including water.34 Finally, close approaches between the pyrrolidine ring segments of the proline molecule occur without facilitating Hbonding interactions. These are considered to be primarily hydrophobically driven contacts (Figure 8d). In order to quantify the importance of the above H-bonded and hydrophobic contact motifs, radial and joint distribution functions, designed such that the major features specific to these motifs can be assigned easily, have been computed. The construction and normalization of appropriate joint distributions have been discussed in detail in ref 35. The basic quantity is a


Troitzsch et al. TABLE 5: Coordination Numbers of Bifurcated and Conventional H-Bonds for Both Proline-Proline and Proline-Water Contacts at Proline’s Carboxyl Groupa HC2 bifurcated HC2 HC3 bifurcated HC3 HW bifurcated HW total HW total bifurcated total nonbifurcated total

Figure 9. Radial and joint radial angular distribution functions of the NP-OC distance and NP-OC distance - cos (NP-OC-OC) angular function respectively.

TABLE 4: Coordination Numbers of the H-Bonding and Backbone (Hydrophobic) Contacts between Proline Molecules HC-OC HC2-OC HC3-OC Cγ-Cγ

TIP3P

TIP4P

SPC/E

EPSR1

rcut [Å]a

0.31 0.12 0.19 0.16

0.20 0.06 0.14 0.16

0.18 0.06 0.12 0.13

0.22 n/a n/a 1.31

2.6 2.6 2.6 4.8

a Here, r cut is defined as the distance to the first minimum in the corresponding g(r).

distribution function of the form g(r,R), where R can be angle or distance, and from which the standard radial distribution g(r) can be recovered by integration via

g(r) )

∫RR

max

J(r,R) g(r,R) d(R)

(5)

min

where J(r,R) is the appropriate Jacobian factor. In Figure 9, both joint and 1d distributions from MD simulation, probing the H-bond dimer configurations given in Figure 8, are shown. This association is characterized by the distance between the nitrogen and the carboxyl oxygens and the cosine of the angle ∠NP-OC-OC. The main peak in the 1d g(r) is split into two main features through the angular component of the 2d function. The features are connected by broad ridges. These can be assigned uniquely, as follows. The H-bond contacts depicted in Figure 8a are confined in the radial component to r < 3.3 Å, and can be assigned the angular range θ ≈ 90-130°, accounting for the relative angular freedom. The lower of the two H-bonds in Figure 8b, as well as the feature contact in Figure 8c, corresponds to the main peak around θ ≈ 66°. Due to the symmetry of the angular quantity, it is clear that this peak must be close to the bifurcated geometry shown in Figure 8c, as discussed further below (Note,

TIP4P

SPC/E

0.01 0.03 0.06 0.05 0.42 2.29 2.71 0.49 2.37 2.86

0.01 0.03 0.01 0.04 0.50 2.28 2.78 0.52 2.35 2.86

EPSR1

a

b 2.1 2.1 2.1 2.1 2.1 2.1

1.67

2.4 2.4 2.4 2.4 2.4 2.4 2.4

1.96

a

Figures 10 and 15 show the distribution function from which the results were computed. Here, a and b (in Å) are the integration limits which define the symmetric shaded areas in Figure 11 used to compute the coordination numbers.

the unique labeling of the two carboxyl oxygens given in Figure 1 was not used to generate the plots given in Figure 9.), e.g. OC17 ) OC16 ) OC. The feature in Figure 9 at θ ≈ 180° or cos θ ≈ -1 on the y-axis, which, combined with the second peak in the 1d distribution around r ≈ 4.9 Å, is due to doubly bonded dimers, which we find from visual inspection (see Figure 8b). In this motif a pair of proline molecules share two H-bonds in which both accept and donate one amine proton, as shown in Figure 8b. More quantitatively, the coordination numbers corresponding to the H-bonded contact can be explored using gHC-OC(r) (see Table 4 top line). The naming conventions are given in Figure 1. It is interesting to note that an inequivalence between the two amine hydrogens is observed in the HC-OC contact. The hydrogen furthest away from the carboxyl group (HC2) forms roughly twice as many inter-proline H-bonds as the other one (HC3), which can be explained by steric hindrance (see Table 4 lines 2 and 3). Let us further consider the H-bond patterns. We observe a high level of bifurcated H-bonds formed with the amine hydrogens in the trajectory snapshots. The relative energetics of conventional two-center and the rarer three-center H-bonds remain a matter of debate.31,32 In principle, two bifurcated bond types are possible: type I involves a single donor and two acceptors (O‚‚‚H‚‚‚O), and type II involves a single acceptor shared by two donors (H‚‚‚O‚‚‚H). In order to explore bifurcation in the solution, we identify corresponding areas in the 2d distributions (Figure 10). Coordination numbers drawn from the relevant distributions are given in Table 5. They are calculated as given in Figure 11, by integration under the shaded areas bounded as indicated. From this analysis, with SPC/E, neither HC3 nor HC2 appear to predominantly donate bifurcated H-bonds,

Figure 10. Radial and joint radial-radial distribution functions of the amine hydrogen-carboxyl oxygen contacts: (left) OC16-HC2 and OC17-HC2 and (right) OC16-HC3 and OC17-HC3.



Figure 11. Schematic representation of the areas in the 2d distributions of Figure 10 defining bifurcated and conventional H-bonds (symmetric around x ) y and marked by dark and light shades, respectively). The integration limits labeled a and b used to obtain the coordination numbers are given in Table 5.

whereas with TlP4P, HC3 donates to an equal population of bifurcated and normal H-bonds. A similar analysis, not depicted, shows that type II bifurcated contacts do not form between prolines. Next, we consider hydrophobic contacts. Defining a coordination number to describe the inter-pyrrolidine ring contact is not as simple as defining a coordination number to describe H-bonding. However, from examination of the molecular crystals, the Cγ-Cγ distance gives a good description of the backbone inter-proline contact. That is, the Cγ carbon is furthest from the polar groups and therefore deemed to be the center of the hydrophobic end (see Figure 1). Figure 12, top shows the gCγ-Cγ(r) radial distribution function. Clearly, the extension of the plateau down to r ) 3.8 Å indicates a presence of a hydrophobic contact, since at this distance no water molecules can act as bridges between the rings. In an effort to pin down orientation of the rings to one another, the joint distribution of the distances of Cγ to the carbons on either side (Cβ, Cδ) is explored in Figure 12 (bottom). The corresponding 1d rdf’s are overlaid. From visual inspection we find that the feature corresponding to the hydrophobic contact occurs chiefly in the region between 3.8 and 4.8 Å on the rCγ-Cβaxis and 5.9-6.8 Å on the rCγ-Cδ-axis. In Figure 8a-c these distance combinations are absent, whereas they are present in 8d. The rest of the plateau in Figure 12 (top and bottom) arises from the backbone-backbone contacts driven by inter-proline H-bond interactions. Finally, it is useful to make at least qualitative comparisons to the structural motifs found in the molecular crystals. In particular, the H-bonded contacts of Figures 2 and 3 are clearly present, as Figures 8 and 9 indicates. A snapshot in Figure 13 illustrates the H-bonding. The crystalline hydrophobic contacts, as outlined in section III.A, are also present in solution (see Figure 12), albeit not through the explicit stacking motif observed in the crystal. Table 6 shows coordination numbers from the solid which can be compared to the liquid results given in Table 5. The large number of Cγ-Cγ contacts in the solid arise from the stacking present in the crystal. In general, since polarization is not included in the model, H-bonded coordination numbers may be somewhat overestimated.30 III.E. Sensitivity of Proline Structure to the Water Model. In order to further assess the proline structure, it is useful to address the dependence of the results on the water model.

Figure 12. Radial and joint radial-radial distribution functions of the backbone contacts. In the top figure, the Cγ-Cγ radial distribution functions for the three models studied are given. Joint distribution functions are presented in terms of distances between Cγ on molecule 1 and the two backbone atoms on either side of Cγ in molecule 2, Cγ-Cδ and Cγ-Cβ in the bottom figure. Close contacts are indicative of the structural motif presented in Figure 8d.

Figure 13. Network structure of aqueous proline. Both direct prolineproline association, and proline contacts mediated by bridging water molecules are present. The bifurcated H-bond connecting the two proline molecules present in the center right portion of the frame is of note.

First, consider the backbone contacts. The coordination number is 0.16 (0.13) for solvation by TIP3P/TIP4P (SPC/E) models. On average this means that one in six prolines (for TIP3P/TIP4P) or one in eight (for SPC/E) form hydrophobic contacts (Table 4).


Troitzsch et al.

TABLE 6: Coordination Numbers of Total (Bifurcated and Conventional) H-Bonds for Proline-Proline and Proline-Water Contacts in the Monohydrate Crystal Cγ-Cγ HW-OC HC2-OC HC3-OC total hydration

TIP3P

TIP4P

SPC/E

rcut

2.30 0.57 0.64 0.25 1.46

2.30 0.25 0.58 0.64 1.47

2.36 0.64 0.57 0.26 1.46

4.86 2.6 2.6 2.6

TABLE 7: Coordination Numbers for OW-OW Contacts in the Solution OW-OW OW-OW

TIP3P

TIP4P

SPC/E

EPSR1

rcut

4.93 4.30

4.63 4.18

4.78 4.21

5.10 4.38

3.54 3.40

Next, consider the case of the H-bonded proline populations, where the sensitivity to empirical water potential is much larger. In the case of the SPC/E solvent, the H-bonded population is only slightly larger than that of the hydrophobic contact. Specifically, the HC-OC coordination number is 0.18 compared to 0.13 for the hydrophobic contact. In contrast, proline in TIP3P water shows a significantly higher population of H-bonded prolines (HC-OC coordination number ≈ 0.3) suggesting that the contact is twice as prevalent in the TIP3P solvent as it is in the SPC/E solvent and also twice as prevalent as the hydrophobic contact. In section III.F, a possible explanation of this dependence based on model water structure is offered. III.F. Hydration and Water Structure. The radial distribution functions of water oxygen atoms, gOW-OW, in aqueous proline solution for the three different empirical water potentials are shown in Figure 14. The corresponding distributions for neat water systems are also shown. Table 7 gives a comparison of the coordination numbers for the OW-OW contact in the solution. It is clear that the TIP3P system is nearly featureless beyond the first coordination shell in both the mixture and the pure TIP3P system. It is well-known that TIP3P water is under structured,36 but the observation may provide some insight into the inter-proline H-bond propensities and solvent sensitivities. Specifically, the formation of an inter-proline H-bond dimer leads to maximum solvent exposure of the pyrrolidine rings. This is energetically unfavorable and leads to a distortion of the surrounding water network. This network is evidently most loosely associated in the TIP3P solvent relative to either TIP4P or SPC/E. The energy penalty for hydrating the proline H-bonded dimer is, therefore, expected to be lowest in TIP3P water and the population of such dimers somewhat higher, as is observed. In the SPC/E and TIP4P water model mixture, the presence of the proline even at a concentration of 1:20, has a clear effect on gOW-OW. Specifically, the second shell is contracted to shorter distance, and the separation between the first and second shell peaks becomes less well defined. The coordination numbers in Table 7 indicate that both the simulations and the experimental inversion (EPSR) are close to one another. In all cases the coordination numbers from EPSR fits to the neutron data are higher than those from any of the models considered here. TIP3P is closest to EPSR, however, the TIP3P g(r)OW-OW is significantly understructured and overall provides a poor description of the water structure. We also note that the SPC/E MD system is under-coordinated compared to EPSR, despite the fact that EPSR uses SPC/E as the reference potential. Next, we investigate the bifurcation of water-carboxyl group H-bonds. The joint distribution describing this contact is shown in Figure 15, and the corresponding coordination numbers are

Figure 14. Oxygen-oxygen (OW-OW) radial distribution functions obtained from the three models of aqueous proline solution studied herein are depicted by the three solid colored lines. The corresponding distribution function predicted by the corresponding three models of neat liquid water are depicted by the three dashed lines.

Figure 15. Radial and joint radial-radial distribution functions for the water hydrogen-carboxyl oxygen contacts, OC16-HW and OC17HW.

included in Table 5. Bifurcation in this contact is not as prevalent as observed in the proline-proline interactions. However, we note that the bifurcation of H-bonds plays a key role in the local structural motifs. Overall, more than one in six H-bonds is bifurcated, which gives rise to a portion of OC’s forming another H-bond, since one hydrogen contact is shared by both carboxyl oxygens. Figure 13 shows how this is incorporated into a selfconsistent network of water hydration and proline contacts that also contains a close approach of the rings. Finally, an investigation of the total carboxyl group solvation predicts that on average 2.86 H-bonds are formed in both TIP4P and SPC/E water respectively. This is notably larger than the value found in the EPSR inversion of the neutron diffraction data.1 Although incorporating polarization will tend to reduce coordination numbers,30 the EPSR prediction of the loss of one H-bond from such a strongly H-bonding center maybe an overestimate. III.G. Long-Range Association. In all the systems studied, local contacts comprising H-bonded dimers, trimers, and rare tetramers are found to coexist with loose apolar associations. These associations/contacts appear to be the only structural motifs that assemble spontaneously on the time scale of the simulations. Aqueous proline is therefore not homogeneously mixed at the molecular level. However, the observed structures are not consistent with the elaborate, stacked pyrrolidine rings proposed by Schobert and others3 to account for the unusual colligative properties of the solutions. In an effort to conclusively



Figure 16. Total structure factors, S(Q), corresponding to the fully deuterated isotopic substitution, Figure 5c. The single-temperature TIP4P-C22 simulation result, as well as the data obtained from the T ) 300 K walker of a parallel tempering simulation are given.

rule out such supramolecular assemblies, we have attempted to introduce them artificially as an initial condition and explore their stability by performing parallel tempering molecular dynamics simulations. Specifically, we start with the TIP4P system and truncate the backbone ring-water interactions by applying an artificially short cutoff of 1.0 Å. This has the effect of biasing the system toward hydrophobic proline-proline contacts. The system is equilibrated under these conditions and spontaneously forms an interconnected network of stacked proline tubes very similar to the proposed structure of Schobert3 (Figure 17, bottom). This is then used as the starting configuration for a first parallel tempering molecular dynamics simulation using the correct CHARMM parametrization. After equilibration, the network of proline was significantly interrupted by water again, especially at higher temperatures. Using the highest temperature configuration of the “melting” PT run, described above, as the starting point for another PT simulation, we find the long-range assemblies are annealed away at all temperatures. Therefore, the structural motif proposed by Schobert2 is unfavorable under C22. Moreover, the S(Q) data obtained from the target (ambient temperature) replica are shown in Figure 16. It is clear that at intermediate Q ( V1/3 in comparison to the result obtained in the thermodynamic limit. The real-space expression for the pair distribution function of a finite system can now be transformed into reciprocal space. Combining eq 7 and the Poisson summation formula,

1

F ˜ (Q ˆ h-1) ∑l F(hl) ) V ∑ Q ˆ F ˜ (Q ˆ h-1) )

∫allspace dr′′ exp(-i2πQˆ h-1r′′)F(r′′)

(8)

Figure 18. Inversion of gOO(r) from SOO(Q) using several reciprocal space cutoffs, Qmax/Å-1 via eq 12. The input SOO(Q) was computed using configurations obtained from a SPC/E model simulation at T ) 300 K and F ) 1 g/mL and hence contains numerical noise.

where 2πQ ˆ h-1 is the reciprocal lattice vector with reciprocal lattice index, Q ˆ , yields

n2(r) )

4πr2 VN

〈|n(Q)|2〉 ∑ Q ˆ

sin Qr Qr

(9)

Here, the atomic (single particle) density function in its reciprocal space form is given by

n(Q) )

∑j eiQ‚r

j

(10)

Taking ΩQˆ to be the number of reciprocal lattice vectors, Q ˆ, with magnitude, Q ) |2πQ ˆ h-1|, in the finite system, it is possible to define the angle and ensemble averaged structure factor,

S(Q) ≡

′

1 NΩQˆ

〈|n(Q′|2〉 ∑ Q ˆ′

(11)

where the prime restricts the sum such that |Q′| ) Q. Hence, in a finite system,

g(r) +

δ(r) 2πr F 2

)

sin Qr

Ω h Qˆ S(Q) ∑ Q ˆ

Qr

(12)

where Ω h Qˆ ) ΩQˆ /N. Note, the sum over lattice vector magnitude is discrete with irregular spacing (Q ) |2πQ ˆ h-1|) as it depends on the integer triplet Q ˆ , and the notation Q ˆ on the sum index is introduced to reinforce this point. Therefore, S(Q) cannot be expressed as a continuous or discrete Fourier sine transform of g(r) (or equivalently a spherical Bessel transform of order zero.). The sum is well behaved because, although the number of allowed Q’s in a given neighborhood increases with system size, the degeneracy decreases correspondingly. It is useful to study the properties of the reciprocal space expression for the radial distribution. First, consider the effect of introducing a cutoff into eq 12 such that Q < Qmax. In Figure 18, the reciprocal space-based expression for the radial distribution function of eq 12 has been computed with several reciprocal space cutoffs, Qmax, and compared to a direct real-space evaluation (e.g., all from simulation data). It is clear that a very large reciprocal space cutoff is required before a direct Fourier



8π3 V

∑ Ωh Qˆ ′ f 4π ∫ dQ′Q′2

Q ˆ ′*0

∆(Q,Q′;σQ) f δ(Q - Q′)

(16)

to generate the standard result

∫0∞ dr4πr2F[g(r) - 1] sinQrQr ) S(Q) - 1 Q > 0

Figure 19. Reciprocal space computation of SOO(Q) compared to an approximate computation based on the Fourier sine transform of gOO(r). The reciprocal space and sine tranform results were generated from configurations obtained from a simulation of 2048 SPC/E water molecules at T ) 300 K and F ) 1 g/mL, and hence, both representations possess numerical noise. The blue curve was obtained by setting the radial distribution function equal to unity beyond rmax ) 10 Å in eq 17. The red curve was generated by setting the distribution equal to unity beyond rmax ) 30 Å. The real space cutoff in the sine transform was kept at rmax ) 30 Å in both cases to preserve the spectral resolution, e.g., the rmax )10 Å result was “zero padded”.

inversion of the structure factor to the radial distribution function becomes possible. These difficulties have bedeviled fluid research for many years, leading to over-interpretations of spurious oscillations and under-estimates of the height of the first peak in the radial distribution function. Next, it is instructive to examine the approach of the formulas appropriate for finite systems to thermodynamic limit. It is a well-defined mathematical operation to compute the Fourier sine transform of eq 12,

∫0R dr4πr2F[g(r) - 1]

{

F

sin Qr

[

∑ Ωh Qˆ ′S(Q′) 4π ∫0 Q ˆ ′*0

R

]}

sin Qr sin Q′r

dr

Q

Q′

- 1 (13)

where R is the prelimit factor that takes the system size to infinity. Clearly, the limit R f ∞ cannot be taken without changing the sum over the discrete Q′ to an integral or equivalently taking the system to size infinity limit in reciprocal space simultaneously. In order to develop a correct limiting formalism, it is useful to define the prelimit δ function,

∆(Q,Q′;σQ) ≡

2 π

∫0R dr sin Qr sin Q′r

1

{

8π

4πQ2 V

3

sin Qr

∑R NR ) N

∑ Ωh Qˆ ′S(Q′) Q ˆ ′*0

N V

F)

gRβ(r) )

[

1

4πFr2cRcβ

[

NR

1

〈

(18)

N

]

×

]

Nβ

∑∑ ∑

δ(|riR - rjβ - hl| - r)〉 (19)

iR)1 jβ)1 l*0iR)jβ

which approaches unity as r f ∞. The notation riR denotes the position of the ith R atom of type R. The associated partial structure factor is defined to be

cRcβSRβ(Q) ≡

′

1 NΩQˆ

nR(Q′)n/β(Q′)〉 ∑ Q ˆ′

〈

(20)

where NR

)

Qr

NR N

cR )

(14)

where σQ ∼ R-1. Inserting eq 14 into eq 13 yields the useful result

∫0R dr4πr2F[g(r) - 1]

It is interesting to consider the approach of the prelimit form of the Fourier sine transformation of the radial distribution function, eq 15, to the thermodynamic limit. Using the EulerMaclaurin summation techniques, it is clear that if the density of states is large around some Q and the prelimit δ-function is strongly peaked at Q′ ) Q, then the continuum limit will be approximately valid. Thus, it can be expected that eq 15 will generate an accurate representation of S(Q) at large Q and fail to generate an accurate representation of S(Q) at small Q. In Figure 19, Fourier sine transformations of radial distribution functions are compared to a direct computation of S(Q) in finite systems of various size. Indeed, the agreement between the two formulas is poor at small Q, but good at large Q. It is therefore quite important when comparing simulation data to experimental results at small Q to avoid the use of the convenient continuum limit and to compute S(Q) directly. The extension to multicomponent systems is fairly straightforward. First, it is useful to define mole fractions and the density

where NR is the number of atoms of type R. The partial radial distribution function is given by

)

Qr

(17)

() Q

Q′

}

∆(Q,Q′;σQ) - 1 (15)

In the thermodynamic limit, it is possible to simplify eq 15 using

nR(Q) )

eiQ‚r ∑ j )1

jR

(21)

R

so that

gRβ(r) +

[

δ(r)

2πr FcRcβ 2

]

1 sin Qr (22) δRβcR ) ΩQˆ SRβ(Q) N Qr


Troitzsch et al.

in a finite system. The quantity SRβ(Q) is real due to the “time reversal symmetry” of the restricted sum. In the thermodynamic limit,

∫0

∞

sin Qr ) Qr cRcβSRβ(Q) - δRβcR

References and Notes

dr4πr2FcRcβ[gRβ(r) - 1]

Q > 0 (23)

as expected. The total cross section in a neutron scattering experiment is

Sneutron(Q) )

′

1 NΩQˆ

〈

|nneutron(Q′)|2〉 ∑ Q ˆ′

(24)

∑j bj eiQ‚r

(25)

where

nneutron(Q) )

j

is the (single particle) atomic density function weighted by the coherent neutron scattering of each atom bj following the first Born approximation. The total cross section can be decomposed into the partials

Sneutron(Q) )

cRcβbRbβSRβ(Q) ∑ Rβ

(26)

Note, Sneutron(Q) approaches ∑RbR2cR as Q f ∞. Therefore, the unitless cross section Sneutron,red(Q) ) Sneutron(Q)/∑RbR2cR which approaches unity asymptotically is presented throughout the paper. In addition, the labels are suppressed for clarity. The associated total distribution function is given by

G(r) )

cRcβbRbβgRβ(r) ∑ Rβ

(27)

However, the unitless quantity, Gred(r) ) G(r)/(∑RcRbR)2, which approaches unity as r f ∞ is typically shown. For simplicity, when dealing with single molecular structure factors, it is useful to define,

S˜ R(Q) )

′

1 NRΩQˆ

|nR(Q′)|2〉 ∑ Q ˆ′

〈

(28)

where NR

nR(Q) )

bi eiQ‚r ∑ i)1

i

time. We thank the U.S.-National Science Foundation for fellowship monies for S.M. under Award OISE-0404938.

(29)

NR is the number of atoms in the R molecule, nR(Q) is the single molecule density, and ri are the constitute atoms in the single molecule. Using this slightly nonstandard definition, it is simple to compute the S˜ R(Q) from single-molecule computations classically and quantum mechanically and employ them in quantum correction procedures. Our definition gives rise to a nonstandard weighting of NR/N ) cR in the correction term. That is, the intensity of the self scattering of the molecule is governed by concentration only, much like the self scattering of single atoms. Acknowledgment. J.C. thanks Scottish Enterprise and the Engineering and Physical Sciences Research Council for support. G.J.M. thanks IBM Research for support and computer

(1) McLain, S.; Soper, A.; Terry, A.; Watts, A. J. Phys. Chem. B 2007, 111, 4668. (2) Schobert, B. Naturwissenschaften 1977, 64, 386. (3) Schobert, B.; Tschesche, H. Biochim Biophys. Acta 1978, 541, 270. (4) Chilson, O.; Chilson, A. Eur. J. Biochem. 2003, 270, 4823. (5) Rudolph, A.; Crowe, J. Biophys. J. 1986, 50, 423. (6) Chattopadhyay, M.; Kern, R.; Mistou, M.-Y.; Dandekar, A.; Uratsu, S.; Richarme, G. J. Bacteriol. 2004, 186, 8149. (7) Kayushina, R.; Vainshtein, B. Kristallografiya 1965, 10, 833. (8) Samuel, D.; Krishnaswamy, T.; Kumar, S.; Ganesh, G.; Jayaraman, G.; Yang, P.; Chang, M.; Trivedi, V.; Wang, S.; Hwang, K.; Chang, D.; Yu, C. Protein Sci. 2000, 9, 344. (9) Jones, G.; Paleg, L.; Winzor, D. Biochim. Biophys. Acta 1994, 1201, 37. (10) Vass, H.; Hossack, W.; Troitzsch, R.; Martyna, G.; Crain, J. Manuscript to be published 2007. (11) Civera, M.; Sironi, M.; Fornili, S. L. Chem. Phys. Lett. 2005, 415, 274. (12) Dixit, S.; Crain, J.; Poon, W.; Finney, J.; Soper, A. Nature 2002, 416, 829. (13) Dougan, L.; Bates, S.; Hargreaves, R.; Fox, J.; Re´at, V.; Finney, J.; Soper, A.; Crain, J. J. Chem. Phys. 2004, 121, 6456. (14) Meirovitch, H.; Rackovsky, S.; Scheraga, H. Macromolecules 1980, 13, 1398. (15) MacKerell, A.; Bashford, D.; Bellott, M.; Dunbrack, R.; Evanseck, J.; Field, M.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D.; Prodhom, B.; Reiher, W.; Roux, B.; Schlenkrich, M.; Smith, J.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (16) Tuckerman, M.; Berne, B.; Martyna, G.; Klein, M. J. Chem. Phys. 1992, 97, 2635. (17) Martyna, G. J.; Tuckermann, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (18) Ryckaert, J.; Ciccotti, G.; Berendsen, H. J. Comput. Phys. 1977, 23, 327. (19) Andersen, H. J. Comput. Phys. 1983, 52, 24. (20) Tobias, D.; Klein, M.; Martyna, G. J. Chem. Phys. 1994, 101, 4177. (21) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (22) Egami, T.; Billinge, S. Underneath the Bragg Peaks: Structural Analysis of Complex Materials; Pergamon Press: Elmsford, NY, 2003. (23) Billinge, S., Thorpe, M., Eds. Local Structure from Diffraction; Plenum Press: New York, 2000. (24) Kuharski, R. A.; Rossky, P. J. J. Chem. Phys. 1985, 82, 5164. (25) Feynman, R.; Hibbs, A. Quantum Physics and Path Integrals; McGraw-Hill: New York, 1965. (26) Tuckerman, M.; Berne, B.; Martyna, G.; Klein, M. J. Chem. Phys. 1993, 99, 2796. (27) Earl, D. J.; Deem, M. W. Phys. Chem. Chem. Phys. 2005, 23, 3910. (28) Koritsa´nszky, T.; R. Flaig, D. Z.; Krane, H.-G.; Morgenroth, W.; Luger, P. Science 1998, 279, 356. (29) Hayashi, Y.; Matsuzawa, M.; Yamaguchi, J.; Yonehara, S.; Matumoto, Y.; Shoji, M.; D, H.; Koshino, H. Angew. Chem., Int. Ed. 2006, 45, 4593. (30) Sprik, M.; Klein, M. L.; Watanabe, K. J. Chem. Phys. 1990, 94, 6483. (31) Yang, J.; Christianson, L.; Gellman, S. Org. Lett. 1999, 1, 11. (32) Yang, J.; Gellman, S. J. Am. Chem. Soc. 1998, 120, 9090. (33) Matsumoto, M.; Tanaka, H.; Nakanishi, K. J. Chem. Phys. 1993, 99, 6935. (34) Walrafen, G. E.; Hokmabadi, M. S.; Yang, W. H.; Chu, Y. C.; Monosmith, B. J. Phys. Chem. 1989, 93, 2909. (35) Whitfield, T. W.; Martyna, G. J.; Allison, S.; Bates, S. P.; Vass, H.; Crain, J. J. Phys. Chem. B 2006, 110, 3624. (36) Mantz, Y. A.; Chen, B.; Martyna, G. J. J. Phys. Chem. B 2006, 110, 3540.

Structure of Aqueous Proline via Parallel Tempering Molecular

Recommend Documents