Introduction

Voltage-dependent gating of ion channels mediates numerous physiological and pathophysiological processes (1–5), by regulating ion flows in response to membrane potential changes in excitable cells such as cardiac myocytes and smooth muscle cells in the lung and blood vessels (6–11). Disruption of normal voltage gating is involved in many cardiovascular and neurological diseases, including epilepsy, mental retardation, chronical pain, hypertension, arrhythmias, stroke, and ischemia (12–20). Among the family of potassium channels, the large-conductance potassium (BK) channels stand out in several ways (21–27). It has the largest single-channel conductance (up to ∼300 ps) and is the only known K⁺ channel to be activated by both intracellular Ca²⁺ and membrane potential. Functional BK channels are homo-tetramers, with each subunit consisting of a seven-helix transmembrane domain (TMD) and a C-terminal Ca²⁺- sensing cytosolic domain (CTD) (Figure 1a). The deactivated state of BK channels contain a physically open central pore (28–30), lacking the classical bundle crossing constriction at the intracellular entrance observed in the voltage-dependent potassium (Kv) channels (31). Instead, BK channels likely follow the hydrophobic gating mechanism, where the pore undergoes hydrophobic dewetting transition to create a vapor barrier to block ion permeation in the deactivated state (32–36). Furthermore, the voltage sensor domains (VSDs) of BK channels are not domain-swapped with respect to the pore-gate domain (PGD), in contrast to the domain-swapped TMD organization of Kv channels (Figure 1b-c). The difference in TMD organization is likely a key factor underlying important but still poorly understood differences in the voltage sensing and gating mechanisms of BK channels (37–41) in comparison to the “canonical” framework established mainly through the study of Kv channels (42–50).

Figure 1.

In Kv channels, the linker connecting transmembrane (TM) helices S4 and S5 forms a α-helix (∼15 aa) that wraps around the S6, with the S4 of one subunit interacting with the S5 of a neighboring subunit and simultaneously connecting to the S5/S6 within the same subunit through the S4-5 linker. Conversely, in BK channels, the S4-S5 linker is a short loop (∼5 aa), and the S4 solely interacts with S5 within the same subunit (Figure 1). Within the canonical voltage activation of Kv channels, each of the first four Arginine residues in the VSD S4 accounts for about 1e gating charge (51–54), giving rise to a total gating charge of ∼ 3 - 3.5 e per monomer. Membrane depolarization drives S4 to move upward by ∼8 Å, which causes the S4-S5 linker helix to pivot upwards and change the interactions with S6 and releases the S6 helix bundle crossing to physically open the gate (55–57). For BK channels, despite sharing many conserved charged residues, the total number of gating charges is only 0.6 e per VSD (22, 58–62). In addition to Arginine residues on BK S4 (R207, R210 and R213), other residues also appear to contribute to voltage sensing, such as D153 and R167 on S2, D186 on S3 and E219 on S4 (39, 63, 64). Furthermore, fluorophore quenching suggested that the movement of BK VSDs during voltage gating was smaller and involved both vertical and lateral motions (65–67). At present, there remains significant ambiguity in the identity of gating charges, details of VSD motions, and how they drive the pore opening in BK channels.

In this study, extensive atomistic molecular dynamic (MD) simulations up to 10 μs in length were performed in explicit solvent and membrane to directly probe the voltage- driven activation of BK channels in the Ca²⁺-free state, using the Core-MT construct that does not include the CTD but retains voltage gating (68, 69). We were able to directly observe spontaneous opening transitions and ion conductance of the channel under 750 mV membrane voltage within the 10 μs simulation timescale. The observed VSD movements appear to be highly consistent with two recent Cryo-EM structural studies of mutant BK channels with constitutively activated VSDs (38, 41). Free energy analysis was preformed to further quantify the total gating charge as well as the contributions of key charged residues. The results are in quantitative agreement with available experimental data (39, 63, 64), providing further support for the predicted voltage sensing and gating mechanism. The analysis reveals central roles of voltage- induced displacement of R210 and R213 on S4 in voltage sensing and the strong interactions at the S4-S5-S6 interface in VSD-pore coupling. This mechanism was further validated using steered MD simulations showing that pulling on the charged tip of R210 and R213 alone could ready drive pore opening in Core-MT BK channels. Taken together, the current work provides for the first time a reliable detailed molecular mechanism of BK voltage activation, which is distinct from the canonical voltage gating mechanism of Kv channels. These novel voltage gating principles are likely relevant in understanding the gating and regulation of other non-domain swapping ion channels.

Results and Discussion

Direct atomistic simulations of voltage activation of BK channels

Starting from a fully equilibrated closed state, multiple atomistic simulations were performed at 0 and 750 mV membrane voltages for up to 10 μs to directly probe voltage-driven activation of Core-MT BK channels using the special purposed supercomputer Anton 2 (70, 71) (see Methods; Figure S1). The simulated voltage of 750 mV is higher than the experimental conditions, where V1/2 of Core-MT BK channels is ∼ 244 mV compared to ∼180 mV for the full-length BK channel (69), to accelerate protein conformational transitions. Similar voltages have been used in atomistic simulations and found to generate realistic transitions (56, 72). We compare the membrane thickness at 300 and 750 mV and the results reveal no significant difference in the membrane thickness (Figure S2). At 0 mV, the Core-MT BK channel remained highly stable in an apparently closed state, with the pore fully dehydrated and not permeable to ions (Figure S1, top row). Similar to what has been observed previously in simulations and cryo-EM maps (30, 32–35, 38), lipid tails can enter the dewetted pore through the fenestration gap between the pore lining S6 helices. They are highly dynamic and likely contribute to the stability of the dewetted state of the pore. Note that the absence of CTDs increases the flexibility of pore-lining S6 helices, and the fully relaxed pore profile (red trace in Figure S1d, top row) shows substantial differences compared to that of the Ca²⁺-free Cryo-EM structure of the full-length channel. Importantly, the VSDs were highly stable and showed minimal movements (∼1 Å or less) for all charged groups as well as the TM helices themselves.

At 750 mV, significant movements were observed with multiple charged groups on S4 (Figure 2a and Figure S1). In particular, the guanidium groups of R210 and R213 moved upward along the membrane normal (z-axis) by up to ∼8 and ∼10 Å, respectively, plateauing after ∼5 μs (e.g., see Figure 2a). Other VSD charges showed much smaller z-displacements. The apparent modest movements of ∼3-4 Å of E219 and R207 seem to be mainly a result of the upshift of S4 itself, which reached a maximum of ∼3 Å during the second half of the simulation (Figure 2b). In contrast, other VSD helices S1-S3 exhibited much smaller z-axis movements (∼1 Å). As illustrated in supplementary Movie S1, the large charge displacements of R210 and R213 charges with modest S4 movement were enabled by side chain snorkeling.

Figure 2.

Voltage-driven VSD movements were apparently coupled with an opening transition of the pore, allowing one to directly observe Core-MT activation and K⁺ conductance during both 10-μs atomistic simulations on Anton 2 (Figure S1). As illustrated in Figure 2 and Movie S2, the backbone root-mean-squared distance (RMSD) of the pore from the open state decreased sharply from ∼4.5 Å around the 4 μs mark to below ∼2.7 Å around the 5 μs mark, as the S4 helix shifted up along z-axis by ∼ 2 Å during the same time span. The pore opening transition mainly involved the increase of the tilt of S6 segments below the glycine hinge (residues 313-324) (Figure 2d). The pore continued to move closer to the open state after the initial rapid response to VSD activation, eventually reaching a state that has similar S6 helix tilt and is only ∼2.2 Å from the Ca²⁺-bound Cryo-EM structure at the end of the 10 μs run (Figure 2e, red vs. cyan cartoons). These pore structural changes roughly doubled the pore diameter, from ∼10 Å to ∼20 Å at the intracellular entrance (z ∼ − 20 Å; Figure 2f). It should be noted that the full-length BK channel in the Ca²⁺ bound state has an even larger intracellular opening (Figure 2f, green trace), suggesting that additional dilation of the pore may occur at longer timescales. On the other hand, the simulation construct does not include the 11-residue Kv mini-tails required for assembly and membrane insertion of Core-MT (68), which could impact the pore conformation. Furthermore, the single- channel conductance of Core-MT is ∼30% lower (68), suggesting that its open pore is more constrictive than that of the full-length BK channels.

The opening transitions were accompanied with rehydration of the pore and the channel became conductive to K⁺ (Figure 3 and Movie S3). Consistent with the smaller pore diameter, the final states of the pore in the simulations only accommodate ∼30 waters, compared to ∼40 waters for the Ca²⁺-bound state of full-length BK (32).

Figure 3.

Furthermore, the single-channel conductance estimated from the last 2 μs of sim 2b (Figure 3a) is only ∼1.5 pS, much lower than the experimental value of ∼220 pS for Core-MT (69). We note that classical force fields, such as CHARMM36m used in the current simulations, are known to underestimate the single channel conductance by about one order of magnitude (73). Indeed, the fully opened state of Core-MT, constructed from the Ca²⁺-bound Cryo-EM full-length BK structure, is predicted to have a conductance of ∼6 pS using the same simulation setup (sim 7; Figure S3), which is only ∼4 fold of that of the final state from the voltage activation simulations. We further note that subconductance open states have been observed in single channel recordings of BK channels (74–77). Nonetheless, the ability of the dilated and hydrated pore to sustain K⁺ permeation throughout the second half of the 10 μs simulation demonstrates that the hydrophobic gate in the closed channel has been broken and the channel enters a conductive state. Consistent with previous simulations (78–80), the conductance follows a multi-ion mechanism, where there are at least two K⁺ ions inside the filter, preferentially occupying positions S1/S3 or S2/S3 (Figure 3b). The two bound ions subsequently move to S0/S1 when an incoming ion takes the S3 position. Notably, between the ions inside the filter, there can be either one or no water molecule, indicating the coexistence of both soft and hard knock-on mechanisms.

Gating charges: how small VSD movements support voltage sensing

The ability of Anton 2 simulations to directly observe voltage-driven activation of Core- MT BK channels allows one to examine the molecular details of VSD activation and VSD-pore coupling. An important observation is that, even though the charged groups of R210 and R213 show z-displacement of 8-10 Å (Figure 4a), the overall movement of S4 along the membrane normal is only ∼3 Å, less than a single α-helical turn (Figure 2b). This is considerably smaller than that observed in canonical Kv channels, which has been estimated to be ∼ 8 Å (55) and up to 15 Å (81). Smaller S4 movements have been suggested in fluorophore quenching experiments (65–67) and molecular simulations (39). Recent Cryo-EM structures of mutant BK channels with constitutively activated VSDs actually reveal minimal S4 movements along the membrane normal (38, 41). Upon voltage activation, R210 switches the salt-bridge partner from D186 on S3 to D153 on S2 and D133 on S1, R213 switches to engage with D153, and R207 switches from engaging with D153 and D133 to becoming exposed to the water/membrane interface (Figure 4b). Similar movements were also observed in recent high-resolution Cryo-EM structures of R207A mutant BK channels with constitutively activated VSDs at 0 mV (41), even though the net z-displacements of R210 and R213 charges in Cryo-EM structures are about ∼2-4 Å smaller. Both the smaller side chain charge movements and a lack of overall S4 z-displacement in the mutant Cryo-EM structures may be due to the absence of membrane voltage. To further evaluate if z-placements observed at 750 mV are an artifact of unphysical voltage, three independent simulations were initiated from the final state of sim 2b at 300 mV for 1.0 μs (sim 9 in Table S1). The results, summarized in Figure S1, show that the activated state of VSD remained stable and the pore stably hydrated in all three simulations.

Figure 4.

Another important feature is that multiple helices in VSD underwent substantial outward movements in the membrane lateral direction in addition to the modest z- displacement of S4. As shown in Figure 4c, both S1 and S4 displayed ∼3 Å outward movements within the x-y plane. These lateral movements are apparently important for driving the opening motion of the pore-lining S6 helices (Figure 4d), likely mediated by direct S4-S5-S6 interactions (see the next section). Importantly, the predicted VSD movements appear highly consistent with observations from fluorophore quenching experiments (65–67), which suggest that S4 likely moves towards S3 and closer to S1/S2 upon activation. Interestingly, the simulations predict that S1 moves during voltage-activation instead of S2, in contrast to the previously proposed model (66). Lateral movements of S4 at the cytosolic end have also been observed in recent Cryo- EM structures of R202Q mutant Aplysia BK channel (equivalent to R213Q in human BK), which presumably locks the VSDs in the activated state (38), even though the overall movements of VSD are much more subtle in the Cryo-EM structures. The later may again be a direct consequence of the absence of membrane voltage.

To understand how the modest movements of VSDs can support effective voltage sensing in BK channels, we calculated the average electrostatic potential maps for the resting and activated states (Figure 5). The results reveal how the protein significantly remodels the local electric field, creating large gradients along both the membrane normal and lateral directions. This is similar the electric field focusing effects proposed for VSDs of both Kv and BK channels (39, 82). In particular, residues R210 and R213 reside in a region of high electrostatic potential (> 1100 mV) in the resting state, but move to a region of much lower electrostatic potential (∼ 600 mV) in the activated state. Therefore, despite the modest S4 movements, charges on R210 and R213 can effectively sense a voltage change of ∼500 mV, compared to the total imposed membrane voltage of 750 mV.

Figure 5.

We performed free energy calculations (83, 84) to further quantify the contributions of key VSD residues to the total gating charge of BK channels (see Methods; Table S1 sim 3-6). As summarized in Table 1, the calculations estimated the total gating charge per VSD to be ∼ 0.45 e, in strong agreement with the range of 0.48 - 0.65 e measured experimentally for the full-length BK channel (39, 63, 64). The results further identify R210 and R213 as the primary contributors to voltage sensing, accounting for approximately 97% of the total gating charge per VSD (Table 1). The calculated residue contributions to gating charge are highly consistent with the latest experimental measurements by Carrasquel-Ursulaez et al (2022) (39), even though earlier experimental studies disagreed on the contributions of several residues, including R207, R167, D153 and D186. It has been recognized that interpretation of gating charge measurements on mutant BK channels is nontrivial (39). For example, neutralizing R213 reduces the total gating charge from 2.62 e for WT to 1.31 e for R213C; yet replacing it with a negative charge (R213E) does not further reduce the net gating charge (63). The simulations reveal that the complication in experimental analysis may be attributed to modest movements of VSD in both membrane normal and lateral directions and a clear interplay of VSD conformation and local electric field (Figure 5). Nonetheless, atomistic simulations together with electrostatic potential and free analyses strongly support a central role of R210 and R213 in BK voltage sensing.

Table 1

The dominant role of R210 and R213 in voltage sensing of BK channels was further tested by steered MD simulations (sim 8), where biasing potentials were applied to steer the movement of the guanidium CZ atoms of R210 and R213 from their positions in the initial resting state to those in the activated state in absence of membrane voltage (see Methods and Figure S4). The simulations show that upward movements of R210 and R213 guanidium tips alone can readily drive pore opening, which were observed in two out of four independent steer MD simulations (Figure S4). In the other two simulations, the S5-S6 packing was ruptured during the first 100 ns, likely due to the rapid rate of pulling. Similar to simulations under 750 mV membrane potential, pore dilation and increase in S6 tilting were observed as the VSD responded to the movement of R210 and R213, along with the breakdown of the hydrophobic barrier, as indicated by an increase in the number of water molecules within the pore. The ability of R210 and R213 charge upward movement alone to drive pore opening in absence of membrane voltage further supports the dominant roles of these two charges in voltage sensing and gating of BK channels.

Central role of the S4-S5-S6 interface in VSD-pore coupling

Examination of the covariance matrices (Figure S5) reveal that all TM helices, particularly, S1-S5, are tightly coupled in both resting and activated states. The lateral movements of S1 and S4 during activation mainly involve the cytosolic ends, concerted with those of S5 and S6 from the pore (Figure 4b). In the Cryo-EM structures, a kink is observed in S4 below R213, resulting in approximately a ∼30° bend in the C-terminal half away from S5 and the pore. During voltage-induced activation, the N-terminal region of S5 maintained strong contacts with the S4 C-terminal segment (Figure S6), such that the upward and outward movement of S4 induces an increase in the tilt and bend of S5 (Fig. 4d). This S4-S5 movement in turn increases the tilt angle of S6 C- terminal half (F315 to E324) by ∼15° to dilate the pore. In particular, the hydrophobic helix-helix contacts between S5 (L235, L239 and F242) and S6 (F315, V319 and I322) are strengthened, while those between N231 and K234 on S5 and E321 and E324 on S6 are weakened during this process (Figure S6). Besides pore dilation, the conformational response of S6 alters the orientation of the conserved E321 and E324 at the cytosolic entrance from membrane interface-facing in the closed state to pore- facing in the open state, as observed in Cryo-EM structures (28–30). The resulting increase in both pore size and surface hydrophilicity promotes pore hydration and disrupts vapor barrier to render the channel conductive.

To further identify structural elements central to BK VSD-pore coupling, we performed unbiased dynamic community and coupling pathway analysis based on the covariance matrices derived from atomistic simulations (see Methods). The results, summarized in Figure 6a, reveal that S4-6 forms a single community distinct from the rest of the VSD (and S0 helix), where the motions of residues are more strongly coupled to each other than to the rest of the channel. Indeed, the optimal and suboptimal pathways of allosteric coupling between S4 R213 (a major voltage sensing charge) and S6 E321 (a key position of pore dilation and hydrophilicity increase) pass exclusively through S5 residues (Figure 6b). The central role of S5 in mediating VSD-pore coupling is further supported by analyzing the information flow betweenness (85), which provides a global measures how conformational perturbation flows through the network from the “source” (e.g., R213) to the “sink” (e.g., E321). The analysis reveals that S5 residues have the largest contributions to the informational flow besides neighboring residues on S4 or S6 (Figure 6c). Other residues on S1 or S3 also contribute significantly to the flow, albeit at much lower levels compared to S5. Taken together, the simulation and dynamic network analysis consistently point to a central role of the S4-S5-S6 interface in VSD-pore coupling of BK channels.

Figure 6:

Concluding Discussion

Long atomistic simulations in explicit water and membrane performed using special purposed supercomputer Anton 2 have allowed direct observation of voltage-induced activation of the Core-MT BK channels for the first time. These simulations provide crucial new insights into how BK VSDs may sense membrane voltage and how the VSD movements may drive the pore opening and activate channel conductance. In particular, the simulations support that the S4 helix remains the major voltage sensing element and that it undergoes modest movements of ∼ 3 Å along both the membrane normal and lateral directions. The ability of BK channels to utilize relatively small VSD movements to sense membrane voltage and drive channel activation is likely attributed to the protein’s ability to modify the local electric field, which leads to large electrostatic potential gradient in both vertical and lateral directions. Further analyses reveal that R210 and R213 are the major voltage sensing residues in the wild-type Core-MT BK channel, contributing ∼97% to the total gating charge of ∼ 0.45 e per VSD. These core features of BK voltage sensing from simulation are largely consistent with a range of existing functional, biophysical and structural studies (38, 39, 41, 66), even though structural studies mutant BK channels with constitutively activated VSDs reveal smaller movements of R210 and R213 guanidinium moieties and minimal vertical movement of S4. The later observations may be a consequence of the complete absence of membrane voltage during structural determination. Importantly, control simulations show that the activated VSD state remains stable at 300 mV, suggesting that the observed vertical S4 movement is not an artifact of high voltage of 750 mV.

Another distinct feature derived from atomistic simulation is that the tightly packed S4- S5-S6 interface in BK channels plays a central role in the VSD-pore coupling, which is further supported by dynamic community, coupling pathway, and informational flow analyses. The modest movements in S4 drive a series of concerted conformational changes in S5 and eventually in S6, with the S4-S5 linker not playing a major role as observed in the canonical Kv channels (55–57). This is supported by the findings that mutations in the S4-S5 linker (N225-K228) do not significantly affect the coupling between VSD and the pore (86). Arguably, the distinct features of voltage-sensing and pore-sensor coupling observed in BK channels are intimately related to the novel non- domain-swapped TM topology. In particular, the VSD and pore domains are more tightly packed compared to the domain swapped configuration (Figure 1), which restricts the VSD motions and provides a direct pathway for their inter-talk. Interestingly, a similar non-canonical pathway for VSD-pore coupling, involving the S4/S1 and S1/S5 interfaces, has recently been proposed for the non-domain-swapped cardiac hERG potassium channel (87). We also note that a similar non-canonical VSD-pore coupling involving S4-S5 interactions between neighboring subunits has been suggested to complement the canonical coupling mode even in domain-swapped Kv channels (57, 88, 89). An increasing number of ion channels have been discovered to adopt non- domain-swapped TM topology besides BK, including hERG (90), KvAP (91), HCN (92), Eag1 (93), and CNG channels (94). It is possible that the mechanistic features observed for BK channels may apply to voltage gating of the important emerging class of non-domain swapped ion channels in general.

Importantly, full-length BK channels can be independently activated by membrane potential and intracellular Ca²⁺ (95–97). The latter requires the C-terminal gating ring that is absent in the Core-MT construct studied in this work. Cryo-EM structures have revealed functional state-dependent interactions between the gating ring and VSD (28–30, 38), strongly supporting that the calcium sensor and VSD are coupled. It is not clear how the VSD-gating ring coupling may affect the nature of VSD movements or how VSDs may drive the pore opening. Nonetheless, mechanistic details revealed from the voltage-gating of the Core-MT construct should provide a solid basis for future computational and experimental studies of BK activation and regulation.

Methods and Materials

Molecular modeling and atomistic simulations

The structure of the Core-MT human BK channel in the closed state was derived from the Cryo-EM structures of the acBK channel in deactivated Ca²⁺-free states (PDB 5tji (29)) as previously described (32, 98). The simulated construct was truncated at R342. The 11-residue C-terminal Kv mini-tails are not involved in voltage gating (99, 100), and they were thus not included. The dynamic loop (C54-V91) and N-terminal tail (M1- N19) were not included either. Residues before and after the missing segments are capped with either an acetyl group (for N-terminus) or a N-methyl amide (for C- terminus). Standard protonation states under neutral pH were assigned for all titratable residues.

The initial structures was first inserted in model POPC lipid bilayers and then solvated in TIP3P water using the CHARMM-GUI web server (101). Even though polarizable water models are probably necessary to capture the precise energetics (and kinetics) of dewetting transitions within a protein cavity, it has also been shown that classical nonpolarizable water models such as TIP3P is sufficient to capture the spontaneous dewetting and rehydration of BK channels (32–34, 36). In fact, a recent analysis using the Drude polarizable force field actually under-estimated the hydration free energy of deactivated BK pore (∼0.5 kcal/mol) (35), which is very likely too small to sustain a dry and nonconductive pore. All systems were neutralized and 150 mM KCl added. The final simulation boxes contain about 593 lipid molecules (POPC) and ∼50,000 water molecules and other solutes, with a total of ∼250,000 atoms and dimensions of ∼160 × 160 × 110 ų. The CHARMM36m all-atom force field (102) and the CHARMM36 lipid force field (103) were used. All simulations were performed using Desmond (104) on Anton 2 (70, 71) or CUDA-enabled versions of Gromacs 2020 (105, 106) on GPU clusters. Electrostatic interactions were described by using the Particle Mesh Ewald (PME) algorithm (107) with a cutoff of 12 Å. Van der Waals interactions were cutoff at 12 Å with a smooth switching function starting at 10 Å. Covalent bonds to hydrogen atoms were constrained by the SHAKE algorithm (108), and the MD time step was set at 2 fs. The temperature was maintained at 298 K using the Nose-Hoover thermostat (109, 110) (in Gromacs). The pressure was maintained semi-isotopically at 1 bar at membrane lateral directions using the Parrinello–Rahman barostat algorithm (111).

All systems were first minimized for 5000 steps using the steepest descent algorithm, followed by a series of equilibration steps where the positions of heavy atoms of the protein and lipid were harmonically restrained as prescribed by the CHARMM-GUI Membrane Builder (112). For Anton simulation, additional 70 ns equilibration step with 0.1 kcal.mol^-1.Å^-2 position restraints on all protein heavy atoms was performed before production runs. Note that the inner pore became dewetted during equilibration simulations. All production simulations were performed under NVT (constant particle number, volume and temperature) conditions at 298 K. The P-loop/filter (T273 to D292) and C-terminus were harmonically restrained with a small force constant of 0.1 kcal.mol^-1.Å^-2 in all production simulations to prevent drift in all production simulations. Using Anton 2, a 2-μs control simulation was first performed without membrane voltage (sim 1, Table S1). Two independent 10-μs simulations were then performed at 750 mV to directly probe voltage-driven activation transitions (sim 2). The voltage was applied as a constant external electric field (E = V/Lz) imposed along the z dimension (84).

Steered molecular dynamics simulations

Steered molecular dynamics (SMD) simulations were performed to test the role of R210 and R213 in voltage sensing and channel activation. In these simulations, a moving reference point was used to simulate the effect of an external electric field by pulling the CZ atoms of R210 and R213 from their initial positions in the resting state to those observed in the activated state. As illustrated in Figure S4, the initial reference points of SMD were positioned at the location of each CZ atom of R210 and R213 in the resting/closed state. From 0 to 100 ns, these reference points were moved along the z-axis at constant speed until they reached the positions corresponding to the predicted activated VSD state at 100 ns. A harmonic positional restraint of 5 kcal/mol·Å² was applied between each CZ atom to its respective reference point along z-axis, allowing the CZ atom to track the reference point’s movement. After 100 ns, the reference points were held stationary at the positions corresponding to the activated VSD state, with the same harmonic restraints imposed to maintain the CZ atoms in their activated position. Four independent simulations (sim 8) were performed, and pore opening was observed in two of the four replicates (replicas 2 and 4). In the other two replicas, the S5-S6 packing was broken during the first 100 ns, likely due to strong steering potentials imposed, and the simulations were terminated.

Free energy analysis of gating charges

The total gating charge ΔQ is the sum of contributions of all charged residues in VSD and can be calculated as (83, 84):

where is q_i is the charge of a specific residue, and the dimensionless quantities f_s(i) for a specific conformational state s, which can be open (o) or closed (c), represent the coupling of charge qi to the transmembrane potential. The value of fs(i) is derived from the difference between charging free energies calculated at two different voltages V1 and V2,

where ΔG_s(V, q_i) is the free energy cost of increasing the charge of residue i from 0 to qi in state s under membrane voltage V. ΔG_s(V, q_i) was calculated using thermodynamic integration (TI), with the charge gradually scaled to the final value in increments of 0.1, i.e., λ = 0, 0.1, 0.2, …., 1.0. In this work, we focused on 7 charged residues on VSD, namely, D153, R167, D186, R207, R210, R213 and E219, and performed four TI free energy calculations for each residue (sim 3-6, Table S1). The last snapshot of simulation 2b was used to represent the activated state of VSD. For each TI window, the system was first equilibrated 200 ps, and then simulated for 4 ns to calculate the mean force <δH/δλ>λ.

Structural and electrostatic potential analysis

The displacement of residue charged groups, as well as movement of helixes (translocation, titling and rotation) is calculated using MDanalysis (113) together with in-house scripts. The TM helices are defined as: T109-S135 (S1), F148-A170 (S2), V181-L199 (S3), G205-N225 (S4), S230-S259 (S5), and A313-E324 (S6 after the glycine hinge). Residue contacts were identified using a minimal heavy atom distance cutoff of 5 Å. Pore water molecules were identified as those occupying the inner pore cavity below the selectivity filter, roughly from L312 to the plane defined by the center of mass of P320. Pore profiles are calculated using program HOLE(114). To examine how the protein modulates the local electric field, average electrostatic potential maps were then calculated from the first (resting) and last (activated) 250 ns of sim 2b using a combination of Python scripts and the VMD’s PMEPot and VolMap plugins (115). The electrostatic potential is obtained by solving the Poisson equation where the sum runs over all atoms, and ρi(r) is the charge electric potential distribution contributed by atom i at position r approximated by a spherical Gaussian: . An Ewald factor β of 0.25 Å^-1 was used for the inverse width of the Gaussian and grid size was set to 1 Å. All molecular illustrations were prepared using VMD (116). External electric potential was added to the PMEPot grid file with in-house scripts (see supplementary files).

Dynamic community and coupling pathways

Dynamic network analysis was performed using the Networkview (117) plugin of VMD . To build the network, each amino acid was represented as a single node at the Cα position, and a contact (edge) was defined between two nodes if the minimal heavy- atom distance between residues was within a cutoff distance (5 Å) during at least 75% of the trajectory. The resulting contact matrix was weighted based on the covariance of dynamic fluctuation (Cij) calculated from the same MD trajectory as wij = - log(|Cij|). The length of a possible pathway Dij between distant nodes i and j is defined as the sum of the edge weights between consecutive nodes along this path. The shortest path, calculated using Floyd-Warshall algorithm (118), is considered the optimal pathway with the strongest dynamic coupling. Suboptimal paths are identified as alternative top-ranked paths with lengths that deviate by less than 50% from the optimal path. We further performed dynamic community analysis to identify groups of residues are more tightly coupled among themselves (119).

Information flow analysis

Information flow provides a global assessment of the contributions of all nodes to the dynamic coupling between selected “source” and “sink” nodes(120, 121). This analysis complements the dynamic pathway analysis to provide additional insights on how different residues may contribute to sensor-pore coupling. For this, a network similar to the one described above was first constructed. Pairwise mutual information was calculated between node i and j as follows: Mij = Hi + Hj − Hij. Hi is calculated as , where ρi(x) is the fluctuation density and x is the distance to the equilibrium position(121). Gaussian mixture model (GMM)(¹²²) is used to estimate the density. The residue network is then defined as Aij = Cij Mij, where Cij is the contact map. To analysis the information flow from the source (S0) to sink (SI) nodes, the network Laplacian, is calculated as L = D - A, where D is diagonal degree matrix: . The information flow through a given node (residue) is defined as . The potentials P is given by P = L̃^-1 b, where L̃^-1 is the inverse reduced Laplacian, and b is the supply vector that corresponds to one unit of current entering at the source node that will exit at sink nodes. The magnitude of fi thus quantifies the contribution of residue i to dynamic coupling between the source and sink nodes.

Data availability statement

All data needed to evaluate the conclusions in the paper are present in the paper and the Supplementary Materials. The PDB structures of Core-MT with resting and activated VSDs, as well as simulation input and analysis script, can be found on GitHub at: https://github.com/mdlab-um/Votage_gating_Core-MT.

Supplementary Information

Supplementary Tables

Table S1.

Supplementary Figures

Figure S1

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Acknowledgements

The authors thank Jianmin Cui and Guohui Zhang for critical discussions. This work is supported by NIH R35 GM144045 (Chen). Anton 2 computer time was provided by the Pittsburgh Supercomputing Center (PSC) through Grant R01 GM116961 from the National Institutes of Health. The Anton 2 machine at PSC was generously made available by D.E. Shaw Research.

Additional information

Author contributions

Chen and Jia conceived and initiated the study. Jia performed simulations and analysis. Jia and Chen wrote the manuscript.