Selective ion binding and transport by membrane proteins – A computational perspective

Accepted Manuscript Selective Ion Binding and Transport by Membrane Proteins – A Computational Perspective Hristina R. Zhekova, Van Ngo, Mauricio Chagas da Silva, Dennis Salahub, Sergei Noskov PII: DOI: Reference:

S0010-8545(16)30503-3 http://dx.doi.org/10.1016/j.ccr.2017.03.019 CCR 112426

To appear in:

Coordination Chemistry Reviews

Received Date: Revised Date: Accepted Date:

5 December 2016 23 March 2017 24 March 2017

Please cite this article as: H.R. Zhekova, V. Ngo, M.C. da Silva, D. Salahub, S. Noskov, Selective Ion Binding and Transport by Membrane Proteins – A Computational Perspective, Coordination Chemistry Reviews (2017), doi: http://dx.doi.org/10.1016/j.ccr.2017.03.019

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Selective Ion Binding and Transport by Membrane Proteins – A Computational Perspective

Hristina R. Zhekovaa,b, Van Ngoa,b, Mauricio Chagas da Silvaa,b,c, Dennis Salahuba,c,d,e*, and Sergei Noskova,b* a

CMS - Centre for Molecular Simulation b

Department of Biological Sciences c

d

Department of Chemistry,

IQST - Institute for Quantum Science and Technology

University of Calgary, 2500 University Drive N.W., Calgary, Alberta, Canada T2N 1N4 e

College of Chemistry and Chemical Engineering, Henan University of Technology, No. 100, Lian Hua Street, High-Tech Development Zone, Zhengzhou 450001, P.R. China

*All authors contributed equally to the best of their abilities. Correspondence should be addressed to Sergei Noskov ([email protected]) or Dennis R. Salahub ([email protected])

1

Abstract Inorganic ions are critical for cellular function and require an efficient mechanism of transport through the cellular membrane. Most often the transport of ions occurs through proteins known as ion channels and transporters. Ion binding and permeation through these proteins is a complicated process that is still under investigation with a wide range of experimental and theoretical methods. Here we present an overview of some of the competing theories of ion transport with special emphasis on the theoretical methods used for the elucidation of the energetics of ion selectivity, coordination and permeation. A large part of the review is dedicated to potassium and sodium channels and transporters, which are among the best studied biological transport systems and provide a frame of reference for all other ion-protein interactions. In addition, we summarize the computational work done on the transport of several other small inorganic ions (calcium, magnesium, chloride, inorganic phosphate). Our aim is to provide a general overview of the current state of knowledge on biological ion-transport phenomena and to evaluate the capabilities of modern computational methods when applied to ion transport. We also strive to draw attention to some underdeveloped areas of ion transport that require further investigation. Keywords Inorganic ions, ion channels and transporters, ion selectivity, ion permeation, quantum mechanics, molecular dynamics

Funding Information HZ was supported by a University of Calgary Eyes High Postdoctoral Fellowship. VN was supported by AIHS and CIHR grants. MCdS was supported by AITF and by NSERC through grants to DRS and SYN. Activities in the SYN and DRS labs were supported by the National Sciences and Engineering Research Council (NSERC) (Discovery Grants RGPIN-315019 to SYN and RGPIN-2014-06606 to DRS) and the Alberta Innovates Technology Futures Strategic Chair in (Bio)Molecular Simulation. Computational resources from Compute Canada/WestGrid are gratefully acknowledged.

2

Abbreviations ABF – Adaptive Biasing Force ADP – Adenosine Diphosphate aHL – alpha-Hemolysin AMBER - Assisted Model Building with Energy Refinement AMP – Adenosine Monophosphate AMOEBA - Atomic Multipole Optimized Energetics for Biomolecular Simulation ATP – Adenosine Triphosphate BD – Brownian Dynamics BetP – Betaine transporter from Corynebacterium glutamicum Cav – Voltage-gated Ca channel CBS - Cystathionine Beta-Synthase CD – Circular Dichroism CHARMM - Chemistry at HARvard Macromolecular Mechanics CLC – Chloride Channels and Transporters CLC-ec1 – Chloride transporter from Escherichia coli CorA – Magnesium transporter from Escherichia coli Cryo-EM – Cryo-Electron Microscopy DFT – Density Functional Theory DFTB – Density Functional Tight Binding DNA – Deoxyribonucleic Acid EAAT – Excitatory Amino Acid Transporters EPR – Electron Paramagnetic Resonance FEP _ Free Energy Perturbation FRET – Förster Resonance Energy Transfer G3P – Glycerol-3-Phosphate gA – Gramicidin A GABA - gamma-aminobutyric acid GGA – Generalized Gradient Approximation GHK - Goldman-Hodgkin-Katz GltPh – Aspartate Transporter from Pyrococcus horikoshii GlpT - Glycerol-3-Phosphate Transporter from Escherichia coli hDAT – Human Dopamine Transporter HF – Hartree-Fock HP - Hairpin IR – Infrared KcsA – Voltage-Gated Potassium Channel from Streptomyces lividans Kv – Voltage-Gated Potassium Channel LDA – Local Density Approximation LeuT – Leucine Transporter from Aquifex aeolicus MC – Monte Carlo MCD – Magnetic Circular Dichroism MD – Molecular Dynamics MgtE – Magnesium Transporter from Thermus thermophilus Mhp1 – Hydantoin Transporter from Microbacterium liquefaciens MM – Molecular Mechanics 3

MscL – Large Mechano-Sensitive Channel MscS – Small Mechano-Sensitive Channel MSM – Markov State Modelling NaChBac – Voltage-Gated Sodium Channel from Bacillus halodurans NaK – Voltage-Gated Sodium Potassium Channel from Bacillus cereus Nav – Voltage-Gated Sodium Channels NavAb – Voltage-Gated Sodium Channel from Arcobacter butzleri NavM - Voltage-Gated Sodium Channel from Magnetococcus marinus NavRh - Voltage-Gated Sodium Channel from alpha proteobacterium HIMB114 NBC – Sodium-coupled bicarbonate cotransporter NCX – Sodium Calcium Exchangers NCX_Mj - Sodium Calcium Exchanger from Methanococcus jannaschii NMR – Nuclear Magnetic Resonance OmpF – Outer Membrane Protein F from Escherichia Coli OprP – Phosphate Transporting Outer Membrane Protein from Pseudomonas aeruginosa PB – Poisson-Boltzmann Pi – Inorganic Phosphate PMF – Potential of Mean Force PNP – Poisson-Nernst-Planck QM – Quantum Mechanics QM/MM – Hybrid Quantum Mechanics/Molecular Mechanics RMSD – Root Mean Square Deviation RNA – Ribonucleic Acid SBP – Sulfate Binding Protein SERCA – Sarco/endoplasmic Reticulum Ca2+-ATPase SF – Selectivity Filter SMD – Steered Molecular Dynamics smFRET – Single Molecule FRET TI – Thermodynamic Integration TM – Transmembrane Segment TMD – Targeted Molecular Dynamics US – Umbrella Sampling UV/VIS – Ultra-Violet/Visible range VDAC – Voltage Dependent Anion Channel vSGLT – Sodium Galactose Transporter from Vibrio parahaemolyticus WHAM – Weighted Histogram Analysis Method

4

Contents Abstract ....................................................................................................................................................................................................2 Keywords ................................................................................................................................................................................................2 Abbreviations ........................................................................................................................................................................................3 1.

Introduction ..................................................................................................................................................................................7 1.1 Ions as key ingredients of life. ...........................................................................................................................................7 1.2 Ion selectivity as a key functionality of many transmembrane proteins. .........................................................7 1.3 Coupled transport of ions and substrates as a key functionality of many transporters. .............................8 1.4 All-atom simulations for ion transport experiments and vice versa. .................................................................9 1.5 Why another review on ion binding to proteins? .......................................................................................................9 1.6 Organization of the review ...............................................................................................................................................10

2.

Potassium Channels: Versatile platform for selective binding and permeation studies ............................11 2.1 Summary..................................................................................................................................................................................11 2.2 Mechanistic models.............................................................................................................................................................11 2.3. Ligands, solvents and ions in the confinement of K+ channel selectivity filter ........................................14 2.4 Challenges in coupling ion binding and conformational dynamics ................................................................16 2.5. Computational models for thermodynamics and kinetics of selective ion binding.................................17 2.6. Ion fluxes and conductance ............................................................................................................................................19 2.7. Challenges and unresolved issues in computational modelling of selective ion transport across K+ channels ...........................................................................................................................................................................................21

3.

Sodium Channels: Charged Ligands, Water and Loose Coordination............................................................24 3.1 Summary..................................................................................................................................................................................24 3.2 Mechanistic models and thermodynamics .................................................................................................................24 3.3 Coordination and importance of polarization effects in the selectivity filter of Nav channels ............26 3.4 Computational Electrophysiology Studies ................................................................................................................27

4.

Sodium and Potassium Transporters: Computationally Demanding ................................................................29 4.1 Summary..................................................................................................................................................................................29 4.2 Active transport of Na+ and K+.......................................................................................................................................29 4.3 Structure and functionality of secondary transporters (GltPh and LeuT): Where ions come to play .29 4.4 Computational challenges ................................................................................................................................................31

5.

Divalent cation channels and transporters: Substantial electronic effects are necessary..........................33 5.1 Summary..................................................................................................................................................................................33 5.2 Transport of Ca2+ and Mg2+ .............................................................................................................................................33 5.3 Cav channels ...........................................................................................................................................................................33 5.4 Primary transport with SERCA......................................................................................................................................34

5

5.5 Secondary transport with NCX ......................................................................................................................................35 5.6 Mg2+ transport by MgtE and CorA ...............................................................................................................................35 5.7 Computational challenges ................................................................................................................................................36 6.

Anion channels and transporters: unexplored territory...........................................................................................37 6.1 Summary..................................................................................................................................................................................37 6.2 Chloride channels and transporters...............................................................................................................................37 6.3 Inorganic Phosphate ............................................................................................................................................................39 6.4. Other ions – F-, Br-, I-, HCO3-, SO42- ..........................................................................................................................40

7.

Outlook and Challenges .......................................................................................................................................................41 7.1 Summary..................................................................................................................................................................................41

APPENDIX - Computational methods used in studies of ion transport ....................................................................44 A1. Quantum Mechanical Methods .....................................................................................................................................44 A2. MM methods .........................................................................................................................................................................45 A3. Computational Electrophysiology ...............................................................................................................................48 A4. QM/MM..................................................................................................................................................................................50 Figure Captions..................................................................................................................................................................................52 References: ..........................................................................................................................................................................................60

6

1. Introduction 1.1 Ions as key ingredients of life. It is common knowledge that the chemistry of life is governed by organic compounds and macromolecules, such as proteins, nucleic acids and lipids. Life, however, would be impossible without a number of small inorganic ions, involved in cell signaling, neural impulse conduction, muscle contraction, pH and osmotic pressure balance, metabolic synthesis and many other biological processes[1]. In humans the physiologically important inorganic ions, including Na+, K+, Ca2+, Mg2+, Cl–, HCO3–, and PO43– (most often in the forms of H2PO4– and HPO42–) are dissolved in the intracellular (cytosol) and extracellular fluid compartments. The constant ionic exchange between the intracellular and extracellular compartments ensures the proper distribution of ions in the human body and is essential for cellular homeostasis. The concentration of these ions in the cytosol and the extracellular fluids is generally different, serving as a source of electrochemical energy for various cellular processes (see Table 1)[2]. For example, K+ ions are found in larger amounts in the cytosol, while Na+ and Ca2+ ions are present at larger concentrations in the extracellular fluids. The plasma membranes which separate the intracellular and extracellular fluid compartments are strongly hydrophobic and act as a barrier, preventing the translocation of ions by passive diffusion mechanisms. Instead, the active and passive transport of ions and other substrates is achieved through specific and highly efficient transmembrane proteins including ion channels and primary and secondary transporters, among others. The small inorganic ions are either the direct targets of transport of these proteins or cotransported with larger organic substrates such as glucose, neurotransmitters, amino acids and so on. The orchestra of these proteins and transporters symphonized with different ionic types and substrates drives and sustains essential functions of cells and cellular organelles including the mitochondria, receptors involved in the immune response[3, 4], endoplasmic reticulum, cardiac and ventricular cells, and fibroblasts. Understanding and controlling ion transport has been a focal point of extensive multidisciplinary scientific studies spanning decades since clinical failures in key functionalities at the molecular level are linked to devastating health issues, known as channelopathies[5]. Among them are the long QT syndrome, various types of epilepsy, cystic fibrosis, hyperkalemic periodic paralysis, neonatal diabetes mellitus, myasthenia gravis, osteoporosis, and many others[5]. Here we will use some of these key functionalities to highlight the progress, but also the drawbacks and the challenges for computational approaches that are widely used to elucidate the transport mechanism of various transmembrane channels and transporters. 1.2 Ion selectivity as a key functionality of many transmembrane proteins. The transport of ions through the cellular membrane is either active or passive, requiring either positive or zero external chemical potential, respectively. Passive transport takes place in ion channels (Fig.1) and porins where the ions translocate rapidly along their electrochemical gradients. Ions may be highly or partially selected, and/or rejected at specific sites within these transmembrane proteins. Ion-selective sites are often composed of highly-conserved chemical motifs (amino-acid residues) that sort ions with high-fidelity. The challenge (compared to various inorganic materials such as ion-exchanging resins or chelating agents) is to achieve selective binding and effective transport at the same time. Ion-selective sites and motifs are usually highly conserved across many families of ion channels, but also commonly adapted to 7

serve various roles in different cells. A selectivity mechanism is a result of the intricate interplay between several key thermodynamic and kinetic properties in (often) high electric-field regions of charged residues. This itself poses a great challenge for the universal understanding of the ion coordination and selectivity mechanisms which often rely on a conservation of multiple selective sites, coupling between ion binding and protein dynamics among other factors. Mapping out these properties from experimental studies in a variety of ion channels has proven difficult due to the complicated evaluation of enthalpic and entropic effects and short-range versus long-range electronic effects. These effects range from polarization, resulting from ion binding to the protein host, to partial charge transfer, to complicated coupling between bound ions and the dynamics of their hydration shells. Computational studies, while very useful for testing various ideas regarding selective transport processes, also add plenty of fuel to intense debates on particular transport mechanisms. For example, classical force-fields, implicitly accounting for polarization effects, have obvious challenges if applied to studies of multiple ions confined in a protein. The amount of polarization response, the stability of various ion arrangements in the selective motif and underestimations or overestimations of solvent shielding effects and coupling to membrane voltage have stirred up debates among researchers in the field of computational chemistry and biology, who often have differing opinions on how these effects should be calculated or even included into the full picture. Adding to the heat of this uneasy debate is the presence of protons impacting ionization states of coordinating residues. Therefore, appropriate computational approaches are needed to address these effects and to integrate them into an optimal and satisfactory mechanism. This requires a good understanding of the performance of various computational approaches when applied to ion transport, which is one of the aims of this review. 1.3 Coupled transport of ions and substrates as a key functionality of many transporters. The elusiveness of adequate sampling of entropic and electronic effects is often amplified in active transport via primary and secondary transporters, which frequently conduct more than one type of substrate and have multiple binding sites. Primary transporters, e.g., the abundant ion pumps[6], use the energy released during common chemical reactions, such as ATP hydrolysis to push a substrate against its concentration gradient. In secondary transport, this typical translocation is coupled to the thermodynamically favorable transport of another substrate along its electrochemical gradient. Both substrates can move in the same (symport) or opposite (antiport) directions. Active transport is generally slower and requires significant conformational changes (Fig. 2) where the protein adopts several distinguishable states. Coordination of ions/substrates to the active sites of the transporters triggers appropriate changes in the protein structures to initiate a transporting cycle. In the so-called alternate access mechanism, a transporter switches consecutively between conformations opened toward the outside and the inside of the cell, during which there might be several stable intermediary states (e.g. a substrateoccluded state)[6]. Some active sites are strongly selective and transport occurs at specific stoichiometry, which is indicative of complex allosteric effects both in the ion/substrate coordination and the resultant structural reorganization[7]. Importantly, small blockers can inhibit the transport in active and passive transporters. This generally hinders the ability of many computational approaches, which can sample the conformational space but cannot include allosteric or long-range effects.

8

1.4 All-atom simulations for ion transport experiments and vice versa. At the macroscopic level, the ion flow to and from the cells can be tracked with various experimental techniques based on spectroscopy, electrophysiology, and mutagenesis, providing information on direction, stoichiometry, thermodynamics, kinetics and structural determinants of transport with a number of profound successes[8, 9]. For example, the atomic-resolution X-ray crystal structure of KcsA solved by the MacKinnon group[10] provided conclusive evidence for the existence of a narrow nanometer-size selectivity filter in the ion channels for the first time after more than thirty years of biochemical and biophysical work mapping possible structural arrangements[11-13]. However, static structures are often insufficient to draw atomistic-level conclusions on the complete mechanism of ion transport. It is obvious that selective binding and transport require substantial conformational changes in the selectivity filter region[14] often coupled with motions in the distant channel domains (e.g. voltage sensing domains in voltagegated channels, ion-induced gating in secondary transporters, etc.)[15-17]. Thus, computer modeling at all levels (from the atomic/molecular level to cell-level models) becomes crucial for the investigation of ion transport, selective binding and ion coordination mechanisms, not only as a complement of the existing experimental methods but also as a predictive tool that can be used to guide experiment[8]. The last 20 years have seen a significant increase of the computer hardware capabilities[18, 19], which in turn has led to the development of new computational methods[20-24] and to the movement of old and well established methods toward the treatment of larger, more complicated and more realistic models despite some pitfalls that will be discussed in this review. In parallel to the advance of computational power, just recently Cryo-Electron Microscope (Cryo-EM) techniques[25] for the determination of molecular structure have reached the stage where X-ray crystallography has dominated for years[26]. This in turn adds new opportunities and challenges for the computational studies since Cryo-EM techniques together with other techniques like solid-state NMR[27] and X-Ray[10, 28] crystallography can afford new and improved structures of transmembrane proteins which would be used as inputs in computational sampling and modeling, or serve as basic tests for computational models. 1.5 Why another review on ion binding to proteins? Our main goals with the review are many-fold. While the improvements in computational and experimental techniques for studies on ion binding to membrane proteins are impressive, there are many obvious challenges and hotly debated topics related to transport mechanisms. Accordingly, (1) We would like to summarize the exciting paths, on which debates have sometimes flared up, to reconsider, reconcile or challenge prior views on ion permeation, selectivity and coupling; (2) We strive to highlight computational hurdles that were fully and partially overcome to afford better agreement with experimental observations and pave the way to better understanding of ion-protein interactions in the contexts of various ion channels and transporters; (3) We also want to outline major computational methods and their domains of applicability in studies of ion-protein interactions underpinning the coordination and transport of various inorganic ions, which include Na+, K+, Ca2+, Mg2+, halides (Cl–, F–, Br–, I–), HCO3–, SO42– and inorganic phosphate Pi (PO43–, H2PO4–, HPO42–);

9

(4) We aim to arrive at a general overview of the capabilities of the current theoretical methods for ion transport; (5) We hope to facilitate the search for more detailed information on these topics, by providing references to several hundred articles of interest in an organized tabulated form; and finally, (6) We suggest a number of critical problems in ion-selective transmembrane proteins that require more attention in future research. With this review, we hope to aid the navigation of the vast sea of available information (especially in the cases of Na+ and K+ channels and transporters), and to draw attention to some less explored topics and research avenues that are rapidly developing with the advent of better experimental methods. 1.6 Organization of the review In Section 2 we recount the forty-year history of studies and debates on different coordination models in potassium channels to highlight key achievements, major computational difficulties and some on-going methodology developments. The choice of the K+ selective channels is not random. The field is, arguably, the most developed with first crystal structures, in-depth electrochemical, spectroscopic and biophysical studies in parallel with theoretical investigations at all levels of detail. Various techniques developed and tested first on K+ channels paved the way for applications to other membrane proteins. In Section 3 we examine how the achievements in studies of potassium channels have aided in the understanding of the transport mechanism of sodium channels, particularly for the selective coordination of sodium over calcium ions. In Section 4, we extend the discussions to the sodium coupled-secondary transporters, where significant conformational changes are essential for ion transport and the ion binding to a protein is coupled to transport events. In Section 5, we review a few available computational studies on calcium and magnesium transport systems, emphasizing the apparent lack of accurate computational treatments for divalent cations interacting with proteins. In Section 6, we introduce the limited research done on anion channels and transporters, which have been largely unexplored. In Section 7, we provide our own visions on the transport studies of the reviewed ions, sketch out research directions currently under development and highlight some largely unexplored areas. Finally, in the Appendix we provide some general information about the different computational approaches used in the studies of ion transport, and how to compute key observables for comparison with experiments. References to relevant books, reviews and papers have been made throughout the review for those readers who are interested to learn more. These references (and more) are organized in four tables in the Supporting Information (Tables S1-S4). Table S1 lists published reviews of experimental and theoretical (or both) results for channels and transporters of all listed ions. In the next three tables, additional relevant papers are grouped based on their major topic (Table S2), some common calculated properties (Table S3) and the used methods (Table S4). The EndNote bibliography file used for these tables is provided in the supplementary information along with an editable document file with the four tables for those readers who would like to continue updating the tables with new data.

10

2. Potassium Channels: Versatile platform for selective binding and permeation studies 2.1 Summary This section discusses several mechanistic models for ion coordination and transport, which have been proposed over the last forty years. These models are based on different sets of available data about crucial interactions in the ion channels. All-atom molecular dynamics (MD) simulations can sample these interactions to map out possible quantities that can give rise to selective coordination of ions and high throughput permeation. Some insightful thermodynamic quantities that can be obtained from MD simulations are the free-energy profiles that contain possible permeation pathways and enable computing observables (Kd, conductance, ion flux) for comparison to experimental data. Some major computational challenges for transport studies in potassium channels remain to this date. Na+ and K+ are the most abundant cations in the extracellular and intracellular fluids,

respectively (Table 1). They are involved in the regulation of cellular volume and osmotic pressure, acid-base balance, and the generation of membrane potentials[1]. Na+ is required in the function of various secondary transporters (vSGLT, LeuT, hDAT, etc.). K+ activates enzymes involved in glycolysis and the respiratory chain and is cotransported by some secondary transporters, similarly to Na+. Both cations are pivotal for the formation and transmission of electrical impulses (Action Potentials) playing roles in neuronal signaling and muscle contractions. Large amounts of Na+ and K+ are moved across the membrane via voltage-gated ion channels (Nav and Kv, respectively) following concentration gradients as a part of the constant membrane polarization/depolarization cycle in excitable membranes[29]. The low intracellular concentration of Na+ and the high intracellular concentration of K+ are maintained by the primary transporter Na,K-ATPase, which pumps Na+ and K+ against their concentration gradient. 2.2 Mechanistic models Several major developments over the last forty years have contributed to our modern understanding of currents and permeability in voltage-gated ion channels. Most famously, the voltage clamp experiment of Hodgkin and Huxley performed in a giant squid axon revealed that the action potential necessary for the conduction of a neural impulse is achieved by concerted movements of Na+ and K+ ions through the axon membrane and that the ionic currents could be related to the time- and voltage-dependent ion conductance according to Ohm’s law. Later, this understanding was expanded to include a relation between the membrane reverse potential (at zero current) and the ion permeabilities, known as the Goldman-Hodgkin-Katz (GHK) voltage equation[29]. The experimentally observed difference in the permeabilities of ions (based on the GHK equation) in some channels meant that these channels are more selective toward a specific ion. For instance, the voltage-gated (Kv) channels select K+ over the very similar (in terms of charge and size) Na+ ion in a 1000:1 ratio. Moreover, they transport the K+ ions at high rates, close to the diffusion limit[30]. It was inferred that the dehydrated K+ ions would likely traverse the selectivity filter (SF) in a single file alternating with water molecules, and at least 2 K+ ions could occupy the SF at a time[31]. The concerted movement of the K+ ions is the essence of the so-called knock-on mechanism, suggested by Hodgkin and Keynes[32] in 1955, in which bound K+ ions are dislodged from their binding sites by the electrostatic repulsion due to an incoming K+ ion. 11

Before 1980, when computational modeling of ion transport was still in its infancy and advanced crystallography techniques were still under development, biophysicists and physiologists had to rely on kinetic criteria to understand the ion permeation and selectivity, arising from the preferential coordination of ions at the protein active sites. The earliest models which tried to explain this selectivity include the rigid pore formulation of Mullins,[33] the fieldstrength model of Eisenman[34], and the snug-fit and selective-exclusion mechanisms of Bezanilla and Armstrong[35]. A pore model dating back to 1935, was first proposed by Daniellli and Davson[36]. Twenty-four years later, Mullins suggested that the cellular membrane must feature a cylindrical pore with diameter smaller than 10 Å in order to exhibit selectivity. Although Mullins could not imagine other kinds of ion-permeation pores, his arguments of hydration shells, hydration free energies, and the relationship with pore size bear the most similarity to the current realm of thoughts in the studies of ion coordination and permeation. Later, the field-strength model was developed to explain the large set of aqueous-glass electrodes that can give rise to different ion selectivity. It serves as the simplest physical model of the thermodynamics of ion selective coordination, simply by considering pKa values in an empirical relationship with an "equivalent-radius" that yields the same electrostatic energy of interaction for each cation. The field-strength model also describes relations between free-energy differences by exchanging ions and corresponding ion selectivity for various glass electrodes[34]. Notably, these ideas were explored further by computational modelling of small ionophores providing critical guidance to experimental studies[11]. In 1972, Bezanilla and Armstrong proposed a more specific pore structure of ion channels with a diameter in the range of 2.6 to 3.0 Å, approximately. Such a binding site would coordinate well the K+ ions providing a high-field environment but would be far too big for optimal coordination of Na+. It is clear now that although these models have some flaws, they captured many crucial quantitative features that are currently considered important for the description of ion selective coordination and permeation in a myriad of ion channels. More detailed reviews offering a historical perspective on the selectivity in voltage-gated ion channels can be found in Refs.[30, 37, 38]. The first X-ray crystal structure of a bacterial potassium channel, namely KcsA, homologous to vertebrate voltage-gated potassium channels was resolved in 1998 and has revolutionized research on excitable membranes[10]. This first structure, albeit of low resolution, provided unprecedented detail about the coordination of K+ ions by the KcsA protein. More high-resolution crystal structures of KcsA and other K+ channels resolved in the following years[28, 39-42] confirmed that the Kv channels consist of four identical alpha subunits with six transmembrane sections (TMs) each. The first four TMs form the voltage sensing domain which responds to changes in the membrane potential[43]. The pore of the ion channel is made of the remaining two TMs and the long loop (P-loop) between them. In the 3D structure of the Kv channels, the P-loops of the 4 subunits form a long and narrow selectivity filter (SF) within the pore region (Fig.1). A highly preserved amino acid sequence (T/S)XXTXGYG (where X is an amino acid) is found in the P-loop and is responsible for the highly selective coordination of K+ ions in these channels[43, 44]. The K+ ions in the SF are coordinated via the carbonyl groups of the GYGXT (X is Val in KcsA) residues in 5 binding sites (S0 - S4, Fig.1). The binding sites positioned relatively deep within the SF (S1-S3) provide 8 carbonyl groups for the K+ coordination, as demonstrated in Fig.3 for the S2 site. The crystal structures offer indisputable evidence that potassium channels have narrow pores that allow ions to permeate, supporting the snug-fit model, but apparently discarding the field-strength hypothesis for the selective coordination of ions in potassium channels. The crystal 12

structures are also consistent with the rigid pore formulation for the selectivity filter and the high throughput permeation of more than a million ions per second. Later on, however, many experimental measurements of the B-factor[45], which is related to root mean square deviation (RMSD) of atoms in crystal structures, determined that the atoms in the selectivity filter of KcsA are often fluctuating at the order of 0.8 Å. Moreover, the fluctuation of X-ray structures was measured at low temperatures, while the ion channel function was measured at room temperatures, suggesting that actual fluctuations may be even larger. The thermal fluctuations, taken together with the difference of only 0.38 Å between the ionic radii of K+ and Na+, shattered the certainty that the snug-fit model can explain how the selective filter in KcsA would coordinate potassium over sodium ions. Early calculations based on the Poisson-Boltzmann formalism for the KcsA structure done by Roux and MacKinnon[46] provided clear evidence for how the electric field of the SF stabilizes K+ with respect to the water cavity and the bulk water, indicating that the cage-like structure of carbonyl oxygen atoms (Fig.3) is a better home for potassium than sodium, thus backing the snug-fit model. Quantitative results based on more advanced molecular dynamics simulations of many crystal structures provided deeper insight into how the large fluctuations or “liquid-like” structure of the selectivity filter in KcsA may actually help in the selective coordination of potassium over sodium ions[30, 47-49]. Evidently, the evaluated Root-MeanSquare fluctuations of the selectivity filter atoms are often larger than 10-15% of the nearest neighbor distance, satisfying the definition of a liquid-like state. The novelty of this finding was that selective ion coordination can naturally arise in reduced models of four to nine oxygen ligands, zero to eight N-methylacetamide, methanol or water molecules, zero to four di-glycines, and zero to two neutral or ionized acetic acids[50], all of which can have relatively large fluctuations in their structural properties[50]. Bostick and Brooks[51] proposed a mechanism of topological control for the selection of K+ ions that put a strong emphasis on the number (=8) of ligands coordinating the ions rather than the chemical type of the ligands as long as the constraints are kept the same for the coordination number of complexes. For example, 8 water molecules and 8 carbonyl ligands with the same applied constraints or forces would yield the same selectivity for potassium over sodium. This model, however, was disputed by calculations in reduced models with 8-coordination states showing that the free-energy difference of about 12 kcal/mol is too small to account for various high ion-selectivity ratios in a vastly large number of potassium-selective channels[48]. By examining the free-energy profiles of all binding sites in the SF of KcsA and computing the alchemical free-energy changes between K+ and Na+, Noskov et. al.[49, 52] found that site S2 (Fig.3) is the most selective with the largest free-energy difference of 4-5 kcal/mol, suggesting a model termed selective-binding. This conclusion is backed by mutagenesis experiments[53, 54], which can switch on and off the selective ion coordination in the potassium channel by mutating site S2. The mutation alters the cage-like structure of S2 into a small water cavity that can coordinate equally K+ and Na+ and the selectivity is lost. This seems like convincing evidence for the S2 selectivity; however, subsequent computations showed that when four sodium ions coordinate at the planes of carbonyl-oxygen atoms in the wildtype selectivity filter in the conducting state of the channel, this binding is almost energetically equal to the binding of multiple potassium ions at the cage-like sites. This implies that the S2 site in the conductive state stops being selective in the presence of multiple cations[55]. Kim and Allen showed that perhaps selective exclusion might happen before the ions reach the selectivity filter to prevent multiple-ion accumulation of sodium in the selectivity filter[55]. By combining MD 13

simulations, electrophysiological experiments and X-ray crystallography, Thompson et. al.[56] showed that the two smaller ions, Na+ and Li+, can block the outward current of K+. As a result, for the selection of K+ over Na+, the entrance of the SF might have a large energy barrier for Na+ entry in the presence of K+, thus supporting the selective exclusion mechanism. Many simulations also demonstrated that in the absence of potassium ions, the selectivity filter would collapse, and that a single potassium ion can stabilize the selectivity filter and then prepare subsequent binding sites[57, 58]. In contrast, a single sodium ion may further destabilize the selectivity, and then hinder the overall ionic currents. Solid-state NMR experiments[59] on KcsA and Kv1.3 support these distinguishable structural rearrangement responses to different ionic types from MD simulations. If sodium ions somehow bind in the SF, they might cause an effect called "punched-through" when high voltages are applied in experiments[60]. These results suggest that a model where multiple Na+ ions enter the SF to form a permeation pathway similar to that of K+ ions is unlikely to occur. Despite all the work done on this topic, the mechanism of the ion selectivity appears unsettled. The permeation pathway of potassium ions in KcsA was once again re-considered recently[61]. The concerted, and widely accepted, knock-on mechanism was disputed based on a different interpretation of electron densities, which suggests that there might be no water molecules sandwiching K+ ions in the SF, thus enabling more efficient permeation pathways. Groot and colleagues reported in their simulations a major portion of pathways with bare potassium ions in the SF within a range of applied physiological voltages in both prokaryotic and eukaryotic potassium channels (KcsA, Kv1.2, Kv2.1) using both CHARMM[62] and AMBER99sb[63] force fields. This became known as the hard-knock model of potassium transport. Although this conclusion appears to contradict the data collected over the course of forty years that seems to support the knock-on mechanism[29], it is solidly based on a large set of simulations that model the high throughput of potassium ions across the SF. To explain such a high throughput in potassium channels, Vaziri and Plenio proposed that there might be shortlived quantum coherence that leads to resonance in transport efficiency[64]. However, this quantum picture of the ion permeation remains untested. 2.3. Ligands, solvents and ions in the confinement of K+ channel selectivity filter Local electrostatic interactions are supposedly the major determinant of the ion selectivity and permeation. These include ion-ion, ion-water, carbonyl-oxygen-ion, ligandligand, and ligand-water interactions, whose magnitudes and dynamic interplays are decisive for the observed selective ion coordination and varying permeation rates. Some of these interactions and their energetics will be addressed in more detail in the following sections. The work done in the last forty years has shown that quantification of the contribution of these interactions to the ion selectivity and permeation in the dynamical sense is not an easy task, even for classical approaches, such as MD. Solving Poisson-Boltzmann equations can yield a simple mechanistic model for how the SF favors coordination of K+ over Na+, but it raises the question of how entropic effects can alter this model[46]. The repulsive interactions between the carbonyl oxygen atoms within the SF are crucial for the maintenance of high selectivity as well[30]. The fluctuations of these carbonyl groups also play an essential role of either stabilizing the coordination of potassium ions and the SF or destabilizing the SF in the absence of K+. This effect depends on multiple-ion binding configurations[55] and the shielding properties of water molecules[47, 48, 65-70]. Water molecules are observed to move behind the SF wall at site S2 to play a critical role in the inactivation processes[16]. Upon mutations in S2, the ion selectivity is 14

lost (see above)[53, 54]. These mutagenesis experiments demonstrate that perhaps the ion channel needs exactly four clear cage-like binding sites (Fig.3) in order to coordinate potassium over sodium, and that there is a correlation between the binding sites. On one hand, it remains inconclusive how this correlation is measured from simulations and experiments, and what differences would be observed if four sodium ions were bound to the selectivity filter; the major culprit behind the difficult resolution of this task is the quality of the MD force fields for sampling. On the other hand, the structure and sequence of KcsA is highly conserved in many families of potassium-selective ion channels[29]. Different types of channels have distinguishable amino acids that are far from the SF, but still contribute to a large range of selectivity ratios and conductance values observed for different ions[29]. These facts raise serious concerns about the capability of the MD simulations to model such distant or long-range effects within affordable timescales. It is evident that the selective coordination is relatively sensitive with respect to different parameters such as van der Waals interactions (NBFIX)[47] and CMAP corrections to torsional interactions[71]. It is clear that the classical force fields may not yield an accurate description for the local electrostatic effects in ion-protein complexes such as induced dipole moments, and possible charge-transfer between divalent ions and negatively charged amino acids[72, 73]. How these effects are involved in the ion permeation and selectivity in other classes of ion channels, discussed in this review, is still unclear. Another critical issue is related to accurate models for describing water dynamics in the protein binding pocket. Water molecules are attracted from the bulk by the polar environment of cation binding sites. The local dynamics, induced dipole moments, hence screening properties of water relevant to ion-ion interactions should depend strongly on the details of the confined environments[74]. To illustrate the relevance of water degrees of freedom and their role in shielding ion-ion interactions, we examine some experimental and computational studies performed for interface systems[74]. In bulk water, both simulations and Ultrafast Infrared (uIR) vibrational echo experiments showed that local fluctuations of hydrogen bonds (in terms of acceptor-donor length) and angles drive the hydrogen bond dynamics on a timescale less than 400 fs, while at larger timescales water networks via hydrogen bonds are disrupted and rearranged. The dynamics of water re-orientations required for shielding effects and the formation or breaking of hydrogen bonds is described as a 60-degree jumping processes. The fluctuations in net dipole of these complicated structural arrangements are related to the large dielectric constant (shielding properties) of water. Therefore, hindrance of water dynamics in the confinement is expected to impact water shielding considerably. For example, for water molecules caught in the confinement created by surfactant [sodium bis(2-ethylhexyl) sulfosuccinate] micelle structures with intermediate and small sizes, the correlation time of hydrogen-bonding networks at the interfaces of the surfactant layers was almost 8 times slower than that in the bulk phase. Of course, the structure of the selectivity filter, in particular, and the entire ion-conducting pore, in general, is markedly different from the surfactant micelles, but the dynamics of water along the conducting pore in closed, transition, and open conformations should not be the same as in the bulk. In one of the simplest cation-selective channels, Gramicidin A, the structure of water molecules was found to be more flexible (unpublished data from SYN lab) than rigid models of water (TIP3P[75] or SPC/E[76]) if we consider polarizabilities of the media explicitly. Furthermore, we may expect ion-dependent water dynamics and hence fine tuning of solvent shielding to a specific ion. This territory is largely unexplored in studies of ion permeation and binding to protein confinements. The most widely used water models in classical MD or MC simulations were originally developed for studies of bulk water reproducing well-defined radial distribution functions, thermodynamics and the 15

diffusion constant. They speed up calculations and shorten the sampling timescales. However, these models prove deficient in producing the expected structural arrangements of water molecules in binary mixtures or ice phases[77, 78]. The quality of solvent models might be critical for resolution of an ongoing debate regarding selective permeation mechanisms such as the knock-on and hard-knock models. The knock-on model suggests that there should be a water molecule shielding the strong repulsive interaction between any two bare K+ ions bound in successive sites. In sharp contrast, the hard-knock model proposed by Groot and colleagues, who used both TIP3P and SPC/E water models in their calculations, suggests that at electrophysiological applied voltages such a shielding water molecule might be absent during high throughput ion permeation. What remains to be established is how both knock-on and hardknock models would be described by polarizable force fields applied to all atoms, which are believed to capture the local and distant electrostatic effects[72, 73]. If the hard-knock model has some merit, then it remains unanswered how it would help in the selection of potassium over sodium, and how water molecules would behave around Na+ in the SF. This task might need thorough quantum mechanical calculations[79]. It is important to note that quantum mechanical (QM) calculations for those interactions must be done and analyzed with care to account for the differences between reduced models and realistic systems. Very often, reduced models use representative ligands found in the ion-channel structures. As systematically examined in CHARMM22 calculations by Roux, out of more than one thousand reduced models, 39% resulted in nontrivial selectivity of K+/Na+ in the confined micro-droplet limit with free-energy differences larger than 2.0 kcal/mol[50]. These models are surprisingly useful in delineation of the two major conditions supporting ion selectivity: substructures around ions and the remaining architectural forces from the surrounding complexes. The free-energy decompositions for these conditions can be used to quantify the differences between proposed ion selectivity mechanisms, and to confirm that the large flexibility of the SF does indeed help to select ions. So far, most substructures exclude the water cavity, which might play a critical role in the overall permeation and selectivity due to its favorable affinity for K+ over Na+[80, 81]. One can systematically build a more sophisticated reduced model for QM calculations using minimum energy structures as well[82]. There appears to be a consensus on the effects of local sub-structures on the selective binding, but a disagreement about what is a "correct" architectural force, which may either come from artificial restraints to keep the quantum systems well defined or actually originate from the correct sampling of the surrounding effects[83, 84]. Either way, the message from using QM calculations with reduced models is clear: architectural forces from the surrounding protein in full systems must be properly taken into account to enforce a precise geometry for the atoms in the substructures. In other words, simply truncating crystal structures without proper sampling of surrounding protein scaffolds for QM calculations may not represent the full mechanistic picture required for understanding of ion selectivity mechanisms. As will be discussed for other channels, most of the QM calculations appear not to include the architectural forces due to the protein matrix outside of the binding pocket. 2.4 Challenges in coupling ion binding and conformational dynamics Proteins, as many other polymers, are highly dynamic molecules. Ion binding impacts protein structure on sub-Ångstrom and often Ångstrom scales[14]. The impacts could potentially be long-ranging and hence it is crucial to describe the vast conformational space sampled by the protein in response to an ion binding[16, 67]. This limits direct applications of QM simulations 16

to analysis of separate frames, often extracted from classical MD simulations. Most of the referenced studies in Tables S2-S4 are based on all-atom MD simulations. Specific examples of all-atom MD simulation results are given in all following sections. Routine information which can be extracted from such MD trajectories includes RMSD values as a measure of mobility and structure deformation, radial distribution functions as a measure of pair probabilities as a function of interatomic distances, dwelling times for ions as a measure of attraction/repulsion to certain portions of the protein, electrostatic interaction energies and free energies of binding, as a measure of binding strength, potentials of mean force (PMF) profiles as a fingerprint of the binding sites and measure of permeability, etc. PMFs and free energies of binding have contributed a lot to our understanding of ion transport systems and have merited their own section (see Thermodynamics below). Occasionally, MD simulations are used for evaluations and predictions of motions of interest. All-atom MD has been used also in the studies of gating in Kv and Nav ion channels[19, 85-91]. Kv channels, although very similar in their selectivity filter, can feature different modes of gating: voltage-, pH- and Ca2+-induced[38]. The concerted movement of positively charged residues in the voltage-gated domains in Kv arises as a response to the change in the membrane potential as suggested by many MD simulations[19, 85-91]. Similar processes were simulated during gating in Nav channels[92-97]. Global charge rearrangement is observed in recent long MD simulations on Kv channels which also afforded glimpses of the activation and inactivation processes in the Kv channel during gating[19, 88, 89]. The structural reorganization in the voltage sensing domains in this case lead to dewetting and collapse of the hydrophobic cavity preceding the selectivity filter[19, 88]. Other inactivation mechanisms were also probed with MD. For instance, the experimentally observed inactivation of the Kv channels at low K+ concentrations was attributed to distortions in the selectivity filter arising from the absence of K+ ions[98-101]. 2.5. Computational models for thermodynamics and kinetics of selective ion binding The MD sampling techniques aim to sample the complete conformational space relevant to the process. If complete sampling is achieved in an unbiased MD simulations one can assess all probabilities of states and construct unperturbed free-energy profiles for the process of interest.[102] The position-dependent free energy metric of the process is often referred to as the Potential of Mean Forces. Although long all-atom MD simulations (hundreds of µs to a few ms) are now feasible with the help of specially designed computer clusters like ANTON[18], only a few computational teams have direct access to such resources for limited applications. Therefore, accelerated sampling techniques and alchemical free energy calculations (umbrella sampling (US)[103, 104], thermodynamic integration (TI)[102, 105], steered molecular dynamics (SMD)[106-109], adaptive biasing force (ABF)[110], Metadynamics[111], Replica Exchange[112], see Appendix), which were discussed elsewhere for ion channels[113-115], would still be the major approaches in the near future. In essence, these techniques allow observation of relevant but rare processes in distant parts of the phase space via biases introduced in much shorter times than in straightforward MD simulations, enabling direct computations of PMFs. PMFs for ion transport/binding contain rich information about the freeenergy barriers, stable binding sites and possible entry and minimum free-energy routes for ions to binding sites in proteins[116]. The obvious challenge in ion binding or transport by a membrane protein is the unambiguous definition of reaction coordinates for describing the process. In the most popular approach, the PMFs of a single ion or multiple ions traversing the channel pore are routinely generated along the Z-axis orthogonal to the membrane plane (Fig.1, 17

Fig.4). This allows one to project out the thermodynamic pathways with free-energy barriers, transition points, and minima on the tractable reaction coordinate (Fig.4). The positions of these free-energy minima indicate possible binding sites in the protein even when such binding sites are unresolved in the X-Ray structures, while free-energy barriers and transition points are important to estimate transition rates. Often, ion binding/transport by ion channels relies on multi-ionic occupancies. In this case, one can use a set of biasing potentials assigned to the individual ions or to group of ions. 1D PMFs for a single ion can be generated from multi-ion PMFs by integrating out the degrees of freedom of the remaining ions (Fig.4, right). The key assumption justifying integrating out some of the degrees of freedom is the uncoupling of two (or more) reaction coordinates and hence allowing projection on effectively 1D reaction coordinates. We were unable to find any methodical studies addressing whether or not this assumption is valid in various ion transport studies. By reducing from 2D to 1D, our new “effective” 1D PMF includes implicitly ion-ion interactions and differs considerably from a standard PMF from simulations done for a single ion permeation case. For instance, single ion PMFs through K+ and Na+ channels often yield energies, which are too low and signify overestimated ion binding. 1D PMFs integrated from multi-ion calculations correct for this problem since they include the ion-ion repulsion energy, which can be as large as 71×(4.5/εr) kcal/mol with r being the distance between the two ions and ε representing the local shielding effect by water molecules and carbonyl oxygen atoms (Fig.4, right). This repulsive interaction can reduce noticeably the free-energy barriers between the ions binding sites in the selectivity filter and propels them to move on following either the knock-on mechanism [55, 65] or the hard-knock mechanism. What the correct local shielding effect should be is, however, nontrivial. Setting aside the issue of accuracy in classical force-fields, it is instrumental to review progress in PMF computations for a variety of model channels. One- and multi-dimensional PMFs have been generated for both potassium[55, 65, 117-119] and sodium[120-124] ion channels. In the case of K+ channels, the PMFs indicated concerted motion of two or three K+ ions through the selectivity filter and a knock-on mode of movement[65, 118]. The 3D PMFs from US identified site S2 (Fig.3) as K+ selective with a 6.6 kcal/mol barrier, in accordance with FEP calculations and Ba2+ blockage experiments[65]. Alternatively, multi-ion PMF studies in the non-selective NaK channel revealed that the lack of selectivity is due to structural changes at the S2 position leading to a larger cavity, which could accommodate hydrated ions (see also discussion on S2 above)[125-127]. The role of the coordination number of K+ and Na+ and in particular the change in coordination from water to the carbonyl groups in the NaK pore was also elucidated in these studies, since MD allows for tracking of ion-water and ion-protein coordination during the ion translocation. The selectivity at the S2 binding site was discussed in terms of the field-strength model, since mutations in NaK reproducing the selectivity filter of K+ channels ultimately lead to dehydration of K+ and its coordination by a rigid carbonyl scaffold[125-127]. A comparative study determined that US, metadynamics and SMD all afford PMFs with similar positions of the binding sites. However, the energy barriers evaluated from SMD were larger than the ones obtained from US[117]. Permeation of Na+ ions in the K+ channel selectivity filter was also examined with multi-ion PMFs and uncovered Na+-specific binding sites indicating that the selectivity in K+ channels may be a result of kinetic factors (i.e. energy barriers) rather than the thermodynamics of binding[55, 119].

18

It is important to note that in the discussion of ion channels, the definitions of selectivity and hence discussions comparing computed and measured properties can be cumbersome[30]. Experimentally, a channel is deemed selective for a certain ion because it preferentially allows the permeation (thus, having higher experimentally measured conductance) of this ion. These experiments involve measurements of reversal potentials that cancel out currents caused by ionic gradients across membranes. Therefore, in this case selectivity is related to the ion permeability depending on both thermodynamic and kinetic properties of the whole channel. That is, during the trans-membrane translocation, ions have to overcome free-energy minima and barriers, which affect the overall permeability often observed indirectly through the signatures of ionic currents and conductance. On the other hand, computed thermodynamic selectivity can be discussed simply as the preferential coordination or the binding affinity of a specific ion with respect to a specific binding site, which is dictated by the free-energy loss upon binding. In this sense, selectivity can be related to the depths of the minima in a free energy surface or to freeenergy penalties when a favored ion in a binding site is replaced with a less favored one. The different free-energy profiles can be integrated to estimate the changes of equilibrium dissociation/association constants for quantitative comparison with experiments (see Appendix). In K+ channels, K+ selectivity was demonstrated from Ba2+ blockage and punch-through experiments[60, 128, 129]. The FEP and TI methods in ion channels are based on stepwise transformation of an ion into another one and evaluation of the relative free energy change ∆∆G from eq.A5 (see Appendix). Such calculations were performed for all binding sites in KcsA and determined that the most selective site was S2 (see above)[47, 65]. Structural changes in the S2 active site of the channel NaK (permeable by both Na and K), have been related to the loss of K+ selective coordination in NaK[127]. FEP calculations in S2 with different force fields yielded similar values of ~6 kcal/mol, in accordance with Ba2+ blockage experiments[127]. The origin of the KcsA selectivity was further elucidated by FEP studies with small “toy-model” systems consisting of 1 ion and 8 freely fluctuating carbonyl groups, which provided similar ∆∆G values (~6 kcal/mol) even in the absence of the rigid protein matrix.86 The selectivity for K+ was lost when the number of coordinating carbonyl ligands was reduced to 6. Further analysis of the enthalpy and entropy components of ∆∆G revealed that selectivity was enthalpy driven and that the enthalpy had two major contributions, which favored smaller and larger ions respectively: the ion-ligand attraction and the ligand-ligand repulsion. This result added another dimension to the already present size, rigidity, field strength and hydration/dehydration formulations for ion selectivity in channels: the number and type of ligands and the ligand-ligand interactions (see above). It must be mentioned that just like standard and biased MD, FEP is sensitive to the length of the simulations. For instance, it was suggested that the Na+ ions actually have binding sites in the SF of potassium channels but these sites are different from S0 - S4 in Fig.1. The short FEP trajectories would not allow the Na+ ions to migrate to their preferred binding sites and would lead to overestimation of ∆∆G[55]. 2.6. Ion fluxes and conductance Another important piece of information that can be extracted from all-atom MD simulations are I-V profiles. I-V profiles have been calculated extensively also with coarsegrained computational models, such as Brownian dynamics and have been used for comparison to experimentally determined ion fluxes and conductance (see below). The use of MD for I-V profile evaluation is often referred to as “computational electrophysiology”. The effect of transmembrane potentials in MD simulations might be modeled either via ion concentration 19

gradients across the membrane or through the introduction of an external electric field as a driving force[89, 130, 131] (see Appendix). The former approach might lead to complications due to the periodic boundary conditions [38] that would practically seal both sides of the membrane, while the latter approach may lead to uncertainty in the experimental mapping of the corresponding applied voltage since either a thickness of lipid bilayers or a simulation dimension can be used[131]. One may want to fix these problems by introducing double layers of an identical system that would separate compartments with different salt concentrations, thus maintaining asymmetrical physiological conditions in ion channels[132-134]. Ionic currents from charge imbalances between the compartments can be sustained by exchanging ions from one compartment with water molecules from the other compartments. This exchanging scheme might introduce artificial Ohmic resistance, which was believed to be negligible. The resulting transmembrane potentials are then estimated through the charge-capacitance relation. Although this approach may double the system size of an original transmembrane protein, it is integrated in Gromacs with fairly inexpensive computational overhead. This computational electrophysiology approach has been recently used to model the high permeation rate of K+ ions in KcsA and has led to the proposal of the hard-knock mechanism highlighted above. In this study calculations were carried out for the simulation of ionic currents in the KcsA channel in both openinactivated and closed-conducting states. Both CHARMM and AMBER99sb force fields together with TIP3P, TIP4P and SPC/E models of water have been used. The ionic current was computed from the imbalanced number of ions between the separated compartments. Several simulations starting with different initial states (KKKKK or KWKWKW) and different saltconcentrations were run for at least one microsecond for determination of the number of events with the hard-knock and knock-on configurations. The resulting conductance was 80 pS, which is in excellent agreement with experiments. These simulations showed that at the simulated transmembrane potential of 90 ± 50 mV, the salt concentration of 100 mM, with SPC/E and AMBER99sb force field[135] there are 30 permeation events, while at 560 ± 120 mV, salt concentration of 400 mM, TIP3P, CHARMM force field there are 29 permeation events. Noticeably, the CHARMM force field used in the same condition as the AMBER99sb force field appears to not capture the same outcomes. These simulations contain more hard-knock events than knock-on events. It is interesting that this computational electrophysiology approach can produce excellent conductance values and challenge electrophysiological data and IR data that support knock-on configurations[67]. At the same time, it exposes many computational challenges that have been present for a long time in the field. Although the PMFs are not directly available from experiments, the properties estimated from PMFs such as ion fluxes, selectivity ratios and overall conductance could be compared with experimental data. There are multiple ways to compute such quantities discussed by Roux et. al.[136, 137] It must be noted that the ion selectivity and permeation arise under applied voltages, although most of the PMF computations exclude the effects of applied forces. The computation of PMFs assumes that all degrees of freedom relevant to the process are properly captured. Therefore, an implicit assumption in connecting inherently non-equilibrium properties measured in experiment, e.g. conductance, to equilibrium PMFs is that the free energy surface remains unchanged across various applied voltages and time-scales of the permeation process. The biasing force due to electrical field could be added post-factum using free energy profiles from PMF computations as input to various methods such as the Brownian Dynamics (BD) and Poisson-Nernst-Planck (PNP) methods (see Appendix) [138-141]. This assumption appears to provide satisfactory correlations to experimental measurements in the potassium channel. It is 20

possible that the simulated structures are the most probable/dominant in the processes of ion selectivity and permeation. Even despite the implicit inclusion of the protein matrix, BD and PNP manage to provide some important structure-activity information, such as ion-channel interactions[142, 143], stability of protonated states[144], potential ion pathways within the pore[143, 145, 146] or effect of the pore radius on the ion fluxes. In working with these models, where the true crystal structure was used instead of an idealized channel, possible binding sites could be identified from PMFs and the ion distribution in the pore[69, 147-149]. Some of the earliest theoretical studies on K+ and Na+ transport in ion channels[150-152] and porins employed continuum electrostatic approaches like PNP and PB or coarse-grained BD simulations[38, 136, 137, 153]. The popularity of these methods was a result of several factors: (1) They were already well established for evaluation of ion fluxes in nanopores and have been applied to model channels like Gramicidin A. (2) They were computationally affordable and allowed simulations at the ms time scales even with the old and limited (from present point of view) computer resources. (3) The lack of crystal structures hindered the modeling of the structural intricacies of the ion channels and generalized implicit models of the protein were often used instead. These attractive qualities of the continuum models and BD make them applicable today – they are still in use in studies of ion transport. To account for membrane potentials and protein dynamics, BD calculations were coupled with Grand Canonical Monte Carlo procedures and incorporated direct account of flexibility of protein fragments[138, 154]. The results were in good qualitative (e.g. selectivity series were reproduced) and even semiquantitative agreement with experiment although routine overestimation with 30-50% is still observed in conductance calculated with PB and PNP[137]. Other systems which were studied with continuum electrostatics methods and BD are the OmpF porin[146] (Fig.5, see below), the α-hemolysin pore[138, 143], model Cav channels[155, 156], and Cl- channels of the CLC family[153, 157]. OmpF and α-hemolysin have a wide modestly selective pore and are a better target for such studies than the narrow, very selective Kv and Cav channels. The performance of these methods has been evaluated also against all-atom MD simulations and experiments and it was determined that BD performs better than PB and PNP. More information on the history, concepts, and applicability of these methods can be found in Refs.[38, 137]. 2.7. Challenges and unresolved issues in computational modelling of selective ion transport across K+ channels One of the structurally simplest cation channels that has been used in MD simulations attempting to reproduce experimental ion fluxes, ion selectivity and conductance is the Gramicidin A channel. Gramicidin A (gA) is an antibiotic forming a narrow channel[4] that transports monovalent cations, is blocked by bivalent cations and does not transport anions. The pore of Gramicidin A is formed by two dimers bound head to head with hydrogen bonds (Fig.5). It is lined with both carbonyl and amide groups and as such would be expected to transport both cations and anions. Attempts at explaining the curious selectivity of gA include protein flexibility and larger mobility of the carbonyl groups, which can reorganize to accommodate cations[158, 159], higher desolvation penalties for the anions, and charge asymmetry leading to larger negative charge at the centre of the pore[137, 160]. Gramicidin was the first ion channel to be simulated with MD[161] and since then has been a target of numerous studies of selectivity and permeation based on almost all known methods in the field[38]. Due to its small size and the abundance of experimental data it is often employed as a test system for newly developed 21

methods[162]. Historical perspective on gA simulations can be found in Ref. [160]. It is also naturally selective for K+ over Na+ providing, arguably, the simplest test ground. The deficiencies of the standard force fields have already been discussed and have contributed also to significant discrepancies between theoretical studies and experimental results on Gramicidin A transport mechanisms[163, 164]. A major conceptual deficiency of MD is its reliance on force fields, which are, by their definition, imperfect. They simplify immensely the real energetics of the system (e.g. excluding the electrons), assume that atoms of the same type would behave the same way even if their surroundings are different and are parameterized for globular proteins in short simulations mostly from calculations in water, which has a much higher dielectric constant than lipid membranes and proteins. In other words, these parameters may not be transferable from bulk water to the diverse protein environments because complicated polarization effects, charge transfer or interactions at the border between two phases (water and protein) are not properly evaluated. The use of partial charges for all atoms in the classical forces may simplify too much the electron densities in high-field environments, such as the SF in ion channels[165, 166]. The assumption of pairwise interactions with additive functions may not hold in high-field and highly correlated environments as well. Density functional theory (DFT) and ab initio MD calculations on models of K+ channels have determined that classical force fields fail to reproduce important protonation effects in the binding sites and tend to overestimate the coordination numbers of the ions both in the pore and in bulk[167-170]. Other QM-based studies point towards the importance of protein backbone polarization by K+ ions and methyl group polarization for selectivity and permeation[166, 171, 172]. Several polarizable force fields are available with Drude[20, 72], and atomic-multipoleoptimized-energetics-for-biomolecular-applications (AMOEBA)[23, 24, 173, 174] being the best known. Other polarizable force fields like Polarizable Charge Equilibration Force Fields[175] were applied to study the permeation of K+ in gA, but are not widely used. In the Drude polarizable force field, massless but charged particles are attached to each atom with a harmonic spring and can reorient with respect to the external electrostatic field. Drude is more computationally feasible and can be applied to larger systems than AMOEBA, which uses a set of multipoles assigned to each atom to reproduce the atomic polarizabilities. Although both methods require additional parameterization and are still under development, Drude has been applied in several recent studies for evaluation of charge transfer effects upon binding, which are especially prominent for polyvalent cations like Ca2+[72, 73, 176]. The agreement for Drude with higher quality QM results was better than the one observed for classical additive force fields. Based on our experience with the Drude force field, one might need a complete parameter set of all different atom pairs (e.g., NBFIX parameters) to account for different environments/interactions, while arithmetic combination rules used for Lennard-Jones interactions are more robust. The possible structural change and dynamics of water molecules and confined conducting pores may require thorough quantum treatment that needs to include architectural forces, multi-ions and applied electric fields. It is critical to address whether the correct “quantum” water molecules and induced dipoles would provide more evidence of the existence of the hard-knock model. This task requires a significant sampling of architectural forces from all-atom MD simulations that may be again force-field dependent. One way to overcome this might be to use an efficient quantum algorithm like density functional tight binding (DFTB), which is ~100 times faster than higher-level DFT theories (see Appendix). However, even a 100-times speed-up might not be sufficient for sampling enough entropic effects in large transmembrane proteins. Another possible choice is to use some QM/MM 22

scheme with the assumption that electronic effects diminish as the distance from the reactive region increases[22, 177, 178]. Nevertheless, the inclusion of the dynamic transport effects in QM/MM calculations would require a constant redefinition of the QM region and may not be feasible. Moreover, QM/MM schemes are still computationally expensive and will remain a major challenge in the studies of ion binding and subsequent transport by transmembrane proteins.

23

3. Sodium Channels: Charged Ligands, Water and Loose Coordination 3.1 Summary Here we discuss the mechanistic models, thermodynamics, essential interactions and allatom sampling for studies of ion selectivity and permeation proposed by QM and QM/MM computations and molecular simulations in sodium channels. This discussion suggests that sodium channels may be even more challenging for computational studies than potassium channels. 3.2 Mechanistic models and thermodynamics The existence of sodium and potassium channels was firmly established in the mid 1950s [32]. However, the first 3D structure of sodium channels was quite elusive until a breakthrough from the Catterall lab in 2014 when the first structure of a bacterial voltage-gated sodium channel, NavAb, was published[179]. The long trans-membrane sequence (alpha-subunit) of a canonical Na+ channel is divided into 4 domains of 6 TMs each and, as in the Kv channels, the SF belongs to the long P-loop between TMs 5 and 6 of each domain. The channel can be formed as a homo-tetramer[179] (in prokaryotic Nav) or as four repeats containing structural elements described above (in eukaryotic Nav). Interestingly, the SF of Nav contains amino acid residues with charged side-chains (Asp, Glu and Lys). In the homo-tetrameric prokaryotic Nav channel 4 Glu residues from the four domains form an EEEE ring in the selectivity filter (the so-called high-field binding site SHFS, Fig.5)[179]. The presence of the charged ring is in line with the idea of “high-field” ligands coordinating smaller ions and compensating for the dehydration penalty, but the reality of selective permeation turns out to me more complicated. The structure of NavAb has surprised many researchers with the size of its channel pore (Fig.5) when compared to the SF of KcsA: its diameter is about 8 Å, more than twice as large as the diameter of the SF of KcsA. The permeating ion is involved in a rather loose coordination by protein residues. At the same time, multi-ion permeation mechanisms appear to play a substantial role in Na+ transport and size-altering mutations (E/D) have a major impact on Na+/K+ selectivity of the Nav filter. At least three binding sites can be distinguished (SHFS, Scen, Sin, see Fig.5) in the crystal structures of NavAb and NavMs[179-184], although a recent study based on long MD simulations in NavMs identified five binding sites in the SF similarly to the K+ channels[185]. The E178D mutation reduces the number of binding sites to only two (as found in the crystal structure)[180]. Another bacterial NaChBac ortholog, NavRh, was also identified with only two binding sites in its SF[186]. Interestingly, in NavRh the ring of Ser180 residues is found at the position corresponding to the ring of EEEE in NavAb, which moves by one amino acid toward the extracellular side. It should be noted that prokaryotic Nav bears a resemblance to some voltagegated Cav channels that also feature an EEEE motif[187]. In the asymmetric eukaryotic Nav, the filter contains a DEKA ring[188]. Considering the larger conducting pore in prokaryotic and eukaryotic sodium channels, the snug-fit model proposed for potassium channels becomes inadequate for the explanation of the ion selectivity. If one recalls the Mullins relationship[33] between channel pore sizes and ion selectivity, i.e., larger pores (but < 10.0 Å) would select ions with larger dehydration energies, then NavAb appears to obey the Mullins rule because the dehydration energy of Na+ is larger than that of K+ by 17-20 kcal/mol[189]. However, NavAb binds and transports selectively Na+ over Ca2+, which has a much larger dehydration energy than Na+ (by about 250-280 kcal/mol)[189]. Similar problems arise when applying the field strength 24

model of Eisenman[34] to these channels. The SF pore of NavAb has a stronger field strength that favors sodium over potassium ions, but cannot be applied to calcium ions simply because Eisenmann did not consider the monovalent-vs-divalent selectivity in his model. Nevertheless, the Mullins and Eisenmann work provided a quantitative description, albeit simplified, for sodium-selective channels based on the strength of the ions or the coordinating ligands. These interdependences have been recently revealed in a number of QM calculations (see next Section). Various MD simulations appear to converge on key features for selective binding and transport mechanisms (see potassium channels discussion) for sodium channels. PMFs in NavAb reveal three favorable binding sites (SHFS, Scen and Sin), at positions consistent with the X-Ray structures (Figs.3-5)[121, 122, 124]. Concerted movement of two Na+ ions was also observed with an overall permeation barrier of ~3 kcal/mol[190]. A third Na+ could also be involved in the permeation process (average occupancy of the selectivity filter was 2.3 ions). Often, the selective binding mechanism implies that a SF has strong and stable coordination sites for one ionic type; while other ionic types are not energetically stable at the same sites. MD simulations by Corry and Thomas[121] suggested that a single potassium ion with the complete first hydration shell would not fit the pore of NavAb at the ring of four negatively charged glutamate residues, while a single sodium ion with the smaller ionic size and large dehydration energy would enable a sufficient space for water molecules to form a molecular bridge with the glutamate carboxyl groups. This conclusion was later confirmed by Furini and Domene[122]. In comparison with the knock-on permeation mechanism, Corry and Thomas also proposed that sodium ions would follow a similar but loosely coupled knock-on conduction, which Boiteux et. al.[191] showed to be more complex. The sodium ions are not just loosely coupled, but concertedly translocate, pass-by each other and correlate with noticeable movements of the glutamate side chains. Importantly, Pomes and colleagues[120] suggested that the permeation in NavAb may involve three to four sodium ions and that a significant number of times three ions are located within the SF of NavAb, which has the three distinct binding sites as identified from the crystal structure. These mechanistic models of ion selectivity and permeation can explain the modest selectivity with free-energy differences of 1-2 kcal/mol favoring sodium ions, and the possibility of potassium ion translocation in any direction within the sodium channels. FinolUrdaneta et. al.[192] suggested a distinct permeation pathway for potassium ions traversing the SF in which a K+ has to wait for another replacement ion before completing a translocation cycle. However, experiments[185, 192] reported that under bi-ionic conditions, potassium ions may block or hinder inward current of ions. Using MD simulations, Ngo et. al.[193] illustrated the possibility that the multiple-ion binding and relatively large pores with negatively charged residues at the extracellular mouth of the SF may give rise to asymmetric translocation for potassium ions in NavAb, resulting in a high-energetic blockage of ionic currents due to 2-3 potassium ions, which are forced to move in the selectivity filter at the same time in the inward direction. This result suggests a selective exclusion mechanism in the inward direction to enable energetically favorable permeation of sodium ions. All of the mechanistic models ascertain the important role of the glutamate residues in the SF. In fact, Finol-Urdaneta et. al.[192] showed that the E/D mutations in NavAb and NaChBac clearly alter the Na+/K+ selectivity and permeation. Those models, however, have not yet been tested in a similar fashion as has been done for the debate on the knock-on vs hard-knock permeation pathways, i.e., by using IR spectroscopy[67]. Much less effort was dedicated to study the Na+/Ca2+ selectivity in NavAb. Both Corry and Ke et. al.[194] showed by means of classical MD simulations that calcium ions could bind strongly to the SF and experience higher barriers than sodium ions during translocation[195]. In 25

the prokaryotic NavAb which is selective for Na+ and conducts very little Ca2+ despite the identical size of the two ions, high desolvation penalties for Ca2+ were discussed as a possible origin of Na+ selectivity based on simulated PMFs. NavAb can select sodium ions over calcium ions simply by allowing calcium bound intact in the SF, while sodium ions can easily pass-by. This is an attractive model for the Na+/Ca2+ selectivity because it may explain some experiments on NaChBac where small concentrations of calcium ions (<1mM) may not alter Na+ currents[196]. However, it sheds little light on the noticeable interference of Ca2+ with the conduction of Na+ in eukaryotic sodium channels, perhaps due to the unavailability of crystal structures of eukaryotic sodium channels. Therefore, the Na+/Ca2+ selectivity and permeation mechanism remains unclear in a large family of sodium channels. Of particular interest here is how multiple ions compete to selectively translocate through the high field-strength environments in sodium channels and some QM calculations (see below) provide meaningful insights into the selective binding. 3.3 Coordination and importance of polarization effects in the selectivity filter of Nav channels All of the key ligating interactions, discussed for potassium channels, are also applicable for description of the Na+/K+ selective coordination and permeation in sodium channels. However, the binding sites and filters we have seen so far are less selective, protein coordination is loose and very dynamic, bound cations retain partial hydration spheres and ion-protein charge interactions are of great importance. In a nutshell, we have a system where polarization effects are expected to play substantial roles along with water and protein dynamics. To support this point, recent QM calculations showed that the strength of ions originates from the ability of the atoms to accept electrons from coordinating negatively charged ligands, e.g., the electronacceptability decreases as Ca2+ > Na+ > K+[197]. The strengths and topologies of coordinating ligands in sodium channels may be more complicated and chemically diverse than in the SF of potassium channels, which contain only carbonyl groups. In both prokaryotic and eukaryotic sodium channels, the resultant strength of the coordinating ligands is a result of complicated interplays between hydroxyl, carbonyl, negatively charged carboxyl, and positively charged amine groups. The flexibility of the protein matrix around the permeation pathway, the packing of ion-coordinating groups are all decisive for the selective coordination of sodium over potassium and, to a lesser degree, calcium ions. QM based on DFT calculations of reduced models of the SF in voltage-gated ion channels have also been done and have described the ion selectivity in Nav as resulting from a complex balance between the number of coordinating ligands, the size and rigidity of the pore[197-200]. That brings forward an important question on how we could possibly construct idealized and yet informative models of the selectivity filters in eukaryotic Nav channels. The selectivity filters of eukaryotic channels with rings of DEKA/DKEA and EEKE/EKEE are sodium selective. When mutating the ring of DEKA to DERA preserving the net charge, the Na+/K+ selectivity is reduced significantly; upon mutation of the positively charged to a negatively charged residue, the eukaryotic channel becomes calcium-selective. Therefore, a minimal model should at least reproduce this feature observed in experimental studies. The presence of Lys in the coordination sphere of a cation is shown to decrease the cation coordination number to three, in contrast to protein coordination numbers found in potassium and calcium channels (eight ligands in KcsA, Fig.3, and ≥ four in calcium channels). This has been suggested to play a role in the favorable binding of sodium ions over others. However, crystal structures we obtained for a bacterial Nav channel show that water may 26

compensate for missing protein coordination efficiently. The preference of the EEEE binding site for Na+ ions arises from the polarization of the water molecules bridging the ion with the protein matrix which is stronger for Na+ than for K+ ions[200]. A plausible role of positively charged Lys in the selectivity filter for Na+ includes net reduction in the negative charge to particularly enhance the selectivity of Na+ over Ca2+ because of the smaller dehydration penalty and attenuated electronic effects. It is also important for the formation of a strong hydrogen-bond network with other negatively charged residues. This salt-bridging network may potentially rigidify the SF of eukaryotic channels with DEKA more than the SF with DKEA further enhancing the Na+/Ca2+ selectivity. For instance, the idealized and rigid DEKA motif used in QM computations suggests that DEKA coordination is more suitable for the Na+ ions than for K+ and Ca2+ ions. It must be noted however that these reduced models do not consider the multipleion effects and influences of applied electric fields, and do not take into account the architectural forces, which might affect the flexibility and rigidity of the SF. 3.4 Computational Electrophysiology Studies Sodium ion fluxes were simulated first via the computational electrophysiology technique [132, 134] (see Appendix) with the AMBER99sb force field and simulation set-up and parameters similar to those used by Kopfer et. al.[61] for studies of ion permeation across the KcsA channel. Although the applied transmembrane voltage of 565±126 mV is larger than physiological voltages, its purpose was the collection of sufficient statistics within a few microseconds. These simulations were able to describe the asymmetry of outward and inward ionic currents, during which the conformations of the glutamate side chains in the SF of openstate NavMs are statistically different. The estimated resulting inward conductance was 27 pS, while the outward conductance was 15 pS[201]. More calculations, similar to those performed in KcsA need to be done for assessment of the force field effects on the permeation pathways of sodium and potassium ions and for probing the flexibility of the side chains in the SF at physiological voltages (< 200 mV), when this flexibility may be subtler. Current all-atom sampling of the crucial interactions in sodium channels is still limited, particularly for the Na+/Ca2+ selectivity and permeation[202]. In physiological conditions, the conductance of Na+ diminishes by 20%-40% due to the interference by extracellular Ca2+ at the resting potential[203205]. In the limited number of classical force fields calculations, the results must be considered with care due to the poor calibration of the force fields for calcium ions (see below). As shown by MD simulations, the occupancy of sodium ions in the negatively charged selectivity filter of NavMs is often three and can go up to four ions[185], creating a high-field environment, where significant polarizable effects might be essential. Thus, it is critical to have more accurate sampling in terms of better polarizable force fields and more efficient QM calculations. While classical parameters for calcium ions might be helpful in some cases, they essentially overestimate the interactions with protein complexes[72, 73]. One of the present challenges in the computational field is the development of reliable parameters for the interactions of calcium ions with all amino acids. Another major challenge for sodium channels is the sequence deviation when going from prokaryotic to eukaryotic sodium-selective channels. All-atom sampling or modeling of eukaryotic sodium-selective channels has not been done, simply because there has been no reliable crystal structure. In this sense, the DFTB approach (see Appendix) might be a suitable computational technique due to its efficiency and ability to describe polarization and small charge-transfer effects. However, DFTB is still computationally inefficient for sampling of the large selectivity filters of many sodium channels with more than 27

5000 atoms. Such sampling is critical for the evaluation of the flexibility-vs-rigidity of the substructures, and for the inclusion of architectural forces by surrounding and distant amino acids. This large number of atoms might pose the same limitation to QM/MM approaches (see Appendix). Very often, in ion selectivity and permeation, there might be very little charge transfer but significant polarization effects, so a particularly promising approach would be using an efficient force field like the Drude polarizable force field. However, to our knowledge, the Drude force field is still under development: there are many missing parameters for interaction of calcium ions with amino-acid side chains, while the parameters for protein structures transferred from the additive CHARMM force field are still being tested. To this end, development of better force fields and more efficient semi-QM or QM approaches should be the main focus for improvement.

28

4. Sodium and Potassium Transporters: Computationally Demanding 4.1 Summary This section discusses some recent computer modelling of sodium coupled transporters to highlight computational difficulties in the selection of reaction coordinates and the description of allosteric effects in primary and secondary active ion transport. 4.2 Active transport of Na+ and K+ Active transport of Na+ and K+ through the membrane is achieved through a number of primary transporters (Na,K-pump, Na,H-pump, K,H-pump) and secondary transporters which couple the translocation of Na+ or K+ along their electrochemical gradients in order to move a substrate against its concentration gradient. A number of MD simulations of pumps transporting Na+ and K+ are available[206-208]. Other pumps like the sarcoplasmic Ca2+-ATPase have been studied more extensively and will be discussed in the section dedicated to Ca2+ transport. But the champions in number of studies are secondary Na+ coupled transporters which are the most extensively studied with all sorts of computational modelling. Like all secondary transporters, they undergo substantial conformational changes during transport (Fig.2) in response to ion and substrate binding. In the most influential model of the field - the alternate access model, the transporter molecules switch between conformations opened toward the extracellular space (outward-open) and conformations opened toward the intracellular space (inward-open)[209, 210]. Intermediate occluded and closed structures can also exist as part of the transport cycle. Most of the computational studies have been performed on two bacterial systems, LeuT and GltPh (Fig.6), both of which are homologues of transporters found in the human brain and critical for its function. LeuT is used as a model of transport of monoamine mediators[210, 211] (serotonin, dopamine and noradrenaline) and GABA in the central nervous system, while Glt Ph serves as a template for translocation of excitatory amino-acids (aspartate, glutamate) in the EAAT systems[209, 212]. A high-resolution crystal structure of LeuT was first resolved in 2005[213]. Since then more than 45 structures of LeuT with different bound substrates and conformational states have been determined[210]. The structure of GltPh was also resolved in several different conformations[214-216]. 4.3 Structure and functionality of secondary transporters (GltPh and LeuT): Where ions come to play LeuT is a homodimer transporter from the thermophilic bacteria Aquifex aeolicus, which cotransports 2 Na+ with 1 small neutral amino acid like Leu or Ala (Fig.6)[210, 213]. Each LeuT domain consists of 12 TMs and features at least two substrate binding sites - S1 and S2, and two Na+ binding sites (Na1 and Na2). The last 2 TMs (purple in Fig.6) are involved in the dimerization of the protein matrix. The first 10 TMs are organized into two sets of 5 TMs each; the second set is inverted with respect to the first, creating a pseudosymmetric structure. This inverse symmetry fold can be found in other secondary transporters (vSGLT, BetP, Mhp1, etc.)[210], even if they have low amino acid homology with LeuT. This has had important implications on modelling of secondary transport: crystal structures have been used as the template in homology modelling studies even when the overall homology between model and template is very low (~20%). This is the case with the homology models of hSERT and hDAT 29

based on LeuT[210, 211, 217]. Homology modelling has been used in Na+ and K+ channels as well, but it is in secondary transporters where it has really shone[217]. The success of LeuT as a model of secondary transport has been explored in several reviews[209-212, 217]. The secondary transporters provide a unique and challenging case for studies of ion coordination and its coupling to thermodynamics/kinetics of transport. Binding of the Na+ ions and the substrates induces the conformational changes necessary for transport. Moreover, the binding processes have an allosteric nature: binding of one species results in structural reorganizations facilitating the binding of a second species. For instance, the Na+ ion, which is not in direct contact with the Leu molecule (in site Na2) is instrumental for the stabilization of the main substrate binding site S1 and the Na1 sodium binding site[7, 210, 211]. Another template for our understanding of selective ion binding is GltPh, a homotrimeric amino acid transporter found in Pyrococcus Horikoshii which couples the translocation of 3 Na+ ions with one Asp residue (Figs.2 and 4)[214, 215]. The GltPh monomers consist of 8 TMs: 4 of them form a rigid trimerization domain (TMs 1, 2, 4 and 5, purple in Fig.6) while the transport domain is formed by the remaining 4 (TMs 3, 6, 7 and 8, orange in Fig.6) and two hairpin loops (HP1 (green) and HP2 (cyan)) which cradle the binding site. HP2 is fairly mobile and tends to open and uncover the binding site spontaneously even in the apo-state[218-221]. Allosteric effects govern the binding and transport in GltPh similarly to those in LeuT. In contrast to LeuT; however, in GltPh the alternate access mechanism is possibly operating via an elevator-like motion during which the transport domain slides along the trimerization scaffold (Fig.2)[214]. Binding and release of Na+ ions and aspartate by Glt Ph is allosterically coupled to the motions of HP1 and HP2 and the consecutive elevator-like translocation as determined from standard and biased MD[214, 218-222]. The binding site of one of the Na+ ions in Glt Ph is shown in Fig.3, to illustrate the type of ligands and the geometry of the coordination shell around the Na+. Unlike K+ in KcsA where the coordination involved uncharged carbonyl groups, the Na+ in GltPh are strongly coordinated to negatively charged residues, such as the ones involved in the EEEE, DEKA, etc. rings in the voltage gated Na+ channels. The crystal structures of LeuT and Glt Ph in inward and outward facing conformations provide important snapshots from their transport cycles and suggest possible binding sites for the ions and the substrates and possible mechanisms of the conformational changes consistent with transport. Yet, they lack the time resolution of transport and in the case of GltPh they do not show the position of the third Na+. MD simulations have therefore been used to probe for a possible binding site of Na3 in GltPh[222]. A putative second substrate binding site (S2) for Leu in LeuT was suggested as well on the basis of SMD calculations (pink sticks in Fig.6) but the determination of its exact location with experiment has been unsuccessful so far[223]. Interaction of monoamine transmitter molecules with protein models based on LeuT have been investigated with docking studies and MD[224-226]. Binding of these compounds in S1 is accomplished primarily via a salt bridge formed by the positively charged monoamine group and an Asp residue at the active site and via a number of secondary amine and hydroxyl group interactions with Tyr, Ser, Thr and Gly residues from the protein matrix. Unlike in channels, the meaning of “selectivity” in the world of secondary transporters is better defined in terms of ions once bound stabilizing a particular state and residing in the binding pocket. That is, transport or ion permeation events across secondary transporters are usually very slow processes and selectivity is essentially defined in terms of binding thermodynamics. It can usually be described as the preferential coordination of a certain ion or 30

substrate based on the energetics of the coordinating sites[68]. However, several recent studies suggested a more dynamic nature of ion binding and coordination to secondary transporters. For example, Zhao et al.[227] and Zomot et al.[228] reported on the dynamical transitions between charged site Na1 (directly coupled to co-transported solute) and tentative site Na1’ formed by one of the negatively-charged residues (E290). The ion occupancy of Na1’ destabilizes the bound substrate and is coupled to the gating dynamics of the protein[229]. The experimental binding affinities are usually averaged over the whole transport cycle, however the discrimination between ions and substrates is a result mostly of the binding thermodynamics since the permeation component, observed in channels is inapplicable. Hence, alchemical free energy methods like FEP have been used successfully to probe the selective coordination at the substrate binding sites both in LeuT and GltPh. Relative free energy calculations of Leu and Asp bound LeuT and GltPh with a different number of Na+ ions present in their Na+ binding sites demonstrated that the substrate binding depends on the ion occupancy in line with the allostericity of transport[209]. Selectivity of the crystallographically determined Na+ binding sites in LeuT and GltPh for monovalent cations other than Na+ (K+ and Li+) has been evaluated from FEP and successfully compared to experimental observations[222, 230]. Water is also of importance for secondary transport, especially in proton coupled transporters – in some cases water pathways remain uninterrupted even in the occluded and closed structures[231]. 4.4 Computational challenges Secondary transport is slow and although µs-long simulations of LeuT are already available[228], the full transport cycles of LeuT and Glt Ph are still beyond the capabilities of allatom MD with present computers. While some phenomena like the coordination of ions and substrates at the binding sites or the gating processes in the very beginning (or end) of the transport cycle can be observed in all-atom MD simulations (see above), the major structural reorganization which leads from an outward facing to an inward facing (or vice versa) structure has to be modelled with more affordable methods. One of the major problems in active transport studies is that the “reaction coordinates” describing the translocation of ions and the protein movements are not as well defined as in the ion channels. Uptake or release of ions and substrates from LeuT and GltPh has been studied with metadynamics and SMD calculations[219, 223], but they represent a very small portion of the full transport cycle. Thus, coarse grained or enhanced sampling approaches like the elastic network model, dynamics importance sampling, two-state anisotropic models, weighted ensemble path sampling or self-guided Langevin dynamics have been used in several different secondary transporters (Mhp1, LacY, LeuT, GltPh) to probe for possible movements in the protein which could lead to conformational changes relevant to transport[7, 212]. Rigid motions consistent with opposite transport of the dynamic and trimerization domain in GltPh were identified from such calculations. The allosteric coupling between the inward and outward gating in LeuT has been tackled with computational methods such as the N-body information theory [232], which identifies higher order correlated motions found in atomic fluctuations from equilibrated MD trajectories. Another option is the construction of an allosteric interaction network from multiple MD trajectories of different mutants and different ion and substrate occupations of the same transporter. Both methods identified residues in TMs 6 and 8 in LeuT as critical for the allosteric coupling between the substrate binding site (S1) and the intracellular

31

gate in agreement with mutagenesis, smFRET and EPR data. More information on these studies, as applied to LeuT and GltPh can be found in Refs.[7, 209, 212]. Most of the aforementioned techniques do not require a predefined “reaction coordinate” like in the PMF simulations and instead try to optimize the transition between an initial and a final state (e.g. the outward facing and inward facing states of a transporter) by following RMSD criteria or by analyzing large combination of higher order correlated motions. A common problem of this treatment is that the extraction of biologically meaningful results from the obtained conformational trajectory is very difficult. The string method with swarm of trajectories (see Appendix) allows for description of the conformational path in terms of a balanced set (neither too many, nor too few) potentially meaningful coordinates, such as bond lengths, valence and dihedral angles, and may be of use in such situations[15, 233]. Nevertheless, we would like to stress that the field of conformational studies in active transporters remains largely unexplored and the performance of the available methods is still insufficiently validated against experiment and benchmark biased all-atom MD methods like umbrella sampling (US, see Appendix).

32

5. Divalent cation channels and transporters: Substantial electronic effects are necessary 5.1 Summary This section summarizes the results from a few computational studies on the transport of Ca2+ and Mg2+ and explains why substantial electronic and polarization effects are required to describe the selective transport of divalent cations. 5.2 Transport of Ca2+ and Mg2+ The vast majority of Ca2+ ions are bound as phosphates in the bones and teeth[1]. Most of the free Ca2+ is found in the extracellular fluid where it is involved in blood clotting, nerve conduction and muscle contraction. Increased intracellular concentration of Ca2+ is linked to DNA replication and cell proliferation. Regulation of the intracellular Ca2+ concentration is critical for the cells and they have developed complex Ca2+-signalling mechanisms to ensure it. A detailed overview of the proteins involved in cell signalling and cell proliferation can be found in Ref.[234]. Transport of Ca2+ through the membrane occurs via Ca2+ channels (e.g. voltagegated Cav or channels, part of the Ca2+ signalling in the cell) and secondary transporters like NCX and the potassium dependent NCKX[235]. A well-characterized serco/endoplasmic Ca2+ATPase (SERCA) is responsible for pumping Ca2+ ions against their concentration gradients in the sarco/endoplasmic reticulum. Mg2+ activates ATP (ATP-Mg2+ complex) and is a cofactor of numerous enzymes responsible for aerobic and anaerobic energy generation[1]. It also participates in the oxidative phosphorylation processes in mitochondria and is vital for the synthesis and stabilization of DNA and RNA. Mg2+ can be involved in the regulation of Ca2+ in the Ca2+ pumps. The systems transporting Mg2+ through the membrane are poorly characterized. Two Mg2+ channels – CorA and MgtE have been subjects of some computational studies and will be discussed here. 5.3 Cav channels Votage-gated Cav channels are numerous and diverse, both in their response to membrane polarization (high-voltage vs. low-voltage channels) and their distribution in cells. A recent review dedicated to the classification of Cav channels and their structural and functional differences can be found in Ref.[236]. Voltage-gated Cav channels share similarities with both Kv and Nav channels: (1) Like Kv, they are very selective for Ca2+ and transport Ca2+ in a 1000:1 ratio with respect to Na+ ions, even though Ca2+ and Na+ ions have identical ionic diameters. (2) Like eukaryotic Nav, they feature a single alpha-subunit made of four non-identical domains. (3) Like prokaryotic Nav, their SF has a characteristic EEEE ring made of Glu residues. X-ray structures of Cav are not available yet. Instead, most of the experimental structural insight comes from sodium bacterial channels, such as NavAb[194, 195]. From a theoretical point of view Ca2+ binding and flux in model ion channels has been studied with reduced models using the integrated Poisson-Nernst, BD and SMD methods, for example. These studies were summarized in Ref.[153, 237]. Generally, it is believed that the ion permeation in the Cav pore follows the same rules as in Kv and Nav: several binding sites are present in the SF, one of them (the EEEE ring) is highly selective for Ca2+ and binds Ca2+ with strong affinity. This Ca2+ ion is 33

then dislodged from the binding site by electrostatic repulsion with other Ca2+ ions (or other divalent cations) as per the knock-on mechanism. Na+ ions (and other monovalent cations) do not generate strong enough electrostatic repulsion with Ca2+, cannot dislodge Ca2+ and are blocked by it. This explains the high selectivity of Cav in the presence of a much higher Na+ concentration. The lack of X-Ray structural data however leaves some unanswered questions, pertaining to the importance of the flexibility of the pore and the Ca2+ dehydration upon entry in the channel. 5.4 Primary transport with SERCA Cytosol Ca2+ is pumped in the sarco/endoplasmic reticulum by the sarco/endoplasmic Ca2+-ATPase (SERCA)[238, 239]. SERCA has a structure typical for the P-type ATPases, such as the Na,K-pump. It consists of 4 domains – N (nucleotide domain featuring the binding sites of ATP), P (phosphorylation domain, where an Asp residue is phosphorylated by ATP during the catalytic cycle), A (actuator domain, involved in opening and closing of the pump and in dephosphorylation of the P domain) and TM (transmembrane domain with binding sites for two Ca2+), see Fig.7. SERCA pumps 2 Ca2+ ions using the hydrolysis of 1 ATP molecule and the anti-transport of 2-4 protons. Several crystal structures of rabbit SERCA, in two different conformations and with different ion (Ca2+, Mg2+) loads exist. The two major conformations of the catalytic cycle are E1 (with Ca2+ binding sites exposed to the cytosol, involved in uptake of Ca2+) and E2 (with Ca2+ binding sites exposed to the sarco/endoplasmic reticulum lumen, involved in release of Ca2+). The transition between E2 and E1 obeys the alternate access mechanism of active transport and occurs via Ca2+/H+ exchange. Mg2+ has been suggested as a catalytically important ion for SERCA as well[239]. Studies of the catalytic cycle of SERCA are fairly challenging due to the large size of the system and the complex allosteric effects during transport. Nevertheless, it has been studied in several impressively long (up to µs)[239-242] all-atom MD simulations as well as a number of shorter MD and coarse-grained simulations[243-248]. Insight was gained into the ion pathways to the Ca2+ binding sites[238, 243-246], the binding process[239, 244, 246], the conformational changes in the cytoplasmic part (P, N, A domains) necessary for ATP hydrolysis and binding induced by Ca2+,[242] the role of Mg2+ and K+ ions[239] and mediator proteins like phospholamban and sarcolipin[240, 241] into stabilization of catalytically relevant states, and the protonation effects[238, 247] of the acidic residues coordinating Ca2+ in good agreement with the limited experimental data from X-Ray structures and FRET spectra. A recent MD study suggested different putative pathways for Ca2+ and H+[238]. Motion Tree diagrams for SERCA revealed that the coupled motion between the TM and cytoplasmic domain is a result of small flexible kinks of the molecule which propel a few large rigid domains into the necessary conformation[249]. Binding affinities for Ca2+ and Mg2+ at the water/protein interface and within TM evaluated from DFT calculations implied that Mg2+ is filtered out of the TM pore by binding to Mg2+-selective sites at the water/protein interface[248]. These studies have already expanded our knowledge about the catalytic cycle of SERCA (and by association, other P-type ATPases) beyond the limits of the experimental data. However, the transport phenomena in SERCA are still far from understood. Longer MD simulations in the future will doubtlessly provide further clarification of the complex process of Ca2+ transport in SERCA and other ion pumps.

34

5.5 Secondary transport with NCX Ca2+ extrusion from the cells is accomplished by secondary transporters, like NCX, which couple the Ca2+ translocation with translocation of 3 Na+ ions along their concentration gradient[250]. The transport follows the alternate access mechanism (e.g. 3 Na+ ions are first transported into the cell from the extracellular space followed by extrusion of 1 Ca2+). Recently, several crystal structures of NCX_Mj from Methanococcus Jannaschi in outward open and outward occluded forms were crystalized[251, 252] (Fig.7) and used for MD and metadynamics studies[252, 253]. Opening and occlusion of NCX were related to the ion load in the active site: fully loaded NCX_Mj with 3 Na+ and 1 water molecule led to occlusion of the structure. This regulation of the catalytic cycle of NCX by the ion occupancy at the active site was observed independently both in the metadynamics simulations (starting from the outward open occluded state, Fig.7) and the X-Ray structures. The coordination of Ca2+ at the NCX_Mj active site is presented in Fig.3. As expected, this coordination is strongly dependent on electrostatic attraction between the Ca2+ ion and a few negatively charged amino acids occupying its first coordination shell. 5.6 Mg2+ transport by MgtE and CorA Crystal structures for two prokaryotic ligand-gated Mg2+ channels, MgtE and CorA, have been resolved and have prompted computational investigations into binding and transport of Mg2+ through the cellular membrane. Their structure and functionality have been reviewed in Ref.[254]. MgtE is a homodimer, featuring a TM part of 2 x 5 TMs and a large regulatory intracellular region consisting of three domains – N (N-terminal domain), and two cystathionine beta-synthase domains, CBS1 and CBS2 (Fig.8)[254, 255]. A narrow channel through the membrane is formed by kinked and twisted Glu and Pro rich TMs 1, 2 and 5. A “plug helix” (blue in Fig.8) connects the TM channel with the cytoplasmic region of MgtE and serves as a hydrophobic gate, controlled by the regulatory intracellular region. The N and CBS domains of the regulatory region are linked together through 12 Mg2+ ions bound at highly conserved acidic residues (Fig.3). The binding of Mg2+ to these sites is related to channel closure. MD simulations of Mg2+ bound and Mg2+ free MgtE models confirm that the channel is closed in the presence of Mg2+ and that this closed structure cannot be maintained by monovalent cations in line with experiment[256]. In particular, binding of Mg2+ at the border between the plug helix and one of the CBS domains has been identified as critical for the channel closure. CorA is a homopentameric channel, shaped like a funnel (Fig.8). Similarly to MgtE, CorA can also be divided into TM and regulatory intracellular domain[254]. Regulatory binding sites for Mg2+ were identified in the intracellular domain at the borders between the monomers. The TM region features several short alpha helices (TM2) and the long helices (TM1) which continue into the intracellular domain and shape the channel pore. The selectivity filter consists of five Gly-MetAsn (GMN) triads at the pore entrance (Fig.8) and controls the entry of hydrated Mg2+ in the lumen. Gating of the channel is induced by Mg2+ binding in the regulatory domain which spreads to the TM region via the long TM1 helices. The lower part of the TM2 helices and the midsection of TM1 are lined with positively charged residues (basic ring) which impose an additional barrier for the Mg2+ entry into the cytoplasm, as determined from MD studies. 35

Beneath the basic ring, in the beginning of the intracellular region is situated an acidic ring of negatively charged amino acids. It is assumed that these two rings interact when the channel is opened and remove the electrostatic barrier for Mg2+ entry. MD results in the absence of Mg2+ indicate asymmetric motions of the intracellular region which lead to dilation and wetting of the pore, consistent with the open Mg2+ permeable structure[257]. In contrast, binding of Mg2+ to the regulatory domain induces a symmetric non-permeable geometry[257, 258]. In addition, hydrated Mg2+ passes through the channel and indirect interactions between the ion and the GMN motifs in CorA occur through the water from the first coordination sphere of Mg2+[257]. 5.7 Computational challenges Considering the abundance and importance of Ca2+ and Mg2+ in the human body, it is somewhat surprising that the number of theoretical studies dedicated to the transport of these two elements is so limited (Tables S1-S4, Supporting Information). One of the major hurdles in the computational investigations is the lack of crystal structures. Moreover, the current force fields commonly used in MD simulations lack polarization and electronic effects which are critical for the proper description of highly charged ions like Ca2+ and Mg2+, involved in polarization and charge-transfer effects beyond the capabilities of the standard classical force fields (see also computational challenges in sodium channels). In some of the more recent MD simulations special corrections in the van der Waals parameters for Ca2+ have been implemented to improve the description of Ca2+ binding to the protein matrix[252, 253]. These corrections however are not validated on diverse sets of model systems. Most MD calculations involving Mg2+ are performed without any corrections. Such corrections are likely necessary for small ions of high charge density like Mg2+. A possible solution to these problems would be the use of polarizable force fields, which are currently under development and validation. For example, calculations of Ca2+ binding with the Drude force field show considerable improvement over results obtained with CHARMM36[73]. As mentioned above, however, the polarizable force fields still lack necessary parameters for the description of ion-amino acid interactions.

36

6. Anion channels and transporters: unexplored territory 6.1 Summary This section highlights the fact that much is unknown about structural mechanisms regulating anion channels and transporters. We review the limited available research in this field and offer some rationale for why we need to investigate these types of transmembrane proteins. 6.2 Chloride channels and transporters Chloride is the most abundant halide in the human body[1]. It is found in large concentrations in the extracellular fluids where it serves as a Na+ counterion and adopts some of the functions of Na+, like nerve impulse conduction or regulation of fluid volume and osmotic pressure and acidity. It is an essential part of gastric acid and participates in the transportation of CO2 via the Cl-/HCO3- exchanger protein, Band 3, found in red blood cells and kidneys[259]. Its homeostasis is achieved with channels and transporters of the CLC family[260]. Other channels which are slightly selective for anions, such as the mechanosensitive channels MscS[261] and the mitochondrial VDAC[262] are also known to transport chloride. VDAC however is mostly associated with the transport of larger organic anions and we will not include it here. Computational studies of VDAC were recently reviewed and the interested readers can find this information in Ref.[262]. The CLC family is divided into Cl- channels, which ensure rapid transport along the concentration gradient of Cl- and Cl-/H+ antiporters which use the translocation of 2 Cl- toward the intracellular space for the extrusion of 1 H+[260]. Chloride coupled antiport of H+ occurs also in the Cl- channels of this family, however to a much lower extent (millions of Cl- per H+). The CLC channels and transporters have very high homology, which implies that they share a similar structure even though their transport mechanisms are so different. In fact, mutations of key gate residues (red sticks in Fig.9) in the CLC-ec1 transporter lead to channel-like conductance of the mutant. This curious behavior of the CLC systems is reflected in a review dedicated to the similarities and differences between channels and transporters of the CLC family, where they have been labeled “proteins with borderline personalities”[260]. Most of the computational studies are based on crystal structures of the CLC-ec1 transporter from Escherichia coli (Fig.9)[263, 264]. CLC-ec1 is a homodimer with two independently operating subunits. Each subunit consists of 18 TMs organized in an inverse repeat of 2 groups of 9 TMs, similarly to other secondary transporters (LeuT, NCX, see above). Three binding sites of low affinity can be identified in each subunit – Sext (external), Scen (central) and Sint (internal) – and form an arched pathway for the Cl- (Fig.9). The selectivity of CLC-ec1 is determined by the Scen site with the following selectivity sequence: SCN- > Cl- > Br> NO3- > I-[260]. The Cl- bound at Scen (Fig.3) is dehydrated and coordinated by OH and NH hydrogens of the surrounding Ser, Tyr, Gly, Ile, Phe residues and by the acidic side chain of a Glu residue, known as the Gluex. Gluex and the Ser and Tyr residue at Scen are likely the extracellular and intracellular gates, involved in the Cl- transport (Fig.9). Double mutations at these residues increase 20-fold the rate of Cl- transport and abolish the conformational rearrangements consistent with transport in the WT transporter[260]. A second Glu residue (Gluin) is situated close to the cytosol and has been regarded as part of the H+ pathway (along 37

with Gluex) through CLC-ec1. The eukaryotic CLC channels (and many of the prokaryotic ones) have in addition cytoplasmic domains (CBS units like in MgtE and CorA, see above), which bind nucleotides like ATP and modulate the transport in both CLC channels and transporters[260]. The putative ion pathways for Cl- and H+ ions in CLC models have been studied with a variety of methods – PB[265], BD[147, 157], MC[266], MD[267-272], biased MD (ABF[273], SMD[269], metadynamics[268], US[268, 272]), QM[274], and QM/MM[275]. The simulated PMFs revealed additional putative Cl- binding sites[265, 268, 272] and indicated that binding of at least two Cl- was necessary for energetically favorable rearrangements along the H+ pathway (between Gluex and Gluin) in line with the observed 2:1 stoichiometry of transport[271]. The binding of Cl- and protonation of Gluin was linked to conformational changes allowing the formation of transient water strings used for proton transfer[267, 273, 274]. The movement of Clions through CLC was found to obey similar principles to the knock-on mechanism in K+ channels (i.e. the electrostatic repulsion of the Cl- decreases the energy barriers of transport)[268, 272]. Exchange of water molecules with the bulk was observed at the Cl- binding sites closest to the periplasm and cytoplasm[275]. The role of the Gluex residue as a competitor of Cl- for Sext and Scen was also examined: it was confirmed that neutralization of Gluex (via H+ or Na+) stabilizes an open structure, in which the side chain of the Glu has moved away from the Clpathway[270, 273]. In addition, ion conduction in model eukaryotic CLC channels, based on the CLC-ec1 structure has been studied with BD, affording good comparison with experimental I-V curves and the presence of a larger number of basic residues along the Cl- pathway has been discussed in conjunction with the higher speed of Cl- transfer in CLC channels[153, 157]. DFT calculations of thermodynamic preferences for Cl-, Br- and F-—binding to the Scen site demonstrate that Scen is weakly selective and its selectivity depends on the presence and orientation of the protonated Gluex[276]. Moreover, charge transfer from the ion to the protein occurs in the cases of Br- and Cl- binding and the observed F- blockage of CLC might be due to interference of F- (which has high H+ affinity) with the transported H+. The MscS are slightly anion selective wide channels which regulate the intracellular pressure in bacterial cells[261]. Another known mechanosensitive channel is MscL which features a larger very conductive and non-selective pore. The smaller pore of MscS prevents the free movement of solutes and as a result MscS has been studied more than MscL[38]. As their name suggests the mechanosensitive channels open under the influence of the surface tension of the membrane and close when the tension decreases. They have been extensively studied with patch-clamp experiments and, understandably, a lot of the existing simulations aim at reproducing ion conductances[277-281]. It was determined computationally that the first known crystal structure of MscS[282] represents a semi-open state[278, 279] which would transform into a more open conformation if external electric field is applied during the simulations[281]. Later, open and closed structures of MscS were obtained from MD simulations with restraints incorporated from electron paramagnetic resonance data and better agreement with experimental conductances was achieved for the opened structure[283, 284]. The Cl- currents have been simulated also with continuum electrostatics methods. These results compare well with more sophisticated BD simulations[280, 285]. The gating process in mechanosensitive channels has generated a considerable interest and has been studied both in all-atom and continuum and coarse-grained MD simulations. 38

Significantly, MscL calculations with the MARTINI[286] model managed to reproduce the gating of the channel by membrane stretching[287, 288]. The theoretical investigations on MscS and MscL gating have been reviewed in more detail in Ref.[38]. 6.3 Inorganic Phosphate Free PO43- in the body is known as inorganic phosphate and at the physiological pH is found in the forms of H2PO4- and HPO42-[1]. Inorganic phosphates (Pi) serve as the main buffering agents in urine. Phosphate is predominantly found in its bound forms – with Ca2+ as a principal component of the skeleton and teeth, in adenosine phosphates (ATP, ADP, AMP) and in the two nucleic acids, DNA and RNA. Free PO43- is released after hydrolysis of the adenosine phosphates during the dephosphorylation processes. Its transport through the cellular membrane is ensured by phosphate selective porins like OprP and Na+ coupled transporters like GlpT. OprP and GlpT are the most common subjects in computational investigations of phosphate transport, since their crystal structures are available[289, 290]. GlpT is a member of the large and diverse Major Facilitator Superfamily, and serves as an antiporter, which exchanges glycerol-3-phosphate (G3P) for inorganic phosphate ions. It consists of 12 TMs, organized in two pseudo-symmetric folds (2 x 6 TMs)[289]. The alternate access in GlpT follows a rocker-switch mechanism, characterized by reorganization of the alpha helices leading to an alternating transition between inward open structure (binding phosphate) and outward open structure (binding G3P)[291]. The available crystal structure for GlpT represents the inward facing state (Fig.10)[289]. A set of 2 Arg and 1 His residues has been discussed as the putative binding site for Pi. Binding of Pi likely brings the two Arg residues together and induces the transition to the outward open structure[291]. A review of modelling studies of GlpT binding and transport can be found in Ref. [291]. Unbiased MD simulations were used for monitoring of Pi movement toward the suggested binding site and initial reorganization of the protein geometry consistent with the rocker-switch model[292]. Spontaneous binding was observed for H2PO4-, HPO42-, and G3P- and G3P2-. Moreover, binding of Pi to the Arg/His binding site led to salt bridge reorganization at the periplasmic site of the protein[293]. Protonation of the His residue was necessary for Pi binding. Large conformational changes, leading to a possibly occluded state were simulated with reduced representation of the protein in the framework of the anisotropic network model[291]. OprP is an inorganic phosphate selective porin, with structure somewhat resembling the OmpF system (see above). It is a trimer consisting of large beta-barrelled pores (Fig.10)[290]. Three loops positioned in the lumen of each pore form a narrow passage for ions and determine the selectivity for inorganic phosphate. The affinity series for OprP is PO43- > SO42- > Cl-[294]. The ion pathway shaped by the pore and the loops is lined with a number of positively charged (Arg and Lys) amino-acids. Interestingly, a single Asp residue is also found on the ion pathway[295]. A snapshot of Pi coordinated to two Arg residues in the OprP protein as determined from X-Ray diffraction is shown in Fig.3. The role of the Arg, Lys and Asp residues in the Pi transport in OprP has been elucidated with biased MD calculations (SMD and ABF)[295, 296]. The simulated PMFs reveal that the Lys residues at the periplasmic site of the pore attract the anions and propel them along their pathway. A central Arg residue (Arg 133, red in Fig.10) was identified as an important site 39

involved in the dehydration of the ion[296]. The negative Asp residue (Asp 94, blue in Fig.10) controls the binding of Pi and ensures its fast translocation[295]. PMFs of other ions (SO42-, Cl-, K+) in the pore are in agreement with the observed experimental selectivity of OprP[294]. 6.4. Other ions – F-, Br-, I-, HCO3-, SO42Other biologically relevant small inorganic ions found in low and high concentration in the human body are also known. Their transport systems are poorly characterized and are investigated in a rather limited number of modelling studies (see Tables S1-S4). We have decided to mention a few of these ions here to draw attention to this relatively unexplored area of ion transport. The halides F-, Br-, I-, are considered essential trace elements and are found in much lower concentrations in the human body than Cl-[1]. F- is involved in ossification processes relevant to the formation of teeth and bones. Br- is a cofactor in enzymes involved in the production of collagen IV and is used by the eosinophil cells in their antiparasitic function. Iodide is necessary for the synthesis of the thyroid hormones and its transfer through the membrane is ensured by the sodium-iodide symporter NIS[297]. NIS is expressed in various other organs – mammary and salivary glands, stomach and intestines, testicles, etc. The role of I- in these tissues is not clear, however iodine deficiency has been linked to extrathyroidal cancers[298]. Transport of the halide ions through the membrane can occur through the CLC channels (except for F-) or other channels which are modestly selective for monovalent anions but transport predominantly Cldue to its high concentration in the body fluids. X-ray structures of prokaryotic F- channels were reported recently and open the grounds for further investigation of their unorthodox binding and selectivity, which is likely a result of interactions between F- and conserved Phe rings[299]. HCO3- is a product of the CO2 metabolism regulated by the carbonic anhydrase enzymes in

the erythrocytes[1]. It serves as a buffer component in the pH regulation of blood. Its abundance and importance are in striking contrast with the lack of information on its regulation. Its transport through the membrane is achieved via Cl-, Na+ and K+ coupled transporters, such as the Band 3 anionic exchanger[259]. Band 3 is considered one of the fastest transporters in the human body and is also one of the most common proteins in the erythrocytes. The crystal structure of its TM anion exchange domain was resolved very recently[259]. This achievement will hopefully stimulate more computational investigations of HCO3- transport. SO42- is predominantly found in the plasma where it is involved in the regulation of ionic homeostasis[1]. It is the most abundant divalent anion in cells and is necessary for the synthesis of a variety of essential compounds such as chondroitin sulfate, dermatan sulfate, etc. Its transport is regulated by sulphate permeases which can be either channels or Na+ coupled transporters[300]. Crystal structures of sulfate binding proteins (SBP) which have been linked to sulfate transport through the membrane have been available for several decades[301]. The binding, selectivity and transport in SBP have been studied with QM[302] and MD[265] calculations. Nevertheless, our understanding of the regulation of these ions is still very limited.

40

7. Outlook and Challenges 7.1 Summary We emphasize some underdeveloped fields of ion transport and the need for more advanced computational approaches that must take into account the subtle polarization and electronic effects and long-time sampling. Computational methods have established themselves as invaluable tools for the elucidation of ion transport through channels and transporters. They provide insight into the structural and dynamic underpinnings of the ion translocation process, which is beyond the capabilities of the experimental methods. Their diversity and affordability have contributed to their increased popularity in the last two decades. In terms of quality, modern computational methods routinely afford qualitative agreement with experimentally observed selectivity and structural trends. Quantitative agreement is also found with those experimental methods, which provide measures for comparison. A plethora of successful simulations have shaped our current understanding of K+ and Na+ transport by selective channels and transporters and attest to both the applicability and power of the theoretical methods. The computational toolkit available for studies of ion transport is already quite impressive and features a large variety of methods based on the QM and MM formalisms (and combinations thereof). Nevertheless, it keeps expanding and improving as new methods are developed and old methods are applied to more demanding calculations. A shift toward bigger, more complicated and realistic models is constantly observed. MD calculations have been used to direct and complement expensive experimental investigations, such as mutagenesis studies, IR, NMR and Cryo-EM. As discussed above, MD simulations can pinpoint crucial interactions that may affect the ion selectivity and permeation; afterwards mutations could be done to validate the roles of these interactions. However, it must be noted that not all mutations identified from MD simulations are achievable in physiological conditions. All-atom sampling allows combination of simulations with 2D IR spectra which can reveal fast dynamics and structural changes in localized environments, transforming the entire field of computational studies. These fused techniques have been used recently to test different ion permeation models such as knock-on and hard-knock[67]. MD simulations have also been implemented to support the tedious task of mapping NMR data into reliable structural features of protein channels. With the recent advancement of Cryo-EM techniques for identifying large complex proteins in physiological environments with near atomic-resolutions (3.8-4.2 Å), allatom MD simulations have become an extremely useful tool for the enhancement of the resolutions obtained with Cryo-EM. Regardless of its capabilities and successes, however, every modelling study needs a starting point: a 3D structure, a protein sequence or information about key residues from experimental measurements and mutagenesis studies. Thus, theory and experiment are always intertwined and deficiency in one leads to a potentially deficient interpretation in the other. It is not surprising then that the computational studies of systems lacking structural insight from experiment are few and far between. The striking contrast between the number of studies on the well-characterized K+ channels and the unknown SO42transporters (Tables S1-S4) conveys how much the computational methods rely on experiment. As such, lack of structural data (e.g. crystal structures) is by far the greatest hindrance for the computational investigations, along with computational hardware limitations. The transport of 41

divalent cations and mono- and polyvalent small inorganic anions remains essentially unexplored from the computational point of view. Ions like Ca2+, Mg2+, HCO3–, SO42– and Pi which are both abundant and involved in some of the most vital metabolic processes, would definitely benefit from further computational investigations. Even the most explored transport systems – those of K+ and Na+ are still far from completely understood, especially in the domain of active transport. Other problems arise from the numerous assumptions and simplifications adopted in the computational methods. As the computers grow more powerful, the MD simulations grow longer and artifacts of improper parameterization, unsuitable for the new timescales start to emerge. The rigid water models are not suitable for the proper description of hydrogen bonding and for interactions at water/other-phase interfaces. The small number of theoretical studies dedicated to anion transport (excluding chloride) has another implication: the force field parameters for these anions are of uncertain quality and may need further validation and re-parameterization. Lack of polarization effects, crucial for the proper description of structures and energetics involving polyvalent ions, has interfered with the studies of Ca2+ and Mg2+ and anion transport. Complex electronic effects are also necessary for the description of coordination and functionality of several other physiologically relevant metal ions such as transition metals or Zn2+, which are almost untreatable with MD at this point. Some strides have been made toward the computational treatment of Zn2+ (e.g recent MD studies of the inhibition of transport in protein channels[303] by Zn2+ or the QM-evaluated thermodynamics of binding of Zn2+ to molecules mimicking potential protein binding sites[304]); however, much more work is necessary for the production and validation of good force fields that can handle metal ions and heavy elements featuring diffuse d-orbitals. Polarizable force fields have been suggested as a possible solution of these issues and have been under development for decades. However, they are computationally demanding, underdeveloped and their performance is still evaluated. QM and QM/MM methods which explicitly include electronic degrees of freedom are still beyond reach for large protein models. In this respect, DFTB is emerging as a promising method for the QM treatment of macromolecules due to its good balance between computational expense and accuracy. It was used recently in a validation study for ion-protein interactions as a benchmark for comparison with Drude and CHARMM36. Generally, QM methods are not very suitable for studies of ion permeation, which is a dynamic and stochastic process. Ab initio molecular dynamics may capture the time evolution and statistical aspects of transport; however, it is extremely demanding and can currently be applied only to very small models. QM/MM may help in the sampling of architectural forces acting on the selectivity filters or binding/transport domains, which are treated with QM. Technically, the bottleneck of QM/MM efficiency lies in the QM domains. Thus, sampling a few hundred nanoseconds using QM/MM is still very demanding. Nevertheless, it is essential to provide more reliable answers for divalent-ion transport and selectivity. The ion-coupled conformational transitions in the transporters are complicated processes in which the definition of reaction coordinates is ambiguous and subjective. Approaches like the string method with swarm of trajectories, targeted molecular dynamics or self-guided Langevin dynamics do not need to bias the simulation along a predefined reaction coordinate and may eventually lead to physically meaningful results, but the analysis of these results in terms of biologically sensible coordinates (interatomic distances, rotations, etc.) can be quite complicated. Most of the reduced representation methods for conformational analysis can give limited information about the general geometry reorganization, but the neglect of the side 42

chains, which may be a critical part of the allosteric networks in secondary transporters, is still a major deficiency. The identification of the allosteric interaction network elements is a difficult task, which requires cross reference of numerous individual simulations of different states and ion loads. New algorithms which automate the search of such elements may need to be developed to standardize and facilitate such work. The conformational changes in transporters (and channels) remain unexplored even with the available methodologies and represent a field of study with unlimited potential for discovery. Presumably, most of these issues can be resolved with all-atom MD simulations which cover the whole relevant time scale of transport. Constant increase in the simulation times is evident in the scientific literature in which, 20 years ago, simulations of several ns were the longest one could possibly do. Today, the longest MD simulation of ion transport spans milliseconds and simulations approaching the microsecond scales are becoming more and more widespread, even for systems of several hundreds of thousands of atoms. It is clear, that this trend will continue with the development of more powerful computers. This, however, may necessitate reparameterization of the present force fields which are optimized for much shorter time scales. Despite the numerous methodological challenges described above, there are no doubts about the significant progress in our modelling and understanding of membrane proteins permeability, ion selectivity, and gating on a structural level. However, the meaningful comparison to experimental observables is often difficult. Even for idealized processes (e.g. model ionic concentration, absence of other reactions in the system, idealized representation of bilayer structure, etc.), one has to develop a comprehensive multi-scale protocol to connect computed properties (gating charge transfer, permeability) from simple models treated at the QM level, single protein simulations common in the MD community, to actual observables, often recorded at the whole cell level. Several attempts were made in the computational electrophysiology community, where I-V relations were computed and compared to the results of single-channel experiments or spectroscopic measurements. The connection between scales often can be achieved by developing cellular models for biological phenomena such as action potential propagation or rates of solute transport based on kinetic and thermodynamic information directly extracted from computational studies. A particularly popular approach is based on the application of Markov-states modelling (MSM)[305-309] to analysis of the conformational space sampled in simulations. There are apparent difficulties in defining relevant metrics of such a space in complex protein environments. The MSM techniques has been used successfully to model various biological phenomena including ion permeation in pores, dynamics of small molecules in lipid bilayers and dynamics of proteins[141, 308]. While MSM requires an exhaustive sampling of the conformational space, new simulation technologies are providing the required tools. Hopefully, additional work in that domain would allow one to sample ion binding and coordination dynamics directly and then organize outputs of computer simulations into an elegant theoretical framework for understanding the thermodynamics and kinetics of selective ion transport in its entirety.

43

APPENDIX - Computational methods used in studies of ion transport The computational methods used in the referenced literature (Tables S1-S4) can be organized into three sections – methods based on quantum mechanics (QM), methods based on classical and statistical mechanics (molecular mechanics, MM) and hybrid QM/MM methods. In this section we provide a brief description of the relevant methods in these three classes with emphasis on their strengths, weaknesses and applicability. Introductory books providing details on the concepts and mathematical expressions behind most of these methods are listed as Ref.[8, 310, 311]. Additional sources are provided for the methods missing from Ref.[8, 310, 311]. A1. Quantum Mechanical Methods

In quantum chemistry a chemical system can be described by a wave function Ψ, , which is a function of the coordinates (spatial and spin) of the particles forming the system (nuclei, R, and electrons, r). The energy of this system can be found when a Hamiltonian , , describing the relevant kinetic and potential energy interactions between the operator,  particles (and existing external electric or magnetic fields if necessary) is applied to Ψ, [8, 310, 311]. , Ψ,  = Ψ,  

(A1)

The form of the real wave function and Hamiltonian are very complicated, so a number of approximations are usually employed for the solution of eq.A1. Nevertheless, once the energy of the system is determined, a variety of experimentally observable properties (IR frequencies, NMR shielding constants, etc.) can be calculated as derivatives of the energy with respect to nuclear coordinates, external electric and magnetic field, nuclear and electron spins, and others. The optimized Ψ yields important information about the electron distribution within the molecular orbitals of the system, which can be related to atomic charges, dipole moments, and magnetic and optical properties of the molecules. QM methods which are based on the wave function Ψ are known as wave function methods. Among them are Hartree-Fock and a variety of post Hartree-Fock (HF) methods, such as Møller-Plesset perturbation theory, configuration interaction, and coupled clusters, which include dynamic electronic correlation. These methods use one reference determinant of the wave function; however, there also are multi configurational and multi reference QM methods using a set of different determinants of Ψ which can yield the limit of the static and dynamic electronic correlation. The wave function methods provide high accuracy in the calculated properties, however they are extremely computationally demanding, due to their multi-dimensional (at least ON, where N is the number of electrons) Ψ. Thus, they are limited to small models, with a small number of nuclei (10-50 atoms) and electrons. Alternatively, one can use the ground state electron density of the system instead of Ψ, since it determines uniquely the ground state properties of the system. This is the basis of the density Functional Theory (DFT) methods. Unlike the wave-function methods where Ψ represents a collection of interacting electrons and the results can be systematically improved by expansion of Ψ with additional terms, representing alternative electron configurations, in Kohn44

Sham DFT the wave function describes a set of non-interacting electrons and the necessary interaction terms are expressed as a functional of the density of this non-interacting system. Currently, a large number of model functionals exist and they can be grouped into several classes (i.e. LDA, GGA, metaGGA, Hybrid, long-range corrected/separated, dispersion corrected) based on the way the density is evaluated at certain points in space. The ground state electron density is only 3 dimensional, therefore DFT allows for QM treatment of much larger systems (~100 atoms). It is the current QM method of choice for protein modelling. A review dedicated to the recent advances in DFT with a focus on biomolecular modelling can be found in Ref.[312]. Both wave-function and DFT methods can be partially simplified by substitution of some terms in their energy ansatz with tabulated or parametrized values. This is done in semiempirical methods based on the HF formalism (AM1, PM3, etc.) and in tight-binding formulations[21] of DFT (DFTB1, DFTB2, DFTB3). These methods are less accurate than postHF and DFT but allow higher computational speed and can be applied to bigger systems. The information, which could be extracted from the QM methods is diverse: optimal geometries, charges, dipole moments, polarizabilities, interaction energies, various spectroscopic parameters observed with experimental techniques (IR, UV/VIS, EPR, NMR, CD, MCD), some thermodynamic properties (enthalpies and entropies calculated from frequency calculation), etc. QM methods allow for studies of charge transfer, bond formation and bond dissociation and can be used in investigations of chemical reaction mechanisms. They have been instrumental for the parametrization of MM force fields. In protein-ion interactions they are usually applied to reduced models of the protein binding sites or to small biomimetic compounds which mimic the actual binding sites. More information on the application of QM methods in biochemistry can be found in a recent review by Cui [177]. A2. MM methods In the MM methods the electrons are not taken into account explicitly. Instead the system is represented as a collection of spherical particles of specific size and charge, connected by harmonically oscillating strings[8, 310, 311]. The energy of the system is expressed in terms of bonded (bond, valence angles, dihedrals) and non-bonded (van der Waals, electrostatic) interactions: =   −   +  



+ 4)* +, 

*2



/

-* . *

 −   +

−,

0

'

-* 3 3* . 1 +

* )* 

 1 + cos $ % − & 

*2

(A2) The bond and angle terms are most often modeled after a harmonic oscillator, the electrostatic term – with the Coulomb law and the van der Waals term – by a Lennard-Jones potential. The terms in eq.A2 rely on pre-existing parameters such as partial charges (3 ),  equilibrium bond ( ) and angle ( ) values, force constants ( ,  ,  ) and )* and -* parameters for the Lennard-Jones potential. The functional form of along with the relevant parameters is known as a force field. Examples of force fields widely used in protein-ion interactions are CHARMM, AMBER and GROMOS. A list of popular force fields and comparison of their performance can be found in Ref.[8, 310]. The quality of the force fields 45

depends heavily on their parameters which are developed from QM calculations and fitting to experimental data. The ad-hoc nature of the parameters makes them non-transferable and applicable only to certain atom types and while the common amino acids, nucleic acids, lipids, sugars and ions have already been parametrized in most force fields, lack of parameters (or insufficient validation of existing parameters) for some less common compounds is still a serious problem in MM calculations. Nevertheless, the MM representation of the energy is computationally inexpensive and allows for calculations of thousands or even a million of atoms. Moreover, the expression in eq.A2 can be easily coupled to randomized algorithms (such as Monte Carlo) or equations of motion, which allows for sampling of ensembles of structures, leading to proper evaluation of thermodynamic properties like entropy and free energy, which require a statistical approach. When a force field is used in conjunction with equations of motion, i.e. for evaluation of the dynamical changes in a system over a period of time, this is known as Molecular Dynamics (MD). MD has been used widely in the study of ion binding and transport phenomena in proteins; most referenced papers in Tables S1-S4 are based on MD studies. The sampling of the possible configurations (phase space) that a large model consisting of a protein, lipid membrane, explicit water and ions can assume, is computationally demanding and time consuming. Transport in secondary transporters often occurs at the µs or ms scales. Allatom MD simulations spanning such long time scales are currently possible for small models of several thousand atoms. In big models, the sampling in a desired direction (one that is consistent with geometry changes of interest, i.e. along a reaction coordinate) can be accelerated with the use of biasing potentials which are imposed on the system in different ways. The reaction coordinate is often divided into windows and numerous calculations are submitted at the same time, each sampling a different portion of the reaction coordinate. Examples of biased MD techniques include Umbrella sampling (US), adaptive biasing force (ABF), metadynamics, and steered molecular dynamics (SMD)[8, 114]. These methods allow the evaluation of one- or poly-dimensional potentials of mean force (PMF) - the free energy with respect to the chosen (one or more) reaction coordinates. Post-processing methods like the Weighted Histogram Analysis Method (WHAM)[313] are employed to remove the biasing potential and to assemble the information from all windows into one continuous PMF. PMFs can be extracted also from non-biased calculations via probability analysis assuming that these calculations have sampled properly the phase space. Biasing potentials are also employed for studies of large conformational changes underlying active ion transport. In such cases, the transition to a certain conformational state is accelerated by pushing the system in this direction, using root mean square deviation (RMSD) criteria (as in Targeted Molecular Dynamics (TMD)[314]) or allowing the system to evolve based on feedback from atomic forces (self-guided Langevin dynamics[315], String method[233]). The calculated PMFs can be related to some important properties measured experimentally. For instance, the maximum ionic conductance gmax in an ion channel can be calculated from a given PMF (W(z)) as follows 7 ?/ ABC/89 : 〉?/ 〈 ?BC/89 : 〉?/ 〈 456 = @ @ (A3) 7 => 89 :;

Here, the brackets represent the average evaluated over the length L of the channel, and the diffusion coefficient => may be computed from the Laplace transformation of the velocity autocorrelation function.

46

The equilibrium dissociation constant per mole for an ion in the cylindrical constraints of the channel can also be evaluated from the PMFs: C

FG?/ = H  IC KLM J>@ ?BC/89 : KNO

(A4)

The selectivity of a protein binding site for a specific ion, which is a measure of how strongly the ions bind to this site, is often estimated from the sign and magnitude of the free energy change ΔΔGA → B calculated upon substitution of a bound ion (A) with another one (B) in the protein and in bulk: ΔΔGA → B = ΔGU'V A → B − ΔGW8 A → B =

XYU'V Z − YU'V [\ − ]YW8 Z − YW8 [^ (A5) Usually, the ion substitutions are achieved in several stages or windows in two series of calculations – one in the protein and one in a box of water. This is the basis of alchemical free energy methods like the Free Energy Perturbation (FEP) and Thermodynamic Integration (TI)[8, 137, 310, 311]. In some cases, the size of the systems prohibits long MD simulations with all-atom force fields. Instead, these simulations can be performed at the desirable time scales with coarse grained force fields (e.g. MARTINI) which describe groups of atoms with specially parametrized pseudo-particles, leading to substantial reduction in the number of relevant particles[8]. Methods based on reduced representation of the protein (elastic network model, two-state anisotropic network model)[316] and importance sampling, such as dynamics importance sampling[317] or weighted ensemble path sampling[318]) are also applied for conformational studies. Other coarse grained approaches involve the Brownian dynamics (BD) method, based on an overdamped Langevin equation in which some components of the system (i.e. solvent, membrane) are included implicitly as random stochastic forces acting on the chemically relevant, explicitly included components (i.e. proteins, ions). Brownian dynamics allows for evaluation of the ion conductance in ion channels which take place at the µs scale and as such is still too expensive for standard all-atom MD. An external electric field mimicking the membrane potential is often applied to the model channel in BD (or MD) in order to induce ion flux through the pore. In these cases, both the ions and their surroundings are accounted for implicitly. The general equation of motion governing dynamics in BD simulations of an ion traversing an ion channel is a specific form of the general Langevin equation developed by Ermak and McCammon [319]: =−

GN 'N  89 :

∇ ` / ,  , ⋯  + ∇ = + b c

(A6)

where W(r1,r2,...) is a many-body PMF that describes interactions between ions and protein, the effect of applied membrane voltage, and the reaction field and static field emerging from protein charges. Di is a position-dependent diffusion coefficient which was kept to corresponding bulk diffusion values for both ions and nucleotides, and ξi(τ) is a term introducing Gaussian random noise to the system dynamics. Detailed comparison of MD, BD, PNP, and PB can be found in Ref.[137]. It is also possible to combine different scales e.g. MD and BD to compute nonequilibrium dynamics of the system. One can obtain a multi-particle PMF W(r1,r2,...) for an ion or solute with comprehensive sampling from MD[30, 139, 141] and then use it as an input for the equation above for the range of environmental conditions. 47

The MM methods cannot give explicit information about the electronic structure of the system, therefore they cannot be used for evaluation of excitations between different electronic levels (UV/VIS spectra) or complex electron redistribution within the system (i.e. atom and bond polarization, bond formation and dissociation). Some of these problems have been addressed with novel force fields which incorporate polarization effects (Drude[20], AMOEBA[320]) or can account for bond breaking/dissociation (e.g. ReaxFF[321]). Nevertheless, classical MM and MD provide abundant information about the structure of binding sites, possible ion pathways, energetics of binding and transport and dynamic structural reorganizations. They incorporate naturally solvent, temperature and pressure effects which makes their results easily relatable to experiment. Large amounts of solvent molecules can be added explicitly and realistic lipid membranes can also be simulated. In contrast, QM methods often need an implicit polarizable continuum description of the surroundings and can handle only a limited number of solvent molecules explicitly. A3. Computational Electrophysiology Ion conductance (see eq.A3) is of great importance in studies of ion channels. Unfortunately, its calculation is not a trivial process, since it requires the application of external electric fields or the introduction of ionic gradients in the presence of periodic boundary conditions. The computational approaches used in such studies (and the challenges associated with them) necessitate a more detailed discussion. As such, they have warranted their own section in this Appendix. Computational electrophysiology in a general sense is the method for simulation of ionic currents induced by voltages applied across transmembrane proteins. This method appears to be trivial if one simply applies a constant electric field across a periodic simulation box to drive ions through a conducting pore. To provide a deeper insight into this problem, Roux[322] discussed a rigorous mathematical framework for addressing simulations with applied electric fields that reproduce common electrophysiological experiments. There are two major issues related to the constant-electric field approach: (1) while the electric fields across realistic cellular membranes induced by ionic gradients may be strongly dependent on the geometries of the conducting pores at the sub-nanometer scale, simulations using constant electric fields impose invariant forces on ions through the entire systems regardless of the geometries; (2) even though the electrostatic potential profiles averaged over various simulation trajectories can capture the changes in the geometries[131], the simulation periodicity should be treated with care. The periodicity implemented in MD simulations serves to enhance sampling, thus reducing finite-size effects, however simulations using the same voltage but different constant electric fields and different system sizes yield non-identical ionic currents through the same conducting pore and afford distinguishable memory functions that underlie or even enhance the inherent finite-size effects, particularly at large voltages. As recommended by Gumbart et al[131], the voltages in the constant electric field simulations should be equal to the product of the electric fields (E) and the length of the simulation box (dC ) along the electric fields (∆f = dC ), rather than depending on the membrane thickness. It should be highlighted that this simulation approach can provide significant insight into the mechanism of ion permeation in voltage-gated potassium channels[19, 88, 89]. The major drawback of this approach is the frequent use of un-physiologically high voltages (large electric fields) to speed up the ion permeation within affordable simulation timescales. This results in an uncertain conductance that may not be reliable for comparison with experimental values. 48

Fig.A1 Schematic setup for computational electrophysiology. The conducting pore is shown as vertical purple lines embedded in the membrane. An alternative computational electrophysiology approach, was designed by Kutzner et al[134]. It can mimic better the electric fields in cellular membranes induced by ionic gradients or charge imbalance at the cost of doubling the size of the simulated system[323]. In the case of ion permeation through an ion channel embedded in a membrane, this single-membrane system is doubled along the z direction (Fig.A1) for simulation of the charge imbalance. As a result, there are two separate compartments of water and ions that are clearly divided by the two identical membranes (see Fig.A1). The charge imbalance between the two compartments is computed at every time step as follows: ∆3g = h$/? g3 ? + $/A g3 A i − h$? g3 ? + $A g3 A i

(A7)

? where $/, g is the number of identical negative charges 3 ? counted at time t in the A compartments (1 or 2); similarly, $/, g is the number of identical positive charges 3 A . In order to constrain the charge imbalance at a certain value during the simulations, excess ions in one compartment are exchanged with water molecules in the other compartment. All exchangeable particles are chosen within the thin central layer in each compartment to minimize the artificial disturbance and resistance caused by the exchange events. One can either exchange ions with water according to the above deterministic protocol or via Metropolis Monte Carlo criteria. As a result, a steady flux of ions driven by the charge imbalance constraint can be simulated with a voltage being computed as the charge imbalance divided by a membrane capacitance.

Note that the computational electrophysiology approach may still suffer from finite-size effects, which are an inherent issue in the MD simulations, particularly at large charge voltages. This method can produce stable voltages as low as 80 mV, but would return as many as 31 permeation events during 1.7 µs[61]. For collection of more statistically relevant data, other simulations[132-134, 194] using this approach have to increase the voltages (charge imbalances) to un-physiological values. The resulting conductances (~ 80 pS KcsA and 27 pS NavMs) appear to agree with experiment. However, as mentioned above most of the applied voltages in this type

49

of simulation remain so unphysiologically high, that the credibility of the results and the method are not convincing. A4. QM/MM QM/MM methods were invented as a way of combining the strong suits of the QM and MM approaches and compensating for their inherent weaknesses. In a QM/MM description a system is divided in (most often) two parts – a small QM part in which the electronic effects are of significance (i.e. a chemical reaction site or a chromophore) and a large MM part which provides a realistic description of the environment[8, 310, 311]. QM/MM methods are indispensable for studies of electron and proton transfer and biochemical reactions. Their computational expense depends strongly on the size of the QM part and the QM method used. In ion transport they are used most often to give an accurate binding description (complemented with charge transfer effects) and in evaluation of relative binding energies (similarly to FEP). The transport of ions, which require exchange of molecules between the MM and QM parts is still beyond their capabilities. Dynamical effects however can be probed with the so-called ab initio molecular dynamics methods, which combine the QM Hamiltonian description of the energetics of a system with equations of motion describing its evolution in time. Such methods (Car-Parinello MD[324]) are very costly and are currently applicable only to very small systems of ~100 atoms.

50

Na+

IC, mM 10

EC, mM 140

Representative channels Nav(NavAb, NavMs, NavRh, NaK, M2-TMD)

K+

140

4

0.0001

2

Kv (KcsA, KiR, SHAKER B, NaK, hERG, MthK), L-type Cav, Orai

1 4 10 7 10

1 104 27 1 1

Ca2+ Mg2+ ClHCO3SO42PO43-

MgtE, CorA CLC-ec1, MscS, MscL OprP

Representative transporters GltPh, LeuT, vSGLT, hSERT, hDAT, NahH, Mhp1, NCX_Mj, NaPi, Na,K-ATPase, EAAT3, vCNT Na,K-ATPase, H,K-ATPase, EAAT3, hSERT SERCA, NCX_Mj, NCKX, BtuB CLC-0, CLC-1, Band 3 Band 3, NBCn1/n2, NBCe1/e2 SBP GlpT

Table 1. Intra- and extracellular concentration (IC and EC, respectively) of the major inorganic ions found in the human body[2]. Included are some ion channels and transporters found in the referenced literature.

51

Figure Captions Fig.1 (Left) Pore and selectivity filter of KcsA (PDB: 1K4C, in orange and green are shown 2 of the 4 monomer components). (Right) Closer look at the SF of KcsA. The K+ binding sites (S0 – S4) and the involved amino acids are also indicated. The z axis shows the reaction coordinate used in PMF simulations of KcsA. Fig.2 Conformational changes (elevator-like motion) during active transport of Na+ ions and aspartate by GltPh (PDB: 2NWX). Fig.3 First coordination shell at the binding sites determined crystallographically for some of the transport systems considered in the review: S2 site in KcsA, one of the Na+ binding sites in GltPh, Ca2+ binding site in NCX_Mj, a peripheral Mg2+ binding site in MgtE, one of the Clbinding sites in CLC_ec1, and Pi binding site in OprP. The used PDB files are 1K4C (KcsA), 2NWX (GltPh), 3V5U (NCX_Mj), 2ZY9 (MgtE), 1OTS (CLC_ec1), 2O4V (OprP). The colors of the unlabeled atoms are as follows: red (oxygen), cyan (carbon), dark blue (nitrogen). Fig.4 (Left) 2D PMF of NavAb with two Na+ in the selectivity filter. (Right) 1D PMF profiles for NavAb calculated with either 1 or 2 (via integration of the 2D PMF profile) Na+ in the selectivity filter. The presence of a second Na+ alters drastically the energy landscape. Some structures corresponding to local minima in the 2D PMF plot (numbered 1-5) are shown as well to illustrate the permeation of the ions through the selectivity filter. Fig.5 Other systems transporting Na+ and K+ ions: a voltage-gated Na+ channel (NavAb, PDB: 4EKW, with two E residues from the EEEE ring displayed at the SHFS site), a monomer of the OmpF porin (PDB: 2OMF), Gramicidin A (gA, PDB: 1JNO), α-Hemolysin pore (αHL, PDB: 7AHL). Fig.6 Na+ coupled secondary transporters. (Left) LeuT dimer (PDB: 2A65) with bound Na+ (blue balls) and Leu (red balls). The orange and green helices represent the “bundle” and the purple helices, the dimerization “scaffold” involved in the rocking bundle mechanism. The yellow (TMs 1, 6) and cyan helices (TMs 3, 8) contain the binding sites Na1, Na2 and S1. In pink are shown residues associated with the putative binding site S2. (Right) A GltPh monomer (PDB: 2NWX) with bound Na+ ions (blue balls) and Asp (red balls). The rigid trimerization scaffold in shown in purple, HP1 – in green, and HP2 – in cyan. Fig.7 Ca2+ transport. (Left) SERCA in its E2 state (PDB: 3N5K). The cytoplasmic domains is color-coded – purple (A), green (P) and red (N). (Right) NCX_Mj (PDB: 3V5U). The positions of the Na+ and Ca2+ in the active site are shown with blue balls and cyan balls, respectively. Fig.8 Mg2+ transport. (Left) MgtE (PDB: 2ZY9) with indicated N (red) and CBS1,2 (green) and plug helices (blue). (Right) CorA (PDB: 2IUB) seen from the side and from the bottom. The positions of the GMN motif, basic and acidic rings are also included in the figure. A pair of

52

short-long TMs (TM2 and TM1) are shown in blue for one of the CorA chains. The Mg2+ in both MgtE and CorA are displayed as purple balls. Fig.9 Systems transporting Cl-. (Left) CLC-ec1 (PDB: 1OTS) transporter with bound 2Cl- (cyan balls). The position of the Cl- binding sites (Sext, Scen, and Sin) are indicated and some relevant residues are highlighted in blue (Gluex and Gluin) and red (Tyr 445 and Ser 107). The curvilinear arrows trace the separate pathways for Cl- and H+ ions. (Right) Mechanosensitive channel MscS (PDB: 2OAV). Fig.10 Inorganic phosphate transport. (Left) OprP porin (PDB: 2O4V). The phosphate ions are shown in yellow and the discussed Arg 133 and Asp 94 residues – in red and blue, respectively. (Right) Inward-facing structure of GlpT (1PW4).

53

Fig.1

Fig.2

54

Fig.3

Fig.4

55

Fig.5

Fig.6

56

Fig.7

57

Fig.8

58

Fig.9

Fig.10

59

References: [1] J.J.R.F. da Silva, R.J.P. Williams, The biological chemistry of the elements : the inorganic chemistry of life 2ed., Oxford University Press, USA, 2001. [2] R.B. Martin, Bioinorganic Chemistry, in: Encyclopedia of Molecular Cell Biology and Molecular Medicine, 2004. [3] J. Brajtburg, W.G. Powderly, G.S. Kobayashi, G. Medoff, Antimicrob Agents Chemother, 34 (1990) 183-188. [4] B.A. Wallace, Annu Rev Biophys Biophys Chem, 19 (1990) 127-157. [5] J.B. Kim, Korean J Pediatr, 57 (2014) 1-18. [6] O. Jardetzky, Nature, 211 (1966) 969-970. [7] M.V. LeVine, M.A. Cuendet, G. Khelashvili, H. Weinstein, Chem Rev, 116 (2016) 6552-6587. [8] E.N. Maldonado, K.L. Sheldon, D.N. DeHart, J. Patnaik, Y. Manevich, D.M. Townsend, S.M. Bezrukov, T.K. Rostovtseva, J.J. Lemasters, J Biol Chem, 288 (2013) 11920-11929. [9] Ion Channels: Methods and Protocols, 2 ed., Springer New York, 2013. [10] D.A. Doyle, J. Morais Cabral, R.A. Pfuetzner, A. Kuo, J.M. Gulbis, S.L. Cohen, B.T. Chait, R. MacKinnon, Science, 280 (1998) 69-77. [11] J. Aqvist, O. Alvarez, G. Eisenman, J Phys Chem-Us, 96 (1992) 10019-10025. [12] G. Eisenman, J. Aqvist, G. Appleby, O. Alvarez, Faseb J, 6 (1992) A144-A144. [13] P. Lauger, W. Stephan, E. Frehland, Biochim Biophys Acta, 602 (1980) 167-180. [14] K. Henzler-Wildman, D. Kern, Nature, 450 (2007) 964-972. [15] C. Zhao, S.Y. Noskov, PLoS Comput Biol, 9 (2013) e1003296. [16] J. Ostmeyer, S. Chakrapani, A.C. Pan, E. Perozo, B. Roux, Nature, 501 (2013). [17] K. Khafizov, C. Perez, C. Koshy, M. Quick, K. Fendler, C. Ziegler, L.R. Forrest, P Natl Acad Sci USA, 109 (2012) E3035-E3044. [18] R.O. Dror, R.M. Dirks, J.P. Grossman, H. Xu, D.E. Shaw, Annu Rev Biophys, 41 (2012) 429-452. [19] M.O. Jensen, V. Jogini, D.W. Borhani, A.E. Leffler, R.O. Dror, D.E. Shaw, Science, 336 (2012) 229-233. [20] J.A. Lemkul, J. Huang, B. Roux, A.D. MacKerell, Jr., Chem Rev, 116 (2016) 4983-5013. [21] M. Elstner, G. Seifert, Philos Trans A Math Phys Eng Sci, 372 (2014) 20120483. [22] H.M. Senn, W. Thiel, Angew Chem Int Ed Engl, 48 (2009) 1198-1229. [23] S. Lindert, D. Bucher, P. Eastman, V. Pande, J.A. McCammon, J Chem Theory Comput, 9 (2013) 46844691. [24] M.L. Laury, L.P. Wang, V.S. Pande, T. Head-Gordon, J.W. Ponder, J Phys Chem B, 119 (2015) 94239437. [25] X.C. Bai, G. McMullan, S.H.W. Scheres, Trends Biochem Sci, 40 (2015) 49-57. [26] E. Callaway, Nature, 525 (2015) 172-174. [27] D. Massiot, F. Fayon, M. Capron, I. King, S. Le Calve, B. Alonso, J.O. Durand, B. Bujoli, Z.H. Gan, G. Hoatson, Magn Reson Chem, 40 (2002) 70-76. [28] Y. Zhou, J.H. Morais-Cabral, A. Kaufman, R. MacKinnon, Nature, 414 (2001) 43-48. [29] B. Hille, Ion Channels of Excitable Membranes, 3rd ed., Sinauer Associates Inc., Sunderland, USA, 2001. [30] S.Y. Noskov, B. Roux, Biophys Chem, 124 (2006) 279-291. [31] Y. Zhou, R. MacKinnon, J Mol Biol, 333 (2003) 965-975. [32] A.L. Hodgkin, R.D. Keynes, J Physiol, 128 (1955) 61-88. [33] L.J. Mullins, J Gen Physiol, 42 (1959) 817-829.

60

[34] G. Eisenman, Biophys J, 2 (1962) 259-323. [35] F. Bezanilla, C.M. Armstrong, J Gen Physiol, 60 (1972) 588-608. [36] J.F. Danielli, H. Davson, J Cell Comp Physiol, 5 (1935) 495-508. [37] P.W. Fowler, K. Tai, M.S. Sansom, Biophys J, 95 (2008) 5062-5072. [38] C. Maffeo, S. Bhattacharya, J. Yoo, D. Wells, A. Aksimentiev, Chem Rev, 112 (2012) 6250-6284. [39] S.B. Long, X. Tao, E.B. Campbell, R. MacKinnon, Nature, 450 (2007) 376-382. [40] S.B. Long, E.B. Campbell, R. Mackinnon, Science, 309 (2005) 897-903. [41] Y. Jiang, A. Lee, J. Chen, V. Ruta, M. Cadene, B.T. Chait, R. MacKinnon, Nature, 423 (2003) 33-41. [42] Y. Jiang, A. Lee, J. Chen, M. Cadene, B.T. Chait, R. MacKinnon, Nature, 417 (2002) 515-522. [43] S. Choe, Nat Rev Neurosci, 3 (2002) 115-121. [44] Z. Lu, A.M. Klem, Y. Ramu, Nature, 413 (2001) 809-813. [45] M.P. Bhate, B.J. Wylie, L. Tian, A.E. McDermott, J Mol Biol, 401 (2010) 155-166. [46] B. Roux, R. MacKinnon, Science, 285 (1999) 100-102. [47] S.Y. Noskov, S. Berneche, B. Roux, Nature, 431 (2004) 830-834. [48] H. Yu, S.Y. Noskov, B. Roux, J Phys Chem B, 113 (2009) 8725-8730. [49] S.Y. Noskov, B. Roux, J Gen Physiol, 129 (2007) 135-143. [50] B. Roux, Biophys J, 98 (2010) 2877-2885. [51] D.L. Bostick, C.L. Brooks, Proc Natl Acad Sci U S A, 104 (2007) 9260-9265. [52] S. Durdagi, S.Y. Noskov, Channels (Austin), 5 (2011) 198-200. [53] A. Alam, Y.X. Jiang, Nat Struct Mol Biol, 16 (2009) 35-41. [54] A. Alam, Y.X. Jiang, Nat Struct Mol Biol, 16 (2009) 30-34. [55] I. Kim, T.W. Allen, Proc Natl Acad Sci U S A, 108 (2011) 17963-17968. [56] A.N. Thompson, I. Kim, T.D. Panosian, T.M. Iverson, T.W. Allen, C.M. Nimigean, Nat Struct Mol Biol, 16 (2009) 1317-1324. [57] V. Ngo, D. Stefanovski, S. Haas, R.A. Farley, PLoS One, 9 (2014) e86079. [58] I.H. Shrivastava, D.P. Tieleman, P.C. Biggin, M.S.P. Sansom, Biophys J, 83 (2002) 633-645. [59] A. Lange, K. Giller, S. Hornig, M.F. Martin-Eauclaire, O. Pongs, S. Becker, M. Baldus, Nature, 440 (2006) 959-962. [60] C.M. Nimigean, C. Miller, J Gen Physiol, 120 (2002) 323-335. [61] D.A. Kopfer, C. Song, T. Gruene, G.M. Sheldrick, U. Zachariae, B.L. de Groot, Science, 346 (2014) 352-355. [62] A.D. MacKerell, D. Bashford, M. Bellott, R.L. Dunbrack, J.D. Evanseck, M.J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F.T. Lau, C. Mattos, S. Michnick, T. Ngo, D.T. Nguyen, B. Prodhom, W.E. Reiher, B. Roux, M. Schlenkrich, J.C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin, M. Karplus, J Phys Chem B, 102 (1998) 3586-3616. [63] V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, C. Simmerling, Proteins, 65 (2006) 712-725. [64] A. Vaziri, M.B. Plenio, New J Phys, 12 (2010). [65] S. Berneche, B. Roux, Nature, 414 (2001) 73-77. [66] B. Egwolf, B. Roux, J Mol Biol, 401 (2010) 831-842. [67] H.T. Kratochvil, J.K. Carr, K. Matulef, A.W. Annen, H. Li, M. Maj, J. Ostmeyer, A.L. Serrano, H. Raghuraman, S.D. Moran, J.L. Skinner, E. Perozo, B. Roux, F.I. Valiyaveetil, M.T. Zanni, Science, 353 (2016) 1040-1044. [68] B. Roux, S. Berneche, B. Egwolf, B. Lev, S.Y. Noskov, C.N. Rowley, H. Yu, J Gen Physiol, 137 (2011) 415-426. [69] B. Roux, S. Berneche, W. Im, Biochemistry, 39 (2000) 13295-13306. [70] H. Yu, S.Y. Noskov, B. Roux, Proc Natl Acad Sci U S A, 107 (2010) 20329-20334. [71] P.W. Fowler, E. Abad, O. Beckstein, M.S. Sansom, J Chem Theory Comput, 9 (2013) 5176-5189.

61

[72] H. Li, V. Ngo, M.C. Da Silva, D.R. Salahub, K. Callahan, B. Roux, S.Y. Noskov, J Phys Chem B, 119 (2015) 9401-9416. [73] V. Ngo, M.C. da Silva, M. Kubillus, H. Li, B. Roux, M. Elstner, Q. Cui, D.R. Salahub, S.Y. Noskov, J Chem Theory Comput, 11 (2015) 4992-5001. [74] M.D. Fayer, N.E. Levinger, Annu Rev Anal Chem, Vol 3, 3 (2010) 89-107. [75] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J Chem Phys, 79 (1983) 926935. [76] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J Phys Chem, 91 (1987) 6269-6271. [77] P.T. Kiss, A. Baranyai, J Chem Phys, 134 (2011). [78] C. Vega, J.L.F. Abascal, M.M. Conde, J.L. Aragones, Farad Discuss, 141 (2009) 251-276. [79] S. Varma, S.B. Rempe, Biophys J, 99 (2010) 3394-3401. [80] Y. Zhou, R. MacKinnon, Biochemistry, 43 (2004) 4978-4982. [81] A. Kariev, M.E. Green, Biochim Biophys Acta Biomembr, 1788 (2009) 1188-1192. [82] D.M. Rogers, S.B. Rempe, J Phys Chem B, 115 (2011) 9116-9129. [83] D.M. Rogers, S.B. Rempe, J Phys Chem B, 116 (2012) 7994-7995. [84] B. Roux, J Phys Chem B, 116 (2012) 7991-7993. [85] L. Delemotte, M.A. Kasimova, M.L. Klein, M. Tarek, V. Carnevale, Proc Natl Acad Sci U S A, 112 (2015) 124-129. [86] L. Delemotte, M. Tarek, M.L. Klein, C. Amaral, W. Treptow, Proc Natl Acad Sci U S A, 108 (2011) 6109-6114. [87] L. Delemotte, W. Treptow, M.L. Klein, M. Tarek, Biophys J, 99 (2010) L72-74. [88] M.O. Jensen, D.W. Borhani, K. Lindorff-Larsen, P. Maragakis, V. Jogini, M.P. Eastwood, R.O. Dror, D.E. Shaw, Proc Natl Acad Sci U S A, 107 (2010) 5833-5838. [89] M.O. Jensen, V. Jogini, M.P. Eastwood, D.E. Shaw, J Gen Physiol, 141 (2013) 619-632. [90] W. Treptow, M. Tarek, Biophys J, 90 (2006) L64-66. [91] F. Khalili-Araghi, V. Jogini, V. Yarov-Yarovoy, E. Tajkhorshid, B. Roux, K. Schulten, Biophys J, 98 (2010) 2189-2198. [92] W.A. Catterall, Neuron, 26 (2000) 13-25. [93] A.F. Barber, V. Carnevale, S.G. Raju, C. Amaral, W. Treptow, M.L. Klein, Biochim Biophys Acta, 1818 (2012) 2120-2125. [94] C. Amaral, V. Carnevale, M.L. Klein, W. Treptow, Proc Natl Acad Sci U S A, 109 (2012) 21336-21341. [95] V. Yarov-Yarovoy, P.G. DeCaen, R.E. Westenbroek, C.Y. Pan, T. Scheuer, D. Baker, W.A. Catterall, Proc Natl Acad Sci U S A, 109 (2012) E93-102. [96] P.G. DeCaen, V. Yarov-Yarovoy, E.M. Sharp, T. Scheuer, W.A. Catterall, Proc Natl Acad Sci U S A, 106 (2009) 22498-22503. [97] P.G. DeCaen, V. Yarov-Yarovoy, Y. Zhao, T. Scheuer, W.A. Catterall, Proc Natl Acad Sci U S A, 105 (2008) 15142-15147. [98] A.C. Pan, L.G. Cuello, E. Perozo, B. Roux, J Gen Physiol, 138 (2011) 571-580. [99] C. Domene, A. Grottesi, M.S. Sansom, Biophys J, 87 (2004) 256-267. [100] C. Domene, M.S. Sansom, Biophys J, 85 (2003) 2787-2800. [101] S. Furini, C. Domene, Proc Natl Acad Sci U S A, 106 (2009) 16074-16077. [102] P. Kollman, Chem Rev, 93 (1993) 2395-2417. [103] G.M. Torrie, J.P. Valleau, J Comput Phys, 23 (1977) 187-199. [104] B. Roux, Comput Phys Commun, 91 (1995) 275-282. [105] S. Kumar, J.M. Rosenberg, D. Bouzida, R.H. Swendsen, P.A. Kollman, J Comput Chem, 13 (1992) 1011-1021. [106] S. Park, M. Jensen, E. Tajkhorshid, K. Schulten, Biophys J, 82 (2002) 487a-487a. [107] S. Park, F. Khalili-Araghi, E. Tajkhorshid, K. Schulten, J Chem Phys, 119 (2003) 3559-3566.

62

[108] S. Park, K. Schulten, J Chem Phys, 120 (2004) 5946-5961. [109] V.A. Ngo, I. Kim, T.W. Allen, S.Y. Noskov, J Chem Theory Comput, 12 (2016) 1000-1010. [110] D. Rodriguez-Gomez, E. Darve, A. Pohorille, J Chem Phys, 120 (2004) 3563-3578. [111] A. Laio, M. Parrinello, Proc Natl Acad Sci U S A, 99 (2002) 12562-12566. [112] Y. Sugita, Y. Okamoto, Chem Phys Lett, 314 (1999) 141-151. [113] C. Domene, S. Furini, Methods Enzymol, Vol 466: Biothermodyn, Pt B, 466 (2009) 155-177. [114] R.C. Bernardi, M.C. Melo, K. Schulten, Biochim Biophys Acta, 1850 (2015) 872-877. [115] B. Roux, Biophys J, 77 (1999) 139-153. [116] G.K. Schenter, B.C. Garrett, D.G. Truhlar, J Chem Phys, 119 (2003) 5828-5833. [117] E. Piccinini, M. Ceccarelli, F. Affinito, R. Brunetti, C. Jacoboni, J Chem Theory Comput, 4 (2008) 173183. [118] J.F. Gwan, A. Baumgaertner, J Chem Phys, 127 (2007) 045103. [119] S. Furini, C. Domene, J Mol Biol, 409 (2011) 867-878. [120] N. Chakrabarti, C. Ing, J. Payandeh, N. Zheng, W.A. Catterall, R. Pomes, Proc Natl Acad Sci U S A, 110 (2013) 11331-11336. [121] B. Corry, M. Thomas, J Am Chem Soc, 134 (2012) 1840-1846. [122] S. Furini, C. Domene, PLoS Comput Biol, 8 (2012) e1002476. [123] S. Furini, C. Domene, Biophys J, 105 (2013) 1737-1745. [124] H. Qiu, R. Shen, W. Guo, Biochim Biophys Acta, 1818 (2012) 2529-2535. [125] S. Furini, C. Domene, Biophys J, 101 (2011) 1623-1631. [126] S. Furini, C. Domene, Biophys J, 103 (2012) 2106-2114. [127] Y. Wang, A.C. Chamberlin, S.Y. Noskov, J Phys Chem B, 118 (2014) 2041-2049. [128] J. Neyton, C. Miller, J Gen Physiol, 92 (1988) 569-586. [129] J. Neyton, C. Miller, J Gen Physiol, 92 (1988) 549-567. [130] A. Aksimentiev, K. Schulten, Biophys J, 88 (2005) 3745-3761. [131] J. Gumbart, F. Khalili-Araghi, M. Sotomayor, B. Roux, Biochim Biophys Acta, 1818 (2012) 294-302. [132] C. Kutzner, D.A. Kopfer, J.P. Machtens, B.L. de Groot, C. Song, U. Zachariae, Biochim Biophys Acta Biomembr, 1858 (2016) 1741-1752. [133] U. Zachariae, C. Kutzner, H. Grubmuller, B. de Groot, Biophys J, 102 (2012) 22a-22a. [134] C. Kutzner, H. Grubmuller, B.L. de Groot, U. Zachariae, Biophys J, 101 (2011) 809-817. [135] I.S. Joung, T.E. Cheatham, J Phys Chem B, 112 (2008) 9020-9041. [136] B. Roux, Curr Opin Struct Biol, 12 (2002) 182-189. [137] B. Roux, T. Allen, S. Berneche, W. Im, Q Rev Biophys, 37 (2004) 15-103. [138] C.J. Solano, K.R. Pothula, J.D. Prajapati, P.M. De Biase, S.Y. Noskov, U. Kleinekathofer, J Chem Theory Comput, 12 (2016) 2401-2417. [139] J. Comer, A. Aksimentiev, J Phys Chem C, 116 (2012) 3376-3393. [140] P.M. De Biase, E.N. Ervin, P. Pal, O. Samoylova, S. Markosyan, M.G. Keehan, G.A. Barrall, S.Y. Noskov, Nanoscale, 8 (2016) 11571-11579. [141] S. Markosyan, P.M. De Biase, L. Czapla, O. Samoylova, G. Singh, J. Cuervo, D.P. Tieleman, S.Y. Noskov, Nanoscale, 6 (2014) 9006-9016. [142] B. Corry, S. Kuyucak, S.H. Chung, Biophys J, 84 (2003) 3594-3606. [143] S.Y. Noskov, W. Im, B. Roux, Biophys J, 87 (2004) 2299-2309. [144] V.M. Aguilella, M. Queralt-Martin, M. Aguilella-Arzo, A. Alcaraz, Integr Biol (Camb), 3 (2011) 159172. [145] T.W. Allen, S. Kuyucak, S.H. Chung, Biophys Chem, 86 (2000) 1-14. [146] W. Im, B. Roux, J Mol Biol, 322 (2002) 851-869. [147] S.H. Chung, T.W. Allen, S. Kuyucak, Biophys J, 83 (2002) 263-277. [148] A. Burykin, M. Kato, A. Warshel, Proteins, 52 (2003) 412-426.

63

[149] S. Berneche, B. Roux, Proc Natl Acad Sci U S A, 100 (2003) 8644-8648. [150] P. Graf, M.G. Kurnikova, R.D. Coalson, A. Nitzan, J Phys Chem B, 108 (2004) 2006-2015. [151] A.B. Mamonov, R.D. Coalson, A. Nitzan, M.G. Kurnikova, Biophys J, 84 (2003) 3646-3661. [152] A.B. Mamonov, M.G. Kurnikova, R.D. Coalson, Biophys Chem, 124 (2006) 268-278. [153] B. Corry, S.H. Chung, Cell Mol Life Sci, 63 (2006) 301-315. [154] W. Im, S. Seefeld, B. Roux, Biophys J, 79 (2000) 788-801. [155] B. Corry, T.W. Allen, S. Kuyucak, S.H. Chung, Biochim Biophys Acta, 1509 (2000) 1-6. [156] B. Corry, T.W. Allen, S. Kuyucak, S.H. Chung, Biophys J, 80 (2001) 195-214. [157] B. Corry, M. O'Mara, S.H. Chung, Biophys J, 86 (2004) 846-860. [158] T. Bastug, A. Gray-Weale, S.M. Patra, S. Kuyucak, Biophys J, 90 (2006) 2285-2296. [159] Y.Z. Tang, W.Z. Chen, C.X. Wang, Eur Biophys J, 29 (2000) 523-534. [160] B. Roux, Acc Chem Res, 35 (2002) 366-375. [161] D.H. Mackay, P.H. Berens, K.R. Wilson, A.T. Hagler, Biophys J, 46 (1984) 229-248. [162] T.W. Allen, T. Bastug, S. Kuyucak, S.H. Chung, Biophys J, 84 (2003) 2159-2168. [163] T.W. Allen, O.S. Andersen, B. Roux, J Gen Physiol, 124 (2004) 679-690. [164] T.W. Allen, O.S. Andersen, B. Roux, Proc Natl Acad Sci U S A, 101 (2004) 117-122. [165] M. Compoint, C. Ramseyer, P. Huetz, Chem Phys Lett, 397 (2004) 510-515. [166] D. Bucher, S. Raugei, L. Guidoni, M. Dal Peraro, U. Rothlisberger, P. Carloni, M.L. Klein, Biophys Chem, 124 (2006) 292-301. [167] D. Bucher, L. Guidoni, P. Carloni, U. Rothlisberger, Biophys J, 98 (2010) L47-49. [168] D. Bucher, L. Guidoni, U. Rothlisberger, Biophys J, 93 (2007) 2315-2324. [169] L. Guidoni, P. Carloni, Biochim Biophys Acta, 1563 (2002) 1-6. [170] L. Guidoni, V. Torre, P. Carloni, FEBS Lett, 477 (2000) 37-42. [171] M. Rossi, A. Tkatchenko, S.B. Rempe, S. Varma, Proc Natl Acad Sci U S A, 110 (2013) 12978-12983. [172] A. Warshel, M. Kato, A.V. Pisliakov, J Chem Theory Comput, 3 (2007) 2034-2045. [173] P.Y. Ren, J.W. Ponder, J Comput Chem, 23 (2002) 1497-1506. [174] J.W. Ponder, C.J. Wu, P.Y. Ren, V.S. Pande, J.D. Chodera, M.J. Schnieders, I. Haque, D.L. Mobley, D.S. Lambrecht, R.A. DiStasio, M. Head-Gordon, G.N.I. Clark, M.E. Johnson, T. Head-Gordon, J Phys Chem B, 114 (2010) 2549-2564. [175] S. Patel, J.E. Davis, B.A. Bauer, J Am Chem Soc, 131 (2009) 13890-13891. [176] N. Manin, M.C. da Silva, I. Zdravkovic, O. Eliseeva, A. Dyshin, O. Yasar, D.R. Salahub, A.M. Kolker, M.G. Kiselev, S.Y. Noskov, Phys Chem Chem Phys, 18 (2016) 4191-4200. [177] Q. Cui, J Chem Phys, 145 (2016) 140901. [178] D. Bucher, U. Rothlisberger, J Gen Physiol, 135 (2010) 549-554. [179] J. Payandeh, T. Scheuer, N. Zheng, W.A. Catterall, Nature, 475 (2011) 353-358. [180] C.E. Naylor, C. Bagneris, P.G. DeCaen, A. Sula, A. Scaglione, D.E. Clapham, B.A. Wallace, EMBO J, 35 (2016) 820-830. [181] X. Zhang, W. Ren, P. DeCaen, C. Yan, X. Tao, L. Tang, J. Wang, K. Hasegawa, T. Kumasaka, J. He, J. Wang, D.E. Clapham, N. Yan, Nature, 486 (2012) 130-134. [182] J. Payandeh, T.M. Gamal El-Din, T. Scheuer, N. Zheng, W.A. Catterall, Nature, 486 (2012) 135-139. [183] E.C. McCusker, C. Bagneris, C.E. Naylor, A.R. Cole, N. D'Avanzo, C.G. Nichols, B.A. Wallace, Nat Commun, 3 (2012) 1102. [184] C. Bagneris, P.G. DeCaen, C.E. Naylor, D.C. Pryde, I. Nobeli, D.E. Clapham, B.A. Wallace, Proc Natl Acad Sci U S A, 111 (2014) 8428-8433. [185] M.B. Ulmschneider, C. Bagneris, E.C. McCusker, P.G. Decaen, M. Delling, D.E. Clapham, J.P. Ulmschneider, B.A. Wallace, Proc Natl Acad Sci U S A, 110 (2013) 6364-6369. [186] X. Zhang, M. Xia, Y. Li, H. Liu, X. Jiang, W. Ren, J. Wu, P. DeCaen, F. Yu, S. Huang, J. He, D.E. Clapham, N. Yan, H. Gong, Cell Res, 23 (2013) 409-422.

64

[187] J.R. Tyson, T.P. Snutch, WIREs Membr Transport Signalling, 2 (2013) 181-225. [188] A. Khan, L. Romantseva, A. Lam, G. Lipkind, H.A. Fozzard, J Physiol, 543 (2002) 71-84. [189] D.W. Smith, J Chem Educ, 54 (1977) 540-542. [190] L. Stock, L. Delemotte, V. Carnevale, W. Treptow, M.L. Klein, J Phys Chem B, 117 (2013) 3782-3789. [191] C. Boiteux, I. Vorobyov, T.W. Allen, Proc Natl Acad Sci U S A, 111 (2014) 3454-3459. [192] R.K. Finol-Urdaneta, Y. Wang, A. Al-Sabi, C. Zhao, S.Y. Noskov, R.J. French, J. Gen. Physiol., 143 (2014) 157-171. [193] V. Ngo, W. Yibo, S. Haas, S.Y. Noskov, R.A. Farley, PLoS Comput Biol, 1004482 (2016) 1-19. [194] S. Ke, E.-M. Zangeri, A. Stary-Weinzinger, Biochem Biophys Res Commun, 430 (2013) 1272-1276. [195] B. Corry, PeerJ, 1 (2013) e16. [196] D.J. Ren, B. Navarro, H.X. Xu, L.X. Yue, Q. Shi, D.E. Clapham, Science, 294 (2001) 2372-2375. [197] T. Dudev, C. Lim, J Am Chem Soc, 136 (2014) 3553-3559. [198] T. Dudev, C. Lim, J Am Chem Soc, 132 (2010) 2321-2332. [199] T. Dudev, C. Lim, Phys Chem Chem Phys, 14 (2012) 12451-12456. [200] T. Dudev, C. Lim, Acc Chem Res, 47 (2014) 3580-3587. [201] S. Ke, E.N. Timin, A. Stary-Weinzinger, PLoS Comput Biol, 10 (2014) e1003746. [202] C. Ing, R. Pomes, Na Channels Phyla Funct, 78 (2016) 215-260. [203] A.M. Woodhull, J Gen Physiol, 61 (1973) 687-708. [204] C.M. Armstrong, G. Cota, Proc Natl Acad Sci U S A, 88 (1991) 6528-6531. [205] C.M. Armstrong, G. Cota, Proc Natl Acad Sci U S A, 96 (1999) 4154-4157. [206] W. Kopec, B. Loubet, H. Poulsen, H. Khandelia, Biochemistry, 53 (2014) 746-754. [207] V. Leone, D. Pogoryelov, T. Meier, J.D. Faraldo-Gomez, Proc Natl Acad Sci U S A, 112 (2015) E10571066. [208] K. Munson, R.J. Law, G. Sachs, Biochemistry, 46 (2007) 5398-5417. [209] I. Zdravkovic, C. Zhao, B. Lev, J.E. Cuervo, S.Y. Noskov, Biochim Biophys Acta, 1818 (2012) 337-347. [210] C.J. Loland, Biochim Biophys Acta, 1850 (2015) 500-510. [211] J. Grouleff, L.K. Ladefoged, H. Koldso, B. Schiott, Front Pharmacol, 6 (2015) 235. [212] J. Li, P.C. Wen, M. Moradi, E. Tajkhorshid, Curr Opin Struct Biol, 31 (2015) 96-105. [213] A. Yamashita, S.K. Singh, T. Kawate, Y. Jin, E. Gouaux, Nature, 437 (2005) 215-223. [214] N. Akyuz, E.R. Georgieva, Z. Zhou, S. Stolzenberg, M.A. Cuendet, G. Khelashvili, R.B. Altman, D.S. Terry, J.H. Freed, H. Weinstein, O. Boudker, S.C. Blanchard, Nature, 518 (2015) 68-73. [215] D. Yernool, O. Boudker, Y. Jin, E. Gouaux, Nature, 431 (2004) 811-818. [216] G. Verdon, S. Oh, R.N. Serio, O. Boudker, Elife, 3 (2014) e02283. [217] S. Manepalli, C.K. Surratt, J.D. Madura, T.L. Nolan, AAPS J, 14 (2012) 820-831. [218] P.J. Focke, X. Wang, H.P. Larsson, Structure, 21 (2013) 694-705. [219] G. Grazioso, V. Limongelli, D. Branduardi, E. Novellino, C. De Micheli, A. Cavalli, M. Parrinello, J Am Chem Soc, 134 (2012) 453-463. [220] Z. Huang, E. Tajkhorshid, Biophys J, 95 (2008) 2292-2300. [221] I.H. Shrivastava, J. Jiang, S.G. Amara, I. Bahar, J Biol Chem, 283 (2008) 28680-28690. [222] H.P. Larsson, X. Wang, B. Lev, I. Baconguis, D.A. Caplan, N.P. Vyleta, H.P. Koch, A. Diez-Sampedro, S.Y. Noskov, Proc Natl Acad Sci U S A, 107 (2010) 13912-13917. [223] L. Shi, M. Quick, Y. Zhao, H. Weinstein, J.A. Javitch, Mol Cell, 30 (2008) 667-677. [224] M. Gabrielsen, A.W. Ravna, K. Kristiansen, I. Sylte, J Mol Model, 18 (2012) 1073-1085. [225] L. Celik, S. Sinning, K. Severinsen, C.G. Hansen, M.S. Moller, M. Bols, O. Wiborg, B. Schiott, J Am Chem Soc, 130 (2008) 3853-3865. [226] A.M. Jorgensen, L. Tagmose, A.M. Jorgensen, K.P. Bogeso, G.H. Peters, ChemMedChem, 2 (2007) 827-840. [227] C.F. Zhao, S. Stolzenberg, L. Gracia, H. Weinstein, S. Noskov, L. Shi, Biophys J, 103 (2012) 878-888.

65

[228] E. Zomot, M. Gur, I. Bahar, J Biol Chem, 290 (2015) 544-555. [229] S. Stolzenberg, M. Quick, C.F. Zhao, K. Gotfryd, G. Khelashvili, U. Gether, C.J. Loland, J.A. Javitch, S. Noskov, H. Weinstein, L. Shi, J Biol Chem, 290 (2015) 13992-14003. [230] S.Y. Noskov, B. Roux, J Mol Biol, 377 (2008) 804-818. [231] J. Li, S.A. Shaikh, G. Enkavi, P.C. Wen, Z. Huang, E. Tajkhorshid, Proc Natl Acad Sci U S A, 110 (2013) 7696-7701. [232] M.V. LeVine, H. Weinstein, PLOS Comput Biol, 10 (2014) e1003603. [233] A.C. Pan, D. Sezer, B. Roux, J Phys Chem B, 112 (2008) 3432-3440. [234] M.C. Pinto, A.H. Kihara, V.A. Goulart, F.M. Tonelli, K.N. Gomes, H. Ulrich, R.R. Resende, Cell Signal, 27 (2015) 2139-2149. [235] H. Zhekova, C. Zhao, P.P. Schnetkamp, S.Y. Noskov, Biochemistry, 55 (2016) 6445-6455. [236] W.A. Catterall, Cold Spring Harb Perspect Biol, 3 (2011) a003947. [237] A. Malasics, D. Gillespie, W. Nonner, D. Henderson, B. Eisenberg, D. Boda, Biochim Biophys Acta, 1788 (2009) 2471-2480. [238] M. Bublitz, M. Musgaard, H. Poulsen, L. Thogersen, C. Olesen, B. Schiott, J.P. Morth, J.V. Moller, P. Nissen, J Biol Chem, 288 (2013) 10759-10765. [239] L.M. Espinoza-Fonseca, J.M. Autry, D.D. Thomas, PLoS One, 9 (2014) e95979. [240] L.M. Espinoza-Fonseca, J.M. Autry, G.L. Ramirez-Salinas, D.D. Thomas, Biophys J, 108 (2015) 16971708. [241] L.M. Espinoza-Fonseca, J.M. Autry, D.D. Thomas, Biochem Biophys Res Commun, 463 (2015) 37-41. [242] L.M. Espinoza-Fonseca, D.D. Thomas, PLoS One, 6 (2011) e26936. [243] J.E. Fonseca, S. Kaya, R.F. Rakowski, Nanotechnology, 18 (2007) 424022. [244] Y. Huang, H. Li, Y. Bu, J Comput Chem, 30 (2009) 2136-2145. [245] M. Musgaard, L. Thogersen, B. Schiott, E. Tajkhorshid, Biophys J, 102 (2012) 268-277. [246] Y. Sugita, M. Ikeguchi, C. Toyoshima, Proc Natl Acad Sci U S A, 107 (2010) 21465-21469. [247] Y. Sugita, N. Miyashita, M. Ikeguchi, A. Kidera, C. Toyoshima, J Am Chem Soc, 127 (2005) 61506151. [248] F. Xiang, R.I. Cukier, Y. Bu, J Phys Chem B, 111 (2007) 12282-12293. [249] C. Kobayashi, R. Koike, M. Ota, Y. Sugita, Proteins, 83 (2015) 746-756. [250] M. Giladi, I. Tal, D. Khananshvili, Front Physiol, 7 (2016) 30. [251] J. Liao, H. Li, W. Zeng, D.B. Sauer, R. Belmares, Y. Jiang, Science, 335 (2012) 686-690. [252] J. Liao, F. Marinelli, C. Lee, Y. Huang, J.D. Faraldo-Gomez, Y. Jiang, Nat Struct Mol Biol, 23 (2016) 590-599. [253] F. Marinelli, L. Almagor, R. Hiller, M. Giladi, D. Khananshvili, J.D. Faraldo-Gomez, Proc Natl Acad Sci U S A, 111 (2014) E5354-5362. [254] J. Payandeh, R. Pfoh, E.F. Pai, Biochim Biophys Acta, 1828 (2013) 2778-2792. [255] J. Sahni, A.M. Scharenberg, Mol Aspects Med, 34 (2013) 620-628. [256] R. Ishitani, Y. Sugita, N. Dohmae, N. Furuya, M. Hattori, O. Nureki, Proc Natl Acad Sci U S A, 105 (2008) 15393-15398. [257] R. Pfoh, A. Li, N. Chakrabarti, J. Payandeh, R. Pomes, E.F. Pai, Proc Natl Acad Sci U S A, 109 (2012) 18809-18814. [258] O. Dalmas, L.G. Cuello, V. Jogini, D.M. Cortes, B. Roux, E. Perozo, Structure, 18 (2010) 868-878. [259] T. Arakawa, T. Kobayashi-Yurugi, Y. Alguel, H. Iwanari, H. Hatae, M. Iwata, Y. Abe, T. Hino, C. IkedaSuno, H. Kuma, D. Kang, T. Murata, T. Hamakubo, A.D. Cameron, T. Kobayashi, N. Hamasaki, S. Iwata, Science, 350 (2015) 680-684. [260] A. Accardi, A. Picollo, Biochim Biophys Acta, 1798 (2010) 1457-1464. [261] R. Peyronnet, D. Tran, T. Girault, J.M. Frachisse, Front Plant Sci, 5 (2014) 558.

66

[262] S.Y. Noskov, T.K. Rostovtseva, A.C. Chamberlin, O. Teijido, W. Jiang, S.M. Bezrukov, Biochim Biophys Acta, 1858 (2016) 1778-1790. [263] R. Dutzler, E.B. Campbell, M. Cadene, B.T. Chait, R. MacKinnon, Nature, 415 (2002) 287-294. [264] R. Dutzler, E.B. Campbell, R. MacKinnon, Science, 300 (2003) 108-112. [265] J.D. Faraldo-Gomez, B. Roux, J Mol Biol, 339 (2004) 981-1000. [266] G.V. Miloshevsky, P.C. Jordan, Biophys J, 86 (2004) 825-835. [267] W. Han, R.C. Cheng, M.C. Maduke, E. Tajkhorshid, Proc Natl Acad Sci U S A, 111 (2014) 1819-1824. [268] F.L. Gervasio, M. Parrinello, M. Ceccarelli, M.L. Klein, J Mol Biol, 361 (2006) 390-398. [269] Y.J. Ko, W.H. Jo, J. Comp. Chem., 31 (2010) 603-611. [270] A. Suenaga, J.Z. Yeh, M. Taiji, A. Toyama, H. Takeuchi, M. Son, K. Takayama, M. Iwamoto, I. Sato, T. Narahashi, A. Konagaya, K. Goto, Biophys Chem, 120 (2006) 36-43. [271] Y. Zhang, G.A. Voth, Biophys J, 101 (2011) L47-49. [272] J. Cohen, K. Schulten, Biophys J, 86 (2004) 836-845. [273] M.H. Cheng, R.D. Coalson, Biophys J, 102 (2012) 1363-1371. [274] G. Kieseritzky, E.W. Knapp, J Biol Chem, 286 (2011) 2976-2986. [275] S. Pezeshki, C. Davis, A. Heyden, H. Lin, J Chem Theory Comput, 10 (2014) 4765-4776. [276] Z. Chen, T.L. Beck, J Phys Chem B, 120 (2016) 3129-3139. [277] A. Anishkin, K. Kamaraju, S. Sukharev, J Gen Physiol, 132 (2008) 67-83. [278] A. Anishkin, S. Sukharev, Biophys J, 86 (2004) 2883-2895. [279] M. Sotomayor, K. Schulten, Biophys J, 87 (2004) 3050-3065. [280] M. Sotomayor, T.A. van der Straaten, U. Ravaioli, K. Schulten, Biophys J, 90 (2006) 3496-3510. [281] M. Sotomayor, V. Vasquez, E. Perozo, K. Schulten, Biophys J, 92 (2007) 886-902. [282] R.B. Bass, P. Strop, M. Barclay, D.C. Rees, Science, 298 (2002) 1582-1587. [283] V. Vasquez, M. Sotomayor, J. Cordero-Morales, K. Schulten, E. Perozo, Science, 321 (2008) 12101214. [284] V. Vasquez, M. Sotomayor, D.M. Cortes, B. Roux, K. Schulten, E. Perozo, J Mol Biol, 378 (2008) 5570. [285] C. Song, B. Corry, PLoS One, 6 (2011) e21204. [286] L. Monticelli, S.K. Kandasamy, X. Periole, R.G. Larson, D.P. Tieleman, S.J. Marrink, J Chem Theory Comput, 4 (2008) 819-834. [287] M. Louhivuori, H.J. Risselada, E. van der Giessen, S.J. Marrink, Proc Natl Acad Sci U S A, 107 (2010) 19856-19860. [288] S. Yefimov, E. van der Giessen, P.R. Onck, S.J. Marrink, Biophys J, 94 (2008) 2994-3002. [289] Y. Huang, M.J. Lemieux, J. Song, M. Auer, D.N. Wang, Science, 301 (2003) 616-620. [290] T.F. Moraes, M. Bains, R.E. Hancock, N.C. Strynadka, Nat Struct Mol Biol, 14 (2007) 85-87. [291] G. Enkavi, J. Li, P. Mahinthichaichan, P.C. Wen, Z. Huang, S.A. Shaikh, E. Tajkhorshid, Methods Mol Biol, 924 (2013) 361-405. [292] G. Enkavi, E. Tajkhorshid, Biochemistry, 49 (2010) 1105-1114. [293] C.J. Law, J. Almqvist, A. Bernstein, R.M. Goetz, Y. Huang, C. Soudant, A. Laaksonen, S. Hovmoller, D.N. Wang, J Mol Biol, 378 (2008) 828-839. [294] N. Modi, R. Benz, R.E. Hancock, U. Kleinekathofer, J Phys Chem Lett, 3 (2012) 3639-3645. [295] N. Modi, I. Barcena-Uribarri, M. Bains, R. Benz, R.E. Hancock, U. Kleinekathofer, ACS Chem Biol, 10 (2015) 441-451. [296] N. Modi, I. Barcena-Uribarri, M. Bains, R. Benz, R.E. Hancock, U. Kleinekathofer, Biochemistry, 52 (2013) 5522-5532. [297] T. Kogai, G.A. Brent, Pharmacol Ther, 135 (2012) 355-370. [298] S. Micali, S. Bulotta, C. Puppin, A. Territo, M. Navarra, G. Bianchi, G. Damante, S. Filetti, D. Russo, BMC Cancer, 14 (2014) 303.

67

[299] R.B. Stockbridge, L. Kolmakova-Partensky, T. Shane, A. Koide, S. Koide, C. Miller, S. Newstead, Nature, 525 (2015) 548-551. [300] A.M. Pajor, Pflugers Arch, 451 (2006) 597-605. [301] J.W. Pflugrath, F.A. Quiocho, Nature, 314 (1985) 257-260. [302] T. Dudev, C. Lim, J Am Chem Soc, 126 (2004) 10296-10305. [303] F. Qiu, A. Chamberlin, B.M. Watkins, A. Ionescu, M.E. Perez, R. Barro-Soria, C. Gonzalez, S.Y. Noskov, H.P. Larsson, Proc Natl Acad Sci U S A, 113 (2016) E5962-E5971. [304] T. Dudev, C. Lim, Chem Rev, 103 (2003) 773-788. [305] N.W. Kelley, V. Vishal, G.A. Krafft, V.S. Pande, J Chem Phys, 129 (2008). [306] J.K. Weber, V.S. Pande, J Chem Theory Comput, 7 (2011) 3405-3411. [307] J.K. Weber, V.S. Pande, J Chem Theory Comput, 11 (2015) 2412-2420. [308] D. Shukla, C.X. Hernandez, J.K. Weber, V.S. Pande, Accounts Chem Res, 48 (2015) 414-422. [309] C.R. Schwantes, D. Shukla, V.S. Pande, Biophys J, 110 (2016) 1716-1719. [310] F. Jensen, Introduction to Computational Chemistry, John Wiley and sons, Ltd. UK 2006. [311] A.R. Leach, Molecular Modelling: Principles and Applications, 2 ed., Pearson Education Ltd, 2003. [312] D. Salahub, A. De la Lande, A. Goursot, R. Zhang, Y. Zhang, Struct Bond, 150 (2013). [313] S. Kumar, J.M. Rosenberg, D. Bouzida, R.H. Swendsen, P.A. Kollman, J Comput Chem, 16 (1995) 1339-1350. [314] J. Schlitter, M. Engels, P. Kruger, J Mol Graph, 12 (1994) 84-89. [315] X. Wu, B.R. Brooks, Chem Phys Lett, 381 (2003) 512-518. [316] A.R. Atilgan, S.R. Durell, R.L. Jernigan, M.C. Demirel, O. Keskin, I. Bahar, Biophys J, 80 (2001) 505515. [317] J.R. Perilla, O. Beckstein, E.J. Denning, T.B. Woolf, J Comput Chem, 32 (2011) 196-209. [318] P.G. Bolhuis, D. Chandler, C. Dellago, P.L. Geissler, Annu Rev Phys Chem, 53 (2002) 291-318. [319] D.L. Ermak, J.A. Mccammon, J Chem Phys, 69 (1978) 1352-1360. [320] Y. Shi, Z. Xia, J. Zhang, R. Best, C. Wu, J.W. Ponder, P. Ren, J Chem Theory Comput, 9 (2013) 40464063. [321] T.P. Senftle, S. Hong, M.M. Islam, S.B. Kylasa, Y. Zheng, Y.K. Shin, C. Junkermeier, R. Engel-Serbert, M.J. Janik, H.M. Aktulga, T. Verstraelen, A. Grama, A.C.T. van Duin, npj Comput Materials, 2 (2016) 15011. [322] B. Roux, Biophys J, 95 (2008) 4205-4216. [323] J.N. Sachs, P.S. Crozier, T.B. Woolf, J Chem Phys, 121 (2004) 10847-10851. [324] J. Hutter, WIREs Membr Transport Signalling, 2 (2012) 604-612.

68

Highlights • • • • •

Comparison of several competing theories for selective ion transport Review of computational studies on more than 10 small inorganic ions Overview and evaluation of available theoretical methods for ion transport studies Suggestion of ion transport areas that need further investigation Tabulated references in Supplementary Information

69

70