|Home | About | Journals | Submit | Contact Us | Français|
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.
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 experimental1–4 and computational5–18 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.3–5,7,8,11–18 From the computational perspective, this task demands computing multidimensional potentials of mean force (PMFs) for which efficient tools have been lacking.19–36
Recently, we have developed and generalized the gradient-augmented Harmonic Fourier Beads (ggaHFB) method21–23 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.42–44 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.45–49 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.1–3,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.52–56 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 B3LYP61–63 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)65–68 with the MM and QM energy functions, respectively.
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.
The generalized gradient-augmented Harmonic Fourier Beads (ggaHFB) method considers an arbitrary system of N atoms described by 3N generalized coordinates = (q1, … ,q3N), and, equivalently, by 3N Cartesian coordinates = (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 = (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 = (qS+1, … ,q3N) comprise the spectator coordinate space (SCS) and do not contribute explicitly to the energy gradient.
Here superscript b indicates the bias, and 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.
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 () = U(,) 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.
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.
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 b,ref in RCS, one can compute the corresponding Cartesian free energy gradient via eq 3a.22,23,39,40
Here Wu is the unbiased free energy, kB is the Boltzmann constant, T is the simulation temperature, and |J()| 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 b = (q1b, … ,qSb) replaces the local minimum 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” b,j = (q1b,j, … ,qSb,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:
Importantly, the averaged configuration provides the complete free energy gradient in RCS and not just one of its components:
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.
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.19–26,72–79 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,80–85 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
where () is the curve tangent and 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:
otherwise different nonlinear coordinate systems will yield different reaction paths for the same stationary points.37,38,80–83,85 Thus, to be invariant the transition path curve in nonlinear coordinates must satisfy the following more complicated condition
where () 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.
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:
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.
To optimize a transition path in Cartesian coordinates, we take K unique configurations that gradually progress from the reactant to the product and assign them to a uniform grid with mesh size of Δα = 1/(K−1). If initial configurations k = (αk) 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
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 , the evolution returns either the minimized or the averaged bead 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.73–75,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.21–23,37–40,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”:
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,89–91 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.
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:
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.
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.
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 model93–97 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.
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 α-strand52–54 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.
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 α-strand52–54 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.
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.
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.98–101
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.
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.
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)65–68 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.
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.
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.98–101 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.
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.
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.