Continuum Gating Current Models Computed with Consistent Interactions.

The action potential of nerve and muscle is produced by voltage-sensitive channels that include a specialized device to sense voltage. The voltage sensor depends on the movement of charges in the changing electric field as suggested by Hodgkin and Huxley. Gating currents of the voltage sensor are now known to depend on the movements of positively charged arginines through the hydrophobic plug of a voltage sensor domain. Transient movements of these permanently charged arginines, caused by the change of transmembrane potential V, further drag the S4 segment and induce opening/closing of the ion conduction pore by moving the S4-S5 linker. This moving permanent charge induces capacitive current flow everywhere. Everything interacts with everything else in the voltage sensor and protein, and so it must also happen in its mathematical model. A Poisson-Nernst-Planck (PNP)-steric model of arginines and a mechanical model for the S4 segment are combined using energy variational methods in which all densities and movements of charge satisfy conservation laws, which are expressed as partial differential equations in space and time. The model computes gating current flowing in the baths produced by arginines moving in the voltage sensor. The model also captures the capacitive pile up of ions in the vestibules that link the bulk solution to the hydrophobic plug. Our model reproduces the signature properties of gating current: 1) equality of ON and OFF charge Q in integrals of gating current, 2) saturating voltage dependence in the Q(charge)-voltage curve, and 3) many (but not all) details of the shape of gating current as a function of voltage. Our results agree qualitatively with experiments and can be improved by adding more details of the structure and its correlated movements. The proposed continuum model is a promising tool to explore the dynamics and mechanism of the voltage sensor.

Introduction

Much of biology depends on the voltage across cell membranes. The voltage across the membrane must be sensed before it can be used by proteins. Permanent charges move in the strong electric fields within membranes, so carriers of sensing charge were proposed as voltage sensors even before membrane proteins were known to span lipid membranes (1). The movement of permanent charges of the voltage sensor is gating current, and the movement is the voltage-sensing mechanism. Permanent charge is our name for a charge or charge density independent of the local electric field (for example, the charge and charge distribution of Na+ but not the charge in a highly polarizable anion like Br− or the nonuniform charge distribution of H2O in the liquid state with its complex time dependent (and perhaps nonlinear) polarization response to the local electric field).

Knowledge of membrane protein structure has allowed us to identify and look at the atoms that make up the voltage sensor. Protein structures do not include the membrane potentials and macroscopic concentrations that power gating currents, and therefore, simulations are needed. Atomic-level simulations like molecular dynamics (MD) do not provide an easy extension from the atomic timescale ∼10−15 s to the biological timescale of gating currents that starts at ∼10−6 s and reaches ∼10−2 s. Calculations of gating currents from simulations must average the trajectories (lasting ∼10−1 s sampled every 10−15 s) of ∼106 atoms, all of which interact through the electric field to conserve charge and current while conserving mass. It is difficult to enforce continuity of current flow in simulations of atomic dynamics because simulations compute only local behavior, whereas continuity of current is global, involving current flow far from the atoms that control the local behavior. It is impossible to enforce continuity of current flow in calculations that assume equilibrium (zero net flow) under all conditions.

A hybrid approach is needed, starting with the essential knowledge of structure but computing only those parts of the structure used by biology to sense voltage. In close-packed (“condensed”) systems like the voltage sensor or ionic solutions, “everything interacts with everything else” because electric fields are long ranged as well as exceedingly strong (2). In ionic solutions, ion channels, even enzyme active sites, steric interactions that prevent the overfilling of space in well-defined protein structures are also of great importance because they produce short-range correlations (3).

Closely packed charged systems are well handled mathematically by energy variational methods. Energy variational methods guarantee that all variables satisfy all equations (and boundary conditions) at all times and under all conditions and are thus always consistent. We use the energy variational approach developed in (4) and (5) to derive a consistent model of gating charge movement, based on the basic features of the structure of crystallized voltage-sensitive channels. A schematic of the model is shown below. The continuum model we use simulates the mechanical dynamics in a single voltage sensor, although the experimental data is from many independent voltage sensors. Ensemble averages of recordings of individual independent voltage sensors are equivalent to macroscopic continuum modeling in a single voltage sensor if correlations are captured correctly in the model of the single voltage sensor.

Materials and Methods

Theory: Mathematical model

The reduced mechanical model for a voltage sensor is shown in Fig. 1 a with four arginines (R, i = 1, 2, 3, 4), each attached to the S4 helix by identical springs with the same spring constant K. The electric field will drag these four arginines because each arginine carries +1 charge. The charged arginines can also move as a group. S4 connects to S3 and S5 at its two ends by identical springs with spring constant KS4/2.

Figure 1

(a) Geometric configuration of gating pore in this model, including the attachments of arginines to the S4 segment. (b) Following (a), an axisymmetric three-zone domain shape is designated in r-z coordinate for the current 1D model. Here, the diameter of the hydrophobic plug is 0.3 nm (arginine’s diameter); L = 0.7 nm; L = 1.5 nm; and the radius of the vestibule is R = 1 nm. BC means boundary condition. To see this figure in color, go online.

Once the membrane is depolarized from, for example, −90 mV inside negative to +10 mV inside positive, arginines together with S4 will be driven toward the extracellular side. A repolarization from +10 to −90 mV moves the arginines back to the intracellular side. This movement is the basic voltage-sensing mechanism. The movement of S4 triggers the opening or closing of the lower gate—consisting mainly of S6 forming the ion permeation channel—by a mechanism widely assumed to be mechanical, although electrical aspects of the linker motion are likely to be involved as well.

When arginines are driven by an electric field, they are forced to move through a hydrophobic plug composed of several nonpolar amino acids from S1, S2, to S3 (6). Arginines reside initially in the hydrated lumen of the intracellular vestibule. They then move though the hydrophobic plug and wind up in the vestibule on the extracellular side. This movement involves dehydration when the arginines move through the hydrophobic plug, in which the arginines encounter a barrier in the potential of mean force (PMF), mainly dominated by the difference of the solvation energy in bulk situation and in the hydrophobic plug (7). Note that Na+ and Cl− (which are the only ions in the bulk solution in this article for simplicity) are found only in vestibules and are not allowed into the hydrophobic plug in our model. The ends of the two vestibules on each side of the hydrophobic plug act as impermeable walls for Na+ and Cl− in our model. When the voltage is turned on and off, these two walls store/release charge (carried by ions) in their electric double layers (EDL) that have many of the properties of capacitors.

In this continuum model, the four arginines (R, i = 1, 2, 3, 4) are described by their individual density distributions (concentrations) (c, i = 1, 2, 3, 4), allowing the arginines to interact with Na+ and Cl− in vestibules. The density (i.e., concentration) distributions represent probability density functions as shown explicitly in the theory of stochastic processes used to derive such equations in (8) using the general methods of (9). The important issue here is how well the correlations are captured in the continuum model. Some are more likely to be faithfully captured in molecular or coarse-grained dynamics simulations (e.g., more or less local hard sphere interactions) (10, 11, 12, 13, 14) and others in continuum models (e.g., correlations induced by far-field boundary conditions like the potentials imposed by bath electrodes to maintain a voltage clamp) (15, 16, 17, 18).

Here, we treat the S4 itself as a rigid body, so we can capture the basic mechanism of a voltage sensor without considering the full details of structure, which might lead to a three-dimensional model difficult to compute in reasonable time. We construct an axisymmetric one-dimensional (1D) model with a three-zone geometric configuration illustrated in Fig. 1 b, following Fig. 1 a. Zone 1 with z ∈ [0, L] is the intracellular vestibule; zone 2 with z ∈ [L, L + L] is the hydrophobic plug; zone 3 with z ∈ [L + L, 2L + L] is the extracellular vestibule. Arginines, Na+, and Cl− can all reside in zone 1 and 3. Zone 2 only allows the residence of arginines, albeit with a severe hydrophobic penalty because of their permanent charge, in a region of low dielectric coefficient, hence called hydrophobic.

Based on Fig. 1 b, the governing 1D dimensionless Poisson-Nernst-Planck (PNP)-steric equations are expressed below with the detailed nondimensionalization process shown in Supporting Materials and Methods, Section S1. The first one is a Poisson equation that shows how charge creates potential:where ϕ is electric potential; c is concentration of species i with valence q = 1, q = −1, q = q = 1, i = 1, 2, 3, 4; with being the Debye length, and the characteristic length (radius of vestibule) R = 1 nm here. A(z) is the channel cross-sectional area at position z. For zones 1 and 3, Γ = 1 by setting NaCl bulk concentration c = 184 mM and ε = 80. For zone 2, we assume a hydrophobic environment with ε = 8 and therefore Γ = 0.1. The value of the dielectric constant inside the hydrophobic plug (zone 2) is not experimentally available; however, the computational result is not sensitive to this value based on our sensitivity analysis.

The second equation is the species transport equation based on conservation laws:with the content of flux J expressed below based on the Nernst-Planck equation for Na+ and Cl−:and for four arginines c, i = 1, 2, 3 and 4 based on the Nernst-Planck equation with steric effect and some imposed potentials:where D is the diffusion coefficient for species i.

The first and second terms in Eqs. 3 and 4 describe diffusion and electromigration, respectively. The third terms in Eq. 4 are external potential terms with V, i = 1, 2, 3, and 4 being the constraint potential for the four arginines c to S4, represented here by a spring connecting each arginine c to S4, as shown in Fig. 1 a. Governing equations Eqs. 1, 2, 3, and 4 were derived by energy variational methods, which is further shown in Supporting Materials and Methods, Section S3.

The elastic system is described bywhere K is the spring constant, z is the fixed anchoring position of the spring for each arginine c on S4, and Z(t) is the center-of-mass z position of S4 by treating S4 as a rigid body. Here, we set z1 = 0.6, z2 = 0.2, z3 = −0.2, and z4 = −0.6 using structural information that gives the arginine anchoring interval on S4 as 0.4 nm. Z(t) follows the motion of equation based on the spring-mass system:where m, b, and K are the mass, damping coefficient, and restraining spring constants for S4. Z is the resting position of Z(t). Here, z is the center of mass for the set of arginines c, which can be calculated by

We assume that the spring-mass system for S4 is overdamped, which means the inertia term in Eq. 6 can be neglected.

The energy barrier V in Eq. 4 is nonzero only in zone 2, which mainly represents the difference in solvation energy, chiefly characterized by the difference of dielectric constants, in the hydrophobic plug and bulk solution. The structure of the energy barrier is actually very complicated. Here, we simply assume a hump shape for PMF (see more in Supporting Materials and Methods, Section S2), although we will seek greater realism in later work.

The last term in Eq. 4 is the steric term that accounts for steric interaction among arginines (5, 19). Here, we set g = 0.5, a reasonable value. Though there is actually no experimental measurement available for g, the computation results have been verified to be insensitive to its value.

Here, we assume quasisteady state for Na+ and Cl−, which means in Eq. 2, and the reasons are elaborated in Supporting Materials and Methods, Section S4. The formulation of boundary and interface conditions is also shown in Supporting Materials and Methods, Section S5.

Besides the main input parameter V, which is the applied voltage bias (corresponding to the command potential in voltage-clamp experiments), other parameters like D (i = 1, 2, 3, 4), K, K, and b are also required. Results are especially sensitive to the values of K, K, and b. We have tried and found D = 50; i = 1, 2, 3, and 4; K = 3; K = 3; and b = 1.5 provide the best fit to the experimental Q(charge)-voltage (QV) curve reported in (20). Some additional explanation on fitting these parameter values is described in Supporting Materials and Methods, Section S6.

Usually, the electric current in the ion channel is treated simply as the flux of charge and is uniform in the z direction when steady in time. This is not so in this nonsteady dynamic situation because the storing and releasing of charge in vestibules is involved. Here, the flux of charge at the middle of hydrophobic plug, z = L + L/2, was computed to estimate the experimentally observed gating current. However, it is actually impossible (so far) to experimentally measure the current at the middle of the hydrophobic plug. In experiments, the voltage-clamp technique is used, and on/off gating current through the membrane is measured, which should be equal to the flux of charge at z = 0 in this framework, as shown in Fig. 1 b. The flux of charges at any z position I(z, t) can be related to the flux of charges at z = 0, I(0, t), simply by charge conservation:whereand flux of charges at any z position I(z, t) is defined by

We identify as the displacement current and denote it as I (z, t) because Eq. 8 is equivalent to Ampere’s law in Maxwell’s equations, and is exactly the displacement current in Ampere’s law. The proof is elaborated on in Supporting Materials and Methods, Section S7. A general discussion about displacement current can be found in (21, 22, 23), which does not involve assumptions concerning the dielectric coefficient ε or polarization properties of matter at all. Hence, Eq. 8 can be simply rewritten aswhere we define the sum of displacement current and flux of charges as the total current I (z, t). The z distribution of the total current should be uniform by Kirchhoff’s law, and we verify this by computations shown in the section under heading Flux of Charges at Different Locations. Note the ionic current I(z, t) changes a great deal with location. The displacement current I(z, t) varies a great deal with location. The total current, the sum I(z, t), does not vary at all with location, although of course it varies a great deal with time. For example, calculations of current in the baths (which are not reported here) would show only ionic current in the time range considered here, but it would equal the total current that flows anywhere in our 1D model of the voltage sensor domain.

We are also interested in observing the net charge at vestibules. Consider, for example, the net charge at the intracellular vestibule, Q(L, t). The net charge consists of arginine charges and their countercharges formed by the EDL of ionic solution in that location. Electroneutrality is approximate but will not be exact there. Flux of charge, displacement current, and net charge at vestibules will be discussed further in the section under heading Flux of Charges at Different Locations.

To evaluate the current theoretical model, it is important to compare our computational results with experimental measurements (20) in the curves of gating current and amount of gating charge moved versus applied voltage (I(current)-voltage [IV] and QV curves). To construct the QV curve, we calculate , , which are the amounts of arginine found in zone 1, 2, and 3, respectively. Usually Q2 ≈ 0 is due to the energy barrier V in zone 2. Arginines tend to jump across zone 2 when driven from zone 1 to zone 3 as the voltage V is turned on. The number of arginines that move and settle at zone 3 depends on the magnitude of V. Besides IV and QV curves, the time course of the movement of arginines and S4, z(t) and Z(t), is important to report here because recording these movements in experiments is becoming feasible nowadays by optical methods. Many qualitative models accounting for the movement of S4 and conformation change of the voltage sensor have been proposed. Readers are referred to review articles (24, 25) for more details.

Numerical method

Equations 1, 2, 3, and 4 are first discretized in space by high-order multiblock Chebyshev pseudospectral methods and then integrated in time under the framework of method of lines. The details of the numerical method are referred to Supporting Materials and Methods, Section S8.

Results and Discussion

Here, numerical results based on the mathematical model described above were calculated and compared with experimental measurements (20). Our 1D continuum model has advantages and disadvantages. The lack of three-dimensional structural detail means that some details of the gating current and charge cannot be reproduced. It should be noted, however, that to reproduce those, one needs more than just static structural detail. One must also know how the structures (particularly their permanent and polarization charge) move and change after a command potential is applied in the experimental ionic conditions. The 1D model has advantages because it computes the actual experimental results on the actual experimental timescale in realistic ionic solutions and with far-field boundary conditions actually used in voltage-clamp experiments. It also conserves total current, as we will demonstrate later. Conservation of current needs to be there and verified in theories and simulations because it is a universal property of the Maxwell equations (21, 22, 23).

QV curve

When the membrane and voltage sensor is held at a large inside negative potential (e.g., hyperpolarized to −90 mV), S4 is in a resting potential position, and all arginines stay in the intracellular vestibule. When the potential is made more positive (e.g., depolarized to +10 mV), S4 is in the active potential position, and all arginines are at the extracellular vestibule.

The voltage dependence of the charge (arginines) transferred from intracellular vestibule to extracellular vestibule is characterized as a QV curve in experimental papers, and it is sigmoidal in shape (20). Fig. 2 a shows that our computed QV curve—the dependence of Q3 on V—is in very good agreement with the experiment (20). This good agreement comes from the fact that our resultant QV curve is also a sigmoidal curve, and, most important of all, the slope of QV curve can be tuned, mainly by the adjustment of K, K, and b, to agree with experiment. Not many theoretical models can achieve this agreement, especially for the slope. Models in (15, 16) show good agreement with experiments, whereas a mismatch of slope was reported in (17, 18). The voltage dependence of activation has been considered a crucial property of the sodium conductance since it was defined (1). Fig. 2 b shows the steady-state distributions of Na+, Cl−, and arginines in the inside negative, hyperpolarized situation (V = −90 mV). As we can see, all the arginines stay in the intracellular vestibule, and none of the arginines move to the extracellular vestibule (Q3 ≈ 0).

Figure 2

(a) QV curve and comparison with (20). Steady-state distributions for Na+, Cl−, and arginines are shown at (b) V = −90 mV, (c) V = −48 mV, and (d) V = −8 mV. Note that the experimental data in (20) were scaled to 4e. To see this figure in color, go online.

Fig. 2 c shows the situation at V = −48 mV, which is the midpoint of the QV curve. As we can see, each vestibule has distributions of c (i = 1, 2, 3, and 4), resulting in half of the arginines staying in it (Q3 = 2). The center-of-mass position for each arginine, presented later in Fig. 6, shows that R1 and R2 are in the extracellular vestibule, and R3 and R4 are in the intracellular vestibule. There are almost no arginines in zone 2 (hydrophobic plug) because of the energy barrier in it. Note that this represents an average because in a single molecule interpretation, half of the sensors will be with all R’s inside and the other half with all R’s outside. The midpoint of −48 mV from (20) requires the resting position of S4, Z, to be biased from L + 0.5L to Z = L + 0.5L + 1.591 nm; otherwise, the midpoint would be 0 mV. Fig. 2 d shows the situation at full depolarization (V = −8 mV), at which time all arginines move to the extracellular vestibule (Q3 ≈ 4) in the fully depolarized, activated state.

Figure 6

(a) and (c) are the time courses of the amount of arginines moved to the extracellular vestibule. (b) and (d) are center-of-mass trajectories of individual arginines and S4. (a) and (b) are the case of large depolarization with V rising from −90 to −8 mV at t = 10, holding on till t = 150, and then dropping back to −90 mV. (c) and (d) are the case of small depolarization with V rising from −90 to −50 mV at t = 10, holding on till t = 150, and then dropping back to −90 mV. To see this figure in color, go online.

Gating current

Fig. 3 shows the time course of gating currents, observed as flux of charge at the middle of hydrophobic plug I(L + L/2, t) because of the movement of arginines when the membrane depolarization is large and when the depolarization is small. In the case of large depolarization, V rises from −90 mV at t = 10 to −8 mV and drops back to −90 mV at t = 150 (Fig. 3 a). The time course of gating current and contributions of individual arginines are shown in Fig. 3 b. As expected, the rising order of each current component follows the moving order of R1, R2, R3, and R4 when depolarized and that order is reversed when repolarized. The area under the gating current is the amount of charge moved. Because arginines move forward and backward in this depolarization/repolarization scenario, the areas under the ON current and the OFF current are same. The areas are equal for each component of current as well. The equality of area is an important signature of gating current that contrasts markedly with the properties of ionic current (26, 27). In the case of small depolarization (V rises from −90 to −50 mV at t = 10 to and drops back to −90 mV at t = 150, Fig. 3 c), the time course of gating current and its four components contributed by each arginine for this situation is shown in Fig. 3 d. Under this small depolarization, not all arginines move past the middle of the hydrophobic plug because of the weaker driving force in the small depolarization compared with the large depolarization case. This can be inferred because the areas under each component current are different (Fig. 3 d).

Figure 3

(a) The time course of V rising from −90 to −8 mV at t = 10 holds on until t = 150 and then drops back to −90 mV. (b) The time course of gating current, I(L + L/2, t), and its components corresponding to (a) are shown. (c) The time course of V rising from −90 to −50 mV at t = 10 holds on until t = 150 and then drops back to −90 mV. (d) The time course of gating current, I(L + L/2, t), and its components corresponding to (c) are shown. To see this figure in color, go online.

The gating currents can be better understood by looking at a sequence of snapshots showing the spatial distribution of electric potential, species concentration, and electric current. The distributions at several times are shown in Fig. 4 a for the case of sudden change in command voltage to a more positive value and a large depolarization, and the distributions are shown in Fig. 4 b for the case of a small depolarization. The electric potential profiles at t = 13 and t = 148 show that the profile of electric potential changes as arginines move from left to right even though the voltage is maintained constant across the sensor. Slight bulges in electric potential profile exist wherever arginines are dense. This can be easily explained by understanding the effect of Eq. 1 on a concave spatial distribution of electric potential.

Figure 4

(a) The top row shows dimensionless species concentration distributions at t = 0, 13 (right after depolarization), and 148 (right before repolarization) for the case of large depolarization with V from −90 mV at t = 10 to −8 mV and dropping back to −90 mV at t = 150. The middle row shows concurrent electric potential profiles. The bottom row shows concurrent electric current profiles with the components of flux of charge, displacement current, and total current. (b) The same as (a) is shown except with V depolarized from −90 to −50 mV. To see this figure in color, go online.

In Fig. 4, the total current defined in Eq. 11, though changing with time, is always constant in z at all times, satisfying Kirchhoff’s law (i.e., conservation of current). At t = 13, when gating current is substantial, as seen from t = 13 in Fig. 3, b and d, we can visualize the z distributions of flux of charges I(z, t), displacement of current I(z, t), and total current I(z, t) individually in Fig. 4.

Flux of charges at different locations

Flux of charges I(z, t), together with displacement current I(z, t) and total current I(z, t), depicted in Fig. 4, deserve more discussion here. Though I(z, t), I(z, t), and I(z, t) are well defined in Eqs. 8, 9, 10, and 11, the actual computation of them takes an indirect path because of the assumption of quasisteady state for Na+ and Cl− in Eq. 2. The details are presented in Supporting Materials and Methods, Section S9. The computed total current I(z, t) does indeed satisfy Kirchhoff’s law by its uniformity in z. This verification is shown in Fig. 4 at several times, and we have checked that this is in fact true at any time.

In the bottom rows of Fig. 4 at t = 13, we observe that I(z, t) is generally nonuniform in z and is accompanied by congestion/decongestion of arginines in between. However, I(z, t) is almost uniform in zone 2 (hydrophobic plug), which means almost no congestion/decongestion of arginines occurs there, and therefore, there is no contribution to the displacement current from zone 2. This is because arginines can hardly reside in zone 2 because of the energy barrier in it.

Several things are worth noting in the time courses of and I(0, t) (equal to uniformly distributed I as depicted by Eq. 11) illustrated in Fig. 5 a under the case of large depolarization. First, is noticeably larger than I(0, t) in the ON period. This is because their difference, exactly the displacement current I, is always negative at zone 2 when depolarized because arginines are leaving zone 1 and make  < 0 for zone 2. It is expected that the area under the time course of would be very close to 4e, as verified by the time courses of Q3 in Fig. 5 b. We use I(0, t) to estimate the experimentally measured voltage-clamp current, whereas the counterpart area of experimentally measurable I(0, t) would be less than 4e because of its smaller magnitude compared with . This may partly explain the experimental observations that at most 13e (25, 28, 29), instead of 16e, are moved during full depolarization in four voltage sensors (for a single ion channel) based on computing the area under voltage-clamp gating current. Therefore, flux of charge at any location of zone 2, though impossible to measure in experiments so far, will give us the amount of arginines moved during depolarization more reliably than the measurable I(0, t).

Figure 5

(a) The time courses of , I(0, t), and despiked I(0, t) for the case of large depolarization with V rising from −90 to −8 mV at t = 10, holding on till t = 150, and then dropping back to −90 mV. The inset plot is a magnification of the ON current to visualize the difference of I(0, t) and despiked I(0, t) more clearly. (b) The time courses of , and are under the same depolarization scenario as (a). To see this figure in color, go online.

Second, we see in Fig. 5 a with magnification in its inset plot that, as in experiments, I(0, t), but not , has contaminating leading spikes in ON and OFF parts of the current. These spikes are capacitive currents from solution EDL of vestibules caused by the sudden rising and dropping of command potential. These spikes need to be removed in voltage-clamp experiments to get rid of the contribution from vestibule solution EDL (and membrane) to the transport of gating charges (arginines) when computing the area under gating current. The technical details of removing these spikes are shown in Supporting Materials and Methods, Section S10, and more details about spikes can be found in Supporting Materials and Methods, Section S11.

Third, in Fig. 5 b, as arginines move from one vestibule to another, the concentrations of Na+ and Cl− also correspondingly change with time at the vestibules. They form countercharges through EDL and balance arginine charges at vestibules. However, these EDL changes only maintain an approximate, not exact, charge balance, as shown in Fig. 5 b. The violation of electroneutrality causes the displacement current, which is not negligible. This further causes the underestimate of arginines that move when the voltage sensor is depolarized if the estimate is made by measuring the area under I(0, t).

As in the previous section, we used flux of charges at the middle of the hydrophobic plug, I(L + L/2, t), instead of experimentally measurable I(0, t) to represent the gating current in discussions. We may as well name I(L + L/2, t) as the arginine current to avoid the confusion with the actual gating current I(0, t) here. This arginine current leaves out its associated displacement current I(L + L/2, t) and serves to represent gating current better for two reasons:

The area under the time course of I(L + L/2, t) gives us the amount of arginines moved during depolarization more faithfully than I(0, t). The fluxes of charge for each arginine shown in Fig. 3, b and d carry important information about how each arginine is moved by the electric field that will be further illustrated in Fig. 6. All these will not be easy to display and comprehend if we use I(0, t) instead.

Using I(0, t) as a definition of gating current would require a decontamination by removing the leading spikes, which is computationally costly. Removing spikes would especially pose a heavy numerical burden when doing parameter fitting in which numerous repeated computations are done.

Time course of arginine and S4 translocation

Fig. 6 shows the time course of Q (amount of arginines moved to extracellular vestibule, equal to Q3 here) and center-of-mass trajectories of individual arginines (z, i = 1, 2, 3, and 4) and S4 segment (Z). Fig. 6, a and b show the case of large depolarization, and Fig. 6, c and d show the case of small depolarization.

In the case of large depolarization (Fig. 6 b), the arginines and S4 z positions quickly reach individual steady states, with almost all arginines transferred to the extracellular vestibule as previously shown in Fig. 4 a. Therefore, Q is close to its saturated value 4 as shown in Fig. 6 a. Arginines and S4 move back to the intracellular vestibule once the voltage drops back to −90 mV. From Fig. 6 b, the forward-moving order of arginines is R1, R2, R3, and R4, and the backward-moving order is the opposite R4, R3, R2, and R1 with agreement with the structure. This agreement might look trivial in molecular dynamics simulations but is not a trivial check here because this model describes arginines not by particles, as in molecular dynamics, but by concentrations. Note that an incorrect order and pace of the movement of arginines would cause disagreement with experiments in the shape of IV curve as well. S4 is initially farthest to the right but lags behind R1 and R2 during movement in depolarization, as shown in Fig. 6 b. This is certainly because S4 is finally relaxed to an almost unforced situation close to its resting position Z during this large depolarization. We can further calculate the displacements of each arginine and S4 during this full-saturating depolarization and find Δz1, ≈ Δz2, ≈ Δz3, ≈ 1.93 nm, Δz4, = 1.76 nm, and ΔZ = 1.51 nm. Besides almost the same displacements for R1, R2, and R3, their average moving velocities are also very close to each other. This seems to suggest a synchronized movement among R1, R2, and R3 that we have not imposed on the arginines in our model. Also, we can see the movements of arginines contribute significantly to the movement of the S4 segment. This can be seen from the steady-state z position of S4 derived from Eq. 6,

Experimental estimates of S4 displacement during full depolarization range from 2 to 20 Å (24, 30), depending on the model of the voltage sensor and its motion, including the transporter model, the helical screw, and the paddle model (24). Our ΔZ = 1.51 nm here is large and seems to agree better with experimental estimates requiring large displacements, such as the paddle model. In contrast, the helical screw model, which is supported by most of the recent data, is known to have shorter displacements. A plausible explanation for our overestimate of ΔZ is that our 1D model uses a straight line perpendicular to the hydrophobic-plug path for the movement of the arginines. In reality, the S4 segment is significantly tilted with respect to the membrane, and the arginines follow a spiral along the helix. Therefore, if the S4 segment rotates and changes its tilt during activation, the total vertical translation needed to cross the hydrophobic plug is significantly reduced, as was shown by Vargas et al. (31). The value obtained in (31) was between 0.7 and 1 nm when comparing the displacement perpendicular to the membrane of the open-relaxed state crystal structure of Kv1.2 (32) and the closed structure that has been derived by consensus from experimental measurements (31).

In the case of small depolarization, the driving force is weaker than in a large-saturating depolarization, so their z positions do not have a chance to reach steady states as they do during a full-saturating depolarization. Rather, in a small depolarization, the motion of the arginines and S4 are aborted. They return to the intracellular vestibule because the depolarization drops (i.e., decreases in magnitude, and the membrane potential becomes more negative) before arginines and S4 have a chance to reach their steady-state positions. This detailed atomic interpretation likely overreaches the resolution of our model. At the single-sensor level, we do not expect partial movements; instead, some sensors will have moved all the way and others not at all, but the distribution of sensors in the two extreme positions should follow what we predict with this model, which is an ensemble average. We look forward to measurements of movements of probes that mimic arginine in its environment that require improvements in the resolution and structural realism of our model.

Fig. 6 c illustrates these aborted motions. Q reaches 1.57 at most, which should be 2 instead if steady-state was reached as it is if time is long enough. See the steady-state behavior shown in the QV curve of Fig. 2 a. Fig. 6 d shows that the S4 segment is initially farthest to the right, lags behind R1 during movement, and is almost caught up by R2. The maximal displacements of arginines and S4 calculated from Fig. 6 d are Δz1, = 1.36 nm, Δz2, = 0.966 nm, Δz3, = 0.459 nm, Δz4, = 0.316 nm, and ΔZ4, = 0.616 nm. The significant difference between Δz1,, Δz2,, Δz3,, and Δz4, may imply that R1 and R2 have jumped across the hydrophobic plug and entered the extracellular vestibule, whereas R3 and R4 still stay at the intracellular vestibule during this small depolarization. This is consistent with the observation from individual gating-current components of arginines in Fig. 3 d.

Family of gating currents for a range of voltages

Though we prefer to I(0, t) for representing gating current as explained in the section under heading Flux of Charges at Different Locations, we here use the actual gating current, despiked I(0, t), to compare with experiment (20). Fig. 7 a shows the time courses of a subtracted gating current (despiked I(0, t)) for a range of voltages V ranging from −62 to −8 mV. The area under gating current, for both ON and OFF parts, increases with V because more arginines are transferred to the extracellular vestibule as V increases. The shapes of this family of gating currents agree well with experiment (20) in both magnitude and time course.

Figure 7

(a) The time courses of subtracted gating current, despiked I(0, t), with the voltage rising from −90 to V mV at t = 10, holds on till t = 150, and then drops back to −90 mV, where V = −62, −50, …, −8 mV. (b) τ2 versus V compared with experiment (20) is shown. To see this figure in color, go online.

We can characterize the time course by fitting the decay part of a subtracted gating current by , τ1 < τ2 as generally done in experiments (20) in which τ1 is the fast time constant and τ2 is the slow time constant. Usually, the movement of arginines is dominated by τ2. Here, τ2 was calculated from simulation and compared with experiment (20) as shown in Fig. 7 b. Because in our computation the time is in arbitrary units, we have scaled the time to have the maximal τ2 to fit with its counterpart in experiment (20). Overall, the trend of τ2 versus V in our result, though not the whole curve, agrees well with experiment (20). To the left of the maximal point in Fig. 7 b, simulation results fit rather well with the experiment compared with the values to the right of the maximal point, at which it overestimates τ2 compared with the experiment. This overestimate is consistent with the observation that the amount of transferred charges Q saturates slightly faster in experimental data than in this simulation as V increases (see QV curve of Fig. 2 a). This phenomenon is related to the cooperativity of movement among arginines, which will be further discussed below.

Effect of voltage pulse duration

Fig. 8 shows the effect of voltage pulse duration with Fig. 8 a for the case of small depolarization and Fig. 8 b for the case of large depolarization. The magnitude and time span of subtracted gating current (despiked I(0, t)) are changed by pulse duration in both cases, but the shape will asymptotically approach the same curve as pulse duration increases, no matter the size of the depolarization. This behavior occurs because it takes time for the command pulse to drive the arginines toward the extracellular vestibule. If the pulse duration is long enough, the time course of Q will approach its steady state for large depolarization as in Fig. 6 a. Small depolarization takes a longer time to reach its steady state, as demonstrated in Fig. 6 c. The shapes of gating currents in Fig. 8 compare favorably with experiment (20) in which the OFF subtracted gating currents for short pulses have very fast decays, whereas for long pulses, the OFF subtracted gating currents have larger rising amplitude and slower decay because of a larger amount of arginines moved.

Figure 8

Subtracted gating currents, despiked I(0, t), showing the effect of voltage pulse duration. (a) V increases from −90 to −35 mV at t = 10 and drops back to −90 mV at various times. (b) V increases from −90 to 0 mV at t = 10 and drops back to −90 mV at various times. To see this figure in color, go online.

Conclusions

Previous work with molecular and coarse-grained simulations have captured some interactions, but they have not yet reproduced the time course and voltage dependence of macroscopic gating currents (10, 11, 12, 13, 14), and previous continuum models have captured only the steady-state properties of charge movement (15, 16, 17, 18).

This 1D continuum mechanical model of the voltage sensor tries to capture the essential structural details of the movement of mass and charge that are necessary to reproduce the basic features of experimentally recorded gating currents. After finding appropriate parameters, we find that the general kinetic and steady-state properties are well represented by the simulations. The good agreement of our numerical results with salient features of gating current measured experimentally would be impossible by simply tuning of parameters if our model had not captured the essence of physics for the voltage sensor. The continuum approach seems to be a good model of voltage sensors, provided that it 1) takes into account all interactions crucial to the movement of gating charges and S4; 2) computes their correlations consistently, so all variables satisfy all equations under all conditions with one set of parameters; and 3) satisfies conservation of current. This last point gave us a new insight: what is measured experimentally does not correspond to the transfer of the arginines because the total current, containing a displacement current, is smaller than the arginine current. It should be noted, however, that the total energy provided by the voltage clamp is qV, where q is the time integral of the measured gating current and V is the applied voltage. This is the total energy that explains the correspondence of charge per channel with the charge estimated by the limiting slope method (33, 34, 35).

We have simplified the profile of the energy barrier in the hydrophobic plug because the PMF in that region, and its variation with potential and conditions, is unknown. There is plenty of detailed information on the amino acid side chains in the plug and how each one of them changes the kinetics and steady-state properties of gating charge movement (6). Therefore, the next step is to model the details of interactions of the moving arginines with the wall of the hydrophobic plug and the contributions from other surrounding charged protein components. Some of the effects to be included are as follows:

Steric and dielectric interactions of the arginines that this model does not include. These include the interaction of arginines with negative charges of the S2 and S3 segments and the negative phospholipids as well as the hydrophobic residues in the plug. These interactions may be responsible for the simultaneous movement of two to three arginines across the plug, which is an experimental result that this model does not reproduce (36, 37).

Time dependence of the plug energy barrier V. Once the first arginine enters the hydrophobic plug by carrying some water with it, this partial wetting of the hydrophobic plug will lower V, chiefly consisting of solvation energy, and enable the next arginine to enter the plug with less difficulty. This might explain the cooperativity of movement among arginines when they jump through the plug. The addition of details in the plug may also produce intermediate states that have been measured experimentally. In this situation, arginines may transiently dwell within the plug.

A very strong electric field might affect the hydration equilibrium of the hydrophobic plug and would lower its hydration energy barrier as well (38). This cooperativity of movement may help explain the quick saturation in the upper right branch of the QV curve (and smaller τ2). It may also explain the experimentally observed translocation of two to three arginines simultaneously (36, 37).

The power of this mathematical modeling is precisely the implementation of interactions and the various effects in a consistent manner. Implementing the various effects listed above is likely to lead to a better prediction of the currents and to the design of experiments to further test and extend the model.

Further work must address the mechanism of coupling between the voltage sensor movements and the conduction pore. For example, the spring constant of the two sides of S4 have been made equal, which does not take into account the structural reality that one side has a linker to S3, whereas the other links to the pore opening. It seems likely that the classical mechanical models of coupling will need to be extended to include coupling through the electrical field. The charges involved are large. The distances are small, so the changes in electric forces that accompany movements of charged mass (and flows of displacement current) are likely to be large and important. It is possible that the voltage sensor modifies the stability of a fundamentally stochastically unstable, nearly bistable, conduction current (of single channels) by triggering sudden transitions from closed to open state in a controlled process reminiscent of Coulomb blockade in a noisy environment (39).

Author Contributions

All authors conceived the research, T-L. H. wrote the code and carried out the computations, and all authors contributed to the interpretation and writing of the manuscript.