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

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

Helix insertion energy of a model polyleucine helix with a central arginine

As the arginine enters the membrane, the upper leaflet bends to allow water penetration. At the upper leaflet, the arginine height is 21 Å, and it is 0 Å at the center. (A) The total electrostatic energy remains nearly constant upon insertion. (B) The nonpolar energy increases linearly to 18 kcal/mol as the membrane bending exposes buried TM residues to water. (C) The membrane deformation energy redrawn from . (D) The total helix insertion energy is the sum of A–C plus ΔG, which is not shown (solid red line). Correcting for the optimal membrane deformation at a given arginine depth, as shown in , produces a noticeably smaller insertion energy (dashed red line). Our continuum computational model matches well with results from fully atomistic MD simulations on the same system (diamonds taken from ). The result from a classical continuum calculation is shown for reference (solid blue line). (E) System geometry when the arginine (green) is positioned at the upper leaflet. This configuration represents the far right position on A–D. Gray surfaces represent the lipid head group–water interface, purple surfaces represent the lipid head group–hydrocarbon core interface. The membrane is not deformed in this instance. (F) System geometry when the arginine is at the center of the membrane. This configuration represents the far left position on A–D. The shape of the upper membrane–water interface (gray) was determined by solving .<p><b>Copyright information:</b></p><p>Taken from "A Continuum Method for Determining Membrane Protein Insertion Energies and the Problem of Charged Residues"</p><p></p><p>The Journal of General Physiology 2008;131(6):563-573.</p><p>Published online Jan 2008</p><p>PMCID:PMC2391250.</p><p></p

The influence of membrane bending on computing the biological hydrophobicity scale and the interplay of electrostatics and nonpolar forces

(A) Amino acid insertion energies for 17 residues calculated using our membrane-bending model (green bars) and compared with the translocon scale () (red bars) and a scale developed from MD simulations of lone amino acids () (blue bars). All three scales were shifted by a constant factor to set the insertion energy of alanine to zero (1.97 kcal/mol green, −0.11 kcal/mol red, 2.02 kcal/mol blue). The insertion energy of charged residues is reduced by 25–30 kcal/mol by permitting membrane bending. Calculations similar to those in indicate that only charged residues result in distorted membranes. (B) The energy difference between very polar amino acids and charged amino acids is quite small; however, our model predicts that the physical scenarios are quite different. The insertion penalty for asparagine is primarily electrostatic, while the nonpolar component stabilizes the amino acid (left bars). Conversely, there is very little electrostatic penalty for inserting arginine and most of the cost is associated with the nonpolar energy required to expose the TM domain to water (right bars). For comparison, we show that the classical flat membrane gives rise to a huge electrostatic penalty and a significant 6 kcal/mol membrane dipole potential penalty for arginine (middle bars).<p><b>Copyright information:</b></p><p>Taken from "A Continuum Method for Determining Membrane Protein Insertion Energies and the Problem of Charged Residues"</p><p></p><p>The Journal of General Physiology 2008;131(6):563-573.</p><p>Published online Jan 2008</p><p>PMCID:PMC2391250.</p><p></p

Determining the optimal membrane shape for a fixed charge in the core of the membrane

(A) Starting with the arginine at the center of the membrane as in , the helix was held fixed, and the point of contact of the upper membrane–water interface with the helix was varied from a height of 0 Å to its equilibrium width of 21 Å. The electrostatic energy decreases by 35 kcal/mol as the arginine residue gains access to the polar head groups and extracellular water. The majority of the decrease in electrostatic energy occurs from 21 to 5 Å, and there is pronounced flattening in the curve between 0 and 5 Å. (B) The total insertion energy exhibits a well-defined energy minimum at a contact membrane height of 5 Å.<p><b>Copyright information:</b></p><p>Taken from "A Continuum Method for Determining Membrane Protein Insertion Energies and the Problem of Charged Residues"</p><p></p><p>The Journal of General Physiology 2008;131(6):563-573.</p><p>Published online Jan 2008</p><p>PMCID:PMC2391250.</p><p></p

Computed biological hydrophobicity scale using a classical view of the membrane

(A) A model peptide with an arginine at the central position (green) is shown spanning a low-dielectric region representing the membrane. Helical segments with this geometry were used to produce the bar graph in B. (B) Insertion energies for 19 of the 20 naturally occurring amino acids. Amino acids are ordered according to the translocon scale (). All molecular drawings were rendered using VMD ().<p><b>Copyright information:</b></p><p>Taken from "A Continuum Method for Determining Membrane Protein Insertion Energies and the Problem of Charged Residues"</p><p></p><p>The Journal of General Physiology 2008;131(6):563-573.</p><p>Published online Jan 2008</p><p>PMCID:PMC2391250.</p><p></p

A screenshot of the user interface.

<p>Parameters pertaining to the calculation are entered in the field on the left, and the molecule and membrane can be viewed in the embedded Jmol viewer on the right. Pictured here is the membrane-embedded single transmembrane helix used for the calculations in CASE I.</p

States used to compute protein solvation energies.

<p>(A) The helix (orange) is pictured embedded in the membrane, which is delineated by the upper blue and lower gray lines. The membrane core between the two red lines is assigned a dielectric value  = 2. A headgroup region of 8 Å is indicated between the water and membrane core. Bulk water above and below the membrane is assigned a dielectric value of  = 80. (B) The helix in the bulk water ( = 80) in the absence of the membrane. The helix carries one charged residue (Arg14) shown in green in (A) and (B). The protein solvation energy is calculated by computing the total electrostatic energy of systems A and B and then calculating the quantity: . Images rendered with VMD <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0012722#pone.0012722-Humphrey1" target="_blank">[63]</a>.</p

Top view of the KcsA channel (green) and the = 2.01 isocontour highlighting the membrane interface (gray).

<p>The K ion in the center of the channel is shown in blue. (A) When the membrane is not excluded from the channel pore, we observe that membrane is added to the pore region. (B) With the exclusion radii set too high at 28 Å, there are large gaps of water between the outer membrane and the protein. (C) The channel should be clear of membrane and the membrane should fit snugly around the outside of the protein as shown here. Membrane exclusion radii are 24 Å and 16 Å at the top and bottom of the channel, respectively.</p

A cartoon representation of the distinct dielectric environments in each calculation.

<p>The orange regions represent protein, the gray membrane, and all white areas indicate water. The inner solution space at the bottom is assigned a voltage of , and correspondingly an effective charge density is assigned and a value of one for the variable . The water in the center of the channel is assigned values for and that correspond to the outer solution space. The lower value of the membrane (dashed line) separates the inner and outer solution spaces. In the gray region, and are set to 0 and , respectively, to mimic the membrane.</p

Parameters for protein solvation CASE I.

<p>Parameters for protein solvation CASE I.</p