Brownian Dynamics Simulations of Ion Transport through the VDAC

GCMC/BD simulation system setup

The GCMC/BD simulation system (Fig. 1) was constructed as follows: A 40-Å-thick implicit membrane with its normal along the Z axis was placed in the middle of an 80 × 80 × 90 Å³ box centered at the origin. Each hVDAC1 NMR structure was oriented with its pore axis parallel to the Z axis, its center of mass at Z = 0 Å, and its N- and C-termini facing the cis-side (cytosol, Z < 0). hVDAC1 spans from −20 Å to 20 Å in the Z-direction, with an outer diameter of ∼35 Å and inner diameter of ∼20 Å. The 5-Å-thick buffer regions were included at the top and bottom of the simulation box to mimic the constant concentration baths for mobile ions (K⁺ and Cl⁻ in this work). The system consisted of three different dielectric regions: an ion-accessible solvent (water) region, a protein region, and a membrane region. The latter two are inaccessible to mobile ions. The dielectric constants for water, protein, and lipid were set to 80, 2, and 2, respectively (15,18,20). The bulk water dielectric constant was assigned in the pore region because of the relatively large pore size of hVDAC1. Although the protein dielectric constant may affect the macroscopic transport properties, as pointed out in other coarse-grained simulation studies (21,22), we used the most common protein dielectric constant of 2 (15,18,20), and there were negligible effects on ion conductance and ion distribution with the protein dielectric constant of 4 (data not shown). The reference of the TM potential (V_mp) is the cis-side (Z < 0), and thus the TM potential in this work is defined as V_mp = V_trans − V_cis.

GCMC/BD simulation

We used the open-source GCMC/BD simulation program (16,17), which is freely available at http://thallium.bsd.uchicago.edu/RouxLab/gcmc.html. In GCMC/BD, only mobile ions are explicitly simulated. Water and membrane are treated as structureless continuum dielectric regions, and the channel protein atoms are modeled as fixed charges with finite atomic radii. Taking into account the spatial variation of the diffusion constant (D_α), the ion movement is described by the BD equation (23):

where r_α, {r_α}, W, and ζ_α are the position of ion α, the position of all ions in the system, the multi-ion potential of mean force (PMF), and the random Gaussian noise representing the frequent bombardment of water molecules, respectively. The multi-ion PMF is expressed as:

where ϕ_sf is the electrostatic potential from the fixed protein charges and V_mp, ϕ_rf is the reaction field potential arising from dielectric heterogeneity, U_core is the core-repulsive steric potential to prevent ions from moving into ion-inaccessible regions, and u_αγ describes the ion-ion interactions between ions α and γ, including the dielectric-screened Coulombic electrostatic interaction, van der Waals interaction, and water-mediated short-range interaction (15). Using the Poisson-Boltzmann (PB)/PNP program (available at http://thallium.bsd.uchicago.edu/RouxLab/pbeq.html), we calculated ϕ_sf, the parameters for ϕ_rf (20), and U_core, and tabulated them on a rectangular grid with a grid spacing of 0.5 Å before running the main GCMC/BD simulation. Thus, we were able to efficiently calculate the multi-ion PMF in Eq. 2 and the corresponding force on each ion at each BD time step of 10 fs. In the buffer regions, we performed one GCMC step at every BD step to maintain the specified ion concentration according to the ion creation and annihilation probability. A detailed description of the GCMC/BD approach and its application to biological systems can be found in elsewhere (16,17). Because the conformation change during the voltage gating process of VDAC cannot be captured with the present GCMC/BD simulation using a fixed channel structure, a study of the voltage gating mechanism is beyond the scope of this work.

Diffusion constant

One of the most important physical parameters in the theoretical study of ion transport is the diffusion constant, which directly affects the ion transport rate in coarse-grained approaches such as BD (see Eq. 1) and PNP calculations. Several different approaches can be used to model the space-dependent diffusion constant of mobile ions in a physically restricted pore, including hydrodynamic ion-solvent modeling (18), fitting of measured current-voltage (I/V) curves (24), and direct calculation from MD trajectory (15).

We used the diffusion constant profile obtained from the MD simulations of NMR01 in 1.0 M KCl solution (see Fig. S1 in the Supporting Material; for further details, see Rui et al. (19)). The average scaling factor of the diffusion constant inside the pore (i.e., the diffusion constant inside the pore normalized by the bulk value) is 0.55 for both K⁺ and Cl⁻. However, the calculated bulk diffusion constants for both K⁺ and Cl⁻ from the MD simulations are 0.30 Å²/ps and 0.29 Å²/ps, which are ∼1.5 times larger than experimental values (0.20 Å²/ps for both K⁺ and Cl⁻). The higher bulk diffusion constants resulted from the modified TIP3P water model used in the MD simulation with the particle mesh Ewald method (25), which is known to underestimate the solvent viscosity and hence overestimate the particle diffusion constants in solution (26). Therefore, we used MD-derived bulk diffusion constants in the comparison of BD and MD simulations, and the experimental bulk diffusion constants in the comparison of the BD simulations and experimental data, assuming the same scaling profile along the Z direction. A third-order polynomial function with a switching length of 17 Å was employed to approximate the diffusion constant between the bulk and pore regions (15).

Calculation of channel conductance and ion selectivity

The channel conductance (G = I/V_mp) can be calculated directly by counting the number of ions that cross the XY plane at Z = 0 at a given TM potential. In this work, we consider two symmetric KCl concentrations for conductance calculation: 1.0 M for comparison of BD and MD simulations (19), and 0.1 M for comparison of the BD simulations and experimental data. The simulation times were 100 ns for 1.0 M KCl and 400 ns for 0.1 M KCl.

Both the ratio of ion currents (e.g., I_Cl/I_K) from a symmetric KCl solution and the permeability ratio from an asymmetric solution can be a measure of ion selectivity. Because it is still technically challenging to maintain an asymmetric bath condition in MD simulations and it is impractical to separate ion current components from the total current in experiments, the current (conductance) ratios from BD simulations were compared with those from MD simulations, and the permeability ratios (or the reversal potential) from BD simulations were compared with experimental data. The reversal potential (V_rev) was measured from an I/V curve as the potential at which the total current becomes zero. With the bath concentrations and reversal potential, the permeability ratio (p) can be estimated by the Goldman-Hodgkin-Katz equation (27):

We used a 10-fold KCl concentration gradient between cis- and trans-side baths (C^cis = 1.0 M and C^trans = 0.1 M) for permeability calculations with the BD simulation time of 200 ns.

Comparison of BD and MD

Our BD simulation uses a hierarchical approach to closely reproduce key microscopic properties of ion transport, such as ion diffusion constant profiles and ion-ion interactions observed in MD simulation. Nonetheless, it is important to assess the efficacy of the coarse-grained BD simulation because various approximations are still used to enhance simulation efficiency in the BD approach, such as a (fixed dielectric) rigid channel structure, implicit solvent-based protein-ion interactions (i.e., grid-based steric and electrostatic potentials, including the reaction field in Eq. 2), and an implicit membrane bilayer. In this work, BD and MD results are compared in terms of the I/V relationship, the ion current (conductance) ratio, and ion number densities inside the pore. For the direct comparison with MD simulations, we used the MD-derived diffusion constants in the BD simulations, which are ∼1.5 times larger than the experimental values (see Materials and Methods).

Fig. 2 compares the I/V curves from the BD and MD simulations of NMR01 in 1.0 M KCl solution with V_mp of 0, ±25, ±50, ±75, and ±100 mV. A 1.0 M KCl solution was used to enhance the sampling of ion transport events in the MD simulation. The I/V curves show excellent overall agreement between BD and MD for K⁺, Cl⁻, and total currents. The I/V curves from BD also show slight asymmetry (|I(−100 mV)| > |I(100 mV)|) as observed in MD simulation results, although the I/V relationship is symmetric at lower TM potentials up to ±50 mV. The average total conductance from BD (G = 4.9 ± 0.1 nS) equals the MD conductance (G = 4.9 ± 0.1 nS; the value after ± corresponds to the standard error throughout this work). The higher diffusion in the BD and MD simulations is responsible for the resultant higher conductance compared with the experimental value of 3.9–4.5 nS for various VDACs (28). Table S1 summarizes the current ratios from both the BD and MD simulations. The current ratios from BD remain almost constant regardless of V_mp and agree with the MD results of the positive V_mp region. The average current ratio from BD is 1.7 ± 0.1, which indicates the anion preference of the channel. The discrepancy between the current ratios from BD and MD—more specifically, the much higher ratios in MD on the negative V_mp side—is mostly due to the limited number of K⁺ crossing in the MD simulations (19).

Table 1 summarizes the average number of K⁺ and Cl⁻ in the pore region (−15 Å ≤ Z ≤ 15 Å) at equilibrium (V_mp = 0 mV) and nonequilibrium (V_mp = 50 mV) conditions. Regardless of V_mp, a higher number of Cl⁻ than K⁺ is found in the pore, which also indicates the anion preference of the channel. The minimal effects of V_mp on the ion distribution are observed inside the pore because the protein-ion interactions are the main determinant for the distribution. Similar results were obtained in the MD simulations under the same conditions, but with slightly elevated numbers for both K⁺ and Cl⁻ in the pore. The ion number profiles along the Z axis (Fig. S2) clearly show that the difference between the pore ion numbers in BD and MD (Table 1) arises mostly from the different number of K⁺ at the center. In the MD simulations, trapping of K⁺ (due to less screened interactions between K⁺ and acidic residues) at the center of the channel was observed (19). The trapped K⁺ ions attract nearby Cl⁻, consequently resulting in the slightly higher number of both ions inside the pore in the MD simulations. Although such detailed protein-ion interactions are approximated by implicit solvent-based protein-ion interactions (Eq. 2), the superposition of the ion trajectories from BD is similar to that from the MD simulations, as shown in Fig. S3. Accumulated K⁺ ions are found close to several charged or polar residues, including Asp-16, Asp-30, Glu-84, and Asn-207. Mutations of several corresponding residues, such as Asp-16 to Lys, and Asp-30 to Lys, in Saccharomyces cerevisiae yeast VDAC (scVDAC) increase the anion selectivity of the channel (29) (see also the mutation study below).

We also performed the BD simulations with explicit lipid molecules taken from a snapshot in the MD simulations (19). The differences in the conductance and number of ions between implicit and explicit membranes are very small (data not shown), indicating that the use of a simple implicit membrane is a reasonable approximation in the BD simulation. Given the aforementioned approximations in BD, this comparison clearly illustrates that the hierarchical approach in BD captures the essential interactions that govern ion transport through hVDAC1, and reproduces the MD results within the limitation of such approximations. Similar observations in terms of the efficacy of the hierarchical BD simulation were reported in previous MD and BD simulation studies of α-Hemolysin (18,30) and α-Hemolysin with a molecular adaptor (20,31).

Ion transport properties of hVDAC1 NMR models

Given the current debate on the validity of the reported hVDAC1 structures (7), it is necessary to thoroughly investigate the structure-channel property relationship of each NMR structure and the structure ensemble by calculating the ion transport properties under physiologically relevant conditions. To compare the BD results directly with available experimental data, we used the same diffusion constant profile (Fig. S1) with the experimental ion diffusion constants for the bulk region (see Materials and Methods).

Table S2 summarizes the single-channel conductance of each hVDAC1 NMR model in 0.1 M KCl solution with V_mp = 50 mV and its permeability ratio (p) in an asymmetric 1.0:0.1 M KCl solution. Fig. 3 shows a typical I/V curve in the asymmetric solution, and p is calculated using Eq. 3 by determining the reversal potential from the I/V curve. The experimental conductance and permeability ratio in corresponding conditions are 0.4–0.6 nS and ∼2.0 (4). From the experimental values, a selection window of conductance and permeability ratio (0.4 nS < G < 0.6 nS and 1.5 < p < 2.5) was chosen to assess the ion transport properties of the NMR structures. As shown in Fig. 4, all NMR structures except NMR02 have permeability ratios > 1.0, indicating anion selectivity of the structures. However, only three NMR models (NMR01, NMR03, and NMR20) yield a conductance and permeability ratio that is not significantly different from the experimentally measured range (considering a 95% confidence interval), and NMR03 best satisfies both selection criteria. No model yields both low conductance and cation selectivity resembling the closed state (4), demonstrating that the NMR structures represent open-state channel structures.

Assuming an equal probability of observing each NMR model, we calculated the distributions of the single-channel conductance and the reversal potentials from the 10 independent simulation runs for all of the structures, and the results are shown in Fig. 5. The reversal potential distribution is determined by linear interpolation of the I/V curves for each structure. The ensemble-averaged conductance is 0.36 ± 0.04 nS, which is not significantly below the lower bound of the conductance selection window considering a 95% confidence interval (0.30 nS < G < 0.43 nS; Table S2). Also, compared to the channel conductance (1.06 nS) from the previous PNP calculation (14), the BD simulations clearly yield a much more reasonable single-channel conductance that is comparable to the range of experimental measurement, demonstrating the importance of including the finite ion size and the direct ion-ion interactions in the pore (15). The average reversal potential is 23.7 ± 2.2 mV (corresponding to p = 3.2), which is higher than the upper bound in the selection criteria of 18.95 mV (p = 2.5). The ensemble properties based on the uniform distribution of the NMR models show good agreement with experimental conductance, whereas the anion selectivity is a bit higher than the experimental data. This observation supports the conclusion that the NMR models represent the anion-selective, open-state hVDAC1 structures.

Structure and ion transport properties for selected structures

Ultimately, one must calculate the channel property as an ensemble property that arises from an ensemble of channel structures. However, since most BD or PNP studies are based on a single, rigid representative structure, it is instructive to investigate the structure-channel property relationship with a few selected NMR structures. For this purpose, three NMR models (NMR02, NMR03, and NMR11) were chosen because they have similar conductance but different permeability ratios. Table 2 summarizes their ion transport properties, including conductance, number of ions in the pore (−15 Å ≤ Z ≤ 15 Å), and permeability ratio. Both the conductance ratio and the permeability ratio show an increasing trend of anion selectivity in the order of NMR02, NMR03, and NMR11. Of interest, variations of the K⁺ number in the pore and the K⁺ conductance among the three models are large, whereas those of Cl⁻ are small. Thus, one can envision that the permeation process of K⁺ would differ among the models, and the structural differences would obviously affect the K⁺ permeation process, resulting in different ion selectivity, whereas the Cl⁻ permeation process would remain similar.

To characterize the microscopic origin of the anion selectivity difference in these models, we calculated the one-dimensional (1D) multi-ion PMFs (32) for K⁺ and Cl⁻ along the Z axis from the average equilibrium ion density profile (C_α):

where α is the ion type and C_ref is the arbitrary reference concentration that is set to the bulk concentration to make ${\mathrm{W}}_{\alpha ,1\text{D}}\left(Z\right)\to 0$ as $\left|Z\right|\to \infty$. ${\mathrm{W}}_{\alpha ,1\text{D}}$ represents the mean potential acting on an ion along the channel axis, including the contributions from protein-ion, ion-ion, and ion-solvent interactions, and the volume exclusion effect (i.e., the variation in the ion-accessible cross-sectional area).

Fig. 6 shows the 1D multi-ion PMFs of the three NMR models with and without protein charges. The PMFs without protein charges include the volume exclusion effect and the reaction field from the membrane-solvent dielectric boundary, and thus serve as the reference PMFs to consider the influence of protein-ion interactions. The volume exclusion effect causes the barriers in all PMF profiles near the pore entrance (Z ∼ ±20 Å). Clearly, the protein charges are the characteristic determinants of the PMFs inside the pore. The Cl⁻ PMFs have lower barriers than the reference PMFs, and the barriers of the K⁺ PMFs are higher, which explains why the Cl⁻ currents in the BD simulations are larger than the K⁺ currents (Table 2). The similar Cl⁻ currents in the three models are also dictated by similar, rather flat (∼1.0 kcal/mol) PMFs inside the pore. The slight differences in the Cl⁻ currents are determined by the PMF barriers in 10 Å ≤ Z ≤ 20 Å, i.e., the order in the Cl⁻ currents (NMR03 > NMR11 > NMR02 in Table 2) is inversely proportional to the barrier heights (NMR02 > NMR11 > NMR03 in Fig. 6). Compared to the Cl⁻ PMF, the variation of the K⁺ PMFs in the pore among the selected models is significant. Clearly, the lowest K⁺ current in NMR11 results from the two large PMF barriers, which are higher than those in NMR02 and NMR03.

A discussion about the relationship between the PMFs and the protein structure, including the position and orientation of the protein residues, can provide detailed insight into the structure-channel property relationship. Of the three models, NMR11 has the smallest ion-accessible area profile, whereas NMR02 and NMR03 have similar area profiles (Fig. S4 A). Hence, one would expect NMR11 to conduct a smaller ionic current because it has the narrowest pore, and the electrostatic effect from protein residues (especially the charged residues) would be amplified in NMR11 due to crowding of the charged residues. This notion is illustrated in Fig. 6, which shows that the reference PMF of NMR11 has a larger barrier than the other models, and the PMF change considering the protein charge effect is the largest in NMR11. The effect of the smaller pore in NMR11 is also clearly captured at the PMF peaks at Z = −8 Å, where pore-facing Lys-61 and Lys-96 are positioned to strongly hinder the K⁺ permeation and attract Cl⁻. The electrostatic effect from these residues is stronger in NMR11 than in NMR02 and NMR03, as the Z position corresponds to the narrowest pore by the N-terminus (see Fig. S4 B for the cross-sectional average ion density in the same Z position).

In addition, the K⁺ PMF barrier in NMR11 is enhanced because Glu-84 (near Z = 0 Å) is buried in the lipid region (Fig. S4 B), which reduces the attraction of K⁺ into the pore, whereas the side chain of Glu-84 in NMR02 points directly into the pore (Fig. S4 C). The side chain of Glu-84 in NMR03 is located inside the channel but its orientation is parallel with the pore wall. Hence, the order of local K⁺ PMF minima near Z = 0 Å (NMR02 < NMR03 < NMR11) can be attributed to the relative orientation of Glu-84. Among the charged residues in the N-terminal region, Arg-15 exemplifies the effect of the different orientation of the residue on the ion permeation. As shown in Fig. S4 D, the side chain of Arg-15 in NMR03 and NMR11 is located at the XY center of the pore and faces the pore, which results in the deep Cl⁻ PMF well near Z = 8 Å. In contrast, the wall-facing Arg-15 in NMR02 results in weaker attraction of Cl⁻ at this position. Consequently, higher densities of Cl⁻ are observed near Arg-15 in NMR03 and NMR11, but not in NMR02. Asp-16 is also found to strongly affect the K⁺ PMF barrier, which causes a decrease of the K⁺ PMF from the reference PMF near the trans-side (Z > 0) pore entrance for NMR02 and NMR03. Localized K⁺ trajectories are observed between Asp-16 and the cluster of the Asn residues (Asn-168, Asn-183, Asn-185, and Asn-207) in NMR02 and NMR03. However, in NMR11, the N-terminal region is positioned between Asp-16 and the Asn cluster to interrupt the K⁺ attraction site and cause the reduction of the K⁺ permeation.

The results of these detailed analyses show that one can explain the microscopic origins of the different ion selectivities among the models by examining the PMF, ion trajectories, and structural information. Specifically, the conformation of the flexible N-terminal region and orientation of its charged residues play important roles in determining the overall conduction properties of the NMR models.

Mutation effects on ion transport properties

Mutation studies have been used to characterize the role of specific residues in ion transport properties. The ability to reproduce or predict mutation effects from computational studies is important because such studies can explain specific mutation-induced changes on a microscopic level. We selected seven mutation sites from the experimental measurements of yeast scVDAC (29,33): D15K, K19E, D30K, K61E, K95E, E220K, and K274E. Previous experimental functional studies indicated a structural similarity between scVDAC and hVDAC1 (34). The corresponding residues in hVDAC1 are D16K, K20E, D30K, K61E, K96E, G220K, and K274E. The five mutations of D15K, K19E, D30K, K61E, and K95E in scVDAC were found to notably affect the ion selectivity, whereas those of E220K and K274E resulted in negligible selectivity changes (29,33). The positions of the mutation sites are depicted in Fig. 7, A and B.

Because NMR03 shows the best agreement with the experimental conductance and reversal potential, we used NMR03 as a template wild-type (WT) structure for the aforementioned mutants. Table 3 summarizes the reversal potentials (V_rev) calculated from the BD simulations and compares them with the experimental data from scVDAC (33). Given the assumption that each mutation does not induce a significant structural change from NMR03, the BD simulation reproduces the experiment reasonably well except for D16K and K274E (see also Fig. 7 C, which compares the permeability ratio change (R_p = p_mutation/p_WT) between the BD simulations and the experimental data).

The changes in V_rev are correlated with those in the 1D multi-ion PMF barriers (Fig. S5). The mutations' impact on the PMF is highly dependent on the residue position. Mutations near the pore entrance are likely to result in a local change of the PMF. For example, the K20E mutation is located near the pore entrance on the trans side, and it has only a slight impact on the Cl⁻ PMF near its location (10 Å ≤ Z ≤ 20 Å). On the other hand, mutations that occur inside the channel pore affect the overall K⁺ and Cl⁻ PMFs. This is illustrated by the K61E mutation on the channel wall inside the pore, which perturbs both the K⁺ and Cl⁻ PMFs over a wide range along the Z axis. Even though these mutations affect the ion PMFs differently, they result in similar shifts in V_rev.

A large discrepancy between the BD simulation and the experiment is found for D16K and K274E (Fig. 7 C and Table 3). In the case of D16K, whereas V_rev changes by 5.5 mV in the experiment, a larger increase by 14.2 mV is observed in the BD simulation. This is because Asp-16 is located in the flexible N-terminal helix, and its mutation to Lys would result in a different conformation and movement of the helix, which is not considered in the current BD simulation of the mutant. The K274E mutation also shows a negligible V_rev change in the experiment, but its BD simulation yields the largest decrease (cation selectivity). This arises from the perturbation of the Cl⁻ PMF by the K274E mutation on β-strand 19 at the trans-side pore entrance; the mutation inverts the deep potential well in the WT to a large barrier near Z = 7.5 Å (Fig. S5 D). This considerable alteration in both the PMF and the reversal potential is again attributed to the rigid-body assumption in BD. In the hVDAC1 structure, the side chain of Lys-274 is toward the channel lumen. Consequently, mutating Lys-274 to Glu affects the electrostatic potential inside the pore and causes a significant shift in V_rev. Considering the location of Lys-274 close to the pore entrance and the flexibility of the native VDAC structure, the side chain of Glu-274 from the mutation can be positioned outside of the pore. Therefore, the K274E mutation may have little effect on V_rev in experimental mutagenesis studies of scVDAC (7,34).

In general, the BD simulations reproduced the experimental ion selectivity change from single-point mutations reasonably well. The discrepancies observed for the D16K and K274E mutations are likely due to possible conformation changes, structural differences of mutated residues (between in hVDAC1 and in scVDAC), and the limitation of the fixed-structure approximation in the BD approach.