J. Phys. Chem. 1995,99, 4855-4859
4855
Nature of Denaturing Agents: Monte Carlo Simulations of Bimolecular Complexes Involving Urea and N-Methylacetamide in Aqueous Solution Julianto Pranata Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701 Received: September 30, 1994; In Final Form: January 5, 1995@
Monte Carlo simulations were utilized to investigate the formation of dimers involving urea and N-methylacetamide (NMA) in aqueous solution. The dimers show much variation in structure, including both singly and doubly hydrogen-bonded as well as stacked configurations, depending on their constitution and intermolecular separation. Relative stabilities of the dimers were found to be urea dimer > mixed dimer > NMA dimer. These results suggest that protein denaturation in the presence of urea involves the participation of aggregates of the denaturant.
Introduction Although urea and guanidinium chloride are widely used as protein denaturants, the mechanism of their action is still not well-understood. One viewpoint is that these denaturing agents act by disrupting water structure, thus weakening the hydrophobic forces that are primarily responsible for protein folding. A variety of evidence has been presented supporting the notion that urea is a water "structure other evidence is consistent with "structure-making" proper tie^;^ yet other evidence says it is neithere5 The opposing viewpoint is that these denaturants interact directly with the protein. Support for this includes Breslow and co-workers' extensive studies on the effects of these agents on the solubilities, binding constants, and reaction rates of a variety of systems;6 other evidence comes from calorimetric measurement^.^,^ Recent developments in computer simulation techniques suggest that additional insights on this phenomenon can be obtained using these methods. Molecular dynamics simulations were utilized ten years ago to study a dilute aqueous solution of ~ r e a . ~More J ~ recent Monte Carlo simulations have resulted in much information about the properties of this system." Simulations of guanidinium cation in water have also been reported.12 Monte Carlo simulations were also used recently to study the interaction of urea and guanidinium cation with a benzene molecule;13 these systems serve as models for the interaction of denaturants with aromatic side chains of protein residues. In the present work, we performed Monte Carlo simulations on bimolecular complexes involving urea and N-methylacetamide (NMA), as a model for the interaction of urea with the main chain atoms in a protein molecule.
Methodology and Computational Details Monte Carlo simulations are employed to compute the structural and thermodynamic properties of a system as statistical mechanical averages.l4 The particular system under consideration consisted of one or two solute molecules (urea and/or NMA) plus 324 water molecules in a rectangular box of dimensions ca. 18.6 x 18.6 x 27.9 A3. Periodic boundary conditions were employed; Le., the central box was surrounded by images of itself in all directions. This allowed the simulation of the bulk properties of a fluid while using a finite number of molecules. Instantaneous geometrical configurations of the system were selected using the Metropolis Monte Carlo @
Abstract published in Advance ACS Abstracts, March 1, 1995.
0022-3654/95/2099-4855$09.00/0
alg~rithm,'~ augmented by preferential sampling16 where the solvent molecules closer to the solute were moved more frequently than those farther away. Interactions between components of the system are represented by empirical potential energy functions. OPLS potent i a l ~were ~ ~ used throughout, including the TIP4P model for water, the conformation-dependent parameters for transNMA,19 and the recently reported parameters for urea.11J3 In the simulations, the trans conformation of NMA was used exclusively, because it is the more stable both intrinsically and in aqueous ~ o l u t i o n and ' ~ also because the majority of amide bonds in proteins are in the trans conformation. The OPLS potential functions represent the intermolecular interaction energies by Coulombic and Lennard-Jones terms between sites primarily centered on the nuclei (eq l).17
Thus, the interaction energy between two molecules a and b is given by the sum of the individual pairwise interactions between sites i on molecule a and sites j on molecule b. Parameters A and C are related to the more familiar u and E by the relation Aii = 4 E ~ c J ~and ' ~ Cii = 4 E~U?, and geometric and CQ = combining rules are utilized; Le., Au = (C.,C.,)l/Z* 18 JJ For the simulations, the geometric structures of urea and NMA were optimized using ab initio calculations at the RHF/ 6-3 1G(d) level,20 with empirical corrections provided by incorporating the offset forces proposed by Fogarasi and coworkers.21 In the simulations the solute molecules were treated as planar, rigid entities with all internal degrees of freedom frozen. Normally, sampling of the torsional degrees of freedom would be included, but the high energy barriers (ca. 18 kcal/ mol' associated with the amide bond rotation should make the rigid body approximation a good one. The particular focus of these simulationswas the computation of free energies. This is made possible by incorporating statistical perturbation theory, which allows the calculation of free energy differences (eq 2).24 1322323)
Gj - Gi= -kT ln(exp[-(Hj - Hi)/kT])i
(2)
Thus, the free energy difference between systems i a n d j is expressed as a statistical mechanical average of a function of their energy difference, with the averaging done on an ensemble generated for system i. In practice, eq 2 is only applicable where 0 1995 American Chemical Society
Pranata
4856 J. Phys. Chem., Vol. 99, No. 13, 1995 Urea
+
Urea
+
Urea dimer
\
- _ _ _ _ _NMA dimer
=‘
> Urea
+ NMA
+
Mixeddimer
I
cu
AG4
r
---_mixed dimer
4.0-
- ureadimer 2.0-
5
0.0Lr.
-
.-Y 2 2
BAG (Mixed dimer vs Urea dimer) = AG2 - AGl AAG (NMA dimer vs Mixed dimer) = AG3 - AGl AAG (NMA dimer vs Urea dimer) = AG, - 2AG,
-.-
I
3.0
Figure 1. Thermodynamic cycles illustrating the relationships among the free energy changes accompanying the processes described in the
text. systems i and j are not too different from one another. However, any two systems can be “connected” via a coupling parameter, and multiple simulations are then performed to gradually transform or “mutate” one system to the other. Another limitation is that only free energy differences can be calculated, rather than absolute free energies, although a judicious choice of the systems involved can sometimes allow the calculation of absolute free energies i n d i r e ~ t l y . ~ ~ Two general types of simulations were performed. In the first one, the free energy profile o(r), or potential of mean force (PMF), was calculated as a function of the distance r between the carbonyl carbon atoms of two solute molecules. Calculations of this type were performed where the solutes were two urea molecules, one urea molecule and one NMA molecule, and two NMA molecules. In each simulation the distance r between the carbonyl carbon atoms was constrained, and a perturbation was performed where r was both increased and decreased by 0.1 A, a procedure known as double-wide sampling.26 Equation 2 was used to calculate the associated free energy changes. In this manner, the PMF was obtained for values of r between 3.3 and 6.9 A. A similar calculation of the PMF for the NMA dimer was reported a few years ago;27 the calculation is repeated here for comparison with the other dimers. From the PMF it is then possible, at least in principle, to calculate the association constant K, (eq 3).28
K , = 4 z K r 2 exp[-o(r)/kr] d r
-2.0-
(3)
The difficulty in applying eq 3 lies in the fact that the absolute position of the PMF is unknown (a consequence of the inability to calculate absolute free energies) and in the choice of the cutoff distance c, defined as the distance beyond which the two molecules are not considered to be “associated”, which appears in eq 3 as the integration limit. The second type of simulation involved the mutation of one solute molecule (e.g., urea) into another (e.g., NMA). These mutations were performed on a single solute molecule in solution or on one of two solute molecules or on both solute molecules. When two solute molecules were involved, each was treated independently; Le., they were not constrained to remain a given distance from each other or by any other type of constraint. Each mutation involved a total of 10 simulations with double-wide sampling.26 The free energy changes associated with these interconversions were calculated using eq 2. A thermodynamic cycle can then be constructed (Figure 1) and used to calculate the differences between the free energies of binding for the various dimers.
I 4.0
I
I
I
5.0
6.0
7.0
C - C distance (angstroms) Figure 2. Potentials of mean force for the urea dimer, mixed dimer, and M A dimer. The absolute vertical displacement of each curve is arbitrary.
Each simulation consisted of the generation of lo6 configurations for equilibration purposes, followed by 4 x lo6 codigurations from which statistical mechanical averages were obtained. Statistical uncertainties were calculated by computing separate averages over blocks of 2 x lo5configurations. The simulations were performed in the isothermal-isobaric (“T) ensemble at 25 “C and 1 atm. The intermolecular interactions were switched off at a distance of 8.5 A, with quadratic feathering over the last 0.5 A. The Monte Carlo simulations were performed using the BOSS program29 installed on a Hewlett-PackardApollo Series 700 computer or on a Silicon Graphics Indigo R4000 computer. The ab initio calculations (for geometry optimizations of the solutes) were performed using the TX-90 program,30 implemented on the Hewlett-Packard computer.
Results and Discussion Potentials of Mean Force. PMF’s for the urea dimer, the NMA dimer, and the mixed dimer are shown in Figure 2. Since eq 2 only gives free energy differences, the absolute vertical position of each curve is arbitrary. Therefore, one cannot be certain as to the depth of the free energy well for each dimer. However, it is apparent from the shape of the curves that the urea dimer has the deepest free energy well. The other dimers have considerably shallower minima, with the mixed dimer having a slightly deeper minimum. The lack of a clearly defined minimum for the NMA dimer is consistent with earlier calculations by J ~ r g e n s e n . ~ ~ Examination of the configurations generated in the Monte Carlo simulations reveals much variation in the structures of the dimers. Figure 3 displays representative configurations of the urea dimer at several values of r. Stacked structures are sometimes observed at r = 3.6 A, but this is likely due to the distance constraint rather than an intrinsic preference for such structures. Hydrogen-bonded structures predominate for r between 4.0 and 5.2 A. Most of these have only one hydrogen bond holding the two urea molecules together. Somewhat unexpectedly, cyclic doubly hydrogen-bonded configurations were only observed at distances of 4.0 and 4.2 A, and even then only occasionally. From the PMF (Figure 2) it can be noted that these distances are near the bottom of the free energy minimum, but it should be noted also that the bottom of the well is very broad. In other words, singly hydrogen-bonded structures can be just as stable. Thus, the presence of a second intersolute hydrogen bond does not appear to contribute much
Monte Carlo Simulations of Bimolecular Complexes
d
J. Phys. Chem., Vol. 99, No. 13, 1995 4857
imposition of distance constraints has some effect on generating these structures. Beyond 5.0 8, the structures are unbound. The PMF for the mixed dimer is essentially flat between 4.0 and 5.0 or even beyond. From these observations it can be inferred that, in this system, hydrogen-bonded and stacked structures are about equally probable, but neither one is particularly stable. Qualitatively, the conclusion can be drawn that the urea dimer is the most stable, followed by the mixed dimer, while the NMA dimer is the least stable of the three. A more quantitative estimate can be obtained by applying eq 3, although as already mentioned this involves making assumptions about the absolute r=5.6 position of the PMF and a somewhat arbitrary choice for the cutoff distance c. A reasonable expectation is that the PMF will eventually approach zero asymptotically as r gets large. It was assumed that this is already happening at the largest r considered in the simulations (6.9 A). As for c, the value chosen was 5.4 A, the largest distance for which the intersolute hydrogen bond in the urea dimer was still observed, at least in some configurations. This choice appears to be reasonable at Figure 3. Typical configurations for the urea dimer at various separations: 3.6 A, stacked; 4.0A, doubly hydrogen-bonded; 4.4 A, least for the urea dimer. With these assumptions, the computed singly hydrogen-bonded; 5.6 A, unbound. For clarity only the water association constant for the urea dimer is 15 M-l, while for the molecules nearest to the solutes are shown. mixed dimer and the NMA dimer they are 1.5 and 0.9 M-l, respectively. The computed Ka for the urea dimer corresponds to a binding free energy of about -1.6 kcaVmo1. This is in reasonable agreement with the notion that the existence of a hydrogen bond contributes about -1 kcaVmol to the stability, a value favored by most protein chemists.31 The dimerization of urea and several derivatives in water has been studied by several group^,^^-^^ and using the data from these experiments the free energy of an amide hydrogen bond was recently estimated to be as high as -6 kcaVm01.~~However, this estimate incorporates a lot of assumptions, including the assumption of a cyclic structure with two hydrogen bonds for the urea dimer. As -5.4 discussed above, this assumption is not valid. Structures with only one hydrogen bond can be at least equally stable in aqueous solution. A subsequent analysis from the same group, using an alternative approach, resulted in a value more in line with the conventional view.36 Q The computed Ka for the NMA dimer (0.9 M-l) is too large. Jorgensen's estimates range from 0.0001 to 0.12 M-1;27 the Figure 4. Ty ical configurations for the NMA dimer at various sepaexperimental estimate is 0.005 M-1.37 This illustrates the rations: 3.6 stacked; 4.6 A, stacked; 5.4 A, unbound. problems associated with the application of eq 3, especially in cases where there is not a well-defined minimum. In such cases to the stability of the dimer. In agreement with this analysis, a the choice of c is not only difficult but can also have a large separate Monte Carlo simulation of the urea dimer with no effect on the calculated Ka. For example, a value of 3.8 8, for distance constraints generated both singly and doubly hydrogenc (which is not an unreasonable distance for a stacked dimer, bonded structures. Obviously this is due to the competition although stacked stuctures are still observed at larger separabetween solute-solute versus solute-solvent hydrogen bonds. tions) resulted in a computed Ka of 0.05 M-I. In any case, the For r values greater than 5.0 A, the structures are mostly conclusion remains that the NMA dimer in water is not as unbound, although singly hydrogen-bonded structures are still strongly bound compared to the other dimers. encountered occasionally at separations as much as 5.4 A. Mutation Calculations. Table 1 summarizes the results of The structures observed for the NMA dimer (Figure 4) show the calculations involving interconversions between urea and less variation. For r values below 5.0 A, stacked structures NMA. The precision of the calculation can be judged by predominate; above 5.0 A the structures are unbound. This comparing the free energy differences obtained from dimer preference for stacked structures reproduces Jorgensen's recalculations, AG2, AG3, and AG4. From the thermodynamic sults.27 cycle illustrated in Figure 1, the sum of AG2 and AG3 should Thus, it appears that the urea dimer prefers hydrogen-bonded equal AG4. In fact the discrepancy is 0.5 kcaymol, which is structures while the NMA dimer prefers stacked structures. What within the level of precision that has come to be expected for about the mixed dimer? F i p e 5 displays the structures. At calculations of this type. small values of r, 3.4-4.0 A, stacked structures predominate. Hydrogen-bonded structures are observed for r between 4.2 and As mentioned above, the mutation calculations on the dimers 4.8 A, except, surprisingly, at r = 4.4 A, where stacked involved unconstrained dimers. Consistent with the PMF structures predominate. At r = 5.0 A, both stacked and calculations, the simulations indicate that systems near the urea hydrogen-bonded structures are observed. Undoubtedly, the dimer end point mostly remain bound to each other (typical
8-.
3
A Q8
l
r=4.0
4858 J. Phys. Chem., Vol. 99, No. 13, 1995
Pranata
r=4.4
r=4.8
4
d-
%
P
r=5.4
Figure 5. Typical configurations for the mixed dimer at various separations: 3.6 A, stacked; 4.2 A, hydrogen-bonded; 4.4 A, stacked; 4.8 A, hydrogen-bonded; 5.4 A, unbound.
TABLE 1: Results of Mutation Calculation9 process
free energy change (kcaVmo1)
urea NMA (AG,) urea dimer mixed dimer (AG2) mixed dimer NMA dimer (AG,) urea dimer NMA dimer (AG4)
6.9 f0.2
- -
-
7.9 f0.2 7.7 f0.2 16.1 f0.2
See also Figure 1. The values shown are the average f 1 standard deviation.
distances from 4.5 to 5 A), while systems closer to the mixed dimer and NMA dimer do not (typical distances from 5 to 8
A).
Some concern may arise from the computed value of AG1, the free energy change for converting urea to NMA in water. This value was calculated earlier to be 4.1 kcal/mol.ll The difference arises primarily because the earlier calculation utilized a different partial charge distribution for NMA.19~23 The partial charges used in the present calculation were developed specifically for trans-NMA, and this difference accounts for 1.4 kcall moll9 or about 50% of the discrepancy. The rest could have arisen from differences in the simulation protocol, such as cutoff distances, number of solvent molecules, length of simulation, and so on. From the thermodynamic cycle (Figure 1) and the data in Table 1 we can obtain the relative stabilities of the dimers. The mixed dimer is 1 kcaVmol less stable than the urea dimer, while the NMA dimer is even less stable by an additional 0.8 kcaY mol. These values may be compared to the results from the PMF calculations. The ratio of the Ka's of the urea dimer and the mixed dimer, as calculated from the PMF, corresponds to a AAG of 1.4 kcallmol. The corresponding AAG between the mixed dimer and the amide dimer is 0.3 kcallmol. The discrepancy between the values obtained from the two different methods is quite small and within the method's limits of precision.
Conclusions Monte Carlo simulations indicate that the urea dimer is the most stable bimolecular complex in water, followed by the mixed dimer (one urea and one NMA molecule), and the least stable complex is the NMA dimer. Under the conditions usually employed for urea-assisted denaturation of proteins, Le., high concentrations (up to 8 M) of the denaturant, it would seem likely that the urea molecules would prefer to aggregate with
each other, rather than interacting with the amide functionalities in the protein. While this appears to be consistent with the notion that the denaturant does not interact directly with the protein in the unfolding process, such a conclusion is not warranted. We did not address the posibility of urea interacting with the side chains of proteins, and more importantly,these simulations are limited to bimolecular complexes. The possibility that aggregates of urea molecules might interact with the protein has not been considered. In our opinion this is the most probable scenario. The simulations indicate that the urea dimer is a fairly stable species, but it is not present exclusively or even primarily as a cyclic doubly hydrogen-bonded structure. By not having the second hydrogen bond, the urea dimer has many additional sites for further hydrogen bonding to other solute and solvent molecules. Additional simulations involving more complex solute systems should provide further understanding of the denaturation process, and we intend to attempt such simulations in the future.
Acknowledgment. This work was partially supported by the National Science Foundation through an EPSCOR grant to the Center for Protein Dynamics at the University of Arkansas. References and Notes (1) Finer, E. G.; Franks, F.; Tait, M. J. J. Am. Chem. SOC.1972, 94, 4424. (2) Hammes, G. G.; Schimmel, P. R. J. Am. Chem. Soc. 1967, 89, 443. (3) Barone, G.; Vitagliano, V. J. Phys. Chem. 1970, 74, 2230. (4) Picker, P.; Ledue, P. A.; Philip, P. R.; Desnoyers, J. E. J. Chem. Thermodyn. 1971, 3, 631. (5) Subramaniam, S.; Balasubramaniam, D.; Auhluwalia, J. C. J. Phys. Chem. 1969, 73, 233. (6) Breslow, R. Acc. Chem. Res. 1991,24, 159. Breslow, R.; Guo, T. Proc. Natl. Acad. Sci. U.S.A. 1990,87, 167. Breslow, R.; Halfon, S. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 6916. (7) Bonner, 0.D.; Bednarek, J. M.; Arisman, R. K. J. Am. Chem. SOC. 1977, 99, 2898. (8) Makhatadze, G. I.; Privalov, P. L. J. Mol. Biol. 1992, 226, 491. (9) Kuharski, R. A.; Rossky, P. J. J. Am. Chem. SOC.1984,106,5786, 5794. (10) Tanaka, H.; Touhara, H.; Nakanishi, K.; Watanabe, N. J. Chem. Phys. 1984, 80, 570. (1 1) Duffy, E. M.; Severance, D. L.; Jorgensen, W. L. Zsr. J. Chem. 1993, 33, 323. (12) Boudon, S.; Wipff, G.; Maigret, B. J. Phys. Clzem. 1990,94,6056. {J-$ D J L € E:~-M,: ~ . _2(nwakyk. P. .I,: Jorgensen. W. 1,. .I. Am. Chem. SOC.1993, 115, 927 1.
Monte Carlo Simulations of Bimolecular Complexes (14) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon: Oxford, U.K., 1987. (15) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (16) Owicki, J. C. ACS Symp. Ser. 1978, 86, 159. (17) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. SOC.1988, 110, 3673. (18) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (19) Jorgensen, W. L.; Gao, J. J. Am. Chem. SOC.1988, 110, 4212. (20) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (21) Fogarasi, G.; Zhou, X.; Taylor, P. W.; Pulay, P. J. Am. Chem. SOC. 1992, 114, 8191. (22) Duffv. E. M.: Severance. D. L.: Jorzensen. W. L. J. Am. Chem. Soc. 1992, li4, 7535.' (23) Jorgensen, W. L.; Swenson, C. J. J. Am. Chem. SOC.1985, 107, 1489. . .. (24) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420. (25) Jorgensen, W. L.; Bucher, J. K.; Boudon, S.; Tirado-Rives, J. J. Chem. Phys. 1988, 89, 3742. (26) Jorgensen, W. L.; Ravimohan, C. J. J. Chem. Phys. 1985,83,3050. ,
Y
,
J. Phys. Chem., Vol. 99, No. 13, 1995 4859 (27) Jorgensen, W. L. J. Am. Chem. SOC.1989, 111, 3770. (28) Rue, J. E. J. Chem. Educ. 1969, 46, 12. Justice, M. C.; Justice, J. C. J. Solution Chem. 1976, 5, 543. (29) Jorgensen, W. L. BOSS Version 3.1; Yale University: New Haven, CT, 1991. (30) Pulay, P. ?X-90; University of Arkansas: Fayetteville, AR, 1990. Pulay, P. Theor. Chim. Acta 1979, 50, 299. (31) See, for example: Dill, K. A. Biochemistry 1990, 29, 7133. Pace, C. N.; Heinemann, U.; Hahn,U.; Saenger, W. Angew. Chem., Znr. Ed. Engl. 1991, 30, 343. (32) Schellman, J. A. C. Trav. Lab. Carlsberg, Ser. Chim. 1955, 29, 223. (33) Susi, H.; Timasheff, S. N.; Ard, J. S. J. Biol. Chem. 1964, 239, 305 1. (34) Gill, S. J.; Noll, L. J. Phys. Chem. 1972, 76, 3065. (35) Doig, A. J.; Williams, D. H. J. Am. Chem. SOC.1992, 114, 338. (36) Searle, M. S.; Williams, D. H.; Gerhard, U. J. Am. Chem. SOC. 1992, 114, 10697. (37) Klotz, I. M.; Franzen, J. S. J. Am. Chem. SOC.1962, 84, 3461. JP942655Q
