Open AccessArticle

Electrophysiological Properties from Computations at a Single Voltage: Testing Theory with Stochastic Simulations

Michael A. Wilson

^1,2 and

Andrew Pohorille

^1,3,*

Exobiology Branch, MS 239-4, NASA Ames Research Center, Moffett Field, CA 94035, USA

SETI Institute, 189 Bernardo Ave, Suite 200, Mountain View, CA 94043, USA

Department of Pharmaceutical Chemistry, University of California, San Francisco, CA 94132, USA

Author to whom correspondence should be addressed.

Entropy 2021, 23(5), 571; https://doi.org/10.3390/e23050571

Submission received: 1 December 2020 / Revised: 24 April 2021 / Accepted: 28 April 2021 / Published: 6 May 2021

(This article belongs to the Special Issue Physics of Ionic Conduction in Narrow Biological and Artificial Channels)

Download

Browse Figures

Review Reports Versions Notes

Abstract

We use stochastic simulations to investigate the performance of two recently developed methods for calculating the free energy profiles of ion channels and their electrophysiological properties, such as current–voltage dependence and reversal potential, from molecular dynamics simulations at a single applied voltage. These methods require neither knowledge of the diffusivity nor simulations at multiple voltages, which greatly reduces the computational effort required to probe the electrophysiological properties of ion channels. They can be used to determine the free energy profiles from either forward or backward one-sided properties of ions in the channel, such as ion fluxes, density profiles, committor probabilities, or from their two-sided combination. By generating large sets of stochastic trajectories, which are individually designed to mimic the molecular dynamics crossing statistics of models of channels of trichotoxin, p7 from hepatitis C and a bacterial homolog of the pentameric ligand-gated ion channel, GLIC, we find that the free energy profiles obtained from stochastic simulations corresponding to molecular dynamics simulations of even a modest length are burdened with statistical errors of only 0.3 kcal/mol. Even with many crossing events, applying two-sided formulas substantially reduces statistical errors compared to one-sided formulas. With a properly chosen reference voltage, the current–voltage curves can be reproduced with good accuracy from simulations at a single voltage in a range extending for over 200 mV. If possible, the reference voltages should be chosen not simply to drive a large current in one direction, but to observe crossing events in both directions.

Keywords:

computational electrophysiology; electrodiffusion model; stochastic simulations; current–voltage dependence; reversal potential; committor probabilities

1. Introduction

Ion channels are ubiquitous in living systems in which they mediate ion transport across cell walls [1,2,3]. Although all confirmed structures of ion channels are either bundles of

α

-helices or

β

-barrels organized around a transmembrane, water-filled pore lined largely with hydrophilic side chains, they markedly differ in their properties. Their activity is regulated by a variety of signals, such as voltage, ligands, pH or mechanical tension. Some channels are made of peptides that barely span the membrane, while others are among the largest protein assemblies in a cell. In terms of ionic conductance, defined as the ratio of ionic current to voltage, channels differ by more than two orders of magnitude and conductance is not correlated with size. For example, the single-channel conductance of a bacterial homolog of pentameric ligand gated ion channels (pLGICs), GLIC, which consists of 317 residues per subunit is 8 pS [4], similar to the lowest conductance level of a channel made of antimicrobial peptide, alamethicin, which is built of 20 amino acids [5]. Another channel-forming peptide trichotoxin (TTX), consisting of 7 helices, each containing 18 residues conducts ions at 850–900 pS [6], which is close to the conductance of mechanosensitve channels MscS containing 250–1100 residues [7], approximately equal to 1 nS [8]. Some channels exhibit exquisite selectivity whereas others are non-selective. Single point mutations can not only markedly affect conductance and selectivity but even render a channel inactive or constitutively open [9,10,11]. How variations in a common, general architecture translate to a variety of electrophysiological behaviors is of great interest not only for understanding regular biological systems but also for explaining a number of diseases associated with improper function of ion channels [12,13,14].

The availability of high-resolution structural models of ion channels has created opportunities to connect structure and function. Molecular dynamics (MD) computer simulations can contribute to this goal by providing mechanistic and thermodynamic descriptions of ion transport that is not readily accessible from experimental studies [15,16,17,18,19,20,21,22,23,24,25,26]. For a recent, comprehensive review, see Flood, et al. [27]. Furthermore, MD simulations can be used to validate experimentally derived structural models, which do not always correspond to the native structures of channels [28], select the native structure among several candidates [29], and predict functional effects of mutations. These simulations, however, have to be validated by demonstrating that they can be used to reproduce measured electrophysiological properties with satisfactory reliability.

Calculating electrophysiological properties from MD simulations with applied voltage can be done simply by way of computing the current across the simulation cell [15,30,31] or counting the number of ions that cross the channel [16,24,25,32]. However, this direct method, especially when applied to obtain I-V curves and reversal potentials, requires significant computational effort, as it involves MD simulations at a number of applied voltages. To obtain the same accuracy, channels with low conductance generally require longer simulations than channels with high conductance. For example, in simulations of TTX, we counted almost 200 K

^{+}

crossing events in 900 ns at 50 mV [33], whereas only 23 Na

^{+}

crossing events were observed in a 7.7

μ

s simulation of GLIC at 100 mV (Wilson and Pohorille, unpublished). Since the ionic currents from both simulations appear to obey Poisson statistics, we expect the relative errors in the conductance of K

^{+}

in TTX and Na

^{+}

in GLIC to be approximately 7% and 20%, respectively. To achieve the same relative errors for currents in GLIC as in TTX, a MD trajectory of over 60

μ

s in length would be required. This means that calculating the I-V curve might present a considerable challenge. To deal with this challenge, it has been common to improve statistics for ion crossing events by applying high voltages, sometimes substantially above their physiological values, or increasing ionic concentration in bulk solution [15,16,17,34,35,36,37,38]. This approach, however, is fraught with dangers, as it might lead to the disruption of membrane structure or saturation effects for ion entry to a channel [37]. Furthermore, if high voltages are used, I-V curves in physiologically relevant ranges are obtained via interpolation or extrapolation procedures of unknown accuracy [35].

If the motion of ions through the channel can be satisfactorily described as diffusion in the applied electric field and the potential of mean force (PMF) exerted by all other components of the system, which is assumed to be independent of voltage, then the computational effort can be markedly reduced. Several approaches not based on MD take advantage of this description. Methods based on Poisson-Nernst-Planck (PNP) theory rely on solving the electrodiffusion (ED) equation for electrical current in which the mobile ions are represented as a mean-field concentration profile whose distribution and motion is determined by electrostatic forces [39,40,41,42,43,44]. In Brownian Dynamics (BD), channel conductance is calculated by way of solving the Langevin equation in which both short-ranged interactions with a static model of the channel and long-ranged, electrostatic interactions are taken into account [40,45,46,47,48]. In both PNP and BD approaches, electrostatic forces are are obtained from the Poisson equation and the medium is represented as a continuum. From this perspective, we take an approach in which atomistic, dynamic information offered by MD is combined with the efficiency of the ED equation. Instead of carrying out a series of extensive MD calculations of a channel over a range of applied voltages, a substantially reduced set of simulations is combined with the one-dimensional ED model in a steady state. In this approach, the actual electrophysiological properties, such as the current, are calculated from the ED equation, whereas the quantities needed to solve this equation are supplied from MD simulations.

Not all channels conform to the ED model. This model cannot be applied directly if ion crossing events are not statistically independent [49,50], ion diffusion is single file rather than Fickian [51], there are strong binding sites for ions in the channel, the channel changes its structure in response to applied voltage in the range of interest [26] or ions experience saturation in the mouth of the channel. Despite these limitations, it appears that ion movement through many channels satisfies the assumptions of the ED model. A number of small, naturally occurring and synthetic channels and pLGICs belong to this category. The channels discussed in this paper were found to be well-described by the model. More generally, a number of different equations that are special cases of the ED equation, such as Goldman-Hodgkin-Katz (GHK) equation, have been extensively and successfully used as basic tools in experimental and computational electrophysiology for nearly 80 years, for example to determine ionic selectivity from the measured reversal potential [1,52]. Further, basic assumptions of the ED model, such as independence of ion crossing events or Fickian nature of diffusion, can be tested without substantial, additional effort. This was previously done for a number of channels [25].

Previously, we developed an approach to calculate electrophysiological properties from the integrated form of the ED equation [19,24,25]. Instead of MD calculations at several voltages, the system is simulated at a single voltage (or with no applied voltage) to obtain the PMF for each ion in the channel. Subsequently, markedly shorter simulations at voltages of interest are required to determine the densities of the ions near the ends of the channel. Calculating the currents from simulations at these voltages is not required. In addition, ionic diffusivity along the channel has to be determined. Both boundary density values and diffusivity obtained by any of these methods are burdened with errors, contributing to inaccuracies in the calculated currents.

Recently, we developed two formalisms for calculations of electrophysiological properties, including I-V curves and reversal potential, from a single MD simulation at one voltage [33]. From this simulation, the PMF, nonequilibrium density profiles and committor probabilities for ions in the channel are obtained and used to calculate currents at different voltages after appropriate transformations of the ED equation. Additional calculations to obtain the density boundary terms at different voltages and diffusivity are no longer needed. These formalisms were tested on a simple model of the TTX channel, comprised of 7 straight

α

-helices, each containing 18 amino acids [6], and were shown to perform very well. The improved efficiency of this novel approach derives from the fact that only one MD simulation instead of multiple ones is needed to obtain the I-V curve or reversal potential.

As is the case for any new approach, it is essential to establish the intrinsic accuracy of our formalism. This is the goal of this paper. Specifically, we focus on the question: how reliable are our approaches to calculating electrophysiological properties, independent of other sources of errors, such as inaccuracies in force fields and insufficient simulation times? Separating errors due to the proposed methods from other error is not simple. It cannot be done through a direct comparison with experiments because, for example, of inaccuracies due to force fields. In principle, it can be done via comparison with accurate MD simulations of the PMF and currents at several applied voltages, but, in practice, it is expensive to obtain sufficiently accurate free energies and currents. Although PMFs for ions in a number of channels were obtained from MD simulations, in most cases no estimates of errors were provided [21,23,26,53,54,55,56,57,58,59,60,61]. In a few cases in which errors are available by way of either direct estimates or comparisons of PMF obtained via different methods [20,24,25,62,63,64,65,66,67,68] they are of the order of 0.2–0.7 kcal/mol, which is similar to what is expected to be the intrinsic accuracy of the formalisms studied here. This means that if there were differences between electrophysiological properties obtained from MD simulations at several applied voltages and reconstructed from simulations at a single voltage it would not be possible to determine whether these differences were due to insufficient accuracy of the simulations or to inaccurate reconstruction from the new methods. For these reasons, we take a different approach.

We assume that the ED model describes ion transport with satisfactory accuracy and that the underlying PMF is known. Then, the ED equation is solved many times by way of stochastic simulations to ascertain how statistical errors depend on the number of stochastic trajectories. Even though the stochastic simulations employed here do not involve time explicitly, the number of trajectories considered in a stochastic dataset can be related to the MD simulation time. In a simple case of TTX, we demonstrate that this can be done consistently. For each set of stochastic trajectories, the PMF and the electrophysiological properties at different applied voltages are reconstructed by way of the new theoretical approaches considered here. Due to the stochastic nature of each solution to the ED equation and the limited number of trajectories in each set, which can be related to a specific simulation time, the quantities of interest obtained from each reconstruction differ between themselves and from the accurate values associated with the underlying PMF. Further, the calculated quantities also depend on the theoretical formalism used. If a sufficient number of trajectories has been generated, statistical errors on the quantities of interest can be estimated as a function of the number of trajectories or equivalently, simulation time and performance of each theoretical approach can be systematically assessed. For sets with a large number of trajectories, which corresponds to long simulation times in MD, the underlying PMFs will be reconstructed accurately. For sets with a smaller number of trajectories, the accuracy will not be as good and is expected to deteriorate as the number of trajectories is reduced. A similar systematic study cannot easily be done in practice by way of MD simulations because the computational effort to generate many MD trajectories of different length would have been prohibitive. Furthermore, no analytical method for error analysis exists for this problem.

In principle, this type of analysis can be carried out for any underlying PMF, even if it is unrelated to real ion channels. This is, however, not the direction that we pursue. Instead, we use the PMFs that we previously obtained from simulations of three actual channel models and, for the purpose of this study, assume that they are accurate. The models were selected such that they differ in size, pore structure, conductance and selectivity. The first model is TTX, which exhibits relatively high conductance, very little structure inside the pore and a weak selectivity for cations.

The second model is the high-resolution NMR structure proposed by OuYang et al. [69] for a hexameric channel p7 from the hepatitis C virus. Each subunit consists of 63 amino acids. The model has an unusual architecture not found in any other channel. The channel does not exist as a bundle of

α

-helices, which is the most common structural motif among membrane proteins, but instead forms an interlocked structure in which each subunit assumes a horseshoe conformation with each side comprised of a short,

α

-helical section. Because of these atypical features there have been concerns about the veracity of this model [28]. Recently, we calculated conductance and ionic selectivity of this model by way of MD simulations and showed that both properties differ significantly from those measured experimentally (Shannon et al., unpublished). Specifically, in contrast to the electrophysiological data, the model exhibits high conductance and strong selectivity for Cl

^{-}

over K

^{+}

. These results strongly suggest that the proposed model does not represent the native structure of the channel, demonstrating that computational electrophysiology can be used not only to support but also to disprove structural models of ion channels.

The third model is based on the crystal structure of a pentameric, cation-selective ion channel, GLIC, from a cyanobacterium Gloeobacter [70]. This channel is a bacterial homolog of receptors belonging to the family of pentameric ligand-gated ion channels. Its main electrophysiological characteristics are low-conductance (9.3 pS) and strong selectivity for cations [70]. Molecular models of all three channels are shown in the Supplementary Materials (SM), Section S4.

Both the p7 and GLIC models have a markedly more varied pore structure than TTX and, consequently, a more complex PMF. Although we will use the names of these three channels further in the text, we do not claim that the underlying PMFs faithfully represent the PMFs for these channels in their native open forms (this does not appear to be the case for p7) and, therefore, we do not compare the electrophysiological properties calculated from stochastic simulations with the same properties obtained experimentally. Instead, we fully concentrate on assessing the accuracy of the underlying theory.

2. Theory and Method

In this section, we briefly outline the theory behind three different approaches to calculating the PMF and electrophysiological characteristics of an ion channel. Two of them require simulations at only one applied voltage. A more detailed derivation of the basic formulas, which follows closely the earlier development [33], is provided in Supplementary Materials, Section S1. Next, we describe how the properties of interest can be obtained from stochastic simulations under the assumptions of the ED model specified in the introduction. Note that while the theory is developed in the context of MD simulations, here, we use the results of the theory to compute the PMF and I-V curves from density profiles and committor probabilities that were obtained from stochastic simulations.

2.1. Calculating the Potential of Mean Force

If the concentrations of ions on both sides of the membrane and the applied voltage remain constant in time, the system is in a steady state, which means that the flux of ions through the channel, J, is also constant in time. These are the conditions most often considered in both experiments and simulations aimed at extracting electrophysiological properties of channels. Then, the one-dimensional ED equation for a given type of ions can be written as

J = - D (z) (\frac{d ρ (z)}{d z} + β ρ (z) \frac{d E (z)}{d z}),

(1)

where

D (z)

is the diffusivity that, in general, depends on position z along the reaction coordinate

z

. For a transmembrane channel embedded in the membrane located in the x,y-plane, a convenient reaction coordinate is the position of an ion along the pore of the channel, which can be measured along the z-coordinate.

ρ (z)

is the line density of ions, which is usually recorded as a histogram in computer simulations.

β = k_{B} T

, where

k_{B}

is the Boltzmann constant and T is temperature.

E (z)

is given by

E (z) = A (z) + q V (z),

(2)

where

A (z)

is the PMF,

V (z)

is the applied voltage and q is ionic charge. In a constant electric field,

E^{e l}

, acting along

z

, which is the most frequent experimental condition,

V (z) = E^{e l} (z - z_{a}) .

(3)

Even though the electric field is applied across the whole system [15,30], it acts only between

z_{a}

and

z_{b}

in the non-polar phase, which has been identified as corresponding to the hydrophobic core of the membrane [25,33]. Thus, electric field is a boxcar function that is equal to

E^{e l}

in the range [

z_{a}, z_{b}

] and zero otherwise. This can be formally written as

E^{e l} [H (z - z_{a}) - H (z - z_{b})]

, where H is the Heaviside function. Although we will not use this notation for simplicity, the range in which

E^{e l}

is non-zero has to be kept in mind.

Integrated with the integrating factor

exp [β E (z)]

and resolved for J, the ED equation takes the form

J = \frac{ρ (z_{m i n}) exp [β E (z_{m i n})] - ρ (z_{m a x}) exp [β E (z_{m a x})]}{\int_{z_{m i n}}^{z_{m a x}} \frac{exp [β E (z)]}{D (z)} d z} .

(4)

For a system in a steady state, J does not formally depend on the limits of integration

z_{m i n}

and

z_{m a x}

. This means that these limits do not have to coincide with the edges of the channel. In practice, the limited precision of MD simulations introduces some dependence on the limits of integration, as analyzed elsewhere [25].

To calculate J from this equation,

E (z)

has to be known, which in turn requires determining

A (z)

. This can be done in equilibrium simulations in the absence of voltage. A host of methods exist for this purpose [71,72,73].

A (z)

can be also calculated from non-equilibrium simulations at an applied voltage. If the ED equation is integrated with the integrating factor

1 / ρ (z)

then

J = \frac{ln \frac{ρ (z_{m i n})}{ρ (z_{m a x})} - β [A (z_{m a x}) - A (z_{m i n}) + q E^{e l} (z_{m a x} - z_{m i n})]}{\int_{z_{m i n}}^{z_{m a x}} \frac{1}{D (z) ρ (z)} d z} .

(5)

Since J is independent of the limits of integration,

z_{m a x}

can be substituted by z. After simple rearrangements, it yields a formula for the PMF relative to its value at

z_{m i n}

Δ A (z, z_{m i n}) = A (z) - A (z_{m i n})

Δ A (z, z_{m i n}) = - k_{B} T [ln \frac{ρ (z)}{ρ (z_{m i n})} + J \int_{z_{m i n}}^{z} \frac{1}{D (z^{'}) ρ (z^{'})} d z^{'}] - q E^{e l} (z - z_{m i n}) .

(6)

We call this method for determining PMF the Integrated Electrodiffusion Equation Method (IEEM).

To solve Equations (5) and (6),

D (z)

has to be known in the full range of z.

D (z)

can be determined by way of calculating the mean square displacement of the ion at several points along the channel obtained from a series of short MD trajectories after subtracting the PMF [19], from the force-force autocorrelation function acting on a stationary ion at different positions in the channel [74], or by way of a Bayesian fitting method [75,76,77]. See Supplementary Materials, Section S5 for a discussion of how diffusivity was computed in our MD simulations.

Once

Δ A (z_{m a x}, z_{m i n})

and

D (z)

, which are both assumed to be independent on voltage, are known, the boundary density terms

ρ (z_{m i n})

and

ρ (z_{m a x})

have to be obtained from either MD or stochastic simulations at each voltage of interest. Since the full knowledge of

ρ (z)

is not needed, these simulations can be markedly shorter than simulations to determine the PMF. Then,

Δ A (z_{m a x}, z_{m i n})

D (z)

ρ (z_{m i n})

and

ρ (z_{m a x})

are used in Equation (5) to calculate J at a given voltage. Previously, we demonstrated that this method performs satisfactorily for simple channels [19,24,25].

Recently, we developed two alternative approaches to calculating the PMF and electrophysiological properties that require markedly less computational effort [33]. Both rely on separating the total ionic current, J, to currents moving in two opposite directions – from

z_{m i n}

z_{m a x}

and from

z_{m a x}

z_{m i n}

. We abbreviate them

J^{f}

and

J^{b}

and call them forward and backward currents, respectively.

J^{f} = \frac{ρ^{f} (z_{m i n}) exp [β E (z_{m i n})] - ρ^{f} (z) exp [β E (z)]}{\int_{z_{m i n}}^{z} \frac{exp [β E (z^{'})]}{D (z^{'})} d z^{'}},

(7)

J^{b} = \frac{ρ^{b} (z) exp [β E (z)] - ρ^{b} (z_{m i n}) exp [β E (z_{m i m})]}{\int_{z_{m i n}}^{z} \frac{exp [β E (z^{'})]}{D (z^{'})} d z^{'}},

(8)

J = J^{f} - J^{b} .

(9)

Here,

ρ^{f} (z)

and

ρ^{b} (z)

are densities of ions that entered the range

[z_{m i n}, z_{m a x}]

z_{m i n}

and

z_{m a x}

, respectively. We assume that both forward and backward currents are in a steady state and, therefore, their values do not depend on the limits of integration. This allows for setting the upper limit to

z_{m i n} < z \leq z_{m a x}

Assume that

J^{f} > 0

and

J^{b} > 0

and take the ratio of Equation (7) to Equation (8). This yields

\frac{J^{f}}{J^{b}} = \frac{ρ^{f} (z_{m i n}) - ρ^{f} (z) exp [β Δ E (z, z_{m i n})]}{ρ^{b} (z) exp [β Δ E (z, z_{m i n})] - ρ^{b} (z_{m i n})},

(10)

where

Δ E (z, z_{m i n}) = E (z) - E (z_{m i n}) .

(11)

Combined with Equations (2) and (3), Equation (10) can be solved for

Δ A (z, z_{m i n})

Δ A (z, z_{m i n}) = - k_{B} T ln \frac{J^{b} ρ^{f} (z) + J^{f} ρ^{b} (z)}{J^{b} ρ^{f} (z_{m i n}) + J^{f} ρ^{b} (z_{m i n})} - q E^{e l} (z - z_{m i n}) .

(12)

From this equation it follows that the PMF can be obtained from non-equilibrium simulations at applied electric field

E^{e l}

simply from an average of ion densities in the forward and backward directions weighed by the backward and forward currents, respectively. We call this method for calculating the PMF the Current-Weighted Density Method (CWDM). Knowledge of diffusivity is not necessary in CWDM. The denominator in the argument of the logarithmic function sets the reference value of the PMF at

z_{m i n}

If we abbreviate the number of crossing events in forward and backward direction as

n^{f}

and

n^{b}

, respectively, then, assuming that crossing events are governed by the Poisson statistics, the corresponding errors will be approximately

1 / \sqrt{(} n^{f})

and

1 / \sqrt{(} n^{b})

. This means that if

n^{f}

n^{b}

is small,

Δ A (z, z_{m i n})

calculated from Equation (12) may become inaccurate. Thus, we developed another, related theoretical approach for determining the PMF from non-equilibrium simulations that does not suffer from this disadvantage. Since it requires calculating committor probability,

P (z)

, we will call it the Committor Probability Method (CPM). For a diffusive process considered here,

P (z)

referenced to the forward direction is defined as the probability that a particle (ion) in position z will reach

z_{m a x}

before it reaches

z_{m i n}

P (z)

can be calculated either directly during computer simulations or in post-processing, as described in Supplementary Materials, Section S3. A general discussion of committor probabilities in more than one dimension and their application to chemical kinetics can be found elsewhere [78,79,80,81].

The PMF can be calculated from ion densities in the forward or the backward direction. The corresponding formulas are

exp [β Δ E (z, z_{m i n})] = \frac{ρ^{f} (z_{m i n}) [1 - P (z)]}{ρ^{f} (z) - ρ^{f} (z_{m a x})},

(13)

exp [β Δ E (z, z_{m i n})] = exp [β Δ E (z_{m a x}, z_{m i n})] \frac{ρ^{b} (z_{m a x}) P (z)}{ρ^{b} (z) - ρ^{b} (z_{m i n})} .

(14)

Their derivation closely follows our earlier work [33] and is given in Supplementary Materials Section S1.

Both equations allow for calculating the same quantity— the PMF. Individually, each of them is not expected to be accurate in the full

[z_{m i n}, z_{m a x}]

range of z, especially away from the entry point. Specifically, as z becomes close to

z_{m a x}

both

ρ^{f} (z) - ρ^{f} (z_{m a x})

and

1 - P (z)

approach zero. Since numerical inaccuracies in Equations (13) and (14) affect mainly the opposite sides of the

[z_{m i n}, z_{m a x}]

range, these two equations can be profitably combined. Then,

ρ^{f} (z) - ρ^{f} (z_{m a x})

and

ρ^{b} (z) - ρ^{b} (z_{m i n})

can be considered as two biased distributions representing the same unbiased distribution

h (z)

. The problem of merging them to reconstruct

h (z)

such that statistical error on

Δ A (z, z_{m i n})

is minimized can be solved by way of the Weighted Histogram Analysis Method (WHAM) [82]. This yields the following formula for reconstructing the PMF from non-equilibrium simulations:

Δ A (z, z_{m i n}) = C - k_{B} T ln [\frac{h (z)}{ρ^{f} (z_{m i n}) (1 - P (z)) P (z)}] - q E^{e l} (z - z_{m i n}),

(15)

where neither C, which is a constant that only shifts the energy scale, nor

ρ^{f} (z_{m i n})

, which is independent of z and is needed to ensure that the PMF at

z = z_{m i n}

is equal to zero, influences the shape of

Δ A (z, z_{m i n})

. Similarly to Equation (12), no knowledge of diffusivity is required.

Typically, MD simulations would be carried out on the channel system of interest at some applied voltage V. From this, the committor probability,

P (z)

and the 1-sided density profiles,

ρ^{f} (z)

and

ρ^{b} (z)

, and the number of forward and backward crossing events would be determined, and used to calculate the forward and backward fluxes,

J^{f}

and

J^{b}

, respectively. The PMF can be determined from either CWDM (Equation (12)) or CPM (Equation (15)). As will be discussed later, we created synthetic data sets of

10^{6}

–

10^{8}

stochastic trajectories. For each data set, we calculate the same quantities that would be calculated in MD,

P (z)

ρ^{f} (z)

and

ρ^{b} (z)

, and then use these to calculate the PMF. Note that the free energy depends on ratios of density profiles, so the absolute normalization of the density profiles is not important. Similarly, the CWDM requires only ratios of the forward and backward currents, so the magnitudes are not required.

2.2. Calculating I-V Dependence from Simulation at a Single Voltage

If the PMF, the current,

J_{μ}

, and the density,

ρ_{μ} (z)

, or the committor probability,

P_{μ} (z)

, are known from simulations at an applied voltage,

Δ V_{μ}

, the current,

J_{ν}

, at a different voltage

Δ V_{ν}

can be obtained without any calculations at this voltage. This allows for reconstructing the I-V curve from simulations at a single voltage.

If Equation (1) is integrated with the same integrating factor,

exp [β E_{ν} (z)]

, for both voltages,

Δ V_{μ}

and

Δ V_{ν}

, we obtain

\begin{matrix} J_{μ} = \\ \frac{ρ_{μ} (z_{m i n}) exp [β E_{ν} (z_{m i n})] - ρ_{μ} (z_{m a x}) exp [β E_{ν} (z_{m a x})] + β q (E_{ν}^{e l} - E_{μ}^{e l}) \int_{z_{m i n}}^{z_{m a x}} ρ_{μ} (z) exp [β E_{ν} (z)] d z}{\int_{z_{m i n}}^{z_{m a x}} \frac{exp [β E_{ν} (z)]}{D (z)} d z}, \end{matrix}

(16)

and

J_{ν} = \frac{ρ_{ν} (z_{m i n}) exp [β E_{ν} (z_{m i n})] - ρ_{ν} (z_{m a x}) exp [β E_{ν} (z_{m a x})]}{\int_{z_{m i n}}^{z_{m a x}} \frac{exp [β E_{ν} (z)]}{D (z)} d z}

(17)

The latter but not the former equation is the standard integrated form of the ED equation, Equation (4).

If we take the ratio of currents

J_{μ} / J_{ν}

then, after some algebra given in Supplementary Materials, Section S2 we obtain

\frac{J_{μ}}{J_{ν}} = 1 + β q (E_{ν}^{e l} - E_{μ}^{e l}) \int_{z_{m i n}}^{z_{m a x}} f_{μ} (z) exp \{β q [V_{ν} (z) - V_{μ} (z)]\} d z,

(18)

where

f_{μ} (z) = exp [β Δ E_{μ} (z, z_{m i n})] \frac{ρ_{μ} (z)}{ρ_{μ} (z_{m i n})} .

(19)

f_{μ} (z) = 1 + exp [β Δ E_{μ} (z, z_{m i n})] \frac{ρ_{μ} (z_{m a x})}{ρ_{μ} (z_{m i n})} - P_{μ} (z),

(20)

depending on whether it is preferred to calculate

J_{ν}

from ion density,

ρ_{μ} (z)

, or committor probability,

P_{μ} (z)

. In both instances, neither diffusivity nor quantities at the applied voltage

Δ V_{ν}

are needed. Equation (19) is expected to be less accurate that Equation (20) because

Δ E_{μ} (z, z_{m i n})

and

P_{μ} (z)

that enter the latter equation are estimated on the basis of both forward and backward simulations, whereas

ρ_{μ} (z)

in the former equation is a one-sided density that looses accuracy away from the entry point.

Unlike the free energy, Equation (18) gives only ratios of forward or backward currents with respect to a reference voltage. Consequently, to calculate the I-V curves, we need currents at this reference voltage.

2.3. Stochastic Simulations

The electrodiffusion equation was solved by generating trajectories on a free energy surface

E (z)

that included the PMF and applied electric field with diffusivity

D (z)

or average diffusion coefficient,

< D >

, at temperature T [83,84]. This allowed us to generate the channel crossing statistics, density profiles and committor probabilities for the ions for this free energy surface. As the crossing events that we have observed in MD simulations appear to obey Poisson statistics, independently for both ions, we consider the ED equation for each ion separately. Then, we calculated statistical errors in recovering the underlying PMF and the I-V curves as functions of the number of trajectories.

As above, we define the channel boundaries as

z_{m i n}

and

z_{m a x}

, and absorbing boundary conditions were located at these points. Trajectories were initiated at a point just inside the boundary at either

z_{m i n}

for forward trajectories or

z_{m a x}

for backward trajectories, and propagated until they reached either of the absorbing boundaries. Forward and backward trajectories are considered separately as we are interested in the 1-sided density profiles and committor probabilities, as well as their 2-sided combination. A trajectory that crossed from

z_{m i n}

z_{m a x}

is said to be a crossing trajectory in the forward direction. Similarly, trajectories that cross from

z_{m a x}

z_{m i n}

are crossing trajectories in the backward direction. For simplicity, these will be referred to as forward or backward crossing events. Since the trajectories are initiated near the absorbing boundaries, the majority of trajectories in either direction do not cross, but they do contribute to the 1-sided density profiles and first-passage statistics that are used to compute the committor probabilities (see Supplementary Materials, Section S3 for further details).

The number of trajectories initiated per data set were typically

10^{6}

10^{7}

10^{8}

, further abbreviated as N6, N7 and N8, respectively. These numbers were chosen because the average number of crossing events for PMFs corresponding to the models of TTX and p7 observed for

N = 10^{6}

is of the same order of magnitude as the numbers of crossing events observed in our MD simulations of 0.5–2

μ

s. For the cation-selective GLIC channel, in which the free energy barrier to permeation of Na

^{+}

is markedly higher, simulations with

10^{6}

trajectories yielded too few crossing events to be useful. In this case, the number of crossing events observed at the N7 level approximately corresponds to the number seen in MD simulations of 8

μ

s. See Supplementary Materials, Section S4 for details of the MD simulations.

The free energy surfaces for the stochastic simulations were obtained by adding the voltage ramp to the PMFs. We use a set of PMFs from our MD calculations. For our problem, the PMFs are the equilibrium free energy surfaces for moving an ion along the 1-dimensional reaction coordinate of the ion with respect to the center-of-mass of the protein channel, at the bulk ion densities of the MD. For this study, we used average diffusion coefficients,

< D >

obtained by averaging the diffusivities estimated from MD (Supplementary Materials, Section S4). Strict matching of diffusivity is not necessary for the primary purpose of this study, but it provides a more realistic connection between statistical errors estimated in stochastic simulations and time scales of MD simulations. Additional details are given in Supplementary Materials, Section S7. The forward and backward ion density profiles were obtained from histograms of either the forward or backward trajectories in each data set. Committor probabilities (Supplementary Materials, Section S3) were calculated from the first passage statistics of the forward and backward trajectories. The density histograms and committor probabilities were computed for each data set, and not as an average over the individual trajectories in the data set. Averages were constructed over multiple data sets.

3. Results and Discussion

3.1. Connection with Molecular Dynamics

To compute the currents for the I-V curves from stochastic simulations, some connection to MD is required. MD simulations provide both forward and backward ion trajectories as part of the simulation, unless the channel is strongly rectifying or a large voltage is applied. The net current due to a particular ion is

J = J^{f} - J^{b}

(Equation (9)). The total current is the sum of these net currents over all types of ions. As mentioned in the introduction, J can be obtained from MD simulations by way of combining the fluxes from forward and backward crossing events or calculating the ionic displacement currents. In stochastic simulations, only the former method can be used. Therefore, we tested whether both method yield the same results for MD and found that this was indeed the case. Both methods and the results of the tests are described in Supplementary Materials, Section S6.

In addition, the detailed balance condition connecting forward and backward crossing events has to be satisfied. In MD simulations this problem is implicitly solved: if there is no external voltage, simulations of transmembrane systems will exhibit no net current to within statistical errors, which means that the number of forward and backward crossing events is equal, again to within statistical errors. In stochastic simulations, detailed balance also has to be satisfied, which means that trajectories in both directions have to be combined with the correct weights. To determine these weights, we carried out sets of

10^{8}

simulations with no applied voltage to obtain the well converged, average numbers of forward and backward crossing events. From these simulations, the ratio of forward to backward trajectories that satisfies the detailed balance condition was established and subsequently used to compute the density profiles, committor probabilities and the PMFs at different voltages.

Once the ratio of forward to backward trajectories needed to satisfy the detailed balance condition is known, the average numbers of crossing events in both directions,

n^{f} (Δ V)

and

n^{b} (Δ V)

, can be obtained from stochastic simulations at a given voltage

Δ V

. This, however, is still insufficient to determine currents; additional information about time scales is required. This can be obtained from a MD simulation of the system. We abbreviate the number of forward and backward crossing events observed in MD simulations at applied voltage,

Δ V_{r e f}

, as

m^{f}

and

m^{b}

. Then, the length of the MD trajectory,

t_{M D}

, can be used to estimate a stochastic time,

t_{S}

, corresponding to the number of stochastic trajectories that produced

n^{f} (Δ V_{r e f})

and

n^{b} (Δ V_{r e f})

crossing events at the voltage

Δ V_{r e f}

. A simple way to make such estimate is to use the number of crossing events in one direction. It is, of course, recommended to choose the direction that provides better statistics. Assuming that there are more forward than backward crossing events in the MD simulations,

t_{S} = t_{M D} \frac{n^{f} (Δ V_{r e f})}{m^{f}} .

If the backward events dominate,

t_{S}

would be estimated using

n^{b} (Δ V_{r e f})

and

m^{b}

. Once

t_{S}

has been determined, the stochastic currents at voltage

Δ V

can be calculated:

J_{S} (Δ V) = J_{S}^{f} (Δ V) - J_{S}^{b} (Δ V) = [n^{f} (Δ V) - n^{b} (Δ V)] / t_{S} .

where

J_{S} (Δ V)

J_{S}^{f} (Δ V)

and

J_{S}^{b} (Δ V)

are the total, forward and backward currents at voltage

Δ V

3.2. Committor Probabilities

The committor probabilities for p7 are shown in Figure 1. The committor probabilities in Figure 1a have been calculated from Supplementary Materials, Equation S30, in which ions arriving at z from both sides are included. The statistical errors associated with

P (z)

at different voltages are small, even at the N6 level. As can be seen in Figure 1b, this is not the case for one-sided

P (z)

, obtained using Supplementary Materials, Equations S28 and S29. This is due to the decreasing number of ions from one direction as they approach the opposite side of the channel. The inset of Figure 1b shows the numbers of first-passage events from the forward and backward calculations as well as the combined number of events. At 140 mV, the statistics are satisfactory only in the forward direction, in which most of crossing events occur, whereas no reliable probabilities are obtained for the backward direction over the full range of z. The opposite is true for

- 140

mV;

P (z)

in the forward direction is unreliable. Thus, combining information about

P (z)

in both directions is preferable whenever possible.

As voltage changes from

- 140

mV to 140 mV, the position of the transition state for K

^{+}

permeation through p7, defined as the x,y-plane at which

P (z) = 0.5

, Ref. [85] shifts substantially and systematically from 7 Å to

- 13

Å with respect to the center of mass of the membrane. Such large shifts, however, are not universal. As we have shown in the example of TTX [33], the position of the transition state changes markedly less with voltage if the underlying PMF is strongly peaked.

Calculating

P (z)

for GLIC is more difficult. This is a slow channel and even at the N7 level, which approximately corresponds to a MD trajectory of 8

μ

s in length (see Supplementary Materials, Section S4, the number of crossing events is small. At 100 mV only an average number of 0.5 forward and 29 backward crossing events were observed. In particular, N7 simulations of forward trajectories frequently produce no crossing events. At the same voltage,

P (z)

in the backward direction is often equal to 1 over a relatively wide range of several Å near

z_{m a x}

, which means that all ions that reached this range exit the channel at

z_{m a x}

. In such circumstances, calculation of

P (z)

from Supplementary Materials, Equation S30 is no longer possible. A different approach is needed.

Direct calculation of the committor probability requires that some number of trajectories successfully cross the channel,

N^{b} (z_{m i n})

> 0 and

N^{f} (z_{m a x})

> 0. If one of these conditions is not met, for example, if

N^{f} (z_{m a x})

= 0, then

P^{f} (z) = N^{f} (z_{m a x}) / N^{f} (z)

= 0. If we consider position

z^{'}

(

z^{'}

z_{m a x}

) at which

N^{f} (z^{'})

> 0, then we can write the committor probability

P (z)

α N^{f} (z^{'}) / N^{f} (z)

, where

α

is unknown, though formally would be equal to

N^{f} (z_{m a x}) / N^{f} (z^{'})

if complete sampling of the forward direction were available.

α

can be determined in a self-consistent manner. Using Equation S30 from Supplementary Materials, Section S3, we can write the total committor probability in the region

z < z^{'}

and the backward committor probability

z > z^{'}

P (z) = \{\begin{matrix} \frac{α N^{f} (z^{'}) + N^{b} (z) - N^{b} (z_{m i n})}{N^{f} (z) + N^{b} (z)} & if z < z^{'} \\ 1 - \frac{N^{b} (z_{m i n})}{N^{b} (z)} & if z > z^{'} . \end{matrix}

If we require that

P (z)

is continuous at

z = z^{'}

, then

α

= 1 −

N^{b} (z_{m i n}) / N^{b} (z^{'})

. Other ways of determining

P (z)

for this problem are possible.

3.3. The Potential of Mean Force

Typically, the PMFs for ions in channels are calculated in simulations in the absence of electric field using enhanced sampling techniques (see the recent review by Flood, et al. [27]). In contrast, the methods outlined here allow for reconstructing PMF from steady-state simulations with an applied electric field. The underlying PMFs for p7 and GLIC used in the present study were obtained by way of this method (see Supplementary Materials, Section S1). Since TTX is a bundle of straight

α

-helices surrounding a featureless water pore, the PMFs for K

^{+}

and Cl

^{-}

are quite generic, which is characteristic of several very simple channels (see Figure 2a) [24,25]. For K

^{+}

, the PMF is fairly flat over a wide range of approximately 18 Å inside the channel, which is reminiscent of classical models of ionic conductance in which it is assumed that the PMF is a step function constant inside the channel [1,86]. For Cl

^{-}

, the PMF is peaked near the center of the bilayer, which can be attributed to the Born barrier experienced by an ion permeating a rigid, featureless non-polar lamella [87]. If an ion is transferred across a membrane through a water-filled pore, the general shape of the PMF remains the same, but the barrier is substantially reduced [87,88]. For TTX, it still remains approximately 1.5 kcal/mol higher than the barrier for K

^{+}

, which is consistent with a weak selectivity of this channel toward cations.

The PMF for permeation of Cl

^{-}

in the OuYang et al. model of p7 [69] is more structured than the PMF for TTX (see Figure 3a). The barriers are low, which explains high chloride current predicted by this model [37]. In contrast to TTX, the barriers to Cl

^{-}

permeation in p7 are located near the mouths of the channel due to the presence of positively charged residues at these locations. Compared to permeation of Cl

^{-}

, the current of K

^{+}

in this model is quite low, which indicates that the channel should be anion-selective. Both predicted selectivity and total currents are at variance with electrophysiological data [89], thus contributing to arguments that the proposed high-resolution structure [69] is not native.

The PMF representing permeation of Na

^{+}

through GLIC is also markedly more structured than the PMF for ions in TTX (see Figure 2b). The barrier is substantially higher than in the other two channels. As a result, the conductance of this channel is relatively low [70]. This presents a challenge because the number of crossing events in both directions is small. The PMF for Cl

^{-}

in this channel is not considered because no crossing events of this ion have been observed in MD simulations.

Taken together, the PMFs considered here are quite different from one another, but are typical of the variety seen in ion channels. In spite of these differences, all three PMFs were successfully reconstructed from non-equilibrium simulations by way of Equations (12) and (15) associated, respectively, with the CWDM and CPM methods. The applied voltages were 50 mV for TTX, 140 and

- 35

mV for p7 and 100 mV for GLIC. For TTX and p7, reconstruction was carried out at the N6, N7 and N8 levels. For GLIC, the number of crossing events at the N6 level was quite small or equal to zero. Thus, only N7 and N8 levels were considered. At the N6 level, 50 and 100 data sets of trajectories were generated for TTX and p7, respectively. At the N7 level, 20 sets of trajectories were generated for TTX and p7, and 50 sets were generated for GLIC. At the N8 level, the number of generated sets was 4, 8 and 25 for TTX, p7 and GLIC, respectively. At each level, all reconstructed PMFs are found to be tightly clustered and their averages at each level are close to the underlying PMF, as shown in Figure 2 and Figure 3a.

In these figures, statistical errors associated with dispersion of the reconstructed PMFs are marked. For

Δ A (z_{m a x}, z_{m i n})

, these errors are approximately ±0.3 kcal/mol at the N6 level for TTX and p7 and at the N7 level for GLIC. As expected, they are reduced by approximately a factor of 3 with each level in which the number of sets increases by an order of magnitude. The mean PMFs obtained by way of CWDM and CPM at different levels are quite close to the underlying PMF and the corresponding statistical errors are very similar, indicating that both methods are successful in reproducing the underlying PMFs. Only for GLIC at the N7 level, does the

Δ A (z_{m a x}, z_{m i n})

reconstructed by way of CWDM appear to be systematically underestimated. At this level, no crossing events in one direction are observed for a considerable fraction of data sets, which makes reconstruction of the PMF from Equation (12) impossible. This systematically biases the sample in favor of sets with higher counts of crossing events and, consequently, lower

Δ A (z_{m a x}, z_{m i n})

. From the comparison between the PMFs reconstructed for p7 from trajectories at 140 and −35 mV, it appears that precision of the reconstruction depends somewhat on applied voltage. If the forward and backward densities are well balanced, precision improves.

In CPM, the PMFs can be calculated from one-sided quantities, Equations (13) and (14), or by combining them. Here, the latter has been done by way of WHAM, Equation (15). As can be seen in Figure 3b, this approach yields improved agreement with the underlying PMF. In one-sided formulas, the densities can become quite low near the exit and, as a result, precision in this range suffers.

In summary, both CWDM and CPM provide a reliable means for reconstructing PMFs from non-equilibrium simulations. However, the relation between statistical errors obtained in stochastic and MD simulations is not straightforward. Even if the assumptions of the ED model are satisfied, precision of stochastic simulations is expected to be higher than precision in MD simulations of equivalent length. Specifically, it is usually uncertain if all degrees of freedom perpendicular to the reaction coordinate have been properly equilibrated on the time scale of the simulations. Torsional angles in the side chains of residues lining the pore or motion of whole helices are examples of degrees of freedom that might undergo slow equilibration and, by doing so, influence the calculated PMF and electrophysiological properties. The same concern applies to all other methods for calculating these quantities.

3.4. Current-Voltage Dependence

Once the PMFs for the ions permeating the channel have been reconstructed and the committor probabilities for these ions have been calculated for a reference voltage, the full current-voltage (I-V) curves can be calculated from Equations (19) and (20) without the need for additional simulations. This is the principal gain in efficiency of the method: the I-V curve can be obtained from a single MD simulation instead of multiple simulations. For example, if constructing the I-V curve required MD simulations at five different voltages in the range of interest the efficiency of our methods would be approximately five-fold. Since numerical results indicate that Equation (19) yields less accurate results than Equation (20), this equation will not be further considered. The results for TTX, p7 and GLIC are shown in Figure 4 and Figure 5. The reference applied voltages are the voltages used for reconstructing PMFs, described in the previous subsection. For comparison, currents calculated directly from stochastic trajectories at several voltages are also shown.

As we can see in Figure 4, the agreement between the I-V curves for both K

^{+}

and Cl

^{-}

in TTX calculated directly and by way of Equations (18) and (20) is excellent, even at the N6 level, for the full range of voltages studied here, which extends from −100 to 100 mV. As shown in Figure 6, the I-V curves obtained for different sets of trajectories are closely clustered and deviate from each other only at the largest absolute applied voltages by no more than a few pA.

For p7, the agreement is not as good if the the reference voltage of 140 mV is used for calculating the I-V curve. The Cl

^{-}

currents calculated directly and from Equation (18) agree well for positive voltages, but diverge for negative voltages, away from the reference state. The corresponding statistical errors also increase and become quite large below −50 mV. The source of this disagreement can be traced to the integrand in Equation (18). As the difference between the reference and the target voltage increases, the exponential term also increases, which magnifies inaccuracies in function

f (z)

. If the reference voltage is chosen to be

- 35

mV, the agreement over the full range of voltages −150 to 150 mV improves markedly, with modest deviations only at high, positive voltages (see Figure 4). A similar situation was observed for GLIC. For the reference voltage of 100 mV, the I-V curves at the N7 level satisfactorily reproduce currents calculated directly for positive voltages. For negative voltages, the performance of the method progressively deteriorates. Again, if one is interested in an I-V curve that extends to both positive and negative voltages, a different choice of reference voltage may yield significant improvements in accuracy.

As pointed out in the introduction, a number of previous studies have used unrealistically large applied voltages to increase the number of crossing events and, by doing so, improved precision of the calculated currents [15,23,35,37]. Furthermore, as discussed earlier, this may lead to electroporation of the membrane, saturation effects during the intake of ions at the mouth of the channel and involves extrapolation or interpolation to the voltages by way of ad hoc procedures of unknown accuracy. The approach developed here is more efficient and accurate and has a substantially stronger theoretical basis than procedures used previously, even though only calculations at the reference applied voltage are necessary. In this approach, the accuracy of the reconstructed I-V curves can be substantially improved through a judicious choice of this reference voltage. This choice depends on the range of voltage that is of interest and on several properties of a channel, in particular its rectification, which characterizes an asymmetry of currents in response to the change in direction of applied voltage. In general, maximizing total ion current through applying high voltage is not the optimal strategy. Instead, it is often better to choose a voltage that yields good statistics in both directions.

3.5. Reversal Potential

The reversal potential,

Δ V_{R}

, is the applied voltage at which there is no net current. If ionic concentrations on both sides of the membrane are equal,

Δ V_{R} = 0

. Experimentally, the reversal potential is measured by maintaining different concentrations on the cis and trans side of the membrane, and then used in conjunction with the GHK equation to estimate channel selectivity [1,6]. In MD, asymmetric concentrations have to be maintained to measure directly the reversal potential [90], which markedly complicates simulations. We have only considered the situation where the concentrations of ions are the same on both sides of the membrane, as this corresponds to the conditions under which we have carried out MD. We wish to expand this to a range of concentrations.

We expect that the net number of crossing events, from which we calculate the I-V curve, depends on this concentration. If the bulk concentrations are low and the channel is not saturated, then we expect the number of crossing events, and hence the currents, to be linearly dependent on the concentrations of ions. For example, if the concentration is doubled on both sides of the membrane, the net currents will also double. Under these assumptions, we can calculate the reversal potential from our formalism. We simply need to scale the fluxes of all ion types on one side of the membrane to match the desired concentration difference.

The K

^{+}

/Cl

^{-}

selectivity of TTX obtained from the MD simulations is 2.2 [24,33]. Using currents scaled by 5:1 in the I-V reconstruction from Equation (20), we obtain a reversal potential of −9 mV. This corresponds to a GHK selectivity of 1.7, which is reasonably close to the selectivity found in MD. Note that experimentally, the reversal potential is −27 mV, corresponding to a K

^{+}

/Cl

^{-}

selectivity of 6 estimated from the GHK Equation [6]. This cannot be compared directly to our results because the actual channel structure is unknown, there are uncertainties due to force fields, and the GHK equation itself is an approximation.

4. Conclusions

Stochastic simulations were used to investigate the reliability of two new methods to calculate PMFs for ion transport across transmembrane ion channels and electrophysiological properties of these channels within the general framework of the electrodiffusion model. Both methods have the desirable features that only simulations at a single voltage are needed and information on the diffusivity is not required. In CPM, knowledge of the committor probability is required. Stochastic simulations containing 10

^{6}

trajectories were shown to have similar numbers of crossing events for models of TTX and p7 in MD simulations of 1–2

μ

s in length. Analysis of 50 or 100 of such simulations indicate that errors in the free energy profiles are approximately ±0.3 kcal/mol. For a model of a slow channel, GLIC, 10

^{7}

trajectories, which approximately corresponds to MD simulations of 10

μ

s in length, are needed to achieve a similar statistical error. For both TTX and for p7 at lower applied voltages, CPM and CWDM yield similar results. In CWDM, one-sided fluxes are used directly, and for cases in which few crossing events are observed in one direction, either due to large applied voltages, such as p7 at ±140 mV, or because the channel is rectifying, such as GLIC, CPM performs better because two-sided quantities are employed in this method. Similarly, even though one-sided CPM calculations are possible, the errors near the end of the channel become substantial because the density becomes quite small, yielding large relative errors.

Stochastic simulations were also used to investigate the reliability of a new expression to calculate the ionic currents at different voltages,

Δ V

, given knowledge of the PMF, committor probabilities and density profiles at a reference voltage

Δ V_{ref}

. We found that the I-V dependence could be reconstructed over a range of ±100 mV, with respect to the reference voltage. Judicious choice of

Δ V_{ref}

can markedly improve the accuracy of the reconstruction. Specifically, the I-V reconstruction for p7 is much better for

Δ V_{ref}

- 35

mV than for

Δ V_{ref}

= 140 mV. Although much of the error can be attributed to the large voltage ramp for voltages away from

Δ V_{ref}

(3.2 kcal/mol at 140 mV), some of the error is due to the poor statistics in the direction against the field. This is also evident in the reconstruction of the I-V curve for GLIC, for which some simulations yielded no crossing events against the field.

Common goals of simulations of ion channels are to obtain the free energy profiles of ions translocating the channel and to determine electrophysiological properties of the channel. In some instances, a reliable estimate of the numbers of crossing events, from which the ionic currents can be calculated, is difficult to obtain from MD even for long simulation times. We have shown that the new methods perform very well both to obtain reliably the free energy profile across the channel and to allow for accurate determination of the I-V curves. In the latter case, it is desirable to use a reference voltage that yields good crossing statistics in both directions rather than a voltage that maximizes the total number of crossing events. In summary, if transport of ions through a channel can be satisfactorily described by the ED model, the new methods offer substantial reductions of computational effort without sacrificing accuracy. Our approach is amenable to extensions in which the advantages of MD and stochastic simulations are further combined on reliable theoretical grounds.

Supplementary Materials

The following are available at https://www.mdpi.com/article/10.3390/e23050571/s1, Figure S1: Pictures from MD simulations of TTX, p7, and GLIC, Figure S2: Total displacement charge calculated from ion crossing statistics and displacement current.

Author Contributions

Conceptualization, A.P.; methodology, A.P. and M.A.W.; software, M.A.W.; validation, M.A.W. and A.P.; formal analysis, A.P. and M.A.W.; writing–original draft preparation, A.P.; writing–review and editing, A.P and M.A.W.; visualization, M.A.W.; project administration, A.P.; funding acquisition, A.P. All authors have read and agreed to the published version of the manuscript.

Funding

Support for this research was provided by NASA’s Planetary Science Division Research Program.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

BD	Brownian Dyanmics
CPM	Committor Probability Method
CWDM	Current-Weighted Density Method
ED	electrodiffusion
GHK	Goldman Hodgkin Katz (equation)
IEEM	Integrated Electrodiffusion Equation Method
MD	molecular dynamics
PNP	Poisson-Nernst-Planck
PMF	potential of mean force
TTX	trichotoxin channel
WHAM	Weighted Histogram Analysis Method

References

Hille, B. Ion Channels of Excitable Membranes, 3rd ed.; Sinauer Associates, Inc.: Sunderland, MA, USA, 2001. [Google Scholar]
Zheng, J.; Trudeau, M.C. Handbook of Ion Channels; CRC Press: Boca Raton, FL, USA, 2015. [Google Scholar]
Kubalski, A.; Martinac, B. Bacterial Ion Channels and Their Eukaryotic Homologs; ASM Press: Washington, DC, USA, 2005. [Google Scholar]
Bocquet, N.; Nury, H.; Baaden, M.; Le Poupon, C.; Changeux, J.P.; Delarue, M.; Corringer, P.J. X-ray structure of a pentameric ligand-gated ion channel in an apparently open conformation. Nature 2009, 457, 111–114. [Google Scholar] [CrossRef] [PubMed]
Sakmann, B.; Boheim, G. Alamethicin-induced single channel conductance fluctuations in biological membranes. Nature 1979, 282, 336–339. [Google Scholar] [CrossRef] [PubMed]
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. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Martinac, B. Mechanosensitive ion channels: An evolutionary and scientific tour de force in mechanobiology. Channels 2012, 211–213. [Google Scholar] [CrossRef] [Green Version]
Bass, R.B.; Strop, P.; Barclay, M.; Rees, D.C. Crystal structure of Escherichia coli MscS, a voltage-modulated and mechanosensitive channel. Science 2002, 298, 1582–1587. [Google Scholar] [CrossRef] [Green Version]
Panaghie, G.; Purtell, K.; Tai, K.K.; Abbott, G.W. Voltage-dependent C-type inactivation in a constitutively open K+ channel. Biophys. J. 2008, 95, 2759–2778. [Google Scholar] [CrossRef] [Green Version]
Anstee, Q.M.; Knapp, S.; Maguire, E.P.; Hosie, A.M.; Thomas, P.; Mortensen, M.; Bhome, R.; Martinez, A.; Walker, S.E.; Dixon, C.I.; et al. Mutations in the Gabrb1 gene promote alcohol consumption through increased tonic inhibition. Nat. Commun. 2013, 4, 2816. [Google Scholar] [CrossRef] [Green Version]
Hou, X.; Outhwaite, I.R.; Pedi, L.; Long, S.B. Cryo-EM structure of the calcium release-activated calcium channel Orai in an open conformation. Elife 2020, 338, 1308–1313. [Google Scholar]
Ashcroft, F.M. Ion Channels and Disease; Academic Press: Cambridge, MA, USA, 1999. [Google Scholar]
Bagal, S.K.; Brown, A.D.; Cox, P.J.; Omoto, K.; Owen, R.M.; Pryde, D.C.; Sidders, B.; Skerratt, S.E.; Stevens, E.B.; Storer, R.I.; et al. Ion channels as therapeutic targets: A drug discovery perspective. J. Med. Chem. 2013, 56, 593–624. [Google Scholar] [CrossRef]
Garcia, M.L.; Kaczorowski, G.J. Ion channels find a pathway for therapeutic success. Proc. Natl. Acad. Sci. USA 2016, 113, 5472–5474. [Google Scholar] [CrossRef] [Green Version]
Aksimentiev, A.; Schulten, K. Imaging α-hemolysin with molecular dynamics: Ionic conductance, osmotic permeability, and the electrostatic potential map. Biophys. J. 2005, 88, 3745–3761. [Google Scholar] [CrossRef] [Green Version]
Sotomayor, M.; Vásquez, V.; Perozo, E.; Schulten, K. Ion conduction through MscS as determined by electrophysiology and simulation. Biophys. J. 2007, 92, 886–902. [Google Scholar] [CrossRef] [Green Version]
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, 1898–1906. [Google Scholar] [CrossRef] [Green Version]
Kutzner, C.; Grubmüller, 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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]
Zhu, F.; Hummer, G. Theory and simulation of ion conduction in the pentameric GLIC channel. J. Chem. Theory Comput. 2012, 8, 3759–3768. [Google Scholar] [CrossRef] [Green Version]
Stock, L.; Delemotte, L.; Carnevale, V.; Treptow, W.; Klein, M.L. Conduction in a biological sodium selective channel. J. Phys. Chem. B 2013, 117, 3782–3789. [Google Scholar] [CrossRef]
Jensen, M.Ø; Jogini, V.; Eastwood, M.P.; Shaw, D.E. Atomic-level simulation of current–voltage relationships in single-file ion channels. J. General Physiol. 2013, 141, 619–632. [Google Scholar] [CrossRef]
Ulmschneider, M.B.; Bagnéris, 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. USA 2013, 110, 6364–6369. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef]
Pohorille, A.; Wilson, M.A.; Wei, C. Validity of the electrodiffusion model for calculating conductance of simple ion channels. J. Phys. Chem. B 2017, 121, 3607–3619. [Google Scholar] [CrossRef] [PubMed]
Callahan, K.M.; Roux, B. Molecular dynamics of ion conduction through the selectivity filter of the NaVAb sodium channel. J. Phys. Chem. B 2018, 122, 10126–10142. [Google Scholar] [CrossRef] [PubMed]
Flood, E.; Boiteux, C.; Lev, B.; Vorobyov, I.; Allen, T.W. Atomistic Simulations of Membrane Ion Channel Conduction, Gating, and Modulation. Chem. Rev. 2019, 119, 7737–7832. [Google Scholar] [CrossRef] [PubMed]
Oestringer, B.P.; Bolivar, J.H.; Hensen, M.; Claridge, J.K.; Chipot, C.; Dehez, F.; Holzmann, N.; Zitzmann, N.; Schnell, J.R. Re-evaluating the p7 viroporin structure. Nature 2018, 562, E8–E18. [Google Scholar] [CrossRef] [PubMed]
Machtens, J.P.; Kortzak, D.; Lansche, C.; Leinenweber, A.; Kilian, P.; Begemann, B.; Zachariae, U.; Ewers, D.; de Groot, B.L.; Briones, R.; et al. Mechanisms of anion conduction by coupled glutamate transporters. Cell 2015, 160, 542–553. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Roux, B. The membrane potential and its representation by a constant electric field in computer simulations. Biophys. J. 2008, 95, 4205–4216. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Gumbart, J.; Khalili-Araghi, F.; Sotomayor, M.; Roux, B. Constant electric field simulations of the membrane potential illustrated with simple systems. Biochim. Biophys. Acta 2012, 1818, 294–302. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]
Pohorille, A.; Wilson, M.A. Computational Electrophysiology from a Single Molecular Dynamics Simulation and the Electrodiffusion Model. J. Phys. Chem. B 2021, 125, 3132–3144. [Google Scholar] [CrossRef]
Biró, I.; Pezeshki, S.; Weingart, H.; Winterhalter, M.; Kleinekathöfer, U. Comparing the temperature-dependent conductance of the two structurally similar E. coli porins OmpC and OmpF. Biophys. J. 2010, 98, 1830–1839. [Google Scholar] [CrossRef] [PubMed] [Green Version]
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. [Google Scholar] [CrossRef] [PubMed]
Zachariae, U.; Schneider, R.; Briones, R.; Gattin, Z.; Demers, J.P.; Giller, K.; Maier, E.; Zweckstetter, M.; Griesinger, C.; Becker, S.; et al. β-Barrel mobility underlies closure of the voltage-dependent anion channel. Structure 2012, 20, 1540–1549. [Google Scholar] [CrossRef] [Green Version]
Holzmann, N.; Chipot, C.; Penin, F.; Dehez, F. Assessing the physiological relevance of alternate architectures of the p7 protein of hepatitis C virus in different environments. Bioorgan. Med. Chem. 2016, 24, 4920–4927. [Google Scholar] [CrossRef] [PubMed]
Wood, M.L.; Freites, J.A.; Tombola, F.; Tobias, D.J. Atomistic modeling of ion conduction through the voltage-sensing domain of the Shaker K+ ion channel. J. Phys. Chem. B 2017, 121, 3804–3812. [Google Scholar] [CrossRef] [Green Version]
Chen, D.P.; Eisenberg, R.S. Flux, coupling, and selectivity in ionic channels of one conformation. Biophys. J. 1993, 65, 727–746. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef]
Coalson, R.D.; Kurnikova, M.G. Poisson–Nernst–Planck theory of ion permeation through biological channels. In Biological Membrane Ion Channels; Springer: Berlin, Germany, 2007; pp. 449–484. [Google Scholar]
Liu, J.L.; Eisenberg, B. Poisson-Nernst-Planck–Fermi theory for modeling biological ion channels. J. Chem. Phys. 2014, 141, 12B640_1. [Google Scholar] [CrossRef] [Green Version]
Liu, J.L.; Eisenberg, B. Molecular mean-field theory of ionic solutions: A Poisson-Nernst-Planck-Bikerman model. Entropy 2020, 22, 550. [Google Scholar] [CrossRef]
Im, W.; Roux, B. Brownian dynamics simulations of ions channels: A general treatment of electrostatic reaction fields for molecular pores of arbitrary geometry. J. Chem. Phys. 2001, 115, 4850–4861. [Google Scholar] [CrossRef]
Chung, S.H.; Kuyucak, S. Recent advances in ion channel research. Biochim. Biophys. Acta 2002, 1565, 267–286. [Google Scholar] [CrossRef] [Green Version]
Chung, S.H.; Allen, T.W.; Kuyucak, S. Conducting-state properties of the KcsA potassium channel from molecular and Brownian dynamics simulations. Biophys. J. 2002, 82, 628–645. [Google Scholar] [CrossRef] [Green Version]
Chung, S.H.; Krishnamurthy, V. Brownian Dynamics: Simulation for Ion Channel Permeation. In Biological Membrane Ion Channels; Springer: Berlin, Germany, 2007; pp. 507–543. [Google Scholar]
Berneche, S.; Roux, B. Molecular dynamics of the KcsA K+ channel in a bilayer membrane. Biophys. J. 2000, 78, 2900–2917. [Google Scholar] [CrossRef] [Green Version]
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. Science 2014, 346, 352–355. [Google Scholar] [CrossRef] [Green Version]
Allen, T.W.; Andersen, O.S.; Roux, B. Energetics of ion conduction through the gramicidin channel. Proc. Natl. Acad. Sci. USA 2004, 101, 117–122. [Google Scholar] [CrossRef] [Green Version]
Keener, J.P.; Sneyd, J. Mathematical Physiology; Springer: Berlin, Germany, 1998; Volume 1. [Google Scholar]
Berneche, S.; Roux, B. Energetics of ion conduction through the K+ channel. Nature 2001, 414, 73–77. [Google Scholar] [CrossRef]
Furini, S.; Domene, C. On conduction in a bacterial sodium channel. PLoS Comput. Biol. 2012, 8, e1002476. [Google Scholar] [CrossRef] [Green Version]
Chakrabarti, N.; Ing, C.; Payandeh, J.; Zheng, N.; Catterall, W.A.; Pomès, R. Catalysis of Na+ permeation in the bacterial sodium channel NaVAb. Proc. Natl. Acad. Sci. USA 2013, 110, 11331–11336. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Boiteux, C.; Vorobyov, I.; Allen, T.W. Ion conduction and conformational flexibility of a bacterial voltage-gated sodium channel. Proc. Natl. Acad. Sci. USA 2014, 111, 3454–3459. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Harpole, T.J.; Grosman, C. Side-chain conformation at the selectivity filter shapes the permeation free-energy landscape of an ion channel. Proc. Natl. Acad. Sci. USA 2014, 111, E3196–E3205. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Finol-Urdaneta, R.K.; Wang, Y.; Al-Sabi, A.; Zhao, C.; Noskov, S.Y.; French, R.J. Sodium channel selectivity and conduction: Prokaryotes have devised their own molecular strategy. J. General Physiol. 2014, 143, 157–171. [Google Scholar] [CrossRef] [Green Version]
Trick, J.L.; Chelvaniththilan, S.; Klesse, G.; Aryal, P.; Wallace, E.J.; Tucker, S.J.; Sansom, M.S. Functional annotation of ion channel structures by molecular simulation. Structure 2016, 24, 2207–2216. [Google Scholar] [CrossRef] [Green Version]
Flood, E.; Boiteux, C.; Allen, T.W. Selective ion permeation involves complexation with carboxylates and lysine in a model human sodium channel. PLoS Comput. Biol. 2018, 14, e1006398. [Google Scholar] [CrossRef] [Green Version]
Cottone, G.; Chiodo, L.; Maragliano, L. Thermodynamics and kinetics of ion permeation in wild-type and mutated open active conformation of the human α7 nicotinic receptor. J. Chem. Inf. Model. 2020, 60, 5045–5056. [Google Scholar] [CrossRef]
Jensen, M.; Borhani, D.W.; Lindorff-Larsen, K.; Maragakis, P.; Jogini, V.; Eastwood, M.P.; Dror, R.O.; Shaw, D.E. Principles of conduction and hydrophobic gating in K+ channels. Proc. Natl. Acad. Sci. USA 2010, 107, 5833–5838. [Google Scholar] [CrossRef] [Green Version]
Corry, B.; Thomas, M. Mechanism of ion permeation and selectivity in a voltage gated sodium channel. J. Am. Chem. Soc. 2012, 134, 1840–1846. [Google Scholar] [CrossRef]
Modi, N.; Bárcena-Uribarri, I.; Bains, M.; Benz, R.; Hancock, R.E.; Kleinekathöfer, U. Role of the central arginine R133 toward the ion selectivity of the phosphate specific channel OprP: Effects of charge and solvation. Biochemistry 2013, 52, 5522–5532. [Google Scholar] [CrossRef]
Jorgensen, C.; Furini, S.; Domene, C. Energetics of ion permeation in an open-activated TRPV1 channel. Biophys. J. 2016, 111, 1214–1222. [Google Scholar] [CrossRef] [Green Version]
Polovinkin, L.; Hassaine, G.; Perot, J.; Neumann, E.; Jensen, A.A.; Lefebvre, S.N.; Corringer, P.J.; Neyton, J.; Chipot, C.; Dehez, F.; et al. Conformational transitions of the serotonin 5-HT 3 receptor. Nature 2018, 563, 275–279. [Google Scholar] [CrossRef]
Klesse, G.; Rao, S.; Tucker, S.J.; Sansom, M.S. Induced polarization in molecular dynamics simulations of the 5-HT3 receptor channel. J. Am. Chem. Soc. 2020, 142, 9415–9427. [Google Scholar] [CrossRef]
Rao, S.; Klesse, G.; Lynch, C.I.; Tucker, S.J.; Sansom, M.S. Molecular Simulations of Hydrophobic Gating of Pentameric Ligand Gated Ion Channels: Insights into Water and Ions. J. Physic. Chem. B 2021, 125, 981–994. [Google Scholar] [CrossRef]
OuYang, B.; Xie, S.; Berardi, M.J.; Zhao, X.; Dev, J.; Yu, W.; Sun, B.; Chou, J.J. Unusual architecture of the p7 channel from hepatitis C virus. Nature 2013, 498, 521–525. [Google Scholar] [CrossRef] [PubMed]
Sauguet, L.; Poitevin, F.; Murail, S.; Van Renterghem, C.; Moraga-Cid, G.; Malherbe, L.; Thompson, A.W.; Koehl, P.; Corringer, P.J.; Baaden, M.; et al. Structural basis for ion permeation mechanism in pentameric ligand-gated ion channels. EMBO J. 2013, 32, 728–741. [Google Scholar] [CrossRef] [Green Version]
Chipot, C.; Pohorille, A. Free Energy Calculations. Theory and Applications to Chemistry and Biology; Springer: Berlin, Germany, 2007. [Google Scholar]
Pohorille, A. Free Energy Calculation for Understanding Membrane Receptors. In Computational Biophysics of Membrane Proteins (Theoretical and Computational Chemistry, Band 10); The Royal Society of Chemistry: London, UK, 2017; pp. 59–106. [Google Scholar]
Furini, S.; Domene, C. Computational studies of transport in ion channels using metadynamics. Biochim. Biophys. Acta 2016, 1858, 1733–1740. [Google Scholar] [CrossRef] [PubMed]
Marrink, S.J.; Berendsen, H.J.C. Simulation of water transport through a lipid membrane. J. Phys. Chem. 1994, 98, 4155–4168. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef]
Peter, C.; Hummer, G. Ion transport through membrane-spanning nanopores studied by molecular dynamics simulations and continuum electrostatics calculations. Biophys. J. 2005, 89, 2222–2234. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Comer, J.; Chipot, C.; González-Nilo, F.D. Calculating position-dependent diffusivity in biased molecular dynamics simulations. J. Chem. Theory Comput. 2013, 9, 876–882. [Google Scholar] [CrossRef]
Hummer, G. From transition paths to transition states and rate coefficients. J. Chem. Phys. 2004, 120, 516–523. [Google Scholar] [CrossRef]
Weinan, E.; Vanden-Eijnden, E. Transition-path theory and path-finding algorithms for the study of rare events. Ann. Rev. Phys. Chem. 2010, 61, 391–420. [Google Scholar]
Elber, R.; Bello-Rivas, J.M.; Ma, P.; Cardenas, A.E.; Fathizadeh, A. Calculating iso-committor surfaces as optimal reaction coordinates with milestoning. Entropy 2017, 19, 219. [Google Scholar] [CrossRef]
Berezhkovskii, A.M.; Szabo, A. Committors, first-passage times, fluxes, Markov states, milestones, and all that. J. Chem. Phys. 2019, 150, 054106. [Google Scholar] [CrossRef]
Kumar, S.; Rosenberg, J.M.; Bouzida, D.; Swendsen, R.H.; Kollman, P.A. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011–1021. [Google Scholar] [CrossRef]
Pieprzyk, S.; Heyes, D.M.; Brańka, A.C. Spatially dependent diffusion coefficient as a model for pH sensitive microgel particles in microchannels. Biomicrofluidics 2016, 10, 054118. [Google Scholar] [CrossRef] [Green Version]
Huber, G.A.; McCammon, J.A. Brownian Dynamics Simulations of Biological Molecules. Trends Chem. 2019, 1, 727–738. [Google Scholar] [CrossRef]
Pratt, L.R. A statistical method for identifying transition states in high dimensional problems. J. Chem. Phys. 1986, 85, 5045–5048. [Google Scholar] [CrossRef]
Gennis, R.B. Biomembranes: Molecular Structure and Function; Springer: New York, NY, USA, 1989. [Google Scholar]
Parsegian, V.A. Energy of an Ion crossing a Low Dielectric Membrane: Solutions to Four Relevant Electrostatic Problems. Nature 1969, 221, 844–846. [Google Scholar] [CrossRef]
Neumcke, B.; Läuger, P. Nonlinear electrical effects in lipid bilayer membranes. II. Biophys. J. 1969, 9, 1160–1170. [Google Scholar] [CrossRef] [Green Version]
Montserret, R.; Saint, N.; Vanbelle, C.; Salvay, A.G.; Simorre, J.P.; Ebel, C.; Sapay, N.; Renisio, J.G.; Bockmann, A.; Steinmann, E.; et al. NMR structure and ion channel activity of the p7 protein from hepatitis C virus. J. Biol. Chem. 2010, 285, 31446–31461. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]

Figure 1. (a) Committor probabilities for Cl

^{-}

in p7 at

- 140

mV (red),

- 70

mV (green),

- 35

mV (black), 0 V (blue), 70 mV (cyan) and 140 mV (magenta). Error bars are shown for the N6 data sets at

- 35

mV and 140 mV; (b) Committor probabilities for p7 at 140 mV from the N6 data set for 1-sided forward (green) and backward (blue) trajectories, respectively, 2-sided data set in the backward direction (red lines), and average in the forward direction with error bars (red symbols). In the inset we show the number of first passage trajectories to reach z for one N6 data set in the forward (green) and backward (magenta) directions and the total (light blue).

Figure 1. (a) Committor probabilities for Cl

^{-}

in p7 at

- 140

mV (red),

- 70

mV (green),

- 35

mV (black), 0 V (blue), 70 mV (cyan) and 140 mV (magenta). Error bars are shown for the N6 data sets at

- 35

Figure 2. (a) PMFs for K

^{+}

(lower curves) and Cl

^{-}

(upper curves) in TTX from stochastic simulations with an applied voltage of 50 mV. The PMFs have been reconstructed by way of CWDM at the N6 (blue) and N7 (gold) levels or by way of CPM at the N6 (green) and N7 (magenta) levels; (b) PMF for Na

^{+}

in GLIC from stochastic simulations with applied voltage of 100 mV. The PMF has been reconstructed by way of CWDM at the N7 (blue) and N8 (gold) level or by way of CPM at the N7 (green) and N8 (magenta) level. In both panels, the underlying PMF is in red.

Figure 2. (a) PMFs for K

^{+}

(lower curves) and Cl

^{-}

^{+}

Figure 3. (a) PMF for Cl

^{-}

in p7 from stochastic simulations with an applied voltage of 140 mV. The PMFs have been reconstructed by way of CWDM at the N6 (blue) and N7 (gold) levels or by way of CPM at the N6 (green) and N7 (magenta) levels. The input PMF (red) is shown for reference. PMFs at the N8 level are not shown, as they coincide with the underlying PMFs and statistical errors associated with this level arequite small and are poorly visible at this scale; (b) PMFs for P7 reconstructed by way of one-sided forward trajectories (green) using Equation (13) and backward trajectories (blue) using Equation (14) from stochastic simulations at the N6 level with applied voltage of 140 mV. Two-sided reconstruction (magenta) and the underlying PMF (red) are shown for comparison. Note that one-sided, but not two-sided reconstructions are burdened with large errors at the ends.

Figure 3. (a) PMF for Cl

^{-}

Figure 4. (a) I-V curves for K

^{+}

(green) and Cl

^{-}

(blue) in TTX reconstructed from simulations at 50 mV at the N6 level. Blue and green dots are currents obtained from direct simulations at specific voltages.; (b) I-V curves for Cl

^{-}

in p7 reconstructed from simulations at 140 mV at the N6 (blue), N7 (green) and N8 (red) level, and for

- 35

mV at the N6 level (magenta). N7 and N8 curves are not shown because they are almost identical to the N6 results. Black dots are currents obtained from direct simulations at specific voltages. All reconstructions were done using the PMFs obtained by way of CPM. The results of reconstructions using the PMFs from CWDM are not displayed because they are nearly identical.

Figure 4. (a) I-V curves for K

^{+}

(green) and Cl

^{-}

(blue) in TTX reconstructed from simulations at 50 mV at the N6 level. Blue and green dots are currents obtained from direct simulations at specific voltages.; (b) I-V curves for Cl

^{-}

in p7 reconstructed from simulations at 140 mV at the N6 (blue), N7 (green) and N8 (red) level, and for

- 35

Figure 5. I-V curves for Na

^{+}

in GLIC reconstructed from simulations at 100 mV at the N7 level with PMF from CPM (blue), at the N7 level with PMF from CWDM (magenta), and N8 with PMF from CWDM (red). N8 with CPM (not shown) is almost identical to N7 CPM. Black dots are currents obtained from direct simulations at specific voltages.

Figure 5. I-V curves for Na

^{+}

Figure 6. Reconstructions of I-V curves in TTX from individual sets of trajectories for K

^{+}

(a) and Cl

^{-}

(b). The PMFs were obtained from CPM (upper panels) or CWDM (lower panels). The curves were calculated by way of Equation (20) (blue) or Equation (18) (green). All reconstructions were carried out from simulations at applied voltage of 50 mV at the N6 level. Note that blue curves, but not green curves, are tightly clustered together indicating that Equation (20) is more accurate than Equation (18).

Figure 6. Reconstructions of I-V curves in TTX from individual sets of trajectories for K

^{+}

(a) and Cl

^{-}

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Wilson, M.A.; Pohorille, A. Electrophysiological Properties from Computations at a Single Voltage: Testing Theory with Stochastic Simulations. Entropy 2021, 23, 571. https://doi.org/10.3390/e23050571

AMA Style

Wilson MA, Pohorille A. Electrophysiological Properties from Computations at a Single Voltage: Testing Theory with Stochastic Simulations. Entropy. 2021; 23(5):571. https://doi.org/10.3390/e23050571

Chicago/Turabian Style

Wilson, Michael A., and Andrew Pohorille. 2021. "Electrophysiological Properties from Computations at a Single Voltage: Testing Theory with Stochastic Simulations" Entropy 23, no. 5: 571. https://doi.org/10.3390/e23050571

APA Style

Wilson, M. A., & Pohorille, A. (2021). Electrophysiological Properties from Computations at a Single Voltage: Testing Theory with Stochastic Simulations. Entropy, 23(5), 571. https://doi.org/10.3390/e23050571

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Electrophysiological Properties from Computations at a Single Voltage: Testing Theory with Stochastic Simulations

Abstract

1. Introduction

2. Theory and Method

2.1. Calculating the Potential of Mean Force

2.2. Calculating I-V Dependence from Simulation at a Single Voltage

2.3. Stochastic Simulations

3. Results and Discussion

3.1. Connection with Molecular Dynamics

3.2. Committor Probabilities

3.3. The Potential of Mean Force

3.4. Current-Voltage Dependence

3.5. Reversal Potential

4. Conclusions

Supplementary Materials

Author Contributions

Funding

Conflicts of Interest

Abbreviations

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI