Selectivity of a Singly Permeating Ion in Nonselective NaK Channel

Sep 16, 2015 - Among all of the possible binding sites, a vestibule is noticed to be nonselective and seen to act as a probable binding site only in t...
0 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCB

Selectivity of a Singly Permeating Ion in Nonselective NaK Channel: Combined QM and MD Based Investigations Biswajit Sadhu,† Mahesh Sundararajan,‡ and Tusar Bandyopadhyay*,‡ †

Radiation Safety Systems Division and ‡Theoretical Chemistry Section, Bhabha Atomic Research Centre, Mumbai 400 085, India S Supporting Information *

ABSTRACT: Ion channels, such as potassium channels are known to discriminate ions to achieve remarkable selective transportation of K+ over Na+ through the membrane. The recently reported NaK ion channel, on the contrary, seems to be an exception, as it is observed to permeate most of the group IA alkali metal cations and hence is suggested to be nonselective in nature. However, does that correspond to a complete annihilation of selectivity inside the selectivity filter (SF) of the channel? What is the origin of such nonselectivity/selectivity, if any? The present computational study is an extensive multiscale modeling approach to find the probable answers to these intriguing questions. Here, we have used density functional theory (DFT) based calculations using a realistic truncated model of SF from the crystal structures of the NaK ion channel to evaluate the binding of various alkali metal ions (Na+, K+ and Cs+), free from “contamination” due to the absence any other “rivalry” cations, in its different binding sites. Among all of the possible binding sites, a vestibule is noticed to be nonselective and seen to act as a probable binding site only in the presence of multiple ions. Binding sites S3 and S4 are found to be selective for K+ and Na+, respectively. As an important observation, we find that calculations on oversimplified models using an isolated ion binding site may lead to an erroneous selectivity trend as it neglects the synergetics of consecutive binding sites on the final outcome. Energy decomposition analysis revealed ion-dipole electrostatics as the major contributing interaction in metal-bound binding sites. Our investigations find that although NaK is permeable to monovalent alkali metal ions, strongly “site specific” selectivity does exist at the three well-defined noncontiguous binding sites of the SF. Different important physicomechanical parameters (such as ligating environment, synergistic influence of binding sites, and topological constraints) are found to be the determining factor to induce the “site specific” selectivity of ions during translocation. Wherever possible, our computed results are compared with the available experimental findings. We finally conduct a detailed umbrella sampling-corrected metadynamics simulation in order to obtain an ion permeation free energy landscape within the SF that corroborates well with the “site specific” selectivity trend. chemically equivalent K+ binding sites, selectively permeating K+ over Na+.9−11 Recently, Shi et al.12 determined the crystal structures of a new class of ion channel, namely, the NaK ion channel from Bacillus cereus. The structure of the NaK channel comprises four subunits arranged in 4-fold symmetric fashion to create a pore at the center. Each subunit is composed of outer and inner transmembrane helices, namely, M1 and M2, which are connected to each other by a re-entrant loop or pore helice. SF with a structural sequence 63TVGDG67 (Figure 1) is a construct of the loop. Note that the functional groups (carbonyl, carboxylate, and hydroxyl groups) of the five residues can possibly create five binding sites. However, out of these amino acids, the carboxylate group of D66 points

1. INTRODUCTION Ion permeating channels are pore forming transmembrane proteins,1 which play critical roles in the transmission of nerve impulse and control several biochemical processes such as regulating cardiac,1 skeletal,2 and muscle contraction.3 The salient feature of various ion channels is their striking ion selectivity. For example, the potassium channel (K+ channel) has a distinct quality of selecting K+ in the presence of Na+ at a ratio of 1000:1,1 in spite of their same valence and miniscule difference in ionic radii (0.36 Å). In recent years, several crystal structures of various ion channels have been reported.4−8 All ion channels have a specific structural sequence in their selectivity filter (SF), which discriminately mediates the translocation of ions. For instance, K+ channels are known to preserve a structural signature sequence of 75TVGYG79 in their SF. The backbone carbonyl ligands of these five residues along with the hydroxyl group from T75 create four contiguous and © 2015 American Chemical Society

Received: June 23, 2015 Revised: September 9, 2015 Published: September 16, 2015 12783

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B

Figure 1. NaK ion channel and its different binding sites in the selectivity filter (SF) are shown. Color schemes used are yellow, C; blue, N; red, O; and gray, H. The left panel shows the ribbon presentation of the tetrameric channel helix. The signature sequence, 63TVGDG67, forming the SF is presented in licorice representation. During the all-atom molecular dynamics (MD) simulation, the entire tetraameric channel is used. The right panel shows the zoomed in view of the SF and its various binding sites. Realistic models of the SF binding sites, structure A, and structure B that has been utilized to perform DFT calculations are shown. Structure B is further truncated (in the truncated model, TM) into stand-alone binding sites, vestibule, S3, and S4, in order to study the ion binding characteristic in isolation. Selected Cα atoms, shown with an asterisk in the model structures, A and B, are physically constrained in order to preserve the optimum pore radius (see methods below). For clarity reasons, hydrogen atoms are not shown in the case of structure A, B, and truncated models of structure B.

recently to explain the “nonselectivity” occurring in the NaK channel. Available experimental structures on wild type and engineered NaK channels suggest that site S3 is most nonselective13 and that three binding sites in place of four (as they are in selective K+ channel)11 cause the loss of its discriminating character. The latter observation is further substantiated by a potential of mean force (PMF) calculation of the ion translocation through the SF of the engineered NaK channel.20 Computational studies on the nonselectivity of the NaK channel entreat a variety of concepts: hydration number (HN), variable metal-protein ligand coordination number (CN), ratio of HN and CN of the permeating cations, water containing vestibule that disrupts the protein SF matrix causing planes of carbonyl oxygen to dislocate, flexibility of the dynamical hydration valve in the vestibular region, and the presence of a hydrogen bond (HB) between D66 and N68.15−19,23,26 Although nonselectivity is the hallmark of the NaK channel, that does not render the same permeability for all alkali metal ions. Structural investigations led by Alam et al.13 revealed that Na+, K+, and Rb+ ions have their affinity toward specific binding sites of SF. In line with this, several computational studies also observed significantly different free energy barriers favoring the higher permeability of K+ than Na+.16,18,19 The above background of our current understanding on monovalent ion permeation raises one intriguing conundrum: the “site-specific” selectivity in the “nonselective” SF of NaK channel. Herein, we address this issue by elucidating the structure, binding and translocation behavior of ions

outside the pore axis, and the carbonyl group of G65 and D66 remains tangential to the pore axis.13 Hence, G65 and D66 practically do not bind to the metal center during ion translocation, which eventually makes the total number of binding sites equal to four including vestibule. The amino acid sequence in the SF of NaK shows strong resemblance to that of the K+ channel (TVGYG) except that tyrosine (Y) is now replaced by the acidic aspartate (D) residue. This led to significant conformational changes in NaK SF, where a wide water vestibule is formed in place of two distinguishable binding sites (cf. sites 1 and 2 of the K+ channel). Reported structural investigations12−14 revealed that the external entrance formed by the carbonyl groups of G67 can act as the probable binding site (Sext) with the involvement of nearby water molecules. The other two sites, S3 and S4 are similar to those in the K+ channel (see Figure 1). Noticeably, in contrast to the K+ channel and similar to cyclic nucleotide gated ion channels (CNG channels), NaK is nonselective in nature and thus can permeate most of the alkali metal cations (Na+, K+, and Cs+) along with few divalent metal cations (e.g., Ca2+).14 Ever since the nonselective nature of NaK is known, a plethora of activities began,11,13,15−23 all aimed at deciphering the root reason behind nonselectivity because that could serve as a testbed for the different mechanisms proposed toward explaining selectivity in ion discriminating channels. While lingering over the years, our understanding of the selectivity issues of tetrameric channels has advanced and became more mature,11,24,25 similar such debetable issues began to crop up 12784

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B

binding site, vestibule, formed by the carbonyl group of V64, is wide enough to accommodate three to four water molecules, where the permeating cation can easily get diffused away. As mentioned earlier, D66 does not play any role in cation binding at the vestibule region. (3) Unlike Sext, the carbonyl groups of V64 and T63 residues together form S3. Although, the structural consideration indicates that binding at the vestibule and S4 requires involvement of only one residue (V64 for vestibule; T63 for S4), the cation binding at both the sites indirectly gets affected by other connected consecutive binding sites, either because of conformational rearrangement or because of the presence of other cations. To this end, in order to study the synergistic effects of contiguous binding sites, structure B was further truncated to form small synthetic models (denoted as TM), each containing only one binding site, i.e., vestibule, S3 and S4. The isolated binding sites of the TM model (Figure 1) are thought to be free from the effects of the presence of other nearby binding sites. 2.1.2. Optimization Protocol. All starting structures were optimized using a pure generalized gradient approximated BP86 functional29,30 with medium sized def2-SV(P) basis set,31 which was earlier shown to produce good structural parameters in accord with the experiments.32−34 Nevertheless, a set of calculations are performed on TM structures using a larger basis set (def2-TZVP)35 to comment on the basis set dependency over the geometry and energetics of metal-bound systems. Neither the alternation of the binding affinity trend nor any significant structural changes were noticed with the use of larger basis sets (see Supporting Information, Table S1). To speed up the calculation, a resolution of identity (RI) approximation is invoked for coulomb exchange integrals using a corresponding auxiliary basis set (def2-SV(P)/J). For energy evaluation, a three-body dispersion corrected36 hybrid B3LYP37,38 functional (B3LYP-D3) with a triple ζ TZVP basis set35,39 was used. For heavy metal Cs+, the core orbitals (containing 46 core electrons) were modeled via the def-ECP pseudopotential. def2-SV(P) (for geometry optimizations) and def-TZVP basis sets (for energy minimizations) were used to describe the valence orbitals. All geometry optimizations were carried out at gas phase with tight SCF convergence criteria (10−8 a.u. on the density). The basis set superposition error (BSSE) for the metal-bound complexes were found to be within the window of 1−3 kcal mol−1, indicating that counter poise corrections will have negligible influence on subsequent energetics. In the process of energy evaluation, a COSMO continuum solvation model40 was used to take care of longrange electrostatic interactions as implemented in TURBOMOLE v6.3.41 In the NaK ion channel, the secondary protein backbone attached to the SF region works as a structural constraint to maintain the vicinity of four nearby subunits. This in turn preserves the optimum pore radius to form appropriate cation binding sites. Optimization of our model structures without such a constraint led to a significant increase of pore radius (data not shown). In most occasions, a maximum three out of four subunits are observed to interact with the attached cations, thus underestimating the cation binding affinities. Previous MD studies42,43 indicated the importance of “topological” constraints on the coordination number in proteins. Structural analysis of ion selectivity in NaK further emphasizes the rigidity of the SF during ion binding.13 The presence of a methylene spacer in the synthetic models used by Dudev and Lim24,25 served as the needed constraint during optimization. Keeping

through the NaK channel, covering the utmost interesting physiological and thermodynamic aspects that influence the site specifity, through a combination of quantum mechanical and classical mechanics simulations. Along with Na+ and K+, we also consider the radiologically important Cs+ ion, as there are reports which suggest that several ion channels can be used as a transport route for Cs+.27 In this work, we have used a density functional theoretical (DFT) approach on a realistic model of NaK SF to determine the origin of site-specific ion selectivity, while metadynamics (MtD) based molecular dynamics (MD) simulations were performed to construct the PMF profiles of the translocating ions in the wild NaK channel.12 The use of synthetic toy models resembling the binding sites of the ion channel often lead to valuable insights.24,25 However, as we show here in the case of the NaK channel, use of a further truncated model (TM) (considering only one binding site at a time), derived from our realistic model, can lead to erroneous conclusions and fail to capture the synergetics of contiguous binding sites (such as treating S3 and S4 together), in accordance with experimental observation.11 Finally, we focus on arriving at a clearer picture of the process of site-specific binding and translocation using the combined results of DFT and MtD-MD based calculations.

2. COMPUTATIONAL DETAILS 2.1. Quantum Mechanical Study. 2.1.1. Model Structures. The published crystal structure of WT NaK contains two subunits,12 which were duplicated to create a functional tetrameric channel and is available from the potassium channel database (KDB).28 Before the QM calculations were initiated, we have performed a 10 ns long equilibrium MD simulation (see below). The last frame of this MD trajectory can be thought to represent an equilibrated structure of the NaK channel. The starting structures for QM calculations were built from an equilibrated structure of the NaK ion channel by extracting the five SF residues (namely, T63, V64, G65, D66, and G67) from each of the four subunits. The resulting structure with all four probable binding sites representing the full SF consists of more than 270 atoms. Performing QM calculations on such a gigantic systems is highly expensive. Thus, to optimize the computational cost, the structure was truncated in between D66 and G65 (Figure 1). This leads to two model structures (A and B) with different binding sites. Structure A is composed of G67 and D66 residues consisting of only one binding site for the cation, i.e., the external entrance or Sext, whereas structure B contains the remaining sites, namely, vestibule, S3, and S4. Further, structures A and B were modified at the peptide junction by substituting the CO-NH- group with −CONH(CH3) in order to avoid unreasonable hydrogen bond formation among the nearby subunits that could have caused the filter’s structural distortion. The truncated structures are not ad hoc and can be justified by the following reasons. (1) Previous structural investigations12−14 predict that the binding of the cation at the external entrance (Sext) is mediated only by the carbonyl group of the G67 residue. Notably, the presence of the hydrogen bond between the negatively charged carboxylate group of D66 and the amide hydrogen of the G67-N68 backbone helps to make the entrance rigid, which in turn maintains structural integrity. As our model structure for Sext contains −CONH(CH3) (representing N68 side chain) at the terminal of the G67 residue, it preserves the “through space” interaction,13 naturally occurring inside the NaK SF. (2) The size of the next possible 12785

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B these facts in mind, physical constraints were applied to a few selected backbone carbon atoms of model structures (for structure A, one and for structure B, two carbon atoms). In Figure 1, these three Cα atoms are marked by *. The coordinates of these carbon atoms are “frozen” in order to restrict their movements in the Cartesian space. The selection of backbone carbon atoms were made in such a way that the applied constraint should not restrict the free movement of directly coordinating protein ligands (carbonyl and hydroxyl groups). Analytical harmonic frequency calculations were performed for all structures to characterize the stationary points. Thermodynamic quantities such as changes in Gibbs free energies (ΔG) and entropy (ΔS) were derived. Natural population analysis (NPA)44 was performed to quantify the charge transfer among metal, SF, and water molecules. We would like to point out that most of the reported optimized structures are at their potential minima, without any imaginary frequency in their respective potential energy surface. Note that only in a few cases the imaginary frequencies of values ranging from 5i−35i were obtained, which correspond to the presence of constraints in the structure. Throughout, only real modes are considered in order to compute the contribution of various components in the partition function calculations. We believe that the error associated with these small number of cases will have a negligible effect45 over the major conclusions drawn on overall binding energetics. 2.1.3. Reaction Scheme. Capture and subsequent permeation of an ion from the bulk to NaK SF is an outcome of the chelating strengths of the water molecules in the bulk and the protein ligands in the SF. First, using the following reaction schemes we have calculated the binding free energies (BFE) at different ion-chelating sites of NaK for each of the three cations (M+, M = Na, K, and Cs), not in competition with other ions.

S4·M1+ + M2(H 2O)6+ → S4·M2+ + M1(H 2O)6+

(5)

ϵ ϵ ΔGexϵ = ΔGexgas + [ΔGsolv ]product − [ΔGsolv ]reactant

(6)

where ΔGgas ex is the free energy change associated with ion exchange reactions in the gas phase, and ΔGϵsolv is the free energy change in transferring a reactant or product partner from gas phase to solution phase with dielectric constant, ϵ. The difference in the sum of electronic energies, thermal energies including zero point energy, and the work term with the entropy change between reactants and products at T = 300 K is used to find ΔGgas. The dielectric constant of the ion channel is known to vary from 5 (for deeply buried binding sites 3 and 4) to 30 (at the external entrance),43,46,47 and thus, the behavior of ion binding is expected to be different from that of the bulk water environment (ϵ = 80). To model the different dielectric media mimicking the channel interior, calculations were performed with three different ϵ values (5, 10, and 20), using gas phase optimized structures. Additionally, another set of calculations were performed at ϵ = 80 to study and compare the fate of these reactions in the bulk. 2.1.4. Energy Decomposition Analysis. To determine the nature of interaction between the alkali metals and SF, energy decomposition analysis (EDA) was performed using the methods of Morokuma48 and Ziegler49 at the GGA/B3LYPD3/TZP level with ZORA,50 as implemented in the program package ADF 2013.51 The total bonding energy (ΔEbond) between the two fragments (metal and protein) was calculated using the following equation:

(1)

ΔE bond = ΔEelstat + ΔEpauli + ΔEorb + ΔEdisp

(2)

(7)

where, ΔEelstat and ΔEpauli represent the electrostatic and repulsive interaction energy contributions between the fragments, respectively. ΔEorb is the stabilizing energy that arises from the orbital contributions and thus represents the strength of covalent bonding between the fragments. ΔEdisp accounts for the dispersion interaction in the process of binding. 2.2. Molecular Dynamics Study. 2.2.1. Structure Modeling. Each subunit of tetrameric NaK protein channel contains 1−104 residues. We have simulated the full protein channel (PDB ID: 2AHY) including the M0 helix (see in contrast, Vora et al.16). The lipid bilayer from the initially equilibrated lipid bilayer-embeded solvated tetrameric channel28 was removed to save computation time. In this study, the presence of the M0 helix was used to maintain the structural integrity of the tetrameric ion channel. However, it is important to mention that although in a real translocation scenario, the M0 helix mediates the opening and closing of the ion channel, the presence of the M0 helix does not cause any structural changes at the SF.14 Four Cα atoms from each of the subunits of M0 helix were tethered using positional constraints with a force constant of 2.4 kcal mol−1 Å−2. The chosen Cα atoms are at the junction where two M0 helices interact with each other due to the proximity of two nearby subunits. The SF region of the channel was thus allowed to reorganize in response to the ion

Crystallographic difference map analysis of the metal-bound NaK complexes13 indicates four water molecules are attached in the first-shell of metal ions during their capture into the Sext site from the bulk. In view of this, we performed calculations with varying HNs (n = 0, 2, 4) at Sext. As the binding sites S3 and S4 at the SF interior are analogues to the corresponding binding sites of the K+ channel, the available carbonyl groups at the respective site can satisfy the overall CN of metal cations. Hence, no water involvement in binding is considered at these sites. Hydrated metal structures used in the above equations are obtained using successive MD and DFT calculations on metal water systems (see below). Site specific selectivity of a given ion in the presence of other copermeant ions can be judged by the free energies of ion exchange reactions. Ion exchange (M1 + → M2 +; M2 substituting M1) free energies between two metal cations, M1+ and M2+ (viz., ΔGK+→ Na+, ΔGCs+→ K+, and ΔGCs+→Na+) at the binding sites are calculated using the following exchange reaction schemes: Sext ·M1+ · (H 2O)n = 0,2,4 + M2(H 2O)6+ → Sext ·M2+·(H 2O)n = 0,2,4 + M1(H 2O)6+

(4)

In the above equations, M+-(H2O)6 represents the hydrated metal ion. The SF bound ions are represented as S3-M+ and S4M+ for sites S3 and S4, respectively, and Sext-M+-(H2O)n=0,2,4 represents the partially hydrated (n = 0,2,4) ion at the external site, Sext. The free energy change associated with the above equations is given by

Sext + M(H 2O)6+ → Sext ·M+(H 2O)n = 0,2,4 + (H 2O)m = 6,4,2

S3/S4 + M(H 2O)6+ → S3/S4·M+ + (H 2O)6

S3·M1+ + M2(H 2O)6+ → S3·M2+ + M1(H 2O)6+

(3) 12786

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B

elegant method for conformational search along a chosen set of CV, and, unlike the umbrella sampling (US) method, is free from the prior assumptions on final state. Here, the MtD method is used to generate the ion permeation pathway with absolute certainty. The MtD discovered PMF is then refined with the help of US correction. The details of this corrective step can be found in refs 62−64. In brief, in the MtD-US correction approach, first a single MtD run along the chosen CV is used to discover a preoptimized biasing potential, VMtD(s(x),tb) up to a certain build-up time, tb. Subsequently, the MtD run is stopped, and a conventional equilibrium US simulation is initiated along the MtD-discovered CV path using the physical force field derived potential energy, Vphys, additionally incremented by the time independent biasing potential (non-negative) that was accumulated over the previous MtD run. The US step allows ergodicity to be satisfied since the metastabilties along the permeation pathway are removed by the elevated biasing potential (obtained through MtD step) and the diffusive dynamics in the chosen CV space can be sampled with ease with a reasonable computational effort in the spirit of a equilibrium run. After the US run is completed, we have calculated the biased probability density, ρUS(s), of finding the system at given COM spearation, s, between the permeating ion and the channel through histogram sorting. Finally, the US corrected PMF is given by

permeation, while the tethered atoms, away from the SF region, fluctuate in a way similar to that in a lipid-embeded system.52 2.2.2. Simulation Strategy. All MD simulations are performed at 300 K in the NVT ensemble with an amder99sb-ildn force field53 using the GROMACS 4.5.4 suite.54 The TIP3P water55 solvated ion channel with the added ions (to manitain electroneutrality) comprised of more than 0.26 million atoms. For short-range electrostatic and van der Waals interactions, a cut off of 14 Å was applied. For the long-range electrostatics interaction, the particle mesh Ewald56 method was used, while the LINCS algorithm57 was implemented to constrain the bond between hydrogen and the heavy atoms at their equilibrium lengths. Finally, a 10 ns equilibration run was conducted by integrating the equation of motion with a 2 fs time step. RMSD of the structure was calculated with respect to the initial structure considering the backbone Cα atoms. Insignificant deformations were noted upon equilibration ( 5 at Sext is unfavorable. This supports the notion that small Na+ can remain favorably bound in a tetra-coordinated protein ligand environment. Subsequent addition of a water molecule led to strongly favorable binding for K+ as compared to Cs+. With the involvement of two or four water molecules, the binding of Na+ and K+ is noted to be very similar (difference of BFE is 1−4 kcal mol−1) (see Table S3). The effect of methyl substitution at Sext increases the cumulative charges over carbonyl due to the lack of charge delocalization and also somewhat enhances the structural flexibility at the entrance, which eventually led to high binding energy at ϵ < 80. However, BFE was found to decrease substantially at ϵ = 80. Notably, obtained BFE does not seem to 12790

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B Table 4. Energy Decomposition Analysisa sites

ion

ΔEpauli

ΔEelstat

ΔEorb

ΔEdisp

total

Sext

Na+ K+ Cs+ Na+ K+ Cs+ Na+ K+ Cs+

28.67 46.94 59.66 13.52 27.82 37.94 25.82 42.9 33.59

−250.87 (83.42) −244.73 (84.30) −239.1 (81.29) −98.59 (66.87) −97.69 (68.71) −88.88 (65.90) −94.66 (67.13) −91.52 (65.02) −69.88 (63.78)

−39.85 (13.25) −33.88 (11.67) −37 (12.58) −37.66 (25.54) −32.97 (23.19) −32.73 (24.27) −38.23 (27.11) −35.93 (25.53) −29.58 (27.00)

−10.02 (3.33) −11.69 (4.03) −18.05 (6.14) −11.18 (7.58) −11.52 (8.10) −13.27 (9.84) −8.11 (5.75) −13.3 (9.45) −10.11 (9.23)

−272.07 −243.36 −234.49 −133.91 −114.36 −96.94 −115.18 −97.85 −75.98

S3

S4

a

Energies are in kcal/mol unit. The values in parentheses refer to the percentage contribution of each energy term in total energy (i.e., sum of ΔEelstat, ΔEorb, and ΔEdisp).

Figure 3. (a) Pore internal surface profile of the NaK ion channel. The dots show the locus of the outer sphere of a flexible sphere with variable radius, squeezing through the pore, which finally sketches the central line of the pore (shown in yellow). (b) Plot of distance along the channel axial position (Z) vs pore radius as obtained from the average structure of the ion-SF trajectories for the three ions using the HOLE suite of programs.74 Axial position is measured from the average center of the SF.

correlate well with the increase of M+-OG67 bond distance. It is to be remembered that our calculated BFE at Sext does not truly correspond to the protein binding but reflects the resultant binding energy of metals in a mixed ligand environment (protein and water molecules). As explained earlier, in the process of single ion translocation through the SF, ions initially at the vestibule are observed to drift toward the more preferable binding site at S3. Hence, binding energy calculations are not performed for the vestibule. For all of the dielectric media (ϵ = 5 to 80), binding strength at S3 is observed to follow the order: K+ > Cs+ > Na+. It is important to note that S3 is a deeply buried binding pocket. Hence, expected ϵ is well below the values of the aqueous medium. Calculated BFE values suggest that the decrease of ϵ substantially enhances affinity toward Cs+ as compared to Na+ and K+. Contrary to S3, the binding strength of K+ at S4 is found to be the least as compared to Na+ and Cs+. Binding of Na+ and Cs+ occurs close to the plane of carbonyls formed by T63, which might cause higher binding strength (more negative ΔG). Calculated thermodynamic parameters strongly suggest that the binding is enthalpy driven and that entropy does not have significant influence on free energies (see Table 2). Nevertheless, as compared to Na+, positive entropic contribution for Cs+ seems to favor its binding at all of the binding sites,

which can be accredited to its low hydration free energy (−59.8 kcal mol−1). 3.2.3. Ion Exchange Reaction at the Selectivity Filter. To evaluate the site-specific preferential binding of alkali ions, free energy change associated with the ion exchange reactions (cf. eqs 3 to 5) were calculated, and the results are presented in Table 3. Gas phase energies are calculated here only for reference purposes. On the basis of this, the implicit solvation model results are calculated at different dielectric media, mimicking the protein matrix in solution. Comparison of ΔGϵex (ϵ = 5, 10, 20, 80) values for the triad of exchange reaction at the external site shows that Na+ has slighly more selectivity as compared to K+, while Cs+ at Sext is noted to have the least selectivity. Interestingly, at low dielectric media (ϵ < 20), in the methyl substituted Sext site, the preferential binding order changes with Cs+ becoming most preferred. This is reminiscent of the fact that in the absence of this HB, the filter loses its compact structural feature, allowing a larger Cs+ to be accommodated more easily. Further, a marginal loss of selectivity is also noted for the Na+-K+ ion exchange couple. Akin to the BFE trend, the selectivity order of ions at S3 is found to be K+ > Cs+ > Na+. Earlier experimental studies on selectivity13 also observed a similar trend, where Cs+ is found to have higher affinity at S3, which however cannot replace K+ at 12791

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B higher concentration. Contrary to S3, K+ has the least selectivity for S4 as compared to Na+ and Cs+. 3.2.4. Energy Decomposition Analysis (EDA). Energy decomposition analysis was performed to distinguish the contribution of various interactions to total bonding energy (see Table 4). In general, favorable ion−dipole interactions (estimated contribution of 60−85% to the total energy) are held primarily responsible for binding at the various sites of the ion channel. Because of the increase of ionic radius from Na+ to Cs+, Pauli repulsion contribution (ΔEpauli) to the total bonding energy is increased, whereas lowering of charge density led to the decrease of favorable electrostatic interaction (ΔEelstat) as we move down the periodic table. We also note a structural change associated with the S4 site after binding to Cs+. The hydroxyl side chain was observed to move away in order to accommodate the large cation. Here, Cs+ interacts with the carbonyl group of T63 and lies below the carbonyl plane, resulting in the lowest ΔEelstat for Cs+. Dispersion also contributes to the total bonding energy (3.33 to 15.36%), and the contribution increases as the ion moves toward the deeply buried binding pockets. Owing to the low charge by radius ratio of Cs+, dispersion contributed more to the binding energy as compared to Na+ and K+. Orbital interactions (ΔEorb) primarily come from the metal−carbonyl binding interactions via covalent interactions. Thus, increase in the number of carbonyl ligands at the binding periphery of the cation increases the orbital energy contributions, which is reflected in the higher value of ΔEorb at the S3 (∼23−26%) and S4 (∼25−27%) as compared to Sext (∼11−13%). 3.3. Metal Ion Permeation through the NaK Channel. One of the intriguing observation made by Alam and Jiang13 while analyzing the structural details of ion-occupied NaK channel is that its SF hardly rearranges in response to the ions. In view of this, we first investigate the filter’s structural integrity during single ion permeation through it, while the channel is embedded in an explicit solvent environment. After the completion of a MtD event, i.e., complete permeation of ions through the filter, we utlized the channel trajectories for this purpose. MtD generated atom trajectories are then subjected to Monte Carlo simulated annealing treatment74 in order to find the best possible route for a size-adjustable sphere to swim through the channel. Graphical presentation of the annealed data is shown in Figure 3. The average channel pore radius is obtained from a simulation trajectory. Upon careful inspection of the ion-permeation trajectories, we have sometimes observed carbonyl flipping at the SF (in the absence of cation). However, frequencies of such events are indeed very low. In the presence of ion, like in experiments,13 here we did not observe any major changes in SF. Clearly, no significant deviations in the pore radius (see Figure 3b) are observed with respect to the size of the permeating ions.13 This observation further adds to the gross nonselective feature of NaK, although, as presented above, the filter possesses the site-specific selectivity characteristics at the microscopic level. Further, we analyzed the MtD trajectories for the calculation of van der Waals and electrostatic interaction energies of the system during the ion translocation process. Favorable ion permeation in the channel depends on the interaction energy contributions from various subgroups (ion−SF, SF−water, and ion−water) involved in the process. We found that at the bulk, ion−water interaction persists. As the ion enters the SF, favorable ion−SF interactions tend to increase (more negative), while partial dehydration of ion

reduces the ion−water interaction (more positive). The average of interaction energy of all the components showed a linear time dependent variation, with little change between the permeating ions, indicating the nonselective nature of the ion channel. However, changes in the variation of the individual components (such as ion−SF) of these interaction energies drastically differ from one ion to the other (data not shown). This result is again a manifestation of site-specific ion−filter interactions, which can be further substantiated by an analysis of the coordination number of the permeating metal ions. For the calculation of coordination numbers, a distance of 3.2 Å was taken as the radius of the first coordination shell for all metal ions. A plot of the total coordination number (with water and protein ligands) of the ions (see Figure S2) reveals that all ions maintain a coordination number of 6−8 at the bulk. During translocation, a coordination number drop of up to four was noticed for all the three ions. In Figure 4, the results of the ion coordination numbers with the SF residues are presented in order to depict their

Figure 4. Change of coordination number of the ions (top, Na+; middle, K+; and bottom, Cs+) with the protein ligand residues of the selectivity filter during their permeation. The different colors repsenting the ligands are G67, blue; D66, green; G65, black; V64, violet; and T63, red. Axes on the left and right sides of plots indicate the contribution of different residues of the filter to the coordination number and COM separation of NaK channel and the permeating ions, respectively. Shaded horizontal strips in each of the plots indicate binding sites of the filter with respect to the COM separation, as indicated.

participation and involvement during ion translocation. In this figure, the MtD CV (COM separation between NaK channel and the ions) is also plotted alongside as a function of MtD time. This CV, along with the participation of various residues in the SF, further helps delineate the locations of the binding sites of the filter. Among the five residues, major contributions from V64 (pink lines) and T63 (red lines) are evident in Figure 4, which is quite expected due to the small pore size at sites S3 and S4. Importantly, due to the tangential alignment of G65 (black) and D66 (green) carbonyls with respect to the pore axis in the vestibular region, no significant contributions due to these residues to total coordination are observed inside the pore. Expectedly, G67 (blue) is seen to contribute at the entrance quite efficiently for all of the ions and especially for the bigger sized Cs+. Note that Cs+ during 6−9 ns of the MtD result in Figure 4 does not coordinate with the 12792

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B

solution. This also provided us a reference unbound state for subsequent estimates of PMF, which is equivalent to having the unbound ion and the channel in water medium. The PMF results in Figure 5 suggest that the permeation of the three ions through the SF is a multistep process with free energy barriers and basins along the pathways. It is highly encouraging to note the resemblence of the selectivity order in each of the binding sites as obtained by the MtD-US method with those of the DFT results pertaining to ion exchange reactions. Moreover, the MtD simulation snapshots of the ion− ligand coordination structures pertaining to the Sext and S3 sites are found to be similar to those obatined by DFT calculations (shown in Figure 2). Entry from bulk to the filter is manifested with a high energy barrier caused by the shading of water molecules of the hydrated ions. Contrary to an earlier structural investigation,13 our calculations involving single ion permeation show that Na+ binds preferentially at Sext over the other two monovalent ions. At the vestibule, the ions remain mainly diffused in the absence of a definite chelating ligand (no minimum is seen in the PMF profile). Note that the DFT results of the TM model representing the vestibule are at variance with this finding pointing toward the weakness of the model. At site S3, all of the three ions exhibit free energy minima, albeit shallow. This indicates that S3 is mostly nonselective, however, as can be inferred from Figure 5 that K+ preferentially binds at S3, corroborating well with the DFT results (Table 4). Na+ is seen to be selective in between the S3 and S4 regions. The selectivity of the ions at S4 remained unresolved (when compared with DFT calculations) in the MtD-US results probably because of the participation of water molecules residing in the immediate cavity of the filter that restores the hydration sphere of ions showing minima in the cavity region.

protein ligands. Rather, during this time this has been found to coordinate with the channel water (cf. Figure S2). The ion-channel COM separation data show that ions starting from the bulk repeatedly enter and exit the channel during the MtD runs. Alongside these repeated events, the chealting ligands also respond in a systematic fashion. For example, G67 coordinates with Na+ at ∼2.5, 14, and 16 ns during the ion’s repeated attempts to re-enter the filter. As mentioned earlier, the accuracy and convergence61 of the MtD method is an outstanding issue, and sufficient sampling that would guarantee a converged PMF profile is difficult to predict a priori. As a result, in the present work, instead of generating the PMF directly from MtD, we use the method to explore the ion permeation pathway. Note that even though the reaction coordinate is defined exactly as normal to the lipid bilayer, a straight line is usually not a proper reaction coordinate for complicated ion transport mechanisms through a confined space such as ion channels.52 This is evidenced from Figure 4 where the time-dependent COM separation between the translocating ions and the channel can be seen to follow a path far from being a straight line. Once the pathway is detected, we use it to discover the converged results by the US correction technique (see Computational Details for the method) using eq 9. The use of the MtD-US approach is an improvement over techniques of either of MtD or the US. MtD improves the possibility of exploring the low probability regions, leading to a more reliable description of the underlying free energy landscape, which is free from the ad hoc choice of a straight path as is often assumed by the users of the US method. However, US allows circumventing the convergence problem associated with the MtD method. Thus, the convergence problem in the MtD method can be avoided, and the resulting PMFs are shown in Figure 5. We have ensured that the external area of the pore is sufficiently filled with the MtD hills bias such that the ions at the maximum distance from Sext are well solvated, like the free ion in dilute

4. DISCUSSION 4.1. Selectivity in the Nonselective NaK Channel. 4.1.1. Does the NaK Channel Have Selectivity? Ever since the SF structure of K+ channel responsible for the selective mediation of ions has been known, the means by which selectivity is achieved raises several useful debates that have enriched the literature. The NaK channel having a similar overall structure as K+, however, has been largely regarded as nonselective, equally permeable to Na+ and K+. A few noteworthy contributions15−22 in the field of nonselectivity in NaK focus on the multiple ion occupancy of the channel, whereas our molecular level knowledge of single ion permeation through NaK, free from all competing effects of rival cations remains incomplete. Present combined investigations using quantum and classical mechanics are an effort toward this direction with an aim to provide greater insight into it. Conclusively, in line with several other studies, we find that NaK is permeable to monovalent metal cations. Nevertheless, strongly site specific selectivity does exist inside the pore (Table 3 and Figure 5). For instance, K+ is found to have higher binding strength at the site S3, whereas Na+ is seen to be selective at the boundary of S3 and S4. However, at least for single ion translocation, no ionic selectivity was observed at the vestibule site. As can be seen in Figure 5, the wide variations in the PMF profiles of the three individual ions, depth of the free energy wells and the height of the barriers, will certainly lead to variations in the conduction of the ions as one would infer from Kramer’s rate.23

Figure 5. PMF profiles of alkali metal ion (Na+, K+, and Cs+) permeation as a function of COM separation between the NaK channel and the permeating ions (nm) are shown. The profiles are obtained after equilibrium US corrections of MtD generated free energy lanscapes. The minimum of the unbinding free energy for each of the ions is used as a reference point to arrive at this plot. Shaded vertical strips indicate the binding sites of NaK. Approaching toward the origin of the plot corresponds to the translocation of ions toward the central cavity. 12793

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B But question remains as to what attributes selectivity to the specific binding sites in an otherwise nonselective channel. Like in K+ channel,42,75−79 the combination of many interdependent factors, such as coordinating environment, conformational constraints, and type of interacting ligands, induce the site specificity. Our computed results suggest that even without the coordinated water molecules, smaller Na+ can favorably occupy the plane of the coordination center formed by the four carbonyls at Sext, while coordinated water molecules are needed for relatively larger K+ and Cs+ ions to make the favorable coordination (Table S3). This finding correlates well with several other studies,43,47,80−82 where octahedron arrangement was found to be an important governing factor for favorable K+ coordination. Stating this, we further emphasize that only a margianl selectivity difference at the Sext between Na+ and K+ was observed (BFE differs only by ∼4 kcal mol−1; see Table 2). This is not surprising as the solvent exposed face of Sext provides enough opportunity for both the ions to procure the suitable coordinating environment. This also stands true for the other two sites (S3 and S4), where K+ is observed to maintain the octahedron coordination arrangement, while Na+ is too small to interact with all of the carbonyls and thus prefer to stay at the junction of two sites with a total coordination number of four. Preservation of such necessary coordination numbers, particularly for Na+ (total coordination number drops to four during translocation) and K+ (maintains overall coordination number of eight), is largely maintained during the MtD analysis (Figure S2). The higher selectivity of K+ at S3 as compared to S4 can be understood by the close inspection of optimized structures (Table 1). At S3, The comparable bond distances of K+−OV64(CO) (2.899−2.932 Å) and K+−OT63(CO) (2.807− 2.816 Å) imply that K+ can occupy the center formed by eight carbonyls, which is certainly not the case at site 4 where the K+OT63(OH) bond length (3.045−3.195 Å) is noted to be significantly longer as compared to that of K+−OT63(CO) (2.577−2.628 Å). NPA analysis at these sites suggests that net charge transfer values of metals are in correlation with the charge density trend of alkali metals (see Table S4). However, the computed trend does not harmonize well with the obtained BFE values. This again emphasizes the importance of proper ligand orientations and other physicomechanical factors for the site-specific selectivity at the binding pockets. In summary, like in the K+ channel,4,6 the binding sites of the NaK channel are also ion-specific. In the K+ channel, however, all of the four binding sites selectively bind to only K+ ion, whereas here we have found that the binding sites of NaK possess varying selectivity toward Na+, K+, and Cs+ ions. For example, Na+ favorably binds at Sext (by ∼2 kcal/mol) when compared with K+. Similarly, sites S3 and S4 preferably bind K+ and Na+ ions, respectively. All of the binding sites are not selective toward Cs+ (in comaprison to other two ions), and the vestibule is the most nonspecific toward any of the three binding ions. This renders overall nonspecificity in an otherwise ion-specific NaK channel. 4.1.2. Importance of Physico-Mechanical Constraints. Our present investigations also shed light on the requirement of constraint for site selectivity. In the NaK channel, the presence of constraints due to the hydrogen bonding interactions between the carboxylate of D66 and G67-N68 peptide linkage restricts the movement of free moving G67. Important observations were reported in previous experimental studies83 on the mutation of D66 with alanine residue, where such changes led to a significantly higher flow rate for Na+ and K+

inside the pore. Our QM study predicts that substitution of the carboxylate group of D66 at Sext by the neutral methyl led to elongation of metal−ligand bond distance and shortening of the metal−Owater bond distance (Table 1). For K+ → Na+, the calculated ΔGϵex value is seen to reduce by 1−2 kcal mol−1. Although this is not significant, it allows us to infer that the presence of aspartate residue has indirect influence over the selectivity. This finding is reminiscent of the findings of Eisenman84 that anionic field strength of the coordinating ligands can impel the selectivity and supports the claim of Thomas et al.82 that “slight constraints” are indeed necessary to induce appreciable selectivity in the binding sites of ion channel. 4.1.3. Interdependency of Consecutive Binding Sites over Ion Selectvity. Previous experimental observations11 indicated that the presence of contiguous binding sites could act as a decisive factor behind ionic selectivity. In view of this, we have also extended our study to find out the direct infleunce of consecutive binding sites upon selectivity trend. In this regard, we have compared the changes of structural parameters, alteration of the binding free energies of ions and subsequent modulation of ion exchange free energies in the truncated models (TM: contains stand-alone binding sites, Figure 1) with respect to structure B (contains more than one binding sites) (see Tables 1, 2, and 3). Contrary to structure B, the vestibule at the TM model is free from the closest S3 site. Such a miniature model now showed a preference order for the exchange reactions among the alkali metals: K+ > Cs+ > Na+ at ϵ ≥ 10, which otherwise does not act as a favorable binding site in structure B (see Tables 2 and 3). For the metal-bound S3 truncated model, upon optimization only marginal structural changes (except for Na+) were noted (by ∼0.15 Å) with respect to metal-bound structure B (see Table 1). This led to similar binding (Table 2) and ion exchange energies (Table 3) as compared to structure B and thus retains the same selectivity trend. On the contrary, significant geometric changes were observed at the S4 truncated model, where lengthening of the M+-OT63 bond length is noticed with respect to the metal-bound structure B (by ∼0.25 Å). Except for K+, significant shortening was visible in the M− OT63(OH) bond distances (by more than 2 Å) for Na+ and Cs+. Owing to the absence of S3 in the truncated model of S4, Na+ and Cs+ ions are enforced to occupy unfavorable octahedron coordination geometry. These have rendered significant changes in ΔG values and altered the observed selectivity order with respect to structure B, making Cs+ at ϵ ≤ 20 to be the most selective ion, while Na+ has the least binding preference. These observations are analogues with the previous findings on the K+ channel,24,42,47 where a greater possibility of attaining suitable coordination geometry was shown to be the determining factor to arrive at a K+ selectivity over Na+. These results further allow us to infer that consecutive binding sites may have their influence over each other during metal binding. Thus, an oversimplified model with stand-alone binding sites may not truly reflect the actual binding and selectivity trend of ions when compared with the more realistic model, where all of the binding sites act in unison inside the permeating SF.

5. CONCLUSIONS In the present article, we elucidate and rationalize the structure, binding, and translocation of various alkali metal ions through the NaK ion channel in a counterion-free noncompetitive environment. Because of the absence of any interionic 12794

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B electrostatic repulsion, the binding and permeation of a single ion can be imagined to depict a clear picture of the microscopic selectivity in an otherwise nonselective NaK channel. Our combined studies on such a system using both density functional theory and molecular dynamics simulations in explicit water medium led us to summarize the following conclusions. (a) Similar to earlier findings, we also found that the channel is permeable to all three monovalent cations studied in this work. However, the filter is shown to possess strong sitespecific selectivity. Channel forming proteins are relatively flexible structures, more so in the presence of water, that undergo rapid thermal atomic fluctuations often larger than the small difference in ionic radius between the three ions. The signature of site-specific selectivity is, however, found to be so pronounced that even with the inclusion of explicit water in the MtD calculaton, the selectivity trend of the permeating ions remained intact throughout the filter, which corroborates well with the implicit solvent model results of DFT calculations. As explained, conformational arrangement, type of coordinating ligand, and constraints present in the system are found to be the main governing factor to fine-tuning of the selectivity. (b) We found that S3 and S4 are more selective toward K+ and Na+, respectively. K+ is observed to prefer an octahedron coordination environment, whereas Na+ is found to remain favorably bound at the carbonyl plane formed at the junction of S3 and S4. However, the presence of cation at S3 can force Na+ toward the S4 site to attain pyramidal conformation involving water from the cavity. Cs+ is noted to be more favorable at S3 as compared to S4. Further, a marginal selectivity is noted for Na+ over K+ at Sext. The vestibule region is completely nonselective in nature, and the presence of interionic repulsion is essential to prompt ion binding at the vestibule. The optimized structural coordinates and selectivity orders are found to be in excellent agreement with the available crystal structures. (c) As an important conclusion, we find that the use of oversimplified synthetic models without the presence of consecutive binding sites may lead to wrong estimations of the selectivity trend. This explains the indirect influences of nearby binding sites over the ionic selectivity of a particular site. Our EDA analysis suggests that the metal binding is dominated by electrostatic interaction. Nevertheless, orbital interactions may also contribute significantly at the deeply buried sites (S3 and S4). Finally, it is important to emphasize that unlike QM methods the MtD-MD methodology is devoid of several electronic effects such as polarization and charge transfer phenomena, which have been known to affect metal binding. Future studies are ongoing in our laboratory in this direction and will be published in due time. Nevertheless, the present combined investigations using QM and MtD-US-MD methods have provided very useful insights over the bottleneck of the ion translocation phenomenon.

■

■

hydroxyl ligand to metal cations, radial distribution functions and the corresponding integrals of the metal cations, and change of total coordination number of the metal ions during translocation (PDF)

AUTHOR INFORMATION

Corresponding Author

*Phone: +91 22 2559 0300. Fax: +91 22 2550 5151. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS M.S. and T.B. thank Dr. B. N. Jagatap for his kind support and constant encouragement and the BARC computer centre for providing the high performance parallel computing facility. B.S. thanks Dr. Anil Kumar, Shri K. Dodiya, Shri R. Singh, and Dr. K.S. Pradeepkumar for their continuous support and encouragement. The work was supported by DAE under project XII-N-R&D-02.04/Theoretical & Computational Chemistry of Complex Systems.

■

REFERENCES

(1) Hille, B. Ion Channels of Excitable Membranes; Sinauer: Sunderland, MA, 2001. (2) Guharay, F.; Sachs, F. Stretch activated single ion channel currents in tissue cultured embryonic chick skeletal muscle. J. Physiol. 1984, 352, 685−701. (3) Nelson, M. T.; Patlak, J. B.; Worley, J. F.; Standen, N. B. Calcium channels, potassium channels, and voltage dependence of arterial smooth muscle tone. Am. J. Physion-cell Ph. 1990, 259, C13−C18. (4) Doyle, D. A.; Cabral, J. M.; Pfuetzner, R. A.; Kuo, A.; Gulbis, J. M.; Cohen, S. L.; Chait, B. T.; MacKinnon, R. The structure of the potassium channel: molecular basis of K+ conduction and selectivity. Science 1998, 280, 69−77. (5) Jiang, Y.; Lee, A.; Chen, J.; Cadene, M.; Chait, B. T.; MacKinnon, R. Crystal structure and mechanism of a calcium-gated potassium channel. Nature 2002, 417, 515−522. (6) Jiang, Y.; Lee, A.; Chen, J.; Ruta, V.; Cadene, M.; Chait, B. T.; MacKinnon, R. X-ray structure of a voltage-dependent K+ channel. Nature 2003, 423, 33−41. (7) Payandeh, J.; Scheuer, T.; Zheng, N.; Catterall, W. A. The crystal structure of a voltage-gated sodium channel. Nature 2011, 475, 353− 358. (8) Zhang, X.; Ren, W.; DeCaen, P.; Yan, C.; Tao, X.; Tang, L.; Wang, J.; Hasegawa, K.; Kumasaka, T.; He, J.; et al. Crystal structure of an orthologue of the NaChBac voltage-gated sodium channel. Nature 2012, 486, 130−134. (9) Heginbotham, L.; Abramson, T.; MacKinnon, R. A functional connection between the pores of distantly related ion channels as revealed by mutant K+ channels. Science 1992, 258, 1152−1155. (10) Heginbotham, L.; Lu, Z.; Abramson, T.; MacKinnon, R. Mutations in the K+ channel signature sequence. Biophys. J. 1994, 66, 1061−1067. (11) Derebe, M. G.; Sauera, D. B.; Zenga, W.; Alam, A.; Shia, N.; Jiang, Y. Tuning the ion selectivity of tetrameric cation channels by changing the number of ion binding sites. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 598−602. (12) Shi, N.; Ye, S.; Alam, A.; Chen, L.; Jiang, Y. Atomic Structure of Na+ - and K+ -conducting channel. Nature 2006, 440, 570−574. (13) Alam, A.; Jiang, Y. Structural analysis of ion selectivity in the NaK channel. Nat. Struct. Mol. Biol. 2009, 16, 35−41. (14) Alam, A.; Jiang, Y. High-resolution structure of the open NaK channel. Nat. Struct. Mol. Biol. 2009, 16, 30−34.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b05996. Choice of basis set, calculated hydration energy of alkali metal ions, effects of water molecules around metal ions in the site, Sext, charge transfer from carboxyl and 12795

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B (15) Noskov, S. Y.; Roux, B. Importance of Hydration and Dynamics on the Selectivity of the KcsA and NaK Channels. J. Gen. Physiol. 2007, 129, 135−143. (16) Vora, T.; Bisset, D.; Chung, S.-H. Conduction of Na+ and K+ through the NaK Channel: Molecular and Brownian Dynamics Studies. Biophys. J. 2008, 95, 1600−1611. (17) Shen, R.; Guo, W. Ion binding properties and structure stability of the NaK channel. Biochim. Biophys. Acta, Biomembr. 2009, 1788, 1024−1032. (18) Shen, R.; Guo, W.; Zhong, W. Hydration valve controlled nonselective conduction of Na+ and K+ in the NaK channel. Biochim. Biophys. Acta, Biomembr. 2010, 1798, 1474−1479. (19) Furini, S.; Domene, C. Gating at the selectivity filter of ion channels that conduct Na+ and K+ Ions. Biophys. J. 2011, 101, 1623− 1631. (20) Furini, S.; Domene, C. Nonselective conduction in a mutated NaK channel with three cation-binding Sites. Biophys. J. 2012, 103, 2106−2114. (21) Furini, S.; Domene, C. K+ and Na+ Conduction in selective and nonselective ion channels via molecular dynamics simulations. Biophys. J. 2013, 105, 1737−1745. (22) Shen, R.; Guo, W. Mechanism for variable selectivity and conductance in mutated NaK channels. J. Phys. Chem. Lett. 2012, 3, 2887−2891. (23) Wang, Y.; Chamberlin, A. C.; Noskov, S. Y. Molecular strategies to achieve selective conductance in NaK channel variants. J. Phys. Chem. B 2014, 118, 2041−2049. (24) Dudev, T.; Lim, C. Determinants of K+ vs Na+ selectivity in potassium channels. J. Am. Chem. Soc. 2009, 131, 8092−8101. (25) Dudev, T.; Lim, C. Factors governing the Na+ vs K+ selectivity in sodium ion channels. J. Am. Chem. Soc. 2010, 132, 2321−2332. (26) Varma, S.; Rempe, S. B. Structural transitions in ion coordination driven by changes in competition for ligand binding. J. Am. Chem. Soc. 2008, 130, 15405−15419. (27) For a review see: Zhu, Y.-G.; Smolders, E. Plant uptake of radiocaesium: a review of mechanisms, regulations and application. J. Exptl. Botany 2000, 51, 1635−1645. (28) Fowler, P. W.; Tai, K.; Sansom, M. S. P. The selectivity of K+ ion channels: testing the hypotheses. Biophys. J. 2008, 95, 5062−5072. (29) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (30) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (31) Schafer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. (32) Sundararajan, M.; Sinha, V.; Bandyopadhyay, T.; Ghosh, S. K. Can functionalized cucurbituril bind actinylcations efficiently? A density functional theory based investigations. J. Phys. Chem. A 2012, 116, 4388−4395. (33) Sundararajan, M.; Ganyushin, D.; Ye, S.; Neese, F. Multireference ab initio studies of zero-field splitting and magnetic circular dichroism spectra of tetrahedral Co(II) complexes. Dalton Trans. 2009, 30, 6021−6036. (34) Sundararajan, M. Quantum chemical challenges for the binding of simple alkanes to supramolecular hosts. J. Phys. Chem. B 2013, 117, 13409−13417. (35) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (36) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104−154119. (37) Becke, A. D. Density functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5682.

(38) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (39) Schafer, A.; Huber, C.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (40) Klamt, A.; Schuurmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (41) Ahlrichs, R.; Bar, M.; Baron, H.-P.; Bauernschmitt, R.; Bocker, S.; Ehrig, M.; Eichkorn, K.; Elliot, S.; Furche, F.; Haase, F.; Haser, M.; Horn, H. et al.; TURBOMOLE, V6.3 2011, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007; TURBOMOLE GmbH, since 2007; available from http://www. turbomole.com. (42) Bostick, D. L.; Brooks, C. L., III Selectivity in K+ channels is due to topological control of the permeant ion’s coordinated state. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 9260−9265. (43) Varma, S.; Sabo, D.; Rempe, S. B. K+/Na+ selectivity in Kchannels and Valinomycin: Over-coordination Vs Cavity-size constraints. J. Mol. Biol. 2008, 376, 13−22. (44) Reed, A. E.; Weinstock, R. B.; Weinhold, F. Natural population analysis. J. Chem. Phys. 1985, 83, 735−746. (45) Liao, R. Z.; Siegbahn, P. E. Mechanism and selectivity of the dinuclear iron benzoyl-coenzyme A epoxidase BoxB. Chemical Science 2015, 6, 2754−2764. (46) Dudev, T. Modeling metal binding sites in proteins by quantum chemical calculations. Computational Chemistry 2014, 2, 19−21. (47) Varma, S.; Rempe, S. B. Tuning ion coordination architectures to enable selective partitioning. Biophys. J. 2007, 93, 1093−1099. (48) Morokuma, K. Molecular Orbital Studies of Hydrogen Bonds. III. CO···H-O Hydrogen Bond in H2CO···H2O and H2CO···2H2O. J. Chem. Phys. 1971, 55, 1236−1244. (49) Ziegler, T.; Rauk, A. On the calculation of bonding energies by the Hartree Fock Slater method. Theor. Chim. Acta 1977, 46, 1−10. (50) Van Lenthe, E.; van Leeuwen, R.; Baerends, E. J.; Snijders, J. G. Regular relativistic two-component Hamiltonians. Int. J. Quantum Chem. 1996, 57, 281−293. (51) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Gisbergen, van S. J. A.; Fonseca Guerra, C.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931−967. (52) Zhang, Y.; Voth, G. A. Combined metadynamics and umbrella sampling method for the aalculation of ion permeation free energy profiles. J. Chem. Theory Comput. 2011, 7, 2277−2283. (53) Lindorff-Larsen, K.; Maragakis, P.; Piana, S.; Eastwood, M. P.; Dror, R. O.; Shaw, D. E. Systematic validation of protein force fields against experimental data. PLoS One 2012, 7, e32131. (54) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular dynamics. J. Chem. Theory Comput. 2008, 4, 435−447. (55) Jorgensen, W. L.; Jenson, C. Temperature dependence of TIP3P, SPC, and TIP4P water from NPT Monte Carlo simulations: seeking temperatures of maximum density. J. Comput. Chem. 1998, 19, 1179−1186. (56) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577−8593. (57) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (58) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (59) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: a portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commun. 2009, 180, 1961− 1972. 12796

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797

Article

The Journal of Physical Chemistry B (60) Moradi, M.; Tajkhorshid, E. Driven Metadynamics: Reconstructing equilibrium free energies from driven adaptive-bias simulations. J. Phys. Chem. Lett. 2013, 4, 1882−1887. (61) Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the accuracy of metadynamics. J. Phys. Chem. B 2005, 109, 6714−6721. (62) Babin, V.; Roland, C.; Darden, T. A.; Sagui, C. The free energy landscape of small peptides as obtained from metadynamics with umbrella sampling corrections. J. Chem. Phys. 2006, 125, 204909. (63) Autieri, E.; Sega, M.; Pederiva, F.; Guella, G. Puckering free energy of pyranoses: a NMR and metadynamics-umbrella sampling investigation. J. Chem. Phys. 2010, 133, 095104. (64) Pathak, A. K.; Bandyopadhyay, T. Unbinding free energy of acetylcholinesterase bound oxime drugs along the gorge pathway from metadynamics-umbrella sampling investigation. Proteins: Struct., Funct., Genet. 2014, 82, 1799−1818. (65) Varma, S.; Rempe, S. B. Coordination numbers of alkali metal ions in aqueous solutions. Biophys. Chem. 2006, 124, 192−199. (66) Rao, J. S.; Dinadayalane, T. C.; Leszczynski, J.; Sastry, G. N. Comprehensive study on the solvation of mono-and divalent metal cations: Li+, Na+, K+, Be2+, Mg2+ and Ca2+. J. Phys. Chem. A 2008, 112, 12944−12953. (67) Soper, A. K.; Weckstrom, K. Ion solvation and water structure in potassium halide aqueous solutions. Biophys. Chem. 2006, 124, 180− 191. (68) Glezakou, V.-A.; Chen, Y.; Fulton, J. L.; Schenter, G k.; Dang, L. X. Electronic structure, statistical mechanical simulations, and EXAFS spectroscopy of aqueous potassium. Theor. Chem. Acc. 2006, 115, 86− 99. (69) Whitfield, T. W.; Verma, S.; Harder, E.; Lamoureux, G.; Rempe, S. B.; Roux, B. Theoretical study of aqueous solvation of K+ comparing ab initio, polarizable, and fixed-charge models. J. Chem. Theory Comput. 2007, 3, 2068−2082. (70) Carrillo-Tripp, M.; Saint-Martin, H.; Ortega-Blake, I. A comparative study of the hydration of Na+ and K+ with refined polarizable model potentials. J. Chem. Phys. 2003, 118, 7062−7073. (71) Sadhu, B.; Sundararajan, M.; Bandyopadhyay, T. Water mediated differential binding of strontium and cesium cations in fulvic acid. J. Phys. Chem. B 2015, 119, 10989−10997. (72) Kim, I.; Allen, T. W. On the selective ion binding hypothesis for potassium channels. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 17963− 17968. (73) Kopfer, D. A.; Song, C.; Gruene, T.; Sheldrick, G. M.; Zachariae, U.; de Groot, B. L. Ion permeation in K+ channels occurs by direct coulomb knock-on. Science 2014, 346, 352−355. (74) Smart, O. S.; Goodfellow, J. M.; Wallace, B. A. The Pore Dimensions of Gramicidin A. Biophys. J. 1993, 65, 2455−2460. (75) Bostick, D. L.; Brooks, C. L., III Selective complexation of K+ and Na+ in simple polarizable ion-ligating systems. J. Am. Chem. Soc. 2010, 132, 13185−13187. (76) Bostick, D. L.; Arora, K.; Brooks, C. L., III K+/Na+ selectivity in toy cation binding site models is determined by the ‘host’. Biophys. J. 2009, 96, 3887−3896. (77) Thomas, M.; Jayatilaka, D.; Corry, B. Mapping the importance of four factors in creating monovalent ion selectivity in biological molecules. Biophys. J. 2011, 100, 60−69. (78) Yu, H.; Noskov, S. Y.; Roux, B. Two mechanisms of ion selectivity in protein binding sites. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 20329−20334. (79) Noskov, S. Y.; Berneche, S.; Roux, B. Control of ion selectivity in potassium channels by electrostatic and dynamic properties of carbonyl ligands. Nature 2004, 431, 830−834. (80) Forman, S. L.; Fettinger, J. C.; Pieraccini, S.; Gottarelli, G.; Davis, J. T. Toward artificial ion channels: A Lipophilic G-Quadruplex. J. Am. Chem. Soc. 2000, 122, 4060−4067. (81) Dudev, T.; Lim, C. Metal binding affinity and selectivity in metalloproteins: insights from computational studies. Annu. Rev. Biophys. 2008, 37, 97−116.

(82) Thomas, M.; Jayatilaka, D.; Corry, B. The predominant role of coordination number in potassium channel selectivity. Biophys. J. 2007, 93, 2635−2643. (83) Alam, A.; Shi, N.; Jiang, Y. Structural insight into Ca2+ specificity in tetrameric cation channels. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 15334−15339. (84) Eisenman, G. Cation selective glass electrodes and their mode of operation. Biophys. J. 1962, 2, 259−323.

12797

DOI: 10.1021/acs.jpcb.5b05996 J. Phys. Chem. B 2015, 119, 12783−12797