Free Energy Landscapes of Alanine Oligopeptides in Rigid-Body and

Article pubs.acs.org/JPCB

Free Energy Landscapes of Alanine Oligopeptides in Rigid-Body and Hybrid Water Models Divya Nayar and Charusita Chakravarty*

Downloaded by TEXAS A&M INTL UNIV on August 31, 2015 | http://pubs.acs.org Publication Date (Web): July 16, 2015 | doi: 10.1021/acs.jpcb.5b02937

Department of Chemistry, Indian Institute of Technology-Delhi, New Delhi 110016, India ABSTRACT: Replica exchange molecular dynamics is used to study the effect of different rigid-body (mTIP3P, TIP4P, SPC/ E) and hybrid (H1.56, H3.00) water models on the conformational free energy landscape of the alanine oligopeptides (acAnme and acA5nme), in conjunction with the CHARMM22 force field. The free energy landscape is mapped out as a function of the Ramachandran angles. In addition, various secondary structure metrics, solvation shell properties, and the number of peptide−solvent hydrogen bonds are monitored. Alanine dipeptide is found to have similar free energy landscapes in different solvent models, an insensitivity which may be due to the absence of possibilities for forming i−(i + 4) or i−(i + 3) intrapeptide hydrogen bonds. The pentapeptide, acA5nme, where there are three intrapeptide backbone hydrogen bonds, shows a conformational free energy landscape with a much greater degree of sensitivity to the choice of solvent model, though the three rigid-body water models differ only quantitatively. The pentapeptide prefers nonhelical, non-native PPII and β-sheet populations as the solvent is changed from SPC/E to the less tetrahedral liquid (H1.56) to an LJ-like liquid (H3.00). The pentapeptide conformational order metrics indicate a preference for open, solvent-exposed, non-native structures in hybrid solvent models at all temperatures of study. The possible correlations between the properties of solvent models and secondary structure preferences of alanine oligopeptides are discussed, and the competition between intrapeptide, peptide−solvent, and solvent−solvent hydrogen bonding is shown to be crucial in the relative free energies of different conformers. fields. Examples of the latter approach include incorporating the statistical mechanics of solvation, including the hydrophobic effect,19−22 and better models of electrostatics that include many-body inductive effects.23 As computational power increases, more complex parametrization of potentials incorporating such theories become feasible. In this context, considerable attention has been focused on the role of water as a solvent that controls the interplay of hydrophilic and hydrophobic hydration, shaping the secondary structure preferences of proteins and peptides. Solvation of small hydrophobic solutes is known to be crucially dependent on the equation of state of bulk solvent.24−26 Though side chains of amino acids fall in the small hydrophobe regime, protein hydration shows a more complex dependence on the solvent model, since the nanoscale size of proteins, the patchy nature of hydrophobicity, the specific interactions between amino acid side chains, and the collective rearrangements associated with hydrophobic core formation are additional factors that must be accounted for.13,21,22,25,27−33 Hydrophilic interactions in peptides originate from electrostatic and inductive interactions

1. INTRODUCTION The free energy landscape of a protein determines the thermodynamic and kinetic aspects of folding and depends on a complex interplay of energetic and entropic factors.1−3 In order for atomistic computer simulations to provide an accurate as well as efficient model for understanding the atomic level changes underlying these processes, standard biomolecular force fields model the dependence of the configurational energy on the atomic position coordinates using relatively simple functional forms.4−8 Force-field parameters are obtained by fitting to thermodynamic, crystallographic, and solution NMR data at or near standard temperature and pressure conditions. This parametrization strategy implies that simulations can sample the free energy landscape accurately in the neighborhood of the native structure. However, the predictive value of the force fields is significantly reduced as deviations from native structure or ambient conditions increase, such as when modeling pressure effects,9 cold denaturation,10 or unfolded ensembles of proteins.9,11,12 One strategy for improving the predictive value of force fields is to enlarge the structural information base used for parametrization, such as NMR information on residual structure in unfolded states and in intrinsically disordered proteins.11,13−18 A complementary approach is to develop an understanding of the physics of inter- and intramolecular forces operating in solvated proteins that can eventually be built into the parametrization of force © 2015 American Chemical Society

Special Issue: Biman Bagchi Festschrift Received: March 27, 2015 Revised: July 1, 2015 Published: July 1, 2015 11106

DOI: 10.1021/acs.jpcb.5b02937 J. Phys. Chem. B 2015, 119, 11106−11120

Article


The Journal of Physical Chemistry B

contrast, the two other three-site models, mTIP3P and TIP3P, provide a very poor description of a water model with TMD values close to 182 K at 1 atm pressure.37,40 Standard protein force fields are parametrized for a given water model, e.g., CHARMM, AMBER, and GROMOS using mTIP3P, TIP3P, and SPC, respectively.6−8 It is only relatively recently that the consequences of varying solvent models, for a given choice of protein force fields, are beginning to be explored. The differences between the water models are quantitative, rather than qualitative, at the level of local hydration shell structure, energetics, and dynamics.41,50−52 These quantitative differences in local hydration behavior translate into significant differences in the secondary structure preferences of biomolecules.9,13,33,53,54 The pressure-induced stabilization of isolated α-helices has been recently examined for the 15-residue Ac(AAQAA)3-NH2 peptide in TIP4P/2005 and TIP3P water models with the Amber03w force field.9 While the TIP4P/2005 model with an accurate equation of state is able to reproduce the pressure-induced stability, the TIP3P model is not. The best model for predicting properties of bulk water need not, however, necessarily provide the best description of aqueous solvation, since the physics underlying solvation of hydrophobic and hydrophilic solutes is substantially different, and can result in a sequence-dependent response of the folding landscapes to changes in the solvent.33,55,56 In the case of the 16-residue βhairpin fragment of 2GB1, the TIP4P model favors a more compact hydrophobic core in the unfolded ensemble compared to the mTIP3P solvent.33 In the case of the Trp-cage miniprotein, the TIP4P-Ew model lowers the free energy for unfolding in comparison to the TIP3P due to differences in solvation of polar groups.56 A recent comparative study of the Trp-cage, GB1 β-hairpin, and AK16 peptide indicates that negative pressure destabilizes the Trp-cage but stabilizes the other two.55 With the increasing computational feasibility of free energy landscape studies, biomolecular simulations have been carried out in modified rigid-body water models in which the magnitude of dispersion−repulsion interactions relative to the electrostatic interactions is increased. These models are useful for understanding the balance between hydrophilic and hydrophobic interactions that underlies bimolecular solvation. For example, hybrid water models (Hn.00) based on modification of the SPC/E water model have been used to study DNA dodecamer and a polyalanine-based peptide, though the free energy landscape was not mapped out.57,58 Pande et al. have modified the LJ parameters of the TIP3P model to study solvation of an alanine-based peptide.59 Shaw et al. have modified the TIP4P water model to examine the role of dispersive interactions in modeling of the disordered states such as those of intrinsically disordered proteins.60 In this paper, we have chosen to use the SPC/E based hybrid models, since the bulk liquid state anomalies and hydrogen-bonding properties of these hybrid water solvents have been well-characterized in previous studies.47,48,61,62 Alanine polypeptides are useful systems to study the effects of solvent on the secondary structure preferences of peptides, since the homopolymeric structure eliminates complexities due to patchy hydrophobicity and residue-specific interactions.9,13,30,31,63−69 Long, polyalanine chains are known to form αR-helical structures where the intrapeptide hydrogen bonds between the carbonyl and amide groups of the ith and (i + 4)th residues hold the helical structure together.1 This study focuses on alanine dipeptide (CH3CO−NHCH(CH3)CO−

between charged and polar amino acid side chains, as well as hydrogen-bonding interactions between backbone CO and NH groups and water. Additional electrostatic effects, such as formation of salt bridges and electrostriction, can play a crucial role for certain proteins.3,34,35 The simulation studies presented here are directed at obtaining greater quantitative insight into the interplay between electrostatic and dispersion− repulsion contributions to peptide−water interactions that shape the free energy landscape. We focus on alanine oligopeptides, so that complexities due to variations in sidechain water interactions are avoided, and we can focus on understanding the competition between intrapeptide, peptide− solvent, and solvent−solvent hydrogen bonds. Biomolecular force fields are constructed from a protein force field and an aqueous solvent model. Water is typically modeled as a liquid bound by rigid-body, effective pair potentials though polarizable models are becoming increasingly more viable computationally.6−8,13,36 In the rigid-body water models, a fixed molecular geometry is assumed and repulsion− dispersion interactions are modeled using a single LennardJones site located on the oxygen atom. Hydrogen bonding and electrostatic interactions are parametrized using a distribution of partial charges. Inductive effects are ignored, except indirectly using an enhanced effective dipole appropriate for bulk, rather than gas- phase, water. Table 1 gives the parameters Table 1. Potential Energy Parameters of Rigid-Body Water Models and Hybrid Solvent Models37,45−48 a ϵOO (kJ/mol) σOO (Å) ϵHH (kJ/mol) σHH (Å) ϵOH (kJ/mol) σOH (Å) rOH (Å) rOM (Å) ∠HOH (deg) qO (e) qH (e) qM (e)

mTIP3P

TIP4P

SPC/E

H1.56

H3.00

0.636 3.1506 0.1924 0.4000 0.3497 1.7755 0.9572

0.649 3.1540

0.649 3.1656

1.012 3.1656

1.949 3.1656

0.9572 0.15 104.52

1.000

1.000

1.000

109.47 −0.8472 0.4238

109.47 −0.8472 0.4238

109.47 −0.8472 0.4238

104.52 −0.8340 0.417

0.52 −1.04

a

The Lennard-Jones sites are associated with a well depth parameter, ϵ, and a length scale parameter, σ, where the subscripts OO, HH, and OH denote distinct pair interactions between atomic sites. The charges on oxygen atoms, hydrogen atoms, and mass-less sites are denoted by qO, qH, and qM, respectively. The molecular geometry and site distribution are determined by the fixed length parameters, rOH and rOM, and H−O−H bond angle.

of commonly used rigid-body water models. These water models yield very similar results for the properties of bulk water under ambient conditions, but they differ very significantly in the temperature regime for liquid state anomalies and phase diagrams.37−44 The four-site models, TIP4P and TIP4P/2005, provide a reasonable description of the equation of state of water and behave as strongly tetrahedral liquids at room temperature. The TIP4P/2005 model is currently probably the most reliable in terms of modeling the thermodynamic properties of bulk water with a temperature of maximum density (TMD) of 280 K at 1 atm pressure.38,49 The closely related TIP4P model has a lower TMD of 253 K at 1 atm. The three-site SPC/E water model has a TMD of 241 K at 1 atm. In 11107


Article



Figure 1. Ramachandran plot of (a) acAnme at 300 K and (b) acA5nme at 298.43 K in mTIP3P. The different types of configurations that the peptides can attain are marked in the figure. The data has been plotted for 6500 and 10 000 configurations (picked at regular intervals) out of a total of 26 000 and 30 000 dumped configurations over a trajectory with a run length of 26 ns for acAnme and 60 ns for acA5nme, respectively, as obtained from REMD simulation. The regions considered for computing average population fractions of the peptides are indicated by the rectangular regions and are also described in Table 6 for acAnme.

conformation).73,74 The second quadrant of the Ramachandran plot of the dipeptide, also known as the β-basin, is shared by three populations: C5, C7eq, and PPII.29,67,74−76 Previous free energy landscape studies with CHARMM22 and CHARMM27 force fields show that the C5 and PPII conformers are not wellseparated by a free energy barrier in typical rigid-body water models and, therefore, are considered as belonging to a single free energy basin called the β-basin.67,75 Therefore, we include the C5, C7eq, and PPII populations in the (C5 + C7eq) population, without any further subdivision of populations, as shown in Figure 1. The pentapeptide can form three intrapeptide backbone hydrogen bonds of which two require participation of CO and NH groups of the capping moieties. The pentapeptide occupies the αR, β-sheet, and PPII conformations in the Ramachandran plot, which we have considered separately and marked in Figure 1. We use the CHARMM22 force field for both acAnme and acA5nme, in keeping with previous studies, since it has identical protein parameters to CHARMM27.29,67 CHARMM22 is

NHCH3, denoted here by acAnme) which is the basic building block of polyalanine chains and alanine pentapeptide (CH3CO−(NHCH(CH3)CO)5−NHCH3, denoted here by acA5nme) which is just large enough for incipient helix formation with a single αR-helical turn. The capping of the free NH2 and CO2H by acetyl (ac) and N-methylamide groups (nme), respectively, ensures that additional charged interactions due to zwitterionic forms are absent. The small size of the peptides implies that they can sample a much wider range of conformation space under ambient conditions, compared to longer chains because of reduced steric hindrance. Experimentally, short alanine peptides show a preference for the PPII configuration in aqueous medium compared to longer chains.32,65,70−72 Dialanine has two peptide bonds and one intrapeptide hydrogen bond, which can be formed either between the carbonyl and amide of the alanine residue (which stabilizes the five-membered ring conformation) or between the carbonyl of the acetyl group and amide of the N-methylamide group (which results in formation of a seven-membered ring 11108


Article

The Journal of Physical Chemistry B Table 2. Some Properties of the Bulk Solvents at 300 K and 1.0 g cm−3 37,45−48 a model

Uconf

Uvdw

Uelec

⟨qtet⟩

EHB

μ (D)

QT (D-Å)

μ/QT (Å−1)

ϵr

mTIP3P TIP4P SPC/E H1.56 H3.00

−40.8 −41.1 −46.3 −41.6 −41.7

4.5 7.6 9.0 4.0 −6.8

−45.3 −48.7 −55.3 −45.6 −34.9

0.56 0.63 0.63 0.56 0.49

−31.56 −30.03 −37.10 −32.40 −25.50

2.350 2.177 2.350 2.350 2.350

1.721 2.147 2.035 2.035 2.035

1.363 1.014 1.155 1.155 1.155

94.0 50.0 68.0 75.9 90.3

a


The total, electrostatic, and Lennard-Jones contributions to the potential energy per molecules are denoted by Uconf, Uelec, and Uvdw and given in kJ mol−1. Hydrogen bonding strength (EHB), tetrahedral order (⟨qtet⟩), dipole moment (μ), quadrupole moment (QT), ratio of dipole−quadrupole moments (μ/QT), and relative dielectric constant (ϵr) are also reported for different solvent models. The H-bond strength has been computed for the mTIP3P and TIP4P water models following the procedure in ref 48 and corresponds to the dimer energy obtained for the O−O distance corresponding to the first peak of the O−O radial distribution function.

biased toward α-helical conformations,65,77 which has been corrrected in the more recent CHARMM36 force field.78 Since our purpose is to map the free energy landscape as a function of the solvent model, for a given protein force field, as well as to compare with previous studies, we use the CHARMM22 force field. Three rigid-body water models discussed above (mTIP3P, TIP4P, SPC/E) and two hybrid water models (H1.56, H3.00) are used as solvents (see Table 1 for parameters). As discussed above, in the biologically relevant temperature regime (approximately 275−325 K), TIP4P and SPC/E behave as tetrahedral liquids with anomalous properties similar to those of bulk water, while mTIP3P must be regarded as essentially a simple liquid.40,41 The hybrid models are designed to interpolate between the Lennard-Jones (LJ) fluid and SPC/E water model, with LJ contribution to the solvent−solvent effective pair interaction enhanced by a factor of λ relative to the SPC/E model.47,48,61,62 While the H1.56 hybrid model is a weakly tetrahedral liquid with residual water-like anomalies, H3.00 behaves essentially as a simple liquid.47 In Table 2, we summarize some of the properties of the bulk model liquids that are relevant to solvation at 300 K and 1 g cm−3.47,48,61,62 For all of the solvents, including the H3.00 model, the configurational energy is dominated by the electrostatic contribution. However, the H3.00 liquid is the only one where the LJ contribution is attractive and this results in distinct changes in structural properties, such as the radial distribution function and the tetrahedral order parameter. As indicated by the ⟨qtet⟩ values, the SPC/E and TIP4P models have relatively high and very similar degrees of tetrahedral order. The mTIP3P and H1.56 liquids have very similar, weak degrees of tetrahedrality. The H3.00 model has very weak tetrahedrality, and its local structure is dominated by dispersion−repulsion interactions. The hydrogen-bonding strengths are strongest for SPC/E and weakest for H3.00; however, hydrogen bonds in H3.00 are about 70% of the strength of SPC/E hydrogen bonds. The dipole moments of all five models are very similar, and the dipole−quadrupole moments ratio is distinctly different only for the mTIP3P model. The variation in local structure and hydrogen-bonding strengths in these model solvents are expected to provide an understanding of the interplay of intrapeptide, solvent−peptide, and solvent−solvent hydrogen bonding in shaping the free energy landscape. Since alanine dipeptide has a well-shielded intrapeptide hydrogen bond while alanine pentapeptide has more exposed backbone hydrogen bonds, we expect different sensitivities to solvent models. We carry out replica exchange molecular dynamics (REMD) simulations of alanine oligopeptides in these five different solvent models to map out the conformational free energy landscape.

The organization of the paper is as follows. Section 2 describes computational details. Section 3 discusses our results, focusing on secondary structure metrics, hydration shell structure, and the competition between peptide−water and water−water hydrogen bonding. Conclusions are summarized in section 4. We note that the terms water−water and peptide− water in this paper refer to the interaction between solvent molecules and between peptide and solvent, respectively, irrespective of whether the solvent molecules belong to rigidbody or the hybrid class of models. Similarly, the term hydration refers to solvation due to solvent molecules of either class of the models.

2. COMPUTATIONAL DETAILS Replica exchange molecular dynamics simulations were carried out for two sets of systems: (a) alanine dipeptide in five different solvent models and (b) penta-alanine in five different solvent models. All simulations were done using the CHARMM22 force field and NAMD package (either version 2.7 or 2.9, depending on the availability of GPUs). Since the CHARMM22 force field is parametrized with the mTIP3P solvent model, the other solvent models are modeled by modifying the Lennard-Jones interaction energy parameters, bond lengths, and bond angles in the parameter file and partial charges on the atoms in the topology file. The potential energy parameters for the three rigid-body water and hybrid models are given in Table 1. The cutoff for nonbonded LJ interactions was chosen to be 9 Å. The particle mesh Ewald method was used for computing long-range electrostatic interactions. System Setup. The initial configurations of alanine dipeptide were prepared using the Molefacture plugin in VMD which allows editing or building of molecules. The Nterminus and C-terminus of the peptides were capped using acetyl and N-methylamide groups, respectively. The structure files for the two peptides were constructed using the psfgen plugin in VMD. The dipeptide (acAnme) and pentapeptide (acA5nme) were then solvated by 544 and 1050 solvent molecules (mTIP3P, TIP4P, SPC/E, H1.56, H3.00), respectively, using the Solvate plugin of VMD in a cubical box (see Table 3). The details of the system setup for acA5nme and acAnme were taken from refs 79−81. Molecular Dynamics Simulations. Energy minimization for each solvated peptide system was performed using the conjugate gradient method. Although it is usually not essential to heat the system to the required temperatures but since the system was initialized with a suitable energy minimized configuration, slow heating of the system in the NPT ensemble was carried out from 0 K to the desired temperature for equilibration, i.e., 300 K for alanine dipeptide and 250 K for 11109


Article


temperature-wise for further analysis. Standard errors on measured quantities were calculated by splitting the trajectory into 26−30 parts each containing 1000 equi-spaced configurations. All the analysis was done using standard or in-house VMD Tcl scripts and in-house Fortran codes. The convergence of the REMD simulations was checked by computing the replica mixing parameter82 and by using Bennett’s histogram overlap method.83,84 The replica mixing parameter (m(T)) allows one to verify whether sufficient exchange between replicas has taken place, by quantifying the time spent by each replica at a given temperature, and is defined as

Table 3. Details of the Solvated Peptide System with Equilibrium Density (ρ) and Box Length (l) Obtained after NPT Equilibration of acAnme and acA5nme at 300 and 250 K, Respectively, at 1 atm Pressure, in Different Solvent Systemsa peptide acAnme

acA5nme

solvent model

ρ (g cm−3)

l (Å)

A

mTIP3P TIP4P SPC/E H1.56 H3.00 mTIP3P TIP4P SPC/E H1.56 H3.00

0.993 0.986 0.863 0.826 0.887 1.059 1.041 1.007 0.953 0.949

25.52 25.58 26.75 27.14 26.5 31.19 31.37 31.71 32.31 32.34

0.18−0.32 0.17−0.38 0.17−0.36 0.20−0.32 0.22−0.38 0.29−0.38 0.23−0.33 0.23−0.31 0.31−0.38 0.39−0.45

m (T ) = 1 −

∑i ti

(1)

where ti is the total time spent by replica i at temperature T as discussed in ref 82. The ideal value of m(T) is obtained by assuming that each replica spends equal time at each temperature and is given by (1 − 1/√N), where N is the number of replicas of the system, i.e., 24 and 20 for acAnme and acA5nme, respectively.82 In all of our simulations, the values of m(T) obtained at all temperatures are within an error of 5% of the expected ideal values of 0.796 and 0.776 for acAnme and acA5nme, respectively, implying that our simulations are well-converged. Bennett’s histogram overlap method checks whether temperature sampling in the simulations obeys the following Boltzmann distribution

a


∑i ti 2

The following densities of the systems were used to carry out REMD simulations. Also reported are the acceptance ratios (A) of exchange moves for REMD simulation of each system.

alanine pentapeptide at a rate of 0.2 K ps−1. The damping parameters for the Langevin thermostat were 10 and 5 ps−1 for acAnme and acA5nme, respectively. A Nosé−Hoover barostat with an oscillation period of 100 fs was used to maintain a pressure of 1 atm. The time step was fixed at 2 fs. Equilibration and production runs were carried out in the NPT and NVT ensembles, respectively, with the same thermostat and barostat parameters as used for heating the system. Production run lengths for the dipeptide were 16 ns at 300 K, while, for the pentapeptide, they were 20 ns long at 250 K. REMD Simulations. In order to carry out REMD simulations, the initial configuration of the system was picked from the MD production trajectory and a suitable number of replicas were initialized over a specified temperature range. In the case of the dipeptide, the initial configuration had the dihedral angles lying in the αR region of the Ramachandran plot and 24 replicas with temperatures ranging from 300 to 600 K were used. In the case of the pentapeptide, a fully helical configuration in which all five dihedral Ramachandran angles were in the αR helical region of the Ramachandran plot (−62°, −41°) and the peptide formed three backbone hydrogen bonds was used as the initial conformation for REMD simulations and 20 replicas with temperatures ranging from 250 to 350 K were used for the pentapeptide. The NAMD package uses the geometric distribution method to generate the set of temperatures of replicas Ti. The box length along with the equilibrium density of each system, used for REMD simulations, are given in Table 3. The parameters for the REMD simulations for both oligopeptides are summarized in Table 4. During the simulation, the trajectory of each replica was stored, as it traversed the temperature space as a function of time. The trajectories, therefore, were sorted

⎡ H(E , Ti + 1) ⎤ ⎛ 1 1 ⎞ ln⎢ − ⎟E + constant ⎥=⎜ kBTi + 1 ⎠ ⎣ H(E , Ti ) ⎦ ⎝ kBTi

(2)

where H(E, Ti) is the potential energy histogram of the system at temperature Ti. Our simulations obey the above linear relationship, with the slope of the plots within an error of 5% of the expected value, for both acAnme and acA5nme simulations.

3. RESULTS AND DISCUSSION 3.1. Free Energy Landscape as a Function of Ramachandran Angles. Figure 1 shows the occupancy of the Ramachandran plot for acAnme and acA5nme solvated in mTIP3P water, using REMD simulations, to clearly define the secondary structure regions of interest with the solvent model appropriate for the CHARMM force field. The conformational flexibility of acAnme implies that it can occupy regions of the Ramachandran plot not accessible to larger peptides or proteins. Figure 1a shows the regions of the Ramachandran plot (ϕ, ψ) occupied by six different conformational structures of the dipeptide, defined by the (ϕ, ψ) dihedral angles as (i) extended β-sheet structure, β2 (−137.8°, 22.9°), (ii) right-handed α-helix, αR (−60°, −45°), (iii) lefthanded α-helix, αL (63.5°, 34.8°), (iv) 5-membered ring

Table 4. Canonical Replica Exchange Molecular Dynamics Simulation Details for the System of acAnme and acA5nme in Different Solvent Modelsa system

Nwat

Nrep

Tmin−Tmax (K)

trep (ns)

ttot (μs)

τLang (ps−1)

τ (fs)

ν (ps)

acAnme acA5nme

544 1050

24 20

300−600 250−350

26 60

0.624 1.200

10 5

2 2

1 2

a The number of replicas is denoted by Nrep. The temperature range of simulations is given by the minimum (Tmin) and maximum (Tmax) temperatures. The simulation run length for each replica (trep), the total simulation run length (ttot), the time step of simulation (τ), the damping coefficient of Langevin thermostat (τLang), and the time interval at which the replicas are exchanged (ν) are reported.

11110


Article



Figure 2. Contour plots of Ramachandran angles of acAnme and acA5nme on a free energy scale, where the units of free energy are J mol−1. The first column compares the Ramachandran plot of acAnme in different solvent models at 300 K. The second and third columns compare the Ramachandran plot of acA5nme in different solvent models at 250 and 298.43 K, respectively.

conformation, C5 (−158.4°, 161.6°), arising due to hydrogen bonding between the carbonyl and amide groups of the alanine residue, (v) 7-membered ring conformation, C7eq (−82.8°, 77.9°), formed due to hydrogen bonding between the carbonyl of the acetyl group and the amide of the N-methylamide group, with the side-chain CH3 group equatorial to the ring, and (vi) 7-membered ring analogous to C7ax (74.4°, −64.1°) with the side-chain CH3 group axial to the ring.74,81 The values of (ϕ, ψ) given in parentheses in the above definitions correspond to the free energy minima in the Ramachandran plot.74 The C7ax conformation occupies the fourth quadrant of the Ramachandran plot, which is a forbidden region for large peptides and proteins.

Figure 1b shows an occupancy plot for acA5nme in the Ramachandran plane. Only a few conformations occupy the αLhelical region. There is significant occupancy of the polyproline II, β-sheet, and αR-helical regions of the Ramachandran plot. In addition, unlike in acAnme, there is minimal occupancy of the C7ax region but significant population in the region with ϕ < −50° and ψ < −50°, which is unpopulated in the case of acAnme. Figure 2 shows the conformational free energy landscape of acAnme and acA5nme in different solvent models as a function of the Ramachandran (ϕ, ψ) angles, i.e., F(ϕ, ψ) = −kBT ln P(ϕ, ψ), where P(ϕ, ψ) is the probability of observing dihedral angle values associated with a peptide bond of any amino acid residue. Results for the dipeptide at 300 K are shown in the first 11111


Article



Table 7 gives the average population fractions of different configurations of alanine dipeptide at 300 K in different solvent

column and for the pentapeptide at 250 and 298.43 K in the second and third columns of the figure, respectively. The acAnme peptide attains all possible conformations in all the solvent models, and there are no qualitative differences in the stability and propensity of different conformations of acAnme between rigid-body water models and the hybrid water models. The conformational free energy landscape of acA5nme, on the other hand, shows that, even though the overall structure of the free energy landscape is not qualitatively altered on changing the solvent model, the relative propensities and stability of populations of different conformers vary in different solvent models at both 250 and 298.43 K. The αR-helical conformation is found to be the most prominent and stable population of this peptide in all rigid-body water models. The stability of this population, however, is seen to decrease as the solvent model is changed from SPC/E to H1.56 to the more LJ-like H3.00 model at both 250 and 298.43 K. 3.2. Secondary Structure Populations. To quantify the relative stability of different conformers of peptides as a function of the temperature and solvent model, we compare the trajectory-averaged population fractions of occupancy of different regions of the Ramachandran plot. For calculating the average fractions of populations of acAnme, configurations of acAnme lying in the (ϕ, ψ) regions described in Table 6

Table 7. Average Fractions of Populations of acAnme Occupying Different Regions of the Ramachandran Plot in Different Solvent Models at 300 K

conformations

ϕ (deg)

ψ (deg)

−139 −119 −57 −49 −57 −83 −78

135 113 −47 −26 −70 158 149

Table 6. Regions in the Ramachandran Plot Used for Calculating Fractions of Different Populations of acAnme (Indicated by Rectangles with Different Colors in Figure 1a) conformations

ϕ (deg)

ψ (deg)

C7eq + C5

−180 < ϕ < 0 −180 < ϕ < 0 130 < ϕ < 180 −180 < ϕ < 0 130 < ϕ < 180 0 < ϕ < 130 0 < ϕ < 130 0 < ϕ < 130

70 < ψ < 180 −180 < ψ < −110 70 < ψ < 180 −110 < ψ < 70 −110 < ψ < 70 −60 < ψ < 90 −180 < ψ < −60 90 < ψ < 180

αR + β2 αL C7ax

C7eq + C5

αR + β2

αL

C7ax

mTIP3P TIP4P SPC/E H1.56 H3.00

0.455 0.529 0.440 0.473 0.437

0.418 0.340 0.405 0.490 0.511

0.081 0.086 0.094 0.035 0.035

0.045 0.045 0.061 0.002 0.017

models. The αL and C7ax regions in all the solvent models have an occupancy of less than 10%, compared to other secondary structures. The adjacent C7eq + C5 regions have a maximum occupancy of ∼53% in TIP4P and a minimum occupancy of ∼44% in SPC/E as well as H3.00 solvent. The occupancy of the αR + β2 region varies from a maximum of 51% in H3.00 to a minimum of 40% in SPC/E. The αR and β2 populations of the dipeptide increase as the model is changed from SPC/E to H1.56 to H3.00, i.e., with increasing strength of repulsion− dispersion interactions, possibly because of the compact nature of the αR structure. Figure 3 shows the temperature-dependent fractions of the population of alanine pentapeptide having αR-helical (Fhel), βsheet (Fβ‑sheet), and polyproline II (FPPII) conformations in different solvent models. In all of the solvent models, Fβ‑sheet and FPPII are much less than Fhel, even though Fhel decreases with an increase in temperature in all solvent models. Fhel is approximately 50% in all rigid-body water models at 300 K and the temperature dependence of the helical fraction in rigidbody water models is identical above 300 K. The Fhel fraction shows significant differences in all solvent models at temperatures below 300 K and increases in the order H3.00 < H1.56 mTIP3P ≈ H1.56 > H3.00. Such a trend in qtet is similar to that observed for the bulk system of these solvent models. mTIP3P and H1.56 show a very similar tetrahedral parameter profile in the solvation shells of both peptides. With H3.00 being a poor tetrahedral liquid, it shows the lowest qtet values in the hydration shells of a peptide. 3.5. Intrapeptide and Peptide−Solvent Hydrogen Bonds. The average numbers of hydrogen bonds have been

Figure 4. Temperature dependence of (a) average radius of gyration, ⟨Rg⟩, (b) average root-mean-square deviation (⟨RMSD⟩), (c) average root-mean-square fluctuation (⟨RMSF⟩), and (d) average solvent accessible surface area (⟨SASA⟩) of acA5nme as a function of temperature in different solvent models. Each inset shows the temperature dependence of the respective order metric of acAnme in different solvent models. The line colors and point symbols in the inset figures represent the same solvent model as those shown in the main figure.

size and fluctuation metrics of the pentapeptide, the RMSF shows the highest variation with temperature, especially in the case of the SPC/E model where the peptide shows 61% variation in RMSF values as the temperature is changed from 250 to 350 K. As summarized in Table 8, the trends in Rg, RMSD, RMSF, and SASA for different solvent models at a given temperature are anticorrelated with the helical propensity of the peptide, given by Fhel. For example, the peptide shows the highest Fhel in SPC/E where the conformations are 11114


Article


Table 8. Trends Observed in Different Conformational Order Metrics and Hydration Shell Metrics of acA5nme in Different Solvent Modelsa observable

ranking

Fhel Fβ‑sheet, FPPII ⟨Rg⟩, ⟨RMSD⟩, ⟨RMSF⟩, ⟨SASA⟩ ⟨nfirst wat ⟩ ⟨qfirst tet ⟩ ⟨NPep HB ⟩ ⟨NPW HB ⟩ ⟨NWW HB ⟩ EHB

H3.00 < H1.56 < mTIP3P ≈ TIP4P < SPC/E SPC/E < mTIP3P ≈ TIP4P < H.156 < H3.00 SPC/E < mTIP3P ≈ TIP4P < H.156 < H3.00 H1.56 < H3.00 < SPC/E < TIP4P < mTIP3P H3.00 < mTIP3P ≈ H1.56 < SPC/E ≈ TIP4P H3.00 < H1.56 < mTIP3P ≈ TIP4P < SPC/E mTIP3P < H1.56 ≈ SPC/E < H3.00 < TIP4P H3.00 < H1.56 ≈ mTIP3P < SPC/E < TIP4P H3.00 < TIP4P < mTIP3P < H1.56 < SPC/E


a The fractions of helical, β-sheet, and PPII populations are denoted by Fhel, Fβ‑sheet, and FPPII, respectively. The trends in average radius of gyration (⟨Rg⟩), average solvent accessible surface area (⟨SASA⟩), average root-mean-square deviation (⟨RMSD⟩), and average root-mean-square fluctuation first (⟨RMSF⟩) are reported. The trends in average number (⟨nfirst wat ⟩) and average tetrahedral order (⟨qtet ⟩) of solvent molecules in the first hydration WW ⟩), peptide−water (⟨NPW shell as well as average number of hydrogen bonds formed between peptide atoms (⟨NPep HB HB ⟩), and water−water (⟨NHB ⟩) are also reported. EHB is the hydrogen bond strength of the solvent dimer as reported in Table 2.

H-bonds around peptide followed by SPC/E. H3.00 forms the least number of H-bonds, as it is a poor tetrahedral liquid. H1.56 and mTIP3P form a similar number of H-bonds, as is also seen in the case of the tetrahedral order parameter, where the two models show very similar qtet values. The conceptual basis of force-field parametrization that forms the background of recent variable solvent studies57−60,92 including this one may be summarized as follows. Simulations based on physically well-motivated parametric potentials are carried out in order to reproduce as much structural, dynamic, and thermodynamic information from the experiments.4 Original force-field parametrization was based on fitting or tuning the parameters to the best extent possible to the available experimental data on biomolecules, typically under ambient conditions.65 Different biomolecular force fields parametrized with different rigid-body water models, therefore, tend to agree on the folded structure of peptides which were used in the parametrization procedure. However, this provides only a limited view of the free energy landscape. The present study demonstrates that, by suitable choice of small alanine oligopeptides and solvent models, considerable insight can be obtained on the competition between intrapeptide, peptide− water, and water−water interactions in determining the conformational free energy landscape. The importance of such a competition between intrapeptide, peptide−water, and water−water hydrogen bonds in determining the α-helical structures was first pointed out by Scheraga.30 Our results show that, in H3.00, intrapeptide hydrogen bonds are disfavored, and helicity is low, because water molecules prefer to form peptide−water hydrogen bonds rather than water−water hydrogen bonds. As electrostatic effects increase in going from H3.00 to H1.56 to SPC/E, the water−water hydrogen bonds become stronger and the peptide−water hydrogen bonds become less favorable energetically, which favors intrapeptide hydrogen bonding and the associated helical structure. The number of intrapeptide hydrogen bonds in the mTIP3P and TIP4P water models is very similar and slightly lower than that of SPC/E at low temperatures, in keeping with trends in the fraction of helical configurations, despite the fact that TIP4P has the largest number of peptide−water and water−water hydrogen bonds. Thus, it is crucial to obtain physical insights and experimental data on hydrogen-bond energetics and dynamics in solutions.

computed for (a) peptide backbone, (b) peptide−water, and (c) water−water interactions for acA5nme in different solvent models. A hydrogen bond is said to be formed if the donor (D)−acceptor (A) distance is less than 3.5 Å and the A−D−H angle is less than 30°.89−91 The peptide backbone hydrogen bonds form between the amide NH group and CO groups of the peptide backbone. The peptide−water hydrogen bonds are formed when water interacts with the NH group of amide or oxygen of the CO group of the peptide backbone. The water−water hydrogen bonds are computed for water molecules lying within the first hydration shell of the peptide, i.e., water molecules lying within 5 Å from any heavy atom of the peptide. Figure 8a shows the average number of peptide backbone Pep hydrogen bonds (⟨NHB ⟩) for acA5nme as a function of temperature. The pentapeptide forms three backbone Hbonds in the native α-helical state. It is found that, on average, the maximum number of H-bonds ranges up to 1.4 as in the case of SPC/E at 250 K. This is because, out of the three Hbonds, two H-bonds require participation of CO and NH groups of the capping residues (acetyl and N-methylamide) present at the termini of the peptide, making the H-bonds more vulnerable to breakage, thereby leading to a low number of total backbone H-bonds. The trend in (⟨NPep HB ⟩) in different solvent models follows the same order as seen for Fhel, which is H3.00 < H1.56 < mTIP3P ≈ TIP4P < SPC/E. This similarity in trends of Fhel and ⟨NPep HB ⟩ is expected, as a higher number of peptide backbone H-bonds results in a higher helicity of the peptide. Figure 8b shows the average number of peptide−water hydrogen bonds for acA5nme as a function of temperature. The number of peptide−water H-bonds decreases with temperature. The trend in solvent models observed for the number of peptide−water hydrogen bonds is mTIP3P < H1.56 ≈ SPC/E < H3.00 < TIP4P. Figure 8c shows the average number of hydrogen bonds formed between water molecules present around the peptide in the first hydration shell in different solvent models, which decrease as a function of temperature. The number of water−water hydrogen bonds increases in the order H3.00 < H1.56 ≈ mTIP3P < SPC/E < TIP4P. This trend is similar to that seen for the qtet(r) profile in different solvent models, indicating that the differences in solvent models in terms of forming a tetrahedral network in a bulk system can be seen in the hydrogen bonding network of hydration shell waters also. TIP4P forms the highest number of 11115


Article



Figure 6. Average number of water molecules in the first hydration shell, ⟨nfirst wat ⟩, of (a) acAnme and (b) acA5nme, as a function of temperature in different solvent models.

in solvent as compared to the i−(i + 4) intrapeptide hydrogen bonds in the larger α-helical peptides. This is in keeping with the earlier metadynamics study of dialanine in various rigidbody water models.67 The probability of occupancy of the second quadrant (C5 + C7eq) and third quadrant (αR + β2) is greater than 44% and less than 50%, respectively, in all solvent models except for the H3.00 model where the occupancy of (αR + β2) is significantly preferred over (C5 + C7eq) population. In H3.00, the repulsion−dispersion interaction is more important than the electrostatic effects, as a result of which packing effects are much more favorable. The increased packing effects of the solvent together with the special nature of the intrapeptide hydrogen bond in the dipeptide give rise to the preference for the (αR + β2) conformations over (C5 + C7eq) conformations (that include the PPII conformation) in the H3.00 model. The pentapeptide, acA5nme, where there are three intrapeptide hydrogen bonds between the ith and the (i + 4)th backbone CO and NH groups, shows a conformational free energy landscape with a much greater degree of sensitivity to the choice of solvent model. The small size and simplicity of the peptide allow us to capture how the interplay between intrapeptide, peptide−water, and water−water hydrogen bonding determines the competition between αR-helical and polyproline (PPII) helical structures, which provides some interesting insights that should be applicable also to more complex peptides. Table 8 summarizes our results by showing the increasing order in various secondary structure and hydration shell metrics as a function of the solvent model below 300 K. In all solvents, with the CHARMM22 force field, the helical structure remains dominant. The helical fraction (Fhel) is maximum in SPC/E, very similar in mTIP3P and TIP4P, and minimum in H3.00; this trend is identical to that observed for the number of intrapeptide bonds (NPep HB ) which hold the helix together. The β-sheet and PPII fractions show trends that are anticorrelated with Fhel and are strongly correlated with the size and positional fluctuation-related

Figure 5. Number of water molecules, n(r), as a function of distance from (a) acAnme at 300 K and (b, c) acA5nme at 250 and 298.43 K, respectively, in different solvent models. The shaded regions demarcate different hydration shells of peptides.

4. CONCLUSIONS We study the effect of different rigid-body (mTIP3P, TIP4P, SPC/E) and hybrid (H1.56, H3.00) water models on the conformational free energy landscape of the alanine oligopeptides (acAnmeand acA5nme) using the CHARMM22 force field. The conformational free energy landscape of these peptides is mapped out as a function of the backbone dihedral angles defining the Ramachandran plot as a function of temperature. Structural metrics related to size and positional fluctuations of the two peptides are computed. The hydration shell (or solvation shell) structure of the peptide in different solvents is studied, focusing in particular on the occupancy, tetrahedral structure, and hydrogen-bonding characteristics within the first hydration shell. The intrapeptide hydrogen bond in the dipeptide, acAnme, forms part of relatively small 5- and 7-membered rings and is, therefore, relatively less exposed and less sensitive to the change 11116


Article



Figure 8. Average number of hydrogen bonds formed between (a) PW peptide backbone atoms (NPep HB ), (b) peptide−water (NHB ), and (c) ) in the first hydration shell of the peptide, as a water molecules (NWW HB function of temperature for acA5nme in different solvent models.

solvents varies significantly and seems to play no major role. Only the H3.00 model for which dispersion−repulsion interactions are important forms the least number of water− water hydrogen bonds, due to which the number of peptide− water hydrogen bonds increases and the helical structure become relatively unfavorable. The SPC/E, H1.56, and H3.00 liquids show a clear trend where better solvent hydrogenbonding propensities weaken peptide−water hydrogen bonds and stabilize the helical structure. In a completely nonpolar, non-hydrogen-bonding solvent (sometimes represented by the vacuum in simulations), only intrapeptide hydrogen bonds can be formed and the α-helical structure is well-known to be the minimum.30,31,69,93 Results presented in ref 58 suggest that this may also be true for longer peptides, such as dodecalanine, in solvents with polarity similar to H3.00, though their conclusions are based on time-dependent RMSDs, rather than free energy landscapes. It is also evident that an adequate representation of the peptide−solvent hydrogen bonds, together with solvent−solvent hydrogen bonds, is required to model protein secondary structures. The tendency of short polyalanine chains to adopt PPII structures in aqueous media may, in fact, reflect that, for these systems which have a large solvent accessible surface area, the competition between different types of hydrogen bonding is particularly pronounced. A general conclusion that emerges from our study is that more insights into the nature of solvent−peptide and solvent− solvent interactions must be built into the design of protein

Figure 7. Average tetrahedral order, qtet(r), of water molecules as a function of distance from (a) acAnme at 300 K and (b, c) acA5nme at 250 and 298.43 K, respectively, in different solvent models. The shaded regions demarcate different hydration shells of peptides.

metrics, reflecting the more extended nature of these secondary structures. The number of water−water hydrogen bonds (NWW HB ) and mean tetrahedral order within the first hydration shell show identical rankings and correlate also with trends in these WW properties for the bulk solvents. Interestingly, NPep HB and NHB PW show very similar trends, while NHB is completely different; neither of these three quantities correlate exactly with trends in hydrogen bond energies from dimer data shown in Table 8. The rigid-body water models all have very similar electrostatic parameters, and show a similar tendency to stabilize the helical structure. SPC/E has the strongest water−water hydrogen bonds, which must reduce the tendency to form peptide−water hydrogen bonds, resulting in helical stablization, especially at low temperatures when energetic effects dominate. The mTIP3P, TIP4P, and H1.56 models have comparable EHB values. TIP4P and mTIP3P show an almost identical temperature dependence of the helical fraction, while H1.56 is somewhat lower. The overall dielectric constant of the 11117


Article

The Journal of Physical Chemistry B force fields. Parameterization approaches that incorporate available experimental data on short peptides that can sample large regions of the Ramachandran plots would be particularly valuable.

■

at Atomic Resolution Using NMR Spectroscopy. Chem. Rev. 2014, 114, 6632−6660. (15) Fink, A. L. Natively Unfolded Proteins. Curr. Opin. Struct. Biol. 2005, 15, 35−41. (16) Tran, H. T.; Mao, A.; Pappu, R. V. Role of Backbone-Solvent Interactions in Determining Conformational Equilibria of Intrinsically Disordered Proteins. J. Am. Chem. Soc. 2008, 130, 7380−7392. (17) Schweitzer-Stenner, R. Conformational Propensities and Residual Structures in Unfolded Peptides and Proteins. Mol. BioSyst. 2012, 8, 122−133. (18) Best, R. B.; Hummer, G. Optimized Molecular Dynamics Force Fields Applied to the Helix-Coil Transition of Polypeptides. J. Phys. Chem. B 2009, 113, 9004−9015. (19) Dill, K. A.; Truskett, T. M.; Vlachy, V.; Hribar-Lee, B. Modelling Water, the Hydrophobic Effect and Ion Solvation. Annu. Rev. Biophys. Biomol. Struct. 2005, 34, 173−199. (20) Southall, N. T.; Dill, K. A.; Haymet, A. D. J. A View of the Hydrophobic Effect. J. Phys. Chem. B 2002, 106, 521−533. (21) Chandler, D. Interfaces and the driving force of hydrophobic assembly. Nature 2005, 437, 640−647. (22) Matysiak, S.; Debenedetti, P. G.; Rossky, P. J. Dissecting the Energetics of Hydrophobic Hydration of Polypeptides. J. Phys. Chem. B 2011, 115, 14859−14865. (23) Freed, K. F. Perturbative Many-Body Expansion for Electrostatic Energy and Field for System of Polarizable Charged Spherical Ions in a Dielectric Medium. J. Chem. Phys. 2014, 141, 034115. (24) Paschek, D. Temperature Dependence of the Hydrophobic Hydration and Interaction of Simple Solutes: An Examination of Five Popular Water Models. J. Chem. Phys. 2004, 120, 6674−6690. (25) Hummer, G.; Garde, S.; García, A. E.; Pratt, L. R. New Perspectives on Hydrophobic Effects. Chem. Phys. 2000, 258, 349− 370. (26) Ashbaugh, H. S.; Collett, N. J.; Hatch, H. W.; Staton, J. A. Assessing the Thermodynamic Signatures of Hydrophobic Hydration for Several Common Water Models. J. Chem. Phys. 2010, 132, 124504. (27) Berne, B. J.; Weeks, J. D.; Zhou, R. Dewetting and Hydrophobic Interaction in Physical and Biological Systems. Annu. Rev. Phys. Chem. 2009, 60, 85−103. (28) Acharya, H.; Vembanur, S.; Jamadagni, S.; Garde, S. Mapping Hydrophobicity at the Nanoscale: Applications to Heterogeneous Surfaces and Proteins. Faraday Discuss. 2010, 146, 353−365. (29) Prabhu, N.; Sharp, K. Protein-Solvent Interactions. Chem. Rev. 2006, 106, 1616−1623. (30) Vila, J. A.; Ripoll, D. R.; Scheraga, H. A. Physical Reasons for the Unusual α-Helix Stabilization Afforded by Charged or Neutral Polar Residues in Alanine-rich Peptides. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 13075−13079. (31) Levy, Y.; Jortner, J.; Becker, O. M. Solvent Effects on the Energy Landscapes and Folding Kinetics of Polyalanine. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 2188−2193. (32) Rhee, Y. M.; Sorin, E. J.; Jayachandran, G.; Lindahl, E.; Pande, V. S. Simulations of the Role of Water in the Protein Folding Mechanism. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 6456−6461. (33) Nayar, D.; Chakravarty, C. Sensitivity of Local Hydration Behaviour and Conformational Preferences of Peptides to Choice of Water Model. Phys. Chem. Chem. Phys. 2014, 16, 10199−10213. (34) Vlachy, V.; Hribar-Lee, B.; Kalyuzhnyi, Y. V.; Dill, K. A. ShortRange Interactions: From Simple Ions to Polyelectrolyte Solutions. Curr. Opin. Colloid Interface Sci. 2004, 9, 128−132. (35) Hribar, B.; Southall, N. T.; Vlachy, V.; Dill, K. A. How Ions Affect the Structure of Water. J. Am. Chem. Soc. 2002, 124, 12302− 12311. (36) Yu, H.; Hansson, T.; van Gunsteren, W. F. Development of a Simple, Self-Consistent Polarizable Model for Liquid Water. J. Chem. Phys. 2003, 118, 221−234. (37) Vega, C.; Abascal, J. L. F.; Conde, M. M.; Aragones, J. L. What Ice can Teach us about Water Interactions: A Critical Comparison of the Performance of Different Water Models. Faraday Discuss. 2009, 141, 251−276.

AUTHOR INFORMATION

Corresponding Author

*Phone: (+)91-11-2659-1510. Fax: (+)91-11-2686-2122. Email: [email protected]. Notes


The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

This work was financially supported by the Department of Science and Technology, New Delhi. Computational resources from the Computer Services Centre of I.I.T.-Delhi are acknowledged. D.N. thanks the Council of Scientific and Industrial Research, New Delhi, for support through a Senior Research Fellowship. It is a pleasure to dedicate this paper to Prof. Biman Bagchi in recognition of his many contributions as a physical chemist.

(1) Creighton, T. E. Proteins: Structures and Molecular Properties; W. H. Freeman and Company: New York, 1992. (2) Cox, M. M.; Nelson, D. L. Lehninger Principles of Biochemistry; W.H. Freeman and Company: New York, 2008. (3) Dill, K. A. Dominant Forces in Protein Folding. Biochemistry 1990, 29, 7133−7155. (4) Leach, A. Molecular Modelling: Principles and Applications; Addison Wesley Longman Limited: Essex, England, 1998. (5) Schlick, T. Molecular Modeling and Simulation: An Interdisciplinary Guide; Springer-Verlag: New York, 2002. (6) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (7) Scott, W. R. P.; Hünenberger, P. H.; Tironi, I. G.; Mark, A. E.; Billeter, S. R.; Fennen, J.; Torda, A. E.; Huber, T.; Kruger, P.; van Gunsteren, W. F. The GROMOS Biomolecular Simulation Program Package. J. Phys. Chem. A 1999, 103, 3596−3607. (8) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (9) Best, R. B.; Miller, C.; Mittal, J. Role of Solvation in PressureInduced Helix Stabilization. J. Chem. Phys. 2014, 141, 22D522. (10) Chatterjee, P.; Bagchi, S.; Sengupta, N. The Non-Uniform Early Structural Response of Globular Proteins to Cold Denaturing Conditions: A Case Study with Yfh1. J. Chem. Phys. 2014, 141, 205103. (11) Piana, S.; Klepeis, J. L.; Shaw, D. E. Assessing the Accuracy of Physical Models used in Protein-Folding Simulations: Quantitative Evidence from Long Molecular Dynamics Simulation. Curr. Opin. Struct. Biol. 2014, 24, 98−105. (12) Lindorff-Larsen, K.; Trbovic, N.; Maragakis, P.; Piana, S.; Shaw, D. E. Structure and Dynamics of an Unfolded Protein Examined by Molecular Dynamics Simulation. J. Am. Chem. Soc. 2012, 134, 3787− 3791. (13) Nerenberg, P. S.; Head-Gordon, T. Optimizing Protein-Solvent Force Fields to Reproduce Intrinsic Conformational Preferences of Model Peptides. J. Chem. Theory Comput. 2011, 7, 1220−1230. (14) Jensen, M. R.; Zweckstetter, M.; Huang, J.-R.; Blackledge, M. Exploring Free-Energy Landscapes of Intrinsically Disordered Proteins 11118


Article



(59) Sorin, E. J.; Rhee, Y. M.; Shirts, M. R.; Pande, V. S. The Solvation Interface is a Determining Factor in Peptide Conformational Preferences. J. Mol. Biol. 2006, 356, 248−256. (60) Piana, S.; Donchev, A. G.; Robustelli, P.; Shaw, D. E. Water Dispersion Interactions Strongly Influence Simulated Structural Properties of Disordered Protein States. J. Phys. Chem. B 2015, 119, 5113−5123. (61) Lynden-Bell, R. M. Towards Understanding Water: Simulation of Modified Water Models. J. Phys.: Condens. Matter 2010, 22, 284107. (62) Prasad, S.; Chakravarty, C. Onset of Simple Liquid Behaviour in Modified Water Models. J. Chem. Phys. 2014, 140, 164501. (63) Drozdov, A. N.; Grossfield, A.; Pappu, R. V. Role of Solvent in Determining Conformational Preferences of Alanine Dipeptide in Water. J. Am. Chem. Soc. 2004, 126, 2574−2581. (64) García, A. E.; Sanbonmatsu, K. Y. α-Helical Stabilization by Side Chain Shielding of Backbone Hydrogen Bonds. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 2782−2787. (65) Gnanakaran, S.; Garcıa, A. E. Validation of an All-Atom Protein Force Field: From Dipeptides to Larger Peptides. J. Phys. Chem. B 2003, 107, 12555−12557. (66) Samanta, A.; Chen, M.; Yu, T.-Q.; Tuckerman, M.; Weinan, E. Sampling Saddle Points on a Free Energy Surface. J. Chem. Phys. 2014, 140, 164109. (67) Vymetal, J.; Vondrasek, J. Metadynamics As a Tool for Mapping the Conformational and Free-Energy Space of Peptides- The Alanine Dipeptide Case Study. J. Phys. Chem. B 2010, 114, 5632−5642. (68) Ishizuka, R.; Huber, G. A.; McCammon, J. A. Solvation Effect on the Conformations of Alanine Dipeptide: Integral Equation Approach. J. Phys. Chem. Lett. 2010, 1, 2279−2283. (69) Hazel, A.; Chipot, C.; Gumbart, J. C. Thermodynamics of Decaalanine Folding in Water. J. Chem. Theory Comput. 2014, 10, 2836− 2844. (70) Ramakrishnan, V.; Ranbhor, R.; Durani, S. Existence of Specific Folds in Polyproline II Ensembles of an Unfolded Alanine Peptide Detected by Molecular Dynamics. J. Am. Chem. Soc. 2004, 126, 16332−16333. (71) Mezei, M.; Fleming, P. J.; Srinivasan, R.; Rose, G. D. Polyproline II Helix Is the Preferred Conformation for Unfolded Polyalanine in Water. Proteins: Struct., Funct., Genet. 2004, 55, 502−507. (72) Creamer, T. P.; Campbell, M. N. Determinants of the Polyproline II Helix from Modeling Studies. Adv. Protein Chem. 2002, 62, 263−282. (73) Cheam, T. C.; Krimm, S. Ab Initio Force Fields of Alanine Dipeptide in C5 and C7 Conformations. J. Mol. Struct.: THEOCHEM 1989, 188, 15−43. (74) Mackerell, A. D., Jr.; Feig, M.; Brooks, C. L., III Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25, 1400−1415. (75) Feig, M. Is Alanine Dipeptide a Good Model for Representing the Torsional Preferences of Protein Backbones? J. Chem. Theory Comput. 2008, 4, 1555−1564. (76) Tobias, D. J.; Brooks, C. L., III Conformational Equilibrium in the Alanine Dipeptide in the Gas Phase and Aqueous Solution: A Comparison of Theoretical Results. J. Phys. Chem. 1992, 96, 3864− 3870. (77) Best, R. B.; Buchete, N. V.; Hummer, G. Are Current Molecular Dynamics Force Fields too Helical? Biophys. J. 2008, 95, L07−L09. (78) Huang, J.; MacKerell, A. D., Jr. CHARMM36 All-Atom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 2013, 34, 2135−2145. (79) Buchete, N.-V.; Hummer, G. Coarse Master Equations for Peptide Folding Dynamics. J. Phys. Chem. B 2008, 112, 6057−6059. (80) Buchete, N.-V.; Hummer, G. Peptide Folding Kinetics from Replica Exchange Molecular Dynamics. Phys. Rev. E 2008, 77, 030902(R). (81) Li, X.; O’Brien, C. P.; Collier, G.; Vellore, N. A.; Wang, F.; Latour, R. A.; Bruce, D. A.; Stuart, S. J. An Improved Replica-Exchange

(38) Vega, C.; Abascal, J. L. F. Relation Between the Melting Temperature and the Temperature of Maximum Density for the Most Common Models of Water. J. Chem. Phys. 2005, 123, 144504. (39) Abascal, J. L. F.; Vega, C. Dipole-Quadrupole Force ratios Determine the Ability of Potential Models to Describe the Phase Diagram of Water. Phys. Rev. Lett. 2007, 98, 237801−237804. (40) Agarwal, M.; Alam, M. P.; Chakravarty, C. Thermodynamic, Diffusional, and Structural Anomalies in Rigid-Body Water Models. J. Phys. Chem. B 2011, 115, 6935−6945. (41) Nayar, D.; Agarwal, M.; Chakravarty, C. Comparison of Tetrahedral Order, Liquid State Anomalies, and Hydration Behavior of mTIP3P and TIP4P Water Models. J. Chem. Theory Comput. 2011, 7, 3354−3367. (42) Nayar, D.; Chakravarty, C. Water and Water-like Liquids: Relationships Between Structure, Entropy and Mobility. Phys. Chem. Chem. Phys. 2013, 15, 14162−14177. (43) Agarwal, M.; Singh, M.; Sharma, R.; Parvez Alam, M.; Chakravarty, C. Relationship between Structure, Entropy, and Diffusivity in Water and Water-Like Liquids. J. Phys. Chem. B 2010, 114, 6995−7001. (44) Shadrack Jabes, B.; Nayar, D.; Dhabal, D.; Molinero, V.; Chakravarty, C. Water and Other Tetrahedral Liquids: Order, Anomalies and Solvation. J. Phys.: Condens. Matter 2012, 24, 284116. (45) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (46) Neria, E.; Fischer, S.; Karplus, M. Simulation of Activation Free Energies in Molecular Systems. J. Chem. Phys. 1996, 105, 1902−1921. (47) Lynden-Bell, R. M.; Debenedetti, P. G. Computational Investigation of Order, Structure, and Dynamics in Modified Water Models. J. Phys. Chem. B 2005, 109, 6527−6534. (48) Lynden-Bell, R. M.; Head-Gordon, T. Solvation in Modified Water Models: Towards Understanding Hydrophobic Effects. Mol. Phys. 2006, 104, 3593−3605. (49) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (50) Nutt, D. R.; Smith, J. C. Molecular Dynamics Simulations of Proteins: Can the Explicit Water Model Be Varied? J. Chem. Theory Comput. 2007, 3, 1550−1560. (51) Glass, D. C.; Krishnan, M.; Nutt, D. R.; Smith, J. C. Temperature Dependence of Protein Dynamics Simulated with Three Different Water Models. J. Chem. Theory Comput. 2010, 6, 1390−1400. (52) Agarwal, M.; Kushwaha, H. R.; Chakravarty, C. Local Order, Energy and Mobility of Water Molecules in the Hydration Shell of Small Peptides. J. Phys. Chem. B 2010, 114, 651−659. (53) Florová, P.; Sklenovský, P.; Banás,̆ P.; Otyepka, M. Explicit Water Models Affect the Specific Solvation and Dynamics of Unfolded Peptides While the Conformational Behavior and Flexibility of Folded Peptides Remain Intact. J. Chem. Theory Comput. 2010, 6, 3569−3579. (54) Best, R. B.; Mittal, J. Microscopic Events in β-hairpin Folding from Alternative Unfolded Ensembles. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 11087−11092. (55) Hatch, H. W.; Stillinger, F. H.; Debenedetti, P. G. Computational Study of the Stability of the Miniprotein Trp-Cage, the GB1 βHairpin, and the AK16 Peptide, under Negative Pressure. J. Phys. Chem. B 2014, 118, 7761−7769. (56) Paschek, D.; Day, R.; García, A. E. Influence of Water-Protein Hydrogen Bonding on the Stability of Trp-cage miniprotein. A Comparison between the TIP3P and TIP4P-Ew Water Models. Phys. Chem. Chem. Phys. 2011, 13, 19840−19847. (57) Gu, B.; Zhang, F. S.; Wang, Z. P.; Zhou, H. Y. Solvent-Induced DNA Conformational Transition. Phys. Rev. Lett. 2008, 100, 088104. (58) Shen, H.; Cheng, W.; Zhang, F.-S. Structural Conservation of the Short α-Helix in Modified Higher and Lower Polarity Water Solutions. RSC Adv. 2015, 5, 9627−9634. 11119


Article


The Journal of Physical Chemistry B Sampling Method: Temperature Intervals with Global Energy Reassignment. J. Chem. Phys. 2007, 127, 164116. (82) Lockhart, C.; Klimov, D. K. Molecular Interactions of Alzheimer’s Biomarker FDDNP with Aβ Peptide. Biophys. J. 2012, 103, 2341−2351. (83) García, A. E.; Herce, H.; Paschek, D. Simulations of Temperature and Pressure Unfolding of Peptides and Proteins with Replica Exchange Molecular Dynamics. Annu. Rep. Comput. Chem. 2006, 2, 83−93. (84) Sindhikara, D.; Meng, Y.; Roitberg, A. E. Exchange frequency in Replica Exchange Molecular Dynamics. J. Chem. Phys. 2008, 128, 024103. (85) Ramachandran, G.; Sasisekharan, V. Conformation of Polypeptides and Proteins. Adv. Protein Chem. 1968, 23, 283−437. (86) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (87) Maiorov, V. N.; Crippen, G. M. Significance of Root-MeanSquare Deviation in Comparing Three-dimensional Structures of Globular Proteins. J. Mol. Biol. 1994, 235, 625−634. (88) Kabsch, W. A Discussion of the Solution for the Best Rotation to Relate Two Sets of Vectors. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1978, 34, 827−828. (89) Bolhuis, P. G. Transition-Path Sampling of β-hairpin folding. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 12129−12134. (90) Tarek, M.; Tobias, D. J. Role of Protein-Water Hydrogen Bond Dynamics in the Protein Dynamical Transition. Phys. Rev. Lett. 2002, 88, 138101. (91) Xu, H.; Berne, B. J. Hydrogen-Bond Kinetics in the Solvation Shell of a Polypeptide. J. Phys. Chem. B 2001, 105, 11929−11932. (92) Best, R. B.; Zheng, W.; Mittal, J. Balanced Protein-Water Interactions Improve Properties of Disordered Proteins and NonSpecific Protein Association. J. Chem. Theory Comput. 2014, 10, 5113− 5124. (93) Vila, J.; Williams, R. L.; Grant, J. A.; Wǿjcik, J.; Scheraga, H. A. The Intrinsic Helix-Forming Tendency of L-Alanine. Proc. Natl. Acad. Sci. U. S. A. 1992, 89, 7821−7825.

11120


Free Energy Landscapes of Alanine Oligopeptides in Rigid-Body and

Recommend Documents