From the gating charge response to pore domain movement: Initial motions of Kv1.2 dynamics under physiological voltage changes

Abstract

Recent structures of the potassium channel provide an essential beginning point for explaining how the pore is gated between open and closed conformations by changes in membrane voltage. Yet, the molecular details of this process and the connections to transmembrane gradients are not understood. To begin addressing how changes within a membrane environment lead to the channel’s ability to sense shifts in membrane voltage and to gate, we performed double-bilayer simulations of the Kv1.2 channel. These double-bilayer simulations enable us to simulate realistic voltage drops from resting potential conditions to depolarized conditions by changes in the bath conditions on each side of the bilayer. Our results show how the voltage sensor domain movement responds to differences in transmembrane potential. The initial voltage sensor domain movement, S4 in particular, is modulated by the gating charge response to changes in voltage and is initially stabilized by the lipid headgroups. We show this response is directly coupled to the initial stages of pore domain motion. Results presented here provide a molecular model for how the pre-gating process occurs in sequential steps: Gating charge response, movement and stabilization of the S4 voltage sensor domain, and movement near the base of the S5 region to close the pore domain.

Keywords::

• Voltage gating
• ionic solution
• double-layer effect
• S4 domain

Introduction

Transmembrane potentials occur based on differences of the electrical charge distributions on the intracellular and extracellular surfaces. In particular, K⁺ channels play a key role in maintaining this ionic balance across the membrane [1–7]. These channels undergo a voltage dependent gating transition when a change in transmembrane voltage produces a structural change within the channel to create a passage of potassium ions across the cellular membrane [7]. One of the most studied K⁺ channels is the Kv1.2 channel. The X-ray structure of this channel was solved by Long et al.[8]. The mammalian Kv1.2 channel is a member of the Shaker K⁺ channel family and is closely related in amino-acid sequence and electrophysiological properties to other K⁺ channels [8]. It consists of a homotetramer with each of the subunits containing six transmembrane segments, S1–S6, and a cytosolic T1 domain [8]. The S5–S6 region forms the ion-conducting pore and determines ion selectivity, whereas the S1–S4 region forms the highly charged voltage sensor domain which includes the ‘paddle domain’. The voltage sensor domain is thought to move in direct response to a change in voltage difference across the membrane. The voltage domain response is mediated by specifically charged amino acids called ‘gating charges’.

Potassium channels cycle between two states, a closed state under a resting potential (−90 mV to −50 mV) and an open state under depolarization (−5 mV to +40 mV), as the strength of the voltage gradient across the bilayer varies [7]. Studies of these channels and their relationship with differing voltage gradients date back to the initial voltage clamp experiments performed by Hodgkin and Huxley in 1952 [6]. In computational simulations, the difficulty in simulating conditions comparable to a voltage clamp experiment lies in determining the spatially and temporally averaged potential on a detailed atomic level. The molecular dynamics timescale is relatively brief compared to an experimental timescale, so the simulations to be discussed can only touch on the initial changes due to a single voltage step, rather than following the full voltage dependent response.

There are several proposed approaches for modeling a transmembrane potential that have been applied to explore the effect of voltage on channel function [9–14]. One method represents the transmembrane potential through the use of a constant electric field [9]. In this approach one creates a constant electric field by adding a force dependent on the partial charge’s z-position relative to the xy-bilayer plane and the magnitude of the voltage drop. In this way a time-averaged constant voltage can be maintained across the bilayer. Recently, a study of the Kv1.2 channel was performed using this method [11]. The study examined how the electric field within the membrane becomes ‘focused’ resulting from interactions of the gating charges with their surroundings. ‘Focusing’ in this context refers directly to these changes brought about by the charge distributions (fixed on the channel and mobile) that create a much larger voltage drop at one location relative to other locations. Another transmembrane potential method was introduced by Delemotte et al. to study the conductance of a nanotube embedded within a bilayer [10]. They considered a bilayer surrounded by electrolyte baths enclosed by an air/water interface slab. The potential is established by an asymmetric distribution of ions on both sides of the membrane. With the constant-field model or the model with the air/water slab, several important characteristics of membrane systems may still not be captured, such as the ionic surface associations and fluctuations of the potential away from the average.

In this study, we investigate fluctuations of the Kv1.2 channel and its environment under different transmembrane potentials. Our approach for simulating the channel under various voltage gradients uses the double-bilayer method, an approach that was first used for pure bilayer simulations [12,13]. For this method, the voltage drops are dependent solely on the explicit charge concentrations felt across the bilayer. The method accounts for a transmembrane potential by simulating more than one bilayer per unit cell thus establishing two baths separated by two bilayers (Figure 1); the total number of atoms between each simulation remains the same. In contrast with the constant-field method, the double-bilayer method has the ability to observe fluctuations and to determine local charge interactions within separate baths. While the air-water method has the ability to observe fluctuations similar to the double-bilayer, the contributions and interactions of the baths may not be entirely separate due to the effect of the air/vacuum slab across the periodic boundary. The advantage of our double-bilayer systems is that one can address both the fluctuations and the averages of the potential distributions.

Figure 1.  All-atom depiction of Kv1.2 channel double-bilayer system. The left and right boundaries are periodically connected. The K⁺ channels are represented as a purple C_α trace. The DOPC lipids (yellow) separate Bath 1 and Bath 2 displayed in pink. The chloride (green) and the potassium (blue) distribution within each bath are described in Table I. This Figure is reproduced in colour in Molecular Membrane Biology online.

Our K⁺ channel calculations, with nearly 813,000 atoms, are the largest systems constructed and simulated using this type of method to date [12–14]. This series of all-atom molecular dynamics (MD) simulations allow us to investigate the movement and interactions of the gating charges and voltage sensor domain behaviour as the transmembrane potential changes. Using these MD simulations in combination with Principal Component Analysis (PCA), we further investigate how motions of the S4 domain are coupled to fluctuations of the pore domain. We comment on the nature and degree of fluctuations for both the voltage sensor domain and the pore domain at different voltage drops across the bilayer (see Table I). Based on these fluctuations, we suggest several pre-gating steps leading up to the full Kv1.2 channel gating process.

Table I.  Bath charge concentrations of the Kv1.2 channel.

System	Bath 2	Bath 1	Voltage drop
RP	694 K^+, 684 Cl^-	706 K^+, 628 Cl^-	−56 mV
DP1	695 K^+, 683 Cl^-	705 K^+, 629 Cl^-	−43 mV
DP2	696 K^+, 682 Cl^-	704 K^+, 630 Cl^-	−30 mV
DP3	697 K^+, 681 Cl^-	703 K^+, 631 Cl^-	−10 mV
DP4	698 K^+, 680 Cl^-	702 K^+, 632 Cl^-	+6 mV
DP5	699 K^+, 679 Cl^-	701 K^+, 633 Cl^-	+15 mV
Protein	147 Arg, 56 Lys, 173 Glu, 48 Asp	53 Arg, 24 Lys, 75 Glu, 72 Asp

RP = resting potential, DP = depolarized potential (The ‘protein’ contribution should be included for each system when summing each bath’s charge distribution as the residue contributions remained consistent for each bath).

Methods

Channel modeling

The structure of the Kv1.2 channel was solved using X-ray crystallography to a resolution of 2.90 Å (PDB ID: 2a79) [8]. The structure represents residues 32 through 421 of the entire Kv1.2 channel and has some missing loop regions. A structural alignment using ClustalW [15] was used to determine the best position of the sequence in comparison with the crystallographically determined regions of the channel, thus filling in the missing loop and sidechain segments and completing the full structure of the channel. The ClustalW results indicated that residues 203–218, 244–254, 276–287 would be loop regions, residues 255 through 275 (the UNK C crystallographic region) aligned with the S3 region, and residues 133 through 202 (the UNK D crystallographic region) aligned with the S1 and T1-S1 linker regions. To confirm the secondary structure of residues 32 through 421 of the entire Kv1.2 channel, the residue sequence was submitted to TMHMM Server v2.0 program [16]. Our secondary transmembrane structure prediction results were validated for that region of the Kv1.2 sequence (See Supplementary Figure 1 in Supplementary information, online version only).

To construct coordinates for the missing regions and place the residue sidechains of the Kv1.2 channel, we used the CHARMM program (version c32b1) with CHARMM27 parameters [17,18]. The titratable residues were simulated using their default charged state in the CHARMM program, where Arg and Lys had a charge of +1, Glu and Asp had a charge of −1, and His had a neutral charge. After the loop coordinates were constructed, a series of short MD simulations were run with constrained non-loop coordinates under GBSW implicit membrane conditions [19–22] to sample for various loop orientations. The GBSW membrane parameters used 30 Å for the thickness of the membrane hydrophobic core, 5 Å for each membrane interface, and a surface tension coefficient of 0.04 kcal/(mol*Ų) to represent the non-polar solvation energy. The model represents a total membrane thickness of 40 Å to mimic a bilayer of dioleoylphosphatidylcholine (DOPC). The simulations were run using a 1 fs timestep at a temperature of 300K for 2 ns. The final loop configuration for each monomer was selected based on the loop configuration that exhibited the lowest energy value.

The most recent Kv1.2-Kv2.1 chimeric structure (PDB ID:2R9R) was released after the completion of our simulations. This replaces the Kv1.2 helix-turn-helix domain (S3b–S4) with that of the Kv2.1 channel [47]. The loops constructed on our Kv1.2 channel do differ from this structure, but not so dramatically that they change channel behaviour. This is because the configurations of the functionally relevant regions of the Kv1.2 channel (2a79) are similar to that of the chimeric channel (2r9r). The backbone root means square difference (RMSD) of the functionally relevant regions (S1, S2, S3, S4, S4–S5 linker, S5, P, and S6) between the two structures is 0.98 Å regions (Supplementary Figure 8, online version only). The full backbone RMS difference between the structures is 5.7 Å, stemming from the differences in loop configurations and a difference in the T1–S1 linker domain (Supplementary Figure 8, online version only). This confirms that our initial construction of loop conformations would not significantly affect the overall channel structure.

Ensemble preparation

We modeled the Kv1.2 channel (with ions in the pore region) in an all-atom double-bilayer system to understand how the channel responds to a transmembrane electrostatic potential and to enhance sampling efficiency (Figure 1). An initial double-bilayer system was constructed using dioleoylphosphatidylcholine (DOPC) lipids. The dimensions of the unit cell was 150 ×150 × 340 Å. The DOPC bilayers each contained a channel and was spaced 137 Å apart to avoid direct surface interactions and interdependence between the two bilayers or channels. The resting-potential (RP = −56 mV) double-bilayer system was created by replacing a water molecule with either a K⁺ or Cl^– ion until the starting concentration was ∼ 100 mM for each bath. The placement of the ions was based on energetically favorable locations for a RP voltage drop across the bilayer. The positions of two K⁺ located within the selectivity filter domain were preserved. The selectivity filter positions one and three, (S₁ and S₃), were used to maintain a conductive pore conformation throughout the simulations. These two ions were initially fixed only during the equilibration of the starting system to prevent collapse of the selectively filter domain. During the double-bilayer production simulations, these ions in the selectivity filter fluctuated within the S₁ and S₃ filter positions by tenths of an Angstrom (Supplementary Figure 11, online version only). The initial system was equilibrated for 11 ns with constant lipid surface area conditions and a Nosé-Hoover thermostat to hold the simulation at a temperature and pressure near 300 K and 1 atm.

Based on the equilibrated RP system, five different depolarized systems (DP) were created by the exchange of ions located in opposite baths. The ion exchange between the two baths facilitated the creation of a transmembrane electrified potential. The net potential for each system can be calculated based on the charge distributions listed in Table I. The potentials used in this paper differ from previous membrane dipole potentials: The applied transmembrane potentials are smaller and remain within the physiological voltage potential conditions. The resting potential system and the five depolarized systems each contained a total of 812,614 atoms.

MD simulations

All-atom molecular dynamics (MD) simulations of all constructed ensembles were performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) version June 2007 and CHARMM27 force field parameters [18,23,24]. LAMMPS has been used for previous biomolecular transmembrane potential simulations with great success, having excellent parallel scaling abilities [13,14]. The simulations of the six different systems were performed at Sandia National Laboratories on the large-scale Computational Plant (Cplant) cluster. An average of 0.8 ns per day on 250 processors was completed.

The unconstrained MD simulations were performed for an equilibration time of 1–2 ns and a simulation run time of 10 ns at a temperature of 300 K for each of the simulations. These unconstrained MD simulations allows for flexibility of the sidechain and backbone of the protein and a free-flow movement for the lipid, water, and ion molecules. The SHAKE algorithm held the bonds of the hydrogen atoms rigid to enable a 2 femtosecond timestep. The TIP3 waters were also held rigid by SHAKE. Each system was simulated under constant lipid surface area conditions and allowing only the simulation cell length in the z-direction to fluctuate in order to maintain isobaric conditions. A Nosé-Hoover thermostat held the simulations at a temperature and pressure near 300 K and 1 atm, values that are consistent with CHARMM lipid force fields. The z-direction box length remained around 340 Å, with minor fluctuations of about 1–1.5 Å, for both the RP and the DP systems. The final dimensions for the periodic boundary conditions of each system were 150 Å × 150 Å × 340 Å. The separation between the bilayers had minor fluctuations of about 1 to 1.5 Å throughout the 10 ns. The particle-particle-particle-mesh (P³M) method was used to compute the long-range electrostatics, which is similar to the commonly used particle mesh Ewald (PME) method. We used a real-space cut-off of 10 Å and smoothly switched the van der Waals (vdW) interaction to zero between 8 and 10 Å.

While a total simulation time of 10 ns is short compared to the timescale of the full gating-transition (microseconds to milliseconds), the use of an ensemble of double-bilayer simulations at different transmembrane potentials gives us some insights into the initial stages of the gating process. At the same time, we want to note that any large molecular dynamics simulation would have sampling problems for those events that are not observed in the simulations. We thus report on averages and analysis of rare events (observed once or twice in a simulation) and attempt to represent our confidence with error bars and discussion in what follows.

Analysis

The analysis for the all-atom calculations was predominately completed using MDAnalysis scripts, an analysis library developed by Michaud-Agrawal [25]. The potential profiles were calculated based on the average of the charge distribution over time and space for particular one-dimensional z-slices. This calculation is similar to the calculation performed by Lee et al., Sachs et al. and Denning and Woolf [12–14]. The phosphate headgroup distance and angle cut-off values used to assess our lipid hydrogen-bond patterns were r_HB = 0.35 nm and α_HB = 60° and have been used to evaluate hydrogen-bond patterns in other studies [14,26–28]. For the number of coordinating water molecules, distance cutoff values of r_HB = 0.35 nm, 0.40 nm, and 0.50 nm were used, similar choices have been made in other studies [14,26–28].

To determine whether equilibrium conformational fluctuations of the pore vary with the voltage drop across the bilayer, additional analysis was performed. Principal component analysis (PCA) determines the effective normal modes observed in the MD trajectory. The method we use is similar to that used to reveal key fluctuations of the NikR protein by Bradley et al. [29]. PCA for α carbons of the channel pore domain used the quasiharmonic vibran module in CHARMM. The vector (Δx) described the displacement of α carbon atoms from the minimized starting open to a final closed structure determined by X-ray crystallography. The Kv1.2 open and closed structures were created through homology modelling using Modeller and threading the Kv1.2 sequence through the K⁺ channel pore X-ray structures. The KcsA structure (1k4c) represented the closed conformation and the Kv1.2 (2a79) structure represented the open conformation. The question of whether the displacement vector from the X-ray structures, Δx, represented in the vector space described by subsets of PCA modes change as a function of voltage was accomplished by calculating the weight factor (α⁽ⁱ⁾) for each PCA mode as follows:

Then a subset S of modes was used to calculate a reconstructed vector () as follows:

The similarity between Δx and each PCA mode (v_i or ) was assessed by their overlap as measured by the angle between the two vectors:

The relative correspondence of the reconstructed vectors was computed:

This relative correspondence provides a measure of how well the displacement due to the reconstructed vector agrees with the observed crystallographic displacement vector. This calculation assesses conformational fluctuations when compared across different voltages and the number of modes required to represent these motions. For each of the different voltages, the relative correspondence, ϵ, as a function of the mode number, x, was found to decay exponentially, i.e., ϵ(x) = exp(-m*x), and thus can be characterized by a single parameter m (or the decay constant 1/m). We calculated m from a linear fit to y = ln(ϵ(x)) via a least-squares linear regression:

where m is the slope, n represents the number of data points (180), x is the mode number, and y is the natural logarithm of the relative correspondence at a particular x value. This fit for the relative correspondence allows us to understand how the conformational fluctuations changed as a function of voltage, averaged over the full trajectory.

Results

Voltage drop profile

The potassium channel undergoes several pre-gating steps in preparation to change conformation between a closed state under a resting potential (−90 mV to −50 mV) and an open state under depolarization (−5 mV to +40 mV) [7]. To simulate physiological changes in voltage and explore steps key to the K⁺ channel gating process, the Kv1.2 channel was simulated using the double-bilayer method. Figure 1 shows the nature of the double-bilayer method. The use of two independent baths allows the explicit variation of salt concentrations between the inner and outer baths which establishes a biological voltage gradient. Placing a single channel in each bilayer supports thorough sampling for changes in channel behavior that occur under the same applied transmembrane potential conditions.

The transmembrane potentials used here reflect physiological conditions necessary to stimulate the channel gating process. Five different voltage systems were used to simulate depolarized states, relative to the resting potential state, for the K⁺ channel. Each of our simulations represents a different physiological transmembrane potential gradient associated with a distinctive voltage profile. The voltage drop across a bilayer for each system established itself over a time period of 1 ns (Figure 2A). The smoothness of the 1 ns voltage profile indicates that the strength of the ionic concentration and the mobile charge distribution is evenly distributed within the two baths for each simulation. While this is not a proof that the voltage drop is constant over each simulation, this suggests that the long-time and spatial average produced within each simulation is that of a Nernst potential determined by the salt baths. The advantage of our double-bilayer systems is that one can address both the fluctuations and the averages of the potential distributions. This benefit allows one to examine questions of how the local charge distributions can contribute to the conformational changes discussed in later sections.

Figure 2.  Establishing K⁺ channel potentials. (A) Transmembrane potential profiles at time specific averages. (Dotted line = 100 ps; dot-dashed = 1 ns; solid = 10 ns). Arrows indicate locations of bilayer. (B) Transmembrane potential profiles. The line colour corresponds to the expected voltage gradient for each system (RP = resting potential, DP = depolarized potential). The arrows along the axis indicate the location of the baths and the bilayer regions. The simulations contain a constant total number of ions. To establish the different voltage gradients ions were exchanged between the two baths. (C) The two-dimensional (xz-plane) transmembrane potential profile of the DP5 system has a 15 mV drop and shows the inhomogeneous nature of the voltage drop across the bilayer (see arrows). This Figure is reproduced in colour in Molecular Membrane Biology online.

The six different voltage profiles, seen in Figure 2B, show how the voltage drop across the bilayers varies as the ionic concentrations within the two baths is changed. These simulations represent the shift from a resting potential (−56 mV) to a depolarized potential (+15 mV). The local charge contributions made from the channel and bilayer distributions are located where the local voltage has the highest values (Figure 2B). These peaks can be attributed to the contributions of the partial charges of the channel and lipids. To further examine the charge distribution details found within the peak regions, a two-dimensional (2D) voltage profile was calculated. Similar to the one-dimensional voltage profile, the 2D voltage profile of the xz-plane establishes itself over a time period of 1 ns (Supplementary Figure 6, online version only). Figure 2C shows a time-averaged 2D voltage profile of the xz-plane for the DP5 system. This Figure shows how the electric field is focused near the S4 domain. Similar behaviour has also been observed in other studies and it has been hypothesized that such a focusing effect contributes to an increased switching sensitivity of the channel as the voltage drop changes [27,28,30,31].

This sensitivity near the S4 domain is a result of the different interaction patterns of the gating charges with their environment as the transmembrane potential changes. Figure 2B shows that within the channel/bilayer regions the peaks shift toward Bath 1, the extracellular bath, by ∼ 2–3 Å as the voltage drop across the membrane becomes depolarized. The shifting of these peaks is due to the re-adjustment of the gating charges within a membrane environment. In other words, as the S4 helix shifts toward the extracellular bath the individual gating charges become more solvated thus shifting the electric field further toward the remaining charges buried within the bilayer (solvation of these charges are later discussed in more detail). In the next section, we show how gating charge interactions with lipid and water molecules in a membrane environment lead to earlier steps of the gating process.

Gating charge interactions

The lipid and water molecule interactions are important for understanding how charged residues behave within a bilayer environment [32–34]. In particular, the atomic interactions of the gating charges, Arg²⁹⁴, Arg²⁹⁷, Arg³⁰⁰, and Arg³⁰³ located on the S4 domain, are of interest because their movements are considered to be the driving force for voltage sensor domain motions [11,14,27,28,30–41]. As no constraints were used during these simulations, the arginine positions are influenced by their interactions with lipid and water molecules and thus contribute to changes within the local electrical field of the membrane [11,14,27,28,30,31]. Note that our understanding of their response to voltage depends on convergence of sampling and thus we are limited to what is seen within the 10 ns simulations at each voltage level. Previous studies of smaller systems, involving only the helix-turn-helix motif of the voltage sensor domain, observed fluctuations in the total number of lipid and water hydrogen-bonds around the gating charges for all proteins when varying the voltage from the resting potential to a depolarized potential [14,30]. Due to the unconstrained nature of our MD simulations, environmental effects such as lipid or water molecule coordination influence arginine sidechain behavior (e.g., flexibility and position within the bilayer). A full understanding of these effects would clearly require several orders of magnitude greater simulation time to understand how frequently different environmental changes occur. In the much shorter time of our simulations, we examined the fluctuations in the total number of lipid hydrogen-bonds and coordinating water molecules with these gating charges using the cutoff values noted in the Methods section. As can be seen in Supplementary Figures 4 and 5 (online version only), the initial fluctuations in the number of lipid headgroups and water molecules with each subunit and each residue are subtle and no major change is detected during the 10 ns of our simulations. To determine how the collective environmental effects are distributed across a single channel as the transmembrane voltage changes, the average number of lipid headgroups and water molecules per subunit were calculated for each channel. These average changes reflect the behavior over a 10 ns time interval and should not be interpreted as changes that occur immediately upon a voltage change.

Our averaged-results show that in the initial stages of gating the number of lipid headgroups in close contact with Arg²⁹⁴ and Arg²⁹⁷ is voltage-dependent (Figure 3). The charged residues coordinate differently with lipids depending on the strength of the applied potential. As the potential becomes more depolarized the number of lipids that form hydrogen-bonds with Arg²⁹⁴ and Arg²⁹⁷ increased, while hydrogen-bonding with Arg³⁰⁰ and Arg³⁰³ remained negligible. This lipid response only occurred for a single channel while the number of lipid hydrogen-bonds for the other channel with the system remained constant independent of changes in voltage. This indicates that the timescales for initial formation of lipid-protein hydrogen bonds to occur may be on the order of nanoseconds.

Figure 3.  The lipid hydrogen bond pattern for the individual channels as a function of voltage. (The square dot represents Arg²⁹⁴; triangle dot represents Arg²⁹⁷; diamond dot represents Arg³⁰⁰; circle dot represents Arg³⁰³) The error bars represent the standard deviation. (A) Channel 1 gating charges and the number of lipid hydrogen bonds. (B) Channel 2 gating charges and the number of lipid hydrogen bonds.

While the number of lipid hydrogen-bonds in our simulations changed relative to the transmembrane potential, the number of water molecules coordinating the gating charges did not change appreciably in the early stages of pre-gating (Figure 4). The number of coordinating water molecules for each of the gating charges was calculated at different cut-off distances r_HB (Figure 4). Arg³⁰⁰ and Arg³⁰³ had an average of 3–5 water molecules coordinating within an r_HB = 3 Å, while an average of 6–8 water molecules were coordinating Arg²⁹⁴ and Arg²⁹⁷. The number of water molecules surrounding the gating charges increased with r_HB of the gating charges. While it is not observed in this study, previous work, in smaller systems, has found that the total number of water molecules around the arginines changed after 10–30 ns when the voltage was varied from −190 mV to +100 mV [14,30]. Due to the short timescales of our simulations, one would expect that on longer timescales the number of water molecules coordinating the gating charges would exhibit a voltage dependence based on the results observed in these previous studies. Both the lipid and the water hydrogen-bonding effects produce changes within the local electrical field of the membrane that influence the position of the gating charges within the membrane and thus cause a shift in the position of the S4 voltage sensor domain.

Figure 4.  Water hydrogen bonding pattern for the gating charges. The distributions show the average number of coordinating water molecules (Nij) as a function of the voltage gradient for 3 of the gating charges at different radial distribution values. (The square dot represents channel 1; the triangle dot represents channel 2). The error bars represent the standard deviation.

Gating charge fluctuations

The gating charge’s positions are influenced by their interactions with lipid and water molecules and shift as the local electrical field of the membrane changes [11,14,28,30]. In the open state or under depolarized conditions, previous studies have shown that Arg²⁹⁴ and Arg²⁹⁷ interact more strongly with both lipid headgroups and water molecules shown by the hydrogen-bond partner analysis [14,30]. Arg³⁰⁰ and Arg³⁰³ were also shown to form salt-bridges with Glu¹⁸³ from the S1 domain as the transmembrane potential is depolarized [11,28,30,38,44]; Arg³⁰⁰ would compete with Arg³⁰³ to form a salt-bridge interaction with Glu²²⁶ of the S2 domain [11,28,30,39].

In our simulations, the shift of the S4 domain in the z-direction (Figure 6A) repositions the first two gating charges, Arg²⁹⁴ and Arg²⁹⁷. Initially, we examine the centre of mass (COM) z-position for each of the four gating charges relative to the COM of the pore domain over the full 10 ns of our simulations. As seen in Supplementary Figures 3 and 10 (online version only), no collective shifts among the charged residues were detected for either channel. This suggests that a longer timescale is needed to observe the full set of changes that underlie gating, a collective shift of all gating charges within the bilayer. The fluctuations for each residue were averaged among individual channels to determine if a subtle trend in the z-position shift for a particular channel occurs as the voltage drop changes. As the voltage changes from a resting potential to a more depolarized potential, they move from pointing toward the S1–S3 domains to pointing towards the extracellular bath region (Figure 5). In addition to Arg²⁹⁴ and Arg²⁹⁷ moving in response to changes in voltage, Arg³⁰³ also showed a response that correlates to a change in the voltage gradient (Figure 5). These initial events suggest these three residues contribute to movements of the voltage sensor domain during the early stages of the gating process.

Figure 5.  The centre of mass position of the sidechain for the gating charges for each channel versus the voltage drop. The error bars represent the time-averaged standard deviation for each z-position. (The dashed line represents Arg²⁹⁴; the solid line represents Arg²⁹⁷; the dash-dot line represents Arg³⁰⁰; and a dot-dot line represents Arg³⁰³) (A) Channel 1 (B) Channel 2.

Figure 6.  The centre of mass position of the voltage sensor domain for each channel versus the voltage drop. The error bars represent the time-averaged standard deviation for each z-position. (The solid line represent channel 1; the dashed line represents channel 2) (A) S4 domain, (B) S2 domain, (C) S3 domain.

Voltage sensor domain fluctuations

As the voltage sensor domain responds to differences in the transmembrane voltage drop, the gating effect can be seen as a type of domain motion. Previous experimental studies have examined the degree and nature of this domain motion [8,31,38–43] as well as interactions within a bilayer environment [35–37]. The S4 domain motion, in particular, is influenced directly by the response of the gating charges. The S4 domain has been shown to move laterally within the bilayer as gating charges move and the strength of the voltage drop changes [14,30,40,43,44]. For example, previous studies have suggested that S4 domain translation within the bilayer is subtle moving only 2–3 Å [11,44,45]. Upon analyzing the z-position for the COM for individual S4 domains relative to the COM of the pore domain, no collective shifts among four S4 domains for a single channel were observed over our simulation period of 10 ns (Supplementary Figures 2 and 10, online version only). To determine if our fluctuations for individual S4 subunits exhibit a general trend for a particular channel, the average S4 COM position was calculated for both channels. When taking the average COM positions of all S4 domains for each channel, a shift in position along the membrane normal is observed as the transmembrane potential decreases. Under depolarization, the S4 domains move towards the extracellular side by 1–2 Å for this channel (Figure 6A). The degree and the direction of the S4 domain shift is similar to previous computational simulations on longer timescales [44,45]. This shift may begin to occur on a timescale of 10 ns as the S4 domain shift only occurred for one of the channels. The observed S4 domain fluctuations and resulting shift are coupled to the degree of pore domain fluctuation.

The motions of the S4 domain, while subtle, do not occur without additional supporting interactions. Previous studies have shown that both lipid headgroups and small S2 and S3 domain motions may be an important part in stabilizing the motions of the S4 domain [11,28,30,35–40,42,43]. On the timescale of our transmembrane voltage simulations, only small COM fluctuations of the S2 and S3 domains were observed (Figures 6B, 6C). The small initial fluctuation of S2 and S3 indicate that these domains do not respond directly to changes in voltage during the steps of the pre-gating process. While this is not a direct assessment that the S2 and S3 domains stabilize the S4 domain, this small motion is consistent with other studies of longer timescales [30,44,45] and may be an initial indication that their small movement contributes to the stabilization of the S4 domain movement.

PCA: Pore fluctuations

As the S4 domain shifts during the gating process, its motions drive the structural transition between the open and closed state of the pore domain by influencing the movement of S5 and S6. Our double-bilayer simulations only explore a small portion of the transition, the pre-gating effects; however, even on this small timescale fluctuations could be present that would eventually lead to the large-scale gating transition motions. To explore the context of the structural transition of the Kv1.2 channel pore domain, Principal Component Analysis (PCA) was used. PCA provides an additional avenue to quantitatively determine if the all-atom equilibrium fluctuations carry within them low frequency modes necessary to complete a specific transition. Since our simulation timescale is short compared to the full K⁺ channel transition, we examine the extent of pore fluctuations to determine which modes are necessary for the transition between the open and closed state, and thus the amount and direction of the motion coupled to different voltage drops across the bilayer.

The initial K⁺ channel transition between the conductive and non-conductive states can be observed through the use of reconstructed PCA modes and a displacement vector. PCA was performed for the α carbons of 404 pore domain residues from the Kv1.2 equilibrium MD simulations yielding a total of 1212 modes based on 3N mode space, N being the total number of atoms. The first 600 PCA modes out of 1212 modes were analyzed and compared to standards used in previous PCA/MD simulations studies [29].

The question becomes what percentage of the total movement represents similarity between Δx and each PCA mode along the open to closed path is captured within these 600 modes. To do this, the relative correspondence, also known as relative error, of subsets of modes was calculated along a displacement vector for the pore domain atoms. The relative correspondence provides an indication of how subsets of modes describe motions that would lead to the closing of the pore domain. Figure 7 shows the relative overlap for a particular subset of modes with the open to closed displacement vector for each voltage drop. A low relative correspondence indicates a higher percentage of the included modes contribute to the motions closing the channel. As the total number of modes increases, the amount of relevant motion closing the channel increases (Figure 7). Figure 7A shows that 180 modes out of the 600 modes analyzed capture most of the movement for each of the individual voltage differences. These 180 modes provided a reasonable representation of the conformational change present at each voltage based on standards of PCA/MD simulations studies [29] and this subset of modes captured over 90% of the closing motions of the pore domain (Figure 7). Since the relative correspondence was found to decay exponentially as a function of the mode number x, Figure 7A was re-plotted (Figure 7B) as a logarithmic plot to show how the natural logarithm of the relative correspondence changes as a function of the mode number. Based on the logarithmic plot, one can see that the natural logarithm of the relative correspondence changes linearly as the number of modes increases within the 180 modes.

Figure 7.  The percent of displacement overlap for the pore domain as a function of mode for each voltage drop. (Yellow = RP; Pink = DP1; Cyan = DP2; Red = DP3; Green = DP4; Blue = DP5). The displacement vector represents an open ◊ closed channel vector. (A) Basic line plot (B) Natural logarithm plot (C) First 180 modes basic line plot (D) First 180 modes natural logarithm plot. This Figure is reproduced in colour in Molecular Membrane Biology online.

While the natural logarithm relative correspondence seems to follow a linear fit, it is difficult to ascertain how these fluctuations differ between the six double-bilayer simulations. The subtle difference between the correspondence profiles for the six systems may indicate a low signal-to-noise ratio. This result though is not surprising because we only sample on a small portion of the full gating transition during our MD simulations. To further analyze the correspondence profiles, we examine how the natural logarithm of the relative correspondence changes as the subset of included modes increases by calculating a linear fit for the single parameter m (or the decay constant 1/m) for each of the different voltage simulations. This calculation allows us to quantitatively analyze the relative correspondence and thus the overall pore domain fluctuations to determine if they would lead to the closing of the pore.

The rate of movement along the subset of reconstructed vectors differs (Figure 8). It changes with the voltage drop across the bilayer and its rate of change, the logarithm of the percent correspondence per mode, is shown in Figure 8. A lower rate of change indicates that the subset of modes involved move more along the open to closed displacement vector. In other words, the lower the rate of change the more the modes at a specific voltage contribute to the global movement toward a more closed state. For example, at −60 mV the rate of change shows the pore domain motions are closing the channel, while at −10 to 15 mV the rate of change values are higher and the modes do not sample along the displacement vector. This result makes sense because the channel is already and should remain in an open state under depolarized conditions (−10 to 15 mV).

Figure 8.  Characterization of motions and displacement based on the PCA modes of the pore domain. (A) The Figure shows the rate of change in the error of overlap as a function of voltage. (represents open ◊ closed). The error bars indicate the error based on subsets of the trajectory data. (B, C) The colours indicate the degree and direction of displacement (Blue = high overlap with Δx; White = low overlap with Δx; Red = minimal overall motion) (B) C_α trace of the pore showing the motions at −56 mV. (C) C_α trace of the pore showing the motions at +15 mV. This Figure is reproduced in colour in Molecular Membrane Biology online.

The depolarized potential simulations exhibit less movement along the open to closed displacement vectors. The movement for this simulation is largely around the loop regions of the pore (Figure 8B). The resting potential simulations, ranging from −56 mV to −30 mV, display greater movement along the displacement vector and movement is localized near the connection of the S5 to the S4–S5 linker (Figure 8C). At depolarizing potentials we observe more fluctuations of the S4–S5 linker and these fluctuations are preferentially directed in the direction towards the closed state structure. This is consistent with simulations of longer timescales [45] where the motions near the S4–S5 linker and S5 region correspond to the opening and closing of the pore domain. Although we do not observe a full gating transition, this result indicates that the pore-closing motions are encoded in the open state structure. Furthermore, these motions are enhanced at transmembrane potentials that lead to channel closure under physiological conditions and eventually facilitate the gating transition at longer time scales [46]. Thus motions of the pore domain are shown to be directly associated with motions of the S4 domain.

Discussion

The focus of K⁺ channel experimental studies has been on understanding the full K⁺ channel gating transition and how gating charge motion and interactions are coupled to the gating of the pore domain in response to changes in the transmembrane potential [6,35–40,42,43]. Computationally, it is difficult to show the full K⁺ channel gating-transition on the timescales of an unbiased all-atom MD simulation as the full transitions occurs on a microsecond to millisecond timescale [1–5]. An unbiased MD simulation would allow one to examine the natural course of interactions and events of a transition without the presence of additional forces. Previous and current results shown in this paper illustrate that full rigid body motion of the S4 sensor alone in unbiased MD simulations begins on the order of ns-μs [14,27,28,30,44].

We hypothesize that the gating transition process occurs in steps: gating charges are influenced by changes in their environment, S4 moves laterally in response, and then this change is transmitted along the S4–S5 linker biasing the pore domain to open or close the channel. The MD double-bilayer simulations together with PCA help to investigate how changes in the voltage drop across the bilayer lead to different motions for the Kv1.2 channel. Our results suggest that the dynamics seen in these simulations contain initial fluctuations relevant to the conformational transition between the conductive and nonconductive states.

The first ‘pre-gating’ phase involves the gating charges of the S4 domain interacting with their environment. The biophysical underpinnings of the local voltage changes around the S4 domain are believed to stem from the ability of the gating charges to interact differently with their environment under various transmembrane potential conditions [14,30]. The suggestion is that the fluctuations of these gating charges can signal the initial pore domain fluctuations as the voltage drop changes from a depolarized to resting-state potential and thus their interactions with their environment influence the degree of pore fluctuations. Previous electron paramagnetic resonance (EPR) measurements of the K⁺ channel show that only residues Arg²⁹⁴ and Arg²⁹⁷ would be exposed to the lipid headgroup region [38]. The EPR study showed that for KvAP the gating charges Arg¹³³ and Arg¹³⁶ were not exposed to the lipid headgroup environment based on their low probe mobility and their inaccessibility to either O₂ or NiEdda [38]. Our results are comparable to the EPR results in that Arg³⁰⁰ and Arg³⁰³ of the Kv1.2 channel, the equivalent of residues Arg¹³³ and Arg¹³⁶ of KvAP channel, do not hydrogen-bond with the lipid headgroups. This result did not change even as the strength of the applied voltage varied. Our results suggest that in the initial stages of gating the number of lipid head groups in close contact with Arg²⁹⁴ and Arg²⁹⁷ is voltage-dependent. The charged residues coordinate differently with lipids depending on the strength of the applied potential.

Water crevices surrounding the voltage sensor domain have been shown to direct the transmembrane electric field around the S4 domain, which leads to an energetic stabilization of the charged domain inside the membrane. It exposes the sensor to strong electrostatic forces so that one would expect the largest motions due to the field to occur at the S4 domain, thus creating a stronger effect within the electric field around the S4 domain [14,27,28]. While this was not analyzed in depth, similar water crevices were observed during our studies (Supplementary Figure 7, online version only). These water crevices are necessary to stabilize and solvate the arginines within the hydrophobic environment of bilayer. Arginines have been found to participate in stable salt-bridges with negatively charged residues from the neighbouring S1 and S2 helices, forming a structural barrier separating the extracellular and intracellular water cavities within the voltage-sensing domains [27,28].

Under physiological conditions, the interactions of the gating charges shift causing a lateral displacement of the S4 domain as the transmembrane potential changes. Experimental studies show the S4 domain moves laterally and rotates within the bilayer plane during gating while the S2 and S3 domains undergo little or no movement [30,31,36,37,42,43]. These fluorescence resonance energy transfer (FRET) and EPR studies showed no variation in signal for these two domains as the potential drop changed indicative of little to no large-scale movement. On the timescale of our transmembrane voltage simulations, the small COM fluctuation of S2 and S3 indicate that these domains do not immediately respond to changes in voltage. This is consistent with other computational studies [11,30,44] and may suggest that their subtle domain motions stabilize the motions of the S4 domain as it rotates within the xy-plane.

The degree and direction of S4 domain shift observed in this study is also similar to the previous experimental and computational simulations on longer timescales [31–40,44,45]. As the voltage drop shifts to a resting-potential, the S4 domain translates in the z-direction within the bilayer by 1–2 Å. The S4 displacement is directly connected to several aspects of the system. This shift affects the local electric field of the bilayer/channel regions and is reflected in a 1–2 Å shift in the voltage profile as the transmembrane potential changes. The S4 domain translation also corresponds to different degrees of pore domain fluctuations.

The shift of the S4 domain is coupled to pore domain fluctuations involved with the closing of the channel gate. As the potential approaches a resting-state value, the S4 domain translates in the z-dimension within the bilayer by 1–2 Å, producing increased fluctuations of the pore domain. The fluctuations of the pore domain increase due to the channel changing conformations to a non-conductive state as the voltage potential relaxes from a depolarized potential to a resting-state potential. In other words, as the strength of the voltage drop decreases so does the degree to which the pore domain fluctuates. The increased pore domain movement occurs near the beginning of the S5 domain. These initial fluctuations may be connected to those interactions between the S4–S5 linker and the C-terminus of the S6 helix responsible for closing the channel observed in other studies [45,46]. However, it is difficult to ascertain a direct connection, as the timescales for these simulations do not extend to the full gating transition. Although the entire transition is not observed these results do indicate that pore-closing motions are encoded in the open state structure and are enhanced during the TM potentials that lead to channel closure under physiological conditions, thus eventually facilitating the gating transition at longer time scales. The coupling between the fluctuations of the pore domain and S4 domain may indicate these are initial pre-gating motions necessary for transitioning between the conductive and non-conductive states.

In summary, the Kv1.2 channel undergoes several steps within the gating process. The results presented here provide evidence for bridging the gap between the initial gating charge response and the motions of the pore domain. Using a combination of MD and PCA, our results suggest how variations in transmembrane potential influence the gating charges and the voltage sensor domain and thus increase pore domain fluctuations. The gating charges are sensitive to changes in the transmembrane potential and it is reflected in their hydrogen-bonding patterns and position. The observed S4 gating domain motions may be indicative that its rigid body motions begin on the order of nanoseconds. The changes in voltage sensor domain position, in particular S4, are associated with pore domain fluctuations that open or close the channel. PCA analysis allows one to provide insights into motions typically observed in longer timescale simulations. This Kv1.2 channel study also contributes to the use of different methods used to study large transmembrane proteins under different voltage gradients and demonstrates the importance of understanding channel behavior and function under physiological conditions.

Acknowledgements

The authors thank Oliver Beckstein and Naveen Michaud-Agrawal for valuable discussions and contributions to MDAnalysis. This research was funded by a grant from National Institutes of Health (GM064746). Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.

Declaration of interest: The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.