Validity of the Electrodiffusion Model for Calculating Conductance of

Article pubs.acs.org/JPCB

Validity of the Electrodiffusion 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, Moffett Field, California 94035, United States Department of Pharmaceutical Chemistry University of California, San Francisco, California 94132, United States § SETI Institute, 189 N Bernardo Ave #200, Mountain View, California 94043, United States ‡

S Supporting Information *

ABSTRACT: We examine the validity and utility of the electrodiffusion (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 diffusivity, obtained from molecular dynamics simulations, the ED equation can be used to determine the full current−voltage dependence with modest or no additional effort. The potential of mean force can be obtained not only from equilibrium simulations, but also, with comparable accuracy, from nonequilibrium 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 sufficient conditions for Ohmic and nonrectifying 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 dynamics2−6 into the realm of all-atom simulations in which system dynamics and nonelectrostatic 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 presence of electric field7,8 and were followed by a number of studies in which ion current through different channels was calculated directly from MD simulations.9−16 The connection between different 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 different channels in the absence of applied potential (for a recent review see Pohorille18). The PMF is the main ingredient of many molecular theories of ion conductance. In particular, it enters the generalized Nernst−Planck equation for the flux of ions in which the gradient of the electrostatic potential is replaced by the average total force experienced by © 2016 American Chemical Society

an ion. This equation is often called the generalized electrodiffusion (ED) equation or, simply, the electrodiffusion equation. From another perspective, the ED equation can be considered a simplified 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,19−22 An important advantage of the ED equation is its simplicity. The only ingredient that is difficult 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 diffusion 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 Special Issue: Klaus Schulten Memorial Issue Received: September 22, 2016 Revised: November 22, 2016 Published: November 23, 2016 3607

DOI: 10.1021/acs.jpcb.6b09598 J. Phys. Chem. B 2017, 121, 3607−3619

Article

The Journal of Physical Chemistry B

diffusivities and the comparison between currents calculated from MD and the ED model. Finally we 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 Supporting Information. The article concludes with a summary.

force acting on a single ion in one dimension. Most of the assumptions underlying the ED equation have not been systematically tested. In this article, we undertake this task with regard to the one-dimensional ED equation, using as models three simple channels. Two of them, are built of natural, nongenomic peptaibols from fungi: trichotoxin23,24 and alamethicin.25−27 The third one is a synthetic channel, LS3, built of two types of amino acids, leucine and serine.28,29 All three channels are known to mediate transport of K+ and Cl−. 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 efficiently and accurately in combination with MD simulations. In this article, we do not address the comparisons between the calculated and experimental structures and electrophysiology of these channels. This has been16 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 fields, are quite different 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 effort. It provides a means to extrapolate currents obtained at one applied voltage to other voltages, a task that is otherwise difficult 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 efficiently 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 differ 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 nonequilibrium ion density profile.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 nonequilibrium 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 profiles, discussion of methods for calculating

■

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 Figure S1 in the Supporting Information. 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.28,29 In an αhelical conformation, the peptide has a hydrophobic, leucinerich 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 ion-conducting pore.28,29,33,34 The initial channel model was built using parameters of an ideal Crick coiled coil bundle.35 The bundle was placed in a 1palmitoyl-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 1 M KCl concentration. The size of the simulation box was 85.283 Å × 85.154 Å × 107.670 Å with periodic boundary conditions applied in all three spatial direction. Long ranged forces were calculated by way of the particle-mesh Ewald method using 64 × 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 field was used to describe the POPC phospholipids,36 and water molecules were represented by the TIP3P model.37 The LS3 peptide was described using CHARMM22 with CMAP corrections.38,39 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 different 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, 2753, and 960 ns. The ion density profiles obtained in the absence of voltage were used to calculate the free energy profiles of K+ and Cl− inside the channel. Note that the electrostatic potential is defined relative to the positive z direction of the simulation box. The orientation of the protein in 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 3608


Article


where A(z) is the PMF at equilibrium, in the absence of the electric potential, V(z), and q is ionic charge. V(z) is assumed to change linearly through the channel between zmin and zmax

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 fluctuations 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 consistent 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 model41 in which six Alm peptides align across the DOPG membrane around a water-filled, ion-conducting pore. After initial equilibration of the system for 30 ns, a nonequilibrium MD trajectory 2400 ns in length was collected with an applied voltage of 200 mV (electric field was pointing from the C-terminal to the N-terminal of the channel), During the simulations, CD-Alm6 underwent structural fluctuations 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 field36 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 field for carbohydrate42 and the newly developed CHARMM general force field for small molecule,43 respectively. CHARMM27 force field with CMAP was used to describe the protein. Other computational details were the same as in the simulations of LS3. Electrodiffusion Equation. Assume that an ion moves through a transmembrane channel diffusively 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. Note that the PMF contains 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-filled channel, and we are assuming that a 1dimensional model can be used.16 The PMF, nonequilibrium free energy profile, and other spatially dependent quantities are considered only within the pore region. Both the electric field 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: ⎛ d ρ (z ) dE(z) ⎞ ⎟ J = −D(z)⎜ + βρ(z) ⎝ dz dz ⎠

V (z) = −,el(z − zmin)

Thus, the total applied voltage, ΔV, is equal to −,el(zmax − zmin), where ,el is the electric field. Since voltage and electric field play the roles of potential and force, respectively, positive applied voltage corresponds to negative electric field. Eq 1 is a first-order, linear ordinary differential 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 g(z) and carry out integration in the limits between z1 and z2. Function g(z) is arbitrary, but nonzero in the integration limits. This yields J

∫z

z2

1

g (z) dz = − D(z)

∫z

z2

1

g (z)

dρ(z) dz − β dz

∫z

z2

g (z)ρ(z)

1

dE(z) dz dz

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

∫z

z2

g (z)

1

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

∫z

= ρ(z 2)g (z 2) − ρ(z1)g (z1) −

∫z

z2

ρ(z)

1

z2

dg (z) dz dz

ρ(z)g (z)

1

d ln(g (z)) dz dz

Then, we obtain ⎧ ⎪ J = z g(z)dz ⎨ρ(z1)g (z1) − ρ(z 2)g (z 2) + 2 ∫z D(z) ⎪ ⎩ 1 1

∫z

z2

1

⎫ ⎪ d ρ(z)g (z) [ln(g (z)) − βE(z)]dz ⎬ dz ⎪ ⎭

(2)

If we choose g(z) = exp[βE(z)], the derivative in the integral in the numerator vanishes identically and the integrated equation for J takes the form: J=

ρ(z1)e βE(z1) − ρ(z 2)e βE(z 2)

∫z

z 2 exp(βE(z)) dz D(z)

1

(3)

This formula is often called the electrodiffusion equation. Since 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 that estimates of the functions that enter the ED equation are burdened with errors, the currents calculated by way of different reweighing schemes will differ. If this difference 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, nonequilibrium density profile. For a system in

(1)

where D(z) is the position-dependent diffusivity, ρ(z) is ion density per unit length along z and β = 1/kBT, where kB is the Boltzmann constant and T is temperature. Here, we use ρ(z) rather than a more customary probability distribution function, which is proportional to ρ(z), because the former quantity is directly available from simulations. E(z) is defined as E(z) = A(z) + qV (z) 3609


Article

The Journal of Physical Chemistry B the steady state, J does not depend, in principle, on the choice of z1 and z2. In practice, this dependence exists, primarily due to inaccuracies in determining A(z) and D(z). If A(z) and D(z) have been already obtained, e.g., from equilibrium simulations, only the boundary densities ρ1 and ρ2 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 zmin and zmax, eq 3 can be simplified. The equilibrium ion density, ρeq(z), is related to A(z)

Jn =

J=

∫z

|Jp |

ρ(z max ) exp(qΔV /kBT ) ρeq (z max )

|Jn |

z max exp(qV (z) / kBT ) dz ρeq (z)D(z)

z max exp(βqV (z)) dz ρeq (z)D(z) min

1 − exp(β|q|ΔV ) z max

exp(β | q | V (z)) dz ρeqn (z max − z)Dn(z max − z)

=

∫0

z max

exp(β | q | V (z)) dz ρeqn (z max − z)Dn(z max − z)

∫0

z max exp(βqV (z)) dz ρeqp (z)Dp(z)

■

(4)

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 voltage gated potassium channel, KcsA, involves concerted motion of several, most likely four, ions.48−53 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 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. If n waiting times between n + 1 crossing events have been ordered in the ascending order according to their length, then the length, ti, of the i-th waiting time is12

From this formula it follows that J can be calculated at different 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 A(z) and D(z) is not required. It is sufficient to know ρ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 flux 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 fluxes 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 sufficient to estimate ionic selectivity. The current of positive ions, Jp, is given as Jp =

∫

Thus, ionic selectivity as a function of voltage can be estimated from equilibrium quantities. This is possible because we can obtain Jp and Jn separately. In experiments, only the total current, Jp + Jn, is known. As a consequence, selectivity is determined from the reversal potential by way of the Goldman−Hodgkin−Katz (GHK) equation.29,44 If a concentration difference 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 sufficiently large that they cannot enter the channel.

min

1 − exp(βqΔV )

∫z

z max exp[β |q| (ΔV − V (z))] dz ρeqn (z)Dn(z) 0

Here, ρpeq(z) and Dp(z) are equilibrium density and diffusivity of positive ions, respectively, and ρneq(z) and Dn(z) are the same quantities for negative ions. The charge on positive ions is q 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

Assuming that the electric field vanishes at zmin and zmax, we assume that ρ(zmin)/ ρeq(zmin) and ρ(zmax)/ ρeq(zmax) are independent of applied voltages (but dependent on concentrations). 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 ρeq then ρ(zmin) = ρeq(zmin) and ρ(zmax) = ρeq(zmax). This yields J=

1 − exp(β|q|ΔV )

=−

∫

∫0

where ρref insures correct unit conversion. Taking advantage of this relation, eq 3 can be rewritten as −

z max exp(−β | q | V (z)) dz ρeqn (z)Dn(z) 0

=−

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

ρ(z min) ρeq (z min)

1 − exp( −β|q|ΔV )

1 ⎛ i ⎞ ti = − ln⎜1 − ⎟ N⎠ λ ⎝

(5)

where λ is the inverse of the average waiting time and 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

1 − exp(βqΔV ) z max exp(βqV (z)) dz p z max ρeq (z)Dp(z)

∫

(

ln 1 −

and the current of negative ions, Jn, is 3610

i N

) with the slope is defined by λ. Then, crossing events DOI: 10.1021/acs.jpcb.6b09598 J. Phys. Chem. B 2017, 121, 3607−3619

Article


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).

Figure 2. Free energy profiles 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 profiles 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 nonequilibrium simulations at several applied voltages.

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, n, is equal to √n. This, in turn, provides error estimate for 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 n for N. This yields a simple, leastsquares method for estimating λ. A number of other methods exist for this purpose, including those based on maximum likelihood and moments.54

As can be seen from Figure 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 waiting−time tail of the distribution, where statistics are 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. 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 sampling55,56 and the adaptive biasing force,57−59 are available18 and have been applied12,16,17,20,31 for this purpose. 3611


Article


Figure 3. Voltage ramps computed from subtraction of the PMF from the nonequilibrium free energy profile 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 shown as colored dots for each applied voltage.

Alternatively, A(z) can be obtained with the aid of the ED equation from nonequilibrium simulations in the presence of electric field. If g(z) in eq 2 is chosen to be equal to 1/ρ(z), integration in the numerator can be carried out analytically. Then, A(z) can be expressed in terms of J, applied field, diffusivity, and nonequilibrium free energy profile. ⎡ A(z) − A(z1) = kBT ⎢ln ρ(z1) − ln ρ(z) − J ⎣ + q ,el(z − z1)

∫z

z

1

dip between 4 and 12 Å. Another difference was found in the PMF for Cl− in LS3 at higher electric fields. 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 fields. At any rate, carrying out simulations at high voltages to improve the statistics for crossing events12,16 may create difficulties 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 zmin and zmax. This assumption has not been so far tested and the values of zmin and zmax have not been determined. It is expected that the electric field 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 eq 6. This time, instead of solving it for A(z), the equation is rearranged such that it is solved for V(z), assuming that A(z) is known, e.g., from equilibrium simulations. This yields

⎤ 1 dz′⎥ ρ(z′)D(z′) ⎦

(6)

Plots of A(z) reconstructed from this equation are shown in Figure 2. In accord with the assumptions underlying the ED equation, the reconstructed A(z) profiles are, in general, independent of voltage and remain in a very good agreement with the PMFs determined in equilibrium simulations. In most cases, A(z) calculated from equilibrium and nonequilibrium simulations differ less than statistical errors of individual profiles. This result is somewhat unexpected. Since equilibrium simulations carried out to calculate A(z) are usually stratified, which means that an ion spends a considerable time in each strata (“window”) along z, even in those corresponding to lowprobability regions, one would expect these approaches to yield more accurate PMF than nonequilibrium simulations in which configurations near the top of the free energy barrier are sampled only during occasional crossing events. Good performance of the procedure in which A(z) is reconstructed from nonequilibrium simulations by way of eq 6 may be in part due to the fact that ion crossings are independent events and, thus, configurations 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 differences between A(z) determined in equilibrium and nonequilibrium simulations appear for K+ 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

⎡ qV (z) − qV (z1) = kBT ⎢ln ρ(z1) − ln ρ(z) − J ⎣ − A(z) − A(z1)

∫z

z

1

⎤ 1 dz′⎥ ρ(z′)D(z′) ⎦

(7)

As can be seen from Figure 3, the voltage across the POPC membrane containing the 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) 3612


Article


Figure 4. Ion density profiles of K+ (left) and Cl− (right) across the LS3 channel (upper) and trichotoxin (TTX) channel (lower panels) at several voltages.

free energy profile (see Figure S2). Alternatively, one can substitute the ED equation, eq 4, for J into eq 8, which yields

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 form a conducting phase and the electric field affects 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 A(z) obtained from equilibrium simulations, which is located near z = 5−10 Å, as discussed in the previous section. This dip is absent in A(z) reconstructed from nonequilibrium simulations. Also, the CDAlm6 channel is markedly more flexible 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. Nonequilibrium Density Profiles. Nonequilibrium density profiles obtained from simulation of trichotoxin and LS3 channels at several voltages are shown in Figure 4. For both channels, the profiles at different voltages, including ΔV = 0, converge to the same values at zmin and zmax. This justifies the validity of eq 4. For a more detailed discussion on calculating density profiles, especially near the ends of the channel, see the Supporting Information. If A(z) and D(z) are known, the density profile for a given type of ions at a specified applied voltage can be calculated if J is known. This is clear after rearranging eq 3 ⎡ ρ(z) = exp(− βE(z))⎢ρ(z1)exp(βE(z1)) − J ⎣

∫z

z

1

⎡ I(zmin , z) ρ (z ) = exp[−βqV (z)]⎢1 − ρeq (z) I(zmin , zmax ) ⎣ + exp(βqΔV )

I(zmin , z) ⎤ ⎥ I(zmin , zmax ) ⎦

(9)

where the integral I(z1,z) is defined as I(z1 , z) =

∫z

z

1

dz′

exp[βqV (z′)] ρeq (z′)

This expression for the density profile is positive definite. This is in contrast to eq 8 where errors in J and in the free energy profile can lead to negative values of density when J > 0. In addition, eq 9 allows for the reconstruction of the density profiles from equilibrium properties alone and guarantees correct behavior of these profiles at the ends of the voltage ramp. This, however, does not mean that the overall agreement with the profiles obtained from MD simulations is better. The results of this reconstruction is shown in Figure S2. If desired, it might be possible to improve the reconstruction of ρ(z), for example by way of combining forward (starting with z1) and backward (starting with z2) calculations with weights chosen to reduce overall error. Diffusivity. There are a number of ways to estimate diffusivity 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 that is independent of the ED equation relies on the Einstein relation. It applies when D 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 diffusivity at position z0 can be obtained from a series of MD trajectories initiated from statistically independent configurations drawn

exp(βE(z′)) ⎤ dz′⎥ ⎦ D(z′)

(8)

If the ED equation holds, ρ(z) obtained from eq 8 and calculated directly from MD simulations should match. As seen in Figure S2, this match is not always perfect, largely because the reconstructed profiles are quite sensitive to errors in the 3613


Article

The Journal of Physical Chemistry B from an equilibrium ensemble in which the ion is at z0. Then, ⟨z2⟩ after time t is equal to 2D(z0)t, as long as D 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 diffuse according to the Fick law, as assumed in the ED equation, or in single file, as is the case in the gramicidin A channel.60 In both types of diffusion, the distribution of ionic positions at fixed time is expected to be Gaussian, but the functional dependence of the width of this distribution on t is different.17 If the distribution does not remain Gaussian with time it might indicate that D 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 diffusion.12 Ions in the trichotoxin, LS3, and CDAlm6 channels exhibit similar behavior. This is shown in Figure S3 in the Supporting Information. For trichotoxin, the diffusivities of both K+ and CL− were calculated from the Einstein relation in the range [−8, 8] Å with respect to the center of the pore.16 No statistical support was found for a hypothesis that diffusivities are positiondependent in this range. Similar results were obtained for the other two channels. For LS3, the ionic diffusivities were calculated at positions z = 0, −4, and 4 Å. They were, respectively, equal to 0.42, 0.31, and 0.40 × 10−5 cm2 s−1 for K+ and 0.25, 0.29, and 0.27 × 10−5 cm2 s−1 for Cl−. The corresponding average values were 0.37 and 0.27 × 10−5 cm2 s−1 for K+ and Cl−. These values are nearly an order of magnitude smaller than the corresponding diffusivities in TIP3P water. They are also smaller than diffusivities in the trichotoxin channel by approximately a factor of 3, consistent with the fact that the transmembrane pore in trichotoxin is over 1 Å wider than the pore in LS3. For the CD-Alm6, the diffusivities were 0.18 (±0.03), 0.14 (±0.04), and 0.19 (±0.04) × 10−5 cm2s−1, for z = −6, 0, and 5 Å, respectively. This means that diffusivity 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 diffusivity is based on the fluctuation−dissipation theorem. Short MD trajectories are generated with an ion at a fixed position along the channel and the diffusivity 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 diffusivity.17,62,63 These methods are based on different approximations or models and some of them are not expected to perform well for ion channels. An elegant approach is to fit simultaneously A(z) and D(z) to the data from MD simulations.64−66 In the same spirit, one could z 1 estimate the running integral ∫ ρ (z)D(z) dz by way of eq 10 z min

Information). For this reason, we have not attempted to extract D(z) from nonequilibrium 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)

ρ(z1)e βE(z1) − ρ(z 2)e βE(z 2) z

∫z 2 exp(βE(z))dz 1

In general, zmv 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 diffusivity is approximately constant in this region it can be substituted for D(zmv) without appreciable loss of accuracy. The value of diffusivity in the regions that provide only small contribution to the integral has only small effect on the calculated current. In other words, the assumption of constant D is not equivalent to the assumption that diffusivity 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 selfconsistently for A(z) and position-independent D to obtain the best fit with A(z) determined from equilibrium simulations. The best values for D of 0.35 and 0.27 × 10−5 cm2s−1 for K+ and Cl−, respectively, are quite close to the diffusivities determined from the Einstein relation. This excellent agreement justifies using constant D for the channels considered in this study. This approximation, however, should be used with care. As shown in Figure S4, in the Supporting Information, the relevant integral increases rapidly in a relatively narrow region of z that is only 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. Ionic Currents. As it has been already mentioned, currents calculated by way of eq 3 for a steady state system should be independent of the limits of integration z1 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 different voltages, offering 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 zmin whereas the upper bound is systematically increased produce stable values of J through LS3 for both ions and at different voltages. This is shown in Figure 5. A single exception is K+ at 200 mV. This might be due to relatively small inaccuracies in the calculated PMF that are amplified at large applied voltages. The effect of integration boundaries on J can be systematically studied by way of changing both the lower and the upper integration limit. As shown in Figure 6,

eq

(below) from the knowledge of A(z), J, and ρ(z), and then extract D(z) (see Figure S4 in the Supporting Information). These approaches, however, share the same problem: the fitted D(z) captures all inaccuracies hidden in other quantities and/or in models used for estimating diffusivity. 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) from a fit in this region would be smaller rather than larger than diffusivity in the middle of the channel, which is clearly unphysical (see Figure S4 in the Supporting 3614


Article


the electric field, 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 = z exp[βA(z)] ⎨ρ(z1)exp[βA(z1)] − ρ(z 2)exp[βA(z 2)] 2 dz ⎪ ∫z D(z) ⎩ 1 1

+

z2

∫z

ρ(z)exp[βA(z)]

1

=

∫z

1

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 zmin = −16 to z = −12 to 16 Å.

⎡

⎢ ρ(z1) 1 dz ⎢⎣ ρeq (z1) ρeq (z)D(z) 1

z2

⎫ d [βA(z) − βA(z) − βqV (z)]dz ⎬ ⎭ dz −

ρ(z 2) + βq ,el ρeq (z 2)

∫z

z2

1

⎤ ρ(z) ⎥ dz ρeq (z) ⎥⎦

(10)

If the limits of integration are set to the edges of the voltage ramp, zmin and zmax, and the densities at these values of z are assumed to be independent of voltage then eq 10 simplifies

integrating different regions along z yields different relative errors. One trend apparent from this figure 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 difficulties in identifying ions that can be considered inside the channel rather than nonspecifically interacting with the headgroup region of the membrane. Relative errors in the currents obtained from the ED equation, compared to MD simulations, are of the order of 15− 50%. These errors should be kept in perspective. It is not unusual that different electrophysiological measurements of conductance for the same channel produce differences 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. Ohmic Behavior and Rectification. If current through the channel is proportional to applied voltage, the channel is Ohmic. If current is asymmetric with respect to the reversal of

J=

βq ,el

∫z

z max min

1 ρeq (z)D(z)

∫ dz z

z max min

ρ (z ) dz ρeq (z)

Nonlinearities in the dependence of current on electric field 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 ρ(z) → ρeq(z) as , → 0 it follows that the integral is equal to zmax − zmin for Ohmic channels. As can be readily verified from eq 4, channels for which ρeq(z)D(z) is constant will exhibit Ohmic behavior. As shown in Appendix A in the Supporting Information, this is not only a sufficient 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

Figure 6. Relative errors in ED currents with respect to the MD values for K+ and Cl− ions in the LS3 channel at 200, 100, and −100 mV. The x-axis is the length of integration range, zmax − zmin and the y-axis, Z, is (zmin + zmax)/2 defining the center of integration for each window. 3615


Article

The Journal of Physical Chemistry B Table 1. Ionic Currents for MD and Electrodiffusion Model IK (pA) V (mV) LS3

TTX

200 100 −100 100 50 −50

ICl (pA)

MD −29.8 −7.8 4.3 −77 −35 29

± ± ± ± ± ±

ED −43 −13 4 −86 ± 13 −43 ± 11 26 ± 8

1.8 0.7 0.5 5 4 4

MD −13.8 −4.6 3.2 −35 −14 14

± ± ± ± ± ±

ED 1.2 0.5 0.4 3 3 3

−26 −7 4 41 ± 12 16 ± 6 −16 ± 6

the shape of A(z), the same trends are observed for the trichotoxin channel, although they are weaker.

admitted by the equation. Which of these behaviors will be observed can be immediately gauged from the shape of the PMF. Specifically, assuming that D(z) is not a quickly changing function of z, the current for positive qΔV 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 Supporting Information. A channel is not rectifying if and only if the free energy profile is symmetric with respect to the midpoint of the voltage ramp. A proof is included in Appendix C in the Supporting Information. Thus, once A(z) is available, simple visual inspection of the profile is sufficient to assess whether the channel is rectifying. How the degree of rectification is related to the 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 fields and the lack of rectification are associated with special properties of the free energy profiles rather than with common features of ion channels. The I−V curves for K+ in LS3 predicted from the ED equation, shown in Figure 7, illustrate the dependencies

■

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 nonequilibrium 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 sufficient 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 different integrating factor to solve the ED equation for A(z). This is an example of “inverse problem”, whereby calculated quantities, in this case J 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, A(z), is recovered from nonequilibrium simulations. In all cases studied here, A(z) obtained from equilibrium and nonequilibrium 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 position-dependent diffusivity might significantly 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 significant loss of precision, it is sufficient to have good estimates of diffusivity 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 field is nonzero only in the nonpolar core of the phospholipid bilayer. As expected, it appears to change linearly along the bilayer normal, although the reconstructed voltage dependence is fairly noisy.

Figure 7. Current−Voltage curves for K+ ions in the LS3 channel calculated using the ED equation with the equilibrium free energy profile, A0(z), and the free energy profile, 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 profile.

discussed above. As expected from the shape of A(z), the channel is both non-Ohmic and rectifying. Since the ion density markedly increases toward the ends of the channel (see Figure 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 is measured, it can be directly compared with the current of K+ because of the channel selectivity for this ion. As expected from 3616



■

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. Diffusion of ions exhibits Fickian behavior. Good agreement between A(z) obtained at different applied voltages indicates that the PMF is independent of electric field, as assumed in the ED equation. Also, electric field appears to be constant across the bilayer. This confirms the validity of the ED equation, at least for simple channels. The results suggest that differences 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 justified 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 sufficient statistics of crossing events. Instead, one can carry out equilibrium simulations or nonequilibrium simulations at a larger electric field 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 diffusion-limited.69 Also, the ED equation provides an efficient 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 different channels. For example, as has been shown here, constant ρeq(z)D(z) across a channel is not only a sufficient, but also the necessary condition for Ohmic behavior. Also, the channel is not rectifying if and only if ρeq(z)D(z) is symmetric with respect to the midpoint of the voltage ramp. Simple, visual inspection of the calculated A(z) is sufficient to assess a number of electrophysiological characteristics of a given channel, such as the direction of deviations from the Ohmic behavior and rectification. 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 difficulty emerges when individual ion crossings cannot be considered as independent events. Another complication arises if ion transport is gated, as diffusion 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.

■

Article

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. *E-mail: [email protected]; Phone: 650-604-5759; Fax: 650-604-1088. ORCID

Andrew Pohorille: 0000-0001-8562-6595 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS 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.

■

REFERENCES

(1) Aksimentiev, A.; Schulten, K. Imaging α-hemolysin with molecular dynamics: Ionic conductance, osmotic permeability, and electrostatic potential map. Biophys. J. 2005, 88, 3745−3761. (2) Chung, S. H.; Allen, T. W.; Kuyucak, S. Modeling diverse range of potassium channels with Brownian dynamics. Biophys. J. 2002, 83, 263−277. (3) Noskov, S. Y.; Im, W.; Roux, B. Ion permeation through the αhemolysin channel: Theoretical studies based on Brownian dynamics and Poisson-Nernst-Plank electrodiffusion theory. Biophys. J. 2004, 87, 2299−2309. (4) Cheng, M. H.; Coalson, R. D. An accurate and efficient 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, 488−498. (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, 81−93. (6) Dyrka, W.; Bartuzel, M. M.; Kotulska, M. Optimization of 3D Poisson-Nernst-Planck model for fast evaluation of diverse protein channels. Proteins: Struct., Funct., Genet. 2013, 81, 1802−1822. (7) Crozier, P. S.; Rowley, R. L.; Holladay, N. B.; Henderson, D.; Busath, D. D. Molecular dynamics simulation of continuous current flow 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, 10847−10851. (9) Biró, I.; Pezeshki, S.; Weingart, H.; Winterhalter, M.; Kleinekathofer, U. Comparing the temperature-dependent conductance of the two structurally similar E. coli porins OmpC and OmpF. Biophys. J. 2010, 98, 1830−1839. (10) Pezeshki, S.; Chimerel, C.; Bessonov, A. N.; Winterhalter, M.; Kleinekathöfer, U. Understanding ion conductance on a molecular level: An all-atom modeling of the bacterial porin OmpF. Biophys. J. 2009, 97, 1898−1906. (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, 2107−2115. (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, 2394−2402. (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, 809−817. (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.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b09598. Pictures of the membrane systems, considerations of density profiles obtained from MD, reconstruction of density profiles, calculation of diffusivities, current integral range, when a channel is Ohmic, shapes of I− V curves, and rectification (PDF) 3617


Article


Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: Validation on six lipid types. J. Phys. Chem. B 2010, 114, 7830− 7843. (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, 926−935. (38) MacKerell, A. D., Jr; Brooks, B.; Brooks, C. L., III; Nilsson, L.; Roux, B.; Won, Y.; Karplus, M. In The Encyclopedia of Computational Chemistry; Schleyer, P. v. R., Allinger, N. L., Clark, T., Gasteiger, J., Kollman, P. A., Schaefer, H. F., III, Schreiner, P. R., Eds.; John Wiley and Sons: New York, 1998; Vol. 1; pp 271−277. (39) MacKerell, A. D.; Feig, M.; Brooks, C. L. Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25, 1400−1415. (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-scaffolded alamethicin with remarkably efficient membrane permeabilizing properties and membrane current conductance. J. Phys. Chem. B 2012, 116, 7652−7659. (41) Duclohier, H.; Wroblewski, H. Voltage-dependent pore formation and antimicrobial activity by alamethicin and analogues. J. Membr. Biol. 2001, 184, 1−12. (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 field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling. J. Chem. Theory Comput. 2011, 7, 3162−3180. (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 field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2009, 31, 671−690. (44) Hille, B. Ion Channels of Excitable Membranes, 3rd ed.; Sinauer Associates Inc.: Sunderland, MA, 2001. (45) Levitt, D. G. Interpretation of biological ion channel flux datareaction-rate versus continuum theory. Annu. Rev. Biophys. Biophys. Chem. 1986, 15, 29−57. (46) Hall, J. E. Access resistance of a small circular pore. J. Gen. Physiol. 1975, 66, 531−532. (47) Aguilella-Arzo, M.; Aguilella, V. M.; Eisenberg, R. Computing numerically the access resistance of a pore. Eur. Biophys. J. 2005, 34, 314−322. (48) Morais-Cabral, J. H.; Zhou, Y.; MacKinnon, R. Energetic optimization of ion conduction rate by the K+ selectivity filter. Nature 2001, 414, 37−42. (49) Roux, B. Ion conduction and selectivity in K+ channels. Annu. Rev. Biophys. Biomol. Struct. 2005, 34, 153−171. (50) Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Dynamics of K+ ion conduction through Kv1.2. Biophys. J. 2006, 91, L72−74. (51) Furini, S.; Domene, C. Atypical mechanism of conduction in potassium channels. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 16074− 16077. (52) Egwolf, B.; Luo, Y.; Walters, D. E.; Roux, B. Ion selectivity of alpha-hemolysin with beta-cyclodextrin adapter. II. Multi-ion effects studied with grand canonical Monte Carlo/Brownian dynamics simulations. J. Phys. Chem. B 2010, 114, 2901−2909. (53) Kö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. Science 2014, 346, 352−355. (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, 187−199. (56) Shell, M. S.; Panagiotopoulosa, A.; Pohorille, A. In Free Energy Calculations. Theory and Applications to Chemistry and Biology; Chipot, C., Pohorille, A., Eds.; Springer-Verlag: Berlin, 2007; pp 79−122.

(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, 6364−6369. (16) Wilson, M. A.; Nguyen, T. H.; Pohorille, A. Combining molecular dynamics and an electrodiffusion model to calculate ion channel conductance. J. Chem. Phys. 2014, 141, 22D519. (17) Roux, B.; Allen, T.; Bernche, S.; Im, W. Theoretical and computational models of biological ion channels. Q. Rev. Biophys. 1999, 37, 15−103. (18) Pohorille, A. RSC Theoretical and Computational Chemistry Series No. 10. In Computational Biophysics of Membrane Proteins; Domene, C., Ed.; Royal Society of Chemistry: Oxford, 2016; pp 59−106. (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, 117−122. (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, 19814−19819. (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, 1907−1918. (22) Zhu, F.; Hummer, G. Theory and simulation of ion conduction in the pentameric GLIC channel. J. Chem. Theory Comput. 2012, 8, 3759−3768. (23) Chugh, J. K.; Brückner, H.; Wallace, B. A. Model for a helical bundle channel based on the high resolution crystal structure of trichotoxin_A50E. Biochemistry 2002, 41, 12934−12941. (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. Biophys. J. 2004, 87, 1705− 1710. (25) Hall, J.; Vodyanoy, I.; Balasubramanian, T.; Marshall, G. Alamethicin. A rich model for channel behavior. Biophys. J. 1984, 45, 233−247. (26) Laver, D. R. The barrel-stave model as applied to alamethicin and its analogs reevaluated. Biophys. J. 1994, 66, 355−359. (27) Duclohier, H. Helical kink and channel behaviour: A comparative study with the peptaibols alamethicin, trichotoxin and antiamoebin. Eur. Biophys. J. 2004, 33, 169−174. (28) Lear, J. D.; Wasserman, Z. R.; DeGrado, W. F. Synthetic amphiphilic peptide models for protein ion channels. Science 1988, 240, 1177−1181. (29) Lear, J. D.; Schneider, J. P.; Kienker, P. K.; DeGrado, W. F. Electrostatic effects on ion selectivity and rectification in designed ion channel peptides. J. Am. Chem. Soc. 1997, 119, 3212−3217. (30) Roux, B. Computational electrophysiology: the molecular dynamics of ion channel permeation and selectivity in atomistic detail. Biophys. J. 2011, 101, 755−756. (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, 465−475. (32) Shaw, D. E.; Deneroff, 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, 91−97. (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, 3−10. (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, 2400−2410. (35) Grigoryan, G.; Degrado, W. F. Probing designability via a generalized model of helical bundle geometry. J. Mol. Biol. 2011, 405, 1079−1100. (36) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; 3618


Article

The Journal of Physical Chemistry B (57) Darve, E.; Pohorille, A. Calculating free energies using average force. J. Chem. Phys. 2001, 115, 9169−9183. (58) Darve, E.; Rodriguez-Gomez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128, 3563−3578. (59) Chipot, C.; Pohorille, A. Free Energy Calculations. Theory and Applications to Chemistry and Biology; Springer-Verlag: Berlin, 2007. (60) Levitt, D. G. Dynamics of a single-file pore: Non-Fickian behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1973, 8, 3050−3054. (61) Zwanzig, R. W. Elementary derivation of time-correlation formulas for transport coefficients. J. Chem. Phys. 1964, 40, 2527− 2533. (62) Berne, B. J.; Borkovec, M.; Straub, J. E. Classical and modern methods in reaction rate theory. J. Phys. Chem. 1988, 92, 3711−3725. (63) Woolf, T. B.; Roux, B. Conformational flexibility of ophosphorylcholine and o-phosphorylethanolamine: A molecular dynamics study of solvation effects. J. Am. Chem. Soc. 1994, 116, 5916−5926. (64) Hummer, G. Position-dependent diffusion coefficients 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 membranespanning nanopores studied by molecular dynamics simulations and continuum electrostatics calculations. Biophys. J. 2005, 89, 2222−2234. (66) Comer, J.; Chipot, C.; González-Nilo, F. D. Calculating position-dependent diffusivity in biased molecular dynamics simulations. J. Chem. Theory Comput. 2013, 9, 876−882. (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. Modifications of alamethicin ion channels by substitution of Glu-7 for Gln-7. Biophys. J. 2002, 83, 219−228. (69) Andersen, O. S. Ion movement through gramicidin A channels. Studies on the diffusion-controlled association step. Biophys. J. 1983, 41, 147.

3619


Validity of the Electrodiffusion Model for Calculating Conductance of

Recommend Documents