Validity of the Electrodiffusion Model for Calculating Conductance of

Subscriber access provided by RYERSON UNIVERSITY

Article

The Validity of the Electrodiffusion Model for Calculating Conductance of Simple Ion Channels Andrew Pohorille, Michael A. Wilson, and Chenyu Wei J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b09598 • Publication Date (Web): 23 Nov 2016 Downloaded from http://pubs.acs.org on December 1, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The Validity of the Electrodiusion Model for Calculating Conductance of Simple Ion Channels Andrew Pohorille,

∗,†,‡,§

Michael A. Wilson,

∗,†,¶

and Chenyu Wei

∗,†,‡

†Exobiology Branch, MS 239-4,

NASA Ames Research Center, Moett Field, CA 94035 ‡Department of Pharmaceutical Chemistry

University of California, San Francisco, San Francisco, CA 94132 ¶SETI Institute,

189 N Bernardo Ave #200, Mountain View, CA 94043 §To whom correspondences should be addressed. E-mail: [email protected]; [email protected]; [email protected] Phone: 650-604-5759. Fax: 650-604-1088

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 39

November 22, 2016

Abstract We examine the validity and utility of the electrodiusion (ED) equation, i.e. the generalized Nernst-Planck equation, to characterize, in combination with molecular dynamics, the electrophysiological behavior of simple ion channels. As models, we consider three systems - two naturally occurring channels formed by α-helical bundles of peptaibols, trichotoxin and alamethicin, and a synthetic, hexameric channel, formed by a peptide that contains only leucine and serine. All these channels mediate transport of potassium and chloride ions. Starting with equilibrium properties, such as the potential of mean force experienced by an ion traversing the channel and diusivity, obtained from molecular dynamics simulations, the ED equation can be used to determine the full current-voltage dependence with modest or no additional eort. The potential of mean force can be obtained not only from equilibrium simulations, but also, with comparable accuracy, from non-equilibrium simulations at a single voltage. The main assumptions underlying the ED equation appear to hold well for the channels and voltages studied here. To expand the utility of the ED equation, we examine what are the necessary and sucient conditions for Ohmic and non-rectifying behavior and relate deviations from this behavior to the shape of the ionic potential of mean force.

Introduction In their pioneering study, Aksimentiev and Schulten carried out explicit, molecular dynamics (MD) simulations of ion transport through a naturally occurring transmembrane ion channel,

α-hemolysin 1 .

This brought computational studies of ion transport outside the traditional

framework of the Poisson-Nernst-Planck equation or Brownian dynamics

26

into the realm

of all-atom simulations in which system dynamics and non-electrostatic forces acting on ions were accounted for. Simulations of ion transport through

α-hemolysin

were preceded

by studies that laid the groundwork for molecular dynamics of membrane systems in the

2


Page 3 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


presence of electric eld

7,8

and were followed by a number of studies in which ion current

through dierent channels was calculated directly from MD simulations

916

. The connection

between dierent theories of ion permeation and the corresponding MD simulations has been reviewed by Roux et al.

17

.

Independently, there has been a considerable body of work on calculating the potential of mean force (PMF) acting on ions along dierent channels in the absence of applied potential (for a recent review see Pohorille

18

). The PMF is the main ingredient of many molecular the-

ories of ion conductance. In particular, it enters the generalized Nernst-Planck equation for the ux of ions in which the gradient of the electrostatic potential is replaced by the average total force experienced by an ion. This equation is often called the generalized electrodiffusion (ED) equation or, simply, the electrodiusion equation.

From another perspective,

the ED equation can be considered a simplied Fokker-Planck equation or a form of the steady-state Smoluchowski equation. The connection between the PMF and ionic current highlighted in the ED equation has been exploited to estimate ionic currents through several channels

12,1922

.

An important advantage of the ED equation is its simplicity. The only ingredient that is dicult to determine and usually requires extensive simulations is the PMF. Other quantities are markedly easier to estimate. There is, however, a price to pay for this simplicity. The equation is based on a number of assumptions, the validity of which is sometimes unclear. For example, as is usually the case in most diusion models, it is assumed that events of interest (here, ions crossing the channel) are statistically independent.

It is also assumed

that the force acting on an ion is a sum of the force due to the applied voltage and the intrinsic, voltage-independent average force acting on a single ion in one dimension. Most of the assumptions underlying the ED equation have not been systematically tested. In this paper, we undertake this task with regard to the one-dimensional ED equation, using as models three simple channels. from fungi: trichotoxin

23,24

Two of them, are built of natural, non-genomic peptaibols

and alamethicin.

2527

The third one is a synthetic channel, LS3,

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

built of two types of amino acids, leucine and serine to mediate transport of K

+

− and Cl .

28,29

.

Page 4 of 39

All three channels are known

For these channels, we carried out extensive MD

simulations in the presence and absence of voltage and analyzed the results in the framework of the ED equation.

Our goal was not only to test the validity of the ED equation but

also to help establish the methodology for using this equation eciently and accurately in combination with MD simulations. In this paper, we do not address the comparisons between the calculated and experimental structures and electrophysiology of these channels. This has been

16

or will be done (Wilson and Pohorille, and Wei and Pohorille, unpublished) separately.

This separation is done on purpose, as problems involved in assessing accuracy of simulations compared with experiments, such as the adequacy of time scales and the reliability of force elds, are quite dierent from problems addressed here and have little bearing on the validity of the ED equation. We only mention that the calculated electrophysiological properties of the channels discussed here are in qualitative and sometimes quantitative agreement with experimental results. The ED equation, if valid, can be of considerable utility. It allows for calculating currents from the quantities obtained in equilibrium simulations with only modest, additional computational eort. It provides a means to extrapolate currents obtained at one applied voltage to other voltages, a task that is otherwise dicult to perform reliably

12,14,30,31

. In fact, the

ED equation, in combination with MD simulations at a single applied voltage, can be used to calculate eciently the full current-voltage (I-V) dependence, which is often used to characterize the electrophysiological properties of channels. This information is quite valuable for evaluating structural models of ion channels. If electrophysiological properties calculated for a given structure dier from measured properties beyond uncertainties expected for MD simulations and the ED equation then the model is most likely incorrect. This approach has been already successfully used to determine the number of helices in a channel formed by a peptaibol, antiamoebin

12

. Another interesting application of the ED equation is to solve an

inverse problem of reconstructing the equilibrium PMF from the current and non-equilibrium

4


Page 5 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ion density prole

16

. Finally, the ED equation can be used to gain general information about

electrophysiological characteristics of channels that behave in accord with this equation. After describing the methods for carrying out MD simulations of the model channels, we give a brief exposition of the ED equation and test whether ion crossing events are statistically independent, as required for the validity of the equation. Next, we systematically discuss all quantities that enter the ED equation. First, we describe the procedure for reconstructing the PMF from non-equilibrium simulations and compare the results with the PMFs obtained directly from the MD simulations.

We apply a similar procedure to reconstruct

voltage across the membrane in order to test the assumption about its linear dependence on the position along the bilayer normal.

Then, we proceed to the reconstruction of the

ion density proles, discussion of methods for calculating diusivities and the comparison between currents calculated from MD and the ED model. Finally discuss relations between general electrophysiological characteristics of ion channels and the intrinsic PMF that can be extracted from the ED equation. The proofs of these relations are included in the Supplementary Material. The paper concludes with a summary.

Methods For all three channels, trichotoxin, LS3 and alamethicin, molecular dynamics simulations were carried out in systems consisting of a planar phospholipid bilayer in contact with aqueous lamella on each side, a protein bundle in the transmembrane orientation and KCl. Most simulations were carried out on the Anton computer at Pittsburgh Supercomputer Center

32

.

Description of model building and MD simulations of the trichotoxin channel has been described previously

16

. Pictures of equilibrated transmembrane systems are shown in Fig. S1

in the Supplementary material. The LS3 peptide contains 21 residues in a heptad repeat (LSSLLSL) 3 blocked with the acetyl and N-methylamide groups at the N- and C-terminus, respectively

5


28,29

.

In an

α-


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 39

helical conformation, the peptide has a hydrophobic, leucine-rich face and a hydrophilic, serine-rich face. When an electric potential is applied across the membrane, LS3 monomers form long-lived, voltage-gated hexameric channels in which the hydrophilic faces line the ionconducting pore

28,29,33,34

Crick coiled coil bundle

. The initial channel model was built using parameters of an ideal

35

.

The bundle was placed in a 1-palmitoyl-2-oleoyl-sn-glycero-3-

phosphocholine (POPC) membrane. Besides the channel, the system contained 16228 water molecules, 218 POPC molecules, 301 K

+

and 301 Cl

−

ions, which corresponded to a

M KCl concentration.The size of the simulation box was 85.283 Å

×

85.154 Å

Å with periodic boundary conditions applied in all three spatial direction. forces were calculated by way of the particle-mesh Ewald method using 64

×

1

107.670

Long ranged

×

64

×

64 k-

vectors. The equations of motion were integrated using the r-RESPA multiple time-stepping algorithm with a base time step of 2 fs, and the long-ranged forces were updated every 6 fs. Temperature of 300 K was maintained using the Berendsen NVT thermostat. The CHARMM36 force eld was used to describe the POPC phospholipids molecules were represented by the TIP3P model CHARMM22 with CMAP corrections

38,39

37

36

, and water

. The LS3 peptide was described using

. All bonds between heavy atoms and hydrogen

atoms were constrained with the aid of the SHAKE algorithm. Initially, a MD trajectory 3122 ns in length was generated with an applied voltage of 200 mV. After approximately 1300 ns, the structure of the bundle shifted into a dierent structure, which remained stable in all simulations thereafter.

Consequently, the initial

1300 ns of the trajectory were not considered in our analyzes. Three more trajectories were generated, at 100, -100 and 0 mV. Their lengths were, respectively, 3026 ns, 2753 ns and 960 ns. The ion density proles obtained in the absence of voltage were used to calculate the free energy proles of K

+

and Cl

is dened relative to the positive

−

z

inside the channel. Note that the electrostatic potential direction of the simulation box. The orientation of the

protein is the MD is opposite to the convention used in experiments on LS bundles, and the sign of the voltage is opposite to that of the experiments

28,29

6


.

Page 7 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The alamethicin hexamer (Alm6) bundle was templated onto alpha-cyclodextrin (CD) at the C-terminal side through an alkyne-containing (bis-(aminoethyl)ethylene glycol) linker. Templating prevents uctuations in the number of helices forming the bundle, which is characteristic of free alamethicin

25,26

.

The templated CD-Alm6 channel was embedded in a

negatively charged 1,2-Dioleoyl-sn-glycero-3-phosphoglycerol (DOPG) membrane in a simulation unit cell of 90.5 Å

×

×

90.5 Å

107.8 Å. The system contained 217 lipid molecules,

20163 water molecules, and KCl ions at 150 mM concentration.

This setup was consis-

tent with the conditions in the experiments on the templated Alm channel carried out by Hjørringgaard

et al. 40 .

The initial structure of the Alm6 bundle was taken from a barrel-stave model

41

in which

six Alm peptides align across the DOPG membrane around a water-lled, ion-conducting pore. After initial equilibration of the system for 30 ns, a non-equilibrium MD trajectory 2400 ns in length was collected with an applied voltage of 200 mV (electric eld was pointing from the C-terminal to the N-terminal of the channel), During the simulations, CD-Alm6 underwent structural uctuations between an open, conducting state and a partly closed, poorly conducting structure. Results reported here are only for the fully open state. The CHARMM36 all-atom force eld

36

was applied to describe DOPG lipids. The CD

molecule and the alkyne linker that is covalently bound to the peptide was described with the updated CHARMM force eld for carbohydrate general force eld for small molecule was used to describe the protein.

43

42

and the newly developed CHARMM

, respectively. CHARMM27 force eld with CMAP

Other computational details were the same as in the

simulations of LS3.

The electrodiusion equation Assume that an ion moves through a transmembrane channel diusively in the presence of a constant applied voltage and the PMF due to all other components of the system, i.e.

the protein assembly, membrane, water and other ions.

7


Note that the PMF contains


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 39

contributions that are often considered separately, such as the dipole contribution from the lipid head groups, the charge layering of the ions and the electrostatics of protein. Additionally, we are only considering the movement of ions through a water-lled channel, and we are assuming that a 1-dimensional model can be used.

16

The PMF, non equilibrium

free energy prole, and other spatially dependent quantities are considered only within the pore region. Both the electric eld and the PMF are acting along the

z -axis perpendicular to

the water-membrane interface. Assume further that the PMF is independent of voltage. If ionic concentrations are constant in time, current,

J , of ions of a given type across the channel

can be described by way of the one-dimensional, steady-state Smoluchowski equation:

J = −D(z) where

D(z)

dρ(z) dE(z) + βρ(z) dz dz

is the position-dependent diusivity,

,

(1)

is ion density per unit length along

β = 1/kB T ,

ρ(z)

rather than a more customary probability distribution function, which is proportional

to

ρ(z),

is the Boltzmann constant and

T

z

and

where

kB

ρ(z)

is temperature. Here, we use

because the former quantity is directly available from simulations.

E(z)

is dened

as:

E(z) = A(z) + qV (z), where

A(z)

is the PMF at equilibrium, in the absence of the electric potential,

is ionic charge.

V (z)

V (z),

is assumed to change linearly through the channel between

and

zmin

q

and

zmax

V (z) = −Eel (z − zmin ). Thus, the total applied voltage, eld.

∆V ,

is equal to -Eel (zmax

− zmin ),

where

Eel

is the electric

Since voltage and electric eld play the roles of potential and force, respectively,

positive applied voltage corresponds to negative electric eld.

8


Page 9 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Eq. 1 is a rst-order, linear ordinary dierential equation which can always be solved formally. In practice, it is advantageous to integrate this equation, as it allows to use data collected for the whole channel. To do so, we multiply both sides of Eq. 1 by a function and carry out integration in the limits between

z1

z2 .

and

Function

g(z)

g(z)

is arbitrary, but

non-zero in the integration limits. This yields

Z

z2

J z1

g(z) dz = − D(z)

z2

Z

z1

dρ(z) g(z) dz − β dz

Z

z2

g(z)ρ(z) z1

dE(z) dz. dz

The rst integral on the right hand side (r.h.s.) can be carried out by parts

Z

z2

z1

dρ(z) g(z) dz = ρ(z2 )g(z2 ) − ρ(z1 )g(z1 ) − dz

Z

z2

dg(z) dz dz z1 Z z2 d ln(g(z)) = ρ(z2 )g(z2 ) − ρ(z1 )g(z1 ) − dz. ρ(z)g(z) dz z1 ρ(z)

Then, we obtain

J = R z2

1

g(z)dz D(z) Z z2

ρ(z1 )g(z1 ) − ρ(z2 )g(z2 )

z1

+ z1

If we choose

d ρ(z)g(z) [ln(g(z)) − βE(z)] dz . dz

g(z) = exp [βE(z)],

(2)

the derivative in the integral in the numerator vanishes

identically and the integrated equation for

J

takes the form:

ρ(z1 )eβE(z1 ) − ρ(z2 )eβE(z2 ) J= . R z2 exp(βE(z)) dz D(z) z1 This formula is often called the electrodiusion equation. Since

(3)

exp [βE(z)]

is known to be

the integrating factor for the ED equation, our approach to integrating the ED equation is slightly more complicated than required. We take this route for two reasons. First, other choices of

g(z)

are possible and, in fact, will be used further in this paper.

Second, this

formulation stresses that integrating Eq. 1 also involves reweighing the data. Considering

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 39

that estimates of the functions that enter the ED equation are burdened with errors, the currents calculated by way of dierent reweighing schemes will dier.

If this dierence is

substantial, it most likely indicates that at least one quantity has been estimated poorly. Calculating

J

from Eq. 3 requires knowledge of

A(z)

and

D(z)

between

z1

and

z2 ,

the

applied voltage and the boundary density terms. This is the only integration scheme that does not involve the full, non-equilibrium density prole. For a system in the steady state,

J

does not depend, in principle, on the choice of

z1

exists, primarily due to inaccuracies in determining

A(z)

If

and

D(z)

and

A(z)

have been already obtained, e.g.

the boundary densities

ρ1

and

ρ2

z2 .

and

In practice, this dependence

D(z).

from equilibrium simulations, only

are required to estimate the current. These densities can

be obtained from relatively short MD simulations at a desired applied voltage. There is no need to accumulate good statistics of crossing events at this voltage. In fact, the boundary terms can be, in principle, calculated without observing any crossing events at all. If

A(z)

and

D(z)

can be estimated in the full range between

simplied. The equilibrium ion density,

ρeq (z),

is related to

zmin

and

zmax ,

Eq. 3 can be

A(z)

ρeq (z) = ρref exp(−βA(z)),

where

ρref

insures correct unit conversion. Taking advantage of this relation, Eq. 3 can be

rewritten as

J=

ρ(zmax ) exp(q∆V /kB T ) ρeq (zmax ) . R zmax exp(qV (z)/kB T ) dz ρeq (z)D(z) zmin

ρ(zmin ) ρeq (zmin )

−

Assuming that the electric eld vanishes at and

ρ(zmax )/ρeq (zmax )

zmin

and

zmax , we assume that ρ(zmin )/ρeq (zmin )

are independent of applied voltages (but dependent on concentra-

tions). If concentrations on both sides of the membrane are equal, which is a common setup in computer simulations, and are the same as those used to obtain

10


ρeq

then

ρ(zmin ) = ρeq (zmin )

Page 11 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and

ρ(zmax ) = ρeq (zmax ).

This yields

1 − exp(βq∆V ) J = R zmax exp(βqV (z)) . dz zmin ρeq (z)D(z) From this formula it follows that

J

can be calculated at dierent applied voltages only from

the knowledge of the equilibrium properties,

A(z) and D(z) between zmin and zmax , as long as

the ED equation is valid. In principle, separate knowledge of It is sucient to know

(4)

A(z) and D(z) is not required.

ρeq (z)D(z).

In most traditional approaches to solving the Nernst-Planck equation, a number of additional approximations about the structure and behavior of the system at the boundaries are involved. Usually, the radius of the pore at the ion entry needs to be known and ion ux needs to be estimated

44,45

. Conventionally, this is done by integrating over a hemisphere

assuming a capture radius, although other schemes have also been used. It has been argued that the resulting ion uxes have to be corrected for electrostatic pore resistance

46,47

. The

corresponding parameters and procedures are burdened with considerable uncertainties. If application of the ED equation is coupled with MD simulations, there is no need to consider these ambiguous issues. The same quantities that are used to calculate current are also sucient to estimate ionic selectivity. The current of positive ions,

Jp ,

is given as

1 − exp(βq∆V ) Jp = R zmax exp(βqV (z)) dz zmax ρpeq (z)Dp (z) and the current of negative ions,

Jn ,

is

1 − exp(β|q|∆V ) 1 − exp(β|q|∆V ) 1 − exp(−β|q|∆V ) . Jn = R zmax exp(−β|q|V (z)) = − R zmax exp[β|q|(∆V −V (z))] = − R zmax exp(β|q|V (z)) dz dz dz ρn (z)Dn (z) ρn (z)Dn (z) ρn (zmax −z)Dn (zmax −z) 0 0 0 eq

Here, and

ρpeq (z)

ρneq (z)

and

and

eq

Dp (z)

Dn (z)

eq

are equilibrium density and diusivity of positive ions, respectively,

are the same quantities for negative ions. The charge on positive ions

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

is

q

Page 12 of 39

and the charge on negative ions is assumed for simplicity to have the same absolute

value. The origin of the coordinate system along

z

is assumed to be at

zmin = 0,

again for

simplicity. The second equality was obtained after multiplying numerator and denominator by

exp(β|q|∆V ),

whereas the third equality follows from change of variables

z → zmax − z .

Then, ionic selectivity is

|Jp | = |Jn |

R zmax 0

exp(β|q|V (z)) dz ρn eq (zmax −z)Dn (zmax −z) . R zmax exp(βqV (z)) dz p 0 ρeq (z)Dp (z)

Thus, ionic selectivity as a function of voltage can be estimated from equilibrium quantities. This is possible because we can obtain current,

Jp + Jn ,

is known.

Jp

and

Jn

separately. In experiments, only the total

As a consequence, selectivity is determined from the reversal

potential by way of the Goldman-Hodgkin-Katz (GHK) equation

29,44

.

If a concentration

dierence is established on both sides of the membrane, the reversal potential is the applied voltage at which there is no net current through the channel. Alternatively, selectivity can be determined from measurements of currents when positive or negative ions are suciently large that they cannot enter the channel.

Results Testing the distribution of crossing events At the core of the ED equation is the assumption that ion-crossing events are statistically independent.

This does not have to be the case.

For example, transport through a volt-

age gated potassium channel, KcsA, involves concerted motion of several, most likely four, ions

4853

.

Therefore, ion-crossing events cannot be considered as independent, at least as

long as a single-ion order parameter is used. In many other instances, the description of ion transport as a series of independent events, in accord with the ED equation, may be quite satisfactory. Even if there are good reasons to believe that this description is correct, it is

12


Page 13 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


still prudent to test whether this is indeed the case. If ions cross a channel independently of each other the distribution of waiting times between consecutive crossing events should follow Poisson statistics. between

n+1

If

n

waiting times

crossing events have been ordered in the ascending order according to their

length, then the length,

ti ,

of the

i-th

waiting time is

12

:

i 1 ti = − ln(1 − ), λ N where

λ

is the inverse of the average waiting time and

(5)

N

is the normalization constant. It

insures that the probability distribution of all waiting times is normalized to 1. From this equation, it follows that

λ.

ti

should be a linear function of

ln(1 − Ni ) with the slope dened by

Then, crossing events can be mapped to a stationary Poisson process, as required in the

ED equation. Furthermore, the slope yields an estimate of the current, which is expected to be very close to the estimate obtained from the direct count of crossing events. Perhaps more importantly, once it has been established that ion transport obeys Poisson statistics, it follows that statistical error in the number of crossing events, turn, provides error estimate for

n,

is equal to

√

n.

This, in

J.

Eq. 5 cannot be used directly because

N

is not known. If the total simulation time is

large, compared to the average waiting time, then we can substitute a simple, least squares method for estimating

λ.

n

for

N.

This yields

A number of other methods exist for this

purpose, including those based on maximum likelihood and moments

54

.

As can be seen from Fig. 1, the assumption that ion transport is a stationary Poisson process holds for all three channels examined here. Deviations are observed only in the long waitingtime tail of the distribution, where statistics is poor. If ion transport follows Poisson statistics, the estimate of the current from

λ

might be more precise than the estimate from

the count of crossing events, as it is less dependent on statistics for long waiting times.

13



0

ln(1-i/N)

0 -3 -6 -9 -12

-2

LS3: K+

20

0

60

40

ln(1-i/N)

200 mV 100 mV -100 mV

-6

0

40

20

60

0

-2

-2

TTX: K+

-4 -6

LS3: Cl-

-4

200 mV 100 mV -100 mV

0

0

5

TTX: Cl-

-4

100 mV 50 mV -50 mV

10

15

-6

20

100 mV 50 mV -50 mV

0

5

0

ln(1-i/N)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 39

10

15

20

t i (ns)

-2 ALM: K+

-4

200 mV

0

10

20

t i (ns)

30

40

Figure 1: Plot of ln(1 − i/N ) vs. ti to test if crossing events follow Poisson statistics for + − K ions (left) and Cl ions (right) for the LS3 channel (top), trichotoxin (TTX) channel (middle) and CD-Alm6 (ALM) channel (bottom).

Determining potential of mean force A(z)

can be determined from equilibrium simulations in the absence of applied voltage in

which ionic concentrations on both sides of the membrane are set to be equal. A number of methods, such as umbrella sampling and have been applied Alternatively,

12,16,17,20,31

55,56

and the adaptive biasing force

18

A(z) can be obtained with the aid of the ED equation from non-equilibrium g(z)

in Eq. 2 is chosen to be equal to

integration in the numerator can be carried out analytically. Then,

J,

, are available

for this purpose.

simulations in the presence of electric eld. If

in terms of

5759

A(z)

1/ρ(z),

can be expressed

applied eld, diusivity and non-equilibrium free energy prole.

Z

z

A(z) − A(z1 ) = kB T ln ρ(z1 ) − ln ρ(z) − J z1

1 0 dz + qEel (z − z1 ). ρ(z 0 )D(z 0 )

14


(6)

Page 15 of 39

A(z) (kcal/mol)

A(z) (kcal/mol)

Plots of

A(z)

reconstructed from this equation are shown in Fig. 2. In accord with the as-

3 2 1

-20

1.6 1.2 0.8 0.4 0 4 3 2 1 0

2

10

20

0

-20

3

TTX K+

-5

0

0

-10

10

20

TTX Cl-

2

0V 50 mV -50 mV 100 mV

-10

0V 200 mV 100 mV -100 mV

1 0

-10

LS3 Cl-

3

0V 200 mV 100 mV -100 mV

0 -1

4

LS3 K+

-15

A(z) (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0V 50 mV -50 mV 100 mV

1

5

10

15

0 -15

-10

-5

0

5

10

15

z (Å)

ALM K+ ABF 200 mV

-20

-10

0

10

20

z (Å) Figure 2: Free energy proles for K

+

ions (left) and Cl

−

ions (right) for the LS3 channel

(top), trichotoxin (TTX) channel (middle) and CD-Alm6 (ALM) channel (bottom).

The

PMF was obtained from equilibrium density proles at 0 V (LS3 and TTX) and by way of ABF for CD-Alm6. The ED equation was used to compute the PMFs in the non-equilibrium simulations at several applied voltages.

sumptions underlying the ED equation, the reconstructed

A(z) proles are, in general, inde-

pendent of voltage and remains in a very good agreement with the PMFs determined in equilibrium simulations. In most cases,

A(z)

calculated from equilibrium and non-equilibrium

simulations dier less than statistical errors of individual proles. This result is somewhat unexpected. Since equilibrium simulations carried out to calculate

A(z)

are usually strati-

ed, which means that an ion spends a considerable time in each strata (window) along

z,

even in those corresponding to low-probability regions, one would expect these approaches to yield more accurate PMF than non-equilibrium simulations in which congurations near the top of the free energy barrier are sampled only during occasional crossing events. Good per-

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

formance of the procedure in which

A(z)

is reconstructed from non-equilibrium simulations

by way of Eq. 6 may be in part due to the fact that ion crossings are independent events and, thus, congurations of the system corresponding to these events are largely uncorrelated. It should be, however, kept in mind that application of Eq. 6 may be less successful for poorly conducting channels for which ion crossing events are quite rare. Some dierences between tions appear for K

+

A(z)

determined in equilibrium and non-equilibrium simula-

in CD-Alm6 ( A(z) for Cl

−

was not considered because no chloride

current was observed). The PMF obtained from the equilibrium simulations, but not from the reconstruction, exhibits a small dip between 4 and 12 Å. Another dierence was found in the PMF for Cl

−

in LS3 at higher electric elds. This might be due to subtle changes in

the shape of the pore or reorientation of some residues caused by the applied voltage. Alternatively, large voltage ramp might amplify small inaccuracies in other quantities that enter the ED equation. Finally, it is also possible that some assumptions regarding the underlying stochastic process used to derive the ED equation start to break down at high electric elds. At any rate, carrying out simulations at high voltages to improve the statistics for crossing events

12,16

may create diculties in extrapolating currents to markedly lower, physiological

voltages within the framework of the ED equation.

Applied voltage In the ED equation, it has been assumed that applied voltage changes linearly across the membrane between of

zmin

and

zmax

zmin and zmax .

This assumption has not been so far tested and the values

have not been determined.

It is expected that the electric eld applied

across the simulation box transfers to the membrane, as water or an aqueous solution of electrolytes are conducting phases. The membrane bilayer, however, is not a homogeneous phase, but instead consists of an interfacial headgroup regions and a nonpolar, hydrocarbon core. Although headgroups are polar, they are not mobile and, therefore, it is not clear to what extent they can be treated as conducting.

To shed light on this issue we return to

16


Page 16 of 39

Page 17 of 39

A(z),

the equation is rearranged such that it is

V (z), assuming that A(z) is known e.g.

from equilibrium simulations. This yields:

Eq. 6. This time, instead of solving it for solved for

Z qV (z) − qV (z1 ) = kB T ln ρ(z1 ) − ln ρ(z) − J

z

z1

1 0 dz − A(z) − A(z1 ). ρ(z 0 )D(z 0 )

(7)

As can be seen from Fig. 3, the voltage across the POPC membrane containing the

|V(z)| (mV)

200

|V(z)| (mV)

200

LS3 K+ 200 mV 100 mV -100 mV

100

100

0

LS3 Cl200 mV 100 mV -100 mV

0

-20

-10

10

0

20

TTX K+ 100 mV 50 mV -50 mV

100 50

-20 100 50

0

-10

0

10

20

0

10

20

TTX Cl100 mV 50 mV -50 mV

0

-20

|V(z)| (mV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-10

0

10

20

-20

-10

z (Å)

ALM K+ 200 mV

100 50 0 -20

-10

0

z (Å)

10

20

Figure 3: Voltage ramps computed from subtraction of the PMF from the non-equilibrium + − free energy prole for K ions (left) and Cl ions (right) for the LS3 channel (top), trichotoxin (TTX) channel (middle) and CD-Alm6 (ALM) channel (bottom). Values for the ramps are shows as colored dots for each applied voltage.

trichotoxin and LS3 channels, reconstructed from Eq. 7, is noisy but exhibits a nearly linear trend with

z

approximately in the range

−15 < z < 15

Å, and the total voltage change

across the membrane is quite close to the applied voltage. Determining the precise shape of

V (z)

would require markedly longer simulations. The limits of the voltage ramp coincide

with the maxima in the densities of the carbon atom adjacent to the phosphorus atom on the hydrocarbon tail side of POPC. This indicates that both water and the headgroup region

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 39

form a conducting phase and the electric eld aects mostly the hydrocarbon core. For the CD-Alm6 system, the reconstructed

V (z) inside the membrane is not as linear as for the other

two systems. This is mainly due to the dip of approximately 0.5 kcal/mol in from equilibrium simulations, which is located near section. This dip is absent in

A(z)

z

A(z)

obtained

= 510 Å, as discussed in the previous

reconstructed from non-equilibrium simulations. Also,

the CD-Alm6 channel is markedly more exible than the other two channels studied here. Fluctuations in the pore may have negative impact on statistical precision in determining

V (z).

Even though the channel was embedded in a charged DOPG rather than POPC

membrane, the limits of the voltage ramp estimated from the slope of

V (z)

are also located

between the phosphorus and the adjacent carbon atom. It is possible that these limits are universal to phospholipid bilayers.

Non-equilibrium density proles Non-equilibrium density proles obtained from simulation of trichotoxin and LS3 channels at several voltages are shown in Fig. 4. For both channels, the proles at dierent voltages, including

∆V = 0,

converge to the same values at

zmin

zmax .

and

This justies the validity

of Eq. 4. For a more detailed discussion on calculating density proles, especially near the ends of the channel see the Supplementary Material. If

A(z)

and

D(z)

are known, the density prole for a given type of ions at a specied

applied voltage can be calculated if

J

is known. This is clear after rearranging Eq. 3

"

Z

z

ρ(z) = exp(−βE(z)) ρ(z1 ) exp(βE(z1 )) − J z1

If the ED equation holds,

ρ(z)

0

#

exp(βE(z )) 0 dz . D(z 0 )

(8)

obtained from Eq. 8 and calculated directly from MD simu-

lations should match. As seen in Fig. S2, this match is not always perfect, largely because the reconstructed proles are quite sensitive to errors in the free energy prole (see Fig. S2).

18


Page 19 of 39

LS3 K+ 0 mV 200 mV 100 mV -100 mV

−1

ρ(z) (Å )

0.4 0.3

0.1

0.1

0

0

0.3

−1

0.3 0.2

-10

10

0

20

-20 0.3

TTX K+ 100 mV 50 mV -50 mV

0.2

0.1

0

0 -10

-5

0

5

10

15

0V 200 mV 100 mV -100 mV

-10

0

-15

-10

-5

z (Å) Figure 4: Ion density proles of K

20

10

TTX Cl- 100 mV 50 mV -50 mV

0.2

0.1

-15

LS3 Cl-

0.4

0.2

-20

ρ(z) (Å )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

5

10

15

z (Å) +

(left) and Cl

−

(right) across the LS3 channel (upper)

and trichotoxin (TTX) channel (lower panels) at several voltages.

Alternatively, one can substitute the ED equation, Eq. 4, for

J

into Eq. 8, which yields:

I(zmin , z) I(zmin , z) ρ(z) = exp[−βqV (z)] 1 − + exp(βq∆V ) , ρeq (z) I(zmin , zmax ) I(zmin , zmax ) where the integral

I(z1 , z)

(9)

is dened as:

Z

z

I(z1 , z) =

dz 0

z1

exp[βqV (z 0 )] . ρeq (z 0 )

This expression for the density prole is positive denite. This is in contrast the Eq. 8 where errors in

0.

J

and in the free energy prole can lead to negative values of density when

J >

In addition, Eq. 9 allows for the reconstruction of the density proles from equilibrium

properties alone and guarantees correct behavior of these proles at the ends of the voltage ramp. This, however, does not mean that the overall agreement with the proles obtained

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

from MD simulations is better.

Page 20 of 39

The results of this reconstruction is shown in Fig. S2.

desired, it might be possible to improve the reconstruction of of combining forward (starting with

z1 )

ρ(z),

and backward (starting with

if

for example by way

z2 )

calculations with

weights chosen to reduce overall error.

Diusivity There are a number of ways to estimate diusivity in ion channels. Some of them require a series of additional, usually relatively short simulations, whereas others rely on extracting

D(z)

from simulations of ion transport with and/or without applied voltage. One approach

D

that is independent of the ED equation relies on the Einstein relation. It applies when is not a rapidly changing function of

z.

If the position-dependent average mean force is

subtracted from the total force acting on the ion then the limiting distribution of ions along

z

in the channel should be uniform. The diusivity at position

z0

can be obtained from a

series of MD trajectories initiated from statistically independent congurations drawn from an equilibrium ensemble in which the ion is at as long as

D

z0 .

Then,

hz 2 i after time t is equal to 2D(z0 )t,

remains approximately constant in the range of

z

traveled by the ion in time

t.

The advantage of this approach is that it allows for testing whether ions diuse according to the Fick law, as assumed in the ED equation, or in single le, as is the case in the gramicidin A channel

60

. In both types of diusion, the distribution of ionic positions at xed time is

expected to be Gaussian, but the functional dependence of the width of this distribution on

t is dierent 17 . D

If the distribution does not remain Gaussian with time it might indicate that

changes as ions move away from their initial position or that

A(z)

near

z0

is inaccurate.

Ions that traverse the antiamoebin ion channel have been shown to follow Fickian diusion

12

.

Ions in the trichotoxin, LS3 and CD-Alm6 channels exhibit similar behavior. This is shown in Fig. S3 in the Supplementary material. For trichotoxin, the diusivities of both K

+

and CL

−

were calculated from the Einstein

relation in the range [-8, 8] Å with respect to the center of the pore

20


16

.

No statistical

Page 21 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


support was found for a hypothesis that diusivities are position-dependent in this range. Similar results were obtained for the other two channels. were calculated at positions and 0.40

×

10

−5

cm

2

s

−1

z

for

For LS3, the ionic diusivities

= 0, -4, and 4 Å. They were, respectively, equal to 0.42, 0.31

K+

and 0.25, 0.29 and 0.27

corresponding average values were 0.37 and 0.27

×

10

−5

×

cm

2

10 s

−5

−1

cm

2

for K

s

+

−1

− for Cl .

− and Cl .

The

These

values are nearly an order of magnitude smaller than the corresponding diusivities in TIP3P water. They are also smaller than diusivities in the trichotoxin channel by approximately a factor of three, consistent with the fact that the transmembrane pore in trichotoxin is over 1 Å wider than the pore in LS3. For the CD-Alm6, the diusivities were 0.18 ( ±0.03), 0.14 (±0.04), and 0.19 ( ±0.04)

×

10

−5

2 −1 cm s , for z = -6 Å, 0 Å, and 5 Å, respectively. This

means that diusivity in the central region of all three channels is approximately constant. This, however, is not expected to be a feature of all channels. Another method for calculating diusivity is based on the uctuation-dissipation theorem. Short MD trajectories are generated with an ion at a xed position along the channel and the diusivity at this position is extracted from the force-force autocorrelation function

61

. Although this method has been used to calculate

D(z) of ions in channels 12 , achieving

high reliability with this approach is a challenge. A number of other methods have been proposed to calculate position-dependent diusivity

17,62,63

.

These methods are based on dierent approximations or models and some

of them are not expected to perform well for ion channels. simultaneously

A(z)

and

D(z)

to the data from MD simulations

one could estimate the running integral knowledge of

A(z), J

and

An elegant approach is to t

ρ(z),

6466

.

In the same spirit,

Rz

1 dz by way of Eq. 10 (below) from the zmin ρeq (z)D(z)

and then extract

D(z).

(See Fig. S4 in the Supplementary

Material.) These approaches, however, share the same problem: the tted

D(z) captures all

inaccuracies hidden in other quantities and/or in models used for estimating diusivity. As a result, it is impossible to determine independently its reliability. For example, to compensate for inaccurate

ρ(z)

at the large-z end of the channel,

D(z)

21


from a t in this region would


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 39

be smaller rather than larger than diusivity in the middle of the channel, which is clearly unphysical.

(See Fig. S4 in the Supplementary Material.)

attempted to extract

D(z)

For this reason, we have not

from non-equilibrium simulations in this study. This, however,

does not mean that such a procedure is generally impossible. Taking advantage of the mean value theorem for integrals, Eq. 3 can be rewritten as:

J = D(zmv )

In general,

zmv

ρ(z1 )eβE(z1 ) − ρ(z2 )eβE(z2 ) R z2 . exp(βE(z))dz z1

is a function of applied voltage. The integrand in the denominator of Eq. 3

is the largest when

E(z) reaches the maximum and decreases exponentially with E(z).

This

means that the value of the integral is largely determined by the integrand in the region of the barrier due to the PMF and applied voltage. If diusivity is approximately constant in this region it can be substituted for

D(zmv )

without appreciable loss of accuracy. The value

of diusivity in the regions that provide only small contribution to the integral has only small eect on the calculated current. In other words, the assumption of constant

D

is not

equivalent to the assumption that diusivity is constant everywhere in the channel, but only in the region of large

E(z).

To test the validity of this approximation, we solved Eq. 6 for

LS3 at voltages -200, -100 and 100 mV self-consistently for

D

to obtain the best t with

for

D

of 0.35 and 0.27

×

10

A(z)

and position-independent

A(z) determined from equilibrium simulations. −5

The best values

2 −1 + − cm s for K and Cl , respectively, are quite close to the

diusivities determined from the Einstein relation. This excellent agreement justies using constant

D

for the channels considered in this study. This approximation, however, should

be used with care. As shown in Fig. S4, in the Supplementary material, the relevant integral increases rapidly in a relatively narrow region of

z

that is only approximately 10-15 Å wide,

but this region shifts depending on the voltage. Thus, there might not be a single choice of

D

that is proper for a wide range of positive and negative voltages.

22


Page 23 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ionic currents As it has been already mentioned, currents calculated by way of Eq. 3 for a steady state

z1

system should be independent of the limits of integration

and

z2 .

In practice, this is

not expected to be the case because of inaccuracies in determining the functions that enter the integral

A(z), D(z), ρ(z1 )

and

ρ(z2 ).

Deviations from constant

J

due to the change of

integration boundaries provide a measure of reliability of the calculated current that includes both statistical and systematic errors. From the practical viewpoint, it is also of interest to establish whether there are preferred integration ranges that yield the best accuracy. At the ends of all channels studied here, ion densities are the highest and, therefore, the associated relative statistical errors are the smallest. This might suggest that integrating in these regions should yield relatively accurate currents. On the other hand, these densities are very similar at dierent voltages, oering very little discriminating power. The opposite is true in the center of the channel. Relative statistical errors increase but so does discriminating power. A simple integration scheme in which the lower bound in Eq. 3 is set to bound is systematically increased produce stable values of

J

zmin

whereas the upper

through LS3 for both ions and at

dierent voltages. This is shown in Fig. 5. A single exception is K

+

at 200 mV. This might be

due to relatively small inaccuracies in the calculated PMF that are amplied at large applied voltages. The eect of integration boundaries on

J

can be systematically studied by way of

changing both the lower and the upper integration limit. As shown in Fig. 6, integrating dierent regions along

z

yields dierent relative errors. One trend apparent from this gure

is that integrating over large ranges of

z

is preferred to integrating over small intervals. On

the other hand, integrating over the full length of the channel is not required to improve accuracy. This is valuable, because obtaining reliable

A(z)

near the ends of the channel is

often challenging due to fraying ends, lipid penetration and diculties in identifying ions that can be considered inside the channel rather than non-specically interacting with the headgroup region of the membrane. Relative errors in the currents obtained from the ED equation, compared to MD simu-

23



60

|J(z - zmin )| (pA)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

K+ 200 mV 100 mV -100 mV

Page 24 of 39

Cl- 200 mV 100 mV -100 mV

40

20

0 0

10

20

30

z - zmin (Å)

Figure 5: Magnitudes of the currents in the ED models of the LS3 channel computed as a function of the range of integration across the channel, from

lations, are of the order of 15-50%.

zmin

= -16 to

z

= -12 to 16 Å.

These errors should be kept in perspective.

It is not

unusual that dierent electrophysiological measurements of conductance for the same channel produce dierences of the same order or larger

40,67,68

. A comparison between currents

calculated from MD simulations and the ED equation is given in Table 1.

Table 1: Ionic currents for MD and the electrodiusion model IK (pA) LS3

V (mV)

MD

ED

± 1.8 ± 0.7 4.3 ± 0.5 -77 ± 5 -35 ± 4 29 ± 4

MD

ED

± 1.2 ± 0.5 3.2 ± 0.4 -35 ± 3 -14 ± 3 14 ± 3

200

-29.8

-43

-13.8

-26

100

-7.8

-13

-4.6

-7

-100 TTX

ICl (pA)

100 50 -50

4

± 13 ± 11 26 ± 8

-86 -43

24


4

± 12 ±6 -16 ± 6

41

16

Page 25 of 39

Z

10

K+ 200 mV

Cl- 200mV

K+ 100 mV

Cl- 100mV

K+ -100 mV

Cl- -100 mV

0

-10

Z

10 0

-10 10 Z

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

1

-10 20 10 Integration Range

0

30

0

20 10 Integration Range

in the LS3 channel at 200 mV, 100 mV and -100 mV. The

zmax − zmin

and the

x-axis

0

30

Figure 6: Relative errors in ED currents with respect to the MD values for K range,

2

+

and Cl

−

ions

is the length of integration

y -axis, Z , is (zmin + zzmax )/2 dening the the center of integration

for each window.

Ohmic behavior and rectication If current through the channel is proportional to applied voltage, the channel is Ohmic. If current is asymmetric with respect to the reversal of the electric eld, the channel is called rectifying. Both, Ohmic and rectifying behavior depend primarily on the shape of the pore and the composition of residues forming the pore. This raises a question: when is a channel Ohmic? Some insight into this problem can be gained from integrating Eq. 2 with

g(z) = exp[βA(z)].

J = R z2

1

exp[βA(z)] dz D(z) Z z2

ρ(z1 ) exp[βA(z1 )] − ρ(z2 ) exp[βA(z2 )]

z1

d + ρ(z) exp[βA(z)] [βA(z) − βA(z) − βqV (z)] dz dz z1 Z z2 1 ρ(z1 ) ρ(z2 ) ρ(z) = R z2 − + βqEel dz . 1 dz ρeq (z1 ) ρeq (z2 ) z1 ρeq (z) z1 ρeq (z)D(z)

25


(10)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

If the limits of integration are set to the edges of the voltage ramp, densities at these values of

z

Page 26 of 39

zmin and zmax , and the

are assumed to be independent of voltage then Eq. 10 simplies

J = R zmax zmin

Z

βqEel 1 dz ρeq (z)D(z)

zmax

zmin

ρ(z) dz. ρeq (z)

Nonlinearities in the dependence of current on electric eld are captured in the last integral in this equation. If it is equal to a constant independent of voltage, the channel is Ohmic. From the requirement that to

zmax − zmin

which

ρ(z) → ρeq (z)

for Ohmic channels.

ρeq (z)D(z)

as

E → 0

it follows that the integral is equal

As can be readily veried from Eq. 4, channels for

is constant will exhibit Ohmic behavior. As shown in Appendix A in the

Supplementary Material, this is not only a sucient but also a necessary condition for Ohmic behavior. If

D is the same everywhere in the channel then ρeq (z) and, therefore, A(z) are also

constant. This assumption is commonly used in several basic equations of electrophysiology, such as the Nernst-Planck and GHK equations

44

.

We further explore possible shapes of I-V curves. As long as the ED equation applies and

ρeq (z)D(z)

is not constant, these curves lie either above or below the straight line

characteristic of the Ohmic behavior for all voltages. No other shapes are admitted by the equation. Which of these behaviors will be observed can be immediately gauged from the shape of the PMF. Specically, assuming that

z,

the current for positive

q∆V

D(z)

is not a quickly changing function of

will be always higher than expected from Ohmic behavior

if the equilibrium ion density at the high-voltage end of the channel is greater than the average density in the channel. The opposite is true if the density is lower. This is derived in Appendix B in the Supplementary Material. A channel is not rectifying if and only if the free energy prole is symmetric with respect to the midpoint of the voltage ramp. A proof is included in Appendix C in the Supplementary Material.

Thus, once

A(z)

is available, simple visual inspection of the prole is sucient

to assess whether the channel is rectifying. How the degree of rectication is related to the

26


Page 27 of 39

asymmetry of

A(z),

is shown in the proof.

In general, in contrast to what is frequently

assumed, both Ohmic behavior in a broad range of electric elds and the lack of rectication are associated with special properties of the free energy proles rather than with common features of ion channels. The I-V curves for K

+

in LS3 predicted from the ED equation, shown in Fig. 7. illustrate

10 0 -10

I (pA)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-20

-30

Ohmic A0 (z) A1 (z) -40 MD

-200

-100

0 V (mV)

Figure 7: Current-Voltage curves for K

+

100

200

ions in the LS3 channel calculated using the ED

equation with the equilibrium free energy prole,

A0 (z),

and the free energy prole,

A1 (z),

reconstructed from the simulations at the applied voltage of 100 mV. The currents calculated from MD simulations are shown for comparison. The Ohmic current is for the equilibrium free energy prole.

the dependencies discussed above. As expected from the shape of non-Ohmic and rectifying.

A(z),

the channel is both

Since the ion density markedly increases toward the ends of

the channel (see Fig. 4), the current increases with increasing positive voltage faster than predicted from the strictly Ohmic behavior. The actual I-V curve for LS3 is quite similar for both positive and negative voltages

28,29

. Although experimentally only the total current

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 39

is measured, it can be directly compared with the current of K selectivity for this ion. As expected from the shape of

A(z),

+

because of the channel

the same trends are observed

for the trichotoxin channel, although they are weaker.

Summary and conclusions The ED equation is a very useful tool to calculate the current-voltage dependence and ion selectivity from a single equilibrium or non-equilibrium simulation. For simple channels, it yields results in good agreement with currents and selectivities calculated from MD simulations. Deviations are typically of a similar magnitude as uncertainties in experimental measurements of the same properties. Since the system is in a steady state, the current is in principle independent of the limits of integration in Eq. 3. This is convenient because many channels are frayed or partially disordered at the ends and, therefore, the PMF determined in these regions is uncertain. If the limits of integration are set to the edges of the voltage ramp then, under mild assumptions, knowledge of equilibrium properties,

A(z)

and

D(z),

is

sucient to calculate the current. If these limits are set closer to the center of the channel the boundary density terms also need to be known. They can be determined from MD simulations at a desired applied voltage that are markedly shorter than simulations required to estimate the current. Accumulating statistics of crossing events is not required. In general, it appears that integrating over large ranges of

z

for which

A(z)

can be reliably determined

yields the best estimates of currents. The Smoluchowski equation can be integrated with a dierent integrating factor to solve the ED equation for

A(z).

quantities, in this case

J

This is an example of "inverse problem", whereby calculated

and

ρ(z),

are used to reconstruct causal functions that produced

them. In statistical mechanics this procedure is often called "inverse Boltzmann". For ion channels, it means that an equilibrium property, simulations. In all cases studied here,

A(z)

A(z),

is recovered from non-equilibrium

obtained from equilibrium and non-equilibrium

28


Page 29 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MD agree with each other very well.

In principle, a similar procedure can be used to

reconstruct the other equilibrium quantity that appears in the ED equation -

D(z).

There

are, however, concerns that inaccuracies in the quantities needed for reconstructing positiondependent diusivity might signicantly impact the results. For this reason, we turned to a method independent of the ED equation that is based on the Einstein relation. Precise knowledge of

D(z)

in the full integration range is usually not required. Without signicant

loss of precision, it is sucient to have good estimates of diusivity in the range near the barrier due to the PMF and applied voltage. Another application of the inverse procedure is to examine how applied voltage changes in the system. It was found that the electric eld is non-zero only in the non-polar core of the phospholipid bilayer. As expected, it appears to change linearly along the bilayer normal, although the reconstructed voltage dependence is fairly noisy. A number of basic assumptions underlying the ED equation have been systematically tested and appear to hold. Waiting times between consecutive ion crossings follow Poisson statistics, indicating that ion crossings are independent events.

This is very convenient

because estimates of statistical errors for the calculated currents immediately follow from this fact. Diusion of ions exhibits Fickian behavior. Good agreement between

A(z) obtained at

dierent applied voltages indicates that the PMF is independent of electric eld, as assumed in the ED equation.

Also, electric eld appears to be constant across the bilayer.

This

conrms the validity of the ED equation, at least for simple channels. The results suggest that dierences between the currents obtained from MD simulations and the ED equation are primarily due to inaccuracies in determining

A(z) rather than to approximations involved

in the ED model. The utility of the ED equation arises largely from the fact that it provides a theoretically justied means for extrapolating simulation results obtained at one applied voltage (including zero voltage) to other voltages.

For example, poorly conducting channels would require

extensive simulations to collect sucient statistics of crossing events.

29


Instead, one can


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 39

carry out equilibrium simulations or non-equilibrium simulations at a larger electric eld and recover the current at the voltage of interest. In the latter case, care should be taken not to step outside the range of applicability of the ED equation. For example, entry to a channel might become diusion-limited

69

. Also, the ED equation provides an ecient route

to calculating I-V curves over a broad range of voltages, a task that would be computationally very intensive if carried out via direct MD simulations.

Finally, from this equation, one

can draw a number of general conclusions about electrophysiological behavior of dierent channels. For example, as has been shown here, constant

ρeq (z)D(z)

across a channel is not

only a sucient, but also the necessary condition for Ohmic behavior. is not rectifying if and only if

ρeq (z)D(z)

Also, the channel

is symmetric with respect to the midpoint of the

voltage ramp. Simple, visual inspection of the calculated

A(z) is sucient to assess a number

of electrophysiological characteristics of a given channel, such as the direction of deviations from the Ohmic behavior and rectication. The conclusions of this work apply directly to simple channels. To what extent the ED equation also holds for large, complex, natural channels or synthetic structures, such as nanotubes, has not been yet determined, even though this issue is of considerable interest. One obvious diculty emerges when individual ion crossings cannot be considered as independent events. Another complication arises if ion transport is gated, as diusion model is not a natural framework for describing gating processes. An interesting and potentially fruitful line of research would be to examine the ED equation in collective variables or combine them with both MD and transition state theories suitable for ion channels.

Acknowledgements This work was supported by the NASA Exobiology and Astrobiology Programs. All simulations were performed at the NASA Advanced Supercomputing (NAS) Division and on the Anton computer at the Pittsburgh Supercomputer Center.

30


Page 31 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Aksimentiev, A.; Schulten, K. Imaging

α-hemolysin

with molecular dynamics: Ionic

conductance, osmotic permeability, and electrostatic potential map.

Biophys. J. 2005,

88, 37453761. (2) Chung, S. H.; Allen, T. W.; Kuyucak, S. Modeling diverse range of potassium channels with Brownian dynamics.

Biophys. J. 2002, 83, 263277.

(3) Noskov, S. Y.; Im, W.; Roux, B. Ion permeation through the

α-hemolysin channel:

The-

oretical studies based on Brownian dynamics and Poisson-Nernst-Plank electrodiusion theory.

Biophys. J. 2004, 87, 22992309.

(4) Cheng, M. H.; Coalson, R. D. An accurate and ecient empirical approach for calculating the dielectric self-energy and ion-ion pair potential in continuum models of biological ion channels.

J. Phys. Chem. B 2005, 109, 488498.

(5) Coalson, R. D.; Kurnikova, M. G. Poisson-Nernst-Planck theory approach to the calculation of current through biological ion channels.

IEEE Trans. Nanobiosci. 2005, 4,

8193.

(6) Dyrka, W.; Bartuzel, M. M.; Kotulska, M. Optimization of 3D Poisson-Nernst-Planck model for fast evaluation of diverse protein channels.

Proteins 2013, 81, 18021822.

(7) Crozier, P. S.; Rowley, R. L.; Holladay, N. B.; Henderson, D.; Busath, D. D. Molecular dynamics simulation of continuous current ow through a model biological membrane channel.

Phys. Rev. Lett. 2001, 86, 2467.

(8) Sachs, J. N.; Crozier, P. S.; Woolf, T. B. Atomistic simulations of biologically realistic transmembrane potential gradients.

J. Chem. Phys. 2004, 121, 1084710851.

(9) Biró, I.; Pezeshki, S.; Weingart, H.; Winterhalter, M.; Kleinekathofer, U. Comparing

31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 39

the temperature-dependent conductance of the two structurally similar E. coli porins OmpC and OmpF.

Biophys. J. 2010, 98, 18301839.

(10) Pezeshki, S.; Chimerel, C.; Bessonov, A. N.; Winterhalter, M.; Kleinekathöfer, U. Understanding ion conductance on a molecular level: An all-atom modeling of the bacterial porin OmpF.

Biophys. J. 2009, 97, 18981906.

(11) Faraudo, J.; Calero, C.; Aguilella-Arzo, M. Ionic partition and transport in multi-ionic channels: a molecular dynamics simulation study of the OmpF bacterial porin.

Biophys.

J. 2010, 99, 21072115. (12) Wilson, M. A.; Wei, C.; Bjelkmar, P.; Wallace, B. A.; Pohorille, A. Molecular dynamics simulation of the antiamoebin ion channel: Linking structure and conductance.

Biophys.

J. 2011, 100, 23942402. (13) Kutzner, C.; Grubmuller, H.; de Groot, B. L.; Zachariae, U. Computational electrophysiology: The molecular dynamics of ion channel permeation and selectivity in atomistic detail.

Biophys. J. 2011, 101, 809817.

(14) Chandler, D. E.; Penin, F.; Schulten, K.; Chipot, C. The p7 protein of hepatitis C virus forms structurally plastic, minimalist ion channels.

PLoS Comput. Biol. 2012, 8,

e1002702.

(15) Ulmschneider, M. B.; Bagneris, C.; McCusker, E. C.; Decaen, P. G.; Delling, M.; Clapham, D. E.; Ulmschneider, J. P.; Wallace, B. A. Molecular dynamics of ion transport through the open conformation of a bacterial voltage-gated sodium channel.

Proc.

Natl. Acad. Sci. U.S.A. 2013, 110, 63646369. (16) Wilson, M. A.; Nguyen, T. H.; Pohorille, A. Combining molecular dynamics and an electrodiusion model to calculate ion channel conductance.

141, 22D519. 32


J. Chem. Phys. 2014,

Page 33 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(17) Roux, B.; Allen, T.; Bernèche, S.; Im, W. Theoretical and computational models of biological ion channels.

(18) Pohorille, A. In

Q. Rev. Biophys. 2004, 37, 15103.

Computational Biophysics of Membrane Proteins ;

Domene, C., Ed.;

Royal Society of Chemistry: Oxford, 2016; in press.

(19) Allen, T. W.; Andersen, O. S.; Roux, B. Energetics of ion conduction through the gramicidin channel.

Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 117122.

(20) Zhu, F.; Hummer, G. Pore opening and closing of a pentameric ligand-gated ion channel.

Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 1981419819. (21) Adelman, J. L.; Grabe, M. Simulating current-voltage relationships for a narrow ion channel ising the weighted ensemble method.

J. Chem. Theory Comput. 2015, 11,

19071918.

(22) Zhu, F.; Hummer, G. Theory and simulation of ion conduction in the pentameric GLIC channel.

J. Chem. Theory Comput. 2012, 8, 37593768.

(23) Chugh, J. K.; Brückner, H.; Wallace, B. A. Model for a helical bundle channel based on the high resolution crystal structure of trichotoxin_A50E.

Biochemistry 2002, 41,

1293412941.

(24) Duclohier, H.; Alder, G. M.; Bashford, C. L.; Bruckner, H.; Chugh, J. K.; Wallace, B. A. Conductance studies on trichotoxin_A50E and implications for channel structure.

Bio-

phys. J. 2004, 87, 17051710. (25) Hall, J.; Vodyanoy, I.; Balasubramanian, T.; Marshall, G. Alamethicin. A rich model for channel behavior.

Biophys. J. 1984, 45, 233247.

(26) Laver, D. R. The barrel-stave model as applied to alamethicin and its analogs reevaluated.

Biophys. J. 1994, 66, 355359.

33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 39

(27) Duclohier, H. Helical kink and channel behaviour: A comparative study with the peptaibols alamethicin, trichotoxin and antiamoebin.

Eur. Biophys. J. 2004, 33, 169174.

(28) Lear, J. D.; Wasserman, Z. R.; DeGrado, W. F. Synthetic amphiphilic peptide models for protein ion channels.

Science 1988, 240, 11771181.

(29) Lear, J. D.; Schneider, J. P.; Kienker, P. K.; DeGrado, W. F. Electrostatic eects on ion selectivity and rectication in designed ion channel peptides.

J. Am. Chem. Soc.

1997, 119, 32123217. (30) Roux, B. Computational electrophysiology:

the molecular dynamics of ion channel

permeation and selectivity in atomistic detail.

Biophys. J. 2011, 101, 755756.

(31) Khalili-Araghi, F.; Ziervogel, B.; Gumbart, J. C.; Roux, B. Molecular dynamics simulations of membrane proteins under asymmetric ionic concentrations.

J. Gen. Physiol.

2013, 142, 465475. (32) Shaw, D. E.; Denero, M. M.; Dror, R. O.; Kuskin, J. S.; Larson, R. H.; Salmon, J. K.; Young, C.; Batson, B.; Bowers, K. J.; Chao, J. C. et al. Anton, a special-purpose machine for molecular dynamics simulation.

Commun. ACM 2008, 51, 9197.

(33) Zhong, Q.; Jiang, Q.; Moore, P. B.; Newns, D. M.; Klein, M. L. Molecular dynamics simulation of a synthetic ion channel.

Biophys. J. 1998, 74, 310.

(34) Randa, H. S.; Forrest, L. R.; Voth, G. A.; Sansom, M. S. P. Molecular dynamics of synthetic leucine-serine ion channels in a phospholipid membrane.

Biophys. J. 1999,

77, 24002410. (35) Grigoryan, G.; Degrado, W. F. Probing designability via a generalized model of helical bundle geometry.

(36) Klauda, J. B.;

J. Mol. Biol. 2011, 405, 10791100.

Venable, R. M.;

Freites, J. A.;

O'Connor, J. W.;

Tobias, D. J.;

Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of

34


Page 35 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the CHARMM all-atom additive force eld for lipids: Validation on six lipid types.

J.

Phys. Chem. B 2010, 114, 78307843. (37) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water.

J. Chem. Phys.

1983, 79, 926935. (38) MacKerell, A. D., Jr.; Brooks, B.; Brooks, C. L., III; Nilsson, L.; Roux, B.; Won, Y.; Karplus, M. In Allinger, N. L.,

The Encyclopedia of Computational Chemistry ; Clark, T.,

Gasteiger, J.,

Kollman, P. A.,

Schleyer, P. v. R.,

Schaefer, H. F., III,

Schreiner, P. R., Eds.; John Wiley and Sons: New York, 1998; Vol. 1; pp 271277.

(39) MacKerell, J., A.D.; Feig, M.; Brooks, I., C.L. Extending the treatment of backbone energetics in protein force elds: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations.

J.

Comp. Chem. 2004, 25, 14001415. (40) Hjørringgaard, C. U.; Vad, B. S.; Matchkov, V. V.; Nielsen, S. B.; Vosegaard, T.; Nielsen, N. C.; Otzen, D. E.; Skrydstrup, T. Cyclodextrin-scaolded alamethicin with remarkably ecient membrane permeabilizing properties and membrane current conductance.

J. Phys. Chem. B 2012, 116, 76527659.

(41) Duclohier, H.; Wroblewski, H. Voltage-dependent pore formation and antimicrobial activity by alamethicin and analogues.

J Membr. Biol. 2001, 184, 112.

(42) Guvench, O.; Mallajosyula, S. S.; Raman, E. P.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T. J.; Jamison, F. W.; Mackerell, A. D. CHARMM additive all-atom force eld for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling.


(43) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I. et al. CHARMM general force eld: A

35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 39

force eld for drug-like molecules compatible with the CHARMM all-atom additive biological force elds.

(44) Hille, B.

J. Comput. Chem. 2010, 31, 671690.

Ion Channels of Excitable Membranes (3rd edition) ;

Sinauer Associates Inc.:

Sunderland, MA, 2001.

(45) Levitt, D. G. Interpretation of biological ion channel ux datareaction-rate versus continuum theory.

Annu. Rev. Biophys. Biophys. Chem. 1986, 15, 2957.

(46) Hall, J. E. Access resistance of a small circular pore.

J. Gen. Physiol. 1975, 66, 531532.

(47) Aguilella-Arzo, M.; Aguilella, V. M.; Eisenberg, R. Computing numerically the access resistance of a pore.

Eur. Biophys. J. 2005, 34, 314322.

(48) Morais-Cabral, J. H.; Zhou, Y.; MacKinnon, R. Energetic optimization of ion conduction rate by the K+ selectivity lter.

Nature 2001, 414, 3742.

(49) Roux, B. Ion conduction and selectivity in K

+

channels.

Annu. Rev. Biophys. Biomol.

Struct. 2005, 34, 153171. (50) Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Dynamics of K+ ion conduction through Kv1.2.

Biophys. J. 2006, 91, L7274.

(51) Furini, S.; Domene, C. Atypical mechanism of conduction in potassium channels.

Proc.

Natl. Acad. Sci. U.S.A. 2009, 106, 1607416077. (52) Egwolf, B.; Luo, Y.; Walters, D. E.; Roux, B. Ion selectivity of alpha-hemolysin with beta-cyclodextrin adapter. II. Multi-ion eects studied with grand canonical Monte Carlo/Brownian dynamics simulations.

J Phys Chem B 2010, 114, 29012909.

(53) Köpfer, D. A.; Song, C.; Gruene, T.; Sheldrick, G. M.; Zachariae, U.; de Groot, B. L. Ion permeation in K+ channels occurs by direct Coulomb knock-on. 352355.

36


Science 2014, 346,

Page 37 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(54) Dudewicz, E.; Mishra, S.

Modern Mathematical Statistics ;

Wiley series in probability

and mathematical statistics; Wiley, 1988.

(55) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free energy estimation: Umbrella sampling.

J. Comput. Phys. 1977, 23, 187199.

(56) Shell, M. S.; Panagiotopoulosa, A.; Pohorille, A. In

and Applications to Chemistry and Biology ;

Free Energy Calculations. Theory

Chipot, C., Pohorille, A., Eds.; Springer-

Verlag: Berlin, 2007; pp 79122.

(57) Darve, E.; Pohorille, A. Calculating free energies using average force.

J. Chem. Phys.

2001, 115, 91699183. (58) Darve, E.; Rodriguez-Gomez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations.

(59) Chipot, C.; Pohorille, A.

J. Chem. Phys. 2008, 128, 35633578.

Free Energy Calculations. Theory and Applications to Chem-

istry and Biology ; Springer-Verlag:

Berlin, 2007.

(60) Levitt, D. G. Dynamics of a single-le pore: Non-Fickian behavior.

Phys. Rev. A 1973,

8, 30503054. (61) Zwanzig, R. W. Elementary derivation of time-correlation formulas for transport coefcients.

J. Chem. Phys. 1964, 40, 25272533.

(62) Berne, B. J.; Borkovec, M.; Straub, J. E. Classical and modern methods in reaction rate theory.

J. Phys. 1988, 92, 37113725.

(63) Woolf, T. B.;

Roux, B. Conformational exibility of o-phosphorylcholine and o-

phosphorylethanolamine:

A molecular dynamics study of solvation eects.

Chem. Soc. 1994, 116, 59165926.

37


J. Am.


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 39

(64) Hummer, G. Position-dependent diusion coecients and free energies from Bayesian analysis of equilibrium and replica molecular dynamics simulations.

New J. Phys. 2005,

7, 34. (65) Peter, C.; Hummer, G. Ion transport through membrane-spanning nanopores studied by molecular dynamics simulations and continuum electrostatics calculations.

Biophys.

J. 2005, 89, 22222234. (66) Comer, J.; Chipot, C.; González-Nilo, F. D. Calculating position-dependent diusivity in biased molecular dynamics simulations.


(67) You, S.; Peng, S.; Lien, L.; Breed, J.; Sansom, M. S.; Woolley, G. A. Engineering stabilized ion channels: covalent dimers of alamethicin.

Biochemistry 1996, 35, 6225

6232.

(68) Asami, K.; Okazaki, T.; Nagai, Y.; Nagaoka, Y. Modications of alamethicin ion channels by substitution of Glu-7 for Gln-7.

Biophys. J. 2002, 83, 219228.

(69) Andersen, O. S. Ion movement through gramicidin A channels. Studies on the diusioncontrolled association step.

Biophys. J. 1983, 41, 147.

38


Page 39 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(left) picture of equilibrated LS3 hexameric bundle in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine membrane and (rigt) current-voltage curve from MD (blue duts), predicted from electrodiffusion model with different free energy profiles (red and green) and ohmic response (dashed line). 40x19mm (300 x 300 DPI)


Validity of the Electrodiffusion Model for Calculating Conductance of

Recommend Documents