Atomistic Simulations of Membrane Ion Channel Conduction, Gating

Review pubs.acs.org/CR

Cite This: Chem. Rev. 2019, 119, 7737−7832

Atomistic Simulations of Membrane Ion Channel Conduction, Gating, and Modulation Emelie Flood,† Ceĺ ine Boiteux,† Bogdan Lev,† Igor Vorobyov,‡ and Toby W. Allen*,† †

School of Science, RMIT University, Melbourne, Victoria 3000, Australia Department of Physiology & Membrane Biology/Department of Pharmacology, University of California, Davis, 95616, United States

Downloaded via GUILFORD COLG on July 17, 2019 at 18:17:57 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: Membrane ion channels are the fundamental electrical components in the nervous system. Recent developments in X-ray crystallography and cryo-EM microscopy have revealed what these proteins look like in atomic detail but do not tell us how they function. Molecular dynamics simulations have progressed to the point that we can now simulate realistic molecular assemblies to produce quantitative calculations of the thermodynamic and kinetic quantities that control function. In this review, we summarize the state of atomistic simulation methods for ion channels to understand their conduction, activation, and drug modulation mechanisms. We are at a crossroads in atomistic simulation, where long time scale observation can provide unbiased exploration of mechanisms, supplemented by biased free energy methodologies. We illustrate the use of these approaches to describe ion conduction and selectivity in voltage-gated sodium and acid-sensing ion channels. Studies of channel gating present a significant challenge, as activation occurs on longer time scales. Enhanced sampling approaches can ensure convergence on minimum free energy pathways for activation, as illustrated here for pentameric ligand-gated ion channels that are principal to nervous system function and the actions of general anesthetics. We also examine recent studies of local anesthetic and antiepileptic drug binding to a sodium channel, revealing sites and pathways that may offer new targets for drug development. Modern simulations thus offer a range of molecular-level insights into ion channel function and modulation as a learning platform for mechanistic discovery and drug development.

CONTENTS 1. Introduction 1.1. Hidden Kinetic Mechanisms of Ion Permeation and Selectivity 1.2. Ion Channel Activation Involves Complex Molecular Processes with Wide-Ranging Time Scales 1.3. Exploration of Ion Channel Drug Binding and Modulation 1.4. Explorative and Quantitative Simulation Approaches to Study Ion Channel Function 2. Molecular Mechanical Models of Membrane Ion Channels 2.1. Molecular Dynamics Simulation for Statistical Ensemble Sampling 2.1.1. The Potential of Mean Force 2.1.2. Thermodynamic Selectivity and Free Energy Perturbation 2.1.3. Enhanced Sampling Approaches 2.2. Development of Empirical Force Fields for Ion Channels 2.2.1. The Atomistic MD Force Field for Membrane Proteins 2.2.2. Developments in Polarizable Force Fields for Membrane Proteins

© 2019 American Chemical Society

2.2.3. Improved Ion Models for Accurate Channel Interactions 2.2.4. Ligand Models, Including Drug Interactions with Membranes and Ion Channels 3. Atomistic Models of Ion Channel Permeation and Selectivity 3.1. Membrane Ion Permeation Controlled by the Potential of Mean Force 3.1.1. Effective Dynamics of Ion Movements 3.1.2. Ion Channel PMF Computation 3.1.3. Making Contact with Experimental Ion Channel Conductance 3.1.4. Calculations to Explain Ion Selectivity 3.2. Case Study. Sodium Channel Conduction and Selectivity: from Bacteria to Human 3.2.1. Mammalian and Bacterial NaV Channels 3.2.2. Sodium-Selective Ion Permeation in Bacterial NaV Channels 3.2.3. Sodium-Selective Ion Permeation in Mammalian NaV Channels 3.2.4. Sodium-Selective Ion Permeation in Acid Sensing Ion Channels

7738 7738

7739 7740 7741 7742 7743 7745 7747 7749 7751 7751 7754

7755

7758 7760 7762 7764 7765 7768 7771 7773 7775 7777 7779 7783

Received: October 18, 2018 Published: June 27, 2019 7737

DOI: 10.1021/acs.chemrev.8b00630 Chem. Rev. 2019, 119, 7737−7832

Chemical Reviews 4. Capturing Protein Conformational Changes Associated with Channel Gating 4.1. Using Atomistic Simulations to Capture Activation Mechanisms 4.1.1. Exploring Conformational Change from a Single State 4.1.2. Approaches to Capture Interconversion of End States 4.1.3. The String Method Approach 4.2. Case Study: Signal Transduction in a Pentameric Ligand-Gated Ion Channel 4.2.1. String Solution for the GLIC Channel 4.2.2. Molecular Events Initiating Activation 4.2.3. PMFs Demonstrate Allosteric Change and Modulation 4.2.4. Communication Across the Gating Interface 4.2.5. Intermediates and the Pore Closing Process 5. Exploring Drug Binding to Ion Channels 5.1. Atomistic Methods for Simulated Docking of Ligands to an Ion Channel 5.1.1. Ligand Docking Approaches 5.1.2. Atomistic MD Simulations of Binding 5.2. Case Study: Local Anesthetic and Antiepileptic Binding to a Sodium Channel 5.2.1. Simulated Flooding to Explore Nav Drug Binding 5.2.2. Biased Sampling Approaches to Study NaV Drug Binding 6. Conclusions and Future Perspectives Associated Content Special Issue Paper Author Information Corresponding Author ORCID Notes Biographies Acknowledgments Abbreviations Used References

Review

7784 7785 7785 7786 7787 7790 7790 7792 7792 7793 7794 7795 7796 7796 7797 7800 7804 7806 7807 7808 7808 7808 7808 7808 7808 7808 7809 7809 7810

Figure 1. (A) Membrane ion channel showing conduction and selectivity. (B) Action potential generated by voltage-gated sodium (NaV) and potassium (KV) channels, with the rising phase caused by Na+ entry into the cell and the falling phase by K+ exit from the cell.

The pioneering studies of Hodgkin and Huxley5,6 started generations of scientists on a path to understanding the fundamental mechanisms of ion channel proteins that control vital processes in all cells and which are central to intelligent life. Ion channels now profit from decades of studies examining mechanisms of activation, ion conduction, and sensitivity to physiological and pharmacological conditions.1,7 The first atomic-resolution structure of an ion channel in 1998, for the bacterial K+ channel, KcsA,8 gave us an initial glimpse of the atoms responsible for channel function, with a collection of such structures now having emerged in the 20 years following this breakthrough.8−13 However, the structure by itself does not explain function. Proteins are highly flexible, and this flexibility not only influences but may determine protein activity. The challenge today is to make use of these new atomic-level structures using computational methods that can observe and quantify molecular mechanisms to make testable predictions of function.

1. INTRODUCTION Membrane ion channels are pore-forming proteins present in the membranes of all excitable cells, controlling the flow of salt ions to generate and shape nerve impulses as well as the regulation of chemical activity and cell volume. These proteins selectively transport ions across cell membranes (Figure 1A) in response to physiological stimuli, including neurotransmitter binding, pH, ions, stress, temperature, and voltage. For example, voltage-gated ion channels (VGICs) are the proteins responsible for propagation of nerve impulses, enabling vital functions like heartbeat, sensation, muscle contraction, and brain activity,1 making them key targets for drug development.2,3 The crucial role of ion channels in the nervous system has been known for the better part of a century,4 but it is only recently that we have begun to even see what they look like. Advances in X-ray crystallography and cryogenic electron microscopy (cryo-EM), highlighted by the 2003 and 2017 Nobel Prizes in Chemistry, have led to several high-resolution ion channel structures, providing an opportunity to understand ion channel function at the molecular level.

1.1. Hidden Kinetic Mechanisms of Ion Permeation and Selectivity

Permeation through a channel pore is an electrodiffusive process down an electrochemical gradient. Specific channels are able to provide the right amino acid functionalities to distinguish different species (providing higher permeability to one species than others), based on valence and size, with high fidelity.1 For example, the firing of neurons via the action potential (Figure 1B) involves voltage-gated sodium (NaV) and potassium (KV) channels responding to precise changes in 7738


Chemical Reviews

Review

Figure 2. (A) Example of a K+ channel pore (KcsA; 2 of 4 subunits shown for clarity, gray ribbons), highlighting the narrow SF (green oval in stick representation) with possible K+ ion positions (S0(top)−S4(bottom) and cavity, shown as blue balls) and the lower hydrophobic gate formed by the crossing of inner transmembrane (TM) helices (gray ribbons). Adapted with permission from ref 14. Copyright 2014. Rockefeller Press. Comparison of soft (B) and hard (C) multi-ion mechanisms of K+ ions (blue balls) permeating through a K+ channel SF (blue/red sticks). Water molecules (red/white sticks) intercalate ions in the soft but not the hard knock-on movements. Reproduced with permission from ref 23. Copyright 2014. American Association for the Advancement of Science.

membrane potential to release stores of Na+ and K+ ions in succession. Given the abundance of Na+ and K+ ions, without both of these channels being selective for their respective ion, their driving forces would become indistinguishable, eliminating nerve impulses. The ability of a protein pore to discern between Na+ and K+ ions has remained an intriguing problem for decades, involving two nearly identical metal ions that differ in radius by rc) that is only updated at regular intervals or heuristically when atomic movements are large. Today, faster grid-based neighbor lists are often used, where the simulation cell is divided up into cells with sides larger than the cutoff distance, such that all neighbors within the cutoff distance are in the same or neighboring cells.84 To overcome the problem of boundary artifacts for a finite system, it can be simulated with implicit solvent surroundings.60 This approach can, however, still have boundary effects on solvent structure that may have impact on the protein’s conformational sampling. One can increase the explicit representation of the bulk solvent or move to new methods that can adjust the system size on the fly.96 The more commonly used approach, however, is to impose periodic boundary conditions (PBC), where as a particle leaves the central image it is translated and enters from the opposite face.84 However, it remains important to simulate a large enough box of the system to avoid artifacts from the boundaries. For instance, it has recently been shown that extending the amount of solvent around the protein 10 times can lead to correct folding.97 For ion channels, the membrane and solvent layer dimensions can play important roles if one aims to compute accurate energetics and/or kinetics of ion permeation or channel conformational changes. One needs to avoid or correct for significant interactions between ion channel images and the effects of those images on ion translocation.98 Electrostatic interactions, generated by a charge distribution, ρ, according to Poisson’s equation, are conditionally convergent in this periodic system and must be treated with Ewald summation,85 where the long-ranged electrostatic contribution is solved (without truncation) via the Fourier transform (although noting that electrostatics still involves switching of a real-space electrostatic component). An efficient variant of this method, widely employed in biomolecular MD simulations, is particle mesh Ewald (PME),99 where the fast Fourier transform is used to compute charge density on a discrete lattice in reciprocal space. Variants of this methodology have been developed over the years, including the Gaussian split Ewald100 approach used on the Anton/Anton 2 supercomputers.55 Recently, PME was extended to Lennard-Jones (LJ) interactions in the method of PME-LJ,921 which better captures long-ranged interactions with minimal computational overhead. This approach has been shown to yield improved properties of liquid alkanes922 and may be an important development in future MD programs. With this basic framework for representing the membrane ion channel system and solving its equations of motion, one

modeling of protein function still provides the best route to sampling protein structure−function relationships, incorporating realistic descriptions of protein interactions with lipids, solvent, ions, and ligands, while covering much of the necessary length and time-scales for function, and enabling calculations of thermodynamic and kinetic measurables.83 Here we describe the current state of empirical MD simulation methods for statistical mechanical calculations and recent developments in MD force fields related to ion channel function. The purpose of this section is not to give a complete treatise of MD simulation methods, but to briefly illustrate its use as a tool for exploring configurations of the membrane ion channel system and computing free energies relevant to ion conduction, protein conformational change, and ligand binding. We then return to describe extensions for specific problems in the subsequent sections. 2.1. Molecular Dynamics Simulation for Statistical Ensemble Sampling

MD simulations can serve several purposes. The first is to follow the time evolution of structure and to visualize the interactions of molecules over time, the second is to sample equilibrium configurations from phase space for the calculation of ensemble averages related to experimental observables, and the third is to force the system to follow a process of interest that may not otherwise occur in an unbiased simulation time frame. We first discuss the basic MD method, as it applies to membrane proteins, although several detailed texts on MD simulation can be found elsewhere.57,84−93 The sole purpose of this section is to establish MD as an engine for generating accurate representations of the equilibrium configurations of the all-atom system to enable statistical averages and free energy estimation pertaining to studies of ion channel function. MD involves the numerical solution of Newton’s equation of motion F = ma, where the forces acting on every atom in the system are calculated based on the potential energy function U defined by the particular force field. This can be formulated as a function of the position vectors, ri, of N atoms i = 1,N, or alternatively in terms of functions of a 3N-dimensional vector of coordinates, q, as F(q1 , q2 , ... q3N ) = −∇U (q1 , q2 , ... q3N ).

(1)

Evaluation on the computer takes the form of numerical solutions of this second-order differential equation in a discretized time coordinate, using algorithms based on Taylor expansions. An example is the Verlet algorithm, which expresses positions at the next time step q(t + δt) in terms of the previous and current positions q(t − δt) and q(t), with velocities, v(t) estimated via finite differences.94,95 The time step, δt, is determined by the fastest motions of the system, such as bond vibrations involving hydrogen atoms, being of the order of 10 fs, requiring a time step of the order 1−2 fs to ensure accurate integration of the equations of motions. The number of time steps to carry out is set by the desired observational window to describe configurational changes related to structure and function, which may occur on the nanosecond, microsecond, or longer scales. For example, if one is interested in studying ion motions across an ion channel, this typically occurs anywhere in the range of ns−μs (or longer), and one must sample a full range of structural fluctuations including bonds, angles, backbone, and side chain torsions spanning fs−μs (Figure 5). The difficulty is that, to simulate 1 μs, one must calculate a billion iterations of the equations of 7743


Chemical Reviews

Review

may use MD simulation to sample the evolution of the system over time. If done at equilibrium, this allows for evaluation of configurational integrals within a chosen statistical ensemble that can allow for analysis of averages and thermodynamic quantities, including free energies that govern the process. In the absence of any treatment of the interaction with the surroundings, the equations of motion used in MD will produce a microcanonical (NVE: fixed number of particles N, volume V, and energy E) ensemble, with the probability of the system existing in a microstate i, Pi, being uniformly distributed, according to the postulate of equal a priori probabilities. However, this situation is not directly relevant to most experiments, although may approximate those results in the macroscopic, large N limit. We must instead use a computational approach to sample from canonical (NVT: fixed temperature T) or isothermal−isobaric (NPT: fixed pressure P) ensemble distributions.101 The probability of the system existing in a microstate i, Pi, with energy Ei, for the canonical distribution, where the system is in thermal equilibrium with a heat bath at temperature T, is given by a normalized Boltzmann factor to the state energy Pi(N , V , T ) =

exp( −Ei(N , V )/kBT ) , Z(N , V , T )

with variance 2

⟨ΔX ⟩ =

∑ exp(−Ei(N , V )/kBT ), i

(2)

⟨X ⟩ =

1 h3N

(3)

(4)

∫ dp3N dq3N e−H(p;q)/k T = h13N PNQ N , B

(5)

where H is the system Hamiltonian, as a function of all 3N positions and 3N momenta, and h is Planck’s constant, in the quantum regime. PN is the N-particle kinetic energy integral, PN =

∫ dp3N e−K(p)/k T ,

Pi =

(6)

N

τ

X (t ) d t =

1 τ ∑ Xi , Nτ i = 1

(11)

e−(Ei + PV )/ kBT , Δ(N , p , T )

(12)

where the isothermal−isobaric partition function is

∫ dq

3N −U (q)/ kBT

e

.

Δ(N , p , T ) = =

(7)

∫ dp3N dq3N X(p; q) e−H(p; q)/ kBT ⟨X ⟩ = , ∫ dp3N dq3N e−H(p; q)/ kBT

QN

e−(Ei + pV )/ kBT

i

∑ Z(N , V , T ) e−pV /k T . B

(13)

V

The natural thermodynamic function for variables N, p, and T is the Gibbs free energy, which is simply related to its partition function G(N , p , T ) = −kBT ln Δ(N , p , T ).

(8)

(14)

To maintain constant pressure, an extended system approach103 can again be used, coupling to a piston via a coordinate and momentum. The mass of the piston defines the decay time for volume fluctuations for the barostat,88 which can be applied isotropically or anisotropically.106−108 The Langevin piston method alternatively uses the Langevin

can be simplified if X is a function of only positional variables,

∫ dq3N X(q) e−U(q)/ kBT

∑∑ V

This is the quantity we compute when we run MD simulations to explore all possible 3N-dimensional vectors q from a canonical ensemble. The corresponding average of a quantity X

⟨X ⟩ =

∫0

B

which evaluates as a product of 3N Gaussian integrals, and QN is the N-particle configurational integral QN =

1 τ

over an interval τ, represented by Nτ discrete steps. To maintain the temperature and ensure a canonical distribution of states, one may employ a Berendsen approach,102 where the system is coupled to an external heat reservoir with weak coupling parameter, although this does not strictly generate a canonical ensemble.85,86 One may also couple to an external heat bath via Langevin dynamics,103 although influencing dynamics of the system.85 Today, extended Lagrangian methods, introduced by Nosé104 and Hoover,105 are most commonly employed. Here the heat reservoir is represented implicitly by an additional degree of freedom with coordinate and momentum in the equations of motion.85 Real-world laboratory experiments are maintained at constant pressure, p, not volume, where we seek the isothermal−isobaric distribution,

is the Helmholtz free energy, and is referred to as the characteristic function for the canonical ensemble, being the thermodynamic potential with natural variables N, V, and T. The partition function can be written as Z=

.

Any scalar, vector, or tensor quantity may be used for statistical analysis, with averaging performed according to eqs 9 and 10 above. We may extend this to calculate means and fluctuations involving any time-dependent variable, including structural (e.g., distribution functions, order parameters), dynamical (e.g., translational and rotational diffusion), or mechanical (e.g., elastic moduli),84−86,88 for example. The mean, according to the Ergodic Hypothesis, is postulated to equal the equivalent time average for that system at equilibrium. This assumes equilibration and sufficiently long production MD to sample the distributions of the variables in order to be ergodic.85 Because configurations are being sampled directly from a canonical ensemble (if N, V, and T are maintained; see below), then the configurations summed in this integral have already been weighted according to the Boltzmann factor appearing in the integrand. This means an average, such as in eq 9 above, becomes simply

from which all thermodynamic state functions can be calculated. For example (omitting independent variables), A = −kBT ln Z

B

Q N2

QN

2

(∫ dq3N X(q) e−U(q)/k T ) −

(10)

where kB is the Boltzmann constant and the denominator defines the canonical partition f unction Z (N , V , T ) =

∫ dq3N X2(q) e−U(q)/ kBT

, (9) 7744


Chemical Reviews

Review

equation (to be described below) to maintain pressure, eliminating unphysical volume fluctuations.109 Such approaches maintain temperature and pressure in ways that the distribution of microstates matches the desired ensemble for sampling configurations. We wish to make use of such ensemble configurational sampling using unbiased and biased MD approaches that allow for accurate calculation of thermodynamics and comparison to experiments for ion channel function. In a microsecond-long simulation, there may be millions of snapshots of the coordinates saved from which analysis of structure, dynamics, thermodynamics, and kinetics can be carried out. For a comprehensive discussion on various analysis methods for a range of properties, we refer to any one of the many simulations textbooks available (e.g., refs 84−86,89). Here we focus on free energy calculation, which is very important in computational biology110−113 and for ion channels is essential for defining the distributions and movements of ions, conformational states, and the binding of drugs, as highlighted in sections 3, 4, and 5. We briefly discuss the standard integrals required for free energies related to the movement and alchemical transformations of molecules related to ion channel function. 2.1.1. The Potential of Mean Force. The free energy defines the reversible work required to move from one region to another in phase space, such as between two conformational states of a protein, or for two states with dissociated and bound ions or ligands, for example. We are interested in the equilibrium distribution, via the probability density ρ(z), of some coordinate z, which may be, for example, the zcomponent of one particle’s position vector, such as that of a permeating ion crossing the membrane. However, one may choose any collective variable (such as a center of mass of a subset of atoms, an angle, dihedral or other function of q), with all other 3N-1 coordinates taking on all possible values. It can be shown that the density and free energy obey a simple relationship, such that we may write a configurational free energy

We can compute an effective force that is the negative derivative of this free energy, which is equivalent to the mean force acting on the coordinate z, Feffective = −

∂

= kBT ∂z

∫ dq3N−1 e−U(z ,q ,q ... q 2 3

3N

)/ kBT

∫ dq3N − 1 e−U(z , q2 , q3 ... q3N )/ kBT

∫ dq3N − 1 e−U(z , q2 , q3 ... q3N )/ kBT

= kBT

=

=

∫ dq3N − 1 ∂∂z e−U(z , q2 , q3 ... q3N )/ kBT ∫ dq3N − 1 e−U(z , q2 , q3 ... q3N )/ kBT

∫ dq3N − 1(− ∂∂Uz ) e−U(z , q2 , q3 ... q3N )/ kBT ∫ dq3N − 1 e−U(z , q2 , q3 ... q3N )/ kBT ∫ dq3N − 1F(z) e−U(z , q2 , q3 ... q3N )/ kBT = ⟨F(z)⟩. ∫ dq3N − 1 e−U(z , q2 , q3 ... q3N )/ kBT

(17)

Integrating both sides, z

W (z) − W (z′) = −

∫z′

⟨F(ξ)⟩ dξ

(18)

such that the free energy profile is given by the potential of mean force (PMF). i.e., it is the amount of reversible work associated with moving along the coordinate z. Therefore, one can calculate the free energy by either computing the (−kBT) natural logarithm of the probability density (eq 16) or by measuring the (negative) mean force and integrating (eq 18). PMF calculations require simulating long enough such that we sample equilibrium distributions of all coordinates, both the chosen order parameter of interest (e.g., z above), and all other (presumed fast) coordinates. However, it can be difficult to

W (z) = − kBT ln ρ(z)

sample all values of a coordinate when it corresponds to a rate

∫ dq3N δ(q1 − z) e−U(q ,q ,q ,... q )/k T + C = − kBT ln ∫ dq3N − 1 e−U (z , q , q ,... q )/ k T + C , (15) = − kBT ln

1 2 3

2 3

3N

3N

B

limiting step and thus involves an energy barrier that is not

B

thermally accessible in a normal simulation time frame. In such cases, one may apply biasing potentials to aid sampling. With

where we have an expression for the reversible work, W(z), determined by an N-particle configurational integral. As the probability density ρ(z) can be normalized or not, the free energy W(z) can be offset arbitrarily by some constant C, i.e., we may be interested in plotting W(z) relative to z = 0, infinity, or to a global or local minimum, for instance. We can eliminate the constant altogether by expressing one position relative to a reference z′ (e.g., a point for a molecule existing far from the ion channel in aqueous solution): W (z) − W (z′) = − kBT ln

∂W (z) ∂ = kBT ln ∂z ∂z

US,61 a series of biasing or “window” potentials, ui, such as harmonic potentials, may be added to the system potential, U, in a set of Nw independent simulations/‘windows’, i. Because it is known how the coordinate, z, was biased, it is possible to remove this bias postsimulation. If the unbiased density, ρ(z), is given by

ρ (z ) ρ(z′)

ρ (z ) =

∫ dq3N − 1 e−U(z , q2 , q3,... q3N )/ kBT . = −kBT ln ∫ dq3N − 1 e−U(z′, q2 , q3,... q3N )/ kBT

∫ dq3N − 1 e−U(z , q2 , q3,... q3N )/ kBT ∫ dq3N e−U(z , q2 , q3,... q3N )/ kBT

(19)

then a simulation under the action of a window biasing potential ui(z) would have density ρbiased(z) given by

(16) 7745


Chemical Reviews ρ biased (z) =

Review

∫ dq3N − 1 e−[U(z , q2 , q3,... q3N ) + ui(z)] / kBT ∫ dq3N e−[U(q1, q2 , q3,... q3N ) + ui(z)] / kBT

l 3N − 1 −[U (z , q2 , q3 ,... q3N ) + ui(z)] / kBT | o o e o o ∫ dq o =o m } o o N − U q q q q k T 3 ( , , ,... )/ o o B 1 2 3 3N o o dq e ∫ n ~ l o o ∫ dq3N e−[U (q1, q2 , q3 ,... q3N ) + ui(z)] / kBT | o o o o m } o o o o ∫ dq3N e−U (q1, q2 , q3 ,... q3N )/ kBT o o n ~ l − u z k T N − − U z q q q ( )/ 3 1 ( , , ,... )/ kBT | B i 2 3 3N o oe e ∫ dq o o o =o m } o o 3N −U (q1, q2 , q3 ,... q3N )/ kBT o o o o dq e ∫ n ~ −ui(z)/ kBT e = −u (z)/ k T ρ(z). B ⟩ ⟨e i

windows to solve for the optimal estimate for the unbiased density, ρ, and thus the PMF using eq 16 above. Here, a weight for each window is determined to minimize the overall variance.62 The unbiased probability distribution in eq 21 can be determined by a linear combination of the unbiased probability functions in each window, which can be shown to yield the WHAM equation N

ρ (z ) = ⟨e−ui(z)/ kBT ⟩

We may invert this relationship to obtain the desired unbiased density from knowledge of the biased distribution and the average of the Boltzmann factor of the biasing potential (21)

from which the PMF can be computed using eq 16, W (z) = −kBT ln ρi biased (z) − ui(z) + Fi + C .

(22)

Here Fi is a free energy constant, which represents the free energy cost of introducing the biasing potential to the system Fi = −kBT ln⟨e−ui(z)/ kBT ⟩ = −kBT ln

∫ ρ(z) e−u (z)/k T dz. i

B

(23)

In the case of multiple simulations with different biasing potentials, used to ensure sampling across an entire reaction coordinate of interest, potentials should be chosen carefully so that sufficient overlap between windows is obtained. If the free energy profile would be known ahead of time, the ideal ui(z) would be −W(z), which would yield a flat surface and good sampling. However, this is usually not the case and a common choice is instead to use a set of harmonic biasing potentials, 1 ui(z) = 2 ki(z − zi)2 , where ki is the force constant and zi is the center of window i.88 In this case, ki should be chosen to allow for overlapping histograms from neighboring windows. For a near-flat underlying PMF, the harmonic bias would yield an approximately Gaussian distribution with standard deviation, σ, of magnitude kBT /ki , based on the equivalence of the 2

Boltzmann factor e−ui(z)/ kBT = e−ki(z − zi) 2

−(z − z ̅ ) /2σ

/2kBT

N

∑ j =w1 nj e−(uj(z) − Fj)/ kBT

, (24)

where ni is the number of data points in window i. Because Fi depends on ρ(z) according to eq 23, these equations are solved iteratively until they are self-consistent, or meeting a selfconsistency criterion, and the optimal probability distribution is found.115 Self-consistency is generally determined by evaluating the change in free energy constants, Fi, between iterations. However, a small change in the free energy constants does not guarantee that the PMF has converged, and it is generally best to test convergence of the PMF W(z) itself, with a small tolerance (such as 0.0001−0.001 kcal/mol) between many iterations (such as testing every 100 iterations).116 Uncertainties in the computed PMFs can be computed by various means, including simple block analysis with calculation of standard error of means or through various bootstrapping approaches, based on the umbrella histograms with consideration given to the autocorrelation time,117 or Bayesian bootstrapping,118 for example (see also discussion on multistate Bennett acceptance ratio (MBAR) below, which includes estimates of statistical uncertainty119). Extensions of multistate analysis methods, such as the variational free energy profile (vFEP) exist to improve on WHAM for optimal calculation of free energy profiles.120 One may alternatively use the independent window simulations to calculate the mean force by creating a time series of instantaneous force acting on the coordinate z, or simply use the magnitude of the constraint force (equal and opposite to the instantaneous force at equilibrium) from the displacement and Hooke’s law, and then obtain the PMF via integration

(20)

ρ(z) = eui(z)/ kBT ⟨e−ui(z)/ kBT ⟩ρ biased (z),

∑i =w1 ni ρi biased (z)

z

W (z) − W (z′) = +

∫z′

⟨Fcons(ξ)⟩ dξ.

(25)

The results will be equal in the limit of complete equilibrium sampling. In the case studies of sections 3−5, each of these approaches will be used. For long multimicrosecond simulations with good sampling of ion or drug coordinates, we will simply apply eq 16. For ion interactions with ligands or channel pores or for drug movements across membranes that experience significant barriers, we will use the US biased approach with WHAM analysis (eq 24). The alternative formulation via the mean force in eq 18 (or via constraint force in eq 25) will prove to be useful in decomposing reversible work to reveal the important interactions controlling ion movements (section 3) and protein conformational changes (section 4). In each case, however, the formalism is only valid if the remaining coordinates, implicit in the configurational integration of eq 16, are thoroughly sampled at equilibrium. The challenge is that some of these coordinates may be slowly varying relative to the simulation time frame. We return later in this section to briefly discuss some modern enhanced sampling approaches to help ensure this is achieved.

and a normal

2

distribution, e , or equally, using the result from equipartition theorem stating that the average of the constraint will equal 1/2kBT. Setting the window spacing and this standard deviation to be of similar size should ensure overlap in most cases. Methods exist to test the overlap of windows, such as, for example, the Eigenvector Method for Umbrella Sampling (EMUS) approach of Dinner and co-workers.114 However, large mean forces may lead to displaced distributions that could prevent overlap without tighter constraints, in which case ki can be increased or additional windows can be added in that region. Once sufficient sampling of those window distributions has been obtained, unbiasing and optimization is needed to yield the best combination of data from all windows. The Weighted Histogram Analysis Method (WHAM)62 can take the biased histograms from multiple 7746


Chemical Reviews

Review

2.1.2. Thermodynamic Selectivity and Free Energy Perturbation. Free energies may also be calculated via perturbations to the Hamiltonian using FEP.111,121 FEP methods allow for chemical as well as physical coordinate changes and may include alchemical transformations such as single site mutations,122−124 ion transformations to explore selectivity of binding or permeation,125−129 protonation states of titratable groups in proteins, 130−132 and chemical modification of drug molecules,133−135 for example. Consider the free energy difference between two states of the system that may differ by positions/collective variables (e.g., molecular shape, ligand binding), application of an external field, or an unphysical change (alchemical transformation). If we perturb the energy by an amount ΔU, leading to a change from potential U0 to U1, corresponding to a change from state 0 to state 1,

absolute binding free energy for an ion or molecule to an ion channel. As λ → 0, a singularity occurs due to the repulsive portion of the LJ potential112,142 that can be alleviated using nonlinear scaling.112,144 An efficient way to solve this is a staged protocol142 based on the Weeks, Chandler, and Andersen scheme to separate repulsive and attractive parts of the LJ potential.145 Here a linear scheme is used for the attractive LJ contribution, but the repulsive part is gradually transformed into a soft-core potential using a nonlinear staging parameter.142 Because any arbitrary number of steps can be taken to yield the same path independent result, we may, in theory, take infinitely n small λ steps, with n → ∞. For infinitesimally small steps in λ, we may write, to first order,

U1(q1, q2 , q3 ... q3N ) = U0(q1, q2 , q3 ... q3N ) + ΔU (q1, q2 , q3 ... q3N ),

then one can show that the corresponding free energy change is given by

ΔG =

ΔG = G1 − G0 = − kBT ln⟨e−ΔU / kBT ⟩0 ≡ + kBT ln⟨e+ΔU / kBT ⟩1 , (27)

ΔG = −

(28)

∫0

e−Fk / kBT =

L

−

Nw

(29)

Nt

∑∑ i

The sum of window free energy changes

dλ . (32)

λ

∂U (z) ∂z

dz = −

∫0

L

F (z ) d z ,

(33)

t

e−Uk(q i(t ))/ kBT N ∑ j w nj

e−(Uj(q i(t )) − Fj)/ kBT

, (34)

where t = 1,Nt time steps, and the energy for the kth window evaluated for the configuration from the ith window

n i − 1 → λi

∂Uλ ∂λ

which is equivalent to the PMF along z in eq 25. Multiple histogram reweighting methods exist to incorporate data from simulations of multiple windows to improve estimates of free energy differences. Estimates of the free energy of perturbation may be improved through a WHAM analysis based on perturbation energy histograms.146,147 Assuming linear switching, for a set of λk values for windows k = 1, Nw, from U0 to U1, the free energy constants Fk corresponding to window k are obtained by iteratively, converging

ΔG λi−1→ λi = G λi − G λi−1 = −kBT ln⟨e−ΔUλi−1→ λi / kBT ⟩λi−1 .

i=1

1

This is the method of thermodynamic integration, where one obtains the path independent free energy change by integrating the force along the chosen coordinate, λ, from state 0 to state 1. This coordinate can be unphysical or physical, as we see for the case where states 0 and 1 are connected simply by a positional coordinate, i.e., if we set λ proportional to coordinate z, ranging from 0 to L, where L can be, for instance, a length of a channel pore. Then

where, for example, λ increases in n steps of 1/n length (e.g., 10 steps of 0.1), called “windows”, such that the free energy change for a single window is given by

∑ ΔG λ

∫0

63

which is simply the −kBT ln of the average of the Boltzmann factor within the ensemble of state 0 (as indicated by the subscript). These equivalent formulae are referred to as FEP or thermodynamic perturbation.111,136,137 The two quantities ⟨e−β ΔU ⟩0 and ⟨e+β ΔU ⟩1−1 should, in theory, be equivalent. In practice they are not because they are sampled from a finite number of configurations for which the system is in state 0 or 1, i.e., the low energy configurations of state 1 may not be sampled in a simulation performed with energy U0, (just as configurations sampled with U1 may not include any low energy configurations for state 0). The calculation will overestimate the free energy cost of changing from state 0 to state 1. To improve configurational sampling between states 0 and 1, we may introduce a coupling parameter, λ, ranging from 0 to 1 in n steps. One may linearly interpolate

ΔG =

(31)

and one may show that this leads to

(26)

Uλ = U0 + λΔU = (1 − λ)U0 + λU1 ,

∂Uλ δλ , ∂λ

Uλi = Uλi−1 +

Uk(q i(t )) = (λk − 1)U0(q i(t )) + λk U1(q i(t )).

(30)

is independent of the intermediate states owing to the path independence of the free energy. If the magnitude of the perturbation is large, then the distribution of configurations at one value of λ and the next could be very different. To ensure good convergence, one ideally should ensure that entropic changes remain within a couple of kBT between windows.138 One can also use nonlinear scaling to avoid rapid changes.139 Difficulties can arise, for example, due to linear coupling for the vdW energy term upon particle annihilation (e.g., see refs 140−143), as is important for calculations of solvation and protein binding free energies; relevant to the calculation of an

(35)

This WHAM approach for perturbation energies is the standard approach used in programs such as CHARMM.148 However, a common analysis technique used by, for example, the Visual Molecular Dynamics (VMD) program,149 is to examine end point energy histograms via the methods of Bennett.150 Equation 27 showed how the free energy associated with a perturbation in energy can be written as a sum over microstates, or alternatively an integral over phase space, weighted according to the unperturbed canonical distribution. One may show that this may be rewritten as an integral over energy levels for the perturbation151 7747


Chemical Reviews

Review

Figure 6. (A) Perturbation energy histograms for the case of a strongly driven harmonic oscillator showing the free energy change at the intersection of end point histograms. Reproduced with permission from ref 152. Copyright 2012 AIP publishing. (B) Possible sites in the KcsA channel for K+ and Na+. (C) Changing binding site locations during the perturbation from K+ (cage of 8 carbonyls) to Na+ (plane of 4 carbonyls) in the KcsA channel (2/4 subunits shown for clarity). Reproduced with permission from ref 153. Copyright 2011 National Academy of Sciences USA.

e−ΔG / kBT =

∫

du e−u / kBT ρ0 (u) =

∫

which is known to introduce bias into calculations of free energy. Moreover, their approach provides direct measures of statistical uncertainty, not available through WHAM (see ref 119 for descriptions of this method). In some cases, poor overlap of histograms from end point simulations is unavoidable even in long simulations. This can also be the case when using thermodynamic perturbation with a number of intermediates, as the states λi−1 and λi can be separated by large barriers. One solution to this problem, when the slow coordinate can be identified, is to combine FEP simulations with US via different approaches.147,152 By ensuring exhaustive sampling of configurational space by using biased sampling along such a slow coordinate, correct equilibrium weights can be imposed on conditional histograms of energy difference to ensure strongly overlapping end states and accurate FEP calculations. This has been demonstrated for cases such as potassium ion channel selectivity, for which end point states of the transformation K+ → Na+ sample kinetically separated regions in phase space but are easily connected by an ion translocation coordinate, z (see Figure 6B,C), that can be mapped out using US.153 Figure 6B shows how K+ and Na+ ions exhibit distinct binding sites where K+ prefers a cage of eight carbonyl groups created by a tetramer of two neighboring amino acids, whereas the smaller Na+ ion prefers the lower coordination offered by a plane of four carbonyl groups created by a single tetramer of amino acids. These distinct coordination complexes are separated by only 1−2 Å, but this has been shown to create a barrier that is difficult to overcome during a K+ → Na+ transformation.153 Figure 6C shows how the free energy minimum for Na+ (λ = 1; near z = −1 Å) does not exist when simulating K+ (λ = 0), and vice versa, and it has been shown that for partial λ values a double well sampling problem appears, with two minima separated by a barrier that is difficult to surmount in typical FEP simulations. Knowing what the slow coordinate is, one can

du eu / kBT ρ1(u), (36)

where ρ0/1(u) defines the energy perturbation histogram obtained in the 0/1 state. According to Bennett,150 we may take the ratio of histograms for end point states and find ln

ρ0 (u) ρ1(u)

= (u − ΔG)/kBT .

(37)

This equation says that at the value of u for which the end point distributions, ρ0 and ρ1, match, i.e., ρ0(u)/ρ1(u) = 1, the free energy change ΔG = u. One therefore simply has to see where the distributions overlap and read off the value of the free energy from the graph, as illustrated in Figure 6A. This Bennett overlapping histogram (BOH) method150 can be useful for cases where end points are well sampled to ensure overlap between states but will be problematic when those states are kinetically separated on the simulation time frame. The free energy difference can also be optimized, according to the Bennett acceptance ratio (BAR) method,150 by a reweighting to minimize the variance in the free energy change, yielding the BAR equation ΔGiest − j − c = − kBT ln

∑ ⟨f ( −(U1 − U0) + c)⟩ , ∑ ⟨f ((U1 − U0) + c)⟩

(38)

where c = ΔG + kBT ln(N1/N0) and f(x) = 1/1(1 + ex) is the Fermi function. This equation is solved iteratively until a selfconsistence criterion is met.150 For two end point windows, BAR is equivalent to WHAM in the limit of large sampling.119 Shirts and Chodera have derived an improved estimator for free energies using multiple state simulations, referred to as multistate BAR (or MBAR).119 This formalism reduces to BAR in the limit of two states and, unlike WHAM, has the advantage of not requiring the binning of data into histograms, 7748


Chemical Reviews

Review

biases, each represented as a product of NCV repulsive Gaussian functions of heights hi and widths σi,

invoke biasing approaches to ensure that coordinate is sampled. In this case, US along the z-direction, spanning only a couple of Å in length, is sufficient to overcome the challenge and ensure accurate free energies.153 Despite the increased computational power available today, sampling limitations remain, where slow coordinates, sometimes unknown to the simulator, can prevent accurate free energy calculation. Many modern simulations make use of rigorous approaches that help overcome these problems through appropriate choices of biasing or replica techniques. Below, we give a short summary of the possible approaches, some of which have been used to sample ion permeation, protein conformational change and drug binding in sections 35. 2.1.3. Enhanced Sampling Approaches. Limitations in resources restrict time scales such as to prevent the ability to observe equilibrium distributions without enhanced techniques; e.g., in section 4, we will face the daunting challenge of simulating ion channel gating, where the sampling of protein conformational changes is complicated by wide-ranging time scales of the thermal fluctuations of the protein.154,155 Here there may potentially exist many coordinates that experience large energetic barriers, with the time scales for crossings comparable to or exceeding simulation times. To achieve improved sampling, biases can be introduced into the system to ensure those slow coordinates are mapped out. An array of enhanced sampling methods, both equilibrium and nonequilibrium, have been developed, as reviewed elsewhere (e.g., refs 90,156). One may employ non-Boltzmann sampling approaches, designed to boost sampling of elevated regions of the underlying free energy surface in a defined configurational space. These approaches require some knowledge about the system and what determines the slowest motions. There is always the risk that hidden variables, which do not necessarily make up the largest movements, may represent the slowest changes and act orthogonally to (or be uncorrelated with) other variables that are controlled. Without inclusion of all slow variables (relative to the simulation time frame), the resulting distributions and free energetics will correspond to a restricted region of phase space, away from thermal equilibrium. Configurational transitions observed in that case may also not be representative of the most probable pathways for ion, drug, or protein movements. This can lead to free energy overestimation, activation barrier underestimation, hysteresis, and potentially irrelevant mechanisms.157 In the simplest form, methods exist that predetermine relevant coordinates and biasing potentials to guide the sampling of changes, such as with US simulations61 above. There are a range of other approaches that also work in a predefined space with set bias potentials to achieve increased sampling,158−161 as well as methods that allow the simulation to learn from its history, via adaptive biasing,162−166 changing their biasing potentials on the fly to ensure the system continues to visit new regions of the space. These methods notably include MTD64 and its variants, such as welltempered65 and multiwalker167 MTD, adaptive biasing force,168 adiabatic bias MD,169 conformational flooding,170 and related approaches.171,172 For example, MTD aims to enhance sampling by discouraging revisiting of previously sampled regions of configurational space through the addition of biasing potentials, ubias, as a function of NCV chosen collective variables (CV; qi), at regular intervals, τ, during MD simulation. The bias is of the form of a sum of the history of

ij ((q (t ) − q (t ′))2 yz zz. i i zz z 2σi2 k {

∑ ∏ hi expjjjjj− t

ubias(q(t ), t ) =

NCV

t ′= 0 i = 1

(39)

When equilibrated, the free energy minima will be flooded with an accumulation of repulsive gaussians, leading to a flat free energy surface allowing free diffusion across the space, with the underlying (unbiased) free energy surface equal to the negative of the total cumulative biasing potential.64 Methods exist to improve convergence, such as well-tempered MTD,65 where the biasing potential is adapted based on the history to speed up early barrier crossing events. As the expected limiting factor in sampling configurational changes is the roughness of the underlying free energy surface, and thus the ability to overcome activation barriers with available thermal energy, the most obvious approach to improve sampling is to raise the temperature. However, results at a much-elevated temperature may not be relevant to those at the physiological temperature due to thermal effects on structure and the entropy. One can use replica exchange (REX),173 or parallel tempering, to simulate replicas of the system at multiple temperatures and regularly exchange configurations between different replicas (i and j, with energies Ui and Uj, and temperatures Ti and Tj, respectively) with Monte Carlo, with probability ij e−Uj / kBTi − Ui / kBTj yz p = minjjj1, −U / k T − U / k T zzz j e i B i j B jz k {

= min(1, e−(Ui − Uj)(1/ kBTi − 1/ kBTj)).

(40)

This way, the low temperature simulation benefits from improved sampling at higher temperatures, requiring little knowledge of the system itself. REX lends itself to common multinode computer clusters because the replicas are run independently.174 However, these methods experience the demand of simulating either enough replicas, or enough time, to ensure high exchange counts via a Metropolis criterion,175 where for a large system the energy change associated with exchange can be high. As the number of degrees of freedom increases, the distribution of energies becomes narrow, reducing transition probabilities between replicas. The number of systems required grows inversely to this narrowing of distribution (as the square root of the number of particles), representing a significant burden, and a waste, given most particles represent solvent. The limitation in REX is therefore that the entire Hamiltonian of the system is scaled by a changing temperature, reducing transition probabilities. To limit the amount of energy change associated with exchange of configurations between replicas, one can instead just change individual terms in the Hamiltonian (rather than scaling all of the Hamiltonian by an inverse temperature inside the Boltzmann factor), with Hamiltonian REX.174 Replicas with weaker interactions have a smoother energy surface, leading to better configurational sampling, and if the chosen terms correspond to the slowly varying coordinates, several-fold increase in efficiency can be achieved. If the process is known to be slowed by the breaking and forming of a set of interactions, then those interaction energies could be scaled down to reduce the activation barrier, isolating and limiting the size of the perturbation between 7749


Chemical Reviews

Review

Figure 7. Configurational sampling methods for slow conformational changes. (A) Swarms of trajectories string method for protein conformational change. A hypothetical contoured free energy surface in two dimensions showing optimization from an initial to final path between two ion channel states “open” and “closed”, defined by a string (thick line) with circles representing the images along that curve. A zoomed-in section reveals a swarm of independent unbiased trajectories used for estimating the drift of a single image. Adapted from ref 68. Copyright 2008 American Chemical Society. (B) Example system used in an MSM approach for ligand−protein binding. Libraries of MD simulations showing a ligand (benzamidine) C7 atom position revealing clustering on the protein (trypsin) surface, with kinetics of interconversions between clusters analyzed with an MSM. Reproduced with permission from ref 209. Copyright 2011 National Academy of Sciences USA.

square displacement (RMSD) constraint,192 offering a guess at the crude path for the protein conformational transition, which may subsequently be followed with separate steered MD simulations.181 Targeted and steered MD,193,194 can thus be used to enforce a chosen conformational change away from equilibrium but, importantly, still allow analysis of equilibrium work via the Jarzynski equality,195−198 relating the equilibrium (≡m) free energy change, ΔG, to the nonequilibrium samples of work, W,

neighboring replicas. In the approach of REX with solute tempering (REST), designed for overcoming the demands in biological systems with large quantities of explicit water, deforming the Hamiltonian for each replica is done such that the exchange probabilities do not depend on the number of explicit water molecules.176 A related approach, where solute− solute (e.g., ligand-protein) interactions are scaled via a Monte Carlo scheme, called simulated tempering,177 shows much promise in sampling drug interactions (section 5), as does a new sequential variant called tempered binding,178 where changes occur via sequential Monte Carlo moves in a long MD simulation, finding application using the new DE Shaw Anton/ Anton 2 supercomputers. Extensions to REX also exist that can enhance sampling of defined coordinates, such as combining it with US, termed replica exchange umbrella sampling, or biasexchange MTD,179 where replicas have their window functions exchanged to obtain enhanced sampling and improved PMFs as a function of known coordinates.180 Techniques that apply bias need not be restricted to equilibrium or near-equilibrium sampling but may move the system far from equilibrium to enforce change in a time that is amenable to simulation (but perhaps still short on the functional time scale). By simulating with enough applied force, rare configurations can be encountered, as has been used to study ion channel permeation,181−185 protein conformational changes and unfolding (e.g., refs 181,186−188) as well as protein−ligand binding.181,189 In steered MD, mimicking experimental techniques such as atomic force microscopy or optical tweezers,181 a time-dependent force pulls the system along a defined direction.181 The movement can be imposed by attaching to a spring with increasing force constant,181 or by pulling on an atom or center of mass with a harmonic spring attached to a dummy atom that is moving at constant velocity.190 However, this method requires a decision to be made about the actual vector in which to propagate the system beforehand. An alternative is targeted molecular dynamics, which employs time-dependent holonomic constraints to move from one state to another,191 such as via a relative root-mean-

e−ΔG≡m / kBT = ⟨e−W / kBT ⟩non ≡m .

(41)

Combined with steered MD, this can serve as a viable alternative to US199 and can be considered equivalent to US in the limit of weak and slowly changing force.181 The approach has been used to analyze simulations185,200,201 with appropriate unbiasing of potentials,202 as well as experiments.198,203 The major concern with extraction of reversible work from irreversible processes via the Jarzynski equation is the need for extensive sampling for large systems.85,200 Although speeding up conformational changes by the action of applied force is useful, it does not eliminate the need for similar levels of MD trajectory sampling. There is, in fact, no shortcut to the calculation of the PMF for the system, with similarly large sampling requirements to converge Jarzynski-estimated PMFs (e.g., ref 204). The above methods require knowledge about the process being studied. In section 4.1, we will discuss methods that are designed to reveal the important coordinates underscoring protein conformational changes that can assist in defining the space in which enhanced sampling and nonequilibrium approaches can then be used more reliably. Another approach is to create an extended space of many possible variables that should encompass those most important slow variables (even if the relative importance of those variables is not known) and search for the most probable or minimum free energy pathway of change connecting two end point states. A prominent example is transition path sampling, where an initial path is perturbed via its atomic positions and momenta to generate a 7750


Chemical Reviews

Review

3, some of these concepts are taken further for the case of ion permeation, including how to carefully define PMFs for multiion pores, how to use a computed PMF to make contact with experimental binding and permeability measurements, and how to decompose the free energy to provide a physical basis for what controls ion channel permeation.

new path, with this random walk accepted or rejected via a Monte Carlo algorithm, with a set of related methods.205−208 An important aspect of such approaches is that they do not exhaustively sample the full 3N-dimensional space but only the optimal path itself, effectively reducing the dimensionality of the problem and making sampling more tractable. A drawback is that there may be multiple paths that may go unseen during the iterative solution. One particular approach that sets itself apart is the string method with swarms of trajectories68 to explore a broad energy surface and multiple paths of change, as illustrated in Figure 7A for the case of ion channel activation and described in detail in section 4. In this approach, extensive libraries of short simulations are sent out to explore configurational space, with their mean drifts guiding the system toward the optimal pathway within the reduced space. In recent years, Markov state models (MSMs) have shown enormous promise for studying processes that may involve any number of collective variables that are slow relative to the observational time frame. The approach has been described in detail in various recent reviews and books.210−214 The purpose of MSMs is to observe the possible states of the system and their interconversions without applying bias along preconceived order parameters. In this approach, libraries of simulations are used to explore the microstates of the system, and clustering methods are used to lump those kinetically associated states together to form longer-lived metastates. These simulations need to be long enough to observe transitions between metastates, even if those simulations may not be long enough to capture the entire process. This kinetic information from MD is used to build a discrete-state stochastic model that can describe the long-time evolution of the complex system using a Markovian/memoryless stochastic model.215 The method involves discretization of system’s state into a number of conformationally distinct sets of states based on defined criteria (e.g., positions, angles, principal components, pairwise RMSD values, etc.), and computation of a transition matrix that defines the statistics of transitions between those sets. The eigenvalues and eigenvectors of this matrix can yield the molecular relaxation time scales and the structural changes that occur at those time scales as well as the stationary distribution of coordinates and metastates.214 Applications of MSM include some of the most challenging problems for MD simulation due to the disparity of real and simulation time scales. MSMs have been used to study folding mechanisms for peptides and small proteins (e.g., ref 216) as well as protein conformational changes for native proteins,214 as discussed further in section 4. MSMs also find themselves applicable to the challenging problem of drug binding. While brute force MD can potentially reach times to observe ligand binding events (e.g., using Anton217), often time scales far exceed MD trajectory time scales. For example, Buch et al.209 have used an MSM to examine the trypsin inhibitor benzamidine using libraries of short MD trajectories for binding and dissociation, revealing long-lived intermediates (Figure 7B), as discussed further in section 5.1, in relation to atomistic drug binding studies. We have therefore covered the theoretical approaches for estimating free energies that can be applied to configurational sampling obtained using unbiased and enhanced MD simulations. In sections 3, 4, and 5, we illustrate how to apply these techniques for the calculation of free energies for ion movements, for conformational changes associated with activation, and for a drug binding to an ion channel. In section

2.2. Development of Empirical Force Fields for Ion Channels

The quality of molecular mechanical force fields, commonly used during MD simulations, has improved to allow stable multimicrosecond simulations of ion channels interacting with their physiological and pharmacological environments. The challenge is to accurately model detailed interactions that control the processes of interest, such as those between an ion and the inside of a channel pore. Ensuring accurate energies of interaction is essential to quantitatively capture functional mechanisms and find correspondence with experiment. We survey the current state of atomic force field for simulation of membrane ion channels and explore calculations that attempt to improve on those models with benchmarking calculations to match experiments for ion−protein and drug−membrane interactions. An empirical force field includes definition of the potential energy function, U(q), describing the molecular interactions within the system as a function of the coordinates of the constituent atoms, along with a set of empirically derived parameters.69 Common force fields for proteins (reviewed, e.g., in ref 218) include Amber,219 CHARMM,220,275 GROMOS,222 and OPLS-AA,223 as well as general force fields for organic and biological molecules.224−228 These force fields can be used in various software packages, such as AMBER, 229−231 CHARMM,148,221,232 GROMOS,233,234 GROMACS,235−237 LAMMPS, 238 NAMD, 239−241 TINKER, 242 and DESMOND.243 We focus on atomistic models for simulations of proteins, but coarse-grained (CG) methods,244,245 such as where each amino acid is represented by one or a few particles, can allow simulations on larger time and length scales246−249 to study protein and membrane systems (e.g., refs 247,248,250) up to large protein complexes (e.g., refs 251,252) and have been used to study conformational changes associated with ion channel gating.253−258 Such coarse-graining is inherently approximate, with reduced accuracy for chemically specific interactions, such as those important for ion channel function and drug interactions covered in this review, and lack internal coordinate changes, including secondary structure and side chain isomerizations that may underscore protein activity. However, CG MD can, in principle, be useful for studying chemically specific processes when combined with atomistic simulations in multiscale approaches (e.g., refs 259−264). Here, we summarize the atomistic MD force field with the purpose of highlighting the current state for membrane protein systems and for ion, water, and lipid interactions that may control ion channel function. We begin with a basic introduction to the force field as it applies to membrane protein simulation and then discuss available models for water, lipids, ions, and drugs and their development. 2.2.1. The Atomistic MD Force Field for Membrane Proteins. The atomistic MD energy function includes bonded terms that involve covalently connected atoms separated by one, two, or three bonds, as well as nonbonded terms. The bonded terms involve harmonic expressions for bond and angle vibrations, while torsional (dihedral) potentials are 7751


Chemical Reviews

Review

gas-phase molecular dipoles.270,275 QM interaction energies and interaction distances are thus scaled, leading to balanced protein−protein, protein−water, and water−water interactions.278,279 Note that only the permittivity of free space, ε0, is used in eq 42 (i.e., the relative permittivity or dielectric constant, ε, is set to be 1) due to the explicit inclusion of the dielectric medium, with molecules, including water, able to polarize via reorientation of their permanent dipoles. However, for condensed phases, it has recently been argued that an underlying electronic component of the static relative permittivity (∼2) could be used to create a more accurate model of interactions.280 More generally, the inability of the charge distribution to change depending on its environment can affect the accuracy of the force field, as discussed below. Nonbonded LJ parameters are obtained by reproducing condensed-phase molecular volumes and heats of vaporization,270,281 although may be obtained by fitting to crystal parameters and heats of sublimation.282,283 Application of QM data for LJ optimization requires a high-level of theory and basis sets to account for dispersion interactions. Because of this, a combined QM−empirical approach can be used to account for model compound interactions with noble gases via high-level QM calculations to get relative LJ contributions for individual atoms and to match experimental solvent properties.284,285 Because the LJ parameters are usually determined for self-interactions (i.e., between atoms of the same type), approximate combining rules are used to calculate the parameters between different atom types. In CHARMM and Amber force fields, the Lorentz−Berthelot (LB) mixing rule is used based on a geometric mean (Berthelot rule286)

modeled with sums of periodic terms to account for multiple energy minima and invariance to 360° rotation. In cases where a molecule must be constrained to be planar, improper dihedrals can be maintained with harmonic potentials. Bonded parameters are determined by fitting model compounds to QM calculations or experimental crystal structures, with force constants fitted to vibrational spectra.221 Additional terms exist, such as a Urey−Bradley angle bending harmonic crossterm term, accounting for 1−3 nonbonded interactions as well as corrections for more accurate protein backbone torsional terms to account for deviations from QM and experimental crystallographic data.265−268 This grid-based coupled torsion correction map (CMAP) has been introduced into the CHARMM force field,269 improving conformational sampling and protein stability,270 and has been shown to be important for modeling ion channel permeation (although in the case of channels like gramicidin A (gA), containing left and righthanded amino acids, CMAPs must be reflected to account for the stereochemistry).271 Recent developments include further adjustments of CMAP potentials to better reproduce various experimental protein properties, such as NMR backbone scalar couplings across hydrogen bonds, residual dipolar couplings, and relaxation order parameters, as done in the CHARMM36 force field,220 as well as conformational ensembles of different secondary structures for intrinsically disordered peptides and proteins, as done in its latest modification, CHARMM36m.272 Nonbonded interactions include Coulomb interactions between partial atomic charges as well as a LJ 12−6 potential to model van der Waals interactions (representing the attractive dispersion forces between instantaneous and induced dipoles and the short-ranged electronic repulsion of atoms). We refer to the form of the LJ potential, where the depth of the interaction between atoms i and j is ϵij, and where rijmin represents the location of the energy minimum such that the total nonbonded energy is given by

U (q)nonbonded

ϵij =

ϵiiϵjj

(43)

and arithmetic mean (Lorentz rule)148,287

É ÅÄÅ min 12 l 6 Ñ| o ÅÅij r yz ij rijmin yz ÑÑÑÑo qiqj o o o ÅÅjj ij zz o o jj zz ÑÑo Å 2 = ∑m + ϵ − }, j z j z Å Ñ ij j z j z o Å Ñ j z j z o o ÅÅ rij rij ÑÑo o i K+ ∼ Cs+ ∼ Rb+,1,558,613 being able to select for Na+ over K+ ions with a PNa+/PK+ ∼ 10−20 1. In contrast, the KV-s with narrow single-file SFs are only permeable to a few ions and have a higher relative permeability ratio, PK+/PNa+ of 100− 2501. Ultimately, the kinetic barriers of the multi-ion mechanism control these relative permeability ratios. For example, recent simulations of potassium channel permeation have demonstrated that ion selectivity can directly arise from differing multi-ion processes for each ion.18,20,153,614 However, each of these mechanisms relies on the thermodynamic stabilities of ions in sites inside the channel pore, with the mechanisms distinct because of their differing affinities or particular spatial arrangements (e.g., K+ and Na+ binding at adjacent sites; see Figure 6B).153 Thus, while ultimately the permeabilities of the ions should be either measured or calculated to directly

ΔΔG K+→ Na+ = ΔG Na+:Bulk → Site − ΔG K +:Bulk → Site = (G Na+:Site − G Na+:Bulk ) − (G K+:Site − G K+:Bulk ) = (G Na+:Site − G K+:Site) − (G Na+:Bulk − G K+:Bulk ) = ΔG K+→ Na+:Site − ΔG K +→ Na+:Bulk .

(95)

Calculation of these transformation free energies is done via FEP within the channel pore and in bulk electrolyte, following the methods discussed in section 2.1. In cases where the ion configurations interconvert on long time scales, calculation of such a free energy can be complicated. In a multi-ion channel, one often has to trap the particular multi-ion configuration, such as, for example, S0/S2/S4 for three K+ ions in the KcsA channel (e.g., Figure 6C) and examine perturbations of one of the particular ions, for example. However, problems arise when the transformation from a state such as K+ − K+ − K+ to K+ − Na+ − K+ leads the system to become unstable. This represents a situation, such as in Figure 6B, where a change in coupling 7771


Chemical Reviews

Review

Figure 18. Controls of K+/Na+ selectivity. (A) Interaction and strain energies in a microdroplet. (B) Effects of number and type of ligand on ΔΔGK+→Na+ using the CHARMM27 and AMBER (parentheses) force fields. Reproduced with permission from ref 622. Copyright 2007 Rockefeller Press.

parameter λ leads to a shift in the position of the minimum in the free energy surface. For K+ → Na+ in a K+ channel, this can mean the ion moving from a “cage” site to the neighboring “plane” site and can happen on long time scales to approach equilibrium. As discussed in section 2.1, one solution is to combine FEP with US to ensure the slowest degrees of freedom are captured in the simulations. See refs 152,153 for more discussion and demonstration of improved convergence of FEP calculations, as well as ref 20 for unbiased observations of such multi-ion effects. There is a competition between the loss in interactions with water molecules, when an ion is forced to dehydrate to enter a narrow SF and gain in interactions with the amino acids in the SF.615 To seek a general understanding of the chemical and physical forces underscoring a relative thermodynamic preference for a particular binding site and enable systematic variations in coordination chemistry and number, one can focus in on models of the ion solvation microenvironment. The chemical properties of the ligands coordinating the ions are important and underscore key differences in KV and NaV channels, for instance. According to Eisenman field strength theory,615 an ion with higher charge density prefers a strong field strength ligand, as it can make up for the large partial dehydration cost of entering the narrow channel, whereas a ligand with low field strength cannot.615 The number of coordinating ligands also play an important role, as increasing the number of coordinating ligands can make up for a weaker ion−ligand interaction.616 Examples of high field strength ligands in cation channels are Asp and Glu side chains, present in Na+-selective NaV channels, with their carboxylate groups expected to be deprotonated in aqueous solution at neutral pH. In contrast, protein backbone carbonyl groups, lining the narrow KV SF, offer weaker fields that favor the larger K+ ion. While larger coordination number can make up for weaker interactions,616 there is competition between the favorable ion−ligand energy and the unfavorable ligand−ligand

repulsions.473 K+ prefers a higher coordination number, as seen in KV channels, lined with sites formed by cages of eight backbone carbonyl oxygen atoms (Figure 6B). Decreased number of coordinating carbonyls creates a Na+ selective site,473 as seen inside potassium channels between the cages where planes of four carbonyls can bind Na+ ions153 (see Figure 18A,B). Reduced microenvironment models have yielded basic insight into the selective binding of ions by molecules of varying number, chemistry, and flexibility.473,617−621 Such an approach is meaningful because long-range interactions do not directly control the relative free energy of similar ions.16 Studies have shown that including a charged carboxylate in the microenvironment tends to favor Na+.620 Observations for ion channels have suggested that such selection involves complexes of ions with carboxylates (Figure 18B). One can systematically vary such complexes in different environments and can also decompose energetics to explain the driving forces for ion selectivity. In the creation of a microenvironment model, one must select a system of molecules representative of the ion’s environment and restrain them in a fashion that mimics the host protein. There are two types of restraining forces: ones that keep the complex associated with allow solvation and ones that influence the ability of molecules to coordinate ions. In the absence of the latter, one obtains a “microdroplet” model, where selectivity is controlled by the number and type of ligands.16 One can examine this “inherent selectivity”, as has been done in several previous models,473,617−621 before imposing protein influences for specific cases. Over the last two decades, studies of K+ selectivity have concentrated on such reduced models, theoretically leading to two limiting possibilities for maintaining selectivity.16,621 In one limit, it is the protein’s constraining forces on the site that maintain a nearly perfect (sub-Å) geometry to thermodynamically select for K+ ions. This traditional view623 is that ion selectivity can arise from a “snug fit”, where the SF fits one ion 7772


Chemical Reviews

Review

perfectly but is unable to adapt to the other.624 However, ion channels are typically highly flexible and therefore able to adapt to sub-Å differences in ion size.21,22,155,473,480 In the other limit, the ligands are unconstrained by protein interactions, albeit confined in a microdroplet, but still able to maintain selectivity. Here the crowded environment around the ion leads to opposing terms associated with attractive ion−ligand and ligand−ligand strain energies473 that maintain a preference for the larger K+ ion, as indicated in Figure 18A. Whereas some proteins, such as the leucine transporter LeuT, rely on the protein’s architectural forces to maintain a selective site geometry to some extent, others, such as KcsA, do not.16 Figure 18B shows the effects of number and type of ligand in the microdroplet model. It has been argued that provided such a cluster of molecules is maintained loosely, a naturally flexible carbonyl-rich site is naturally K+-selective, as a result of the ligand−ligand strains imposed on the crowded complex. However, it has been shown that this selectivity can be lost by reducing the number of carbonyls (as observed within K+ channels153), replacing them with water molecules, and is consistently reversed by introducing a higher field strength charged carboxylate ligand.622 Another way to offset the ligand−ligand cost of formation of multiple high field ligand complexes, so as to favor Na+ ions, is to introduce additional ions to the site, increasing the attractive ion−ligand energy and outweighing the ligand−ligand and ion−ion terms.625 This is what has recently been seen in NaV channels with multiple Na+ ions binding concurrently to carboxylate groups.626 Figure 22 (described below), for example, shows how such a particular ion complex can be critical to permeation events and can explain ion selectivity for the channel (section 3.2). FEP can be used to study the relative stabilities of different ions, such as Na+ and K+, in such defined complexes. In the absence of dedicated reduced models, one can instead isolate the complex from simulations of the full system through identification of dominant configurations (such as via cluster analysis of trajectories478). To ensure the complex is maintained, flat-bottom/square well potentials can be used in collective coordinates that best define the complex, such as distance between ion and ligand center of mass, for example. The constrained FEP becomes, from eq 27, ΔGcons = − kT ln⟨e−ΔU / kBT ⟩0;cons

ΔΔGcorrection ≈ − kBT ln

⟨e ucomplex / kBT ⟩1 − cons ⟨e ucomplex / kBT ⟩0 − cons

(97)

and simply entails averaging inverse Boltzmann factors of terms from a time series of the collective variable obtained from the constrained FEP simulations. However, the constraint could be corrected through analysis of all simulations (using all values of the coupling parameter λ). One way to achieve this is to rebuild histograms from individual simulations by reweighting configurations according to the constraint energy. Instead of binning values of 0 or 1 inside each bin, as is normally done to form a histogram, the value binned is the inverse Boltzmann factor to the constraint energy, eucomplex / kBT . By doing so, configurations that were disfavored because of that constraint in the constrained FEP simulation are reweighted to be more favored. One can show that ρλ − corrected (u) = ⟨e ucomplex / kBT δ(u − ΔU (q))⟩λ − cons × e−Fcomplex,λ / kBT ,

(98)

where Fcomplex,λ is the free energy cost of introducing the constraint to that window, with the multiplier in the form of a Boltzmann factor on the right-hand side representing a normalization constant for each window, obtained over the whole trajectory. Such analysis requires side-by-side time series of energy perturbations and collective variable values for the complex constraint, modifying histograms that could then be fed into WHAM analysis (section 2.1). However, in general, a prudent choice of constraint can lead to small corrections that can be adequately captured by eq 97. The above sets up a framework for calculation of alchemical transformation free energies that can report on the relative free energies of binding to sites or identified ion complexes related to ion permeation. As physiologically the relative permeabilities of channels may correspond to energy differences of the order of 1 or a few kcal/mol, it is clear that the value in such calculations hinge on accurate models for ion−ligand interactions. In the case study to follow, it will be seen that sodium channel permeation involves highly specific multi-ion/ multicarboxylate complex formations that were not originally designed for within most force fields and depend critically on the interaction of an ion with a carboxylate group. Efforts to test and improve such parameters were discussed in section 2.2. However, the quality of the force field remains a critical factor in studies of relative free energies governing ion channel selectivity.

(96)

simulated in the presence of a flat-bottom constraint, ucomplex(q), applied to some set of coordinates, within q, that maintains the complex. However, this constraint may have to be tight because the complex may become unstable in the other end state (e.g., K+ may not exist in a minimum in that complex, whereas Na+ does). In the case study of section 3.2, we will see how Na+ forms a stable two-ion/two-carboxylate complex, but that this complex is quite unstable for K+. In that case, even a very sharply defined square well, based on steep half-harmonic boundaries, ucomplex(q), will be in action for a good fraction of the trajectory. This means that non-negligible free energy was involved in introducing the constraint into the simulation. As it potentially has a different cost in each end point state, it must therefore be unbiased postsimulation. To a first approximation, this amounts to computing a postsimulation FEP correction involving the turning off of the constraint term for states 1 and 0 and finding the difference

3.2. Case Study. Sodium Channel Conduction and Selectivity: from Bacteria to Human

NaV channels play an essential role in electrical signaling in the body and are necessary for functions such as heart beat, muscle contraction, and brain activity.1 These functions depend on the ability of these channels to select for Na+ while discriminating against other ions with high fidelity. The current understanding of Na+ channel selectivity is limited, but recent high-resolution structures present an opportunity for elucidating these mechanisms using MD simulations and resultant free energy calculations. In this section, we summarize recent simulations that employ either long unbiased simulations or follow biased PMF and FEP calculations to examine ion binding and permeation (section 3.1), with each approach yielding basic insight into the mechanisms of Na+ conduction and selectivity. 7773


Chemical Reviews

Review

Figure 19. (A) Structure of a prototypical BacNaV channel NaVAb, showing two of the four subunits for clarity, with each segment (S1−S6), S4− S5 linker, P1 and P2 helices, the SF, and intracellular gate marked. (B) Sequence alignment of NaChBac, NaVAb, CaVAb, and NaV1.2 channel SF regions. The amino acids are colored according to their properties using the Zappo coloring scheme.630

Figure 20. Sketch of the topology of a bacterial (top) and mammalian (bottom) NaV. The mammalian sodium channel has four domains (DI−IV) linked into one long polypeptide chain, whereas the bacterial sodium channel has four identical subunits. Each domain/subunit consists of six helical TM segments, S1−S6, and loops in between. The VSD, S1−S4, is shown in blue, the PD, S5−S6, in purple, the P-loops in orange, and the DIII−DIV loop critical for channel inactivation with a signature IFM motif in red.486

The first high-resolution structure of a Na+ channel (the bacterial NavAb channel from Arcobater Butzleri (6); Figure 19) is reminiscent of the voltage gated K+ (Kv) channels, with four identical monomers comprised of a voltage sensor domain (VSD) and a pore domain (PD), symmetrically organized around a central pore. In contrast to the long narrow, ion-sized SF of the K+ channel (compared in Figure 12), where the backbone carbonyls of the amino acids directly interact with multiple K+ ions in single-file,1 the Na+ channel SF is slightly wider and shorter and lined not only by protein backbone (residues 175−178) but also by a ring of glutamate side chains (E177). The NaV SF is thought to allow partly hydrated Na+ ions to cross627 and is wide enough to allow protein side chains to be directed into the pore to assist in coordinating the permeating ions.1,628,629 Importantly, we will discuss observations which reveal considerable flexibility of these side chains, which is not only critical to the translocation of ions but enables nonsingle file permeation events that are distinct from those seen in K+ channels. In contrast, early simulations had suggested that Na+/K+ selectivity may arise from small structural differences,557 where Na+ (but not K+) is “snugly fitted” in the narrow SF, implying a traditional mechanism

relying on size-based discrimination. Such a mechanism of ion discrimination appears inconsistent with the high flexibility of those coordinating side chains.155 Moreover, as discussed in section 2.2, the intimate involvement of these charged side chains means that ion permeation hinges on the precise energetics of ion complex formation that is dependent on the choice of force field, as was the case in K+ channels. As we learnt from studies of K+ channels, permeation is a multi-ion process and mechanisms of discrimination must capture those differing processes for each ion species, even if the underlying causes may involve the thermodynamics of binding to particular ion microenvironments within the pore. Because of the combination of side chain conformational changes and slow interconversion of multi-ion configurations, this case study presents a significant challenge to MD simulations as they attempt to observe permeation events. Yet brute force MD, such as available on the Anton/Anton 2 supercomputers, offers an advantage in reaching sufficient time scales to begin to cover the relevant processes as we seek to understand the ability of a Na+ channel’s SF to favor Na+ over K+ ions. These mechanisms are, however, currently obscured by seemingly unrelated chemistries of different Na+-selective 7774


Chemical Reviews

Review

Figure 21. (A) Radial−axial densities showing Na+ (blue), water (red), and carboxylate/carbonyl oxygen atoms (purple) confirm proposed locations of binding sites in NaV. Reproduced from ref 576. Copyright 2011 American Chemical Society. (B) 2-ion PMF from US showing the pathway for conduction. Reproduced with permission from ref 557. Copyright 2012 American Chemical Society.

make up the PD that is symmetrically positioned to create an ion-conducting pore. S5−S6 are linked by a pore-loop (Ploop), made up by two short helices, P1 and P2, connected by the narrow SF and extracellular loops. The SF is made up a signature Asp-Glu-Lys-Ala (DEKA) ring, shown to be crucial for Na+ selectivity, and a vestibular negatively charged ring EEDD.482 The domains are linked together by large intracellular loops,634 known to control channel inactivation process, which renders the channel in a state unresponsive to further stimuli (especially, DIII-DIV; Figure 20, red loop)634 as well as drug binding635 and self-regulation.636 Recent near-atomic resolution structures of eukaryotic NaVs, NaVPaS from cockroach,637 NaV1.4 from electric eel,638 and

channels, which can be explored by examining bacterial and human NaV channels and evolutionarily, structurally, and functionally unrelated ASIC channels that all exhibit similar preference for Na+ ions. 3.2.1. Mammalian and Bacterial NaV Channels. There are nine members of the mammalian Nav family, Nav1.1Nav1.9,631 sharing several features structurally and functionally.632 These channels consist of four domains, DI-DIV, linked together by intracellular loops, to form one long polypeptide chain (Figure 20).486 Each of the domains consists of six helical TM spanning segments (S1−S6). S1−S4 form the voltage-sensing domain (VSD) responsible for gating of the channel, with S4 important for voltage sensing.1,633 S5−S6 7775


Chemical Reviews

Review

Figure 22. (A) Crystal structure showing static SF conformation with presumed ion binding sites SHFS, SCEN, and SIN and E177 side chains reaching upward. (B) Structure from MD simulations showing an SF conformation with a multi-ion complex formed by flexible E177 side chains. (C) Network representation revealing diverse Na+ binding modes, with different colored 1−3 ion states 1′, 2, 2′, and 3 corresponding to a large number of substates with different numbers of carboxylate and carbonyl O atoms coordinating the ions, shown with circle size representing abundance, lines representing transitions during simulations and with substates clustered to form the different states illustrated in the insets. Reproduced with permission from ref 538. Copyright 2013 National Academy of Sciences USA.

the new NaV1.4487 and Nav1.7639,640 channels from human, solved using cryo-EM, have sparked a new phase of investigation into mammalian NaV-s. However, in the past decade, studies have made use of high-resolution X-ray structures of bacterial channels (BacNaV), such as the first one, NaVAb,482 followed by a set of structures including NaVRh,483 NaVAep1484 and the proposed open NaVMs30,485 and NaVAb.31 The BacNaV topology is simpler than that of the mammalian Navs, being made up of four identical subunits (Figure 20) and otherwise sharing several features with mammalian NaVs.486 This includes their ion selectivity fingerprint (Li+ ∼ Na+ > K+ ∼ Cs+ ∼ Rb+), with relative Na+/K+ permeabilities, PNa+/PK+ of ∼5−170558,613,641 and ∼10−201, and similar single channel conductances of ∼12 pS613 and ∼4−24 pS,1 for bacterial and mammalian Navs, respectively. This is a surprising result given the distinct sequences in bacterial and mammalian channels, as we will discuss here. The structures of bacterial channels, as illustrated in Figure 19A,486 have told us much about how this family of channels function. In particular, they have revealed the structure of the narrow SF, made up of a loop between P1 and P2, which is

known to control ion permeation and discrimination. One of the most pronounced differences between bacterial and mammalian channels is the amino acids in the SF (Figure 19B). Instead of a DEKA ring, the bacterial channel has a ring of four Glu side chains, EEEE. These residues are conserved throughout all BacNaVs (except for NaVRh, which in its place has a ring of Ser at this location, replaced by a ring of Glu three residues higher in the SF) and have been shown to be responsible for Na+ selectivity in BacNaVs.486,628,629 This EEEE ring creates a high-field strength binding site (SHFS) that is thought to be ideal for binding small cations such as Na+,482 although it is also a characteristic of Ca2+-selective channels.642 The SF of the bacterial channel NaVAb contains these four Glu177 side chains as well as the backbone carbonyls of Thr175 and Leu176 beneath.482 These residues form three sites for binding Na+ ions; the EEEE ring forming site SHFS, the backbone carbonyls of Leu forming central site, SCEN, and the carbonyl and hydroxyls of Thr forming an inner site, SIN,480,482,558,643 as indicated below in Figures 21a and 22a. An additional external site in the outer vestibule has also been identified in some BacNavs.484,557 Because of the shorter and wider SF in NaVs, being lined by amino acid side chains, rather 7776


Chemical Reviews

Review

Figure 23. Na+ in the SF of NaVAb. 2D free energy projections showing (A) two-ion occupancy and (B) three-ion occupancy. z1 corresponds to the z position of the bottom ion, z23 is the z position of COM of the top two ions, and z3 is the z position of the top ion. Snapshots with the EEEE ring in purple indicate the corresponding Na+ ion (orange balls) movements. Reproduced with permission from ref 478. Copyright 2018 PLoS.

tional changes, to experience the activation barriers that control conduction as well as the flexible response of the protein that responds to the presence of those ions. Despite the bacterial channel NaVAb being in a closed state, with narrow lower pore (Figure 19), one may sample at equilibrium the movements of multiple ions to describe the permeation mechanism. The mechanism may, however, be altered by the application of significant membrane potentials, e.g, simulations of a proposed-open structure (NaVMs; with pore held open with strong constraints) under the action of large driving membrane potentials (up to 665 mV), observed reduced ion occupancy compared to equilibrium simulations.559 The situation with a closed pore can be considered relevant to the limit of low physiological membrane potentials. 3.2.2. Sodium-Selective Ion Permeation in Bacterial NaV Channels. Early MD simulations of bacterial NaVs

than solely by backbone, NaVs are expected to conduct ions and select for Na+ with a different mechanism. BacNaVs select for Na+ over K+ with relative permeabilities of PNa+/PK+ ∼ 5− 170 and Na+ over Ca2+ with PNa+/PCa2+ ∼ 7−70.558,613,641 Site directed mutagenesis has provided some clues by demonstrating that mutation of the Glu-s to Asp-s decreases Na+ over K+ selectivity,556 suggesting the role of precise location of charged side chains, while addition of Glu or Asp to the SF increases Ca2+ selectivity,491,641 indicating that net charge is critical. These ideas are based on a thermodynamic preference for one ion over the other, however, it is important to also take the kinetics of the multi-ion movement into account, as demonstrated for K+ channels.14,18,19,21,153,480,614,643 To understand ion conduction and selectivity, we must study ion channels on time scales where ions are allowed to move through the SF, capturing their preferred multi-ion configura7777


Chemical Reviews

Review

Figure 24. Comparison of Na+ and K+ movements in the SF of NaVAb. (A,B) Time series of ion movements for Na+ and K+ (different ions of the same species as different colors), revealing predominance of three and two ions, respectively, with better defined binding sites for Na+. Insets on the right show representative configurations (Na+, pink; K+, green balls).

revealed the role of the SF in permeation.481,557,559,576,644−646 Carnevale et al. described the ion binding sites and showed that the Na+ permeation differs greatly from that previously described for KV channels, with Na+ ions permeating the NaV SF off axis while surrounded by water molecules that are decoupled from ion conduction, Figure 21A,576 clearly distinct from the single-file permeation familiar in K+ channels. It was further shown that the chemistry of the EEEE ring and the particular SF geometry create an environment where Na+ ions can permeate in a partly hydrated state.481,557,559 Studies using US free energy simulations (sections 2.1 and 3.1), allowing the mapping out of the positions of one or two ions across the SF, have shed light in the multi-ion permeation mechanism for Na+. They have shown that, while a single ion encounters a high barrier to permeation, a second ion can help promote ion movement via a loosely coupled knock-on mechanism,481,557 by concentrating in the SHFS site, as seen in Figure 21B.481 FEP simulations (sections 2.1 and 3.1) have been used to investigate the origin of the two-ion knock-on mechanism by Zhang et al., who showed that the entrance of a second ion reduces the binding affinity and destabilizes the resident ion.483 Several simulations used enhanced sampling methods to provide sampling of two-ion movements on short (multinanosecond) scales.481,557,558 However, on such time scales, the channel structure does not deviate considerably from the X-ray crystal structure to observe Glu side chain involvement and cannot capture the mix of pore occupancy states. It was suggested, however, that the barrier encountered might be reduced by the entrance of a third ion to enable more efficient conduction,558 requiring more extensive investigation of the multi-ion free energy surface. It was not until longer (microsecond-order) simulations were performed that conformational isomerization of the Glu side chains forming the SHFS site were observed,480,538,647,648 dependent on the varying pore occupancy by ions, as shown in Figure 22. Here the initial E177 upward configuration, Hbonded to Ser178 (Figure 22A), breaks away and the Glu side

chains can then take on different conformations to coordinate 1−3 ions (two-ion example in Figure 22B), exploring a range of substates with varying carboxylate coordination, with an example illustrated in Figure 22C.558 This figure demonstrates the complexity of protein involvement during the conduction process, with states characterized by a number of ions having a large number of interconnected microstates, distinguishable by differing carboxyl and carbonyl coordination sampled at equilibrium (see caption and ref 558 for full methodological details). These changing rotamers allowed for the dynamic formation and breaking of multi-ion/multicarboxylate complexes that facilitate Na+ conduction,480,538 also proposed to be central to Na+ selectivity, owing to their inherent preference for Na+ over K+.615 Furthermore, these unbiased simulations uncovered a low energy barrier for conduction, where occupancy fluctuates between two and three ions480,538,649 leading to a knock-on mechanism, where the Glu side chains facilitate rapid translocation of the ions through the SF.480,538,647,648 This emphasizes the benefits of long equilibrium sampling of multi-ion configurations and flexible protein response for computing meaningful free energy surfaces to understand ion permeation. With good equilibrium sampling of varying channel occupancy, one can divide the simulation up to produce PMFs with defined occupancy numbers, as described by eq 58. In the multimicrosecond simulations of NaVAb by Boiteux et al.,480 two ions were often seen occupying the SF, but the ions appeared trapped, bounded within a few Å range (Figure 23A). The two ions move between SHFS, SCEN, and SIN, and occupy one of the states A, B2, or C2. The SHFS site is occupied by at least one ion all of the time, with the bottom ion in SIN occupied in state A, SCEN in state B2, and with both ions in SHFS in state C2 (Figure 23A; insets). It is not until a third ion enters the SF that complete permeation can be observed, with an ion translocating the whole length of the SF (Figure 23B). Permeation in this 3-ion case can happen via two separate mechanisms, involving four distinct states: A, B3, C3 and D3, as 7778


Chemical Reviews

Review

Figure 25. Models of human NaVs from (A) Xia et al., showing how the SF is permeable only when the DEKA Lys is bound to the side of the SF. Reproduced with permission from ref 643. Copyright 2013 Elsevier. (B) Mahdavi et al., showing one ion binding in the lower SF and two ions in the vestibule. Reproduced with permission from ref 659. Copyright 2015 PLoS. (C) Flood et al., showing how Na+ permeation relies on the active participation of the DEKA Lys, which acts as an additional ion and forms a complex with Na+ and the SF Glu-s and Asp-s. Reproduced with permission from ref 478. Copyright 2018 PLoS. Inset above shows which residues were included from the human NaVs in the respective models.

comparison of ion permeabilities has not yet been achieved. Moreover, it has been shown that the propensity for multiple ion binding to carboxylates is sensitive to the chosen parameters.479 Given the lack of tight constraints imposed on these parameters by available experimental and QM calculations (section 2.2), there is room for future force field developments in improved ion−carboxylate interactions to distinguish between the efficiencies of multiple Na+ and K+ mechanisms. 3.2.3. Sodium-Selective Ion Permeation in Mammalian NaV Channels. Mammalian channels benefit from a long history of experimental characterization because of their relevance to human health but until recently had lacked structural definition to enable studies of function at the atomic level. The solution of multiple bacterial structures in the last eight years created a platform from which one can extrapolate to the mammalian channels with models to explain NaV conduction and selectivity.650,651 Here we discuss how models to date have benefited from large-scale MD simulation to shed light on the fundamentals of NaV conduction and selectivity. However, several recent structures of eukaryotic Na V channels487,637−640 in the last months, now including human channel structures,487,639,640 represent exciting breakthroughs to improve on this knowledge in the years to come.

seen in the insets of Figure 23B (where A is the same state in the two and threeion states; with the third ion entering from above). As the third ion enters the SF, it binds just above SHFS and the middle ion, while the bottom ion occupies SIN, forming state B3. In state C3, the top and middle ion are both bound by the SHFS EEEE ring in a multi-ion/multicarboxylate complex. This repels the bottom ion, which enters the cavity, forming state D3 (solid line). Alternatively, the top ion can pass by the middle ion and push the bottom ion into the cavity (dashed line). Both of these paths encounter low (

Atomistic Simulations of Membrane Ion Channel Conduction, Gating

Recommend Documents