The Free Energy Barrier for Arginine Gating Charge Translation Is Altered by Mutations in the Voltage Sensor Domain

The gating of voltage-gated ion channels is controlled by the arginine-rich S4 helix of the voltage-sensor domain moving in response to an external potential. Recent studies have suggested that S4 moves in three to four steps to open the conducting pore, thus visiting several intermediate conformations during gating. However, the exact conformational changes are not known in detail. For instance, it has been suggested that there is a local rotation in the helix corresponding to short segments of a 3-helix moving along S4 during opening and closing. Here, we have explored the energetics of the transition between the fully open state (based on the X-ray structure) and the first intermediate state towards channel closing (C), modeled from experimental constraints. We show that conformations within 3 Å of the X-ray structure are obtained in simulations starting from the C model, and directly observe the previously suggested sliding 3-helix region in S4. Through systematic free energy calculations, we show that the C state is a stable intermediate conformation and determine free energy profiles for moving between the states without constraints. Mutations indicate several residues in a narrow hydrophobic band in the voltage sensor contribute to the barrier between the open and C states, with F233 in the S2 helix having the largest influence. Substitution for smaller amino acids reduces the transition cost, while introduction of a larger ring increases it, largely confirming experimental activation shift results. There is a systematic correlation between the local aromatic ring rotation, the arginine barrier crossing, and the corresponding relative free energy. In particular, it appears to be more advantageous for the F233 side chain to rotate towards the extracellular side when arginines cross the hydrophobic region.


Introduction
Voltage-gated ion channels are membrane proteins that conduct ions, regulated by the electrostatic potential across the membrane. These channels play fundamental roles for instance in the generation and propagation of nerve impulses and in cell homeostasis. They are made up of four homologous domains, each of which contains six transmembrane helices. The first four make up a voltage-sensing domain (VSD), with the helices labeled S1 through S4. The final two helices (S5, S6) line the ion conducting pore together with the same helices from the other three subunits.
Within the VSD, the S4 helix is remarkably atypical for a transmembrane helix since it contains several charged side-chains, primarily arginine but also lysine and occasionally polar substituents. As the electrostatic potential across the membrane is changed, these charges will be subject to large forces that cause S4 to move in the voltage sensor. This in turn initiates a conformation change in the pore domain that opens or closes the channel to control the flow of K z ions from the intracellular to the extracellular side [1][2][3][4][5][6].
Significant progress has been made in understanding the structure and function of voltage-gated ion channels since the early days when it was not even clear that S4 was a transmembrane helix. In particular the X-ray structure of K v 1.2 [7,8] and later a higher-resolution structure of a K v 1.2/2.1 chimera (K v chim) [9] have been instrumental in this process. However, both these structures are for the open state of the channel, and it is only with recent models [10][11][12] as well as functional characterizations [13,14] and simulations [15] it has become possible to characterize the gating and resting state on atomic level. How are the S4 gating charges stabilized as they effectively move from one side of the bilayer to another, what is the character of the free energy barrier in this process, and are there stable intermediate conformations? What is known is that the gating charges bias the protein to flip between different functional states as a function of the applied voltage. Experimental studies have suggested at least three different functional and conformational states of the voltage sensor [16]; the resting (R) state corresponding to a closed channel when the membrane is polarized, the active (A) state where the channel is open due to depolarization, and the relaxed (AL) state which is an openinactivated channel reached after prolonged depolarization, which then can require higher and/or longer hyperpolarization to get back to the A and R states. These individual states are separated by free energy barriers [17], and by applying an external potential the relative free energy of the states change, which causes the voltage sensor to make a transition to a different state [18].
In practice, this motion occurs when the additional force on the arginines in S4 cause them to cross the hydrophobic core of the voltage sensor, breaking the salt bridges of the extracellular cluster (with E183/E226) [9]. During this process, each arginine must surmount the free energy barrier that comes from transiently losing the salt bridges, and potentially rotating or distorting the S4 structure [17]. The existence of these intermediate states is well supported experimentally [12].
As a consequence, the activation gating currents can have kinetics in the millisecond up to probably second range kinetics. Functionally, one typically separates between deactivation (A to R), activation (R to A) and inactivation (A to AL). Since the kinetics of these steps differ, there could be important differences either in the barrier, or at or at least differences related to the direction in which the free energy barrier has to be surmounted. The activation and deactivation barriers are likely to be caused by the arginines crossing a hydrophobic zone in the center of the VSD, while the inactivation barrier could also have components related to conformational change (e.g. ato 3 10 -helix transition [16,17,19,20]) in the S4 helix.
The number of charged residues that effectively move across the membrane has been the subject of some debate. A study by Khalili-Araghi et al. [20] found the highest-located arginine (R1) in S4 to interact with E183 and E226 in the resting state. This is consistent with the consensus model by Vargas et al. [11] and confirmed by a recent gating simulation [15] where three arginines, R2-R4, cross the barrier and R1 is located above the barrier. However, there are experimental indications that R1 too under some circumstances might be able to pass the barrier [12,14,21], although this might not be the physiologically most populated down state.
Since there are at least three arginines that cross, there should be intermediate states in the process, with S4 moving down in multiple steps towards the resting state as suggested already by Armstrong [22]. Tao et al. [21] suggested a model with five different voltage-sensor states. This would mean that each state indicating one charged residue passing the barrier, i.e. four separate energy barrier steps to overcome. These intermediate states appear to be well supported both by recent models based on metal bridges [12] and simulations of gating [15], although it is unclear whether the intermediate states are also metastable conformations.
The mixture of conserved positive charges in S4 and a hydrophobic region in the center of the VSD supports a switchlike transition between clearly defined arginine interactions above and below the hydrophobic core [9,23]. We have shown this quantitatively in recent simulations, and suggested it is easier for S4 to make this transition as a 3 10 -helix [17] based on a somewhat extreme model where the entire segment was forced into 3 10 helix.
While there are several residues in the core of the VSD that might affect the barrier (Fig. 1), Long et al. suggested the main candidate to be F233 (F290 in Shaker) in the S2 helix, which is almost universally conserved in voltage-gated potassium channels [9]. Tao et al. have conducted experiments on F233 mutations [21]. Only F233Y and F233W produced similar results to wildtype Shaker, but they were also able to get functional channels with a non-aromatic cyclohexylalanine side-chain, which suggests a rigid ring is important, but not aromaticity. This appears to be confirmed by our recent simulations that found F233 to act as a structural lock [17] as well as the clear separator for different states [12]. In contrast, Upadhyay et al. [24] found F233L to behave almost similar to the wild-type, while it produced a shift of 240-50 mV in the hands of Tao et al.. Campos and co-workers [13] found I230C to form a proton pore, which indicates this residue too contributes to the hydrophobic plug. Lacroix and Bezanilla [25] used gating current measurements of several F233 mutants to determine that a hydrophobic residue in position only controls the transfer of R4, the lowest gating charge, while the movement of R1-R3 is not affected. However the nature of the interaction between F233 and the gating charges remains unclear.
Notably, a substitution of the Phe side-chain to a Trp leads to a pronounced relative open state stabilization and consequently a massive left-shift in G/V (pdb: 3LMN) [21]. The replacement with Trp at this position has striking effects on voltage sensor function by a about 40 mV hyperpolarizing shift in activation [21], raising the possibility of an induced intramolecular cation-p interaction with S4. The Trp substitution was also qualitatively observed to result in significantly slower deactivation. An aromatic side-chain could catalyze the transmembrane passage of S4 charges through cation-p interaction and has been suggested to assist the S4 movement via this mechanism based on refined structural data of the K v 1.2 channel [8]. Pless et al. [26] have investigated the possibility that introducing tryptophan at position 233 could form an induced cation-p interaction in Shaker and found that the electronegative surface potential of the Trp contributes to a possible such interaction with K5 (Lys302).
To characterize the transition between states, we have created and equilibrated structural models of the first intermediate state during deactivation (one-down; referred to as closed-1 or C 1 ) of a voltage sensor based on experimental constraints [12], and used steered molecular dynamics simulations to slowly force it back across the barrier in the direction of the S4 helix axis. We believe this provides a much more natural transition pathway than our (and other) previous constrained studies, and in particular it includes a strict assessment since we can measure how closely the open end state reproduces the known X-ray structure. WIth an unconstrained high-quality reaction coordinate for the S4 translation we have been able to carry out potential of mean force free energy calculations at multiple equilibrium positions rather than nonequilibrium pulling (which leads to hysteresis effects). This has enabled us to perform systematic in-silico mutations in the hydrophobic core of the voltage sensor and directly estimate how much mutations affect the free energy barrier for arginine translation. Based on free energy landscapes for the F233 side-chain conformations we also quantify the molecular properties of the lock; the F233 side-chain dynamics appears to be highly important for the kinetics of the gating, and this residue might play the role to keep the the gating kinetics of potassium channels sufficiently slow compared to sodium channels.

A reaction coordinate for R4 crossing the hydrophobic core
The conformation of the VSD intermediate C 1 structure after free relaxation for 100 ns is shown in Fig. 2A (left). Helices S1-S3 are virtually identical to the X-ray VSD structure, while S4 has been translated one step down, i.e. the R4 arginine side-chain is located below F233 instead of above it. No restraints were applied during simulation, but the S4 helix spontaneously adopts 3 10 -helix structure between residues R3 and R4, exactly where F233 is located, further up in the sequence compared to the X-ray VSD. Limited-range umbrella sampling in both directions confirmed this initial model is very close to the local free energy minimum structure, which appears to confirm that these intermediate states are indeed metastable structures in the free energy sense. The steered molecular dynamics simulation forces S4 to move up, which in turn forces R4 to cross the hydrophobic core and instead adopt a position above F233. The advantage of this limited approach is that it specifically targets a single arginine translation and that the end state of the simulation can be critically compared to experiments. During the simulation, there is a spontaneous change of secondary structure in S4, with the 3 10 -helix region sliding to be most pronounced between R4 and K5. At the end of the upward steered simulation, both the side-chain positions and 3 10 -helix location correspond exceptionally well to the features observed in the X-ray structure of the VSD, as illustrated in Fig. 2A (right). Note that S4 is merely pulled along its own helix axis -no information about the X-ray structure is used -but the final simulation conformation reaches an RMSD of only 3.0 Å from the crystal structure (Fig. 2B). This strengthens our confidence in the quality of the obtained reaction coordinate compared e.g. to our own previous qualitative studies where we forced the entire S4 segment into a 3 10 -helix as a proof-of-concept in lack of detailed experimental structural information. In particular, these simulations strongly support the model with a conformational shifts between a-helix and 3 10 -helix, but that this shift occurs through gradual sliding of a constant-length 3 10 -region, rather than a net growth of either structure. In the simulation, the 3 10 -helix region is always located roughly opposite E226 and E236, at the same height as the F233 phenyl ring. Interestingly, this results in a stepwise local rotation of the arginine side-chains around S4 which is very similar to that reported in Fig. 1C of the recent gating simulation of Jensen et al. [15], although the secondary structure is not specifically discussed there.

F233 clearly contributes to the energy barrier
We have previously observed significant free energy barriers to S4 translation [17] at the point where charged side-chains enter the hydrophobic core, but this alone is not sufficient to uniquely identify the source of the free energy cost. To enable this, singlepoint mutations to alanine were introduced at all 21 residues within 1 nm ( Fig. 1) of the hydrophobic core (V170, M171, V172,  I173, L174, I175, I177, V178, F180, V225, C229, I230, I231,  W232, F233, I260, V261, I263, I264, Y267, V268), and PMF curves calculated as described in the methods. The relative difference in the peak free energy value along the PMF for each residue compared to that of the wild-type is listed in Table 1, and illustrated in Fig. 3 for a subset of positions. Most of the positions (15 out of 21) have little or no effect on the free energy barrier, and there are only two positions where there is a significant decrease when the side-chain is mutated to alanine: F180 and F233. Thus out of all tested residues F233 has the highest effect on the barrier.The shift for F233 is over 10 kJ/mol, more than twice as high as for F180. Both of these phenyl rings face the inside of the cavity, with F233 from helix S2 directly blocking the narrowest part of the VSD pore. In contrast, F180 is located in helix S1 above R4 close to R3 (Fig. 1). The barrier effect from this residue is more likely to come from changes in its interaction with R3 parallel to F233 is interfering with R4. In that case R3 is suspected to rotate slightly as well to overcome F180, but not as much compared to the rotational effect R4 has when crossing F233. In addition, F180 is not as conserved in voltage sensors as F233 is; this position can also contain other hydrophobic residues such as valine or leucine (Fig. 1A). There could be additional contributions from positions such as I230, but these effects are significantly smaller than for the two phenylalanines.

Effects of additional mutations on F233
Given the clear influence on the free energy barrier when the F233 side-chain is removed, it is highly interesting to consider the effect of introducing other amino acids in this position. Fig. 4 displays the PMF profiles for the other available aromatic sidechains, tyrosine and tryptophan, compared to the wild-type (WT). Neither of these mutations have any significant structural effect on the voltage sensor. Both F233Y and F233W exhibit higher free energy barriers, with in particular tryptophan making S4 motion significantly harder. Judging from the simulations, the main reason for this barrier is again simple steric hindrance, where all rigid rings provide obstacles, but the double ring of the tryptophan is particularly difficult to rotate away. This would in particular make it harder to deactivate the channel, which appears to agree with experimental results of Tao et al. [21] (in contrast to Upadhyay [24]) who found that F233W closing is significantly slowed down due to the mutation. The experiments also found a left voltage shift for F233W, but not F233Y (which instead behaved similar to WT). This is still compatible with the present study since the voltage shift describes the free energy difference between fully open O and fully closed C 3 states rather than the O-C 1 barrier kinetics. (Fig. 2C). Since there is no kinetics data available for F233Y it is not yet possible to directly compare the barrier to experiments for this mutation. The polar groups present both in tyrosine and tryptophane do not appear to play any significant role for their interactions with arginine. On the contrary, despite phenylalanine being more hydrophobic, the fundamental process for all three side-chains is still that the rings have to rotate away, which agrees very well with the observation that a rigid ring at the F233 position is important, but not its aromaticity.
In the steric obstacle model for the barrier, mutations of F233 to smaller residues should have the opposite effect of the large aromatic chains. This is confirmed not only by the F233A mutation, but also by F233L that replaces the phenyl with the most hydrophobic non-ring side-chain. As shown in Fig. 5, the leucine mutation reduces the free energy barrier by roughly 10 kJ/ mol, almost as much as the mutation to alanine. Thus, here too the main effect appears to be the bulkiness of the side-chain rather than its hydrophobicity. There is in fact one voltage sensor (K v AP) that has a leucine instead of a phenylalanine in the position corresponding to F233, but for virtually all others the phenylalanine appears to be conserved. This could have important implications for gating kinetics, but the relationship between K v AP other K v channels is not entirely certain [27].

F233 conformations during S4 translation
Based on the simulations, the phenyl ring appears to alternate between two preferred positions; it either has a largely horizontal orientation (with respect to the membrane frame-of-reference) where it effectively blocks the narrowest part of the voltage sensor cavity completely, or it rotates away to a largely vertical orientation when the R4 arginine side-chain needs to pass it. From our earlier simulations where the S4 helix was pulled down slowly from the X-ray structure on microsecond scale [17] it is possible to characterize the phenyl ring conformations in terms of the F233 side-chain torsions x 1 and x 2 as a function of time (Fig. 6). With the arginine side-chain entering from above (the arginine moves further down to the right in the plot), the ring is initially pushed even harder into a blocked position by R4 (the first 1-1.5 Å of translation), but after a while this lock opens by making a jump in the x 2 torsion. The x 1 torsion exhibits an initial fluctuation to enable the ring to rotate around x 2 and adopt a largely vertical orientation, but most of the change is observed in x 2 . Once R4 has translated down, both x 1 and x 2 go back to their initial values, and in particular for x 2 this is a very rapid, switch-like, process. This flipping behavior appears similar to that observed in other studies of aromatic ring motion [28][29][30][31][32] that supports the observation that In the C 1 confirmation the R4 arginine side-chain is located below F233, and there is a 3 10 -helix region between residues R3 and R4. After forcing S4 (see B) upwards along the helix axis the R4 arginine is located above F233, and the 3 10 -helix has to slide to be located between R4 and K5. This coincides remarkably well with the X-ray VSD structure (O state); the structures are within 3.0 Å RMSD of each other. The X-ray structure is superimposed in gray with arginines in transparent blue. Only the R0 side-chain has a slightly different conformation, which is due to interactions with a lipid. (B) C a RMSD of the VSD relative to the X-ray structure when starting steered molecular simulation from the C 1 model. Notably, the VSD was not specifically pulled towards any coordinates of the O state; the only external component was a force on the charged residues along the S4 helix axis to mimic a potential change. (C) Model for the first step of deactivation barrier (O to C 1 state) caused by R4 moving across the hydrophobic lock. While this study does not characterize the subsequent barriers, they are likely to be smaller, and other simulations too find the first step to be slowest [15]. Having one barrier of greater magnitude could help avoid trapping the S4 helix in intermediate states, while the population of open vs. closed channels will depend on the relative free energies of the O and C 3 states. doi:10.1371/journal.pone.0045880.g002 a phenylalanine primarily tends to rotate around the x 2 axis. For phenylalanine and tyrosine rings, these motions typically consist of restricted rotational diffusion in x 1 and x 2 degrees of freedom, combined with 180 0 ring flips around the symmetry axis in x 2 [33].
To further study the dynamics of the phenyl ring we performed 2D PMF calculations (see methods) to obtain the relative free energies of its conformations when the R4 side-chain is either above or directly opposite the F233 residue. As illustrated in Fig. 7A, the ring initially has relatively large free energy minima corresponding to the pore-blocking conformations where it has favorable hydrophobic interactions with the VSD core (x 1 * 270 0 , x 2 * 20 0 , conformation A 3 in Fig. 7A). The phenyl ring is pushed by R4 moving down into an even more blocked orientation. This effect has already been characterized in Fig. 6 (initial fluctuations). From these conformations the phenyl ring can rotate away, but this motion is not favorable. As the R4 arginine side-chain is forced down and relocated opposite F233 (Fig. 7B), steric clashes between R4 and the phenyl ring (Fig. 7B, configuration B 1 ) increase the free energy of the pore-blocking state, and eventually make this conformation disadvantageous, at which point the phenyl ring rotates by 100-120 0 to the vertical open state (x 1 * 270 0 , x 2 * 150-170 0 , conformation B 2 in Fig. 7B). This is followed by a flip-like motion of the arginine across the hydrophobic zone, after which the phenyl ring can rotate back into the lower free energy, pore-blocking, conformation.

Discussion
There is now increasing evidence -both from experiments and simulation -that the voltage sensor activation/deactivation process is closely coupled to secondary structure alterations in the S4 helix [16,17,19,20] or, equivalently, local rotation of the side-chain adjacent to F233 [15] rather than merely rotation, translation and tilting as a rigid element. The present study shows that our recent models of the intermediate C 1 state from experiments is a metastable conformation with a local free energy minimum, and it appears to directly confirm a 3 10 -helix region that slides along the S4 sequence as the helix moves from one intermediate conformation to the next. Rather than blindly moving between unknown states, we have shown the end conformation of the present steered molecular dynamics simulations to be remarkably close to the Xray structure of the chimera voltage sensor, both in terms of sidechain interactions, 3 10 -helix contents, and not least a 3 Å RMSD. In relation to our earlier studies qualitatively comparing barriers for constrained a-helix vs. 3 10 -helix translation, this confirms the advantage of 3 10 -helix structure close to the F233 lock [17] and the gradually changing secondary structure that largely maintains the total fraction of each secondary structure type [12] (thereby removing any net cost for 3 10 -helix formation). The standard error for the WT is the absolute standard error. A decrease or increase is defined by the cutoff + 5 kJ/mol. Compare to Fig. 3. In this context, it is worth mentioning the relative differences for Hanatoxin dissociation constants in the closed state of the drk1 VSD, as reported by Swartz [34] and Li-Smerin [35]. While the Hanatoxin data measures the relative stability of O vs C 3 states rather than the kinetics, their study found a 25-fold change upon mutating away the phenylalanine (F274 in drk1), which corresponds to a free energy stabilization of 8 Table 1). The difference is statistically significant both for F180 in the S1 helix and F233 in S2, with the latter showing the greatest shift. doi:10.1371/journal.pone.0045880.g003 Our potential of mean force calculations largely support the experimental results of Tao [21]. Based on the simulations, it is quite clear that F233 is the functionally single most important residue in the hydrophobic core and provides the main part of the R4 translation barrier. It also suggests that F233 is primarily a steric barrier for the arginine side-chains in the most narrow part of the VSD cavity rather than general hydrophobic effect. Larger aromatic rings increase the barrier despite containing some polar groups, while smaller but almost equally hydrophobic side-chains decrease it. The present PMF-based free energy barriers are also lower than the previous estimates, which makes sense since the S4 helix is now completely unconstrained and allowed to adopt whatever conformation is most advantageous rather than forced to assume either the X-ray or all-3 10 helix structure. However, the agreement is not perfect, in part because the studies are measuring slightly different properties: The simulations presently only cover the path from O through C 1 states, which means they can study the deactivation kinetics and also specific stabilization of the O state (relative to C 1 ). Tao et al. too observed a significant deactivation slowdown for F233W (no kinetics was reported for F233Y), but the left-shifted G/V curves are more related to relative stabilization of O vs. C 3 states. In particular, Tao found a left-shift for F233W but not F233Y, and since the present study does not observe any significant differences between O and C 1 for these mutants that difference would have to occur in the fully closed C 3 state. Alternatively, it could of course reflect shortcomings in the simulations, for instance different mutants having slight differences in the transition pathway. This is not realistic to sample with the present approach since it would require orders of magnitude more computer resources, but we believe this effect is limited since we do observe clearly lower barriers for some mutants (while they would be higher if we were only distorting the system). However, we believe this question can partly be resolved by considering the mutation data by Swartz [34] and Li-Smerin [35]. These authors measured the relative change in Hanatoxin (known to stabilize closed channels) dissociation constant from voltage sensors of the drk1 channel through alanine scanning, and found large relative changes both for the phenylalanine, the preceding residue, and the residue three positions later. This corresponds to the open state being stabilized (which we do not observe) -or the closed state destabilized -by 6-8 kJ/mol. However, there is also a smaller change in the same direction when mutating the arginines in S4. It is largest for R0 at the top of S4 (5.5-fold, 4.2 kJ/mol), and gradually drops for positions further down in the helix (R1 is glutamine in drk1; R2 5-fold, 4 kJ/mol; R3 3-fold, 2.7 kJ/mol; R4 2.3-fold, 2 kJ/mol). Since R0 & R1 do not interacting much with the VSD in the open state but they are centrally located in the closed one (and the opposite is true for R3 & R4), this pattern strongly suggests the mutations cause a destabilization of the closed (C 3 ) state of the channel, which agrees very well with the results of the present study.
In particular, it needs to be stressed that F233 is not merely a simple steric barrier since a substitutions by smaller hydrophilic residues do not cause omega currents, while it does occur for similar substitutions in R1 or E1. One could speculate that alterations in the salt bridge patterns are particularly important to allow omega currents, while F233's primary role is to regulate the gating kinetics, and the high conservation of this residue could be a key factor in keeping the gating kinetics of voltage-gated potassium channels sufficiently slow compared to sodium channels.
It is highly interesting that all the potential of mean force curves indicate a significantly lower free energy in the open (X-ray) structure to the left in the figures, while the intermediate C 1 state is clearly higher. Since the present results are based on equilibrium free energy calculations this is unlikely to be an hysteresis effect, and even less so considering we deliberately constructed the reaction coordinate starting from the intermediate state further down; if there were hysteresis effects, we would expect them to overestimate -not underestimate -the relative free energy in the final conformation. These results indicate a hypothetical model for the complete gating process could have a larger step for the first open to C 1 transition, followed by smaller barriers for the remaining steps. Not only would this be necessary to have a low enough free energy difference between the two end states that it is possible to alter it with the externally applied membrane potential, but it is likely also important to avoid having multiple deep minima where the voltage sensor could otherwise be trapped.  There are two limitations in the present model. First, we focus on isolated voltage sensors rather than complete voltage-gated ion channels. This is a slight limitation for more general studies of ion channels since pore conformations likely have some influence at least on the initial stage of deactivation, but it is a good model for the dynamics for isolated voltage-sensing-domains and the core concept of structural transitions in the domain. Second, it has not been possible to include the direct effect of an applied voltage in the present study since we want to study the free energy dependence on conformations during constant external conditions, which makes it difficult to model the O state as depolarized and C1 with some sort of (unknown) intermediate polarization.
Both these are important challenges for further studies, and combination with more experimental results for the kinetics of the voltage sensor activation should help address the questions whether the activation and deactivation processes are mirrors of each other, or have important differences. In particular the F233 properties found in the present study poses interesting questions about the connection between mutations, free energy barriers, and gating kinetics that might be possible to study with simulations. Similarly, a better knowledge of the energetics of the multiple steps should eventually help us understand the important differences between the open vs. open-inactivated conformations observed experimentally [16].

Experimental derived constraints
The strategy to achieve experimental constraints for the modelling of intermediate states, including the closed-1 (C 1 ) state, has recently been described in Henrion et al. [12]. Briefly, residue 325 in helix S3 of the Shaker H4 channel (Acc No NM_167595.3) [36] with the fast, N-type, inactivation removed [37] was probed for possible interactions with residues in a long stretch of S4 (355-369) by constructing double-cysteine mutants (one cysteine in S3 and one in S4). If the pair of cysteines are close to each other (v6.5 Å [12]) at some point during channel gating, those cysteines could coordinate a cadmium ion. If so, addition of CdCl 2 will induce the formation of a metal-ion bridge which alters the behavior of the channel. The effect of cadmium on each mutant was tested electrophysiologically in two-electrode voltage clamped Xenopus oocytes expressing the Shaker mutants. One double-cysteine mutant significantly affected by cadmium was 325C/361C showing altered voltage dependence and pronounced current reduction with slow recovery upon cadmium washout. The effect was attributed to the double-cysteine mutation, as the singlecysteine mutants responded differently to cadmium, and was interpreted as a strong metal-ion bridge.

Isolated voltage sensor simulation model
Helices S1 through S4 from the K v 1.2/2.1 chimera (pdb: 2R9R) [9] were used to model an isolated voltage sensor, since this is presently the highest-resolution voltage-gated ion channel available. Many experimental studies on voltage-gated channels [38][39][40][41][42] have instead been carried out on the Shaker K z channel, but the two have very high sequence similarity and their voltage sensors should be functionally equivalent. The simulation system was set up as previously described [17], with the modification of the Amber99SB-ILDN force field being used to model the protein [43] (Fig. 1B). All simulations were carried out with Gromacs 4.5.3 [44] using 5 fs time steps with virtual interaction sites for hydrogens.

Modeling the intermediate C 1 state & transitions
To overcome the lack of structural intermediates, a model of the state where S4 has translated one step down (C 1 state) was built using Rosetta using the experimental derived constraints in the same way as previously described [12]. The residues starting from S275 at the end of S3 through S307 at the end of S4 were rebuilt using the loop modeling protocol followed by all-atom refinement in Rosetta-membrane [45]. A cadmium-linker constraint between sulfurs of C268 in S3 and C289 in S4 (residues 325 and 361 in Shaker, Henrion et al. [12]) was applied in the form of a harmonic constraint centered at 6 Å . No particular secondary structure was enforced in the modelling, but the resulting C 1 model spontaneously adopted 3 10 -helix between residues 296 and 302, i.e. the region between R3 and K5 facing F233. Comparing with the other models, the 3 10 -part is conserved in this spatial region, and slides along the amino acid sequence when S4 moves between states. Steered molecular dynamics (SMD) was used to generate a trajectory between the two states by applying a harmonic potential to the center-of-mass of the R0-R6 charged side-chains and the center of mass of S1-S3 helices as described in [17], without constraining any parts of the S4 helix. The steered simulations were started from the C 1 state, pulling upwards in the direction along the helix, which makes it possible to assess to what extent the end state of the SMD simulation coincides with the the original VSD X-ray structure. Additionally, the X-ray structure state was also pulled downwards towards to obtain a transition coordinate for the deactivation process and get information whether C 1 is a metastable state.

In-silico mutagenesis & potential of mean force calculations
To investigate the barrier, frames spaced 0.02 nm of the S4 center-of-mass apart were selected from the wild-type SMD trajectories, and each hydrophobic residue in the core region of S1 through S3 (S4 is the moving part) mutated to alanine. F233 was additionally mutated to a number of other residues to assess the influence of the phenyl ring and shed light on differences between experiments. All single-point mutations were introduced with PyMol [46] using the mutagenesis plugin. Each mutated system was then energy minimized for 1,000 steps with steepest descent to resolve steric clashes.
Potential of mean force (PMF) curves were calculated by using umbrella sampling on the relative distance between the charged side-chains in S4 (R0-R6) and the reference group formed by helices S1-S3. The umbrella harmonic force constant was set to 10,000 kJ/mol/nm 2 , and the reference distance in each window was kept constant during the run. Each umbrella point was sampled for 50 ns, with the first 20 ns (40%) being reserved for relaxation and the remaining 30 ns used for production. Standard errors were estimated by splitting the data into three parts and using jack-knife statistics. PMF curves were calculated from the ensemble of simulations using the Gromacs g_wham program [47] Figure 7. Relative free energy of F233 phenyl ring orientations. 2D PMFplots of relative free energy as a function of side-chain torsions x 1 and x 2 , with the R4 arginine (blue) above (A) and at the same level (z-wise) as (B) F233 (green). E226 is shown in red. There appears to be two possible pathways for F233 to rotate away and open the lock, the first involving a switch-like transition in x 1 , and the second only rotating around x 2 . The lowest-barrier path is shown in a bold dashed line. As R4 is located opposite of F233 (B1), the ring's conformational space is restricted by steric clashes (red circle), and it can only rotate around x 2 until the arginine has moved further down. doi:10.1371/journal.pone.0045880.g007 based on the Kumar et al. Weighted Histogram Analysis Method [48].

2D free energy profile of the phenyl ring rotation
The conformations of the F233 phenyl ring and its interactions with arginine side-chains during S4 translation was investigated through 2D umbrella sampling on the side-chain torsions in F233. This was performed for two different S4 positions, first with the R4 arginine located above F233, and second when directly opposing F233, which corresponds to the peak in the PMF curves. The umbrella term was introduced by modifying parameters of sidechain torsions x 1 (N-Ca-Cb-Cc), x 2 (Ca-Cb-Cc-CD 1 ) to force the phenyl ring into different conformations. The biasing potential was set to 100 kJ/mol/rad. For x 1 , sampling was performed in the range [2180 0 , 180 0 ] with a spacing of 10 0 , while x 2 was limited to [0, 180 0 ] because of the ring symmetry. Each umbrella simulation was relaxed for 100 ps followed by 1 ns of production. 2D PMF profiles were computed from the simulations using Alan Grossfield's 2D-WHAM program [49].

Author Contributions
Conceived and designed the experiments: CS BH FE EL. Performed the experiments: CS BW. Analyzed the data: CS SB BH BW FE EL. Contributed reagents/materials/analysis tools: CS SB BW. Wrote the paper: CS SB FE EL.