Multiscale modeling shows that dielectric differences make Na_V channels faster than K_V channels

The most important role of voltage-gated sodium (Na_V) and potassium (K_V) channels is to generate action potentials, all-or-none (i.e., binary) rapid and transient changes of the electric potential difference across the plasma membrane, that propagate signals long distances in the nervous system. This ability is essentially due to the different timing with which these two channels respond to a potential change (Fig. 1 A); Na_V channels activate rapidly upon an initial depolarization of the membrane, allowing a Na current to enter the cell and induce a further membrane depolarization along the nerve during the first millisecond or so after the initial trigger. By contrast, K_V channels remain essentially closed during this initial interval and open only later, causing the repolarization of the membrane to the resting condition (Bezanilla et al., 1970; Hodgkin and Huxley, 1952; Rojas et al., 1970).

The molecular and physical basis for the different kinetics of Na_V and K_V channels has long been sought. The recently determined crystal structure of the voltage sensor domain (VSD) of K_V and Na_V channels has allowed investigation of the structure-function relationship of the voltage sensor, its environment, and the forces controlling its movement. In both channels, the VSD is formed by the transmembrane segments 1–4 (S1–S4), with the S4 segment containing positively charged residues every third amino acid and forming the voltage sensor, the structure that senses the electric field applied across the membrane and moves upon depolarization (Long et al., 2007; Long et al., 2005; Pan et al., 2018; Payandeh et al., 2011).

The S4 segment is surrounded for most of its length by wide intracellular and extracellular vestibules, where water and ions can easily access and balance the gating charges (Fig. 1 B). The two vestibules are separated by a series of mostly nonpolar residues that form a region, the gating pore, thought to be hardly accessible to water (green and red residues in Fig. 1 B). Several mutagenesis experiments suggest that some of these residues are in direct contact with the moving gating charges and modulate the voltage-dependent properties of the channels (Lacroix et al., 2012, 2013, 2014; Lacroix and Bezanilla, 2011; Tao et al., 2010).

One of the residues of the gating pore that determine the slow kinetics of K_V channels, compared with faster Na_V channels, is isoleucine at position 287, located in the middle of the S2 segment (red in Fig. 1 B; notice that in this paper we will always use the residue numbering of the Shaker channel). This nonpolar, highly hydrophobic residue is found in many K_V channels. It is conserved by evolution presumably because of its essential role. By contrast, in the S2 segment of the first three domains of fast Na_V channels, at the corresponding position of isoleucine in K_V channels, we consistently find the uncharged but much more hydrophilic residue threonine (Fig. 1 C; see also Fig. S1 for a more extended sequence alignment; Lacroix et al., 2013). This seems to be an example of divergent side-chain evolution, in which two related proteins (Na_V and K_V channels) have selected two different side chains as an adaptation to create different dielectric polarization and different activation kinetics, and thus to create a propagating action potential. Notably, when isoleucine at position 287 is mutated to threonine (or to other polar amino acids), the Shaker gating currents become much faster and similar to the activation rate of Na_V channel gating currents (Lacroix et al., 2013, 2014). These results strongly suggest that the physical structure and the underlying mechanism controlling the activation rate of voltage-gated ion channels involve the hydrophobic strength (i.e., the polarizability) of the gating pore, especially of the residue at position 287. They suggest that a highly nonpolar, hydrophobic gating pore creates an energetically unfavorable condition, a high-energetic barrier (mostly a dielectric barrier) for the gating charges to move through upon depolarization, resulting in a slowing of the voltage sensor activation (Fig. 1 D). In this view, increasing the polarizability of this region by replacing a highly nonpolar side-chain residue with a more polar amino acid lowers the dielectric energy barrier, provides counter charge (in the wall of the gating pore) for the gating charges, and ultimately increases the activation rate of the voltage sensor.

We test this mechanistic hypothesis using a multiscale modeling approach (i.e., an approach that investigates the system at two different time and space scales) in a series of hierarchical computational methods whereby the calculated quantities from a computational simulation at an atomistic scale are used to define the parameters of the model at a coarser, macroscopic scale that satisfies conservations laws (Eisenberg, 2018). In the specific case studied here (voltage-dependent gating), the idea translates into using Brownian modeling to macroscopically describe the system and obtain gating currents and then using MD simulations to determine the needed free parameters present in the Brownian model. It is difficult to determine these parameters experimentally, but there is hope from the exciting work of the Boxer group (Fried et al., 2014; Fried and Boxer, 2017; Suydam et al., 2006). The multiscale approach appears feasible, because assessing free parameters through MD simulations is much less susceptible to artifact and is less computationally demanding than simulating the overall movement of the voltage sensor. It also has the advantage of using a Brownian model that automatically satisfies conservation laws of mass, current, and charge and Poisson’s equations. It is difficult to satisfy Poisson’s equations and conservation of current when using the periodic boundary conditions common in MD.

The VSD used in the simulation was taken from the Shaker model of Henrion et al. (2012) in the activated state. For MD simulations, the protein was embedded on a membrane environment made of POPC, with dimensions of 80 × 80 Å, and solvated in two steps. The distance between the maximum and minimum z coordinates and the water box edges in the z axis was set to 15 Å. The system was first solvated below the membrane plane with a water box of dimensions 83.46 Å, 80.05 Å, and 16.19 Å (x, y, and z, respectively), where 9,462 TIP3 water molecules (Jorgensen et al., 1983) were placed. Then, the system was solvated above the membrane plane with a water box of dimensions 83.46 Å, 80.05 Å, and 20.02 Å (x, y, and z, respectively), where a total of 9,462 TIP3 water molecules were placed. After the addition of the water molecules, the system presented a total net charge of −2. The system was balanced, and a total of 27 Cl ions plus 29 Na ions were added to the system, making up a salt concentration of 0.15 mol/liter. The MD simulations in the present study were performed using the NAMD MD package (Phillips et al., 2005). The CHARMM36 force field (Best et al., 2012; MacKerell et al., 1998) was used in all MD simulations. A distance cutoff of 12.0 Å was applied to short-range, nonbonded interactions and a cutoff of 10.0 Å was used for the smothering functions. Long-range electrostatic interactions were treated using the particle-mesh Ewald (Darden et al., 1993) method. Before the MD simulations, all systems were submitted to an energy minimization protocol for 2,000 steps. After minimization, an annealing of 0.29 ns was performed, with a temperature ramp of 0.24 ns of simulation, where the temperature was raised from 60 K to 300 K. The pressure was maintained at 1 atm using a Nosé-Hoover Langevin piston (Feller et al., 1995; Martyna et al., 1994). The equations of motion were integrated using the r-RESPA (reversible reference system propagator algorithm) multiple time-step scheme (Phillips et al., 2005) to update the short-range interactions every one step and long-range electrostatics interactions every two steps. The time step of integration was chosen to be 2 fs for all simulations. During annealing, the motion of the backbone of the protein was restrained. After annealing, a 10-ns equilibration was performed while maintaining the restriction of the backbone motion, during which the temperature was maintained at 300 K using Langevin dynamics. The pressure was maintained at 1 atm using a Nosé-Hoover Langevin piston. A distance cutoff of 12.0 Å was applied to short-range, nonbonded interactions and a cutoff of 10.0 Å was used for the smoothing functions. Long-range electrostatic interactions were treated using the particle-mesh Ewald method. Finally, we performed unrestricted MD simulations by applying different values of the electric field in the z direction. These simulations were used to assess the running average of total dipole moment, which was used to estimate the dielectric constant as described in the Results section.

A full description of the Brownian model of voltage gating can be found elsewhere (Catacuzzeno and Franciolini, 2019; Catacuzzeno et al., 2020). Briefly, in our model, the VSD was approximated by an hourglass structure consisting of a cylindrical gating pore with a length of 7 Å and a diameter of 10 Å flanked by internal and external water-accessible (WA) vestibules with a length of 3.15 nm each and a conical shape opening with a half angle of 15° into two hemispherical subdomains of bath solution, both having a radius of 1 µm. The S4 charge profile (Z_S4) was built by considering the six positive charges, whose mean distance between the charged atoms was determined from a 3-D Shaker channel model, with each giving rise to a charge profile normally distributed with a standard deviation of 0.1 nm. The permanently charged profile (Z_F), which in this case deserves the name “fixed” because it does not move, was similarly built by considering the vertical dimension of the C_α of all the charged residues on segments S1–S3 of the VSD.

The energy profiles are assessed self-consistently by computing (and not assuming) electrostatics. More specifically, for each point of the reaction coordinate (i.e., for each position of the voltage sensor along its activation pathway), the corresponding associated energy was assessed by solving iteratively the flux conservative equation for the ionic distribution in the baths and vestibules and the Poisson's equation for the voltage profile until reaching a steady-state (or equilibrium, for equal bath ion concentrations) solution. Once steady state (or equilibrium) is reached, the resulting voltage profile is used to assess the electrostatic energy by integrating the charge on the voltage sensor multiplied by voltage over the entire space domain (cf. the mathematical relationship in the legends of Figs. 5 and 7; note that Fig. 7 reports the sole electrostatic energy, while Fig. 5 reports the total energy experienced by the voltage sensor, including the mechanical energy due to the attached spring). Obviously the calculation of the energy profile as a function of the voltage sensor position does not involve the solution of the Langevin (or the Fokker–Planck) equation, since in our model these equations are used to assess the movement of the voltage sensor (e.g., during the evolution of the gating currents), and because we choose the position of the voltage sensor while assessing the energetic profile, as it represents the independent variable of the plot. Thus, the electrostatic energy profile shown in Fig. 7 essentially comes from the solution of the Poisson–Nernst–Planck equations while varying the position of the voltage sensor. It is important to note that the reliable output of our calculations is the shape of the gating current and other associated variables. The potential profiles are computed to make our results more understandable and to relate to traditional views of channels. The variable nature of the potential profile is not included in this classical view or in our calculation of potential profiles. Much of the variable nature of the potential profile is included in our solution of the Langevin equations and gating currents.

In our model, we assume ionic steady-state equilibrium (i.e., the overdamped system established mathematically by, for example, Eisenberg et al., 1995) that was derived rather than assumed using the overdamped high friction limit. This assumption was also used in Horng et al. (2019), justified by the fact that the voltage sensor moves much more slowly than the ionic re-equilibration in the baths and vestibules. Under this assumption, the energetic profile shown in Fig. 5 (i.e., total energy = electrostatic + mechanical) is exactly the energy encountered by the voltage sensor during its movement and the energy used to drive the voltage sensor movement through the Langevin (or Fokker–Planck) equation. Mathematically speaking, the spatial gradient of the same energy profile represents the force acting on the voltage sensor, F_ex, present in Eq. 3 (the Langevin equation) and Eq. 4 (the Fokker–Planck equation).

We emphasize that the main outputs of this paper are the actual calculations of charge movements, potential fields, effective dielectric constants, and most importantly the functional output of the system, the gating currents themselves. The potential profiles were used to provide a comforting link with the past in Figs. 7 and 8. They are most helpful in a human sense but are likely to be misleading when the more complex interactions described in the last paragraphs are important.

Fig. S1 shows conservation of the speed-control residue in the S2 segment in different types of fast and slowly activating Na_V channels and slowly activating K_V channels. Fig. S2 shows an extended sequence alignment for Na_V and K_V channel voltage sensors and additional information regarding the water accessibility of the VSD in the I287T Shaker channel.

MD simulations performed on available crystal structures can be used to predict the polarizability of a protein or part of it and thus its effective dielectric constant (Guest et al., 2011; King et al., 1991; Patargias et al., 2010; Simonson and Perahia, 1996; Smith et al., 1993; Voges and Karshikoff, 1998; Warshel et al., 2006), a parameter needed in the coarser-scale Brownian modeling of voltage gating because it determines the amount of polarization charge available to balance movements of permanent charges, like those of the gating particle. We thus used MD to assess whether the substitution of threonine (found in the Na_V channel) in place of isoleucine (found in the K_V channel) at position 287 can sensibly change the polarizability and consequently the dielectric constant, ε, of the gating pore and ultimately the movement and stability of the gating charges inside it.

As a first step, we determined the water accessibility of the VSD in the voltage-gated Shaker K channel by performing an MD simulation using the activated state of the Shaker VSD model (Fig. 2). After inserting the protein in a phospholipid bilayer and water box and performing the minimization, annealing, and equilibration of the overall system, we ran a 100-ns MD simulation to look for water inside the vestibules and the gating pore. At the end of the simulation, we looked at the position of water molecules inside the VSD to determine whether water could also access the gating pore. The structure on the left of Fig. 2 A shows a frame at the end of the simulation, while the one on the right shows the WA space of the VSD obtained by the superposition of the water molecules from 60 consecutive frames (one every 0.1 ns). These pictures show that water readily enters both intracellular and extracellular vestibules. The pictures also identify two distinct regions inside the gating pore with different water accessibility. One region, in direct contact with the intracellular vestibule, is totally inaccessible to water, thus forming an effective water-inaccessible (WI) region. This region is very narrow (∼2–3 Å), as judged by the distance between the closest water molecules in the two vestibules (Fig. 2 D). It is located at the level of the gating pore in adjacency to residues F290 on S2 and V236 and I237 on S1 (Fig. 2 E). The second region, located between the WI region and the extracellular vestibule, including the other seven hydrophobic residues forming the gating pore (including residue I287; Fig. 2 E), is WA. As shown in Fig. 2 B, the equilibrated simulation shows three water molecules on average in this region (defined as a cylinder of 5 Å radius and extending from −1 Å to 6 Å with respect to the center of mass of the F290 aromatic ring; cf. inset to Fig. 2 C). These results show the gating pore has two regions with different accessibility to water, a WI region totally inaccessible to water and an adjacent WA region containing three water molecules most of the time. Due to the different water accessibility, we expect these two regions to have a quite different effective polarizability, with a relative effective dielectric constant close to 4 (the value typical for the interior of proteins) for the WI region but sensibly higher for the WA region, since here, water dipoles could significantly contribute to the stabilization of charges.

The dielectric constant has two main components, one given by the polarization of the electronic clouds, $εel,$ and the other due to the reorientation of permanent dipoles along the applied electric field, $εpd.$ Of these, $εel$ can be measured experimentally (in bulk solutions) by impedance spectroscopy (Macdonald, 1992) or by the refractive index of the material (Gussoni et al., 1998), which depends strictly on the dielectric constant at the appropriate frequency. In addition, the $εel$ of a protein region may simply be crudely estimated by summing up the individual polarizability levels of the various amino acids present that can be obtained experimentally (Voges and Karshikoff, 1998). There is regrettably no experimentally recognized method to assess $εpd$ inside a protein region, although it is known to depend on the reorientation of dipoles along the electric field. Indeed, the assignment of polarizability in an atomic scale system has serious ambiguities, as described in Purcell and Morin (2012). For a small spherical volume of material surrounded by an infinite and homogeneous medium, $εpd$ can be estimated by running extended MD simulations and analyzing the equilibrium fluctuation of the total dipole moment of the sample volume by applying the Kirkwood–Frohlich theory of dielectrics (Frohlich 1949; Simonson and Perahia, 1996). Unfortunately, the structure we are studying—the voltage sensor—is immersed in a highly inhomogeneous medium (water and protein residues) on which this approach becomes impractical. In this case, the dielectric constant in the small sample volume in the protein interior will also depend on the polarizability of the surrounding medium, whose value is neither known nor can be assumed to be homogeneous. The assumption of homogeneity is crucial because inhomogeneity creates dielectric boundary charges that can easily dominate the properties of a system (Nadler et al., 2004; Nadler et al., 2003). One can hardly imagine a more inhomogeneous environment than that of a gating pore embedded in a dielectric membrane immersed in an ionic solution.

We then applied this method to estimate the effective dielectric constant of the WI and WA regions of the gating pore, both in the WT and in the mutant structure in which isoleucine 287 (found in K_V channels) had been replaced by threonine (mutant I287T) found in Na_V channels. In the calculation of the dielectric constant of the WI region, we considered all three residues present in this region, namely residues F290, V236, and I237, since they all surround the gating charges while they move through this region of the gating pore (Fig. 3 D, green residues). As expected for a very hydrophobic, nonpolar region, the polarizability of the WI region was electronic in nature, with a value ∼4.0, irrespective of the residue present at position 287 (Fig. 3, C and E).

The situation is quite different for the WA region, where the gating pore widens, and not all residues are close to the moving gating charges. In the assessment of the effective dielectric constant, we thus considered only those residues adjacent to the gating charge (i.e., placed at a distance less than 3 Å), namely, residues I287 and S240, plus few water molecules (Fig. 3 D, blue residues and red/white water molecules). The dielectric constant estimated in this region was much higher than in the WI region, but, most relevantly, it was quite different in the two WA structures being considered (∼20 in the WT and ∼60 in the I287T mutant; Fig. 3, C and E). These results indicate that the nature of the residue at position 287 critically contributes to the polarizability of the WA region around the gating charge, and a mutation of hydrophobic isoleucine at this position introducing a more hydrophilic amino acid like threonine markedly increases the polarization charge in this region.

Two main reasons explain why the WA region has a higher polarizability in the mutant compared with WT. First, the threonine residue has a much higher permanent dipole because of its hydroxyl group, and thus stronger hydrophilic power than isoleucine. As a result, the dipole moment contributed by this one residue is much higher in the I287T mutant than in the WT structure (Fig. 4 A). Second, the region close to the gating charge R371 in the I287T mutant appears to have a slightly higher water content than in WT (Fig. 4 B), presumably because threonine stabilizes nearby water molecules with hydrogen bonds. However, the overall water content inside the gating pore was the same in WT and mutant structures (Fig. S2).

The macroscopic scale description of voltage gating is based on the Brownian model that we have recently developed and is a good predictor of the major signature features of experimental gating currents of Shaker channels (Catacuzzeno and Franciolini, 2019; Catacuzzeno et al., 2020; see Materials and methods). Other papers (Horng et al., 2019; Kim and Warshel, 2014; Peyser et al., 2014; Peyser and Nonner, 2012b, 2012a) give other views of the voltage sensor and resulting currents at different resolutions. Only Horng et al. (2019) and Catacuzzeno and Franciolini (2019) calculate gating currents comparable to those measured experimentally under a range of conditions, with Horng et al. (2019) having less atomic detail. Our model takes the relevant structural and electrostatic parameters directly from the 3-D structural model of the Shaker VSD (Fig. 5, A and B). The electrodiffusion of ions in the bath is described through a flux conservation equation that assumes a quasi-equilibrium, on the time scale of gating currents, and the motion of S4 segments is described by Brownian dynamics. The electric potential profile is modeled using the Poisson equation and includes all the charges present in the system, including the gating charges on the S4 segment and the fixed charges on segments S1–S3 of the VSD. Notably, the Poisson equation that is used in our model to assess the electric potential profile contains the relative dielectric constant, ε, as an undetermined parameter that we have estimated above for both the WT and mutant gating pores using MD simulations. It should clearly be understood that this form of coarse graining has the significant advantage that the coarser, lower-resolution treatment satisfies conservation laws of mass, charge, and current. Coarse graining procedures that do not explicitly introduce the Poisson equation and boundary conditions are unlikely to satisfy conservation laws of charge and current. Periodic boundary conditions in traditional MD make it difficult to satisfy Poisson’s equations and conservation of current in a spatially inhomogeneous, nonperiodic system like an ion channel protein or gating pore.

In our previous model, we assumed a 2-Å-long gating pore localized at the level of F290 that separated electrically the intracellular and extracellular solutions. The effective dielectric constant was set to 4 in this short region and to 80 elsewhere, including the adjacent region containing residue I287 (Catacuzzeno and Franciolini, 2019). The gating pore used in the present study has been changed to account for our MD results shown above. More specifically, the gating pore is now pictured as made of two distinct regions. One region has all the properties previously described, with 2 Å length comprising F290, and is fully impermeant to water (green bar in Fig. 5, A and B). The second region (5 Å long), comprising several hydrophobic residues, including I287, is accessible to water and has now been assigned a dielectric constant of 20, in accordance with the MD results (Fig. 5, A and B). As shown in Fig. 5 C, the updated model continues to predict the major features of experimental gating currents obtained in response to a depolarizing pulse. First, a very fast gating current component is present at the beginning of the depolarizing step, rising instantaneously and then falling very rapidly. Second, a later slower current component develops, starting with a plateau/rising phase and proceeding with a slow decay. Both these features have been observed in experiments (Bezanilla, 2000; Sigg et al., 2003). The plot in Fig. 5 C represents the total energy profile assessed for the various allowed positions of the voltage sensor. It can be seen that the model predicts five energy wells corresponding to the five different positions the S4 segment can assume in accordance with many studies showing that the voltage sensor moves in multiple steps during the activation process (Bezanilla et al., 1994; Henrion et al., 2012; McCormack et al., 1994; Schoppa et al., 1992; Tytgat and Hess, 1992; Zagotta et al., 1994). In our model, each step represents the passage of a different gating charge through the high-resistance and WI region of the gating pore (compare drawings in Fig. 5 C).

We verified that a different dielectric constant of the WA region can explain the different kinetics of the macroscopic gating current of WT and I287T mutant found experimentally (Lacroix et al., 2013). We used our macroscopic Brownian model of voltage gating to see whether changing the dielectric constant of the WA region from 20 to 60 reproduces the increased rate of gating current seen in mutagenesis experiments.

As shown in Fig. 6, B and C, there was a significant increase in the rate of change of the gating currents evoked by depolarizing pulses (ON gating currents), especially at intermediate levels of depolarization (Fig. 6 C), as observed experimentally. The model also predicted an increased rate of the OFF gating currents induced by repolarization after a short depolarizing pulse (Fig. 6, E and F). No change was observed in the charge versus voltage relationship (Fig. 6 D). Overall, these results are very similar to those observed experimentally with the mutation I287T, suggesting that a change in the dielectric constant within the gating pore is a major determinant of the functional effects induced by the mutation.

We then examined the effects of changing the polarizability of the WA region on the electrostatic energy profile and found that the mutation selectively lowers the first of the four peaks present in the profile (Fig. 7). This result appears quite reasonable, since the first energetic peak essentially represents the dielectric boundary energy commonly (but incorrectly) called the Born energy for the addition of one (gating) charge to the WA region of the gating pore (Nadler et al., 2004, 2003), and we expect this quantity to depend on the dielectric constant present there (compare inset to Fig. 7)¹. The remaining energy peaks are instead essentially unmodified as the voltage sensor continues to move. This presumably occurs because the charge residing in the WA region is simply replaced by the next charge of the S4 segment, as the voltage sensor moves in response to an applied voltage, keeping the overall charge content in the WA region constant.

The energy profile presented in Fig. 7 shows that the I287T mutation, besides lowering the first energy barrier, also reduces (by approximately the same amount) the energy of the neighboring well, corresponding to the first gating charge sitting close to the WA region of the gating pore. If, intuitively, a reduction of an energy barrier leads to a faster activation time course while a reduction of a local minimum is instead associated with a slowing of the voltage sensor activation, then the overall effect of simultaneous lowering an energy barrier and nearby energy well is difficult to foresee and could result, in principle, in either speeding or slowing the voltage sensor movement, depending on which of the two effects will prevail.

Based on our modeling results, we conclude that the decrease of the first energy barrier is by far the most important effect of the I287T mutation. To verify this conclusion, we modeled the effect on gating currents kinetics of altering selectively either the first peak or the second well of the energetic profile by amounts comparable to those observed following the I287T mutation (Fig. 8). The results of these computations show that the effect of altering the energy barrier is much more pronounced than the effect of altering the neighboring local minimum and thus explain why the mutation results in a net increase in the gating current kinetics.

Several conclusions may be drawn from the results of our multiscale analysis. First, using MD simulations, we have shown that the gating pore of the Shaker channel, defined in terms of the 10 residues proposed by Lacroix et al. (2014), does not have uniform accessibility to water. More specifically, we have found that the region totally inaccessible to water is actually very short, with a length of 2–2.5 Å, and includes the three most intracellular residues, namely F290 on the S2 segment and V236 and I237 on S1. The remaining part of the gating pore, covering a length of ∼5 Å, is by contrast accessible to water, with three water molecules often found in this region, as shown by our MD simulations. These results appear in agreement with previous MD simulations performed on the VSD of other voltage-gated ion channels. More specifically, an almost continuous arrangement of water molecules was found in the VSD of the KvAP channel (a bacterial K_V channel from Aeropyrum pernix), with a very short interruption at the level of the R133–D62 salt bridge (Freites et al., 2006). A similar conclusion was reached in MD simulations of the Shaker homologue K_V1.2, where a very short constriction, localized at the K5–D259 salt bridge and forming the so-called gating charge transfer center together with F233 (the homologue of F290 in Shaker), prevented communication between the upper and lower crevices (Delemotte et al., 2010; Tao et al., 2010). Finally, a short 2-Å WI region centered on a conserved phenylalanine was also found in the VSDs of the first three domains of Na_V channels and thought to be functionally similar to the VSD of the K_V channel (Gosselin-Badaroudine et al., 2012). The potential difference applied to the membrane will then drop mostly over the very short WI region, because this part of the gating pore is the only region with very high resistance to water and ions. The notion of a very narrow high-resistance region that focuses virtually the entire transmembrane voltage drop has already been suggested by several studies (Ahern and Horn, 2005; Asamoah et al., 2003; Islas and Sigworth, 2001). The focusing of the electric field is of course incompatible with the idea of “constant fields” that has been a mainstay of biophysics for more than 50 years and remains so in teaching materials to this day. Electric fields change shape as conditions change because of the fundamental properties of the electric field, as pointed out in the context of channels some years ago (Eisenberg, 1996a). F290 has been proposed to act as a gating charge transfer center (Tao et al., 2010), although its role has been later refined to stabilizing the activated VSD conformation (Lacroix et al., 2014; Lacroix and Bezanilla, 2011). In our view, this residue is important for creating the only very-high-resistance region inside the VSD, where most of the voltage drops. It is in the 2-Å-wide region that includes F290 that the moving gating charges experience the electrostatic force that stabilizes either the closed or open conformation, depending on the applied voltage.

A second result of this study is the quantitative demonstration that the polarizability (i.e., the hydrophobicity) of the residue at position 287 of Shaker channels can greatly modulate the field strength acting on the gating charge as it enters the gating pore. We showed this by estimating the dielectric constant of the region surrounding the gating charge, using MD simulations at variable applied electric fields, and assessing the change in the total dipole moment of the region. In contrast to the computational methods usually employed to assess the dielectric constant (Guest et al., 2011; King et al., 1991; Simonson and Perahia, 1996; Smith et al., 1993; Warshel et al., 2006), which are based on the Kirkwood–Frohlich theory of spatially uniform dielectrics (Frohlich, 1949; Simonson and Perahia, 1996), our method does not assume a homogeneous environment of known polarizability surrounding the sample volume of interest. For this reason, our approach can be applied also to restricted regions inside a protein containing a small number of side chains. This method has been already successfully applied to estimate the dielectric constant of a homogeneous liquid, as well as of water confined inside the pore of a synthetic ion channel (Kolafa and Viererblová, 2014; Riniker et al., 2011; Sansom et al., 1997; Xu et al., 1996; Yeh and Berkowitz, 1999). However, to the best of our knowledge, it has not been applied before to assess the local dielectric constant inside a (channel) protein.

We found that changing the residue at position 287 significantly changes the effective dielectric constant acting on the gating charge. More specifically, when the isoleucine found at this position in Shaker (as well as in other K_V channels) is changed to threonine, the residue found at the corresponding position in Na_V channels, the estimated dielectric constant appears to increase by as much as three times. The observed increase may have two components: (1) the much higher permanent dipole present on the threonine side chain (as compared with isoleucine, mostly because of the hydroxyl group there), which contributes directly to increase the polarizability, and (2) the possibility that threonine, as a polarizable dipole, attracts more waters inside the WA region of the gating pore, electrostatically and by establishing weak hydrogen bonds (Desiraju and Steiner, 1999). In summary, our data suggest that the WA region of the Na_V channel has a higher effective dielectric constant, and this possibly allows the first gating charge to enter the gating pore more easily than in K_V channels. The change in dielectric constant lowers the energetic barrier gating charge encounters in entering the gating pore.

We are aware that our model has limitations. One is related to the estimation of the local dielectric constant inside the gating pore, which we have obtained with classical MD. As there is no electronic polarization with MD, it might be possible that the resulting dipole moment responses are overestimated and so would be the dipolar component of the dielectric constant. The estimation of the dielectric constant could also be biased by our assumption of a linear relationship between the applied field and the total dipole moment, which we know to be valid in water for moderate fields. We do not know whether this assumption is valid in the complex dielectric environment of the voltage sensor. A second limitation relates to our MD simulations in two respects. They were run in a limited spatial environment using periodic boundary conditions and might have produced time-varying electric fields as the dimensions of the box changed. In addition, their limited time span (100 ns) might not be sufficient to cover wetting/dewetting transitions, which also have an impact on the estimation of the dielectric constant. Nonetheless, in spite of these limitations, we believe that our study provides an important testable hypothesis and opens new avenues for the field of voltage-dependent gating.

We also used our Brownian model (Catacuzzeno and Franciolini, 2019) to see if the predicted change in the effective dielectric constant found by MD explains the functional effect of the I287T mutation found in experiments (i.e., a noticeable increase in both the ON and OFF gating currents kinetics but only a small change in the charge–voltage relationship). The output of our Brownian model was quite similar to the experimental results, suggesting that the polarizability change caused by the I287T mutation is a major determinant for the effects observed on the macroscopic gating currents. More specifically, we found that the mutation caused a maximum decrease in the time constant of the ON gating currents of ∼57% of its WT value, not very far from the ∼48% found in experiments (Lacroix et al., 2013). Thus, the real effect of the I287T mutation on the gating current kinetics appears slightly higher than that predicted by our multiscale model. In this regard, it should be pointed out that mutagenesis experiments have shown that the size of the residue is also important (Lacroix et al., 2013, 2014). Considering the critical location of this residue (very close to the moving gating charges), it is not surprising that steric effects may also be involved, perhaps by changing the friction experienced by the moving voltage sensor and thus the rate of the gating currents.

We then looked at the mechanism linking the effective dielectric constant of the WA region to the rate of movement of the voltage sensor by analyzing the effects of the mutation on the energy profile encountered by the voltage sensor as it moves into the gating pore. We found that the reduction in the dielectric constant of the WA region produces a selective decrease in the first of the four energy peaks encountered by the voltage sensor, corresponding to the entry of the first gating charge (R1) into the gating pore. This result is not unexpected since this step involves the entry of a charge into the WA region, a process modulated by its polarizability. It is important to realize that what we show is a prediction of gating currents performed without using adjustable parameters, not a fitting routine where free parameters are adjusted to match the model output with the experimental data (Horng et al., 2019).

In this study, we have analyzed in detail the functional properties of Shaker channel gating, with isoleucine or threonine at location 287, since the type of residue present at this position (or equivalent position in Na_V channels) is thought to create the kinetic differences between Na_V and K_V channels (Lacroix et al., 2013). The results obtained suggest that different activation kinetics of voltage-gated Na and K channels come from different dielectric forces exerted by threonine compared with isoleucine, as the first gating charge enters the WA region of the gating pore. It seems that evolution uses a single side chain to control the speed of the K_V and Na_V channels and ultimately allow the action potential to exist.

Other processes could contribute to the time course of Na and K channels, as pointed out by Francisco Bezanilla. First, residue 363 in S4 is a highly conserved threonine in domains I, II, and III of Na_V channels, while it is a highly conserved isoleucine in K_V channels. Mutation of this isoleucine to threonine makes the voltage sensor of K_V channels move faster, as inferred from the time course of the gating currents (Lacroix et al., 2013). From the position of this residue relatively far from the gating charge R1, we suspect that its contribution to the polarizability would not be significant, suggesting that a different mechanism is likely responsible for this effect. Second, the presence of the β1 subunit of Na_V channels also accelerates the voltage sensor by an unknown mechanism (Lacroix et al., 2013). Finally, the sensors in Na_V channels exhibit positive cooperativity, which also accelerates the overall speed of channel activation (Chanda et al., 2004).

We point out that, within the limitations of our models, what we propose is both logically sufficient and logically necessary to account for the difference in speed and the structural movements of the voltage sensor. It is sufficient because we reproduce the properties of gating current under a realistic range of conditions. It is necessary because the movement of charges that we calculate must produce the currents we calculate given the universal and exact nature of the Maxwell equations that link charge and current. Of course, evolution is not logical. It often provides redundant mechanisms that are beyond the necessary and sufficient, that is to say, beyond the barely sufficient. These redundant mechanisms may provide properties not glimpsed in gating currents studied with step functions in the way they are usually measured. These mechanisms may make the gating system controllable and/or robust in conditions we have not explored. After all, most devices designed by engineers have multiple controls of the same output. It is not unreasonable to expect the same in devices designed by evolution, like the VSD.

A similar multiscale approach may also be used in the future to understand the physics of the functional effects of the many other mutations that have been performed on the VSD, thereby giving a deeper understanding of voltage-dependent gating.

Jeanne M. Nerbonne served as editor.

We are grateful to Prof. Bezanilla for teaching us the importance of his results in a discussion that motivated the models of this paper. We thank the referees for suggesting significant improvements to our calculations and presentation of results.

This work was supported by the University of Perugia (grant Ricerca di Base 2017).

The authors declare no competing financial interests.

Author contributions: L. Catacuzzeno conceived the modeling strategy and performed the theoretical simulations. L. Catacuzzeno, L. Sforna, F. Franciolini, and R.S. Eisenberg discussed the model results. L. Catacuzzeno, F. Franciolini, and R.S. Eisenberg wrote the paper. All authors approved the final version of the paper.