Simulating Current-Voltage Relationships for a Narrow Ion Channel Using the Weighted Ensemble Method

Abstract

Ion channels are responsible for a myriad of fundamental biological processes via their role in controlling the flow of ions through water-filled membrane-spanning pores in response to environmental cues. Molecular simulation has played an important role in elucidating the mechanism of ion conduction, but connecting atomistically detailed structural models of the protein to electrophysiological measurements remains a broad challenge due to the computational cost of reaching the necessary time scales. Here, we introduce an enhanced sampling method for simulating the conduction properties of narrow ion channels using the Weighted ensemble (WE) sampling approach. We demonstrate the application of this method to calculate the current–voltage relationship as well as the nonequilibrium ion distribution at steady-state of a simple model ion channel. By direct comparisons with long brute force simulations, we show that the WE simulations rigorously reproduce the correct long-time scale kinetics of the system and are capable of determining these quantities using significantly less aggregate simulation time under conditions where permeation events are rare.

Figure 1

A model ion channel. (A) A cut-away view of the solvated model channel showing the solvated pore. The particles comprising the hydrophobic sheets mimicking the membrane are shown as gray spheres, whereas the charged atoms of the intervening pore are blue, gold, gray and light blue and carry charges of −0.5 e, +0.5 e, −0.35 e, and +0.35 e, respectively. Mobile Na⁺ (green) and Cl⁻ (cyan) ions are shown in the bath of SPC/E water molecules. (B) For the WE simulations, the simulation box in (A) is discretized into bins, shown here for the x and z dimensions. A cartoon of the membrane and pore atoms are overlaid to highlight the fine discretization through the pore and the coarser binning in the bulk regions. The bin spacing in the y dimension (not shown in this projection), is identical to binning in x.

Figure 2

Effective steady-state energy profiles for Na⁺ and Cl⁻ as a function of ion position at 0.55 V calculated from the projection of the steady-state ion distribution along the z axis of the simulation box. Profiles calculated from the brute force simulations are shown as solid lines for the low (left) and high (right) concentration systems. The profiles for the corresponding WE simulations are calculated from the final 400 iterations and are shown as open circles.

Figure 3

Effective steady-state energy profiles for the pairwise steady-state distribution of Na⁺ at 0.55 V for the high concentration system along the z axis of the simulation box. The probability of finding a pair of Na⁺ ions at positions z₁ and z₂, P (z₁, z₂), is calculated over all 8 cations in the box using either all of the brute force data (left panel) or the final 400 iterations of the WE simulation (right panel).

Figure 4

I–V relationship of the model ion channel. At each applied voltage, the current carried by the Na⁺ and Cl⁻ is shown for the brute force trajectories as well as the corresponding WE simulations. The latter are analyzed using both the direct and reweighting procedures. Simulations carried out using the low concentration system are shown in the upper panel and simulations using the high concentration system are shown in the lower panel. Error bars indicate the 95% confidence interval. The solid lines (Na⁺) and dashed lines (Cl⁻) are spline fits to the data and carry no theoretical meaning.

Figure 5

Convergence of the current estimates from WE simulations. For the simulations carried out at 0.55 V, the estimates of the current as a function of aggregate simulation time are shown for the Na⁺ current (A and B) and the Cl⁻ current (C and D). Panels B and D show the same data as in A and C, respectively, but are plotted on a log-log scale to highlight the estimated currents early in the simulations. The solid black line is the current estimated from brute force trajectories, with dashed lines corresponding to the evolution of the 95% confidence interval calculated using Eq. 6 assuming n = k_ref T, where k_ref is the estimate from all of the available brute force simulation data. The direct and reweighted curves for each ion type are calculated from a single simulation discarding the first 50% of the WE iterations up until a given time point. The restart simulations (direct - red and reweighted - purple) are initiated from a converged WE simulation with an applied voltage of 0.6 V.

Figure 6

Efficiency of the WE method as a function of observed current. Efficiency is defined as the ratio of the total aggregate simulation time required to reach and maintain an error ≤1, T₁, to the MFPT of the system at a particular current. The error is defined in Eq. 7. Points represent all of the WE simulations for both Na⁺ and Cl⁻ at the two ion concentrations. The dashed line at T₁/MFPT = 1 corresponds to the ratio expected from a conventional brute force simulations over the same range of currents.