Destabilization of Hydrophobic Core of Chicken Villin Headpiece in

Article pubs.acs.org/JPCB

Destabilization of Hydrophobic Core of Chicken Villin Headpiece in Guanidinium Chloride Induced Denaturation: Hint of π‑Cation Interaction Sridip Parui, Rabindra Nath Manna, and Biman Jana* Department of Physical Chemistry, Indian Association for the Cultivation of Science, Jadavpur, Kolkata 700 032, India S Supporting Information *

ABSTRACT: Despite their routine use as protein denaturants, the comprehensive understanding of the molecular mechanisms by which urea and guanidinium chloride (GdmCl) disrupts proteins’ structure is still lacking. Here, we use steered molecular dynamics simulations along with the umbrella sampling technique to elucidate the mechanism of unfolding of chicken villin headpiece (HP-36) in these two denaturants. We find that while urea denatures protein predominantly by forming hydrogen bonds with the protein backbone, GdmCl commences unfolding by weakening of the hydrophobic interactions present in the core. The potential of mean force calculation indicates the reduction of hydrophobic interactions between two benzene moieties in 6 M GdmCl as compared to 6 M urea. We observe a near parallel orientation between the guanidinium cation and aromatic side chains of the HP-36 suggesting π-cation type stacking interactions which play a crucial role in weakening of the hydrophobic interaction. We use QM/MM optimization calculations to estimate the energetics of this π-cation interaction. Additionally, the consistency of the unfolding paths between high temperature (400 K) unfolding simulations and steered molecular dynamics simulations strengthens the proposed molecular mechanism of unfolding further.

■

the mechanisms in action.10−12,18−20 The sequences of events in terms of secondary structure melting of several proteins in these two solvents have produced diverse results.21−26 Acid and base catalyzed peptide hydrogen exchange study suggested that urea is able to form an H-bond with the peptide, but GdmCl cannot.27While probing the thermodynamic aspects of the process, Thirumalai and co-workers have shown that the potential of mean force (PMF) between two methane molecules does not get too much affected in the presence of these two denaturants when compared to bulk water. They have, however, shown that the contact minima of the PMF between two charge molecules get destabilized in GdmCl.11 Recently, Zangi et al. have shown the effect of urea on large hydrophobes.28The solubility measurement of different hydrocarbons suggests a profound ability of GdmCl for weakening hydrophobic interactions.29Despite these huge efforts, a comprehensive understanding of their mechanisms of denaturation is still lacking. In this article, we choose a simple model protein, chicken vilin headpiece (HP-36), consisting of three α helices (stabilized by hydrogen bonding interactions) and a clearly defined hydrophobic core consisting of three phenylalanine (PHE) residues to gain insight on the molecular mechanism of denaturation by urea and GdmCl (please see Supporting

INTRODUCTION According to the concept of the protein folding funnel which assumes the principle of minimal frustration, the process of folding proceeds through an ensemble of pathways toward its native folded state. Few pathways emerge as dominant over others due to the topological frustrations present in the protein.1−3The concept of multiple pathways is naturally entangled with multiple intermediates along those pathways. Therefore, uncovering the structural details of these intermediates is of huge interest in the field of protein folding and unfolding.1−6 Besides other techniques, the use of chemical denaturants to explore this fascinating problem has become popular in recent years.7−9The mechanism of protein unfolding can be classified broadly in two classes. In one such mechanism, denaturants directly interact with the protein to weaken the protein− protein interactions (both hydrophilic and/or hydrophobic) responsible for maintaining its three-dimensional structure.10−12 In a different school of thought, it is assumed that denaturant molecules can cause extensive damage of the threedimensional hydrogen bond network of water around protein that plays a crucial role in the stability of proteins’ threedimensional structure13−17and, hence, causes instability of the protein structure.18−20 Among several denaturants, urea and guanidinium chloride are the two most potent denaturants used in unfolding experiments quite routinely. For both the denaturants, there exist both theoretical and experimental studies claiming both © XXXX American Chemical Society

Received: June 22, 2016 Revised: August 19, 2016

A

DOI: 10.1021/acs.jpcb.6b06325 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Unfolding of HP-36 in 6 M urea and 6 M GdmCl. Time evolution of (a) backbone−backbone H-bonds in protein, (b) number of native contacts in protein, and (c) backbone-denaturants H-bonds. (d) Snapshot of representative protein structure showing the H-bonds between protein backbone and urea.

The results show that the numbers of backbone−backbone hydrogen bonds and native contacts are decreasing with time as expected. We also note that the numbers of backbone− backbone hydrogen bonds and native contacts are more in GdmCl compared to that in urea at the end of simulations. These results indicate the possibility of extensive disruption of both tertiary and secondary structures of HP-36 in 6 M urea. The results also infer the presence of considerable residual structure in the unfolded ensemble in 6 M GdmCl solution. Next, we monitor the number of backbone−denaturant hydrogen bonds (Figure 1c) to examine the direct interactions between denaturant and protein during unfolding. In 6 M urea, we find a drastic increase of the number of backbone− denaturant hydrogen bonds from 10 to 35 in the course of unfolding. It is also interesting to note here that the decrease of the number of backbone−backbone hydrogen bonds (Figure 1a) follows closely with the increase of the number of backbone−denaturant hydrogen bonds (Figure 1c) in 6 M urea. Therefore, from a mechanistic point of view, urea forms hydrogen bonds with the backbone of the protein which in turn causes the breaking of backbone−backbone hydrogen bonds in HP-36 leading to the melting of the helices. In contrast, the number of backbone−denaturant hydrogen bonds in 6 M GdmCl hardly shows any significant change during unfolding. This is also in agreement with the smaller change in backbone− backbone hydrogen bond number in 6 M GdmCl compared to 6 M urea. These results also confirm the suggestions of Lim et al. that urea forms an H-bond with the peptide backbone, but GdmCl cannot.27In Figure 1d, it is clear that urea comes close to the backbone and forms extensive hydrogen bonds. Both oxygen and hydrogen of urea participate in hydrogen bonding. Oxygen from urea makes a hydrogen bond with a hydrogen of the backbone amide group of protein, and the hydrogen of urea forms a hydrogen bond with the oxygen of a backbone carbonyl group. We have also calculated the average work profile with

Information Figure S1 for structural classifications). The scenario provides us a unique opportunity to address both hydrophilic and hydrophobic interactions at the same time. We perform standard, steered molecular dynamics (SMD)30−33 and umbrella sampling simulation techniques of chicken vilin headpiece (HP-36) in 6 M urea and 6 M GdmCl solution to explore the unfolding pathways and relevant intermediates. To provide a thermodynamic picture, we also discuss the effect of urea and GdmCl on attractive hydrophobic association of larger hydrophobes (propane−propane, benzene−benzene) by the umbrella sampling method.34 The details of the MD simulation technique are described in the System Setup section.

■

RESULTS AND DISCUSSION Evolution of Intraprotein and Protein Denaturant Interactions During Progressive Unfolding. Unfolding of protein using denaturants generally requires a time scale of milliseconds to seconds, and it can be computationally challenging to explore. Therefore, we have used the strategy to steer the simulations30−33 from native to unfolded state of the protein in denaturating conditions. However, it is important to note that the steering strengths have been kept quite low to reveal the effect of denaturant on the unfolding process of the protein. The root mean square deviation (RMSD) from the native state has been chosen as the steering order parameter. We have performed eight SMD runs (four for urea and four for GdmCl), and the progression of RMSD of HP-36 during all the SMD runs has been shown in Supporting Information Figure S2a. The RMSD started its journey from a very low initial value to 2.5 nm in 20 ns. To monitor the unfolding of protein in terms of both secondary and tertiary structures, we have calculated backbone−backbone hydrogen bonds (Figure 1a and Supporting Information Figure S2b) and the total number of native contacts (Figure 1b) during the progress of simulations. B


Article


Figure 2. Unfolding pathways of HP-36 in 6 M urea and 6 M GdmCl. The unfolding paths are shown in two-dimensional space defined by hydrophilic and hydrophobic SASA. The representative structures of HP-36 along the unfolding pathways are shown.

Figure 3. Evolution of SASA of PHEs and backbone−backbone hydrogen bonding as a function of RMSD during unfolding in (a) 6 M urea and (b) 6 M GdmCl. Contour plot of SASA of PHEs and number of backbone−backbone H-bonds of the protein at different stages of unfolding. The color codes represent different RMSD states.

mechanism of denaturation in these two solvents further, we have calculated the evolution of solvent accessible surface area (SASA) (both hydrophobic and hydrophilic) during the course of unfolding (Supporting Information Figure S3b). In Figure 2, we have shown the path of unfolding on a two-dimensional plot with order parameters and hydrophilic and hydrophobic SASA, along with representative structures at different stages of unfolding. In urea, at the initial stage of unfolding, we find a decrease in hydrophobic SASA while hydrophilic SASA does not show any significant changes (Figure 2 and Supporting Information Figure S4a) indicating a compact structure of HP36. Structural characterization of this compact state in terms of a contact map is presented in Supporting Information Figure S4b. Following this compaction event, we find an increase in both the SASAs with the gradual unfolding of different helices. Unfolding of a second helix is followed by the melting of the first helix, and finally the third helix is melted leading to the

standard deviation as a function of RMSD for two solvent systems (Supporting Information Figure S5a) and the corresponding free energy profile (Supporting Information Figure S5b) using the second order cumulant expansion of the irreversible work (Jarzinsky’s equality).31,32 Please note that a much larger number of trajectories are needed for a better estimation of the free energy. Now, to analyze the behavior of water in both the solvents, we also calculate backbone-water Hbonds during the course of the unfolding (Supporting Information Figure S3a). Unlike a guanidium cation, water is forming a significant number of H-bonds with the backbone of protein in GdmCl as well as in urea in the unfolding ensemble. As described in earlier studies, water entropy change can play an important role in determining the pathways of unfolding.35,36 In our future work, we will investigate in this line. Denaturant Dependent Unfolding Pathways in 6 M Urea and 6 M GdmCl from SMD. To elucidate the C


Article

The Journal of Physical Chemistry B final denaturated state. The melting of all the helices supports that the unfolding of urea proceeds via disruption of backbone−backbone hydrogen bonds which are responsible for helix formation. Two representative structures near the unfolded state also reveal the presence of few hydrophobic interactions. The pathways for all four SMD runs in urea are found to be consistent and have been shown in Supporting Information Figure S2c. In 6 M GdmCl, however, we find a steep increase of the hydrophobic SASA at the initial stage with very little change in the hydrophilic SASA. We find the disruption of the hydrophobic core of HP-36 while structures of all the helices are still preserved. Following this event, a gradual disruption of hydrophobic contacts is observed leading to a comparatively less compact unfolded state. Interestingly, we find that a significant amount of residual secondary structure of helices is present in the unfolded state. A comparatively lesser hydrophilic SASA of the unfolded state in 6 M GdmCl also indicates the lesser disruption of helices. These observations suggest that the mechanism of unfolding by GdmCl proceeds through the disruption of hydrophobic interactions. The pathways for all four SMD runs in GdmCl are found to be consistent and have been shown in Supporting Information Figure S2c. Since HP-36 is a helical protein and its core is mainly made up of aromatic residues, it is important to look at the process of unfolding in terms of backbone−backbone H-bonds in protein and the change of SASA of aromatic residues during the course of unfolding. In urea at the initial stage of simulation we get a decrease in SASA of aromatic residues (Supporting Information Figure S4a). In the contour plot (Figure 3a), for this particular change, we obtain an increase in RMSD (deep blue → light blue) with no significant change in native hydrogen bonds. So this increase in RMSD is only because of the decrease in SASA of aromatic residues. The structural inspection of this state through the contact map (Supporting Information Figure S4b) reveals that PHE-36 forms contacts with the hydrophobic core of HP-36 which is composed of three other PHE residues (7, 11, 18). A representative structure of this state is shown in Supporting Information Figure S4c. Such a type of state is sometimes referred to as a molten globule state.37With an increase in RMSD, both SASA of PHEs and backbone− backbone H-bonds change significantly in urea. However, in GdmCl, only SASA, not backbone−backbone H-bonds, changes significantly during unfolding (Figure 3b). In urea, at RMSD around 2.0 nm, protein unfolds completely with no native H-bonds while in GdmCl at this RMSD value, half of the secondary structure still does not melt though the value of SASA is much higher. This only happens if urea attacks the backbone of the protein and guanidinium cation disrupts the aromatic core of the protein. Validation of Unfolding Pathways Using Umbrella Sampling. To check whether the small biasing force used in our SMD simulations is affecting the pathway or not, we perform umbrella sampling simulation at different values of RMSD in these two solvent systems. The initial structures at different RMSD values are taken from SMD trajectories. Five such configurations have been taken for each denaturant with increasing RMSDs ranging from 0 to 1.5 nm. From each of these umbrella simulation trajectories, we have calculated average values of SASA of PHEs and number of backbone− backbone hydrogen bonds in the protein. The RMSD dependent values of these two quantities are presented in Figure 4. As in SMD simulations, we find that initial unfolding

in 6 M GdmCl proceeds via the disruption of the hydrophobic core made up with PHEs.

Figure 4. Average values of SASA of PHEs and number of backbone− backbone hydrogen bonds of HP-36 at different stages of unfolding in 6 M urea and 6 M GdmCl calculated from umbrella sampling simulations. Different colors are used for different RMSD values. Solid line is used for urea, and dotted line is used for GdmCl.

Effect of Denaturant on the Hydrophobicity Involving Larger Hydrophobes. To provide a thermodynamics view, we have calculated the potential of mean force (PMF) between two hydrophobic solutes in aqueous urea and GdmCl solution. We have calculated the PMF between two methane molecules (not shown here) and find no significant difference in these two solutions compared to bulk water. Such observations are in agreement with the previous results obtained by Thirumalai and co-workers.11 We realize the fact that the hydrophobic cores of the proteins are mostly composed of bigger aliphatic (leucine, isolucine, etc.) or aromatic (phenylalanine, tryptophan, etc.) amino acids rather than alanine. Therefore, we have calculated PMF between two propane molecules and two benzene molecules in bulk water, 6 M urea, and 6 M GdmCl solutions. The results are displayed in Figure 5a,b. For two propane molecules, we find a contact minimum and a solvent separated minimum in water. The contact minimum in 6 M urea becomes comparatively more stable indicating stronger association. In 6 M GdmCl, the depth of the contact minimum (CM) decreases slightly, and that of solvent separated minimum (SSM) increases. This result indicates the increase in stability of the solvent separated pair of propane in GdmCl when compared to bulk water or 6 M urea solution. The desolvation barrier (barrier between CM and SSM), however, increases in comparison to the other two solvent environments. In the PMF between two benzenes, the hydrophobic weakening by GdmCl is more evident. In GdmCl, the depth of CM is drastically decreased, and that of SSM is increased further as a global minimum. This suggests a huge disruption of hydrophobic association between aromatic moieties in 6 M GdmCl. The desolvation barrier is also less here. However, the PMF in urea is found to be negligibly affected when compared to bulk water. We have also calculated a radial distribution function (RDF) between the hydrophobes and the carbon atoms of the denaturants (Figure 5c,d). It shows that more guanidinium cations are housed around the hydrophobes than urea. This result indicates that the water network around hydrophobic residues is affected strongly in 6 M GdmCl rather than in 6 M urea. The corresponding water network for these two solvent systems is shown in Supporting Information Figure S6. We have also probed the local composition fluctuations relative to the average composition ⟨Δχi⟩ = ⟨(χi − χi )/ χi ⟩ ; χi = instantaneous composition of ith species, χi = average composition of the ith species, and ⟨⟩ denotes the ensemble D


Article


Figure 5. Potential of mean force between (a) two benzenes and (b) two propanes in water, 6 M urea and 6 M GdmCl. Radial distribution of carbon atoms of denaturants (c) from benzene and (d) from propane. (e) Schematic description of angle (θ) between two planes. Here, two planes are the guanidinium cation and the phenyl ring. (f) Probability distribution of the angle (θ) in 6 M GdmCl showing the parallel stacking interactions.

∼13°. This indicates that the guanidinium cation is stacked with the benzene in near parallel orientation. Such a stacking interaction can be accounted for by a π-cation type of interaction between the guanidinium cation and benzene. Such π-cation types of interactions are also found to be present in folded protein structures between aromatic amino acids and arginine (PDB ID 3EA1,43 PDB ID 2OR244). In contrast, we do not find any such specific stacking interactions between urea and the benzene molecule (Supporting Information Figure S7) suggesting the absence of π-cation types of interactions. This observation is particularly important because while an internal π-cation type of interactions can stabilize the protein structure, the similar interaction involving a foreign molecule can destabilize the native state. Comparison with High Temperature Unbiased Unfolding Simulations. We have used the SMD30−33 technique to obtain unfolded configurations by applying moderate force. Although this technique is very useful, the unfolded structures generated by this approach can be different from that obtained from equilibrium simulation. However, we use low force to have a minimal effect of mechanical unfolding and a maximum effect of solvents. Yet, to vindicate the solvent effect found in SMD, a further set of simulations has been performed at 400 K. Such a high temperature is so chosen such that all physical processes become faster. Although at this temperature the

average of water and denaturant near the benzene ring in both the solvent systems. We find larger composition fluctuations (⟨ΔχWater⟩ = 0.19 and ⟨ΔχUrea⟩ = 0.34) for constituent species in aqueous urea than in aqueous GdmCl (⟨ΔχWater⟩ = 0.12 and ⟨ΔχGdm+⟩ = 0.20). An enhanced density fluctuation near the hydrophobic solute indicates enhanced hydrophobicity.38,39 Thus, our results of PMF calculation and composition fluctuation analysis are in agreement with each other. Therefore, weakening of the hydrophobic interaction, particularly between two aromatic hydrophobes by GdmCl, is much stronger than that by urea. This weakened hydrophobic interaction in GdmCl correlates well with the disruption of the hydrophobic core of HP-36 by GdmCl. Stacking Interaction between GdmCl and Benzene Destabilizes the Hydrophobic Interaction. To explore the reason behind this stabilization of solvent separated minima as compared to contact minima for two benzene molecules in 6 M GdmCl, we have analyzed the orientation of denaturant around benzenes. Such orientation dependence between different amino acid side chains in the PMF has been studied earlier.40−42 We have defined an angle (θ) between the plane of denaturant and plane of benzene, and it is pictorially described in Figure 5e. We have calculated θ for the guanidinium cation present in the first solvation layer, and the corresponding distribution of θ (Figure 5f) shows a peak at E


Article


Figure 6. (a) Radial distribution function between backbone amide hydrogen (HB) and urea oxygen (OU) (red) as well as between backbone carbonyl oxygen (OB) and guanidinium hydrogen (HG) (green) in the first and last 10 ns of the simulation. (b) Radial distribution function between side chain phenyl group (PHE) and urea carbon (CU) (red) and guanidinium carbon (CG) (green)) in the first and last 10 ns of the simulation. (c) Probability distribution of angle (θ) between plane of guanidinium cation and plane of phenyl ring of PHEs. (d) Representative snapshot from GdmCl induced unfolding trajectory showing π-cation interaction between the aromatic residues of the protein and guanidinium cation. In QM/MM calculation, phenylalanine (PHE18) and guanidinium cation are included in the QM part as highlighted by the circle. Electrostatic potential mapped on the electron density plot on the QM system is shown inside the rectangular box. Red to blue color indicates a negative to positive electrostatic potential.

Figure 7. Comparison of the result of standard molecular dynamics at 400 K with that obtained from steered molecular dynamics. The unfolding paths are presented as a function of hydrophobic and hydrophilic SASA.

the hydrophobic effect of urea and the guanidinium cation, we calculated a time dependent radial distribution of urea and a guanidinium cation around the phenylalanine groups of the protein (Figure 6b). A sharp increase in the first peak of RDF between the phenyl and the guanidinium cation has been found. On the other hand, in the case of urea, no such change is observed. The distribution of θ (as defined in Figure 5c) shows a peak at ∼13° indicating a π-cation type of interaction between aromatic residues and the guanidinium cation in parallel orientation (Figure 6c). The representative protein structure showing the stacking interaction between aromatic amino acids and guanidinium cation is presented in Figure 6d. We have also employed quantum mechanics (QM)/molecular mechanics (MM) optimization calculations to have a better understanding on the electronic nature of the π-cation interaction. In QM/ MM calculation, phenylalanine (PHE18) and one adjacent

unfolded state of the protein is the stable state, it is still interesting to look at how urea and GdmCl affect the course of the unfolding. Supporting Information Figure S6a shows the change in RMSD, and the number of backbone−backbone and backbone−denaturant hydrogen bonds during the unfolding of HP-36 in urea and in GdmCl. As in the case of SMD, we find that urea unfolds the protein by forming a hydrogen bond with the backbone. This can further be supported by the radial distribution function (RDF) between urea oxygen (CU) and backbone amide (HB) which increases significantly during the unfolding (Figure 6a). In contrast, no noticeable change in the RDF is observed between the amide hydrogens (HG) of the guanidinium cation and the carbonyl oxygen (OB) of the protein backbone. This indicates that the urea backbone interaction increases during unfolding, but guanidinium cation−backbone interaction does not alter much. To check F


Article


(iii) standard MD at high temperature (400 K); all of these simulation techniques accelerate the processes in order to explore the structural information on a time scale, easily accessible to MD simulation. SMD was carried out using GROMACS51 with PLUMBED software.56 SMDs were performed, in NPT ensembles, at 300 K and 1 atm along the order parameter RMSD which takes only the heavy atoms of protein for alignment. The temperature was kept constant by a Nose−Hoover thermostat,57,58 and a Parinello−Rahaman barostat59 was applied for pressure coupling. We have also performed an umbrella sampling of the intermediates found in SMD using GROMACS51 with PLUMBED software.56 In each denaturant, five such intermediates have been taken. Like SMD, samplings were carried out with respect to RMSD. In each window, we perform a production run for 2 ns followed by an equilibration for 1 ns. The production run for high temperature MD was performed in NVT ensembles at 400 K. Before the production run was executed, forceful NPT simulations were carried out at 400 K by restraining the position of the protein while allowing the solvents to move. PMF between Propane and Benzene by Umbrella Sampling. To quantify how denaturants affect hydrophobic association at the pair level, we calculated the potential of mean force for propane−propane and benzene−benzene systems along their separation in water, 6 M urea, and 6 M GdmCl. Having known that urea and GdmCl cause stabilization of hydrophobic interaction for methane at pair level,11,60−62 we have taken propane and benzene with increased hydrophobicity (Table 1). Umbrella samplings34 were performed along the axis

guanidium cation are included in the QM part. The rest of the protein and the solvents are treated in the framework of MM. The π-cation interaction energy between PHE18 and adjacent stacked guanidinium cation is found to be −26.3 kJ mol−1 as obtained at the M06-2X/6-31+G(d,p)/OPLS-AA level of theory.45−47 The optimized mean distance between the carbon of the guanidinium cation and the center of mass of the benzene ring of PHE18 is 3.68 Å. The computed π-cation interaction energy is in agreement with the value obtained by Davis et al.48 Electrostatic potential mapped on the electron density plot on the QM system is also shown in the rectangular box of Figure 6d. The red to blue color indicates a negative to positive electrostatic potential. Finally, we have projected the unfolding trajectories onto the plane with hydrophobic SASA and hydrophilic SASA to show the pathways of unfolding (Figure 7). These pathways are compared with those obtained from SMD (Figure 7), and it is evident that the paths followed by SMD and standard MD at 400 K are similar in both the solvent environments. This confirms that the mechanism by which urea denatures protein is mainly by attacking the backbone via formation of a hydrogen bond whereas GdmCl denatures it by attacking the hydrophobic side chain by weakening the hydrophobic interaction through π-cation type interactions.

■

CONCLUSION In this study, we show that urea unfolds HP-36 by primarily making hydrogen bonds with the backbone of the protein. In contrast, in GdmCl, the initial unfolding of protein is governed by the disruption of the hydrophobic interaction present in the aromatic core of the HP-36. Using PMF calculations, we have shown that GdmCl can weaken the hydrophobic interaction between two aromatic molecules extensively. Analysis of organization of denaturants around benzene molecules and around the aromatic side chains of HP-36 along with QM/MM optimization calculations suggest that a π-cation interaction is playing an important role in the unfolding of HP-36 in GdmCl.

Table 1. Force Field Parameters for Propane and Benzene molecule

atomic center

σ

ϵ

Q

propane

CH3 CH2 C H

0.375 0.407 0.358 0.237

0.867 0.410 0.277 0.118

0.00 0.00 −0.14 0.14

benzene

■

of separation between two hydrophobes to obtain the free energy profile. There were 24 separate windows (each 6 ns long) with sampling every 0.05 nm with a span from 0.3 to 1.5 nm. All simulations were performed in NPT ensembles at 1 atm and 300 K. The potential of mean force (PMF) was calculated by the weighted histrogram analysis method (WHAM).63 QM/MM Optimization. The initial guess structure of HP36 for the QM/MM optimization methodology was taken from the MD simulations. All the QM/MM optimization calculations were performed using the micro−macro-iteration scheme64 implemented in the f-DYNAMO library.65 For HP36, the phenyl side chain of a phenylalanine residue (residue id: PHE18) and the adjacent stacked guanidinium cation (residue id: GND1150) were included in the QM part while the rest of the system was considered in the MM part. On the basis of previous reported data on model studies of the π-cation systems, the M06-2X functional45 in combination with the split valence 6-31+G(d,p) basis set46 was employed to optimize the QM region in our QM/MM calculations. On the other side, the rest of system was described with the OPLS-AA force field.47 Herein, π-cation interaction energy (E) was predicted using eq 1

SYSTEM SETUP, SIMULATION, AND QM/MM OPTIMIZATION DETAILS The simulation began with the crystal structure of chicken villin headpiece (HP-36) (protein data bank (PDB) accession code: 1VII).49,50 All of the computer simulations were performed with MD package GROMACS.51 Force field gromos53a652,53 and SPC/E water model were used.54 The parameters for urea and guanidinium chloride were taken from refs 43 and 21, respectively. Long range electrostatic interactions were calculated with the PME method55 with a grid spacing of 0.16 nm (with quadratic interpolation), and a cutoff of 1.0 nm was applied to evaluate van der Waals interaction. The protein was immersed in a 343 nm3 cubic box containing well-equilibrated 6 M urea and 6 M GdmCl. The system was minimized with 1000 steps of steepest descent minimization. We performed position restrained molecular dynamics simulation for 10 ns to hold the protein fixed and to allow the solvents to relax while keeping temperature constant. Protein Simulation. Although molecular dynamics simulations can provide structural information during unfolding, it is a daunting task to accomplish unfolded state results in the time scale of a conventional MD approach at room temperature. We therefore used three ways to overcome this: (i) steered molecular dynamics (SMD),30−33 (ii) umbrella sampling,34and

E = E(total) − [E(cation) + E(π )] G

(1)


Article


(9) Shortle, D. The denatured state (the other half of the folding equation) and its role in protein stability. FASEB J. 1996, 10, 27−34. (10) Hua, L.; Zhou, R.; Thirumalai, D.; Berne, B. J. Urea denaturation by stronger dispersion interactions with proteins than water implies a 2-stage unfolding. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 16928−16933. (11) O’Brien, E. P.; Dima, R. I.; Brooks, B.; Thirumalai, D. Interactions between Hydrophobic and Ionic Solutes in Aqueous Guanidinium Chloride and Urea Solutions: Lessons for Protein Denaturation Mechanism. J. Am. Chem. Soc. 2007, 129, 7346−7353. (12) Robinson, D. R.; Jencks, W. P. The Effect of Compounds of the Urea-Guanidinium Class on the Activity Coefficient of Acetyltetraglycine Ethyl Ester and Related Compounds1. J. Am. Chem. Soc. 1965, 87, 2462−2470. (13) Bagchi, B. Water in Biological and Chemical Processes: From Structure and Dynamics to Function; Cambridge University Press: London, U.K., 2013. (14) Bagchi, B. Water Dynamics in the Hydration Layer around Proteins and Micelles. Chem. Rev. 2005, 105, 3197−3219. (15) Bagchi, B.; Jana, B. Solvation dynamics in dipolar liquids. Chem. Soc. Rev. 2010, 39, 1936−1954. (16) Nandi, N.; Bhattacharyya, K.; Bagchi, B. Dielectric Relaxation and Solvation Dynamics of Water in Complex Chemical and Biological Systems. Chem. Rev. 2000, 100, 2013−2046. (17) Bhattacharyya, K. Nature of biological water: a femtosecond study. Chem. Commun. 2008, 2848−2857. (18) Hammes, G. G.; Schimmel, P. R. An Investigation of WaterUrea and Water-Urea-Polyethylene Glycol Interactions. J. Am. Chem. Soc. 1967, 89, 442−446. (19) Finer, E. G.; Franks, F.; Tait, M. J. Nuclear magnetic resonance studies of aqueous urea solutions. J. Am. Chem. Soc. 1972, 94, 4424− 4429. (20) Bennion, B. J.; Daggett, V. The molecular basis for the chemical denaturation of proteins by urea. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 5142−5147. (21) Camilloni, C.; Rocco, A. G.; Eberini, I.; Gianazza, E.; Broglia, R. A.; Tiana, G. Urea and Guanidinium Chloride Denature Protein L in Different Ways in Molecular Dynamics Simulations. Biophys. J. 2008, 94, 4654−4661. (22) Rocco, A. G.; Mollica, L.; Ricchiuto, P.; Baptista, A. M.; Gianazza, E.; Eberini, I. Characterization of the Protein Unfolding Processes Induced by Urea and Temperature. Biophys. J. 2008, 94, 2241−2251. (23) Roy, S.; Bagchi, B. Comparative Study of Protein Unfolding in Aqueous Urea and Dimethyl Sulfoxide Solutions: Surface Polarity, Solvent Specificity, and Sequence of Secondary Structure Melting. J. Phys. Chem. B 2014, 118, 5691−5697. (24) Huerta-Viga, A.; Woutersen, S. Protein Denaturation with Guanidinium: A 2D-IR Study. J. Phys. Chem. Lett. 2013, 4, 3397−3401. (25) Heyda, J.; Kožíšek, M.; Bednárova, L.; Thompson, G.; Konvalinka, J.; Vondrásě k, J.; Jungwirth, P. Urea and Guanidinium Induced Denaturation of a Trp-Cage Miniprotein. J. Phys. Chem. B 2011, 115, 8910−8924. (26) Canchi, D. R.; García; Angel, E. Backbone and Side-Chain Contributions in Protein Denaturation by Urea. Biophys. J. 2011, 100, 1526−1533. (27) Lim, W. K.; Rösgen, J.; Englander, S. W. Urea, but not guanidinium, destabilizes proteins by forming hydrogen bonds to the peptide group. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 2595−2600. (28) Zangi, R.; Zhou, R.; Berne, B. J. Urea’s Action on Hydrophobic Interactions. J. Am. Chem. Soc. 2009, 131, 1535−1541. (29) Watlaufer, D. B.; Malik, S. K.; Stoller, L.; Coffin, R. L. Nonpolar Group Participation in the Denaturation of Proteins by Urea and Guanidinium Salts. Model Compound Studies. J. Am. Chem. Soc. 1964, 86, 508−514. (30) Mai Suan, L.; Binh Khanh, M. Steered Molecular Dynamics-A Promising Tool for Drug Design. Curr. Bioinf. 2012, 7, 342−351.

E(total) is the total MM perturbed QM energy of the system where both the GND1150 as a cation and phenyl side chain of PHE18 as a π-system were included in the QM system. Moreover, E(cation) is the MM perturbed QM energy of the system where GND1150 was considered in the QM system. E(π) is the MM perturbed QM energy of the system where the phenyl side chain of PHE18 was included in the QM system.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b06325. H-bond and native contact calculation, native structure of HP-36, additional figures obtained from steered molecular dynamics (SMD) simulation, umbrella samplings, and standard molecular dynamics at high temperature (400 K) (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +91 33 2473 4971. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS S.P., R.N.M., and B.J. thank the supercomputing facility CRAY at IACS for computational support. We also thank Dr. Sandipan Chakraborty, for useful discussion. S.P. thanks CSIR for awarding fellowships.

■

ABBREVIATIONS GdmCl, guanidinium chloride; MD, molecular dynamics; SMD, steered molecular dynamics; RMSD, root mean square deviation; SASA, solvent accessible surface area; PHE, phenylalanine; PMF, potential of mean force; CM, contact minima; SSM, solvent separated minima; RDF, radial distribution function; PDB, protein data bank; PME, particle mesh Ewald; WHAM, weighted histogram analysis method; QM/MM, quantum mechanics/molecular mechanics

■

REFERENCES

(1) Onuchic, J. N.; Wolynes, P. G. Theory of protein folding. Curr. Opin. Struct. Biol. 2004, 14, 70−75. (2) Wolynes, P.; Onuchic, J.; Thirumalai, D. Navigating the folding routes. Science 1995, 267, 1619−1620. (3) Bryngelson, J. D.; Onuchic, J. N.; Socci, N. D.; Wolynes, P. G. Funnels, pathways, and the energy landscape of protein folding: A synthesis. Proteins: Struct., Funct., Genet. 1995, 21, 167−195. (4) Fersht, A. R.; Itzhaki, L. S.; elMasry, N. F.; Matthews, J. M.; Otzen, D. E. Single versus parallel pathways of protein folding and fractional formation of structure in the transition state. Proc. Natl. Acad. Sci. U. S. A. 1994, 91, 10426−10429. (5) Pande, V. S.; Grosberg, A. Y.; Tanaka, T.; Rokhsar, D. S. Pathways for protein folding: is a new view needed? Curr. Opin. Struct. Biol. 1998, 8, 68−79. (6) Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How FastFolding Proteins Fold. Science 2011, 334, 517−520. (7) Tanford, C.; Kawahara, K.; Lapanje, S. Proteins in 6-M guanidine hydrochloride. Demonstration of random coil behavior. J. Biol. Chem. 1966, 241, 1921−1923. (8) Greene, R. F., Jr.; Pace, C. N. Urea and guanidine hydrochloride denaturation of ribonuclease, lysozyme, alpha-chymotrypsin, and betalactoglobulin. J. Biol. Chem. 1974, 249, 5388−5393. H


Article


(52) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656−1676. (53) Smith, L. J.; Berendsen, H. J. C.; van Gunsteren, W. F. Computer Simulation of Urea−Water Mixtures: A Test of Force Field Parameters for Use in Biomolecular Simulation. J. Phys. Chem. B 2004, 108, 1065−1071. (54) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269−6271. (55) Understanding Molecular Simulation: From Algorithms to Applications; Academic Press, Inc., 1996. (56) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (57) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (58) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (59) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (60) Godawat, R.; Jamadagni, S. N.; Garde, S. Unfolding of Hydrophobic Polymers in Guanidinium Chloride Solutions. J. Phys. Chem. B 2010, 114, 2246−2254. (61) Shimizu, S.; Chan, H. S. Origins of protein denatured state compactness and hydrophobic clustering in aqueous urea: Inferences from nonpolar potentials of mean force. Proteins: Struct., Funct., Genet. 2002, 49, 560−566. (62) Ikeguchi, M.; Nakamura, S.; Shimizu, K. Molecular Dynamics Study on Hydrophobic Effects in Aqueous Urea Solutions. J. Am. Chem. Soc. 2001, 123, 677−682. (63) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. THE weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011−1021. (64) Turner, A. J.; Moliner, V.; Williams, I. H. Transition-state structural refinement with GRACE and CHARMM: Flexible QM/MM modelling for lactate dehydrogenase. Phys. Chem. Chem. Phys. 1999, 1, 1323−1331. (65) Field, M. J.; Albe, M.; Bret, C.; Proust-De Martin, F.; Thomas, A. The dynamo library for molecular simulations using hybrid quantum mechanical and molecular mechanical potentials. J. Comput. Chem. 2000, 21, 1088−1100.

(31) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free energy calculation from steered molecular dynamics simulations using Jarzynski’s equality. J. Chem. Phys. 2003, 119, 3559−3566. (32) Park, S.; Schulten, K. Calculating potentials of mean force from steered molecular dynamics simulations. J. Chem. Phys. 2004, 120, 5946−5961. (33) Grubmüller, H.; Heymann, B.; Tavan, P. Ligand Binding: Molecular Mechanics Calculation of the Streptavidin-Biotin Rupture Force. Science 1996, 271, 997−999. (34) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187−199. (35) Sasikala, W. D.; Mukherjee, A. Single Water Entropy: Hydrophobic Crossover and Application to Drug Binding. J. Phys. Chem. B 2014, 118, 10553−10564. (36) Mukherjee, A. Entropy Balance in the Intercalation Process of an Anti-Cancer Drug Daunomycin. J. Phys. Chem. Lett. 2011, 2, 3021− 3026. (37) Dasgupta, A.; Udgaonkar, J. B.; Das, P. Multistage Unfolding of an SH3 Domain: An Initial Urea-Filled Dry Molten Globule Precedes a Wet Molten Globule with Non-Native Structure. J. Phys. Chem. B 2014, 118, 6380−6392. (38) Evans, R.; Wilding, N. B. Quantifying Density Fluctuations in Water at a Hydrophobic Surface: Evidence for Critical Drying. Phys. Rev. Lett. 2015, 115, 016103. (39) Chandler, D. Physical chemistry: Oil on troubled waters. Nature 2007, 445, 831−832. (40) Mukherjee, A.; Bhimalapuram, P.; Bagchi, B. Orientationdependent potential of mean force for protein folding. J. Chem. Phys. 2005, 123, 014901. (41) Buchete, N. V.; Straub, J. E.; Thirumalai, D. Orientationdependent coarse-grained potentials derived by statistical analysis of molecular structural databases. Polymer 2004, 45, 597−608. (42) Vaitheeswaran, S.; Chen, J.; Thirumalai, D. Hydrophobic and Ionic-Interactions in Bulk and Confined Water with Implications for Collapse and Folding of Proteins. J. Stat. Phys. 2011, 145, 276−292. (43) Shi, X.; Shao, C.; Zhang, X.; Zambonelli, C.; Redfield, A. G.; Head, J. F.; Seaton, B. A.; Roberts, M. F. Modulation of Bacillus thuringiensis phosphatidylinositol-specific phospholipase C activity by mutations in the putative dimerization interface. J. Biol. Chem. 2009, 284, 15607−15618. (44) Shao, C.; Shi, X.; Wehbi, H.; Zambonelli, C.; Head, J. F.; Seaton, B. A.; Roberts, M. F. Dimer structure of an interfacially impaired phosphatidylinositol-specific phospholipase C. J. Biol. Chem. 2007, 282, 9228−9235. (45) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (46) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent molecular orbital methods 25. Supplementary functions for Gaussian basis sets. J. Chem. Phys. 1984, 80, 3265−3269. (47) Jorgensen, W. L.; Tirado-Rives, J. The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (48) Davis, M. R.; Dougherty, D. A. Cation-[small pi] interactions: computational analyses of the aromatic box motif and the fluorination strategy for experimental evaluation. Phys. Chem. Chem. Phys. 2015, 17, 29262−29270. (49) McKnight, J. C.; Doering, D. S.; Matsudaira, P. T.; Kim, P. S. A Thermostable 35-Residue Subdomain within Villin Headpiece. J. Mol. Biol. 1996, 260, 126−134. (50) McKnight, C. J.; Matsudaira, P. T.; Kim, P. S. NMR structure of the 35-residue villin headpiece subdomain. Nat. Struct. Biol. 1997, 4, 180−184. (51) Gromacs User Manual; www.gromacs.org. I


Destabilization of Hydrophobic Core of Chicken Villin Headpiece in

Recommend Documents