Reordering Hydrogen Bonds Using Hamiltonian Replica Exchange

Apr 9, 2009 - significantly accelerates the folding process (Wolf, M. G.; de Leeuw, S. W. Biophys. J. 2008, 94, 3742). As the biasing potentials are o...
0 downloads 0 Views 4MB Size
Download PDF
6484

J. Phys. Chem. B 2009, 113, 6484–6494

Reordering Hydrogen Bonds Using Hamiltonian Replica Exchange Enhances Sampling of Conformational Changes in Biomolecular Systems Jocelyne Vreede,*,† Maarten G. Wolf,‡ Simon W. de Leeuw,§ and Peter G. Bolhuis† Van’t Hoff Institute for Molecular Sciences, UniVersity of Amsterdam, Amsterdam, The Netherlands, DelftChemTech, UniVersity of Technology Delft, Delft, The Netherlands, and Theoretical Chemistry, UniVersity Leiden, Leiden, The Netherlands ReceiVed: October 31, 2008; ReVised Manuscript ReceiVed: February 16, 2009

Hydrogen bonds play an important role in stabilizing (meta-)stable states in protein folding. Hence, they can potentially be used as a way to bias these states in molecular simulation methods. Previously, Wolf et al. showed that applying repulsive and attractive hydrogen bond biasing potentials in an alternating way significantly accelerates the folding process (Wolf, M. G.; de Leeuw, S. W. Biophys. J. 2008, 94, 3742). As the biasing potentials are only active during a fixed time interval, this alternating scheme does not represent a thermodynamic equilibrium. In this work, we present a Hamiltonian replica exchange molecular dynamics (REMD) scheme that aims to shuffle and reorder hydrogen bonds in the protein backbone. We therefore apply adapted hydrogen bond potentials in a Hamiltonian REMD scheme, which we call hydrogen bond switching (HS). To compare the performance of the HS to a standard REMD method, we performed HS and temperature REMD simulations of a β-heptapeptide in methanol. Both methods sample the conformational space to a similar extent. As the HS simulation required only five replicas, while the REMD simulation required 20 replicas, the HS method is significantly more efficient. We tested the HS method also on a larger system, 16-residue polyalanine in water. Both of the simulations starting from a completely unfolded and a folded conformation resulted in an ensemble with, apart from the starting structure, similar conformational minima. We can conclude that the HS method provides an efficient way to sample the conformational space of a protein, without requiring knowledge of the folded states beforehand. In addition, these simulations revealed that convergence was hampered by replicas having a preference for specific biasing potentials. As this sorting effect is inherent to any Hamiltonian REMD method, finding a solution will result in an additional increase in the efficiency of Hamiltonian REMD methods in general. Introduction Molecular simulation is a powerful tool to perform detailed investigations of molecular processes in biological systems, such as protein conformational changes. The study of a protein folding reaction at atomic resolution poses a problem, as it involves time scales that exceed the range of conventional atomistic methods, such as molecular dynamics (MD). The folding of a protein requires traversing a rugged free energy landscape with many high barriers. Consequently, many (local) free energy minima exist including unfolded states, partially folded states, and the native state, which are not trivial to escape from. Nevertheless, the folding of a protein is overall a guided process,1 in which several types of interactions play an important role. Hydrophobic interactions are the driving force behind the collapse of a polypeptide chain,2 partially caused by the entropy gain of the solvent molecules. During this collapse, the formation of intramolecular hydrogen bonds can compensate the high free energy cost associated with burying unsatisfied (partially) charged groups and loss of solvent interactions.3 The free energy barrier associated with breaking a single intramolecular hydrogen bond is approximately 3 kBT in a protein in water, which is not a high barrier in itself.4-6 The formation of one non-native hydrogen bond, however, may lead to additional non-native * To whom correspondence should be addressed. E-mail: [email protected]. † University of Amsterdam. ‡ University of Technology Delft. § University Leiden.

contacts, eventually resulting in a misfolded state. Such a misfolded structure may have to unfold again for a new attempt to find the native state. As such, even one incorrectly formed hydrogen bond may lead to kinetic traps that slow down the folding process enormously. Introducing potentials that target the breaking and formation of backbone hydrogen bonds can therefore achieve a significant speed-up when simulating protein folding in explicit solvent, as was recently shown.7 When applying hydrogen bond forming and hydrogen bond breaking potentials in an alternate way, fast reordering of the backbone hydrogen bonds occurs, leading to accelerated folding of the protein.7 In this alternating hydrogen bond potential scheme, the time interval during which the hydrogen bond breaking/forming potential acts on the system is tuned to optimize folding times. Also, alternating hydrogen bond potentials do not conserve the Hamiltonian. These two aspects prevent obtaining thermodynamic properties directly from such simulations. Alternatively, replica exchange MD (REMD), also known as parallel tempering,8 is a method often used to escape local minima in biomolecular systems.9 In REMD, a number of MD simulations run simultaneously at different temperatures. By letting the replicas exchange temperatures, according to a Metropolis criterion, a system can diffuse through temperature space and overcome barriers at high temperature, while sampling the stable states at the temperature of interest. Afterward, the free energy can be obtained by constructing a histogram along

10.1021/jp809641j CCC: $40.75  2009 American Chemical Society Published on Web 04/09/2009

Reordering Hydrogen Bonds Using Hamiltonian REMD any order parameter. Several new insights resulted from the application of this method to protein folding processes.10-12 When using an implicit solvent model, the gain in efficiency in comparison to conventional MD is high.12 In order to investigate the role of solvent molecules in protein conformational changes, the use of explicit solvent models is required, resulting in systems containing many particles. As neighboring replicas must have sufficient overlap in potential energy for an efficient exchange, many replicas are needed between the temperature of interest and high temperature that facilitate fast conformational changes.13,14 Improving the efficiency of parallel tempering is an active area of research and focuses on adjusting the temperature distribution dynamically15 or identifying entropic bottlenecks.16-19 Another possible solution to this problem is letting the replicas switch (part of) their Hamiltonians instead of their temperatures.20-22 A smooth interpolation between the regular Hamiltonian and one that enables swift conformational changes of the protein allows for efficient sampling of the conformational space. Such Hamiltonians target, for example, hydrophobic interactions,23 Lennard-Jones interactions in general,24 or torsional interactions.25 In this work, we present continuous hydrogen bond potentials that we use in a Hamiltonian replica exchange scheme: hydrogen bond switching (HS). Hydrogen bond switching is in principle applicable to any system containing hydrogen bonds. The potentials are derived in such a way that only the atoms in the weakest hydrogen bond will feel a repulsive interaction, creating a less deep free energy minimum, leading to the breaking of that hydrogen bond. Similarly, the atoms closest to forming a hydrogen bond will undergo an attractive interaction, creating a wider free energy minimum, leading to the formation of that hydrogen bond. In this work, we only targeted the hydrogen bonds in the protein backbone. Before setting up an HS simulation, no knowledge of hydrogen bonds in stable states or the order in which hydrogen bonds will break is required. As a test system, we use a β-heptapeptide solvated in methanol. This system is small enough to enable a direct comparison between conventional MD and temperature REMD, in favor of REMD, as Periole et al. have shown in great detail.26 We show that there is a significant gain in efficiency, compared to REMD, when using HS simulations to extensively sample the conformational space of this system. We also study 16-residue polyalanine, A16, a system with more hydrogen bonds. As shown by solid state NMR, polyalanine tends to adopt an R-helical conformation.27 In solution, polyalanine forms fibrillar aggregates when containing more than 15 alanine residues,28 indicating that the R-helical conformation is not stable. Using HS simulations, we have investigated the conformational space of the (un)folding of polyalanine, thus elucidating that A16 favors conformations with R-helical content, but also collapsed coil conformations occur. The HS simulations of polyalanine also uncovered a problem in Hamiltonian REMD methods that has not been explicitly described up to now. As replicas have a preference for a specific biasing potential, diffusion through λ-space is hampered. This problem is inherent to all Hamiltonian REMD methods that involve entropic barriers, presenting a challenge for future research. Methods Continuous Hydrogen Bond Biasing Potentials. Molecular simulations can provide thermodynamic properties, if the underlying methods, i.e., molecular dynamics (MD), sample the canonical ensemble correctly. As shown in previous work, applying attractive and repulsive hydrogen bond biasing po-

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6485 tentials in an alternating way induces fast conformational changes in protein systems.7 As these simulations are out of equilibrium, they cannot provide correct thermodynamic properties. In this work, we aim to apply attractive and repulsive hydrogen bond biasing potentials in a replica exchange scheme. To integrate correctly the equations of motion, using the standard Leap-frog algorithm, all potentials including the hydrogen bond biasing potentials have to be continuous. The biasing potentials as described in Wolf et al. are discontinuous7 and require modification. In addition, to reorder hydrogen bonds successfully, the biasing potentials must fulfill two other criteria. The first criterion is that the energy and the forces of the hydrogen bond biasing potentials must result in efficient hydrogen bond reordering. Reordering requires a large change in the hydrogen bond distance rather than in the hydrogen bond angle. Large changes in the hydrogen bond angle would result in the occurrence of unlikely structures. Second, as a hydrogen bond is a very specific interaction between a donor and an acceptor, excluding other interaction partners, the number of additional biasing potentials acting on the donor/acceptor must be restricted to one. For general applicability, the selection of these hydrogen bond interaction pairs must be based on the conformation in the simulation. The continuous hydrogen bond biasing potential Vhb(r) (eq 1) matches these criteria: D,A

Vhb(r) ) fc ·

∑ max[uhb,X(rD, rH , rA)] D

k

(1)

Here fc is the force constant, r represents the coordinates of all particles in the system, and uhb,X (rD,rHD,rA) is a positive definite hydrogen bond interaction function as a function of the positions of the hydrogen bond donor (rD), the hydrogen (rHD), and the hydrogen bond acceptor (rA). X represents either the donor D or the hydrogen bonded to the donor HD. This potential can be either attracive, aimed at forming hydrogen bonds, or repulsive, aimed at breaking hydrogen bonds. The repulsive (rep) and attractive (att) potentials are D,A

Vrep(r) ) fc ·

∑ max[uhb,D(rD, rH , rA)] D

k

(2)

D,A

Vatt(r) ) fc ·

∑ max[uhb,H (rD, rH , rA)] k

D

D

(3)

The hydrogen bond biasing potential Vhb(r) is a summation over all donors (D) and acceptors (A) in the protein backbone. The max function indicates that for each donor (acceptor) only the hydrogen bond pair with the largest value for the interaction function is selected to contribute to the potential. Apart from providing specificity, selecting only one interaction per donor (acceptor) minimizes the bias force introduced in the system when compared to including all possible donor-acceptor interactions. As the system evolves, these interacting pairs vary; i.e., one particular hydrogen bond interaction is lost while another is activated as the structure of the protein evolves during the simulation. At such a crossing point, the energy of both interactions is equal and therefore swapping the interaction potential is continuous. Also, the total momentum is conserved due to the use of pair interactions. The hydrogen bond interaction function uhb,X(rD,rHD,rA) comprises a distance function ud(dXA) and an angle function uθ(θDHDA):

uhb,X(rD, rHD, rA) ) ud(dXA) · uθ(θDHDA)

(4)

The distance function ud(dXA) depends on the distance dXA ) |rX - rA|.

6486

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

Vreede et al.

Figure 1. (a) Schematic representation of a donor-acceptor (D-A) pair including the hydrogen (HD). The labels indicate parameters used in the hydrogen bond potentials. In the repulsive potential, the distance between the donor D and the acceptor A (dDA) is used, and in the attractive potential, the distance between the hydrogen HD and A (dHDA) is used. θ is the angle between D, HD, and A. (b) Hydrogen bond interaction function corresponding to the repulsive hydrogen bond interaction (uhb,D). The interaction function uhb,D is a function of the acceptor (A) position in the plane around a fixed donor (D) and hydrogen (H). The color indicates the value of uhb, and the contour lines connect similar values of uhb,D at intervals of 0.2.

X represents either the donor D or the hydrogen bonded to the donor HD and A represents the acceptor. For the repulsive interaction, dXA is the distance between D and A (dDA). To ensure that hydrogen bonds are formed correctly for the attractive biasing potential, dXA is the distance between HD and A (dHDA). This way the hydrogen atom is correctly oriented during formation of a hydrogen bond. The explicit expression for the distance function ud(dXA):

ud(dXA) )

{

1 2(dXA - dmin)3 (dmax - dmin) 0

3

dXA < dmin 3(dXA - dmin)

2

-

(dmax - dmin)2

+ 1 dmin e dXA < dmax (5) dmax e dXA

This cubic spline function smoothly varies from the maximum value 1 for distance dXA smaller than a minimum cutoff distance dmin to the minimum value 0 for distance dXA larger than a maximum cutoff distance dmax. Using this particular spline function ensures that the derivative of this function, the bias force, is continuous, preventing the occurrence of spurious impulses. The angle function uθ(θDHDA) in eq 4 depends on the angle θDHDA between D, HD, and A. This is schematically shown in Figure 1a. The angle function uθ(θDHDA) is given by

uθ(θDHDA) )

{

1 (1 - 0.6) · 0.6

[

2(θDHDA - θmin)3 (θmax - θmin)3

-

3(θDHDA - θmin)2 (θmax - θmin)2

]

+1

θDHDA < θmin θmin e θDHDA < θmax (6) θmax e θDHDA

This spline function affects the D-HD-A angle and varies smoothly from its maximum value of 1 for angle θDHDA smaller than a minimum cutoff angle θmin to its minimum value of 0.6 for angle θDHDA larger than a maximum cutoff angle θmax. When using a minimum value of zero, rather than 0.6, the angle potential would generate very high forces. This could lead to

the combination of a very small dDA with a large θDA, resulting in the occurrence of unlikely or even nonphysical conformations. Again, with this function, the force is continuous. The values for the cutoff distances dmin and dmax should be chosen such that the repulsive biasing potentials will act on weak hydrogen bonds, flattening the free energy surface, and that the attractive biasing potentials act on donor-acceptor pairs that are close to forming a hydrogen bond, creating deeper and, more importantly, wider free energy minima. The cutoff angles θmin and θmax should ensure the direction in which the donor-acceptor pairs are pulled apart or drawn together. Figure 1b displays a graphical representation of the repulsive hydrogen bond interaction function uhb,D as a function of the position of the acceptor A for fixed donor D and hydrogen H positions. For a proper hydrogen bond, the acceptor will be at the red plateau (uhb,D ) 1), and consequently, no force resulting from the hydrogen bond potential acts on A and D. When the hydrogen bond is not in an optimal configuration, uhb,D will decrease. Then, the repulsive hydrogen bond potential results in forces exerted on both A and D, causing further disruption of the hydrogen bond. As the angle potential has a lower bound of 0.6, there is an intermediate configuration in which no forces act on the donor and acceptor. This plateau has no physical meaning, and the probability of finding a hydrogen bond in this area is very low, due to steric hindrance in the backbone. The purple plateau indicates the conformation with the lowest value for the biasing potential, in which the hydrogen bond is completely broken. Hydrogen Bond Switching. The continuous hydrogen bond biasing potentials we have developed can be used in an MD simulation of any (bio) molecular system to enhance hydrogen bond reordering. We can also use these potentials in a Hamiltonian REMD scheme in which the replicas exchange their Hamiltonians. As the Hamiltonian of a replica also includes the hydrogen bond biasing potential, exchanging two replicas i and j has the following acceptance probability pacc:20

pacc(i, j) )

{

1 ∆(i, j) e 0 -∆(i,j) ∆(i, j) > 0 e

∆(i, j) ) β[(Hj(ri) - Hj(rj)) - (Hi(rj) - Hi(ri))]

(7) (8)

with H the Hamiltonian, dependent on the coordinates r, and β equal to 1/kBT. The Hamiltonians in replicas i and j differ only

Reordering Hydrogen Bonds Using Hamiltonian REMD

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6487

in the contribution from the hydrogen bond biasing potential energy Ehb, reducing eq 8 to

∆(i, j) ) β[(Ehb,j(ri) - Ehb,j(rj)) - (Ehb,i(rj) - Ehb,i(ri))] (9) The hydrogen bond biasing potential energy Ehb can have two forms: repulsive, aimed at breaking hydrogen bonds (rep), and attractive, aimed at forming hydrogen bonds (att). In our Hamiltonian REMD scheme, a replica can have either a repulsive potential, an attractive potential, or no hydrogen bond biasing potential at all. The hydrogen bond biasing potential energy for replica i, Ehb, is given by

Ehb,i(ri) ) λrep,iVhb,rep(ri) + λatt,iVhb,att(ri)

(10)

with Vhb,rep and Vhb,att the repulsive and attractive hydrogen bond biasing potentials, respectively. The weighting constants λrep,i and λatt,i tune the strength of the repulsive and attractive biasing potentials in the replica. Equation 9 then becomes

∆(i, j) ) β[(λrep,i - λrep,j)(Vhb,rep(rj) - Vhb,rep(ri)) + (λatt,i - λatt,j)(Vhb,att(rj) - Vhb,att(ri))] (11) In our Hamiltonian REMD scheme, a replica can have an additional repulsive hydrogen bond interaction (0 < λrep,i e 1, λatt,i ) 0) or an additional attractive hydrogen bond interaction (λrep,i ) 0, 0 < λatt,i e 1). Also, there is a neutral replica in which both biasing potentials are switched off (λrep,i ) 0, λatt,i ) 0). As no biasing potentials are active, all configurations that visit this neutral replica sample the correct canonical distribution. Using the configurations from the neutral replica, we can calculate order parameters Q for which it is possible to construct a probability histogram P(Q) and compute the corresponding (relative) Landau free energy:

F(Q) ) -kBT ln P(Q)

(12)

Including properties of the rejected Monte Carlo moves in a replica exchange scheme significantly improves the accuracy of the histograms.29,30 This means that the contributions of all other replicas to the histogram of order parameters Q in the neutral replica are scaled with a Metropolis factor. The probability distribution P(Q) can then be estimated as N

P(Q) )

∑

[1 - min(1, exp(-∆))]δ(Qi - Q) +

j)1,j*i N

∑

[min(1, exp(-∆))]δ(Qj - Q) (13)

j)1,j*i

with ∆ from eq 11, Qi the value of the order parameter of replica i, δ(x) the Dirac delta function, and N the total number of replicas. The min function returns the smaller of its arguments. Note that this approach only improves the quality of the histograms by including all replicas, but it does not improve the sampling itself.29,30 Computational Details Systems. In all simulations, we used the GROMACS software package,31 in combination with the GROMOS96 43a1 force field.32,33 We investigated two systems: (1) a β-heptapeptide in βVAL-βALA-βLEU-βALA(Rmethyl)-βVAL-βALA-βLEU methanol and (2) a 16-residue polyalanine, A16, in water. The β-heptapeptide system was provided by Periole and required no further preparation.26 The system was 3028 atoms in size in a cubic periodic box with an edge of 3.96 nm. Nonbonded

TABLE 1: Parameters for the Hydrogen Bond Biasing Potentials parameter

attractive

dmina (nm) dmaxa (nm) θmina (rad) θmaxa (rad)

0.32 0.40 0.85 (145°) 0.50 (60°)

repulsive

fcrep fcatt simulation (kJ/mol) (kJ/mol)

0.23 NVE 0.40 alternating 0.85 (145°) HS 0.50 (60°)

10 50 1

10 50 1

a

The variables dmin, dmax, θmin, and θmax indicate the range in which the potential introduces a force to the system.

interactions were treated with a twin-range cutoff (1.0-1.4 nm). Interactions within the short-range cutoff were considered every time step, whereas interactions within the longer-range cutoff were evaluated every 10 steps, together with the updating of the pair list. Long-range electrostatic interactions beyond the cutoff were treated with a reaction-field correction.26 The polyalanine in SPC water was initialized in an R-helical f u , and a fully extended conformation, A16 . The conformation, A16 R-helical conformation was generated manually by taking residues 129-145 from PDB entry 101M and removing all side chains up to Cβ. This conformation was solvated in a periodic cubic box of 2950 SPC water molecules,34 resulting in a system of 8952 atoms. This system was equilibrated in the NpT ensemble for 100 ps at 300 K and 1 bar, resulting in a box with an average edge length of 4.53 nm. To obtain an unfolded conformation for A16, we performed an MD simulation of the f system at a temperature of 600 K and constant volume. A16 Polyalanine unfolded into a fully extended chain after 100 ps, after which we equilibrated the system again at 300 K and 1 bar for 100 ps. We used a time step of 2 fs throughout the preparation procedure. The van der Waals interaction cutoff was 1.2 nm. The electrostatic interactions were treated with the PME method with a cutoff of 0.9 nm.35 For both the β-heptapeptide and polyalanine, all bonds were constrained using LINCS.36 Temperature was kept constant at 300 K using the Nose´-Hoover thermostat.37,38 Pressure was kept constant at 1 bar using the Parrinello-Rahman barostat.39,40 NVE Simulations. To test if swapping the hydrogen bond interactions does not affect the continuity of the energy functions, we performed three simulations in the NVE ensemble, f system, of 50 ps, with a time step of 0.1 fs. The of the A16 three simulations had, respectively, no biasing potentials, repulsive, and attractive hydrogen bond biasing potentials. Both the repulsive and attractive potentials had a force constant of 10 kJ/mol. All parameters used for the hydrogen bond potentials are listed in Table 1. Alternating Hydrogen Bond Simulations. We performed a simulation of Au16 using the alternating hydrogen bond scheme.7 For 2.5 ns, the system underwent alternating cycles of 2 ps. During such a cycle, repulsive hydrogen bond biasing potentials were active for 0.5 ps and were switched off for the next 0.5 ps. Subsequently, attractive hydrogen bond biasing potentials were switched on for 0.5 ps, that were switched off for the last 0.5 ps of the alternating cycle. Both repulsive and attractive potentials had a force constant of 50 kJ/mol. All parameters used for the hydrogen bond potentials are listed in Table 1. Hydrogen Bond Switching Simulations. We performed hydrogen bond switching (HS) simulations for the β-heptapeptide, the Af16 and the Au16 systems. The β-heptapeptide simulation consisted of 5 replicas, while the polyalanine simulations contained 11 replicas. For all three, we used maximum repulsive and attractive force constants of 1 kJ/mol. The scaling factor λ jump between neighboring replicas was 0.5 for the β-heptapep-

6488

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

Vreede et al.

f u tide and 0.2 for both A16 and A16 . During each HS simulation, 41 all-pair exchange attempts were performed every 1 ps. These 1 ps MD runs were performed using GROMACS31 extended with the continuous hydrogen bond biasing potentials, using a time step of 1 fs. The constant volume MD simulations were kept at a temperature of 298 K with the Nose´-Hoover thermostat.37,38 Table 1 lists the parameters used for the repulsive and attractive biasing potentials at their maximum values. All HS simulations also contained a neutral replica, in which no biasing potentials were active. Temperature REMD simulations. To compare the HS method to temperature REMD, we performed a temperature REMD simulation of the β-heptapeptide using 20 replicas. Each replica started from a different conformation, ranging in temperature from 275.34 to 399.58 K.The imposed temperature Ti of each replica i is given by Ti ) T0 exp(ic), where T0 and c can be varied to give the desired exchange ratio.9 We used T0 ) 270 K and c ) 0.0196, similar to the settings in ref 26. We did not perform any further optimization. Exchanges between all replicas were attempted every 1 ps using a homemade Perl script.11 The 1 ps MD runs were performed using GROMACS,31 with a time step of 2 fs. The Nose´-Hoover thermostat imposed a temperature on the replicas.37,38 Analysis. Conformations at λ ) 0 for the HS simulations, or at T ) 298 K for the temperature REMD simulation, were combined into trajectory files. Analysis of these files consisted of tracking several order parameters using the tools included in the GROMACS software package, in combination with Perl scripts. For the β-heptapeptide, we used the number of backbone hydrogen bonds and the positional rmsd of the backbone atoms of residues 2-6 from the NMR structure.42 We also performed a cluster analysis on the trajectories,43 starting by calculating a matrix of positional rmsd between all conformations. Then, the conformation with the most neighbors within a cutoff of 0.07 nm was determined, and together with the neighbors removed from the ensemble. This procedure was repeated until all structures were part of a cluster.43 The order parameters used to investigate the conformational free energy barriers of the A16 system were the backbone R-helical hydrogen bonds between residue n and n + 4, and the contact order:

N

CO )

∑

1 ∆Sij LN i 0 representing λrep.

simulations of 16-residue polyalanine, A16, starting from two initial conformations. The first conformation is a folded, f , and the second R-helical structure, further referred to as A16 structure is an unfolded, extended structure, further denoted as u . We used 11 replicas with a spacing of 0.2 in λ, using a A16 maximum force constant of 1 kJ/mol. With these conditions, the probability of exchanging two neighboring replicas at room temperature is around 0.6, if the difference in the number of hydrogen bonds is 12. u Figure 5a displays the exchange of the replicas for the A16 simulation. The value for λ for each replica is plotted as a function of the simulation time, as a running average of 100 ps to enhance visibility. Each line represents one replica, and most replicas visit each value of λ often. There are two replicas that remain for a large part of the simulation time at extreme values f system. In this of λ. Figure 5b shows the exchange in the A16 simulation, the replicas spend most of the time in either the repulsive or attractive λ-region, as indicated by the white gaps between λ ) -0.2 and λ ) 0.6. There are replicas that cross from repulsive to attractive and vice versa, but most of these replicas quickly return. There is one replica that remains in the repulsive region after crossing for longer than 100 ps; this replica crosses from attractive to repulsive at around 7.5 ns. In developing hydrogen bond switching, we aimed for fast diffusion through the entire range of biasing potentials. In the HS simulations, we observe that every replica visits each value

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6491 of the biasing potential, although on average it prefers a specific bias. Because all replicas traverse the whole range of biasing potentials, they also visit the neutral replica without any biasing potential. Using the conformations from the neutral replica, we calculated free energy profiles, using eq 13, as a function of the number of helical hydrogen bonds and the contact order (see Computational Details, eq 14). Figure 6a and b shows the free energy profiles obtained from, respectively, the Au16 and Af16 f simulations. In the A16 free energy profile, we identified three regions: F(HBn-n+4 ≈ 11, CO ≈ 4.4), P(HBn-n+4 ≈ 7, CO ≈ 6.6), and C(HBn-n+4 ≈ 0, CO ≈ 6.4). We found the u free energy profile, as well as free latter two also in the A16 energy minimum U(HBn-n+4 ≈ 0, CO ≈ 2.4). F represents the fully folded, R-helical conformation, whereas U represents the fully unfolded, extended chain. C and P contain intermediate conformations and respectively represent a collapsed coil and a partially formed R-helix. Figure 6c shows representative conformations of these four regions. It is clear that after 20 ns the two simulations have not converged to identical free energy profiles. Nevertheless, the profiles do overlap substantially. In both profiles, free energy minimum C occurs, although, for Au16, this region covers a larger contact order interval. Partially folded f free energy profile, intermediate P is not a minimum in the A16 u but in the free energy profile of A16, it is. Both the partial helical conformations of state P and the collapsed coils in state C may be part of the structural ensemble of polyalanine A16. Experimentally, the structure of polyalanine in solution, and A16 in particular, could not be determined, as it readily forms aggregates, prohibiting structural determination at high resolution.28 Solid state NMR experiments on a polyalanine stretch, capped by glycine residues, show that a polyalanine stretch can adopt a helical conformation.27 From a physicochemical point of view, the C and P states are favorable compared to the extended conformations in state U and the fully formed R-helical conformations in state F. In both C and P, the solvent accessibility of the hydrophobic regions, comprising R and β carbon atoms, are minimized, and as such, these states maximize the solvent entropy. Several research groups have investigated polyalanine in varying compositions using temperature REMD. Polyalanine A21 required at least 32 replicas ranging from 275 to 475 K.49 Using the temperature distribution webtool,14 polyalanine A16 as used in this work would require 60 replicas ranging from 275 to 475 K for an acceptance probability of 0.25. In the HS scheme, only 11 replicas are required for A16. This is a significant improvement. Sorting of Replicas. Although the hydrogen bond switching scheme requires fewer replicas than temperature REMD, we have encountered a problem in the exchange between replicas that slows down convergence. Once a replica has adjusted to a certain biasing potential, it becomes significantly more difficult to switch and adjust to another biasing potential. This is evident from the fraction of simulation time, the residence, each replica spends at each value of λ, as shown in Figure 7. Ideally, each replica would spend 1/11 (≈0.09) of the simulation time at each u simulation (Figure 7a), 9 of the 11 value of λ. For the A16 replicas show this ideal behavior. Two replicas spend significantly more time at extreme values of λ. The inset in Figure 7a shows the exchange of these two replicas for a short time interval. Surprisingly, the two replicas that remain at high λ-values do exchange with other replicas, even with replicas at opposite λ. Nevertheless, in the next exchange step, they return back to a high value of λ but not necessarily the same value. f simulation (Figure 7b), This effect is even stronger for the A16

6492

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

Vreede et al.

Figure 6. Free energy profiles of HS simulations of (a) the Au16 and (b) the Af16 systems. The profiles are calculated as a function of the number of helical hydrogen bonds and the contact order (eq 14). The contour lines indicate levels of kBT, decreasing with 1 with darker shading. The labels indicate free energy minima. (c) Snapshots of conformations of A16 representing the various free energy minima, shown as a ribbon representation.

where none of the replicas spend equal amounts of time at all λ-values. Instead, each replica favors specific values of λ. For extreme λ-values, the replicas are more trapped than replicas at λ-values close to zero. The inset in Figure 7b shows the exchange steps of the replicas that spend most of their time at extreme λ-values. Most exchanges occur between replicas that are close in λ, but there are exchanges with replicas at λ-values u system, close to zero and even opposite sign. Similar to the A16 the replica returns to an extreme λ-value after such an exchange. Examination of the conformations in the replicas that are stuck at extreme λ-values shows that at high λ (strong repulsive hydrogen bond biasing potential) most conformations are unfolded. At low λ (strong attractive hydrogen bond biasing potential), the conformations are folded. A folded conformation with many hydrogen bonds favors λ-values corresponding to attractive biasing potentials, as its energy is lowest in this region. Similarly, unfolded conformations with little to no hydrogen bonds favor λ-values corresponding to repulsive biasing potentials. There are exchanges that result in unfolded structures at attractive hydrogen bond potentials and vice versa, but in the next exchange attempt, the unfolded replica goes back to a repulsive hydrogen bond potential (and vice versa). After a successful exchange attempt to an opposite biasing potential, the time interval until the next exchange attempt is not long enough to induce conformational changes, i.e., adapt to the new potential. Consequently, the replica is still in a configuration that agrees better with the old biasing potential than with the new biasing potential and therefore the return exchange is highly favorable. Furthermore, a conformational change leading to a structure with a preference for a different biasing potential depends on random “lucky” fluctuations, and are considered rare events. As a result, replicas have a preference for a specific value of the biasing potential, leading to sorting of replicas. Also, in other Hamiltonian REMD schemes, sorting of replicas

occurs. The distribution of replicas visiting the different force fields in ref 23 is skewed, and as such, shows similar behavior as the HS scheme. In solute tempering, a replica exchange method that tunes solvent-solvent and protein-solvent interactions as a function of the temperature, replicas tend to segregate, resulting in exchange occurring only between folded conformations or unfolded conformations.50 Also when using replica exchange umbrella sampling sorting occurs, as the residence time of replicas is not evenly distributed over all umbrella windows.51,52 Convergence is slowed down, as the diffusion of the replicas through the λ values is hampered by sorting, which is a common trait in Hamiltonian REMD schemes that target protein conformational changes. We expect a much faster convergence if this sorting can be prevented, which requires the replicas to adapt to the new potentials within the time between exchange attempts. However, adjusting to the new potential might take a very long time, depending on the system and the force constant. There are two straightforward ways to deal with this problem. The first solution to facilitate faster adaptation to the new potential is increasing the time between exchange attempts. Diffusion through the replicas is then severely slowed by the lower frequency of exchange attempts, and as a consequence, the increase in efficiency is less. The second solution is increasing the force constant. To illustrate this solution, we performed a hydrogen bond switching simulation using a value of 10 kJ/mol for the force constant of the A16unf system. As the acceptance ratio for the hydrogen bond switching simulations with fc ) 1 kJ/mol was very high (80%), we did not include additional replicas, so the λ-spacing remained 0.2. The acceptance ratio in this simulation was 50%. In Figure 8, the migration of replicas through λ-space is plotted as a function of the simulation time. Replicas exchange only in a subset of the values of the biasing potentials, and for this simulation, the

Reordering Hydrogen Bonds Using Hamiltonian REMD

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6493 to very high forces. Correctly addressing those forces requires a shorter time step, leading to a higher computational cost for the simulations. As gaining efficiency is the motivation behind developing these Hamiltonian REMD methods, clearly these straightforward solutions are not very satisfactory. Improving the efficiency of parallel tempering is an active area of research and focuses on adjusting the temperature distribution dynamically,15 identifying entropic bottlenecks16,18,19 or minimizing the first passage time through replica space.17 Using these strategies to optimize the λ-spacing of replicas in hydrogen bond switching might prove very useful to prevent this sorting effect. With an optimized distribution of biasing potentials, Hamiltonian REMD methods may gain an additional increase in efficiency. Another interesting approach would be to use the hydrogen bond biasing potentials in a replica exchange umbrella sampling scheme.52 In contrast to the attractive and repulsive hydrogen bonds described in this work, this alternative formulation aims to restrain the average distance between hydrogen bond donors and acceptors by a harmonic or other suitable potential. The umbrella windows can then be positioned at different values of the average hydrogen bond distance. Short, intermediate, and long average distances will induce folded, partly folded, and unfolded conformations, respectively. As the transition from folded to unfolded will be much smoother in this formulation of hydrogen bond switching, we expect better migration behavior through the range of replicas. Conclusions

Figure 7. Fraction of simulation time, the residence, of each replica at each value of λ, for (a) Au16 and (b) Af16. The insets show each exchange step for replicas that are indicated by the lines. The λ used here is a condensed version of λatt and λrep, with -1 e λ < 0 representing λatt and 0 < λ e 1 representing λrep.

Figure 8. Exchange in the HS simulation of A16fol at a force constant of 10 kJ/mol. For four replicas, λ is plotted as a function of the simulation time. The λ used here is a condensed version of λatt and λrep, with -1 e λ < 0 representing λatt and 0 < λ e 1 representing λrep.

sorting is much more apparent than in the simulations using a force constant of 1 kJ/mol. In general, this solution would require many more replicas to maintain reasonable acceptance probability between the subsets of replicas, thus losing the gain in efficiency. In addition, using very strong potentials will lead

In this work, we present a Hamiltonian REMD scheme, using biasing potentials that either break or form hydrogen bonds in the protein backbone, called hydrogen bond switching (HS). To this end, we developed continuous hydrogen bond biasing potentials based on potentials used in a nonequilibrium alternating scheme aimed at fast folding.7 These potentials still result in conservation of the total energy in the NVE ensemble, and when applied in the alternating scheme, rapid sampling of the conformational space of polyalanine occurs. To address the efficiency of the HS method, we performed a HS and a temperature REMD simulation of a β-heptapeptide in methanol. Both simulation methods resulted in very similar free energy profiles and required 5 ns or less to sample most of the conformational space of the β-heptapeptide. As the HS simulation required only 5 replicas, while the temperature REMD required 20 replicas, the HS method is significantly more efficient. In addition, we performed HS simulations of polyalanine A16 in explicit water, starting from both a fully folded R-helical structure and an extended chain. Both simulations sampled similar intermediate states, representing collapsed coil-like structures and partially folded R-helices. As the HS simulations require only 11 replicas to sample the conformational space of A16, our method presents a significant improvement in efficiency over temperature REMD. When setting up an HS simulation, no knowledge of folded states is required, similar to temperature REMD. To conclude, HS presents an efficient way of sampling protein conformational space by reordering hydrogen bonds in the protein backbone. Our simulations also revealed that replicas have a preference for specific biasing potentials, hampering the diffusion of replicas through λ-space, thereby slowing down convergence. This sorting of replicas also occurs in other Hamiltonian REMD schemes that target protein conformational changes. As straightforward solutions require a transition to much less efficient

6494

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

REMD protocols, future research should focus on a more sophisticated solution to the sorting effect. Preventing this sorting will presumably further increase the efficiency and utilization of Hamiltonian REMD methods in general. Acknowledgment. The authors would like to thank Xavier Periole for providing the β-heptapeptide test system. This work is financially supported by The Netherlands Organisation for Scientific Research. M.G.W. acknowledges support from the CLS Research Project, grant 635.100.012. The authors acknowledge the Stichting Nationale Computerfaciliteiten (NCF) for computer time on the national supercomputer. References and Notes (1) Levinthal, C. How to Fold Graciously. Mossbauer Spectroscopy in Biological Systems: Proceedings of a meeting held at Allerton House, Monticello, Illinois, 1969; pp 22-24. (2) Tanford, C. J. Am. Chem. Soc. 1962, 84, 4240–4247. (3) Fernandez, A.; Kardos, J.; Goto, Y. FEBS Lett. 2003, 536, 187– 192. (4) Fersht, A. R.; Shi, J. P.; Knill-Jones, J.; Lowe, D. M.; Wilkinson, A. J.; Blow, D. M.; Brick, P.; Carter, P.; Waye, M. M. Y.; Winter, G. Nature (London) 1985, 314, 235–238. (5) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Timneanu, N. J. Phys. Chem. B 2006, 110, 4393–4398. (6) Williams, D. H.; Searle, M. S.; Mackay, J. P.; Gerhard, U.; Maplestone, R. A. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 1172–1178. (7) Wolf, M. G.; de Leeuw, S. W. Biophys. J. 2008, 94, 3742–3747. (8) Hansmann, U. H. E. Chem. Phys. Lett. 1997, 281, 140–150. (9) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141–151. (10) Baumketner, A.; Shea, J. E. J. Mol. Biol. 2006, 366, 275–285. (11) Vreede, J.; Hellingwerf, K. J.; Bolhuis, P. G. Proteins 2008, 72, 136–149. (12) Mohanty, S.; Meinke, J.; Zimmerman, O.; Hansmann, U. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 8004–8007. (13) Rathore, N.; Chopra, M.; dePablo, J. J. Chem. Phys. 2005, 122, 02411. (14) Patriksson, A.; van der Spoel, D. Phys. Chem. Chem. Phys. 2008, 10, 2073–2077. (15) Rick, S. J. Chem. Phys. 2007, 126, 054102. (16) Trebst, S.; Troyer, M.; Hansmann, U. J. Chem. Phys. 2006, 124, 174903. (17) Nadler, W.; Hansmann, U. Phys. ReV. E 2007, 76, 065701. (18) Sabo, D.; Meuwly, M.; Freeman, D.; Doll, J. J. Chem. Phys. 2008, 128, 174109. (19) Zheng, W.; Andrec, M.; Gallicchio, E.; Levy, R. J. Phys. Chem. B 2008, 112, 6083–6093. (20) Fukunishi, H.; Watanabe, O.; Takada, S. J. Chem. Phys. 2002, 116, 9058–9067. (21) Mitsutake, A.; Sugita, Y.; Okamoto, Y. J. Chem. Phys. 2003, 118, 6676–6688. (22) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 2000, 329, 261–270.

Vreede et al. (23) Affentranger, R.; Tavernelli, I.; Di Iorio, E. E. J. Chem. Theory Comput. 2006, 2, 217–228. (24) Kwak, W.; Hansmann, U. Phys. ReV. Lett. 2005, 95, 138102. (25) Kannan, S.; Zacharias, M. Proteins 2007, 66, 697–706. (26) Periole, X.; Mark, A. E. J. Chem. Phys. 2007, 126, 014903. (27) Nakazawa, Y.; Asakura, T. J. Am. Chem. Soc. 2003, 125, 7230– 7237. (28) Shinchuk, L. M.; Sharma, D.; Blondelle, S. E.; Reixach, N.; Inouye, H.; Kirschner, D. A. Proteins 2005, 61, 579–589. (29) Frenkel, D. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 17571–17575. (30) Coluzza, I.; Frenkel, D. ChemPhysChem 2005, 6, 1779–1783. (31) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701–1718. (32) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular simulation: The GROMOS96 manual and user guide; Hochschulverlag AG an der ETH Zu¨rich: Zu¨rich, Switzerland, 1996; 1-1042. (33) Daura, X.; Mark, A.; van Gunsteren, W. J. Comput. Chem. 1998, 19, 535–547. (34) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Interaction models for water in relation to protein hydration; D. Reidel Publishing Company: Dordrecht, The Netherlands, 1981; pp 331342. (35) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089– 10092. (36) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463–1472. (37) Hoover, W. G. Phys. ReV. A 1985, 31, 1695–1697. (38) Nose´, S. Mol. Phys. 1986, 52, 187–191. (39) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182–7190. (40) Nose´, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055–1076. (41) Brenner, P.; Sweet, C. R.; von Handorf, D.; Izaguirre, J. A. J. Chem. Phys. 2007, 126, 74103–74106. (42) Seebach, D.; Ciceri, P. E.; Overhand, M.; B., J.; Rigo, D.; Oberer, L.; Hommel, U.; Amstutz, R.; Widmer, H. HelV. Chim. Acta 1996, 79, 2043– 2066. (43) Daura, X.; van Gunsteren, W.; Mark, A. Proteins 1999, 34, 269– 280. (44) Plaxco, K. W.; Simons, K. T.; Baker, D. J. Mol. Biol. 1998, 277, 985–994. (45) Liu, P.; Huang, X.; Zhou, R.; Berne, B. J. J. Phys. Chem. B 2006, 110, 19018–19022. (46) Faraldo-G/’omez, J. D.; Roux, B. J. Comput. Chem. 2007, 28, 1634– 1647. (47) Hritz, J.; Oostenbrink, C. J. Chem. Phys. 2008, 128, 144121. (48) Min, D.; Wei, Y. J. Chem. Phys. 2008, 128, 094106. (49) Garcı´a, A. E.; Sanbonmatsu, K. Y. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 2782–2787. (50) Huang, X.; Hagen, M.; Kim, B.; Friesner, R.; Zhou, R.; Berne, B. J. Phys. Chem. B 2007, 111, 5405–5410. (51) Lou, H.; Cukier, R. J. Phys. Chem. B 2006, 110, 24121–24137. (52) Wolf, M.; Jongejan, J.; Laman, J.; de Leeuw, S. J. Phys. Chem. B 2008, 112, 13493–13498.

JP809641J