Unifying Single-Channel Permeability From Rare-Event Sampling and Steady-State Flux

Various all-atom molecular dynamics (MD) simulation methods have been developed to compute free energies and crossing rates of ions and small molecules through ion channels. However, a systemic comparison across different methods is scarce. Using a carbon nanotube as a model of small conductance ion channel, we computed the single-channel permeability for potassium ion using umbrella sampling, Markovian milestoning, and steady-state flux under applied voltage. We show that a slightly modified inhomogeneous solubility-diffusion equation yields a single-channel permeability consistent with the mean first passage time (MFPT) based method. For milestoning, applying cylindrical and spherical bulk boundary conditions yield consistent MFPT if factoring in the effective bulk concentration. The sensitivity of the MFPT to the output frequency of collective variables is highlighted using the convergence and symmetricity of the inward and outward MFPT profiles. The consistent transport kinetic results from all three methods demonstrated the robustness of MD-based methods in computing ion channel permeation. The advantages and disadvantages of each technique are discussed, focusing on the future applications of milestoning in more complex systems.

Introduction

Ion channels are complex biological nanopores that perform vital physiological functions with high sensitivity and precision. Over the decades, molecular dynamics (MD) simulation has become an indispensable tool for computing the functional properties of ion channels directly from their dynamic structures. Various MD-based methods were developed for investigating the thermodynamics and kinetics of ions or small-molecules permeation at the single-channel level. Many pioneering atomistic MD simulations on ion channels have focused on computing ion permeation from equilibrium free energy profiles, or potential of mean force (PMF), in conjunction with electro-diffusion theory (Bernèche and Roux, 2001; Allen et al., 2004; Domene et al., 2008; Fowler et al., 2013). The increased computing power and performance of MD engines have also enabled researchers to simulate single-channel conduction explicitly under a constant external electric field (Khalili-Araghi et al., 2006; Roux, 2008) or an asymmetric ionic concentration across the channel (Kutzner et al., 2011; Khalili-Araghi et al., 2013). If the system reaches a steady-state under voltage or concentration gradient, a mean flux rate and a steady-state density profile can be obtained from the ensemble of nonequilibrium processes. These equilibrium and nonequilibrium MD simulations have significantly deepened our understanding of the ion channel permeation process at the high temporal and spatial resolution (Roux, 1998; Zheng and Trudeau, 2015; Flood et al., 2019; Carnevale et al., 2021).

Unlike the steady-state flux under voltage, the equilibrium MD approaches can be generally applied to any small-molecule permeation (neutral or charged). Several theoretical frameworks can be used to compute crossing rates from PMF profiles obtained from enhanced sampling simulations. Particulary, if the PMF is dominated by a single large barrier and the permeant diffusion is constant at the barrier region, the crossing rate can be estimated via Kramer’s theory or transition state theory (TST) borrowed from reaction kinetics. However, for complex biological ion channels, the aforementioned assumptions may be far from satisfied. Alternatively, molecule permeation may be considered a one-dimensional nonreactive diffusive process that can be described using the fluctuation-dissipation theorem. For instance, PMF can be used together with the position-dependent diffusion coefficient to estimate permeability using the inhomogeneous solubility-diffusion (ISD) equation (Diamond and Katz, 1974). ISD equation has been applied successfully in studying solute permeation across the membrane (Awoonor-Williams and Rowley, 2016; Venable et al., 2019). Herein, we show that a slightly modified ISD equation yields a single-channel permeability consistent with a mean first passage time (MFPT) based method, which extracts detailed kinetics along the molecular permeation pathway directly from rare-event sampling methods. Recent examples of such rare-event sampling applied on ion channels include milestoning (Alberini et al., 2018; Cottone et al., 2020; Jiang et al., 2021a), weighted ensemble sampling (Adelman and Grabe, 2015), and Markov state models (Teo and Schulten, 2013; Choudhary et al., 2014; Domene et al., 2021; Hempel et al., 2021). In theory, the PMF-based method, MFPT-based method, and steady-state flux under voltage should converge to the same single-channel permeability for the same studied system. However, a systemic comparison between different methods is still lacking.

In this work, we use a carbon nanotube (CNT) as a model (Figure 1A) of small conductance (∼2 pS) ion channel to compare K⁺ permeability from milestoning, umbrella sampling (US), and voltage simulations. This CNT system has been used to compute K⁺ permeability using a transition path approach similar to the reactive flux method (Zhou and Zhu, 2019). We chose this system because its free energy barrier height (∼4 kcal/mol) and microsecond-timescale crossing rate are physiologically relevant. Such a system requires nontrivial sampling (beyond the capability of brute-force MD), but the rigidity of the CNT still allows good convergence and unambiguous comparison of all methods tested here.

FIGURE 1

FIGURE 1. (A) Simulated CNT system, consisting of two fixed layers of carbon atoms as water impermeable membrane. 6 K⁺ and 6 Cl⁻ ions are shown in tan and pink, CNT in orange color, and water oxygen in cyan. Membrane carbon atoms are shown as blue spheres (B). z-coordinate distribution of a tagged K⁺ from each umbrella sampling window (corresponding data from milestoning are shown in Figure 2B). (C) PMFs from two blocks of umbrella sampling trajectories, generated from MBAR with 80,000 data points per window, with error bars in shaded color (D). The position-dependent diffusion constant of K⁺ computed from umbrella sampling. The transparent lines represent data from the first 20 ns and second 20 ns per window. A thick green line represents averaged and symmetrized values. Inset is the plot of the correlation function used to compute correlation time $i$ and D(z) (see methods) (E). Local resistance (the integrand of Eq. 3) for permeating K⁺ (F). Integration of the permeation resistance, 1/ $P$, as a function of z.

The original milestoning simulation requires running short trajectories in each milestone until they reach another milestone (Faradjian and Elber, 2004). Here, we use the “soft-walls” Voronoi-tessellated Markovian milestoning, which confines the sampling within the Voronoi cells using flat-bottom harmonic restraining potentials (Maragliano et al., 2009). The implementation of this “soft-walls” version (referred to as milestoning thereafter) resembles, to a large degree, the conventional umbrella sampling setup. A detailed comparison of the sampling, PMF, and MFPT results from milestoning and umbrella sampling is the focus of this study. In addition, we tested two bulk boundary conditions, namely the cylindrical and spherical boundaries, that are particularly useful for conducting milestoning on ion channels.

The CNT system chosen here is designed to satisfy the symmetric single barrier requirement, thus allowing us to check the robustness of the milestoning method by computing both inward and outward permeation rates. We also show that the MFPT from milestoning is extremely sensitive to the frequency of recording the relevant collective variables (e.g., the coordinates of the tagged ion). All physical quantities and the obtained results are summarized in Table 1. The overall consistent single-channel permeability demonstrated the robustness of the theoretical and computational framework tested here. The limitation and strengths of each method are discussed and compared.

TABLE 1

TABLE 1. Summary of the three methods for computing single-channel permeability.

Theory and Methods

Relation Between Single-Channel Permeability, Mean First Passage Time, and Conductance

Assuming permeating molecules do not interact under sufficient low concentration, under physiological conditions, single-channel permeability $P$ (cm³/s) can be related linearly to the rate of crossing $k$ or mean first passage time (MFPT) $<t>$ under equilibrium, $P=k/c=1/c<t>$ , in which $c$ is the symmetric solute concentrations. Here, we use the number of molecules per second for $k$, seconds per molecule for MFPT, and molecule/cm³ for $c$.

For ionic permeation under voltage and/or concentration gradient, Goldman–Hodgkin–Katz (GHK) flux equation describes the ionic flux across a homogenous membrane as a function of a constant electric field (voltage) and an ionic concentrations gradient. Under symmetric concentration and constant voltage, the current ($I$) and the permeability ($P)$ can be related by the GHK flux equation, $I=P{z}^2\frac{V_m{F}^2}{RT}c,$ where $z$ is the charge of the permeant, $V_m$ is the voltage, F is the Faraday constant, R is the gas constant, T is the absolute temperature, and $c$ is the concentration (Hille, 2001). When applied to a membrane-embedded single-channel model and ions only cross the membrane through the channel, $P$ in the GHK equation corresponds to the single-channel permeability. GHK flux equation thus relates single-channel conductance $\gamma$, a nonequilibrium property, to the equilibrium property $P$.

It can be seen from the equations above that the crossing rate $k$, MFPT, and conductance $\gamma$ are all concentration-dependent; only $P$ is independent of concentration. Experimentally, $P$ is usually measured relative to the potassium ion permeability, thus representing an intrinsic property of each channel. It is hence an ideal quantity for rigorous comparison between different computational methods.

System Setup and Equilibrium Protocols

The coordinates of the carbon tube (CNT) were taken from Zhou and Zhu (2019). Briefly, it is an uncapped armchair CNT with 13.5 Å in length and 5.4 Å in radius. Two carbon sheets form an artificial membrane to separate the solution (Figure 1A). Constraints were applied to all carbon atoms to keep the system rigid. The CHARMM36 force field was used (MacKerell et al., 1998; Mackerell et al., 2004). After solvation, the box size was 38 × 38 × 75 Å³, which contained 2503 TIP3P water molecules, six K⁺, and six Cl⁻. All MD simulations were performed using 1 fs time step using NAMD2.13 package under NVT ensemble, with 1 atm and 300 K temperature using Langevin thermostat (Hoover et al., 1982; Evans, 1983). Cutoff for calculating vdW interaction and short-range electrostatic interaction was set at 12 Å and force-switched at 10 Å. Long-range electrostatic interactions were calculated using the particle mesh Ewald algorithm (Darden et al., 1993). The system was equilibrium for 100 ns before conducting umbrella sampling, milestoning, and voltage simulations.

Umbrella Sampling Simulations

A total of 28 windows (−27Å < z < +27 Å) were sampled. Each window was separated by 2 Å apart. The tagged K⁺ was restraint by a harmonic restraint on z and a flat-bottom harmonic cylindrical restraint. The force constant for harmonic restraint was 2.5 kcal/mol/Å, and that for cylindrical restraint was 10 kcal/mol/Å within 6 Å on the X-Y radius plane. The reference dummy atom to pull the K⁺ was set at (0, 0, 0). The harmonic distance restraint was determined by the projected vector along z between the dummy atom and tagged K⁺. The cylindrical restraint was determined by the center of mass of all carbon atoms from CNT. Each window was run for 40 ns (Figure 1B). The PMF (Figure 1C) was computed using Pymbar 3.0.3 (Chodera et al., 2007; Shirts and Chodera, 2008). The output frequency was 0.5 ps per frame.

Position-Dependent Diffusion Coefficient D(z)

D(z) of K⁺ inside the CNT was calculated from umbrella sampling trajectories (Figure 1D). The correlation time was extracted from each umbrella window $i$ using ${\tau }_i=\underset{0}{\overset{\infty }{{\displaystyle \int }}}\frac{\langle \delta z\left(t\right)\delta z\left(0\right){\rangle }_i}{\langle \delta z{\rangle }_i^2}\text{d}t$ , where $\delta z\left(t\right)=z\left(t\right)-\langle z\rangle$ is the deviation of the z-position of the ion at time $t$, $z\left(t\right)$, from the time-averaged position $z$ in each window. $\langle \delta z^2{\rangle }_i=\langle z^2{\rangle }_i-\langle z{\rangle }_i^2$ is the variance. Following the formulation of Berne et al. (1988), Woolf and Roux (1994), Hummer (2005), the Laplace transformation of the velocity autocorrelation function along the reaction coordinate $z$ in the harmonically restrained umbrella sampling gives $D\left(z_i\right)=\frac{\langle \delta z^2\rangle }{{\tau }_i}$ .

Markovian Milestoning With Cylindrical Bulk Restraint

The same as in umbrella sampling, the z-coordinate of the tagged ion was used to define a set of Voronoi cells along the channel pore and identify the milestones as the boundaries between the cells. To facilitate the comparison in analysis, we kept the Voronoi cell setup (28 cells and 2 Å apart) and the bulk cylindrical restraint (6 Å radius) identical to our umbrella sampling windows (Figure 2A). The only difference is that a flat-bottom harmonic restraint of force constant 100 kcal/mol/Å, instead of the weak harmonic restraint, was used to confine the sampling within each cell (Figure 2B vs. Figure 1B). We then ran 28 local simulations confined in each cell and collected the kinetics of transitions between milestones. More specifically, let us introduce a set of M Voronoi cells B_i, i = 1, … ,M. Since the total flux in and out of each cell is zero at statistical equilibrium, the rate of attempted escape from cells B_j to B_i, $k_{j\to i}$ , and the equilibrium probability π_i for the tagged ion to be in cell B_i satisfies a balance equation:

${\displaystyle \sum }_{j=1,j\ne i}{\pi }_j{k}_{j\to i}={\displaystyle \sum }_{j=1,j\ne i}{\pi }_i{k}_{i\to j} , {\displaystyle \sum }_{i=1}{\pi }_i=1(1)$

FIGURE 2

FIGURE 2. (A). Raw data plotted along the channel z-axis and radial distance $R=\sqrt{x^2+y^2}$ from channel center axis (x, y = 0, 0) from 28 milestoning sampling cells. (B) Distribution of the milestoning data along the z-axis (same plot from the US is shown in Figure 1A). (C) PMF from Milestoning sampling. Different color represents different Colvars frequency (i.e., the frequency of recording the z-coordinates of the tagged ion). The bold purple (0.5 ps) is the data used for the final comparison. (D) MFPT plot with the same color representation as PMF plot. The curves with dot solid lines are outward directions of tagged ion, and dash curves present inward direction. (E) Maximum waiting time between successive transitions (F). Minimum waiting time between successive transitions (G). z-position decorrelation time in each milestoning cell. (H) Convergence of the variables in Eq. 2 for computing the rate matrix. Milestoning cell 11 is chosen as an example here.

The free energy of each cell can be obtained from the solution of Eq. 1 as −k_BTln(π_i) (Figure 2C). By defining a milestone S_ij as the boundary between two adjacent Voronoi cells B_i and B_j, the dynamics of the system is reduced to that of a Markov chain in the state space of the milestone indices (Vanden-Eijnden and Venturoli, 2009). The MFPT between any pair of milestones S_ij and S_ik can hence be calculated from the rate matrix whose elements $q_{ij,ik}$, the rate of moving from milestone S_ij to S_ik, are given by

$q_{ij,ik}=\frac{{\pi }_i{n}_{ij,ik}^i}{{\pi }_i{r}_{ij}^i+{\pi }_j{r}_{i\dot{J}}^j}(2)$

where $n_{ij,ik}^i$ is the number of transitions from S_ij to S_ik, normalized by the time spent in cell B_i . $r_{ij}^i$ is the time passed in cell B_i after having hit S_ij before hitting any other milestone, normalized by the total time spent in cell B_i . The inward and outward MFPT profiles were obtained by reversing the milestone indices when constructing the rate matrix (Figure 2D). The $n_{ij,ik}^i$ and $r_{ij}^i$ can be used to monitor the convergence of the rate matrix (Figure 2H).

The total sampling time of all 28 cells was 4.9 μs, where each cell was sampled between 150 and 300 ns. The NAMD Colvars output frequency was 0.5 ps. The PMF and MFPT were computed using a set of in-house python scripts https://github.com/yichunlin79/CNT_milestoning_method with different frame sizes. In order to check whether the Colvars output frequency has any effect on MFPT, additional milestoning simulations were conducted with the Colvars output frequency of 0.2 ps and a total sampling time of 2.74 μs.

Markovian Milestoning With Spherical Bulk Restraint

For spherical bulk restrained milestoning, a total of 14 Voronoi cells were used, including eight cells inside the channel (identical to the milestoning above) and three layers of spherical shell on each side of the channel (Figure 3A). The distance between the tagged K⁺ ion and two dummy atoms fixed at the Cartesian coordinates of (0, 0, −8) and (0, 0, 8) are used to set up the spherical shells in bulk with radius increments of 3, 3, and 4 Å. Additional z > |8|Å restraint is applied to keep the ion outside the channel. All restraint force constant is 100 kcal/mol/Å. The length of each bulk window is 150 ns with Colvars output frequency of 0.5 ps⁻¹.

FIGURE 3

FIGURE 3. (A) Sampling plot with spherical restraint at two bulk ends. The three spherical radius intervals are 3 Å, 3 Å, and 4 Å from channel to bulk (B). PMF of the spherical restraint system. The x-axis represents the cell index number. (C) MFPT of the spherical restraint system. x-axis represents the cell index numbers.

Voltage Simulations

After 100 ns equilibrium simulation, constant electric fields corresponding to the transmembrane potential of ±0.3 and ±0.4 V were applied perpendicular to the membrane to all the atoms using NAMD2.13. In order to be consistent with the umbrella sampling and milestoning, a cylinder restraint of 6 Å radius was applied to a tagged K⁺ in bulk with 10 kcal/mol/Å force constant for ±0.4 V systems. All other ions moved freely beyond the cylinder restraint region. The K⁺ conductance was computed by counting the total number of crossing events and computing the charge displacement along the z-axis (Figure 4). Error bar was computed from three independent replicas of 200 ns. For ±0.3 V systems, all six K⁺ were restrained inside the same cylindrical bulk boundary, and a single replica of 60 ns was conducted. The time step was 1 fs, and the output frequency was 5 ps for all voltage systems.

FIGURE 4

FIGURE 4. Current-voltage (I–V) plot for a single K+ with a cylinder restraint of radius 6 Å in bulk. The slope of the linear fitting defines the conductance (A). K⁺ current computed using direct counting method. (B). K⁺ current computed using charge displacement method. The error bars are calculated from three replicas of 200 ns at each voltage.

Results and Discussion

Single-Channel Permeability From Inhomogeneous Solubility-Diffusion (ISD) Equation

Molecular permeation through ion channel can be described by ISD if the one-dimensional free energy and diffusion along channel normal is sufficient to describe the diffusion process (relaxation of orthogonal degrees of freedom is fast relative to the reaction coordinate) and the permeant velocity relaxation time is instantaneous (on the scale of integration time step). The only difference with the ISD equation used for membrane permeation is that a flat-bottom lateral potential $u(x,y)$ is often used to confine a single tagged ion in a cylindrical bulk region outside the ion channel. The effective cross-sectional area due to the lateral restraint is thus $\mathrm{\iint }{e}^{-\beta u\left(x,y\right)}dxdy$, which can be approximated to $\pi r^2$ in a homogenous bulk, where $r$ is the radius of the cylinder. Hence, $\pi r^2$ defines the effective bulk concentration in the simulated region, which led the probability of the ion inside the channel $p\left(z\right)$ over the true bulk density $\rho$ to be $\frac{p\left(z\right)}{\rho }= \pi r^2{\text{e}}^{-\beta w\left(z\right)}$ , where $w\left(z\right)$ is the PMF with the bulk value set to zero at the cylindrical region and $\beta =1/k_{B}T$ (Allen et al., 2004; Zhu and Hummer, 2012). $T$ is the temperature, and $k_B$ is Boltzmann’s constant. Therefore, at low ionic concentration, single-channel permeability can be estimated using a slightly modified ISD equation:

$P=\pi r^2{\left(\underset{{\mathrm{z}}_1}{\overset{{\mathrm{z}}_2}{{\displaystyle \int }}}\frac{e^{\beta \mathrm{w}\left(\mathrm{z}\right)}}{D\left(\mathrm{z}\right)}d\mathrm{z}\right)}^{-1}(3)$

In Eq. 3, $D\left(z\right)$ is the position-dependent diffusion constant of the studied permeant along the z-axis (Figure 1D). The interval of the integration, $\left[z_1,z_2\right]$, is the lower and upper boundaries of the channel pore, beyond which PMF reaches the bulk value. $r$ is the radius of the cylindrical restraint when the ion is outside of the $\left[z_1,z_2\right]$ interval. It is necessary to set $r$ larger than the maximum pore radius so that it has no energetic contribution inside the pore. The radius of the cylindrical restraint defines the effective bulk concentration. Thus, it offsets the bulk PMF value and ensures that the single-channel permeability from Eq. 3 is concentration-independent.

In both umbrella sampling and milestoning, the same cylindrical restraint with $r$ = 6 Å is applied in the bulk region, and the same window size of 2 Å was used. The only difference is that a weak harmonic restraint with a force constant of 2.5 kcal/mol is applied for all umbrella windows to ensure sufficient overlapping between neighboring samplings, but a strong flat-bottom harmonic restraint with a force constant of 100 kcal/mol is applied for all milestoning cells to confine the sampling within each cell. Figure 1B and Figure 2B illustrate the biased sampling distribution imposed by these two types of restraints. The PMFs from the US are shown in Figure 1C. With bulk value offset to zero, a broad energy barrier of 3.8 kcal/mol located inside the channel region is consistent with a previously reported PMF (Zhou and Zhu, 2019).

Using the PMF or $w\left(z\right)$ in Figure 1C and D(z) in Figure 1D, the permeability estimated from Eq. 3 is (8.96 ± 0.02) × 10⁻¹⁶ cm³/s. We can also plot local resistance (the integrand of Eq. 3) for permeating K⁺ through CNT (Figure 1E) and the integration of the permeation resistance, 1/ $\mathrm{P}\left(\mathrm{z}\right)$, as a function of the z-axis (Figure 1F). It is not surprising that the 1/ $P\left(z\right)$ bears the same feature as the outward MFPT in Figure 2D.

Permeability Computed From Mean First Passage Time

The MFPT of a single K⁺ crossing the CNT is computed from Voronoi-tessellated Markovian milestoning simulations (see Methods). The distributions of the tagged K⁺ confined in each 2 Å cell by flat-bottom harmonic restraint are shown in Figures 2A,B. Milestoning simulation yields a consistent PMF profile with the highest energy barrier of 4.1 kcal/mol at the center of the CNT (Figure 2C). As the CNT used here is symmetric by design, a rigorous check of sampling convergence is the perfect symmetricity (mirror image) of the inward and outward MFPT profiles (Figure 2D). The inward and outward MFPT profiles can be obtained by reversing the milestone indices when constructing the transition rate matrix.

We found that while PMF is nearly insensitive to the Colvars frequency (i.e., the frequency of recording the z-coordinates of the tagged ion) tested here, MFPT is susceptible to this frequency. In the current study, the frequency of 5 ps⁻¹ severely overestimates the MFPT due to the missing transition events. Lower frequency also yields less data, which leads to asymmetric MFPTs. Here, the MFPT from the sampling saved per 5 ps has ten times fewer data points than the one from 0.5 ps. Thus, it failed to converge even after 18 μs of sampling. The ideal frequency has to be system-dependent (local diffusion and shape of the underlying free energy landscape). For our CNT system, the Colvars frequencies of 0.2 ps⁻¹ and 0.5 ps⁻¹ yield an identical and symmetric MFPT of 2.6 ± 0.03 μs for K⁺ permeation.

Using the PMF, $w\left(z\right),$ and MFPT, $t,$ from the milestoning (with 0.5 ps⁻¹ Colvars frequency), the single-channel permeability computed from Eq. 4, derived from Eq. 3 and Eq. 5, (Votapka et al., 2016) is 2.92 × 10⁻¹⁶ cm³ s⁻¹, in fairly good agreement with the permeability of 8.96 × 10⁻¹⁶ cm³s⁻¹ computed from ISD equation (Eq. 3) using umbrella sampling data:

$P=\frac{\pi r^2{{\displaystyle \int }}_{\mathrm{z}1}^{\mathrm{z}2}{e}^{-\beta \mathrm{w}\left(\mathrm{z}\right)}\text{d}\mathrm{z}}{2\langle t\rangle }(4)$

Decorrelation Time Versus Waiting Time in Milestoning

Vanden-Eijnden et al. have shown that Markovian milestoning yields exact MFPTs if the milestones are chosen such that successive transitions between them are statistically independent (Vanden-Eijnden et al., 2008; Vanden-Eijnden and Venturoli, 2009) and thus no definition of lag time is needed. To check this assumption for transitions between two neighboring milestones, in each cell, the maximum and minimum waiting time between two neighboring milestones is extracted and plotted in Figures 2E,F. The mirror image-like relation between these two plots manifests that the transition down the slope of PMF is the fastest, and the one against the slope is the slowest. Hence, the longest waiting time of 94.0 ps and the shortest waiting time of 9.0 ps are in the same cell where the PMF is steepest. The velocity decorrelation time is less than the smallest frame size (0.2 ps). The z-position decorrelation time for the tagged K⁺ in each cell is plotted in Figure 2G. The maximum positional decorrelation time (7.6 ps) is also located at the steepest PMF region. Hence, for all pairs of milestones, the velocity decorrelation time and position decorrelation time are both less than the minimum waiting time for successive transitions between milestones.

Mean First Passage Time Computed From Spherical Boundary Condition

Laterally confining the ion in the bulk region (cylindrical restraint) is convenient for describing the thermodynamics and kinetics of the ions across the channel along the channel pore axis (z-axis). However, the geometries of ion channels are diverse. For funnel-shaped channel pores [e.g., connexin hemichannel (Jiang et al., 2021a)] or pores connected with lateral fenestration [e.g., Piezo1 channel (Jiang et al., 2021b)], a spherical boundary may be a better choice to capture the distribution and dynamics of ions near the channel entrance. Hence, we further tested milestoning simulation using spherical restraint for ions in the bulk region (Figure 3A). Unlike the cylindrical restraint, which yields constant ionic concentration along the z-axis, the effective ionic concentration in the current spherical bulk cells decreases as the radius of the sphere increases. Thus, the PMF and MFPT are plotted against the milestoning cell index rather than the z-axis (Figures 3B,C).

At low concentration, single-channel permeability (cm³/s) can be related linearly to the mean first passage time (MFPT) $<t>$ under equilibrium p = 1/c⟨t⟩, in which c is the symmetric solute concentrations. Because single-channel permeability is an intrinsic property of a channel, independent of solute concentration or the shape of the bulk cells, the ratio of MFPT from spherical restraint over cylindrical restraint should be equal to the reciprocal bulk concentration ratio. The concentration of a single ion in a hemisphere with a radius of 10 Å is 1.26 M and in a cylinder with a radius of 6 Å and length of 10 Å is 0.68 M. Therefore, the concentration ratio of ∼2 is indeed consistent with the MFPT of 2.6 μs from cylindrical restraint (Figure 2E) and 1.3 μs from spherical restraint (Figure 3C).

Mean First Passage Time Computed From Umbrella Sampling

In the high diffusion limit, the MFPT of diffusive motion of K⁺ from the lower to upper boundaries of the channel pore, $\left[z_1,z_2\right]$ , can be written as

$MFPT=\underset{\mathrm{Z}1}{\overset{\mathrm{Z}2}{{\displaystyle \int }}}\frac{e^{\beta \mathrm{w}\left(\mathrm{Z}\right)}}{D\left(\mathrm{z}\right)}\text{d}\mathrm{z} \underset{\mathrm{Z}1}{\overset{\mathrm{Z}}{{\displaystyle \int }}}{e}^{-\beta w\left({\mathrm{Z}}^{\prime }\right)}\text{d}{\mathrm{z}}^{\prime }(5)$

Eq. 5 was originally developed for computing the average reaction time for diffusion processes governed by a Smoluchowski-type diffusion equation (Szabo et al., 1980). It is also used to derive Eq. 4 from Eq. 3 (Votapka et al., 2016). To cross-validate our results, we apply the $w\left(z\right)$ and D(z) from umbrella sampling to Eq. 5 and obtain an MFPT of 2.2 $\pm$ 0.02 $\mu$s, fairly similar to the 2.6 $\mu$s MFPT computed directly from milestoning. This consistency further demonstrated the robustness of the tested MD methods for computing transport kinetics.

Permeability Computed From Steady-State Flux

As mentioned above, under low concentration and constant electric field, we can simplify the GHK flux equation for computing the single-channel permeability $P$ from conductance measurement. Under symmetric concentration, the GHK flux equation can be written as

$P=\frac{\gamma RT}{q^2{F}^{2}C}=\frac{\gamma k_{B}T}{q^{2}C}(6)$

where $\gamma$ is the unitary conductance of a single channel, $R$ is the gas constant, $q$ is the charge of the permeating ion, $F$ is the Faraday constant, $C$ is the bulk concentration of the ion, and $k_{B}T$ has the same meaning above, except that the unit is eV here (0.026 eV at 300 K). At sufficiently small voltage, the current-voltage (I-V) relation is expected to be linear, and the slope defines the conductance.

Ionic conductance from MD simulations can be computed using two approaches. The most commonly used is a direct counting method, in which the currents were computed from the number of permeation events (N) over a simulation period (τ), I = N/τ. In our code (see Github link in Methods), the channel is split into upper, inner, and lower regions. A positive permeation event of K⁺ is counted if the time evolution of ion coordinates follows a lower-inner-upper sequence, and a negative permeation event is in reverse order. The carbon nanotube is applied with positive and negative 0.4 V voltage with 200 ns per replica. Each voltage simulation was repeated three times with different initial velocities. A least-square fitting of the I/V curve gave the conductance of 2.7 ± 0.94 ps by the direct counting method (Figure 4A).

A more efficient approach that does not rely on completed permeation events is to compute the instantaneous ionic current from charge displacement along the z-axis, $I\left(t\right)={\displaystyle \sum }_{i=1}^n{q}_i\left[z_1\left(t+\Delta t\right)-z_i\left(t\right)\right]\mathrm{∕}\Delta tL$ , in which $q_i$ and $z_i$ are the charge and z coordinate of ion i and L is the length of the channel pore (Aksimentiev and Schulten, 2005). This charge displacement method yielded a similar conductance of 2.3 ± 1.2 pS, indicating a good convergence of the voltage simulations (Figure 4B).

With an ionic concentration of 0.52 M (one single K⁺ in a cylinder bulk of radius 6 Å and length 28 Å), the permeability is (11.9 ± 6.36) × 10⁻¹⁶ cm³/s by the charge displacement method, and (13.6 ± 4.81) × 10⁻¹⁶ cm³/s by the direct counting method. In addition to ± 0.4 V simulations, ionic current calculated from 60 ns of ±0.3 V with sixfold higher concentration yields a similar permeability of 10.9 × 10⁻¹⁶ cm³/s. Only the results from low concentration are reported in Table 1.

Discussions

In this study, we used a carbon nanotube (CNT) as a toy model of a small conductance ion channel and computed the single-channel permeability from umbrella sampling, Markovian milestoning, and steady-state flux under voltage (Table 1). The PMF and MFPT for a single K⁺ permeating through the CNT produced from Markovian milestoning and umbrella sampling are in good agreement. Milestoning with cylindrical bulk restraint and spherical bulk restraint were tested and yielded consistent MFPTs when the effective bulk concentration is accounted. The single-channel permeability from voltage simulation is also within the same order-of-magnitude as those obtained from PMF-based and MFPT-based methods. These results are also in the same range as the previously reported K⁺ permeability of (25 ± 7) ×10⁻¹⁶ cm³/s computed using a transition path approach with the CHARMM22 force field (Zhou and Zhu, 2019).

It should also be noted that the current CNT model is chosen because it reproduces macroscopic properties (e.g., free energy, conductance) similar to common small-conductance ion channels. However, two carbon sheets were used as an artificial membrane to separate the solution. In the absence of a dielectric medium surrounding the channel transmembrane region, this toy model is unsuitable for investigating detailed electrostatic interaction between the ions and the channel.

In terms of computational resources, umbrella sampling has the advantage because PMF converges much faster than MFPT. However, MFPT allows extracting the kinetics directly from sampling. Thus, it does not rely on the assumption of ISD formulism and does not require additional calculations of position-dependent diffusion coefficient. Steady-state flux is straightforward to apply if the permeant is charged and sufficient sampling is achievable under reasonable voltages. However, if the permeant is not charged, a constant concentration gradient needs to be applied. Furthermore, the effect of an unphysiologically large voltage bias or electrochemical gradient on channel property is likely system-dependent and difficult to predict over a long simulation time. Compared with the steady-state flux approaches, the milestoning approach does not depend on the charge of the permeant. Thus, it can be used to study any type of small molecular permeation, such as the transport of a second messenger, cAMP, through a connexin26 hemichannel (Jiang et al., 2021a). Therefore, the use of milestoning has a significant promise for future applications on complex systems that are challenging to extract kinetics from unbiased MD or PMF-based enhanced sampling approaches.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

Y-CL performed all MD simulations. YL designed and supervised the project. Y-CL and YL analyzed data and wrote the paper together.

Funding

This work was supported by NIH Grants GM130834. Computational resources were provided via the Extreme Science and Engineering Discovery Environment (XSEDE) allocation TG-MCB160119, supported by NSF grant number ACI-154862.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are grateful to Dr. Fangqiang Zhu for sharing the coordinates of the CNT model and Drs. Luca Maragliano and Andrew Harris for inspiring discussions.

References

Adelman, J. L., and Grabe, M. (2015). Simulating Current-Voltage Relationships for a Narrow Ion Channel Using the Weighted Ensemble Method. J. Chem. Theor. Comput. 11, 1907–1918. doi:10.1021/ct501134s

PubMed Abstract | CrossRef Full Text | Google Scholar

Aksimentiev, A., and Schulten, K. (2005). Imaging α-Hemolysin with Molecular Dynamics: Ionic Conductance, Osmotic Permeability, and the Electrostatic Potential Map. Biophysical J. 88, 3745–3761. doi:10.1529/biophysj.104.058727

CrossRef Full Text | Google Scholar

Alberini, G., Benfenati, F., and Maragliano, L. (2018). Molecular Dynamics Simulations of Ion Selectivity in a Claudin-15 Paracellular Channel. J. Phys. Chem. B 122, 10783–10792. doi:10.1021/acs.jpcb.8b06484

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, T. W., Andersen, O. S., and Roux, B. (2004). Energetics of Ion Conduction through the Gramicidin Channel. Proc. Natl. Acad. Sci. U.S.A. 101, 117–122. doi:10.1073/pnas.2635314100

PubMed Abstract | CrossRef Full Text | Google Scholar

Awoonor-Williams, E., and Rowley, C. N. (2016). Molecular Simulation of Nonfacilitated Membrane Permeation. Biochim. Biophys. Acta (Bba) - Biomembranes 1858, 1672–1687. doi:10.1016/j.bbamem.2015.12.014

CrossRef Full Text | Google Scholar

Berne, B. J., Borkovec, M., and Straub, J. E. (1988). Classical and Modern Methods in Reaction Rate Theory. J. Phys. Chem. 92, 3711–3725. doi:10.1021/j100324a007

CrossRef Full Text | Google Scholar

Bernèche, S., and Roux, B. (2001). Energetics of Ion Conduction through the K+ Channel. Nature 414, 73–77. doi:10.1038/35102067

PubMed Abstract | CrossRef Full Text | Google Scholar

Carnevale, V., Delemotte, L., and Howard, R. J. (2021). Molecular Dynamics Simulations of Ion Channels. Trends Biochem. Sci. 46, 621–622. doi:10.1016/j.tibs.2021.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Chodera, J. D., Swope, W. C., Pitera, J. W., Seok, C., and Dill, K. A. (2007). Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations. J. Chem. Theor. Comput. 3, 26–41. doi:10.1021/ct0502864

PubMed Abstract | CrossRef Full Text | Google Scholar

Choudhary, O. P., Paz, A., Adelman, J. L., Colletier, J.-P., Abramson, J., and Grabe, M. (2014). Structure-Guided Simulations Illuminate the Mechanism of Atp Transport through Vdac1. Nat. Struct. Mol. Biol. 21, 626–632. doi:10.1038/nsmb.2841

PubMed Abstract | CrossRef Full Text | Google Scholar

Cottone, G., Chiodo, L., and Maragliano, L. (2020). Thermodynamics and Kinetics of Ion Permeation in Wild-type and Mutated Open Active Conformation of the Human α7 Nicotinic Receptor. J. Chem. Inf. Model. 60, 5045–5056. doi:10.1021/acs.jcim.0c00549

PubMed Abstract | CrossRef Full Text | Google Scholar

Darden, T., York, D., and Pedersen, L. (1993). Particle Mesh Ewald: AnN Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 98, 10089–10092. doi:10.1063/1.464397

CrossRef Full Text | Google Scholar

Diamond, J. M., and Katz, Y. (1974). Interpretation of Nonelectrolyte Partition Coefficients between Dimyristoyl Lecithin and Water. J. Membrain Biol. 17, 121–154. doi:10.1007/bf01870176

CrossRef Full Text | Google Scholar

Domene, C., Klein, M. L., Branduardi, D., Gervasio, F. L., and Parrinello, M. (2008). Conformational Changes and Gating at the Selectivity Filter of Potassium Channels. J. Am. Chem. Soc. 130, 9474–9480. doi:10.1021/ja801792g

PubMed Abstract | CrossRef Full Text | Google Scholar

Domene, C., Ocello, R., Masetti, M., and Furini, S. (2021). Ion Conduction Mechanism as a Fingerprint of Potassium Channels. J. Am. Chem. Soc. 143, 12181–12193. doi:10.1021/jacs.1c04802

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, D. J. (1983). Computer ''experiment'' for Nonlinear Thermodynamics of Couette Flow. J. Chem. Phys. 78, 3297–3302. doi:10.1063/1.445195

CrossRef Full Text | Google Scholar

Faradjian, A. K., and Elber, R. (2004). Computing Time Scales from Reaction Coordinates by Milestoning. J. Chem. Phys. 120, 10880–10889. doi:10.1063/1.1738640

PubMed Abstract | CrossRef Full Text | Google Scholar

Flood, E., Boiteux, C., Lev, B., Vorobyov, I., and Allen, T. W. (2019). Atomistic Simulations of Membrane Ion Channel Conduction, Gating, and Modulation. Chem. Rev. 119, 7737–7832. doi:10.1021/acs.chemrev.8b00630

PubMed Abstract | CrossRef Full Text | Google Scholar

Fowler, P. W., Abad, E., Beckstein, O., and Sansom, M. S. P. (2013). Energetics of Multi-Ion Conduction Pathways in Potassium Ion Channels. J. Chem. Theor. Comput. 9, 5176–5189. doi:10.1021/ct4005933

PubMed Abstract | CrossRef Full Text | Google Scholar

Hempel, T., del Razo, M. J., Lee, C. T., Taylor, B. C., Amaro, R. E., and Noe, F. (2021). Independent Markov Decomposition: Toward Modeling Kinetics of Biomolecular Complexes. P Natl. Acad. Sci. USA 118, e2105230118. doi:10.1073/pnas.2105230118

CrossRef Full Text | Google Scholar

Hille, B. (2001). Ion Channels of Excitable Membranes. 3rd ed. Sunderland, MA: Sinauer, p xviii, 814 p.

Google Scholar

Hoover, W. G., Ladd, A. J. C., and Moran, B. (1982). High-Strain-Rate Plastic Flow Studied via Nonequilibrium Molecular Dynamics. Phys. Rev. Lett. 48, 1818–1820. doi:10.1103/physrevlett.48.1818

CrossRef Full Text | Google Scholar

Hummer, G. (2005). Position-Dependent Diffusion Coefficients and Free Energies from Bayesian Analysis of Equilibrium and Replica Molecular Dynamics Simulations. New J. Phys. 7, 34. doi:10.1088/1367-2630/7/1/034

CrossRef Full Text | Google Scholar

Jiang, W., Del Rosario, J. S., Botello-Smith, W., Zhao, S., Lin, Y.-c., Zhang, H., et al. (2021b). Crowding-Induced Opening of the Mechanosensitive Piezo1 Channel In Silico. Commun. Biol. 4, 84. doi:10.1038/s42003-020-01600-1

CrossRef Full Text | Google Scholar

Jiang, W., Lin, Y.-C., Botello-Smith, W., Contreras, J. E., Harris, A. L., Maragliano, L., et al. (2021a). Free Energy and Kinetics of Camp Permeation through Connexin26 via Applied Voltage and Milestoning. Biophysical J. 120, 2969–2983. doi:10.1016/j.bpj.2021.06.024

CrossRef Full Text | Google Scholar

Khalili-Araghi, F., Tajkhorshid, E., and Schulten, K. (2006). Dynamics of K+ Ion Conduction through Kv1.2. Biophys. J. 91, L72–L74. doi:10.1529/biophysj.106.091926

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalili-Araghi, F., Ziervogel, B., Gumbart, J. C., and Roux, B. (2013). Molecular Dynamics Simulations of Membrane Proteins under Asymmetric Ionic Concentrations. J. Gen. Physiol. 142, 465–475. doi:10.1085/jgp.201311014

PubMed Abstract | CrossRef Full Text | Google Scholar

Kutzner, C., Grubmüller, H., de Groot, B. L., and Zachariae, U. (2011). Computational Electrophysiology: The Molecular Dynamics of Ion Channel Permeation and Selectivity in Atomistic Detail. Biophysical J. 101, 809–817. doi:10.1016/j.bpj.2011.06.010

CrossRef Full Text | Google Scholar

MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. J., et al. (1998). All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 102, 3586–3616. doi:10.1021/jp973084f

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackerell, A. D., Feig, M., and Brooks, C. L. (2004). 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. 25, 1400–1415. doi:10.1002/jcc.20065

PubMed Abstract | CrossRef Full Text | Google Scholar

Maragliano, L., Vanden-Eijnden, E., and Roux, B. (2009). Free Energy and Kinetics of Conformational Transitions from Voronoi Tessellated Milestoning with Restraining Potentials. J. Chem. Theor. Comput. 5, 2589–2594. doi:10.1021/ct900279z

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, B. (1998). Molecular Dynamics Simulations of Ion Channels: How Far Have We Gone and where Are We Heading? Biophysical J. 74, 2744–2745. doi:10.1016/s0006-3495(98)77981-0

CrossRef Full Text | Google Scholar

Roux, B. (2008). The Membrane Potential and its Representation by a Constant Electric Field in Computer Simulations. Biophysical J. 95, 4205–4216. doi:10.1529/biophysj.108.136499

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirts, M. R., and Chodera, J. D. (2008). Statistically Optimal Analysis of Samples from Multiple Equilibrium States. J. Chem. Phys. 129, 124105. doi:10.1063/1.2978177

PubMed Abstract | CrossRef Full Text | Google Scholar

Szabo, A., Schulten, K., and Schulten, Z. (1980). First Passage Time Approach to Diffusion Controlled Reactions. J. Chem. Phys. 72, 4350–4357. doi:10.1063/1.439715

CrossRef Full Text | Google Scholar

Teo, I., and Schulten, K. (2013). A Computational Kinetic Model of Diffusion for Molecular Systems. J. Chem. Phys. 139, 121929. doi:10.1063/1.4820876

PubMed Abstract | CrossRef Full Text | Google Scholar

Vanden-Eijnden, E., Venturoli, M., Ciccotti, G., and Elber, R. (2008). On the Assumptions Underlying Milestoning. J. Chem. Phys. 129, 174102. doi:10.1063/1.2996509

PubMed Abstract | CrossRef Full Text | Google Scholar

Vanden-Eijnden, E., and Venturoli, M. (2009). Markovian Milestoning with Voronoi Tessellations. J. Chem. Phys. 130, 194101. doi:10.1063/1.3129843

PubMed Abstract | CrossRef Full Text | Google Scholar

Venable, R. M., Krämer, A., and Pastor, R. W. (2019). Molecular Dynamics Simulations of Membrane Permeability. Chem. Rev. 119, 5954–5997. doi:10.1021/acs.chemrev.8b00486

PubMed Abstract | CrossRef Full Text | Google Scholar

Votapka, L. W., Lee, C. T., and Amaro, R. E. (2016). Two Relations to Estimate Membrane Permeability Using Milestoning. J. Phys. Chem. B 120, 8606–8616. doi:10.1021/acs.jpcb.6b02814

PubMed Abstract | CrossRef Full Text | Google Scholar

Woolf, T. B., and Roux, B. (1994). Conformational Flexibility of O-Phosphorylcholine and O-Phosphorylethanolamine: A Molecular Dynamics Study of Solvation Effects. J. Am. Chem. Soc. 116, 5916–5926. doi:10.1021/ja00092a048

CrossRef Full Text | Google Scholar

Zheng, J., and Trudeau, M. C. (2015). Handbook of Ion Channels. 2nd Edition. Boca Raton: CRC Press.

Google Scholar

Zhou, X., and Zhu, F. (2019). Calculating Single-Channel Permeability and Conductance from Transition Paths. J. Chem. Inf. Model. 59, 777–785. doi:10.1021/acs.jcim.8b00914

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, F., and Hummer, G. (2012). Theory and Simulation of Ion Conduction in the Pentameric Glic Channel. J. Chem. Theor. Comput. 8, 3759–3768. doi:10.1021/ct2009279

CrossRef Full Text | Google Scholar