Accurate Evaluation of Ion Conductivity of the Gramicidin A Channel

May 12, 2016 - Na+ permeating through the gramicidin A channel are characterized by using the AMOEBA polarizable force field with a total sampling tim...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/JCTC

Accurate Evaluation of Ion Conductivity of the Gramicidin A Channel Using a Polarizable Force Field without Any Corrections Xiangda Peng,†,‡ Yuebin Zhang,† Huiying Chu,† Yan Li,† Dinglin Zhang,† Liaoran Cao,† and Guohui Li*,† †

Laboratory of Molecular Modeling and Design, State Key Laboratory of Molecular Reaction Dynamics, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian 116023, P. R. China ‡ Chinese Academy of Science, University of Chinese Academy Sciences, Beijing 100049, P. R. China S Supporting Information *

ABSTRACT: Classical molecular dynamic (MD) simulation of membrane proteins faces significant challenges in accurately reproducing and predicting experimental observables such as ion conductance and permeability due to its incapability of precisely describing the electronic interactions in heterogeneous systems. In this work, the free energy profiles of K+ and Na+ permeating through the gramicidin A channel are characterized by using the AMOEBA polarizable force field with a total sampling time of 1 μs. Our results indicated that by explicitly introducing the multipole terms and polarization into the electrostatic potentials, the permeation free energy barrier of K+ through the gA channel is considerably reduced compared to the overestimated results obtained from the fixed-charge model. Moreover, the estimated maximum conductance, without any corrections, for both K+ and Na+ passing through the gA channel are much closer to the experimental results than any classical MD simulations, demonstrating the power of AMOEBA in investigating the membrane proteins.



the fixed-charge force field was found to have problems in the heterogeneous system, especially, in the membrane protein system such as ion transduction, in which the electronic polarization interactions play an essential role in modulating the dynamical properties and reproducing the experimental observations. The next generation of polarizable force fields provides a systematic improvement over the classical force fields by taking into account the electrostatic interaction changes of many-body systems in response to the environment.11 For example, the AMOEBA (Atomic Multipole Optimized Energetic for Biomolecular Applications) developed by Ponder and coworkers12−15 is one of the most widely tested polarizable force fields, which introduces the permanent dipoles and quadrupoles to accurately reproduce the molecular electrostatic potentials. Although there already are lines of evidence indicating significant improvements of the AMOEBA over the nonpolarizable models have been achieved in the case of solvation free energy calculations for ions16−18 and small molecules19,20 and binding free energy calculations for protein−ligand complexes,21,22 the implementation of AMOEBA in membrane proteins has been hindered due to the absence of the polarizable force field parameters for lipids. Recently, we have developed the polarizable force field parameters of lipid molecules,23 for the first time, based on the

INTRODUCTION Membrane proteins play an essential role in cellular metabolism, transportation, and signal transduction across the biological membrane, which also account for more than 50% of current medicinal targets.1 However, due to the experimental challenges in conducting biochemical and structural studies on membrane proteins, computational methods have become an indispensable tool to investigate their structural and functional relationship and to unveil the underlying mechanisms, such as ion conduction, drug recognition, material transporting, and conformational changes of the membrane protein in single transduction.2−5 In particular, molecular dynamic simulations have been used extensively to study the membrane protein systems such as the gA channel,6−9 the K+ channel,6,8,9 and the AM2 channel10 and gained considerable biophysical insights of these complexes at an atomic level. In molecular mechanics, the dynamic behaviors of atoms in a MD simulation system are described by solving the Lagrangian equations of motion of each individual atom, and their interactions are treated by well-parametrized empirical potential energy functions, namely the force field, into the bonded and nonbonded terms, respectively. Unlike the high level ab inito quantum mechanical (QM) methods, the classical force field adopts a fixed-charge model to account for the electronic interactions in the many-body system, therefore, reducing the computational cost and allowing the simulations of biological systems of larger sizes and longer time scales. However, despite its huge achievements in homogeneous and bulk-like systems, © XXXX American Chemical Society

Received: February 3, 2016

A

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. Comparison of the MD simulations results of the pure DMPC bilayer using the AMOEBA based polarizable force field for the experiments: (A) the area per DMPC lipid value; (B) the volume per DMPC lipid; (C) the electron density profile of each segment of DMPC; (D) the order parameters of the tail carbon atoms.

improve the classical MD simulation results, such as by taking into account the membrane size and polarization effects.25−27,35 On the other hand, Sandeep and co-workers used the charge equilibration polarizable force field to evaluate the PMF for K+ permeation in the gA channel within a small 20 DMPC bilayer. By assuming a constant K+ diffusion coefficient in the channel (about one-tenth of the bulk value),39 they reported a maximum conductance of 57 pS for K+ permeation, which is in semiquantitative agreement with the experimental value (21pS).39 In this work, the AMOBEA polarizable force field was employed to systematically characterize the PMF profiles for ions permeation through the gA channel using extensive umbrella sampling simulation (a total of 1 μs) for K+ and Na+ in a larger 96-DMPC bilayer, without any correction. The diffusion coefficient of K+ and Na+ along the gA channel was calculated using the GLE-HO method to estimate the maximum conductance of K+/Na+ ions with the 1D NernstPlank equation. Our calculated PMFs suggested that, by explicitly introducing the multipole terms and polarization into the electrostatic potentials, the permeation free energy barrier of K+ through the gA channel is considerably reduced compared to the fixed-

same framework as the AMOEBA model, which opens an avenue to study the membrane proteins using the AMOEBA polarizable force field. As a benchmark model to test the MD simulations of membrane proteins, the gramicidin A (gA) channel has been extensively used to study the mechanisms of ion permeation owing to its simplicity.6−8,24 However, it has long been recognized that classical force fields face significant challenges in accurately reproducing and predicting experimental observables such as ion conductance and permeability even in the small gA channel25−28 due to the dramatic changes of the ion environment during permeation and the lack of an accurate description of the electrostatic interactions in response to the environment changes in the classical force field. The single channel conductance is overestimated by several orders of magnitude compared to the experiment observation, depending on the classical force field and MD simulation protocol used. For example, previous calculated potential mean forces (PMF) for K+ permeation using the classical force field demonstrated an unexpectedly high free energy barrier of 10−20 kcal/ mol,24−36 whereas the experimental data indicated that the K+ permeates through the gA channel with a small free energy barrier.37,38 Therefore, there are many corrections utilized to B

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

the previous set (amoeba09.prm) derived from liquid simulations of small organic molecules containing similar functional groups.13,19 PMF Calculations. An effective and appropriate method to deal with insufficient sampling is to use the umbrella sampling (US) algorithm51,52 by adding a simple spring potential to restrain the motions of ions at some dispersed positions along a predefined path. After sampling, the biased probability distribution can be reweighted using some methods such as the Weighted Histogram Analysis Method (WHAM),53 the Umbrella Integration method (UI),54 the multistate Bennett acceptance ratio method (MBAR),55,56 the variational Free Energy Profile (vFEP) method,57 and others.58−61 WHAM is widely used and works well when the knowledge of overlap between umbrella windows is accessible. So we used the umbrella sampling method and the WHAM algorithm to get the PMFs of K+/Na+ permeation. For the US simulations of the gA channel, the Andersen Thermostat and the Monte Carlo Barostat were applied to couple the system’s temperature and pressure to 330 K with a collision frequency of 1.0 ps−1 and 1 bar with a trial frequency to change the box every 1000 MD steps, respectively. All hydrogen bond lengths and the H−O−H angles of water were fixed using the RATTLE algorithm with an integration step size of 2 fs. The convergence threshold of the induced dipole in the SCF iteration was set to 0.0005D. Other MD parameters are the same with that used in pure DMPC bilayer MD simulation. 10 ns NPT equilibration of the gA channel embedded in the DMPC bilayer was performed. The average Root Mean Square Deviation (RMSD) between the simulation structure and the NMR structure is 0.632 Å (Figure S1.A in the Supporting Information). The tilt angle of the gA channel fluctuated within 0 to 12 degrees (Figure S1.B in the Supporting Information). Then, the last structure of equilibration was used to carry out the umbrella sampling (US) simulations. 6 ns US simulations were conducted for each window, and the last 5.0 ns trajectories were used for analysis. In the potassium simulation, samples were restrained by a force constant of 10 kcal/mol/Å2 around 72 windows of 0.5 Å each to cover the total −18 Å to 18 Å range of the reaction coordinate. On the other hand, for the sodium simulation, 28 windows of 0.5 Å each were used to cover the region ranging from −18 Å to −11 Å as well as 11 Å to 18 Å with the force constant set to be 10 kcal/mol/Å2, while the region between −11 and 11 Å was covered by 73 windows of 0.3 Å each with the force constant set to 15 kcal/mol/Å2. Besides, a cylindrical restraint with a 5 kcal/mol/Å2 force constant was applied on the same ion to keep it within 5 Å from the centerline of the gA channel. All simulations were employed using the PLUMED free energy calculation library,62,63 which has been incorporated into the OpenMM package42 in our previous work64 to take advantage of GPU-accelerated MD simulation and enhanced samplings. The g_wham module65 in GROMACS66,67 was used to reweight the average distribution function of the permeation ions in the gA channel. A sampling bin size of 0.18 Å and tolerance of 2.4 × 10−7 kcal/mol were used to check the convergence. Free energies at a distance of 18 Å from the gA dimer center were chosen as reference. Symmetrized PMF profiles were obtained by averaging the density of states between both sides of the profile. For calculation of the 2D PMF with respect to the ion position along channel axis z and the distance r from the pore

charge model. Moreover, the calculated single channel conductance for K+ using AMOEBA reaches the same magnitude compared to the experimental observation, demonstrating the power of AMOEBA in predicting properties of the membrane protein system. In addition, the selectivity of the gA channel for K+ over Na+ is also discussed. Our results indicated that the Na+ exhibits a stronger interaction with carbonyl groups lining on the gA channel than that of K+, which damps the permeation of Na+.



METHODS Pure DMPC Bilayer MD Simulation Using the AMOEBA Based Polarizable Force Field for Lipid. Recently, we have developed the AMOEBA based polarizable force field for the POPC molecule.23 Here, the parametrization for the DMPC lipid adopts the same protocol as our previous work. To further testify the validation of our parameters for the DMPC molecule, 100 ns MD simulation of the pure DMPC bilayer was performed using our AMOEBA based polarizable force field for lipid. The simulation with the polarizable force field accurately reproduces the experimental observations, such as area per lipid, volume per lipid, order parameters, and the electron density profile, indicating the validation of our parameters for DMPC (Figure 1). The initial conformation of the pure DMPC bilayer containing 72 DMPC molecules and 2824 AMOEBA03 water molecules was built using the CHARMM-GUI40 Web server. The MD simulation was performed under the NPT ensemble at 303.5 K and 1 bar with an Andersen Thermostat41 and a Monte Carlo anisotropic barostat implemented in the OpenMM package.42 The Rattle algorithm43 was also adopted in the MD to constrain all bonds involving the hydrogen, which ensures the stability of the system with a 2 fs integration time step. The Particle-Mesh Ewald (PME)44,45 method for longrange electrostatic calculations was employed, and the cutoff was set to 10.0 Å. The cutoff for the nonbonded van der Waals interactions was set to 12.0 Å, and the dispersion correction algorithm46 was used to estimate the contribution of all van der Waals interactions beyond the cutoff. Mutual polarization was used, which employs the self-consistent field (SCF) with the Direct Inversion in Iterative Subspace (DIIS) method47 to calculate the induced dipole moment in every integration. The convergence threshold of the induced dipole in SCF iteration was set to 0.00001D. MD Simulation of the gA Channel Using the AMOEBA Polarizable Force Field. The initial head-to-head gA dimer conformation was taken from the Protein Data Bank (1JNO). Then, the gA dimer was embedded into a 96 DMPC bilayer and solvated with 3025 AMOEBA03 waters using the CHRMM-GUI Web server. 1.0 M KCl/NaCl were used to balance the charges of the whole system to neutral. The force field parameters for water, K+, Na+, and Cl− ions were taken from AMOEBA potential (amoebapro13.prm)14 in the TINKER package.48 The parameters for the D-amino acids of gA were obtained from the same potential by referencing the respective L-amino acid. The electrostatic parameters of formyl and ethanolamide groups in gA were derived from quantum mechanical calculation (geometry optimization at HF/6-31G*, single point calculation at MP2/6-311++G(2d,2p) with the Gaussian0949 package; multipoles calculation with the GDMA program50) through a similar process reported in ref 14. The van der Waals (vdW), bond, angle, and atomic polarizability parameters of formyl and ethanolamide groups were taken from C

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 2. Asymmetrical and symmetrical PMF profiles (upper panel) and the time evolution of the symmetrical profiles (lower panel).

center, the histograms of r at the same z position were computed, and the probabilities ρ(r, z) were used to obtain the free energy according to −kBT ln ρ(r, z). Then, the 2D PMFs were estimated using the following equation W2d(r , z) = W1d(z) − kBT ln ρ(r , z)

In the GLE-HO method, the memory function is associated with the velocity autocorrelation function in terms of a Generalized Langevin Equation (GLE) kBT 2

(1)

⟨z ̇ ⟩i

where W is the PMFs along reaction coordinates, kB is the Boltzmann constant, T is the thermodynamic temperature, and ρ(z) is the unbiased probability density at reaction coordinate (r, z). Ion Coordination Number. The coordination numbers were determined by integrating the Radial Distribution function (RDF),83 g(r)



nc = 4πρ

∫o

R

r 2g (r )dr

C ̇ (t , zi ) = −

∫0

kBT

∫0

2

⟨δz ⟩i

t

dt ′C(t ′, zi)

t

dt ′M(t − t ′, zi)C(t ′, zi)

(3)

where C(t, zi) is the normalized velocity autocorrelation function of the motion of ions under a harmonic oscillator in our umbrella sampling simulations at each zi bin. Through Einstein’ relation, one can use the Laplace transform of memory function to get the diffusion coefficient at zi D(zi) = lim

(2)

s→0

kBT ̂ M (s , zi )

(4) 76

where nc is the coordination number, ρ is the density of ligand atoms, and R is the first minimum determined from the ionoxygen radial RDF in bulk. The value of R is 3.2 Å for Na+ and 3.6 Å for K+ (Figure S2 in the Supporting Information). Diffusion Coefficients. According to the 1D Nernst-Plank equation, the local diffusion constant of K+/Na+ inside the gA channel is required to estimate the channel conductance. A variety of methods can be used to handle this task, such as Mean Square Displacement (MSD),68−71 Velocity Autocorrelation Function (VACF),68,71,72 Second Fluctuation Dissipation Theorem (SFDT),73,74 and analysis of the Generalized Langevin Equation for a Harmonic Oscillator (GLEHO),25,75,76 etc. All the methods were tested by Artem and partners for K+ inside gA. They found that MSD and VACF methods were not reliable owing to the effect of the systematic force from the membrane-channel system to the ion.77 Therefore, in this work, the GLE-HO method was employed to calculate the local diffusion constant of ions.

Another convenient form

of this method can be expressed

as D(zi) = lim

s→0

−Ĉ(s , zi)⟨δz 2⟩i ⟨z 2̇ ⟩i Ĉ(s , zi)[s⟨δz 2⟩i + ⟨z 2̇ ⟩i /s] − ⟨δz 2⟩i ⟨z 2̇ ⟩i (5)

Ĉ (s, zi) is the Laplace transformation of C(t, zi) C ̂ (s , zi ) =

∫0



dt e−st C(t , zi)

(6)

The Maximum Conductance. Previously, Roux and coworkers25 have provided an elegant protocol to calculate the maximum conductance of ion passing through the gA channel using the 1D Nernst-Plank equation78 J = −D(z) D

dP(z) D(z) dW (z) − P(z) dz kBT dz

(7)

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 3. Theoretical and experimental PMFs for (A) K+ ion (top) and Na+ ion (bottom) permeation in the gA channel. The results of CHARMM27 for K+ and Na+ were obtained from refs 25 and 27. (B) A snapshot shows the K+ ion located at the center between two monomers of gA, which is surrounded by several carbonyl oxygen atoms and two water molecules. The figure was prepared using the VMD program.89

where J is the net stationary flux, D(z) is the ion diffusion coefficient as a function of channel axis z, P(z) is the probability density per unit length of finding a ion, and W(z) is the potential mean force as a function of z. Based on the small-field approximation,79 the maximum conductance of the single-file channel is calculated using the following equation78−80 gmax =

e2 ⟨D(z)−1e+W (z)/ kBT ⟩−1⟨e−W (z)/ kBT ⟩−1 kBTL2

(8)

where L is the length (28 Å) that the 1D PMF of a single ion region is meaningful in the channel (|z| < 15 Å), and e is the charge quantity of one K+/Na+ ion.



RESULTS AND DISCUSSION Potential Mean Forces (PMF) for K+/Na+ Permeation through the gA Channel Using the AMOEBA Polarizable Force Field. In order to ensure the convergence of the PMF profiles calculated using the AMOEBA polarizable force field, extensive samplings were conducted in this work, and a total of 0.4 μs for K+ and 0.6 μs for Na+ umbrella sampling simulations were carried out. However, the asymmetric of the PMF profiles was still observed due to the simulation artifacts caused by the finite size of the periodic system,25,26 the slow fluctuations of the membrane,81,82 and gA titling during MD simulations.30,83,84 At the end of our 6 ns umbrella sampling simulations of each bin, the asymmetric PMF exhibits a deviation less than 2 kcal/mol on the two sides (Figure 2 A-B). Following previous work,26 the PMF profiles were symmetrized by averaging the density of states between both sides of the profile. The Root Mean Square Deviation (RMSD) between the symmetrized PMFs of 5 and 6 ns US simulations is less

Figure 4. 2D PMFs for (A) K+ and (B) Na+. The second dimension of 2D PMFs is the vertical distance from the central axis of the channel.

than 0.2 kcal/mol for both K+ and Na+ (Figure 2 C−D and Figure S3 in the Supporting Information), indicating the convergence of the PMFs. In Figure 3, the symmetrized PMF profiles obtained using the AMOEBA polarizable force field were compared to those calculated using the CHARMM27 classical force field.25,27 As can be seen from Figure 3, the PMF profiles for both K+ and Na+ permeation through the gA channel are significantly reduced compared to the classical force field because of the explicitly defined polarization effects into the electrostatic potentials in AMOEBA. The calculated free E

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 5. Mean force decomposition of PMFs. The total permeation PMFs (black lines) for (A) K+ and (B) Na+ were decomposed into bulk water (pink lines), DMPC + bulk ions (green lines), the gA channel (red lines), and water within a single-file region (blue lines).

energy barrier of K+ permeation through the gA channel using AMOEBA is 4.8 kcal/mol with respect to the binding site, whereas the free energy barrier for Na+ is 5.6 kcal/mol. These results are in agreement with the selectivity of the gA channel.38 2D-PMF Indicates a Closer Interaction of Na+ with the Carbonyl Groups Lining along the Channel. In order to provide the energetics of the ion permeation perpendicular to the channel axis, we also depicted the 2D-PMFs as a function of the position of the ion along the channel axis and the radial distance of the ion from the center of the pore. It can be clearly seen from Figure 4 that the entrance of the channel is located at |z| = 11 Å on the landscapes, above which (|z| > 11 Å), the 2DPMF reveals nearly a flat landscape in the bulk. In the case of K+, two outer and inner cation-binding sites are observed near the pore entrance at |z| = 12 Å and |z| = 9 Å, respectively. The positions of the cation-binding sites observed in AMOEBA are consistent with the experimental observation and previous calculated results.29,79,85 In the case of Na+, the 2D-PMF displays a deeper binding well at |z| = 9 Å, while at |z| = 7 Å, Na+ spans a wider region perpendicular to the pore axis, indicating the closer interactions with the carbonyl groups. Mean Force Decomposition. Since the free energy change is equivalent to the thermodynamic reversible work, the total PMF profile can be further decomposed into individual components of forces38

Figure 6. Distribution of water dipole moments in bulk water (top) and gA (bottom). The red lines show the dipole moments in the first hydration shell of a K+ ion (green line for Na+ ion), and the dashed lines show those of water in the absence of an ion.

Wα(z) − Wα(z 0) = −∑ α

∫z

z

dz′⟨Fα(z′)⟩ 0

(9)

Fα can be any microscopic force which is independent of each other, and the sum is equal to the total force on the reaction coordinate. Generally, the PMF decomposition can be accomplished in many ways. In this work, to gain a better understanding of the mechanism of ion permeation through the gA channel, the total PMF profile was decomposed into the contributions from bulk water, lipid + bulk ions, gA, and water within the single-file region, respectively. As shown in Figure 5, by decomposing the PMF profile in this manner, the ion movement through the gA channel could be clearly divided into a series of discrete steps. (I) Ion entry: This process is facilitated by the attractive contributions from both DMPC and gA. As mentioned above, the total PMF profile exhibits two cation-binding wells at the outer and inner sites of the entrance. Here, we can easily identify that the outer binding site at |z| =

Figure 7. Coordination number changes of (A) K+ and (B) Na+ with oxygen atoms. The blue and green lines represent the coordination number from water oxygen atoms and carbonyl oxygen. The red line refers to the total coordination numbers.

F

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 8. Exchanges of coordinating carbonyl oxygen during ion permeation. O0 refers to the oxygen atom on the formyl group of gA. O2, O4, O5, O6, O7, O9, and O11 refer to the carbonyl oxygen of Gly 2, D-Leu 4, Ala 5, D-Val 6, Val 7, Trp 9, and Trp11 on the backbone of gA, respectively.

cation facilitate the ion diffusion to the channel entrance; meanwhile, the gA channel provides a stronger trapping effect to the alkali metal cation in the deeper channel at |z| = 9 Å. (II) Ion translocation through the single-file region: In this region, the ions were totally shielded from the bulk water. Thus, the bulk water arises a free energy plateau about 71.2 kcal/mol for K+ and 89.1 kcal/mol for Na+. These two values correspond to the total dehydration free energy of K+ (72.6 kcal/mol) and Na+ (89.9 kcal/mol) in bulk water calculated using the AMOEBA polarizable force field as previously reported by Ponder and co-workers.16 Accordingly, the ion permeation could be regarded as the compensation of the ion dehydration free energy contributed by the rest of each individual component during the transfer. Particularly, the single-file water molecules contribute a significant effect in stabilization of the ions, which compensates nearly half of the dehydration free energy of ions. (III) Ion exit: It can be considered as a reverse process of ion entry, corresponding to the dissociation from the pore and rehydration of ions into the bulk solution. Another decomposition was employed to identify the effect of each individual nonbonded term to the interactions between permeation ions and gA, water, DMPC, and ions in bulk (See Figure S4 in the Supporting Information for details.). Single-File Water Polarization Stabilizes the Ions in the Pore. The explicit description of polarization effects in the AMOEBA polarizable force field significantly improves the accuracy of PMF calculations in the case of ion permeation through the gA channel. One of the major reasons for the improvement may be attributed to the flexible water dipole moment changes along with the environmental variations. As shown in Figure 6, when the ions are present in the gA channel, the neighboring water dipole moments will be enhanced from 2.59D to 2.70D for K+ and 2.65D to 2.76D for Na+, respectively. Since the ion−water interactions are proportional to the dipole moments, the ions would be further stabilized by

Figure 9. Diffusion coefficients for K+ (blue) and Na+ (red) passing through the gA channel.

Table 1. Comparison of the Calculated Ion Conductances Using Various Force Fields with Experiments K+

conductance (pS) experimental AMOEBA CHARMM27 CHARMM27 with correction AMBER AMBER with correction GROMOS GROMOS with correction

27

87

23.8, 21 36.8 4.7e-3,25 1.54,26 0.012727 9.27,26 1.1727 6.926 28.126 0.006926 0.0526

Na+ 12.4,27 1588 14.7 6.0e-427 0.07427

12 Å is contributed from the DMPC headgroup, and the inner binding site at |z| = 9 Å is provided by gA protein. The longrange charge attractions between the DMPC head groups and G

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

through the gA channel, our calculated PMF using AMOEBA provides an encouraging agreement with experiments, indicating the accuracy of the AMOEBA polarizable force field and its power in predicting experiment observation. Moreover, the calculated, without any corrections, single-channel conductances of the gA channel for both K+ and Na+ reach the same order of magnitude with experiment results, which are closer to the experimental measurements than any previous calculations. The results demonstrated here would provide a promising avenue to improve the accuracy of MD simulations in membrane proteins.

the coordinating water molecules in the channel compared to the bulk solution, thus, reducing the permeation free energy barrier. However, in the nonpolarizable force field, this effect is not taken into account. Therefore, free energy barriers obtained using the fixed-charge force fields would be significantly overestimated. Ion Hydration/Dehydration Rearrangement during Permeation. Besides the single-file water, the protein matrix of gA, particularly for the backbone carbonyls, also plays an essential role in stabilization of the ions during the ion permeation along a narrow channel. As shown in Figure 7, the average water coordination numbers are 6.89 for K+ and 5.69 for Na+ in the bulk water, whereas the water coordination numbers of both K+ and Na+ decrease dramatically to 2 when the ions enter the single-file channel. However, the total coordination numbers would be still maintained to an average value of 6.11 for K+ and 5.39 for Na+ during the permeation because of the interactions of ions and backbone carbonyl oxygen atoms of gA. Hence, the protein matrix would provide an average of four carbonyls to stabilize the ions, which account for about ∼30−40 kcal/mol to compensate the ion dehydration barrier inside the channel. Moreover, the uneven total PMF profile during the ion permeation could also be attributed to the dynamic interaction changes with the carbonyls. During the ion permeating across the series of small free energy barriers which was predicted by Roux and Karplus,86 two carbonyls at one site are gradually replaced by the other two belonging to the successive free energy well (Figure 8), therefore, preventing the disruption of the hydration state of ions during permeation and prohibiting the hopping of ions inside the channel, which excellently agrees with the results reported by Roux and Karplus.38 The K+/Na+ Conductance Calculations Using the AMOEBA Polarizable Force Field. Using the GLE-HO method, the diffusion coefficient D(z) along the gA channel was calculated and then used to obtain the gmax through the 1D Nernst-Plank equation. In Figure 9, we show the diffusion coefficient of K+ and Na+ passing through the gA channel, respectively. The lower diffusion rate of Na+ at each local minimum explains the selectivity of gA of K+ over Na+. At each maximum, the diffusion coefficients of K+ and Na+ are equivalent to each other and just slightly lower than that in the bulk region. Due to the accurate description of the permanent multipole moments and polarization effects in AMOEBA, the K+/Na+ ions are identified to transit faster between two successive binding sites than the fixed point charge models.80 The calculated ion conductance for K+ and Na+ using AMOEBA are listed in Table 1 and compared to those obtained from experiment observations27,87,88 and classical force fields such as CHARMM,25−27 AMBER,26 and GROMOS.26 The maximum conductance calculated using AMOBEA are 36.8 pS for K+ and 14.7 pS for Na+, which are closer to the experiment than any previous classical MD calculations without any corrections.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00128. Figures S1−S4 and AMOEBA based polarizable force field data for the DMPC molecule (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail:[email protected]. Author Contributions

X.P., Y.Z., and H.C. contributed equally to this work. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS



REFERENCES

This work is supported by grants from the National Natural Science Foundation of China under Contracts No. 21573217, No. 31370714, and No. 91430110.

(1) Lindahl, E.; Sansom, M. S. P. Membrane proteins: molecular dynamics simulations. Curr. Opin. Struct. Biol. 2008, 18 (4), 425−431. (2) Karplus, M.; McCammon, J. A. Molecular dynamics simulations of biomolecules. Nat. Struct. Biol. 2002, 9 (9), 646−652. (3) Biggin, P. C.; Bond, P. J. Molecular dynamics simulations of membrane proteins. Methods Mol. Biol. (N. Y., NY, U. S.) 2008, 443, 147−60. (4) Stansfeld, P. J.; Sansom, M. S. P. Molecular Simulation Approaches to Membrane Proteins. Structure 2011, 19 (11), 1562− 1572. (5) Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. In Annual Review of Biophysics; Rees, D. C., Ed.; Palo Alto, 2012; Vol. 41, pp 429−452.10.1146/annurev-biophys-042910-155245 (6) Maffeo, C.; Bhattacharya, S.; Yoo, J.; Wells, D.; Aksimentiev, A. Modeling and Simulation of Ion Channels. Chem. Rev. 2012, 112 (12), 6250−6284. (7) Kelkar, D. A.; Chattopadhyay, A. The gramicidin ion channel: A model membrane protein. Biochim. Biophys. Acta, Biomembr. 2007, 1768 (9), 2011−2025. (8) Roux, B. Computational studies of the gramicidin channel. Acc. Chem. Res. 2002, 35 (6), 366−375. (9) Sigg, D. Modeling ion channels: Past, present, and future. J. Gen. Physiol. 2014, 144 (1), 7−26. (10) Pielak, R. M.; Chou, J. J. Influenza M2 proton channels. Biochim. Biophys. Acta, Biomembr. 2011, 1808 (2), 522−529. (11) Timko, J.; Kuyucak, S. Investigation of polarization effects in the gramicidin A channel from ab initio molecular dynamics simulations. J. Chem. Phys. 2012, 137 (20), 205106.



CONCLUSION The free energy profiles of K+/Na+ permeating through the gA channel were accurately estimated using the AMOEBA polarizable force field. The calculated potential of mean force (PMFs) profiles of both K+/Na+ exhibits a significant improvement over the conventional force field with the fixedcharge electrostatic model. In the case of K+ permeating H

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation (12) Ren, P. Y.; Ponder, J. W. Polarizable atomic multipole water model for molecular mechanics simulation. J. Phys. Chem. B 2003, 107 (24), 5933−5947. (13) Ren, P.; Wu, C.; Ponder, J. W. Polarizable Atomic MultipoleBased Molecular Mechanics for Organic Molecules. J. Chem. Theory Comput. 2011, 7 (10), 3143−3161. (14) Shi, Y.; Xia, Z.; Zhang, J. J.; Best, R.; Wu, C. J.; Ponder, J. W.; Ren, P. Y. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9 (9), 4046−4063. (15) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A., Jr.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; HeadGordon, T. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114 (8), 2549−2564. (16) Grossfield, A.; Ren, P. Y.; Ponder, J. W. Ion solvation thermodynamics from simulation with a polarizable force field. J. Am. Chem. Soc. 2003, 125 (50), 15671−15682. (17) Piquemal, J.-P.; Perera, L.; Cisneros, G. A.; Ren, P.; Pedersen, L. G.; Darden, T. A. Towards accurate solvation dynamics of divalent cations in water using the polarizable amoeba force field: From energetics to structure. J. Chem. Phys. 2006, 125 (5), 054511. (18) Wu, J. C.; Piquemal, J.-P.; Chaudret, R.; Reinhardt, P.; Ren, P. Polarizable Molecular Dynamics Simulation of Zn(II) in Water Using the AMOEBA Force Field. J. Chem. Theory Comput. 2010, 6 (7), 2059−2070. (19) Shi, Y.; Wu, C.; Ponder, J. W.; Ren, P. Multipole Electrostatics in Hydration Free Energy Calculations. J. Comput. Chem. 2011, 32 (5), 967−977. (20) Manzoni, F.; Soderhjelm, P. Prediction of hydration free energies for the SAMPL4 data set with the AMOEBA polarizable force field. J. Comput.-Aided Mol. Des. 2014, 28 (3), 235−244. (21) Shi, Y.; Jiao, D. A.; Schnieders, M. J.; Ren, P. Y. In TrypsinLigand Binding Free Energy Calculation with AMOEBA, Annual International Conference of the Ieee Engineering in Medicine and Biology Society, Vols 1−20, Minneapolis, MN, Sep 03−06, 2009; IEEE: Minneapolis, MN, 2009; pp 2328−2331. (22) Satpati, P.; Clavaguera, C.; Ohanessian, G.; Simonson, T. Free Energy Simulations of a GTPase: GTP and GDP Binding to Archaeal Initiation Factor 2. J. Phys. Chem. B 2011, 115 (20), 6749−6763. (23) Chu, H.; Peng, X.; Li, Y.; Zhang, Y.; Zhang, C.; Ren, P.; Li, G. Polarizable Atomic Multipole-Based Force Field for POPC Membrane Lipids. J. Comput. Chem. 2016, submitted for publication. (24) Allen, T. W.; Bastug, T.; Kuyucak, S.; Chung, S. H. Gramicidin A channel as a test ground for molecular dynamics force fields. Biophys. J. 2003, 84 (4), 2159−2168. (25) Allen, T. W.; Andersen, O. S.; Roux, B. Energetics of ion conduction through the gramicidin channel. Proc. Natl. Acad. Sci. U. S. A. 2004, 101 (1), 117−122. (26) Allen, T. W.; Andersen, O. S.; Roux, B. Ion permeation through a narrow channel: Using gramicidin to ascertain all-atom molecular dynamics potential of mean force methodology and biomolecular force fields. Biophys. J. 2006, 90 (10), 3447−3468. (27) Ingolfsson, H. I.; Li, Y.; Vostrikov, V. V.; Gu, H.; Hinton, J. F.; Koeppe, R. E., II; Roux, B.; Andersen, O. S. Gramicidin A Backbone and Side Chain Dynamics Evaluated by Molecular Dynamics Simulations and Nuclear Magnetic Resonance Experiments. I: Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115 (22), 7417−7426. (28) Jensen, M. O.; Jogini, V.; Eastwood, M. P.; Shaw, D. E. Atomiclevel simulation of current-voltage relationships in single-file ion channels. J. Gen. Physiol. 2013, 141 (5), 619−632. (29) Song, H. D.; Beck, T. L. Temperature Dependence of Gramicidin Channel Transport and Structure. J. Phys. Chem. C 2013, 117 (8), 3701−3712. (30) Siu, S. W. I.; Boeckmann, R. A. Low Free Energy Barrier for Ion Permeation Through Double-Helical Gramicidin. J. Phys. Chem. B 2009, 113 (10), 3195−3202.

(31) Liu, Z. W.; Xu, Y.; Tang, P. Steered molecular dynamics simulations of Na+ permeation across the gramicidin a channel. J. Phys. Chem. B 2006, 110 (25), 12789−12795. (32) De Fabritiis, G.; Coveney, P. V.; Villa-Freixa, J. Energetics of K+ permeability through Gramicidin A by forward-reverse steered molecular dynamics. Proteins: Struct., Funct., Genet. 2008, 73 (1), 185−94. (33) Forney, M. W.; Janosi, L.; Kosztin, I. Calculating free-energy profiles in biomolecular systems from fast nonequilibrium processes. Phys. Rev. E. Stat. Nonlin. Soft. Matter. Phys. 2008, 78 (5 Pt 1), 051913. (34) Mustafa, M.; Busath, D. D. The gramicidin channel ion permeation free-energy profile: direct and indirect effects of CHARMM force field improvements. Interdiscip. Sci.: Comput. Life Sci. 2009, 1 (2), 113−27. (35) Vorobyov, I.; Bekker, B.; Allen, T. W. Electrostatics of deformable lipid membranes. Biophys. J. 2010, 98 (12), 2904−13. (36) Giorgino, T.; De Fabritiis, G. A High-Throughput Steered Molecular Dynamics Study on the Free Energy Profile of Ion Permeation through Gramicidin A. J. Chem. Theory Comput. 2011, 7 (6), 1943−50. (37) Olah, G. A.; Huang, H. W.; Liu, W. H.; Wu, Y. L. Location of ion-binding sites in the gramicidin channel by x-ray-diffraction. J. Mol. Biol. 1991, 218 (4), 847−858. (38) Roux, B.; Karplus, M. Ion-transport in a model gramicidin channel - structure and thermodynamics. Biophys. J. 1991, 59 (5), 961−981. (39) Patel, S.; Davis, J. E.; Bauer, B. A. Exploring Ion Permeation Energetics in Gramicidin A Using Polarizable Charge Equilibration Force Fields. J. Am. Chem. Soc. 2009, 131 (39), 13890. (40) Kumar, R.; Iyer, V. G.; Im, W. CHARMM-GUI: A graphical user interface for the CHARMM users. Abstr. Pap. Am. Chem. Soc. 2007, 233, 273−273. (41) Andersen, H. C. Molecular-dynamics simulations at constant pressure and-or temperature. J. Chem. Phys. 1980, 72 (4), 2384−2393. (42) Eastman, P.; Friedrichs, M. S.; Chodera, J. D.; Radmer, R. J.; Bruns, C. M.; Ku, J. P.; Beauchamp, K. A.; Lane, T. J.; Wang, L.-P.; Shukla, D.; Tye, T.; Houston, M.; Stich, T.; Klein, C.; Shirts, M. R.; Pande, V. S. OpenMM 4: A Reusable, Extensible, Hardware Independent Library for High Performance Molecular Simulation. J. Chem. Theory Comput. 2013, 9 (1), 461−469. (43) Andersen, H. C. RATTLE: A velocity version of the shake algorithm for molecular-dynamics calculations. J. Comput. Phys. 1983, 52 (1), 24−34. (44) 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 (19), 8577−8593. (45) Lagardere, L.; Lipparini, F.; Polack, E.; Stamm, B.; Cances, E.; Schnieders, M.; Ren, P. Y.; Maday, Y.; Piquemal, J. P. Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Toward Massively Parallel Computations Using Smooth Particle Mesh Ewald. J. Chem. Theory Comput. 2015, 11 (6), 2589−2599. (46) Shirts, M. R.; Mobley, D. L.; Chodera, J. D.; Pande, V. S. Accurate and efficient corrections for missing dispersion interactions in molecular Simulations. J. Phys. Chem. B 2007, 111 (45), 13052−13063. (47) Goings, J. J.; Ding, F.; Li, X. Self-Consistent Field using Direct Inversion in Iterative Subspace Method and Quasi-Newton Vectors. In Proceedings of Mest 2012: Electronic Structure Methods with Applications to Experimental Chemistry; Hoggan, P. E., Ed.; 2014; Vol. 68, pp 77− 86.10.1016/B978-0-12-800536-1.00004-6 (48) Ren, P.; Ponder, J. W. Tinker polarizable atomic multipole force field for proteins. Abstr. Pap. Am. Chem. Soc. 2002, 224, U473−U473. (49) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. I

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (50) Stone, A. J. Distributed multipole analysis: Stability for large basis sets. J. Chem. Theory Comput. 2005, 1 (6), 1128−1132. (51) Kastner, J. Umbrella sampling. Wiley Interdisciplinary ReviewsComputational Molecular Science 2011, 1 (6), 932−942. (52) Torrie, G. M.; Valleau, J. P. Non-physical sampling distributions in Monte-Carlo free-energy estimation - umbrella sampling. J. Comput. Phys. 1977, 23 (2), 187−199. (53) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. THE weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13 (8), 1011−1021. (54) Kastner, J.; Thiel, W. Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: ″umbrella integration″. J. Chem. Phys. 2005, 123 (14), 144104. (55) Shirts, M. R.; Bair, E.; Hooker, G.; Pande, V. S. Equilibrium free energies from nonequilibrium measurements using maximum-likelihood methods. Phys. Rev. Lett. 2003, 91 (14), 140601. (56) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 2008, 129 (12), 124105. (57) Lee, T. S.; Radak, B. K.; Pabis, A.; York, D. M. A New Maximum Likelihood Approach for Free Energy Profile Construction from Molecular Simulations. J. Chem. Theory Comput. 2013, 9 (1), 153− 164. (58) Bartels, C. Analyzing biased Monte Carlo and molecular dynamics simulations. Chem. Phys. Lett. 2000, 331 (5−6), 446−454. (59) Hansen, N.; Dolenc, J.; Knecht, M.; Riniker, S.; van Gunsteren, W. F. Assessment of enveloping distribution sampling to calculate relative free enthalpies of binding for eight netropsin-DNA duplex complexes in aqueous solution. J. Comput. Chem. 2012, 33 (6), 640− 651. (60) Maragakis, P.; van der Vaart, A.; Karplus, M. Gaussian-Mixture Umbrella Sampling. J. Phys. Chem. B 2009, 113 (14), 4664−4673. (61) Sprik, M.; Ciccotti, G. Free energy from constrained molecular dynamics. J. Chem. Phys. 1998, 109 (18), 7737−7744. (62) 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 (10), 1961−1972. (63) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185 (2), 604−613. (64) Peng, X.; Zhang, Y.; Chu, H.; Li, G. Free energy simulations with the AMOEBA polarizable force field and metadynamics on GPU platform. J. Comput. Chem. 2016, 37 (6), 614. (65) Hub, J. S.; de Groot, B. L.; van der Spoel, D. g_wham-A Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation Estimates. J. Chem. Theory Comput. 2010, 6 (12), 3713−3720. (66) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4 (3), 435−447. (67) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26 (16), 1701−1718. (68) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976.

(69) LyndenBell, R. M.; Rasaiah, J. C. Mobility and solvation of ions in channels. J. Chem. Phys. 1996, 105 (20), 9266−9280. (70) Allen, T. W.; Kuyucak, S.; Chung, S. H. Molecular dynamics estimates of ion diffusion in model hydrophobic and KcsA potassium channels. Biophys. Chem. 2000, 86 (1), 1−14. (71) Smith, G. R.; Sansom, M. S. P. Effective diffusion coefficients of K+ and Cl- ions in ion channel models. Biophys. Chem. 1999, 79 (2), 129−151. (72) Burykin, A.; Kato, M.; Warshel, A. Exploring the origin of the ion selectivity of the KcsA potassium channel. Proteins: Struct., Funct., Genet. 2003, 52 (3), 412−426. (73) Kubo, R. Many-Body Theory; Syokabo and Benjaming: Tokyo, 1966. (74) Koneshan, S.; Lynden-Bell, R. M.; Rasaiah, J. C. Friction coefficients of ions in aqueous solution at 25 degrees C. J. Am. Chem. Soc. 1998, 120 (46), 12041−12050. (75) Straub, J. E.; Borkovec, M.; Berne, B. J. Calculation of dynamic friction on intramolecular degrees of freedom. J. Phys. Chem. 1987, 91 (19), 4995−4998. (76) Crouzy, S.; Woolf, T. B.; Roux, B. A molecular-dynamics study of gating in dioxolane-linked gramicidin-a channels. Biophys. J. 1994, 67 (4), 1370−1386. (77) Mamonov, A. B.; Kurnikova, M. G.; Coalson, R. D. Diffusion constant of K+ inside Gramicidin A: A comparative study of four computational methods. Biophys. Chem. 2006, 124 (3), 268−278. (78) Levitt, D. G. Interpretation of biological ion channel flux data reaction-rate versus continuum theory. Annu. Rev. Biophys. Biophys. Chem. 1986, 15, 29−57. (79) Allen, T. W.; Andersen, O. S.; Roux, B. Molecular dynamics potential of mean force calculations as a tool for understanding ion permeation and selectivity in narrow channels. Biophys. Chem. 2006, 124 (3), 251−267. (80) Roux, B.; Allen, T.; Berneche, S.; Im, W. Theoretical and computational models of biological ion channels. Q. Rev. Biophys. 1999, 37 (1), 15−103. (81) Marrink, S. J.; Berendsen, H. J. C. Simulation of water transport through a lipid-membrane. J. Phys. Chem. 1994, 98 (15), 4155−4168. (82) Tepper, H. L.; Voth, G. A. Mechanisms of passive ion permeation through lipid bilayers: Insights from simulations. J. Phys. Chem. B 2006, 110 (42), 21327−21337. (83) Bastug, T.; Patra, S. M.; Kuyucak, S. Molecular dynamics simulations of gramicidin A in a lipid bilayer: From structure-function relations to force fields. Chem. Phys. Lipids 2006, 141 (1−2), 197−204. (84) Kato, M.; Warshel, A. Through the channel and around the channel: Validating and comparing microscopic approaches for the evaluation of free energy profiles for ion penetration through ion channels. J. Phys. Chem. B 2005, 109 (41), 19516−19522. (85) Tian, F.; Wang, J. F.; Cross, T. A. Structural basis for the gramicidin a channel specificity probed by solid state NMR. Biophys. J. 1998, 74 (2), A233−A233. (86) Roux, B.; Karplus, M. Ion Transport in the Gramicidin Channel: Free Energy of the Solvated Right-Handed Dimer in a Model Membrane. J. Am. Chem. Soc. 1993, 115 (8), 3250−3262. (87) Busath, D. D.; Thulin, C. D.; Hendershot, R. W.; Phillips, L. R.; Maughan, P.; Cole, C. D.; Bingham, N. C.; Morrison, S.; Baird, L. C.; Hendershot, R. J.; Cotten, M.; Cross, T. A. Noncontact dipole effects on channel permeation. I. Experiments with (5F-Indole)Trp(13) gramicidin A channels. Biophys. J. 1998, 75 (6), 2830−2844. (88) Finkelstein, A.; Andersen, O. S. The gramicidin-a channel: A review of its permeability characteristics with special reference to the single-file aspect of transport. J. Membr. Biol. 1981, 59 (3), 155−171. (89) Halgren, T. A. The Representation of van der Waals (vdW) Interactions in Molecular Mechanics Force Fields: Potential Form, Combination Rules, and vdW Parameters. J. Am. Chem. Soc. 1992, 114 (20), 7827−7843.

J

DOI: 10.1021/acs.jctc.6b00128 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX