Letter pubs.acs.org/JPCL
Simulating the Activation of Voltage Sensing Domain for a VoltageGated Sodium Channel Using Polarizable Force Field Rui-Ning Sun and Haipeng Gong* MOE Key Laboratory of Bioinformatics, School of Life Sciences, Tsinghua University, Beijing 100084, China S Supporting Information *
ABSTRACT: Voltage-gated sodium (NaV) channels play vital roles in the signal transduction of excitable cells. Upon activation of a NaV channel, the change of transmembrane voltage triggers conformational change of the voltage sensing domain, which then elicits opening of the pore domain and thus allows an influx of Na+ ions. Description of this process with atomistic details is in urgent demand. In this work, we simulated the partial activation process of the voltage sensing domain of a prokaryotic NaV channel using a polarizable force field. We not only observed the conformational change of the voltage sensing domain from resting to preactive state, but also rigorously estimated the free energy profile along the identified reaction pathway. Comparison with the control simulation using an additive force field indicates that voltage-gating thermodynamics of NaV channels may be inaccurately described without considering the electrostatic polarization effect.
V
partial remedies, due to their considerable sequence similarities to NaV channels.35,36 Molecular dynamics (MD) simulations have been employed to study the large-scale conformational change of VSD during voltage gating for NaV channels19,31 as well as other VGICs.37−39 However, all previous simulations adopted additive force fields and thus ignored charge redistribution,40−42 a phenomenon that is non-negligible, especially for the gating process where multiple positive charges uniformly move in one specific direction. In this work, we investigated the VSD activation of a prokaryotic NaVAb channel (from Arcobacter butzleri),7 using a polarizable force field for the first time to our knowledge. Different from previous studies, we first generated a molecular model of the resting state using Ci-VSP35 as template, and then simulated its spontaneous conformational change at the depolarization TM voltage of 0 mV. The simulation finally reached a conformation that was strikingly close to the crystal structure of NaVAb, which allowed us to identify the reaction pathway between the resting and preactive states and to rigorously evaluate the potential of mean force (PMF) along the pathway. Comparison with the PMF calculation using an additive force field suggested that electrostatic polarization should be considered for the quantitative description of VSD gating. Compared with AtTPC1,36 Ci-VSP35 shows a higher sequence similarity to the VSD of NaVAb, especially at the highly conserved residues (Figure S1A). Therefore, we chose Ci-VSP as the template to model the resting state of NaVAb. In contrast to the crystal structure, the resting-state model has the
oltage-gated sodium (NaV) channels control the transmembrane (TM) influx of Na+ ions and are thus responsible for initiation of action potentials in excitable cells including neurons, muscles, and neuroendocrine cells.1,2 NaV channels are therefore essential for proper functioning of organisms, and their malfunctions are linked to many diseases.3 Unlike the complicated architecture of eukaryotic NaV channels,4−6 the homotetrameric organization of prokaryotic NaV channels simplifies their structure determination. As suggested by the reported prokaryotic NaV structures,7−12 each subunit consists of six TM helices (denoted as S1−S6), where S5−S6 jointly constitute the pore domain (PD) while S1−S4 compose four individual voltage sensing domains (VSDs). Upon membrane depolarization, the S4 helix that carries multiple positive charges is proposed to move in the outward direction, eliciting gating currents2 and subsequent channel opening. Despite extensive studies on the gating mechanism,4,12−32 general consensus has not been reached for NaV channels, nor is the extent of conservation of mechanistic details among the voltage-gated ion channel (VGIC) family known.33 However, a recent study implies that electrostatic properties governing the electric component of free-energy landscapes are likely conserved in distinct VGICs,34 which lays the foundation of mechanistic investigation on voltage gating using individual channels. Moreover, since available NaV structures are determined at either preactive or inactive states,7−9,12 structural properties of the resting state is unknown, which further hinders atomistic illustration of the gating process. Fortunately, recently determined VSD structures from the Ciona intestinalis voltage-sensing phosphatases (Ci-VSP) and the two-pore channel from Arabidopsis thaliana (AtTPC1) may provide © XXXX American Chemical Society
Received: January 4, 2017 Accepted: February 7, 2017 Published: February 7, 2017 901
DOI: 10.1021/acs.jpclett.7b00023 J. Phys. Chem. Lett. 2017, 8, 901−908
Letter
The Journal of Physical Chemistry Letters
Figure 1. (A) Sequence alignment for VSDs of NaVAb and Ci-VSP. Only S3 and S4 helices are compared, with highly conserved residues highlighted. (B) Time-dependent RMSD profiles for all Cα atoms in S1−S4 helices and Cβ atoms of gating charge residues (R1−R4) during the 100 ns cMD and 190 ns aMD simulations. Comparisons with the NaVAb crystal structure and resting-state model are colored in blue and green, respectively. (C,D) The resting-state model after cMD equilibration (C) and the final snapshot from aMD simulation (D) aligned against the NaVAb crystal structure. Cβ atoms of gating charge residues (R1−R4) are shown in spheres, with the crystal structure and simulated structures colored in blue and green, respectively.
Figure 2. (A) The PMF profile along the reaction coordinate, with five metastable states labeled from A to E. (B) The free energy profile plotted as a function of gating charge. (C−F) Relative vertical positions of geometrical centers of different groups with respect to the center of membrane bilayer (red line) in the transition from state A to E: backbone atoms of R1−R4 (C), side-chain guanidyl groups of R1−R4 (D), backbone atoms of S4 helix (E), and side chains of conserved residues (Asn25, Glu32, Asn49, Glu59, and Asp80) in S1−S3 helices (F). R1, R2, R3, and R4 are colored in orange, blue, magenta, and green, respectively.
force field that shows significant improvement in accuracy40−43 and successful application for biological macromolecules,44−47 to simulate the resting-state structure model of NaVAb VSD alone in a membrane (Figure S2A). The conventional MD (cMD) simulation (Sim #1 of Table S1) supports the structural stability of this model, with root-mean-square-deviation (RMSD) and root-mean-square-fluctuation (RMSF) of the
four conserved Arg residues on the S4 helix (denoted as R1− R4 and also named as gating charge residues; see Figure 1A) remarkably relocated toward the cytoplasmic side, with R4 positioned at a similar level to AtTPC1 (Figure S1B). The widely used additive force fields have limitations in simulation accuracy, because of the neglect of charge redistribution.40 In this work, we adopted the DRUDE force field, a polarizable 902
DOI: 10.1021/acs.jpclett.7b00023 J. Phys. Chem. Lett. 2017, 8, 901−908
Letter
The Journal of Physical Chemistry Letters
Figure 3. (A) The hydrogen or ionic bonds formed between gating charge residues (R1−R4) and conserved residues on S1−S3 helices as well as the lipid headgroup of lower and upper leaflets. Partners of gating charge residues are arranged according to their relative position along z-axis (identical coloring style to Figure 2C). (B) Gating charge residues and their conserved interacting partners are labeled in the VSD structure. (C) Rotation of gating charge residues (R1−R4) around the axis of S4 helix in the transition from state A to state E (identical coloring style to Figure 2C). (D) Distribution of relative angles of conserved residues in S1−S3 helices that interact with gating charge residues. The inset panel shows a structural top view of these residues, where the S4 helix is represented as a blue cylinder.
the tendency of spontaneous activation observed in the aMD simulation (Figure 1B). The presence of a free energy barrier of 3.358 kcal/mol between states B and C agrees with the observation that only states A and B could be sampled in the cMD simulation (Figure S7). Similarly, the free energy barrier of 5.473 kcal/mol in the reverse direction explains the rare observation of VSD deactivation in previous simulations on the NaVAb crystal structure. Notably, due to the uncertainties in free energy calculation, the PMF profile could also be explained by two stable states (B vs D+E separated by an energy barrier). Following the protocol of Treptow et al.52 and Delemotte et al.,39 we estimated the gating charge for each individual BEUS window (Figure S8). The resting and preactive states show a difference of ∼2 e in gating charge, which is less than previous measurements of 3−3.75 e for the full activation process (12− 15 e for 4 VSDs).26,53 The underestimation in gating charge suggests that we only observed a partial activation process. On one hand, our calculations stop at the preactive state and therefore neglect the ongoing transition from preactive to fully active state that accounts for ∼20% of the overall change of gating charge.19 On the other hand, our resting-state structure model might be an early intermediate formed upon activation (see Discussion for details). Nevertheless, the free energy could be plotted as a function of gating charge, and the profile thus generated exhibits a similar shape to the PMF profile (Figure 2B). By neglecting the coupling between external TM voltage and charge polarization, we could roughly predict the free energy profiles at TM voltages of −100 mV and 100 mV respectively, following ΔG(Q,V) = ΔG(Q,0) − QV,39 where Q and V are the gating charge and applied TM voltage, respectively. As shown in Figure S9, at the TM voltage of −100 mV that is close to resting potential, state B is greatly
TM domain steadily 10 ns (Figure S3). Furthermore, polarizable π electrons of Phe could form favorable electrostatic interactions with the positively charged Arg and Lys,64,65 which may stabilize the protein structure when gating charge residues pass the HCS during activation. Notably, such favorable interactions cannot be properly handled without considering the polarization effect. Consistently, representative structures extracted from the transition state of the PMF profile derived using the polarizable force field suggest nearly parallel stacking between R3 and Phe56 at short distance, when compared with similar structures obtained using the additive force field (Figure 4CD). We compared the electrostatic interaction energy between R3 and Phe56 in simulations using the two different force fields (Figure 4B). Clearly, the interaction is significantly more favorable when simulated using the polarizable force field, especially at the reaction coordinates where the PMFs show great difference. Therefore, the inaccuracy of additive force fields at least partially arises from their failure in handling electrostatics involving electrically neutral residues that are highly polarizable (e.g., aromatic residues). Although our calculations support the importance of favorable interaction between R3 and Phe56, this interaction may contribute at most ∼2 kcal/mol in respect of free energy, considering the contribution of the other polarizable residues (Figure 4A). Notably, an increase of energy barrier by
