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

(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 and D(z) (see methods) (E). Local resistance (the integrand of Eq. 3) for permeating K+ (F). Integration of the permeation resistance, 1/ , as a function of z.

FIGURE 2

(A). Raw data plotted along the channel z-axis and radial distance 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 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

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

Method	Umbrella sampling	Markovian milestoning	Steady-state flux
Permeability equation	P=πr^2(∫z_1z_2e^βw(z)D(z)dz)^−1	P=πr^2∫Z1Z2e^−w(Z)k_BTdz2×MFPT	P=γk_BTq^2C
Effective bulk conc.	0.77 M	0.77 M	0.52 M
Biasing potential	Harmonic	Flat-bottom harmonic	Voltage ±0.4 V
Biasing force constant	2.5 kcal/mol	100 kcal/mol	—
Lateral restraint in bulk	6 Å	6 Å	—
PMF barrier (kcal/mol)	3.8 ± 0.17	4.1 ± 1.7	—
Key parameters	D(z) 0.004 ∼0.028 Å2/ps	MFPT: 2.6 μs	Conductance 2.3 ± 1.2 pS
Permeability (cm3/s)	(8.96 ± 0.02) ×10−16	2.92 × 10–16	(11.7 ± 6.3) ×10−16
Total sampling time (μs)	1.12	4.98	0.60

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 (cm3/s) can be related linearly to the rate of crossing or mean first passage time (MFPT) under equilibrium, , in which is the symmetric solute concentrations. Here, we use the number of molecules per second for , seconds per molecule for MFPT, and molecule/cm3 for .

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 ( ) and the permeability ( can be related by the GHK flux equation, where is the charge of the permeant, is the voltage, F is the Faraday constant, R is the gas constant, T is the absolute temperature, and is the concentration (Hille, 2001). When applied to a membrane-embedded single-channel model and ions only cross the membrane through the channel, in the GHK equation corresponds to the single-channel permeability. GHK flux equation thus relates single-channel conductance , a nonequilibrium property, to the equilibrium property .

It can be seen from the equations above that the crossing rate , MFPT, and conductance are all concentration-dependent; only is independent of concentration. Experimentally, 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 Å3, 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 using , where is the deviation of the z-position of the ion at time , , from the time-averaged position in each window. 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 in the harmonically restrained umbrella sampling gives .

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, , and the equilibrium probability π for the tagged ion to be in cell B i satisfies a balance equation:

The free energy of each cell can be obtained from the solution of Eq. 1 as −k Tln(π ) (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 , the rate of moving from milestone S ij to S ik, are given by where is the number of transitions from S ij to S ik, normalized by the time spent in cell B 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 and 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−1.

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

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 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 , which can be approximated to in a homogenous bulk, where is the radius of the cylinder. Hence, defines the effective bulk concentration in the simulated region, which led the probability of the ion inside the channel over the true bulk density to be , where is the PMF with the bulk value set to zero at the cylindrical region and (Allen et al., 2004; Zhu and Hummer, 2012). is the temperature, and is Boltzmann’s constant. Therefore, at low ionic concentration, single-channel permeability can be estimated using a slightly modified ISD equation:

In Eq. 3, is the position-dependent diffusion constant of the studied permeant along the z-axis (Figure 1D). The interval of the integration, , is the lower and upper boundaries of the channel pore, beyond which PMF reaches the bulk value. is the radius of the cylindrical restraint when the ion is outside of the interval. It is necessary to set 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 = 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 in Figure 1C and D(z) in Figure 1D, the permeability estimated from Eq. 3 is (8.96 ± 0.02) × 10−16 cm3/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/ , as a function of the z-axis (Figure 1F). It is not surprising that the 1/ 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−1 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−1 and 0.5 ps−1 yield an identical and symmetric MFPT of 2.6 ± 0.03 μs for K+ permeation.

Using the PMF, and MFPT, from the milestoning (with 0.5 ps−1 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−16 cm3 s−1, in fairly good agreement with the permeability of 8.96 × 10−16 cm3s−1 computed from ISD equation (Eq. 3) using umbrella sampling data:

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 (cm3/s) can be related linearly to the mean first passage time (MFPT) 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, , can be written as

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 and D(z) from umbrella sampling to Eq. 5 and obtain an MFPT of 2.2 0.02 s, fairly similar to the 2.6  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 from conductance measurement. Under symmetric concentration, the GHK flux equation can be written as where is the unitary conductance of a single channel, is the gas constant, is the charge of the permeating ion, is the Faraday constant, is the bulk concentration of the ion, and 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, , in which and 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−16 cm3/s by the charge displacement method, and (13.6 ± 4.81) × 10−16 cm3/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−16 cm3/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−16 cm3/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.