Determining the molecular basis of voltage sensitivity in membrane proteins

The membrane potential (MP) in living cells results from the uneven distribution of ions between the two sides of the cell membrane (Hille, 2001). It regulates several critical physiological processes, such as the formation and propagation of the action potential in excitable cells (Hille, 2001) and the progression along the cell cycle in nonexcitable ones (Cone, 1969, 1971; Yang and Brackenbury, 2013). Voltage-sensitive membrane proteins are able to detect changes in the MP and in some cases respond directly to it (Bezanilla, 2008). The ability to understand and eventually modulate the function of these proteins provides an interesting strategy to diagnose and treat neurological diseases by changing either the MP or the proteins’ response to it.

The group of voltage-sensitive membrane proteins (VSMPs) includes representatives from transporters (Nakao and Gadsby, 1986; Catterall, 1988; Weer et al., 1988; Colombini, 1989; Schnetkamp and Szerencsei, 1991; Bernardi, 1992; Parent et al., 1992; Kavanaugh, 1993; Petronilli et al., 1994; Jentsch et al., 1995; Wadiche et al., 1995; Caterina et al., 1997; Halestrap et al., 1997; Scorrano et al., 1997; Kaim and Dimroth, 1998, 1999; Künkele et al., 1998; Lostao et al., 2000; Sugawara et al., 2000; Yao et al., 2000; Zheng et al., 2000; Jutabha et al., 2003, 2011; Mackenzie et al., 2003; Melzer et al., 2003; Bezanilla, 2005, 2008; Anzai et al., 2008; Hub et al., 2010; Zhang et al., 2011; Malhotra et al., 2013; van der Laan et al., 2013; Zander et al., 2013), enzymes (Reddy et al., 1995; Beltrán et al., 1996; Cooper et al., 1998; DeCoursey et al., 2003; Valdez and Boveris, 2007; Rosasco et al., 2015; Vorburger et al., 2016), receptors (Ben-Chaim et al., 2003; Rinne et al., 2013; Vickery et al., 2016b), and proteins playing a role in structure or adhesion (Bargiello et al., 2012; Fig. 1). These proteins are extremely diverse in their structure, but all of them have a common ability to convert the electrical energy into a mechanical response. To achieve this, VSMPs have developed one or more voltage-sensing elements, which detect changes in the MP and alter their conformational state accordingly (Bezanilla, 2008; Swartz, 2008). In voltage-gated potassium and sodium channels, the protein itself plays the role of the voltage sensor: the S4 helix of the S1–S4 helical bundle carries several positively charged residues, which on application of an electric field, move in the direction of this field (Long et al., 2005; DeCaen et al., 2009; Khalili-Araghi et al., 2010; Tao et al., 2010; Wu et al., 2010; Delemotte et al., 2011; Payandeh et al., 2011; Henrion et al., 2012; Jensen et al., 2012; Vargas et al., 2012; Lacroix et al., 2014; Machtens et al., 2017). In most other VSMPs, the voltage-sensing elements are nowhere near as well characterized. Several studies suggest different origins of their voltage sensitivity, including specific protein residues and/or ions trapped inside protein cavities (Nakao and Gadsby, 1986; Weer et al., 1988; Mirzabekov and Ermishkin, 1989; Thomas et al., 1993; Verselis et al., 1994; Pusch et al., 1995; Reddy et al., 1995; Popp et al., 1996; Boukalova et al., 2010; Kwon et al., 2012; Pinto et al., 2016; Vickery et al., 2016a). Ions moving in an electric field have, for instance, been suggested to trigger conformational rearrangements in voltage-gated chloride channels (Pusch et al., 1995; Chen and Miller, 1996; Smith and Lippiat, 2010; Zifarelli et al., 2012; Grieschat and Alekov, 2014; De Jesús-Pérez et al., 2016), solute carriers (Parent et al., 1992), active transporters (Nakao and Gadsby, 1986, 1989; Weer et al., 1988; Holmgren et al., 2000), and G-protein–coupled receptors (Vickery et al., 2016a).

In this work, we characterize the voltage sensitivity of several selected membrane proteins using a newly developed computational approach. It evaluates the propensity of the molecular system to detect changes in the local electric field and to respond to them. The system’s response propensity is directly connected to the gating charge, a well-known characteristic of voltage sensitivity that corresponds to the amount of charge transferred during protein activation. Our approach also estimates the local electric field response, which reflects how sensitive the local electric field is to the application of an external electric field. We show that, in all tested voltage-sensitive proteins, the system’s response propensity is large compared with voltage-insensitive ones, which suggests this could be an efficient tool to probe voltage sensitivity in other membrane proteins. We also show how the local electric field response can be large in selected regions of voltage-insensitive proteins. However, the lack of charges in these regions results in their inability to respond to changes in the MP. Finally, we investigate six VSMPs for which the voltage sensors have not yet been characterized. For each of them, we depict putative voltage-sensitive residues and/or voltage-sensitive ions trapped in protein cavities.

Our approach only requires knowledge of a protein structure and modest computational resources, which makes it possible to apply before mutagenesis experiments in order to reduce the number of potential voltage-sensitive elements to be tested.

The Charmm-GUI server was used to prepare the systems for molecular dynamics simulations (Jo et al., 2008). Briefly, every protein of interest was embedded into a 1-palmytoyl-2-oleoylphosphatidylcholine bilayer and solvated with 150 mM of either KCl or NaCl solution. The CHARMM36 force field (Mackerell et al., 2004) was used to describe proteins and lipids with TIP3P as a water model (Jorgensen et al., 1983). In the case of NavMs, we also considered a system with a TIP4P water model to test whether the results depend on the water model used. Fig. 2 A shows that our method is indeed robust to changes in the water model. However, the choice of force field parameters, and in particular of partial charge definition, can in principle influence the results. Note that the Na⁺ ion bound inside the protein cavities of the M₂ receptor was not placed there initially; instead, it reached the binding site during the molecular dynamics simulations. For the detailed description of the systems’ composition and properties, see Table 1.

The molecular dynamics simulations were performed by using GROMACS 2016.1 (Abraham et al., 2015). Each system was equilibrated by following a multistep protocol. During the first 2.5 ns, the protein and lipid headgroups were restrained to their initial positions to allow for the rearrangements of lipid tails and solution. The lipid headgroups were subsequently released, and the simulations continued for 2.5 ns. In the next 35 ns, the restraints applied to the protein backbone and sidechains were gradually decreased from 400 kJ/mol/Å and 200 kJ/mol/Å, respectively, to 0 kJ/mol/Å. Finally, the protein was fully relaxed without restraints until the root mean square deviation with respect to the initial structure reached a plateau value (for the duration of this last step of the equilibration, see Table 1).

A Nosé-Hoover thermostat (Nosé, 1984) and Parrinello-Rahman barostat (Parrinello and Rahman, 1981) were used to keep the temperature (300°K) and pressure (1 atm) constant. A cutoff of 12 Å was applied for short-range electrostatics and Van der Waals interactions. For the latter, a switching function between 10 and 12 Å was applied to smoothly bring the forces to 0 at 12 Å. The particle mesh Ewald summation (Darden et al., 1993) was used for the long-range component of electrostatics. A 1-fs time step was used for the first two steps of the equilibration and 2 fs for the rest of the equilibration and the simulations under an electric field.

For each protein of interest, 10 conformational states were extracted from the equilibration trajectory with a stride of 20 ns. For every state, 8 short (2-ns) molecular dynamics simulations were performed with restraints applied to all heavy atoms of the protein and under an electric field; the following values of the field were used: −0.020, −0.015, −0.010, −0.005, 0.005, 0.010, 0.015, and 0.020 V/Å. Note that these values correspond to an MP up to ∼0.5–3.5 V (depending on the system’s size). Smaller electric fields can be used as well but, in such cases, the uncertainty in the estimated variables is larger and the results are more qualitative (Fig. 2 A). Five frames were extracted from each of the trajectories to compute an average local electrostatic potential map $φ$ by using the PMEpot plugin of VMD (Aksimentiev and Schulten, 2005):

Here, $ρ i (r)$ corresponds to a point charge approximated by a spherical Gaussian with an Ewald factor of 0.25. Eq. 1 was solved on a grid with a resolution of 1 Å. R, $R j$ (per-element local electric field response: $R j =〈 R(r i j ) 〉$⁠; the average is computed over all atoms of the element j), and $ζ j$ were calculated by using an in-house python code (see Results and Discussion for more details). Finally, the average and the standard deviations of $R j$ and $ζ j$ were computed based on all 10 conformational states and all subunits in the case of homomultimers.

The detection threshold was identified based on the calculations of $R j$ (Fig. 2 B). For every protein of interest, $R j$ was computed for the two dimensions orthogonal to that of the external electric field (i.e., in the definition of $R j$ the components of the local field orthogonal to the external field were used). The largest $R j$ value, 0.02, was further considered as the detection threshold.

Video 1 shows the local electric field response estimated for the voltage-gated potassium channel Shaker, projected onto a ribbon representation of the protein backbone. Video 2 shows the local electric field response estimated for the voltage-gated sodium channel NavMs, projected onto a ribbon representation of the protein backbone. Video 3 shows the local electric field response estimated for Cx26, projected onto a ribbon representation of the protein backbone. Video 4 shows the local electric field response estimated for TRPV1, projected onto a ribbon representation of the protein backbone. Video 5 shows the local electric field response estimated for VDAC1, projected onto a ribbon representation of the protein backbone. Video 6 shows the local electric field response estimated for ClC1, projected onto a ribbon representation of the protein backbone. Video 7 shows the local electric field response estimated for the muscarinic acetylcholine receptor M₂, projected onto a ribbon representation of the protein backbone. Video 8 shows the local electric field response estimated for the Na⁺/K⁺ ATPase, projected onto a ribbon representation of the protein backbone. Video 9 shows the local electric field response estimated for GLIC, projected onto a ribbon representation of the protein backbone. Video 10 shows the local electric field response estimated for the two-pore domain potassium channel TWIK-1, projected onto a ribbon representation of the protein backbone.

Voltage sensitivity occurs in many different membrane proteins with diverse structure and function (Nakao and Gadsby, 1986; Catterall, 1988; Weer et al., 1988; Colombini, 1989; Schnetkamp and Szerencsei, 1991; Bernardi, 1992; Parent et al., 1992; Kavanaugh, 1993; Petronilli et al., 1994; Jentsch et al., 1995; Reddy et al., 1995; Wadiche et al., 1995; Beltrán et al., 1996; Caterina et al., 1997; Halestrap et al., 1997; Scorrano et al., 1997; Cooper et al., 1998; Kaim and Dimroth, 1998, 1999; Künkele et al., 1998; Lostao et al., 2000; Sugawara et al., 2000; Yao et al., 2000; Zheng et al., 2000; Ben-Chaim et al., 2003; DeCoursey et al., 2003; Jutabha et al., 2003, 2011; Mackenzie et al., 2003; Melzer et al., 2003; Bezanilla, 2005, 2008; Valdez and Boveris, 2007; Anzai et al., 2008; Hub et al., 2010; Zhang et al., 2011; Bargiello et al., 2012; Malhotra et al., 2013; Rinne et al., 2013; van der Laan et al., 2013; Zander et al., 2013; Rosasco et al., 2015; Vorburger et al., 2016); the unifying factor is how they are all able to detect changes in the MP and convert these changes into mechanical work (Bezanilla, 2008). The gating charge is a characteristic variable of voltage sensitivity. It was initially proposed to report the amount of charge transferred upon gating in voltage-gated ion channels (Armstrong and Bezanilla, 1973; Aggarwal and MacKinnon, 1996; Noceti et al., 1996) and then generalized to characterize any voltage-sensitive transition between two states of a given membrane protein. Computationally, the gating charge Q can be estimated via the excess free energies $ΔG( V m )$ caused by the application of the MP $V m$ (Roux, 2008; Treptow et al., 2009) in Eq. 2:

As such, it measures how the MP perturbs the equilibrium population of the activated versus deactivated states of a protein, referred to as a and d, respectively. It can be further decomposed into individual contributions from different protein elements (residues and ions trapped in cavities; Roux, 2008; Treptow et al., 2009) as $Q=∑ j Q j$⁠, where

In Eq. 3, the summation runs over all charges $q i j$ of the element j; $f a ( r i j )$ and $f d ( r i j )$ represent dimensionless coupling of the charge $q i j$ to $V m$ in the deactivated and activated states (also known as electrical distance). The coupling function can be approximated as the rate of change of the local electrostatic potential $φ ( r i j )$ with respect to $V m$⁠, $∂φrij∂Vm.$ In practice, when structures of the deactivated and activated states are available, $Q j$ can be computed by using this approximation and Eq. 3. However, for many membrane proteins, the structure of only one state is known, and therefore, direct estimation of $Q j$ is not possible. Instead, one can estimate the fraction of the gating charge, which is transferred by element j upon an infinitesimally small displacement $dr$ along the direction of the applied electric field. Hereafter, we will consider a scenario in which the applied electric field is aligned with the normal to the membrane, which coincides with the $z$ axis. Assuming $f( r i j )$ changes monotonically on $dz$⁠, one can obtain:

where $E z (r i j )$ corresponds to the local electric field. By using Eq. 4, it is possible to estimate the amount of charge transferred by the element j upon 1 Å displacement along the direction of an applied electric field,

We further show that $ζ j$ can be considered as a good predictor for $Q j$ and used for the detection of voltage-sensing elements.

Importantly, $ζ j$ reflects the propensity of the system to respond to changes in the local electric field upon application of $V m$⁠. It does not fully correspond to $Q j$ but to a fraction of it, assuming also that the element j is free to displace by 1 Å. In membrane proteins, however, conformational changes triggered by application of voltage can be substantially larger than 1 Å. In voltage-gated ion channels, for instance, S4 is displaced by ∼10 Å between the activated and resting states (Delemotte et al., 2011; Vargas et al., 2012). Such a large displacement is not well described by the first-order derivative (Eq. 4) assumed in our method. Therefore, the identified elements are state-specific; this means that other voltage-sensing elements may remain unidentified if they were not subject to a large local electric field in the examined conformation. Also, the list of identified elements may include false positives if not all of these elements are free to displace along the direction of the applied electric field during a voltage-sensing transition. Because of the potential presence of false-positives, it is important to further validate the predictions made by the suggested approach experimentally.

In addition, we estimate the contribution to $ζ j$ caused by the change in the local electric field, $∂Ezrij∂Vm.$ This property, which we call the local electric field response R, shows how sensitive the local electric field $E z (r i j )$ is to the application of an external electric field, and it is, therefore, useful for detection of protein regions subject to voltage sensitivity.

To test our approach, we applied it to a set of voltage-sensitive and voltage-insensitive membrane proteins. In the former case, we considered two voltage-gated ion channels, whose voltage-sensing elements have been well characterized, and six other membrane proteins, for which these elements remain largely unknown (Table 2). For each of these proteins, we calculated the local electric field response R and the system’s response propensity $ζ j$⁠. Fig. 3 shows the workflow of the suggested approach. A structure of the protein is embedded into a bilayer-solution system and equilibrated. Several protein conformations (N in total) are then extracted from the molecular dynamics trajectory, and each of them is further submitted to eight 2-ns-long runs under different values of an applied electric field (Fig. 3). After that, for every produced run, the local electrostatic potential φ is computed by using the PMEpot plugin of VMD (Aksimentiev and Schulten, 2005). R and $ζ j$ are finally calculated by using their definitions. Importantly, overall, we performed only N× 2 × 2 ns of molecular dynamics simulations (without considering the equilibration step), which only require modest computational resources.

In many voltage-sensitive membrane proteins, we detected regions with a large R value and therefore with a potential to respond to changes in $V m$ (Fig. 4 and Videos 1–10). The center of the voltage sensors in Shaker and NavMs (Long et al., 2005; DeCaen et al., 2009; Khalili-Araghi et al., 2010; Tao et al., 2010; Wu et al., 2010; Delemotte et al., 2011; Payandeh et al., 2011; Henrion et al., 2012; Jensen et al., 2012; Vargas et al., 2012; Lacroix et al., 2014; Machtens et al., 2017), the N terminus in voltage-dependent anion channel 1 (VDAC1) (Mirzabekov and Ermishkin, 1989; Peng et al., 1992; Thomas et al., 1993; Popp et al., 1996; Song et al., 1998; Teijido et al., 2014; Briones et al., 2016), and the Cl⁻ binding site in ClC1 (Pusch et al., 1995; Chen and Miller, 1996; Smith and Lippiat, 2010; Zifarelli et al., 2012; Grieschat and Alekov, 2014; De Jesús-Pérez et al., 2016) all show R values larger than 0.08. However, in the M₂ receptor (Ben-Chaim et al., 2003, 2006; Navarro-Polanco et al., 2011; Rinne et al., 2013; Vickery et al., 2016a,b), which is also known to be voltage-sensitive, R is small and reaches a maximum of only 0.04 in the Na⁺ binding site. Moreover, in voltage-insensitive membrane proteins, large R values were detected along the conductive pore, with the selectivity filter in the two-pore domain potassium channel (TWIK-1; Schewe et al., 2016) and the gate in the ligand-gated ion channel (GLIC; Bocquet et al., 2007), showing R values as large as 0.08. Therefore, R does not directly correlate with voltage sensitivity and cannot be used to discriminate between voltage-sensitive and voltage-insensitive proteins.

$ζ j$⁠, on the other hand, correlates rather well with voltage sensitivity: for all voltage-sensitive membrane proteins, $ζ j$ is large, while it is close to or below the detection threshold for all voltage-insensitive ones. Together with the previous observation, this suggests that voltage-sensitive proteins have evolved to place their charges in the regions with a large R value and therefore maximize their potential to respond to changes in $V m$⁠. In contrast, in voltage-insensitive proteins, charges are located far from the regions with large R values and thus do not sense changes in the local electric field. This also suggests our approach can be used to probe voltage sensitivity in other membrane proteins in which this property has not yet been tested.

We further describe the results obtained for the two voltage-gated ion channels, Shaker and NavMs. For voltage-gated ion channels in general, the mechanism of voltage sensitivity has been well characterized based on numerous experimental and computational studies (Long et al., 2005; DeCaen et al., 2009; Khalili-Araghi et al., 2010; Tao et al., 2010; Wu et al., 2010; Delemotte et al., 2011; Payandeh et al., 2011; Henrion et al., 2012; Jensen et al., 2012; Vargas et al., 2012; Lacroix et al., 2014; Machtens et al., 2017). These channels have four-helix bundle domains that sense changes in the MP–voltage sensor domains. One of the four helices (S4) carries several positively charged residues, while the other three (S1–S3) include a few negative countercharges (Long et al., 2005; Payandeh et al., 2011). Upon application of an electric field, the positive charges of S4 are displaced along the direction of this field to trigger conformational rearrangements in the pore domain (DeCaen et al., 2009; Khalili-Araghi et al., 2010; Tao et al., 2010; Wu et al., 2010; Delemotte et al., 2011; Henrion et al., 2012; Jensen et al., 2012; Vargas et al., 2012).

Our approach correctly detects the positively charged residues on S4 and their negative countercharges on S1–S3 as voltage-sensing elements (Fig. 5). Among them, residues located in the center of the voltage sensor domain show large $ζ j$ values, while those exposed to the extracellular or cytosolic solutions show $ζ j$ values close to the detection threshold (see Materials and methods). This indicates that, in the conformational states of Shaker and NavMs chosen for the analysis, residues such as R4 and K5/R5 contribute most to the gating charge, while the contribution from R1 and R2 is small.

A few other residues show $ζ j$ above the detection threshold, including R118 on the S4–S5 linker and E178 in the selectivity filter of NavMs, and a positive residue on S2 in both Shaker and NavMs. To our knowledge, the role of these residues in voltage sensitivity has not yet been assessed, and thus experimental verification will be required to confirm whether they are true- or false-positive signals.

The correct identification of voltage sensors in Shaker and NavMs strengthened our confidence in the predictive ability of our approach and suggested applying it to other membrane proteins, whose voltage sensors remain unknown. Overall, we applied it to six membrane proteins, including connexin-26 (Cx26; Maeda et al., 2009), transient receptor potential channel 1 (TRPV1; Cao et al., 2013; Liao et al., 2013; Kasimova et al., 2018), VDAC1 (Ujwal et al., 2008), voltage-gated chloride channel 1 (ClC1; Park and MacKinnon, 2018), muscarinic acetylcholine receptor M₂ (Haga et al., 2012), and Na⁺/K⁺ ATPase (Kanai et al., 2013; Figs. 6 and 7).

For Cx26, we detected two regions with large $ζ j$ values, including the extracellular opening (E42, D46, D50, K41, R75, R184, and K188), and the center of the protein subunit (E147, R32, and R143; Fig. 6). Many of these residues have already been suggested to play crucial roles in voltage sensitivity (Chen et al., 2005; Deng et al., 2006; Kwon et al., 2011, 2012; Pinto et al., 2016). For instance, molecular dynamics simulations (Kwon et al., 2011, 2012) suggested an electrostatic network among E42, D46, R75, R184, E187, and K188 to be a part of the voltage sensor. Experimental evidence indicates that mutagenesis of R75, K41, and E42 significantly alters Cx26 voltage sensitivity (Chen et al., 2005; Deng et al., 2006; Pinto et al., 2016). In particular, K41 neutralization results in an increase of the apparent gating charge, while the double mutant K41E/E42S was shown to be more sensitive to voltage than the wild type (Pinto et al., 2016). Although experimental evidence suggests D2 also contributes to voltage sensitivity (Verselis et al., 1994), this residue did not exhibit any large $ζ j$ value in our analysis. We believe this discrepancy might be caused by the low resolution of the N terminus, where D2 is located (Maeda et al., 2009): shortly after the start of molecular dynamics simulations, the N terminus lost its secondary structure and significantly deviated from its initial position, which limits the predictive power in this region.

TRPV1 has a four-helical domain, which is similar to the voltage-sensor domain of Shaker and NavMs (Cao et al., 2013; Liao et al., 2013). However, because this channel does not include several crucial charges on S4 and the corresponding countercharges on S1–S3 (Cao et al., 2013; Liao et al., 2013), it likely responds to changes in the MP through a different mechanism. Our analysis reveals that, in addition to the lack of charges, the four-helical domain also shows small R values, emphasizing that this region is not able to sense changes in the MP (Fig. 4). On the other hand, we found $ζ j$ values above the detection threshold for residues of the extra- and intracellular openings of the four-helical domain (R455, R474, E478, R491, and E513) and of the S4–S5 linker (R557, E570, D576, and R579; Fig. 6). Mutagenesis of several of these residues (R557, E570, D576, and R579) has already been shown to significantly affect TRPV1 voltage sensitivity (Boukalova et al., 2010).

In VDAC1, we observed multiple residues spread over the entire protein with large $ζ j$ values (Fig. 6). These include 5 residues on the α-helical N terminus and 27 residues on the β-barrel. Experimental and computational evidence also indicates that both the N-terminal α-helix and the β-barrel contribute to voltage sensitivity (Peng et al., 1992; Thomas et al., 1993; Popp et al., 1996; Song et al., 1998; Briones et al., 2016). Together with our data, this suggests a scenario in which the entire protein is involved in the response to changes in the MP rather than a single domain playing the role of a voltage sensor. Among the residues pinpointed by our analysis, mutagenesis of D16, K20, K61, and E84 is known to modulate the steepness of the conductance-voltage relationship (Thomas et al., 1993), while that of D32 and K96 does not significantly affect voltage sensitivity (Thomas et al., 1993). This indicates there are likely a few false-positive signals in our analysis; therefore, it has to be combined with experiments to confirm voltage sensitivity of the detected elements.

Previous experimental studies performed on ClC channels and transporters have revealed two sources of voltage sensitivity: protein residues and chloride ions (Pusch et al., 1995; Chen and Miller, 1996; Smith and Lippiat, 2010; Zifarelli et al., 2012; Grieschat and Alekov, 2014; De Jesús-Pérez et al., 2016). Reducing extracellular Cl⁻ concentration in ClC-0 or a substitution of intracellular Cl⁻ by impermeant anions in ClC-2 shifts their opening to more positive voltages (Pusch et al., 1995; Chen and Miller, 1996; De Jesús-Pérez et al., 2016). Measurements of the gating charge in a permeation-deficient mutant of ClC-5 revealed that Q is reduced upon extracellular ion depletion but does not go to zero, indicating that protein residues contribute to voltage sensitivity along with ions (Smith and Lippiat, 2010; Zifarelli et al., 2012; Grieschat and Alekov, 2014). One of these residues could be the so-called glutamate gate (E232 in ClC1), whose mutagenesis completely abolishes voltage sensitivity (Dutzler et al., 2003; Estévez et al., 2003; Melzer et al., 2003). Our analysis shows large $ζ j$ values for both E232 and bound Cl⁻ ions (Fig. 7). It also pinpoints several residues along the conduction path (D136, K215, K231, R421, and K467) and a region at the interface between the two ClC1 subunits (E291, R317, E347, and E548); these residues all show $ζ j$ values above the detection threshold, suggesting they are likely significant contributors to the ClC1 gating charge.

In the case of the M₂ receptor, the identified voltage-sensing element is complex as well and composed of the protein residues and a Na⁺ ion bound inside the internal hydrophilic pocket. The role of a Na⁺ ion in M₂ receptor voltage sensitivity was in fact recently suggested based on molecular dynamics simulations (Vickery et al., 2016a). Our data agree with these findings and show large $ζ j$ values for a Na⁺ ion and for D69 and D107 (Fig. 7). In addition, our results support the recent experimental finding that the DRY motif on the TM3 segment does not contribute to voltage sensitivity (Navarro-Polanco et al., 2011): direct measurements of the gating charge revealed that the Q value of the wild-type and the D120N-R121N mutant is similar.

Finally, for the Na⁺/K⁺ ATPase, we found three regions where $ζ j$ is above the detection threshold, including the ionic binding site and the extracellular and cytosolic openings (Fig. 7). The ionic binding site, which is composed of E327, D808, D926, E953, and three Na⁺ ions, shows the largest contribution to the gating charge. This agrees with the known data reporting that changes in the extracellular Na⁺ concentration has a strong effect on voltage dependence of the Na⁺/K⁺ ATPase ionic current (Nakao and Gadsby, 1989) and the gating charge (Holmgren et al., 2000).

For the set of membrane proteins investigated here, we were able to rediscover several elements that were previously shown to play a role in voltage sensitivity. It is important to note that in ClC1 the M₂ receptor and Na⁺/K⁺ ATPase voltage sensitivity was accessed directly by estimating the change in the gating charge upon mutagenesis of residues or upon variation of the ionic concentrations (Nakao and Gadsby, 1989; Holmgren et al., 2000; Ben-Chaim et al., 2006; Smith and Lippiat, 2010; Navarro-Polanco et al., 2011; Zifarelli et al., 2012). In Cx26, VDAC1, and TRPV1, on the other hand, it was accessed indirectly: as a reporter for voltage sensitivity, the effect of a given mutation on the conductance-voltage relationship (⁠$V 1/2$ and the voltage-independent component of gating) was used. Such an effect, however, can be possibly caused by structural perturbations of the voltage sensor itself or of other protein domains, and therefore, the experiments indirectly accessing voltage sensitivity must be interpreted with caution.

Traditional mutagenesis tends to introduce perturbations that simultaneously affect electrostatic and steric properties. One can disentangle these effects on protein function by using nonnatural amino acids. For example, citrulline is a neutral and nearly isosteric analogue of arginine that has recently been used to quantify electrostatic contribution of the voltage-sensing residues to gating (Infield et al., 2018). Because estimation of the gating charge is often challenging if at all possible, the usage of nonnatural amino acids can emerge as an alternative strategy to characterize voltage sensitivity.

One of the reasons voltage sensitivity in membrane proteins is hard to detect is that it depends on a combination of two factors: the ability of a protein to focus an electric field in some of its regions and its ability to respond to the electric field based on the presence of charged elements in these regions. We show that it is possible to use a simple and cheap computational tool to distinguish between voltage-sensitive and voltage-insensitive proteins and to detect true voltage-sensitive residues in voltage-gated ion channels. The application of this tool to six membrane proteins without well-characterized voltage sensors leads to the prediction of several candidates for voltage-sensing elements, some of which have already been shown to contribute to voltage sensitivity while others remain to be tested experimentally. The suggested approach is general, does not require extensive computational resources, and can be applied to any other membrane protein regardless of structure or function. This application is important to unravel the variety of ways voltage sensitivity has been manifested in biological molecules, and it should be highly useful both to identify how different factors such as mutations, membrane composition, or ligand binding might alter voltage sensitivity and to engineer this property into voltage-insensitive membrane proteins.

We acknowledge Kawai Lee, Marie Lycksell, and Rebecca J. Howard for providing the GLIC trajectory, Asghar Razavi for providing the Na⁺/K⁺ ATPase trajectory, and Samira Yazdi for providing the Shaker trajectory. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing at PDC Centre for High Performance Computing.

This work was supported by grants from the Vetenskapsrådet (2017-04641) and the Science for Life Laboratory.

The authors declare no competing financial interests.

Author contributions: M.A. Kasimova, E. Lindahl, and L. Delemotte designed the research. M.A. Kasimova performed the simulations, analyzed the data, and wrote the original draft. M.A. Kasimova, E. Lindahl, and L. Delemotte edited the paper.

José D. Faraldo-Gómez served as editor.