PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Chem Theory Comput. Author manuscript; available in PMC 2010 March 30.
Published in final edited form as:
J Chem Theory Comput. 2008 September 9; 4(9): 1541–1554.
doi:  10.1021/ct800086s
PMCID: PMC2847312
NIHMSID: NIHMS184077

Intrinsic Free Energy of the Conformational Transition of the KcsA Signature Peptide from Conducting to Nonconducting State

Abstract

We explore a conformational transition of the TATTVGYG signature peptide of the KcsA ion selectivity filter and its GYG to AYA mutant from the conducting α-strand state into the nonconducting pII-like state using a novel technique for multidimensional optimization of transition path ensembles and free energy calculations. We find that the wild type peptide, unlike the mutant, intrinsically favors the conducting state due to G77 backbone propensities and additional hydrophobic interaction between the V76 and Y78 side chains in water. The molecular mechanical free energy profiles in explicit water are in very good agreement with the corresponding adiabatic energies from the Generalized Born Molecular Volume (GBMV) implicit solvent model. However comparisons of the energies to higher level B3LYP/6–31G(d) Density Functional Theory calculations with Polarizable Continuum Model (PCM) suggest that the nonconducting state might be more favorable than predicted by molecular mechanics simulations. By extrapolating the single peptide results to the tetrameric channel, we propose a novel hypothesis for the ion selectivity mechanism.

Introduction

Organisms transmit electric impulses by means of cellular membrane polarization that critically depends on the work of ion channels. These channels permit passage of specific ion types across the membrane. Ion channels selective for potassium such as KcsA1,2 are particularly interesting as they solve a nontrivial problem of selecting larger K+ over smaller Na+ ions. Despite the wealth of information derived from both experimental14 and computational518 studies of potassium channels, the mechanism of selectivity in these biological machines remains too difficult to tackle as it requires probing the multi-ion permeation transition states.35,7,8,1118 From the computational perspective, this task demands computing multidimensional potentials of mean force (PMFs) for which efficient tools have been lacking.1936

Recently, we have developed and generalized the gradient-augmented Harmonic Fourier Beads (ggaHFB) method2123 that allows studying rare events in complex molecular systems by extending Fukui’s intrinsic reaction coordinate (IRC) approach37,38 with the help of the multidimensional free-energy gradient.22,23,39,40

In the present paper we apply the ggaHFB methodology to study an important functional transition of the signature peptide TATTVGYG of the KcsA selectivity filter that pinches the filter shut by flipping its V76 carbonyl group away from the channel axis coupled with the V76 side chain rotation in response to lowering the K+ concentration.1,3,4 The V76 carbonyl group flip in the KcsA channel is associated with the αL to pII backbone conformational transition at the G77 residue of the signature peptide and is believed to switch the selectivity filter from a conducting (αL) to a nonconducting (pII) state. This transition has been alluded to by X-ray crystallography that detected a partial flip of the V76 backbone carbonyls in the wild type KcsA upon lowering K+ concentration1,2 and recently a more pronounced flip in the E71A mutant.41 Similar transitions have been observed in numerous molecular dynamics (MD) simulations of KcsA7,8 and other related channels.4244 Interestingly, X-ray crystallographic studies indicate that the αL to pII backbone transition is accompanied by rotation of the V76 side chain. However, to the best of our knowledge, previous MD simulations of the KcsA and related potassium channels did not report such a rotation. Furthermore, while the carbonyl flip into pII state observed by X-ray crystallography preserved the 4-fold symmetry of the channel, the MD simulations reported only a single strand out of four identical strands to undergo the αL to pII transition, thus breaking the symmetry of the channel.

It is possible that averaging over the four strands of the filter might artificially diminish the extent of the transition seen by the X-ray crystallography, thus masking the symmetry breaking. However, unambiguous demonstration of the symmetry breaking requires assessing the free energy of the conformational transition in the full tetrameric channel. Although possible to accomplish with the help of the ggaHFB method, this task is computationally intensive as it requires free energy optimization of a transition path ensemble for a relatively large system. On the other hand, exploring the same transition using a single peptide might provide useful insights into the function of the tetrameric channel with reduced computational burden. In particular, the intrinsic free energy profile should provide relative free energies of the αL and pII states along with the corresponding free energy barrier outside the channel environment and thus suggest whether multiple transitions inside the channel are likely.

We define the intrinsic free energy profile of the peptide as that of a single peptide in water. Our choice of water medium has been motivated by the following observations. The distributions of the Ramachandran dihedral angles of various residues in the existing protein structures resemble those from the corresponding adiabatic maps in water but differ markedly from those in gas phase.4549 Even though KcsA is a trans-membrane protein, when fully assembled and in conducting state, water molecules can access the back of the selectivity filter, where they participate in hydrogen bonding with E71 and D80 residues (not present in our model).1,2,15,50 Additional water molecules reach behind the selectivity filter to interact with other residues of the signature peptide in the nonconducting state.2,3 Furthermore, the filter is known to conduct water with and without the ions and hence has a water accessible interior.13,51 Therefore, we feel that the study of the behavior of a single selectivity peptide in water will provide useful insights for understanding the behavior of the same peptide in the tetrameric channel.

This paper is organized as follows. First, we review the ggaHFB methodology for finding minimum adiabatic potential energy paths and minimum free energy transition path ensembles and computing corresponding energy profiles. Combining the ggaHFB transition path ensemble optimization and free-energy evaluation capabilities with the available X-ray structural information, we then explore the intrinsic free energy profile of the signature peptide underlying the flip of the V76 carbonyl from conducting into the nonconducting state.5256 Furthermore, we evaluate the effect of the V76 side chain rotation on the backbone transition. To derive additional support for the functional importance of the specified transition to the KcsA channel, we compare the free energy profile of the wild type peptide to that of the GYG to AYA mutant. Note that a closely related G77A mutant either abolishes the selectivity57 or abrogates the activity of the channel.58 To diffuse any doubts regarding the choice of the water environment for our study, we examine the changes to the functional transition upon removing the peptide from water and placing it into gas phase. Here we fully utilize the ggaHFB capabilities in finding minimum adiabatic potential energy pathways and computing the corresponding energy profiles via the generalized line integral formalism. Finally, we provide some benchmarks to lend credence to the computed energy profiles in water. In particular, we gauge the molecular mechanical (MM) CHARMM22 force field59,60 against a popular Quantum Mechanical (QM) Density Functional Theory model, namely B3LYP6163 with a 6-31G(d) basis set. To account for the solvent contribution, we employ the Generalized Born Molecular Volume (GBMV)47,64 and Polarizable Continuum Model (PCM)6568 with the MM and QM energy functions, respectively.

Methodology

Given the novelty of the employed transition path and path ensemble optimization technique–the generalized gradient augmented Harmonic Fourier Beads method–that makes this study possible, we briefly describe the main points of the method in the following paragraphs.

1. Reactive Coordinate Space (RCS) and Biasing Potential

The generalized gradient-augmented Harmonic Fourier Beads (ggaHFB) method considers an arbitrary system of N atoms described by 3N generalized coordinates Q = (q1, … ,q3N), and, equivalently, by 3N Cartesian coordinates X = (x1, … ,x3N). The method derives the gradient of either adiabatic potential energy or the free energy of the system with respect to a selected subset of S ≤ 3N coordinates q = (q1, … ,qS) that comprise the reactive coordinate space (RCS) by employing either biased optimization or biased molecular dynamics (MD) or Monte Carlo (MC) simulations, correspondingly. The remaining 3N-S degrees of freedom [r with macron] = (qS+1, … ,q3N) comprise the spectator coordinate space (SCS) and do not contribute explicitly to the energy gradient.

The biasing potential is a linear combination of relatively stiff harmonic restraints and applies only to the RCS degrees of freedom centered at a reference configuration q¯b,ref=(q1b,ref,,qSb,ref) :22,23

Vb(q1,,qs;q1b,ref,,qsb,ref)=i=1skib(qiqib,ref)2
(1)

Here superscript b indicates the bias, and kib is the ith coordinate bias force constant. This biasing potential allows deriving the desired energy gradients using a very simple idea described in the following section.

2. Adiabatic Potential Energy Gradient from Biased Optimization

The key idea for computing the energy gradients is most clearly demonstrated on the example of the adiabatic potential energy. Let us add the biasing potential (1) to the total energy of the system U (Q) = U(q,[r with macron]) and then perform potential energy optimization on the modified potential energy surface. Such optimization should reach an equilibrium point at which the forces from the biasing potential that apply only to the S degrees of freedom balance those from the potential energy. Because the forces on the remaining 3N-S degrees of freedom become identically zero due to optimization, the equilibrium point provides the gradient not of the full potential energy but instead of the adiabatic potential energy. Therefore, the biased optimization yields the gradient of the adiabatic potential energy in RCS. The following equations summarize the above.

U(q¯,r¯)qi|q¯=[q¯]b,r¯=[r¯]=V(q¯)qi|q¯=[q¯]b=2kib([qi]bqib,ref),fori=1,,SU(q¯,r¯)qi|q¯=[q¯]b,r¯=[r¯]=0,fori=S+1,,3N
(2)

The square brackets indicate the local minimum on the modified potential energy surface. This procedure effectively reduces the full potential energy surface of 3N degrees of freedom to the adiabatic potential energy surface of S ≤ 3N degrees of freedom.

It is worth noting that in order to compute the adiabatic potential energy gradient on steep slopes in the vicinity of transition states one has to use somewhat stiff springs. Otherwise the minimum on the modified energy surface will slide downhill close to the corresponding minimum on the full energy surface providing little or no information about the transition state region. This remark also applies to the free energy gradient discussed in the next paragraph.

3. Free Energy Gradient From Biased Simulations

The idea used to derive the gradient of the adiabatic potential energy can be applied to derive the gradient of the free energy from biased simulations. For the proof of this statement we refer the reader to the previous work22,23,39,40 and only summarize the results here. It has been demonstrated that for somewhat stiff Cartesian restraint (1) with reference configuration [x with macron]b,ref in RCS, one can compute the corresponding Cartesian free energy gradient via eq 3a.22,23,39,40

Wu(x¯)xi|x¯=x¯b2kib(x¯ibxib,ref)
(3a)

Similarly, for the restraint (1) in generalized coordinates centered at qb,ref the corresponding free energy gradient is given by eq 3b.23

Wu(q¯)qi|q¯=q¯b2kib(q¯ibqib,ref)+kBTln|J(q¯)|qi|q¯=q¯b
(3b)

Here Wu is the unbiased free energy, kB is the Boltzmann constant, T is the simulation temperature, and |J(q)| is the ensemble-reduced Jacobian for the transformation from Cartesian to the generalized coordinates. Note that eq 3a is practically identical to eq 2 for the adiabatic potential energy gradient, where the biased ensemble average left angle bracketqright angle bracketb = (left angle bracketq1right angle bracketb, … ,left angle bracketqSright angle bracketb) replaces the local minimum [q]b configuration. The additional logarithmic Jacobian term on the right-hand side of the generalized gradient expression 3b is the consequence of using Cartesian MD or MC propagators with the nonlinear restraints.23,69 Unlike the case for the adiabatic potential energy gradient, the free energy gradient expression is approximate.

The quality of the free energy gradient depends on the stiffness of the harmonic restraint22,23 and on the quality of the corresponding configuration averages. To achieve the highest quality, one can either run a single very long simulation or run several short simulations and then combine the results into the cumulative average. We prefer the latter approach for accurate free energy calculations as it allows monitoring convergence of the gradient. Specifically, running P batches of short MD or MC simulations of equal length subject to the restraint (1) provides P sets of averaged coordinates or “evolved beads” left angle bracketqright angle bracketb,j = (left angle bracketq1right angle bracketb,j, … ,left angle bracketqSright angle bracketb,j) for a given reference bead, where j is the batch number. These averages could then be easily combined to yield the higher quality cumulative average:

q¯b=1Pj=1Pq¯b,j
(4)

Importantly, the averaged configuration provides the complete free energy gradient in RCS and not just one of its components:

Wu(q¯)|q¯b=(Wu(q¯)q1|q¯b,,Wu(q¯)qs|q¯b)
(5)

This property of the ggaHFB method is a great advantage over the histogram-based free energy estimates that require much larger arrays of simulations to populate multidimensional histograms.22,23,39,40,70,71 Therefore, the ggaHFB method offers a practical alternative to the conventional umbrella sampling simulations with weighted histogram analysis method (WHAM).70,71

The ability to compute the free energy gradient efficiently makes it possible to perform gradient-driven optimization on free energy surfaces and ultimately to find minimum free energy transition path ensembles.

4. Minimum Adiabatic Energy Transition Path

The ggaHFB method as a path finding tool belongs to the class of double-ended reaction path methods that require a reactant and a product state to describe a transition of interest.1926,7279 Importantly, the ggaHFB method finds reaction or transition paths that are invariant with respect to coordinate transformations. The concept of invariant reaction paths, called “intrinsic reaction coordinate” (IRC), has been developed by Fukui for the full potential energy surfaces37,38 and has been further elaborated by many authors since.26,75,8085 In simple terms IRC represents the center curve of the reaction path region that follows the invariant energy gradient.

In particular, in Cartesian coordinates the IRC curve satisfies the following simple condition

U(x¯)=U(x¯)n(x¯)n(x¯)U(x¯)n(x¯)n(x¯)=0
(6)

where n(X) is the curve tangent and [0 with right arrow above] is the null vector.

Importantly, for nonlinear coordinates the direction of the gradient vector has to be corrected using the corresponding contravariant metric tensor G that potentially depends on all 3N degrees of freedom:

G=(gij)=(k=13Nqixkqjxk)
(7)

otherwise different nonlinear coordinate systems will yield different reaction paths for the same stationary points.37,38,8083,85 Thus, to be invariant the transition path curve in nonlinear coordinates must satisfy the following more complicated condition

(GU(Q¯))=GU(Q¯)n(Q¯)n(Q¯)(GU(Q¯))n(Q¯)n(Q¯)=0
(8)

where n(Q) is the curve tangent in the nonlinear coordinates.

Both eqs 6 and 8 apply also to the adiabatic energy surfaces. Because the system of eq 8 is somewhat complicated by the need to compute the metric tensor, the ggaHFB method employs Cartesian coordinates for the path curve optimization instead of the generalized coordinates.

5. Minimum Free Energy Transition Path Ensemble

Using the free energy gradient, the ggaHFB method generalizes the concept of the Fukui’s IRC37,38 to free energy surfaces. In deriving the free energy gradient the SCS degrees of freedom orthogonal to RCS are averaged over, which results in each point in the RCS representing an ensemble. Thus, the ggaHFB method finds continuous curves that connect the provided reactant and product ensembles through a series of transition and intermediate state ensembles. These curves must satisfy the condition that the invariant free energy gradient be tangential to the path curve at any point. In particular, the ggaHFB method uses the straightforward generalization of eq 6 to free energy surfaces in Cartesian coordinates:

Wu(x¯b)=Wu(x¯b)n(x¯)n(x¯)Wu(x¯b)n(x¯)n(x¯)=0
(9)

As noted above, working with nonlinear coordinates requires computing logarithmic Jacobian corrections to the free energy gradient. Furthermore, finding invariant paths requires additional metric tensor corrections.19,74 No such complications arise in Cartesian coordinates, which is why the ggaHFB method employs these coordinates to optimize transition path ensembles.

6. Transition Path Optimization

To optimize a transition path in Cartesian coordinates, we take K unique configurations {Q¯k}k=1,K¯ that gradually progress from the reactant to the product and assign them to a uniform grid {αk=(k1)(K1)}k=1,K¯ with mesh size of Δα = 1/(K−1). If initial configurations Qk = Qk) are unavailable they could be derived via a linear interpolation or by the activated evolution procedure21 that is similar to the growing string method.86 Using these K configurations, we obtain up to K corresponding Fourier amplitudes for each degree of freedom by applying the standard Fourier transform integration with the trapezoidal rule on the grid87

bni=k=1K1(fni,k+fni,k+1)Δα
(10)

where fni,k=[qi(αk)qi(0)(qi(1)qi(0))αk]sin(nπαk).

This procedure globally interpolates between all the K points, yielding a continuous Fourier curve21,88 which is an analytical function of a progress variable α ε [0;1]:

qi(α)=qi(0)+(qi(1)qi(0))α+n=1Kbnisin(nπα)
(11)

We then redistribute the K beads along the path curve such that they conform to a particular metric. Usually, we reposition the beads to make the arc lengths between adjacent beads of equal length in the RCS.

The newly redistributed beads serve as reference beads to compute the corresponding adiabatic potential energy gradients or the free energy gradients via the evolution procedures described in sections 2 and 3. Thus, for each reference bead q¯kref=q¯(αkref) , the evolution returns either the minimized [q¯]kb=([q1]kb,,[qS]kb) or the averaged bead q¯kb=(q1kb,,qSkb) also called the “raw evolved bead”.

The ggaHFB method borrows the idea of redistributing beads along the curve and reparametrizing the curve given the redistributed beads from the string method.7375,86 All the other essential ingredients of the ggaHFB method, such as the multidimensional energy gradient derived on the fly from the harmonic biasing potential and the Fourier representation of both the path and of the corresponding energy gradient (see below) employed in the energy profile integration via generalized line integral as well as optimization strategies, have been obtained from sources independent of the string method despite apparent similarity.2123,3740,79,88,104

In the following discussion, we omit the complementary SCS coordinates for clarity. These coordinates are assumed to be either completely minimized or averaged over and do not explicitly affect either the path or its energy. Optimization implies that the SCS coordinates are passed along either through dynamics restart files or through the complete coordinate files. In addition, it is assumed that the changes in the SCS coordinates between the beads are continuous.

For brevity, we only discuss how to drive optimization of the transition path ensembles and compute the corresponding free energy profiles. The same strategies apply to finding the transition paths on adiabatic potential energy surfaces and computing the corresponding energy profiles. In this case, the adiabatic potential energies could also be calculated exactly for all the points along the path and compared to those computed using the ggaHFB’s generalized line integral formalism.

Substituting the raw evolved beads into eq 3a gives estimates of the free energy gradients for each bead. These gradients are then used in the steepest descent step to generate the “enhanced evolved beads”:

q¯kSD=q¯kb+γkWu(q¯kb)
(12)

Here γk is the parameter that controls the SD step size for the kth bead. In the present paper we use the uniform step size parameter γ for all the beads for simplicity.

Following the Fourier transform of the enhanced beads to obtain new Fourier amplitudes, redistribution of the beads along the resulting curve provides new reference beads. These reference beads are realigned to maintain the coordinate system. For this purpose we invoke a mass-weighted best-fit procedure in a suitable space, usually RCS, to enforce the Eckart conditions on the beads.22,79,8991 In cases where only a few coordinates are available for the best fit or if their geometric arrangement breaks down the standard best fit procedure, simpler alignment methods could be used. The final realigned beads then replace the previous reference beads in the next round of evolution. This procedure is repeated until convergence of the path, i.e. until the path curve changes cease. The final optimized curve represents an invariant minimum free energy transition path ensemble that satisfies the Fukui’s IRC criteria.37,38

The convergence rate of the ggaHFB method depends to some extent on the employed bias force constant and step size parameter. Therefore, devising an optimization strategy to achieve the fastest convergence possible is desirable and is an active area of research in our laboratory.

7. Computing the Free Energy Along the Fourier Path

Given a Fourier path in the generalized multidimensional coordinate space and the corresponding free energy gradients, we can compute the free energy profile along that path via the generalized line integral formalism. To achieve the highest accuracy, we Fourier transform both the evolved beads (4) and the corresponding free energy gradients (5) along the path. With the continuous Fourier representations of the forces and the path, we could then analytically evaluate the corresponding reversible work line integral passing through the evolved beads:

Wu(α)=i=1S0α[Wu(α)qiqi(α)]dα
(13)

In practice, we evaluate the generalized line integral of the second order in eq 13 on a fine uniform grid with L»K quadrature points.

This procedure provides the free energy or the potential of mean force (PMF) profile as an analytical function of the progress variable. Unlike umbrella sampling with WHAM, the interpolation-based ggaHFB free energy integration procedure does not require overlap between the windows. Furthermore, the ggaHFB integration procedure allows natural decomposition of the free energy into contributions from individual coordinates.

The analytical form of the energy profile and that of the corresponding path provided by the ggaHFB method renders pinpointing the energy extrema and their accurate RCS coordinates particularly trivial. One can easily find the values of the progress variable α corresponding to extrema on the energy profile and then substitute these values into eq 11 to get the matching structures.

8. Summary of the ggaHFB Methodology

In summary, the ggaHFB method finds the Fukui’s IRC curves on the adiabatic potential energy surfaces and further generalizes this approach to Cartesian free energy surfaces. Thus, the ggaHFB method finds either minimum adiabatic potential energy paths or minimum free energy transition path ensembles via a gradient driven optimization procedure. Optimizing the transition paths and path ensembles in Cartesian coordinates bypasses the need to calculate the corresponding metric tensors. The optimized transition paths provide structural and energetic information about all the intermediates and transition states connecting given reactants and products at once. Furthermore, the global Fourier representation of the path and the forces provide useful means to control various aspects of the path optimization and ultimately makes the ggaHFB optimization extremely robust.

Independent from the path optimization, the ggaHFB method is a practical alternative to the conventional approach to free energy calculations via umbrella sampling with WHAM. Advantageously, the ggaHFB method is histogram-free, which makes it applicable to cases with arbitrary many dimensions. Even though ggaHFB uses somewhat stiff springs, it does not require the overlap between the windows to integrate the free energy profile. Additionally, the Cartesian version of the ggaHFB method avoids the need to compute the logarithmic Jacobian correction that is required if either WHAM or ggaHFB is used with nonlinear coordinates such as bond distances, angles, dihedrals, etc. to compute free energy profiles.23,92 Finally, the energy profiles can be straightforwardly decomposed into contributions from the individual degrees of freedom that could be useful for analysis and design purposes.

Results

1. Minimum Free Energy Transition Path Ensembles

To explore the free energy of the αL to pII backbone transition of a TATTVGYG signature peptide in water and the effect of the V76 side chain rotation, we use the ggaHFB method with two reactive coordinate spaces (RCSs) of different dimensionalities. Specifically, we include all heavy atoms of the peptide into RCS1 and derive RCS2 from RCS1 by excluding side chain atoms. The RCS1 surface provides the free energy of the backbone configuration subject to a particular side chain orientation. In contrast, the RCS2 surface provides the free energy of the backbone configuration irrespective of the side chains. Unless otherwise stated, throughout this work we employ molecular dynamics in the isothermal isobaric NPT ensemble at 298 K and 1 atm using the CHARMM22 molecular mechanical force field59,60 with the CHARMM-modified TIP3P explicit water model9397 to derive all the required free energy gradients.

Our preliminary free-energy optimization runs revealed that both the αL and pII states at position 77 are local free-energy minima of the isolated peptide in water. Interestingly, the partially flipped, nonconducting conformation observed by X-ray crystallography at low K+ concentration (PDB code 1R3K) is unstable by itself in water despite the rotation of the V76 side chain away from the high K+ concentration, conducting conformation (PDB code 1R3J). During optimization of the peptide from the partially flipped state its backbone, but not the V76 side chain, collapses to the conducting state conformation.

Therefore, to study the full range of the peptide flip, we have constructed an initial path that includes both αL and pII states of G77 backbone. To assess the effect of the V76 side chain rotation, we have included two such rotations by requiring that the end points have the same V76 side chain orientation, matching that of the conducting state. Furthermore, we have inserted the crystallographic nonconducting state with partially flipped backbone and rotated V76 in the middle of the path (refer to the Supporting Information for details).

Performing thorough optimization on the RCS1 free energy surface (see the Supporting Information), we obtain an intrinsic transition path ensemble for the G77 backbone conformational transition from αL to pII state in the TATTVGYG peptide in water. Using the RCS1-optimized path ensemble as the reference we then compute the final free-energy profile, both coupled with (RCS1) and uncoupled from (RCS2) the V76 side chain rotation. To elaborate on the energetics of the backbone transition further, we compute an analogous optimal path and the free-energy information for a GYG to AYA mutant. The resulting cumulative PMFs at different collection times are depicted in Figure 1, and the representative structures of the wild type peptide are shown in Figure 2.

Figure 1
Cumulative PMFs for the conformational transition of the signature TATTVGYG KcsA peptide and its AYA mutant from the αL to pII state in explicit water on RCS1 and RCS2 free-energy surfaces at different collection times. The ggaHFB method employed ...
Figure 2
Representative structures from the free energy transition path ensemble of the wild type TATTVGYG signature peptide of the KcsA selectivity filter in explicit water. The values of the progress variable α provided relate structures to the free ...

The PMFs at the RCS1 level display two events involving the V76 rotation as two sharp peaks with barrier heights ranging from 5 to 8 kcal/mol in both directions. Interestingly, the V76 side chain rotation from the conducting into nonconducting orientation destabilizes the αL state by 2 to 3 kcal/mol, indicating that the hydrophobic interaction between the V76 and Y78 side chains provides additional stabilization of the αL state in the conducting conformation. The αL to pII backbone transition at position 77 follows the second, restoring V76 side chain rotation. In the wild type, the free energy barrier for the backbone transition given a specific orientation of the side chains (the RCS1 PMF) has a forward barrier of 6.0 kcal/mol. The pII state is 2.2 kcal/mol less favorable than the αL state and converts back with a barrier of 3.8 kcal/mol. In sharp contrast, in the mutant the forward barrier is only 0.7 kcal/mol and the pII state is 6.3 kcal/mol more favorable than the α-strand. Restoring the α-strand5254 state in the mutant requires surmounting a high 7.0 kcal/mol free-energy barrier.

Setting the side chains free (the RCS2 PMF) permits evaluating the free energy of the backbone transition alone. As is seen from Figure 1, switching to RCS2 space collapses the sharp peaks (labeled with arrows) corresponding to the V76 side chain rotations but leaves the portion of the PMF underlying the backbone transition from the αL to the pII state virtually unchanged. In particular, for the wild type peptide the forward activation barrier is 5.9 kcal/mol, and the pII state is still less stable than the αL by slightly smaller 1.7 kcal/mol. Restoring the conducting state requires over-coming a slightly higher barrier of 4.2 kcal/mol. In contrast, the mutant exhibits a forward barrier of 0.9 kcal/mol and the relative pII state stabilization energy of 7.0 kcal/mol that makes the reverse barrier increase to 7.9 kcal/mol.

2. Minimum Adiabatic Potential Energy Paths

To explicitly evaluate the effect water has on the conformational transitions of the signature peptide and to further demonstrate the capabilities of the ggaHFB method, we have computed the minimum adiabatic potential energy paths for the wild type TATTVGYG peptide and its GYG to AYA mutant in gas phase. Both peptides have three threonine and one tyrosine residues with rotatable OH bonds that were averaged over in the minimum free energy transition path ensembles computed in water. The orientation of these hydrogens significantly perturbs the overall potential energy; therefore, we include these four hydrogen atoms in the reactive coordinate space. Thus, by adding the polar hydrogen atoms to the all-heavy-atom RCS1 we derive the RCS1h for adiabatic potential energy path optimization.

We have to assign some initial values to the tyrosine and threonine OH groups, which have two and three rotameric states, respectively, in order to compute the adiabatic potential energy paths. The total number of possible initial path configurations is therefore 21 × 33 = 54. To control the configurations, we follow the Protein Data Bank atom naming convention and use dihedral angles Cε1-Cζ-Oη-Hη and Cα-Cζ-Oγ-Hγ for the tyrosine and the threonines, respectively. Here, we arbitrarily choose dihedral angles of 180, 180, −30, and 0 degrees for the T72, T74, T75, and Y78, respectively, as the initial conditions for the path optimization. To prepare the initial path with these conditions, we fix the RCS1 coordinates and apply stiff harmonic restraints of 1000 kcal/(mol · rad2) on the corresponding dihedral angles during an optimization of the hydrogen positions.

Because the optimization on adiabatic potential energy surfaces with the bare CHARMM22 molecular mechanical force field is relatively inexpensive, we have initiated the ggaHFB optimization using 89 beads. Using 89 beads is sufficient to integrate the adiabatic potential energy, given the initial orientation of the four OH bonds. Nevertheless, optimization of the OH groups requires increasing the number of beads further to correctly integrate the adiabatic potential energy. The increase reflects the fact that rotations of the OH groups correspond to small changes in the RCS1h, resulting in very sharp transitions along the path. Although proper integration could be achieved by locally increasing the number of beads at the sharp transitions leading to nonuniform bead distributions,23 in the present work we use uniform grid for simplicity. Thorough path optimization increases the overall path length dramatically, in the end requiring 705 beads to properly integrate the adiabatic potential energy along the path.

The final paths in the gas phase have little if any resemblance with the paths optimized in water and exhibit a greater number of local minima and transition states. For the wild type peptide, the α-strand5254 disappears almost completely. First, the G79 residue spontaneously flips into the C5 conformation and then converts into the C7ax conformation. In the flipped configuration on the reactant side of the path G79 forms a hydrogen bond with the OH group of T71 using its carbonyl oxygen. The G79 residue flip significantly perturbs the rest of the α-strand, which quickly collapses further residue-by-residue along the path.

In the mutant, the α-strand is annihilated completely in the reactant basin, where the A79 along with the A77 residues flip into the C7ax conformation. The four residues V76, A77, Y78, and A79 surround the T71 residue like a belt, with alternating axial and equatorial configurations, namely C7ax, C7eq, C7ax, and C7eq, respectively. During the optimization the mutant pathway deviates substantially from that of the wild type.

Figure 3 depicts the corresponding adiabatic energy profiles that underscore the complexity of the changes in the gas phase. It also provides the benchmarks for the adiabatic potential energy integration via the generalized line integral formalism. In particular, comparison of the line integral energies with the exact adiabatic potential energies from the CHARMM22 force field shows the accumulated errors of 0.07 and 0.12 kcal/mol for the wild type and mutant adiabatic energy profiles, respectively. We consider this a very good agreement between the generalized line integral energy and the exact energy profiles.

Figure 3
Adiabatic potential energy profiles along the optimized reaction paths of the signature TATTVGYG KcsA peptide and its AYA mutant in gas phase. The ggaHFB method employed 705 beads to integrate the potential energy with the RCS1h (see text for details) ...

The V76 side chain rotations (labeled with arrows) have been preserved in both wild type and mutant paths, although in some cases they have been coupled with other structural rearrangement as seen in Figure 3. The forward and reverse barrier heights for the V76 side chain rotation vary but are similar to those in water.

Overall the gas phase structures are more compact than the ones in water and establish as many intramolecular hydrogen bonds as possible. Given the complexity of the adiabatic paths and their divergence from the structures obtained by either the X-ray crystallography or by the free energy optimization in water, we omit a detailed description of the structural changes along the path and simply provide the corresponding trajectories in Supporting Information.

3. Comparison of the MM and QM Energy Profiles

A. Gas Phase

Because the present paper investigates an important conformational transition of the TATTVGYG signature peptide from the KcsA potassium channel, it might be useful to assess the molecular mechanical (MM) force field employed. Of particular interest is evaluating the energetics of the signature peptide and its mutant along the path optimized in water. To establish useful benchmarks, we first compute the gas phase adiabatic energy profiles along the minimum free energy transition path ensembles in the RCS1. In particular, we compare the MM energy profiles with one of the most popular density functional theory models, namely B3LYP, with the 6-31G(d) basis set as a high-level quantum mechanical (QM) model (see the Supporting Information for details). This model is not expected to produce accurate energy profiles when it comes to dispersion interactions between the Y78 residue and the V76 side chains and hence should be used with caution.98101

Figure 4 shows the corresponding adiabatic potential energy profiles for the peptides with the rotatable OH bonds fixed at the conformation used as initial condition for the path reoptimization in gas phase. Interestingly, the energy profiles obtained with MM and the QM models for the same path differ substantially. The V76 side chain rotation barriers appear reduced in the QM model.

Figure 4
Gas phase single point energy profile along the α-strand to pII state conformational transition of the signature TATTVGYG KcsA peptide and its AYA mutant in water. MM - uses the bare CHARMM22 force field, QM - B3LYP/6-31G(d) Density Functional ...

Both the MM and QM models favor the pII state over the α-strand. The QM model predicts the α-strand to be much less stable than the pII state in the wild type but relatively more stable in the mutant peptide. In contrast, the MM model suggests that the α-strand is much more stable in the wild type peptide, not the mutant. It is likely that the adiabatic energy surfaces of the MM and QM models are significantly different in the gas phase, and such single point energy profile comparisons should be taken with caution.

B. Implicit Solvent

As mentioned above, we are primarily interested in the energetics of the peptides in water and not in the gas phase. After all, the transition path ensembles for the functional conformational transitions of the peptides have been optimized using the MM model in the explicit water. To compare the MM with QM models in water, we choose the Generalized Born Molecular Volume (GBMV)47,64 and Polarizable Continuum Model (PCM)6568 implemented in Gaussian 03102 implicit solvent models, respectively. Using the free-energy optimized transition path ensembles in explicit water as reference profiles, we have performed optimization of all the degrees of freedom orthogonal to the RCS1 and subject to the same dihedral restraints on the rotatable OH bonds as discussed above with the MM-GBMV model and then computed single point QM-PCM energies for all the beads along the path (see the Supporting Information for details). The results are provided in Figure 5.

Figure 5
Implicit solvent single point energy profile along the α-strand to pII state conformational transition of the signature TATTVGYG KcsA peptide and its AYA mutant in water. MM-GBMV - uses the CHARMM22 force field with the Generalized Born Molecular ...

The MM-GBMV model yields the free-energy profiles that are in very good agreement with the explicit solvent calculations (Figure 1), thus further validating the intrinsic free energy profiles. We cannot expect better agreement between the two profiles given the fixed conformations of the rotatable OH bonds necessary to compute the adiabatic energy profile with GBMV.

The QM-PCM model produces energy profiles very different from those of the MM-GBMV model, and most importantly does not favor the αL state over the pII state in the wild type peptide. In the mutant, on the other hand, the αL state remains unstable with the QM-PCM model.

Because we have performed optimization on the water modified RCS1 free energy surfaces of the peptides with the MM model and explicit water any comparison with other surfaces that have not been optimized do not warrant good agreement, unless the surfaces are exactly the same. This conflict could in principle be resolved by the QM-PCM optimization of the product, reactant, and a few key intermediates, which unfortunately presents a significant challenge at present.

Discussion

The finding that the signature peptide taken from the partially flipped nonconducting state (PDB code 1R3K) collapses back into the conducting state in water despite the V76 side chain rotation shows the intrinsic width of the peptide αL basin. Furthermore, it suggests that either the channel provides additional interactions to stabilize the partially flipped backbone structure or that only one of the four strands of the tetrameric channel undergoes the full transition. If the latter symmetry breaking were to occur, the apparent configuration observed by X-ray crystallography would correspond to the average over four strands, thus artificially reducing the extent of the transition in a single peptide.

The fact that the free energy profiles for the αL to the pII transition are relatively insensitive to the position of the side chains (see Figure 1) reflects the robustness of the backbone transition. The forward and reverse free energy activation barriers in the wild type peptide are 5.9 and 4.2 kcal/mol. Interestingly, previous calculations of an even lower dimensional PMF for a similar transition inside the wild type KcsA channel in the presence of two ions gave a rough estimate of the free energy barrier between 0.5 and 4.0 kcal/mol.8 In sharp contrast, the mutant exhibits a forward barrier of 0.9 kcal/mol and the reverse barrier of 7.9 kcal/mol.

The width of the αL basin in the intrinsic free energy profile of a single peptide might determine the range of a local dilation/contraction of the tetrameric KcsA channel at the V76 carbonyl ring. If the V76 carbonyl were pushed away from the channel axis beyond the limits of the αL basin, the peptide would go over the transition barrier and into the nonconducting pII state. We emphasize that by local dilation/contraction of the channel we imply the change in the distance between the V76 carbonyls associated with the backbone motions within the bounds of the αL basin and not with the transition from αL to pII or back.

The full αL to pII transition has been demonstrated to be unnecessary for the ion selectivity, at least in a synthetic channel with the d-Ala residue in place of G77.3 Note that the wild type KcsA channel, in addition to K+, permits ions of larger size, namely Cs+ and Rb+, which are expected to pass the V76 carbonyl ring without triggering the transition from αL to pII.4 Such a wide range of dilation/contraction would not have been possible if G77 was substituted for regular Ala, as the width and the depth of the αL basin would have been dramatically reduced as seen from the PMFs for the mutant depicted in Figure 1. Note however, that the AYA mutant would be sterically prevented from forming the conducting α-strand conformation in the tetrameric channel.3

In an effort to validate the results obtained with the MM force field in explicit water we have profiled the energetics along the paths using MM and QM methods both in gas phase and in implicit solvent. The results of these calculations are summarized in Figure 4 and Figure 5 that highlight the stark disagreement between the MM and QM models. Although QM models usually have higher fidelity than MM models, the particular DFT method used in this work, namely the B3LYP functional, is well-known to fail to account for short-range dispersion interactions necessary to properly describe the energetics of the hydrophobic interactions such as those between V76 and Y78 side chains.98101 Higher level, more expensive ab initio methods that properly account for the dispersion suggest that the interactions between the CH bonds of the V76 and the phenol ring of the Y78 could favor the α-strand by about 1 kcal/mol.101,103 Additional discrepancies between the MM and QM in this work may arise due to the fact that no optimization has been performed at the QM level of theory. Therefore, the differences between the QM and MM models should be interpreted with caution.

It appears that the stability of the pII state in the wild type peptide might be overestimated by the QM model with implicit water, because it would require at least 20 kcal/mol to assume the conducting conformation in the tetrameric channel. On the other hand the MM model with implicit water predicts the αL and pII configurations to be nearly degenerate. If the QM-PCM model more accurately reproduced the energetics of the solvated peptide even without optimization, the ground or resting state of the wild type KcsA channel would be the nonconducting pII state. Thus the channel would have to be activated by a conformational change from the pII resting state into the α-strand state to conduct ions. This would only be possible due to a strong perturbation such as strong attraction of the ions in the lumen of the tetrameric channel to its carbonyl oxygens.

The switching between the nonconducting and conducting state and the functional contraction/dilation of the tetrameric KcsA channel would require a certain balance between the electrostatic repulsion of the V76 carbonyls and the free energy of the backbone rotation of the residue at position 77. Because the electrostatic repulsion can be relaxed by transiently flipping one or more of the carbonyls out from the αL into the pII state, the filter must also ensure to favor the αL over the pII state at least in the presence of ions in the lumen of the filter. The potassium channel seems to have achieved the α-strand stabilization by using a G residue that has a high propensity for the αL configuration at position 77 and in addition by the hydrophobic interaction between V76 and Y78 side chains. The importance of the hydrophobic interaction is supported by the experimental observation that the V76A mutant abrogates tetrameric assembly of the channel.58 Taking the above into consideration, it appears that the free energy profiles computed with the MM model in explicit water agree better with the proposal than the corresponding single point energy profiles obtained with B3LYP/6-31G(d)-PCM QM model.

In the absence of the actual tetrameric channel in our model, the bulk water better reproduces the environment of the KcsA selectivity filter than the gas phase and therefore provides useful insights into the channel function. In particular, the differences between the adiabatic energy maps of the peptide residues in water and gas phase suggest that the gas phase transition pathways must deviate strongly from those of the transition path ensembles optimized in water. This is particularly true of the RL region that is forbidden in gas phase.48,49 The ggaHFB optimization of the adiabatic paths in gas phase explicitly demonstrates that water plays an active and important role in defining the intrinsic path and the energetics of the peptide backbone transition.

The outcome of the gas phase optimization can be predicted based on the previous studies of the glycine and alanine dipeptides.48,49 In particular, the referred work demonstrated that the αL configurations collapse into the C7ax, while pII configurations collapse into the C7eq.48,49 These are the exact changes we observe upon the adiabatic potential energy optimization in gas phase. The final optimized paths in gas phase are rather complex (see Figure 3) and seem irrelevant for the functional transition of the selectivity peptide in the KcsA channel. On the other hand the pathways in water show very good qualitative agreement with the peptide conformations observed in the tetrameric channel.

Finally, based on the present findings, we are able to propose a novel hypothesis for the mechanism of ion selectivity in the tetrameric KcsA channel. Specifically, we conjecture that in its conducting α-strand state the carbonyl rings should contract around the ion entering the channel and that this contraction would propagate to the nearby carbonyl rings along the channel axis (see Figure 1S in the Supporting Information for an illustration). Because ions are believed to pass the KcsA filter stripped of all but two water molecules that cotranslate with the ion while hydrogen bonding to the carbonyls, these water molecules will experience greater difficulty to pass neighboring carbonyl rings due to the contraction, in turn impeding the ion movements along the channel.

This hypothesis could explain why the channel selects larger K+ over smaller Na+ ions. Specifically, we anticipate that smaller Na+ ions would contract the carbonyl rings to a greater extent than the larger K+ ions thus impeding the passage of the cotranslating water molecules to a greater degree. With water passage impeded the ions themselves must in turn slow down.

Our hypothesis suggests that in the absence of ions the KcsA channel should stay open to water permeation unless one of the four V76 carbonyls flips out pinch-shutting the channel. Indeed, the KcsA channel has been experimentally demonstrated to conduct water in the absence of permeating ions.51 The partial flipping of the carbonyls (while still within the αL basin) might serve a selectivity purpose, whereas a complete flip (transition into the pII basin) can be used to gate the channel.8 To provide further support of this hypothesis we are currently performing an optimization of several transition path ensembles for ion–water copermeation through the tetrameric KcsA selectivity filter. The results of this work will be reported in a forthcoming publication.

To conclude, we have explored an intrinsic free energy landscape of an important functional transition of the signature peptide from the KcsA selectivity filter that is responsible for locally dilating/contracting the channel at the V76 carbonyl ring, in addition to switching the channel between conducting and nonconducting states. We have found that the wild type peptide intrinsically favors the conducting state due to the combination of the high G77 backbone propensity for the αL configuration and the stabilizing hydrophobic interaction between V76 and Y78 side chains. In sharp contrast, the mutant strongly favors the nonconducting state. However, additional steric effects in the tetrameric channel that are absent in the present study are expected to prevent formation of the conducting conformation in the mutant.

We have found the αL to pII transition to be exceptionally robust and intrinsically funneled toward the conducting state in the wild type KcsA peptide at the MM level with explicit water. Although the intrinsic free energy profiles have been validated using the MM with an implicit water model, efforts to gauge performance of the MM model against the QM model indicated that our results should be interpreted with caution. Based on the QM-PCM model it may be possible that the ground state of the channel in the absence of ions could in fact be the nonconducting state and that the conducting state would only form upon ion entrance into the lumen of the channel. Nevertheless, the present study has allowed us to propose a novel hypothesis for the ion selectivity within the KcsA channel in which local contraction of the channel interior in response to the ion presence regulates copermeation of water through the channel to a degree that is inversly proportional to the ion size. Work is currently underway in our laboratory to test the proposed hypothesis in the tetrameric KcsA channel model. We hope that the present work will stimulate future transition path ensemble studies of rare events in complex molecular systems.

Supplementary Material

supporting information

Acknowledgment

This work was in part supported by grants from NIH, NSF, NBCR, CTBP, and HHMI. Computer time on clusters at CTBP, NBCR, and the Department of Bio-Engineering at UCSD is greatly appreciated. M.F. was supported in part by an NIH Molecular Biophysics Training Grant: GM-08326, and subsequently by an NSF Graduate Research Fellowship.

Footnotes

Supporting Information Available: Peptide coordinates from the optimized transition path ensembles together with details of the transition path optimization and PMF calculations, a cartoon illustration of the ion selectivity hypothesis, and details of the B3LYP/6-31G(d)-PCM calculations. This material is available free of charge via the Internet at http://pubs.acs.org.

References

1. Zhou Y, MacKinnon R. The occupancy of ions in the K+ selectivity filter: charge balance and coupling of ion binding to a protein conformational change underlie high conduction rates. J. Mol. Biol. 2003;333(5):965–975. [PubMed]
2. Zhou Y, Morais-Cabral JH, Kaufman A, MacKinnon R. Chemistry of ion coordination and hydration revealed by a K+ channel-Fab complex at 2.0 Å resolution. Nature. 2001;414:43–48. [PubMed]
3. Valiyaveetil FI, Leonetti M, Muir TW, MacKinnon R. Ion Selectivity in a Semisynthetic K+ Channel Locked in the Conductive Conformation. Science. 2006;314(5801):1004–1007. [PubMed]
4. Lockless SW, Zhou M, MacKinnon R. Structural and Thermodynamic Properties of Selective Ion Binding in a K+ Channel. PLoS Biol. 2007;5(5):1079–1088. [PMC free article] [PubMed]
5. Berneche S, Roux B. Molecular dynamics of the KcsA K+ channel in a bilayer membrane. Biophys. J. 2000;78(6):2900–2917. [PubMed]
6. Aqvist J, Luzhkov V. Ion permeation mechanism of the potassium channel. Nature. 2000;404:881–884. [PubMed]
7. Berneche S, Roux B. Energetics of ion conduction through the K+ channel. Nature. 2001;414:73–77. [PubMed]
8. Berneche S, Roux B. A gate in the selectivity filter of potassium channels. Structure. 2005;13:591–600. [PubMed]
9. Mashl RJ, Tang Y, Schnitzer J, Jacobsson E. Hierarchical Approach to Predicting Permeation in Ion Channels. Biophys. J. 2001;81(5):2473–2483. [PubMed]
10. Guidoni L, Carloni P. Potassium permeation through the KcsA channel: a density functional study. Biochim. Biophys. Acta. 2002;1563(1–2):1–6. [PubMed]
11. Compoint M, Boiteux C, Huetz P, Ramseyer C, Girardet C. Role of water molecules in the KcsA protein channel by molecular dynamics calculations. Phys. Chem. Chem. Phys. 2005;7(24):4138–4145. [PubMed]
12. de Haan HW, Tolokh IS, Gray CG. Nonequilibrium molecular dynamics calculation of the conductance of the KcsA potassium ion channel. Phys. Rev. E. 2006;74(3):030905. [PubMed]
13. Grabe M, Bichet D, Qian X, Jan YN, Jan LY. K+ channel selectivity depends on kinetic as well as thermodynamic factors. Proc. Natl. Acad. Sci. U.S.A. 2006;103(39):14361–14366. [PubMed]
14. Kraszewski S, Boiteux C, Langner M, Ramseyer C. Insights into the origins of the barrier-less knock-on conduction in the KcsA channel: molecular dynamics simulations and ab initio calculations. Phys. Chem. Chem. Phys. 2007;9(10):1219–1225. [PubMed]
15. Domene C, Vemparala S, Furini S, Sharp K, Klein ML. The role of conformation in ion permeation in a K+ channel. J. Am. Chem. Soc. 2008;130(11):3389–3398. [PubMed]
16. Gwan J-F, Baumgaertner A. Cooperative transport in a potassium ion channel. J. Chem. Phys. 2007;127(4):045103:1–045103:10. [PubMed]
17. Noskov SY, Berneche S, Roux B. Control of ion selectivity in potassium channels by electrostatic and dynamic properties of carbonyl ligands. Nature. 2004;431:830–834. [PubMed]
18. Noskov SY, Roux B. Ion selectivity in potassium channels. Biophys. Chem. 2006;124(3):279–291. [PubMed]
19. Maragliano L, Vanden-Eijnden E. On-the-fly string method for minimum free energy paths calculation. Chem. Phys. Lett. 2007;446(1–3):182–190.
20. Wienan E, Ren W, Vanden-Eijnden E. Simplified and improved string method for computing the minimum energy paths in barrier-crossing events. J. Chem. Phys. 2007;126(16):164103. [PubMed]
21. Khavrutskii IV, Arora K, Brooks CL., III Harmonic Fourier Beads Method for Studying Rare Events on Rugged Energy Surfaces. J. Chem. Phys. 2006;125(17):174108. [PubMed]
22. Khavrutskii IV, McCammon JA. Generalized gradient-augmented harmonic Fourier beads method with multiple atomic and/or center-of-mass positional restraints. J. Chem. Phys. 2007;127(12):124901. [PubMed]
23. Khavrutskii IV, Dzubiella J, McCammon JA. Computing Accurate Potentials of Mean Force in Electrolyte Solutions with the Generalized Gradient-Augmented Harmonic Fourier Beads Method. J. Chem. Phys. 2008;128(4):044106. [PubMed]
24. Hu H, Lu Z, Parks JM, Burger SK, Yang W. Quantum mechanics/molecular mechanics minimum free-energy path for accurate reaction energetics in solution and enzymes: Sequential sampling and optimization on the potential of mean force surface. J. Chem. Phys. 2008;128(3):034105. [PubMed]
25. Burger SK, Yang W. Sequential quadratic programming method for determining the minimum energy path. J. Chem. Phys. 2007;127(16):164107. [PubMed]
26. Burger SK, Yang W. Quadratic string method for determining the minimum-energy path based on multiobjective optimization. J. Chem. Phys. 2006;124(5):054109. [PubMed]
27. Barducci A, Bussi G, Parrinello M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008;100(2):020603/1–020603/4. [PubMed]
28. Branduardi D, Gervasio FL, Parrinello M. From A to B in free energy space. J. Chem. Phys. 2007;126(5):054103/1–054103/10. [PubMed]
29. Laio A, Parrinello M. Computing free energies and accelerating rare events with metadynamics. Lecture Notes in Physics. 2006;703:315–347. (Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology, Volume 1)
30. Bussi G, Laio A, Parrinello M. Equilibrium Free Energies from Nonequilibrium Metadynamics. Phys. Rev. Lett. 2006;96(9):090601/1–090601/4. [PubMed]
31. Laio A, Rodriguez-Fortea A, Gervasio FL, Ceccarelli M, Parrinello M. Assessing the Accuracy of Metadynamics. J. Phys. Chem. B. 2005;109(14):6714–6721. [PubMed]
32. Laio A, Parrinello M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002;99(20):12562–12566. [PubMed]
33. van der Vaart A, Karplus M. Minimum free energy pathways and free energy profiles for conformational transitions based on atomistic molecular dynamics simulations. J. Chem. Phys. 2007;126(16):164106. [PubMed]
34. Babin V, Roland C, Darden TA, Sagui C. The free energy landscape of small peptides as obtained from metadynamics with umbrella sampling corrections. J. Chem. Phys. 2006;125(20):204909. [PMC free article] [PubMed]
35. Babin V, Roland C, Sagui C. Adaptively biased molecular dynamics for free energy calculations. J. Chem. Phys. 2008;128(13):134101. [PubMed]
36. Huber T, Torda AE, van Gunsteren WF. Local elevation: A method for improving the searching properties of molecular dynamics simulation. J. Comput.-Aided Mol. Des. 1994;8(6):695–708. [PubMed]
37. Fukui K. The Path of Chemical Reactions - The IRC Approach. Acc. Chem. Res. 1981;14(12):363–368.
38. Fukui K, Kato S, Fujimoto H. Constituent Analysis of the Potential Gradient Along a Reaction Coordinate. Method and an Application to CH4+T Reaction. J. Am. Chem. Soc. 1975;97(1):1–7.
39. Hummer G, Szabo A. Free Energy Surfaces from Single-Molecule Force Spectroscopy. Acc. Chem. Res. 2005;38(7):504–513. [PubMed]
40. Kaestner J, Thiel W. Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: “Umbrella integration” J. Chem. Phys. 2005;123(14):144104. [PubMed]
41. Cordero-Morales JF, Cuello LG, Zhao Y, Jogini V, Cortes DM, Roux B, Perozo E. Molecular determinants of gating at the potassium-channel selectivity filter. Nat. Struct. Mol. Biol. 2006;13(4):311–318. [PubMed]
42. Domene C, Grottesi A, Sansom MS. Filter Flexibility and Distortion in a Bacterial Inward Rectifier K+ Channel: Simulation Studies of KirBac1.1. Biophys. J. 2004;87(1):256–267. [PubMed]
43. Capener CE, Proks P, Ashcroft FM, Sansom MS. Filter Flexibility in a Mammalian K+ Channel: Models and SImulations of Kir6.2 Mutants. Biophys. J. 2003;84(4):2345–2356. [PubMed]
44. Khalili-Araghi F, Tajkhorshid E, Schulten K. Dyanmics of K+ Ion Conduction through Kv1.2. Biophys. J. 2006;91(6):L72–L74. [PubMed]
45. Hovmoeller S, Zhou T, Ohlson T. Conformations of aminoacids in proteins. Acta Crystallogr. 2002;D58:768–776. [PubMed]
46. MacKerell AD, Jr, Feig M, Brooks CL., 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(11):1400–1415. [PubMed]
47. Lee MS, Feig M, Salsbury FR, Jr, Brooks CL., III New analytic approximation to the standard molecular volume definition and its application to generalized Born calculations. J. Comput. Chem. 2003;24(11):1348–1356. [PubMed]
48. Pettitt BM, Karplus M. The potential of mean force surface for the alanine dipeptide in aqueous solution: a theoretical approach. Chem. Phys. Lett. 1985;121(3):194–201.
49. Lau WF, Pettitt BM. Conformations of the Glycine Dipeptide. Biopolymers. 1987;26(11):1817–1831.
50. Long SB, Tao X, Campbell EB, MacKinnon R. Atomic structure of a voltage-dependent K+ channel in a lipid membrane-like environment. Nature. 2007;450(7168):376–383. [PubMed]
51. Saparov SM, Pohl P. Beyond the diffusion limit: Water flow through the empty bacterial potassium channel. Proc. Natl. Acad. Sci. U.S.A. 2004;101(14):4805–4809. [PubMed]
52. Milner-White EJ, Watson JD, Qi G, Hayward S. Amyloid Formation May Involve alpha- to beta-Sheet Interconversion via Peptide Plane Flipping. Structure. 2006;14(9):1369–1376. [PubMed]
53. Armen RS, DeMarco ML, Alonso DOV, Daggett V. Pauling and Corey’s {alpha}-pleated sheet structure may define the prefibrillar amyloidogenic intermediate in amyloid disease. Proc. Natl. Acad. Sci. U.S.A. 2004;101(32):11622–11627. [PubMed]
54. Armen RS, Alonso DOV, Daggett V. Anatomy of an Amyloidogenic Intermediate: Conversion of beta-Sheet to alpha-Sheet Structure in Transthyretin at Acidic pH. Structure. 2004;12(10):1847–1863. [PubMed]
55. JiJi RD, Balakrishnan G, Hu Y, Spiro TG. Intermediacy of Poly(L-proline) II and beta-Strand Conformations in Poly(L-lysine) beta-Sheet Formation by Temperature-Jump/UV Resonance Raman Spectroscopy. Biochemistry. 2006;45(1):34–41. [PMC free article] [PubMed]
56. Asher SA, Mikhonin AV, Bykov S. UV Raman Demonstrates that α-helical Polyalanine Peptides Melt to Polyproline II Conformations. J. Am. Chem. Soc. 2004;126(27):8433–8440. [PubMed]
57. Heginbotham L, Lu Z, Abramson T, MacKinnon R. Mutations in the K+ Channel Signature Sequence. Biophys. J. 1994;66:1061–1067. [PubMed]
58. Splitt H, Meuser D, Borovok I, Betzler M, Schrempf H. Pore mutations affecting tetrameric assembly and functioning of the potassium channel KcsA from Streptomyces lividans. FEBS Lett. 2000;472:83–87. [PubMed]
59. MacKerell AD, Jr, Bashford D, Bellott M, Jr, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Reiher WE, III, Roux B, Schlenkrich M, Smith JC, Stote R, Straub J, Watanabe M, Wiorkiewicz-Kuczera J, Yin D, Karplus M. All-atom empirical potential for molecular modeling and dynamics Studies of proteins. J. Phys. Chem. B. 1998;102(8):3586–3616. [PubMed]
60. MacKerell AD, Jr, Brooks BR, Brooks CL, III, Nilsson L, Roux B, Won Y, Karplus M. CHARMM: The Energy Function and Its Parameterization with an Overview of the Program. In: Schleyer PvR, Schreiner PR, Allinger NL, Clark T, Gasteiger J, Kollman P, Henry F, Schaefer I., editors. The Encyclopedia of Computational Chemistry. Chichester: John Wiley & Sons; 1998.
61. Becke AD. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A. 1988;38(6):3098–3100. [PubMed]
62. Becke AD., III The role of exact exchange. J. Chem. Phys. 1993;98(7):5648–5652.
63. Lee C, Yang W, Parr RG. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B. 1988;37(2):785–789. [PubMed]
64. Lee MS, Salsbury FR, Jr, Brooks CL., III Novel generalized Born methods. J. Chem. Phys. 2004;116(24):10606–10614.
65. Cossi M, Barone V, Mennucci B, Tomasi J. Ab initio study of ionic solutions by a polarizable continuum dielectric model. Chem. Phys. Lett. 1998;286(3–4):253–260.
66. Cances E, Mennucci B, Tomasi J. A new integral equation formalism for the polarizable continuum model: theoretical background and applications to isotropic and anisotropic dielectrics. J. Chem. Phys. 1997;107(8):3032–3041.
67. Mennucci B, Tomasi J. Continuum solvation models: A new approach to the problem of solute’s charge distribution and cavity boundaries. J. Chem. Phys. 1997;106(12):5151–5158.
68. Cossi M, Scalmani G, Rega N, Barone V. New developments in the polarizable continuum model for quantum mechanical and classical calculations on molecules in solution. J. Chem. Phys. 2002;117(1):43–54.
69. Ciccotti G, Ferrario M, Hynes JT, Kapral R. Constrained molecular dynamics and the mean potential for an ion pair in a polar solvent. Chem. Phys. 1989;129:241–251.
70. Kumar S, Bouzida D, Swendsen RH, Kollman PA, Rosenberg JM. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem. 1992;13(8):1011–1021.
71. Boczko EM, Brooks CL., III Constant-temperature free energy surfaces for physical and chemical processes. J. Phys. Chem. 1993;97(17):4509–4513.
72. Koslover EF, Wales DJ. Comparison of double-ended transition state search methods. J. Chem. Phys. 2007;127(13):134102. [PubMed]
73. Weinan E, Ren W, Vanden-Eijnden E. Finite Teperature String Method for the Study of Rare Events. J. Phys. Chem. B. 2005;109(14):6688–6693. [PubMed]
74. Maragliano L, Fischer A, Vanden-Eijnden E, Ciccotti G. String method in collective variables: Minimum free energy paths and isocommittor surfaces. J. Chem. Phys. 2006;125(2):024106. [PubMed]
75. Weinan E, Ren W, Vanden-Eijnden E. String method for the study of rare events. Phys. Rev. B. 2002;66:052301:1–052301:4.
76. Elber R. Long-timescale simulation methods. Curr. Opin. Struct. Biol. 2005;15:151–156. [PubMed]
77. Elber R, Cárdenas A, Ghosh A, Stern H. Bridging the gap between long time trajectories and reaction pathways. Adv. Chem. Phys. 2003;126:93–129.
78. Olender R, Elber R. Calculation of classical trajectories with very large time step: Formalism and numerical examples. J. Chem. Phys. 1996;105(20):9299–9315.
79. Khavrutskii IV, Byrd RH, Brooks CL., III A line integral reaction path approximation for large systems via nonlinear constrained optimization: Application to alanine dipeptide and beta-haripin of protein G. J. Chem. Phys. 2006;124(19):194903. [PubMed]
80. Quapp W, Heidrich D. Analysis of the concept of minimum energy path on the potential energy surface of chemically reacting systems. Theor. Chim. Acta. 1984;66(3–4):245–260.
81. Sana M, Reckinger G, Leroy G. An Internal Coordinate Invariant Reaction Pathway. Theor. Chem. Acc. 1981;58(2):145–153.
82. Basilevsky MV. Modern Development of the Reaction Coordiante Concept. J. Mol. Struct. (Theochem) 1983;103(12):139–152.
83. Lazaridis T, Tobias DJ, Brooks CL, III, Paulaitis ME. Reaction paths and free energy profiles for conformational transitions: an internal coordinate approach. J. Chem. Phys. 1991;95(10):7612–7625.
84. Hirsch M, Quapp W. Reaction pathways and convexity of the potential energy surface: applicaitons of Newton trajectories. J. Math. Chem. 2004;36(4):307–340.
85. Maragliano L, Vanden-Eijnden E. On-the-fly string method for minimum free energy paths calculation. Chem. Phys. Lett. 2007;446:182–190.
86. Peters B, Heyden A, Bell AT, Chakraborty A. A growing string method for determining transition states: Comparison to the nudged elastic band and string methods. J. Chem. Phys. 2004;120(17):7877–7886. [PubMed]
87. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in Fortran 77: The Art of Scientific Computing. Vol. 1. Port Chester, NY: Cambridge University Press; 2001.
88. Cho AE, Doll JD, Freeman DL. The construction of double-ended classical trajectories. Chem. Phys. Lett. 1994;229(3):218–224.
89. Eckart C. Some studies concerning rotating axes and polyatomic molecules. Phys. Rev. 1935;47:552–558.
90. Kudin KN, Dymarsky AY. Eckart axis conditions and the minimization of the root-mean-square deviation: Two closely related problems. J. Chem. Phys. 2005;122(22):224105:1–224105:2. [PubMed]
91. Kabsch WA. solution for the best rotation to relate two sets of vectors. Acta Crystallogr. 1976;A32:922–923.
92. Trzesniak D, Kunz A-PE, van Gunsteren WF. A Comparison of Methods to Compute the Potential of Mean Force. ChemPhysChem. 2007;8:162–169. [PubMed]
93. Yu H-A, Roux B, Karplus M. Solvation thermodynamics: An approach from analytic temperature derivatives. J. Chem. Phys. 1990;92(8):5020–5033.
94. Mark P, Nilsson L. Structure and Dynamics of Liquid Water with Different Long-Range Interaction Truncation and Temperature Control Methods in Molecular Dynamics Simulations. J. Comput. Chem. 2002;23(13):1211–1219. [PubMed]
95. Mark P, Nilsson L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A. 2001;105(43):9954–9960.
96. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983;79(2):926–935.
97. Lamoureux G, Jr, Roux B. A simple polarizable model of water based on classical Drude oscillators. J. Chem. Phys. 2003;119(10):5185–5197.
98. Ishimura K, Pulay P, Nagase S. A new parallel algorithm of MP2 energy calculations. J. Comput. Chem. 2006;27(4):407–413. [PubMed]
99. Johnson ER, Becke AD. A unified density-functional treatment of dynamical, nondynamical, and dispersion correlations. II. Thermochemical and kinetic benchmarks. J. Chem. Phys. 2008;128(12):124105/1–124105/3. [PubMed]
100. Becke AD, Johnson ER. A unified density-functional treatment of dynamical, nondynamical, and dispersion correlations. J. Chem. Phys. 2007;127(12):124108/1–124108/8. [PubMed]
101. Becke AD, Johnson ER. A density-functional model of the dispersion interaction. J. Chem. Phys. 2005;123(15):154101/1–154101/9. [PubMed]
102. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Montgomery JA, Jr, Vreven T, Kudin KN, Burant JC, Millam JM, Iyengar SS, Tomasi J, Barone V, Mennucci B, Cossi M, Scalmani G, Rega N, Petersson GA, Nakatsuji H, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Klene M, Li X, Knox JE, Hratchian HP, Cross JB, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Ayala PY, Morokuma K, Voth GA, Salvador P, Dannenberg JJ, Zakrzewski VG, Dapprich S, Daniels AD, Strain MC, Farkas O, Malick DK, Rabuck AD, Raghavachari K, Foresman JB, Ortiz JV, Cui Q, Baboul AG, Clifford S, Cioslowski J, Stefanov BB, Liu G, Liashenko A, Piskorz P, Komaromi I, Martin RL, Fox DJ, Keith T, Al-Laham MA, Peng CY, Nanayakkara A, Challacombe M, Gill PMW, Johnson B, Chen W, Wong MW, Gonzalez C, Pople JA. Gaussian 03, ReVision B.1. Pittsburgh, PA: Gaussian, Inc; 2003.
103. Tsuzuki S, Honda K, Uchimaru T, Mikami M, Tanabe K. The magnitude of the CH/pi interaction between benzene and some model hydrocarbons. J. Am. Chem. Soc. 2000;122(15):3746–3753.
104. Ulitsky A, Elber R. A new technique to calculate steepest descent paths in flexible polyatomic systems. J. Chem. Phys. 1990;92(2):1510. CT800086S.