Binding of the general anesthetic sevoflurane to ion channels

The direct-site hypothesis assumes general anesthetics bind ion channels to impact protein equilibrium and function, inducing anesthesia. Despite advancements in the field, a first principle all-atom demonstration of this structure-function premise is still missing. We focus on the clinically used sevoflurane interaction to anesthetic-sensitive Kv1.2 mammalian channel to resolve if sevoflurane binds protein’s well-characterized open and closed structures in a conformation-dependent manner to shift channel equilibrium. We employ an innovative approach relying on extensive docking calculations and free-energy perturbation of all potential binding sites revealed by the latter, and find sevoflurane binds open and closed structures at multiple sites under complex saturation and concentration effects. Results point to a non-trivial interplay of site and conformation-dependent modes of action involving distinct binding sites that increase channel open-probability at diluted ligand concentrations. Given the challenge in exploring more complex processes potentially impacting channel-anesthetic interaction, the result is revealing as it demonstrates the process of multiple anesthetic binding events alone may account for open-probability shifts recorded in measurements.


Introduction
Volatile and injected general anesthetics encompass a diverse array of small and uncharged chemotypes including haloalkanes, haloethers and alkylphenols. Despite efforts reaching back over a century, clarification of their microscopic mechanism in general anesthesia has proven difficult and wanting. A favored hypothesis proposes that ion channels in the brain are implicated, among which members of ionotropic neurotransmitter receptors, voltage-gated and non-gated ion channels are best-known players [1][2][3]. Primary exemplars are the Cysloop nicotinic acetylcholine and γ-aminobutyric acid class A receptors, the voltage-gated sodium and potassium channels, and the tandem pore potassium channels. An extensive series of electrophysiological studies corroborate the hypothesis by demonstrating a range of effects, from inhibition to potentiation, of general anesthetics on the various receptor targets. Beyond these electrophysiological studies of reductionist systems, the current view has gained additional support from gene knockout experiments demonstrating for some of these channels the in vivo role on a clinically-relevant anesthetic outcome. For instance, the knockout of the nongated tandem pore potassium channel trek-1 produces an animal model (Trek1-/-) resistant to anesthesia by inhalational anesthetics [4]. How general anesthetics modulate ion channels to account for endpoints of anesthesia must at some point build on understanding electrophysiological data in the context of ligand binding, a reasoning that has driven mounting efforts in the field. Currently, though not refuting other molecular processes likely to contribute to anesthetic action [5][6][7], crystallography and molecular dynamics studies support that anesthetics bind ion channels at clinical concentrations [8][9][10][11][12][13][14][15][16]. Binding interactions have been evidenced in anesthetic containing systems of mammalian voltage and ligand-gated channels, as well as bacterial channel analogs. Specifically, partitioning of anesthetics in the membrane core allows it to access and bind multiple transmembrane (TM) protein sites, featuring single or multiple occupancy states-a process that might depend further on chemotypes, channel types and conformations. Although some progress has been made in validating one or more aspects of the direct-site hypothesis, a firstprinciple demonstration that anesthetics bind ion channels to affect protein equilibrium and function as recorded in experiments is still unaccounted for.
Here, we focus our efforts on the haloether sevoflurane and its molecular interaction to Kv1.2, a mammalian voltage-gated potassium channel. Experimental work demonstrates that sevoflurane potentiates the channel in a dose-dependent manner [3,17,18]. Effects on current tracings include a leftward shift in the channel's conductance-voltage relationship and an increased maximum conductance. As extensively discussed in these past publications, at least two molecular mechanisms are expected to be involved in Kv channels potentiation by sevoflurane. One mechanism (i) might involve sites allosterically coupled to the electromechanical transduction directly responsible for controlling voltage-dependent gating. The other (ii) might involve distinct sites, which could modulate the channel's pore region and influence the stability of the conductive state and/or the unitary conductance. Here, we are interested in the investigation of mechanism (i) and its underlying structural hypothesis that sevoflurane binds the channel's open and closed states to impact protein equilibrium and therefore its voltage dependence. Among all other aspects that might impact channel-anesthetic interactions in general, our specific goal is to determine if sevoflurane binds the well-characterized open-conductive (O) and resting-closed (C) structures of Kv1.2 [19,20] in a conformation-dependent manner to impact its voltage-dependent open probability as recorded experimentally. Very recently, we have put forth an innovative structure-based study [21] dealing with the concentration-dependent binding of small ligands to multiple saturable sites in proteins to show that sevoflurane binds the open-pore structure of Kv1.2 at the S4S5 linker and the S6P-helix interface-a result largely supported by independent photolabeling experiments [22,23]. To our current goal, we aim therefore at extending these calculations to investigate sevoflurane interactions with the entire channel TM-domain and, more importantly, to resolve any conformational dependence in its binding process to channel structures C and O. In the following sections, we first provide the theoretical framework to study sevoflurane binding to a specific channel conformation under equilibrium conditions. A state-dependent strategy is put forward to describe anesthetic binding in terms of occupancy states of all identified potential channel binding sites, embodying both concentration and multiple sites saturation effects. The strategy is then generalized to account for ligand impact in the C-O equilibrium, allowing for reconstruction of voltage-dependent open probabilities of the channel at various ligand concentrations. Anticipating our results, we find that sevoflurane binds Kv1.2 structures at multiple sites under saturation and concentration effects. Despite a similar pattern of molecular interactions, binding of sevoflurane is primarily driven towards the open-conductive state shifting leftward the open probability of the channel at diluted ligand concentrations.

Binding of anesthetics to multiple channel sites
We applied large-scale and flexible docking calculations to solve sevoflurane interactions to Kv1.2 structures X � {C, O} (Fig 1). A total of~6,000 docking solutions was generated per channel conformation and clustered into 21 ligand interaction sites. The interaction sites spread over the transmembrane region of the channel at the S4S5 linker, S6P-helix interface and at the extracellular face, next to the selectivity filter. Further docking sites were resolved within the voltage-sensor, at the S4Pore interface and within the channel central cavity. Redocking of sevoflurane generated in turn a total of~13,000 solutions per channel conformation, solving the interaction of two ligands for all sites but the extracellular face.
From the docking ensembles, there are up to 2 × 3 21 channel occupancy states that might contribute to sevoflurane binding and functional effects. To quantitatively evaluate this, we performed an extensive series of decoupling FEP calculations to estimate the per-site binding affinity for one and two bound ligands against the channel structures (cf. Computational Methods for details). Here, FEP calculations started from equilibrium ligand-bound channel structures, embedded in an explicit water-membrane environment (cf. RMSD analysis in S1 Under these technical details, per-site equilibrium binding constants were quantified relative to a homogeneous and diluted aqueous solution occupied by ligands, with an excess chemical potential of m ¼ 0:10 � 0:09kcal:mol À 1 . As shown in S1 and S2 Tables, per-site binding constants are heterogeneous and take place over a diverse range, i.e. 10 −8 (mM -1 ) -10 +2 (mM -2 ). There is however a decreasing trend of affinities involving sites respectively at the S4S5 linker, S4Pore and S6P-helix interfaces, voltage sensor, central cavity and extracellular face.
To determine if sevoflurane binds channel structures X � {C, O} at clinically relevant concentrations, we computed binding probabilities ρ X (n 1 ,. . .,n s ) for dilute concentrations of the ligand in solution, i.e. 1mM, 10mM and 100mM. Equilibrium constants K X (n 1 ,. . .,n s ) for every occupancy state of the channel were then reconstructed from the per-site affinities to determine state probabilities via eq (2). Here, estimates of K X (n 1 ,. . .,n s ) were determined for the condition of independent binding sites, as the minimum site-to-site distances of~15 Å demonstrated their non-overlap distributions in each of the channel structures (cf. Computational Methods for details). At low 1mM concentration, ρ X (n 1 ,. . .,n s ) is largely dominated by the empty state probability ρ X (0 1 ,. . .,0 s ) indicating only a small fraction of bound states with nonnegligible occurrences (S4 Fig). Within this fraction, the most likely states involve single occupancy of the S4S5 linker or the S4Pore interface as shown by the marginal probabilities ρ X (n j ) of individual sites (Fig 2). At higher concentrations, there is a clear shift of ρ X (n 1 ,. . .,n s ) towards channel occupancy states that significantly enhance the average number of bound ligands. Careful inspection of ρ X (n j ) confirms the major relevance of sites at the S4S5 linker and S4Pore interface over the entire concentration range, accompanied by an increasing importance of binding regions at the S6P-helix interface. In contrast, ρ X (n j ) for sites within the voltage-sensor, central cavity and nearby the extracellular face of the channel remain negligible over all concentrations. For completeness, note in S1 Table that equilibrium constants for doubly-occupied sites are comparable to or even higher than estimates for one-bound molecule thus revealing important saturation effects in which one or two sevoflurane molecules can stably bind the channel structures at individual sites. The result is especially true for spots at the S4S5 linker and S4Pore interface.
The complex distributions of the multiple occupied states of structures X � {C, O} were described in three dimensions by mapping ρ X (n 1 ,. . .,n s ) into the position-dependent density r j X ðRÞ of sevoflurane in each binding site j (cf. Computational Methods for details). As shown in Fig 3 and supplementary S1 and S2 Movies, the density of sevoflurane better convey the results by showing the spatially-mapped concentration dependent population of bound ligands. Projection of r j X ðRÞ along the transmembrane direction z of the system, r j X ðzÞ, stresses further the results. Note from r j X ðRÞ that sevoflurane binds channel structures in a concentration dependent manner, binding preferentially the S4S5 linker and the interfaces S4Pore and S6P-helix over a range of concentrations.
So far, our calculations demonstrate that sevoflurane binds Kv1.2 structures over a spectrum of concentrations, preferentially at the linker S4S5 and at the segment interfaces S4Pore and S6P-helix. From a physical-chemical point of view, spots at these channel regions are primarily dehydrated, lipid accessible, amphiphilic pockets providing with favorable interaction sites for the polar lipophilic sevoflurane molecule (S5 Fig). It is worth mentioning that these findings recapitulate recent photolabeling experiments demonstrating that photoactive analogs of sevoflurane do interact to the S4S5 linker and at the S6P-helix interface of the open-  [22,23]. In detail, Leu317 and Thr384 were found to be protected from photoactive analogs, with the former being more protected than the latter. As shown in S6 Fig, atomic distances of bound sevoflurane to these amino-acid side chains are found here to be respectively 7.28±2.5 Å and 10.44±3.66 Å, in average more or less standard deviation. Such intermolecular distances are consistent with direct molecular interactions and therefore consistent with the measured protective reactions-similar conclusions hold for the closed channel as well. Besides that, our calculations recapitulate the stronger protection of Leu317 in the sense that, relative to sites at S6P-helix, the affinity of sevoflurane is found to be higher at the S4S5 linker considering its stable occupancy either by one or two ligands. The stable occupancy of the linker by one or two ligands as computed here, is also consistent with recent flooding-MD simulations of the homologous sodium channel NaChBac [14,24] and more importantly, with previous Ala/Val-scanning mutagenesis showing a significant impact of S4S5 mutations on the effect of general anesthetics on members of the K + channel family [10]. Densities pertaining to sites S4S5 linker, S6P-helix and S4Pore are indicated with yellow symbols. As described in Computational Methods, the determination of r j X ðRÞ involved reweighing the marginal probability ρ X (n j ) at the binding site j by the local equilibrium density of sevoflurane ρ X (R|n j ). The marginal ρ X (n j ) was computed by coarse-graining over state probabilities in S2  In particular, a single residue (Gly329) at a critical pivot point between the S4S5 linker and the S5 segment underlies potentiation of Kv1.2 by sevoflurane [18]. Sevoflurane is found to be close to that amino acid when bound to the S4S5 linker.
In contrast to the aforementioned spots, sites within the voltage-sensor, within the main pore and nearby the extracellular face of the Kv1.2 structures are primarily hydrated, lipidinaccessible, amphiphilic pockets (S5 Fig) that weaken sevoflurane interaction as reflected in the state-and space-dependent densities shown in Figs 2 and 3. The binding probabilities at these sites thus support that a non-negligible fraction of poses determined from docking ( Fig  1D) corresponds to low affinity or false positives. In particular, because sevoflurane induces potentiation rather than blocking of Kv1.2 [17,18], we read the negligible or absent density of the ligand in the channel central-cavity as a self-consistent result of the study-especially for the open-conductive state. Supporting that conclusion, note that binding constants as computed here are upper bounds for the affinity of sevoflurane under ionic flux conditions in which potentiation takes place. Accordingly, as shown in S7 Fig, the binding affinity of a potassium ion at the central cavity overcomes that of sevoflurane due its binding and excess free-energies under applied voltages. Once bound, the ion destabilizes sevoflurane interactions and the molecule is not expected to bind the channel cavity at low concentrations. As also shown in S7 Fig and supplementary S3 Movie, even under the occurrence of rare binding events, sevoflurane appears unable to block the instantaneous conduction of potassium which is also consistent with its potentiating action.
Weak interactions at the main pore and nearby the selectivity filter of Kv1.2 contrasts with sevoflurane binding at analogous regions of NaChBac [14,24], likely due major structural differences between Na + and K + channels. Specifically, the pore of potassium channels lacks lipid-accessible open-fenestrations of the sodium relatives and K + -selective filters are sharply distinct from Na + -selective ones.

Anesthetic binding impacts channel energetics
Despite a comparable pattern of molecular interactions, careful inspection of ρ X (n j ) or r j X ðRÞ reveals that for most sites there is an obvious differential affinity of sevoflurane across Kv1.2 structures (Figs 2 and 3). The overall consequence for sevoflurane binding is then clear: the average number of ligands bound to the open-conductive channel systematically exceeds that for the resting-closed channel over the entire concentration range. There is therefore a remarkable conformational dependence for the anesthetic interaction, such that sevoflurane preferentially binds the open-conductive structure.
Implications for Kv1.2 energetics were then investigated by quantifying changes to the channel open probability ρ O (V) induced by sevoflurane at concentrations of 1mM-100mM (Fig 4). Specifically, from the partition functions Z C (n 1 ,. . .,n s ) and Z O (n 1 ,. . .,n s ) across the entire ensemble of occupancy states of the channel, solution of eqs (5) and (9) show that sevoflurane shifts leftwards the open probability of Kv1.2 in a concentration-dependent manner-voltage shifts amount from -1.0 mV to -30.0 mV with concentration increase of the ligand in solution. The result is particularly interesting, supporting that the approximately 3 mV probability shift recorded experimentally at 1 mM sevoflurane concentration can be explained by our structurebased probability predictions in the concentration range of 1-10 mM. The latter thus provides a theoretical basis to predict sevoflurane impact on channel energetics in a larger, not yet experimentally probed, range of concentrations. Additionally, for a fixed ligand concentration (100 mM), decomposition analysis reveals further that ratio values for the partition functions at individual sites j can be smaller, equal or larger than unity, implying a non-trivial interplay of conformation-dependent modes of action involving distinct sites (cf. Computational Methods for  (5) and (9) by taking into consideration parameters, V m = −21.9mV and ΔQ = 3.85e o , for best two-state Boltzmann fit of measured data for Kv1.2 free of ligands [18]. A reference experimental curve (blue) is shown for sevoflurane at 1 mM concentration, with best two-state Boltzmann parameters V m = −25.1mV and ΔQ = 4.00e o [18]. Relative to the ligand- . In detail, binding of sevoflurane at low affinity sites within the voltage-sensor, central cavity and next to the extracellular face of the channel are mostly conformation-independent and do not impact open probability (ratio � 1). On the other hand, conformation-dependent binding of sevoflurane to sites at the S4S5 linker and the S4Pore interface accounts for the overall stabilization of the open channel (ratio < 1). That effect contrasts with the mild stabilization of the closed conformation of Kv1.2 induced by binding of sevoflurane at S6P-helix and reflected in rightward shifts of ρ O (V) (ratio > 1). The overall conformation-dependent binding process is therefore differentially encoded across distinct channel regions.
As extensively discussed in past publications [3,17,18], potentiation of Kv1.2 by sevoflurane has been attributed to stabilization of the open-conductive state of the channel via at least two molecular mechanisms. One mechanism (i) likely involving sites allosterically coupled to the electromechanical transduction responsible for controlling voltage-dependent gating; and another (ii) implicating distinct sites, which could influence the pore conductive state stability and/or unitary conductance. Here, our structural calculations based on the open-activated and resting-closed states of Kv1.2 were consistently designed to investigate the first (i) of these mechanisms. Given the critical role of S4 and S4S5 linker on the channel gating mechanism [19], it is reasonable that sevoflurane interactions with these segments, as found here, are at the origins of the experimentally measured voltage-dependent component of anesthetic action. While restricted to sevoflurane interactions with the resting-closed and open-conductive structures, the presented two-state binding model only embodies left or rightward shifts in the open probability of the channel, therefore it cannot clarify any molecular process accounting for the maximum conductance increase recorded experimentally. As supported by a recent kinetic modeling study [17], generalization of eqs (9) to include a third non-conducting open state yet structurally unknown is needed to account for such conductance effects and for that reason, the investigation of mechanism (ii) is beyond the scope of our study. We speculate however that binding of sevoflurane at the S4Pore and S6P-helix interfaces could allosterically interfere with pore domain operation, thus affecting channel's maximum conductance. A working hypothesis also raised in the context of anesthetic action on bacterial sodium channels [12,14], assumes indeed that non-conducting states of the selectivity filter are implicated. Corroboration of such an assumption from a molecular perspective is however not trivial and will necessarily involve further structural studies to demonstrate how ligand binding might impact non-conducting open states of the channel to affect maximum conductance.

Discussion
Here, we carried out extensive structure-based calculations to study conformation-dependent binding of sevoflurane to multiple saturable sites of Kv1.2 structures X �{C, O} under equilibrium conditions-the total MD simulation time was~2.0 μs. Binding of sevoflurane was studied for ligand concentrations in the range of 1mM-100mM and saturation conditions up to n max j ¼ 2. Our study relied on the assumption that molecular docking calculations performed binding sites j. Per-site ratio values can be equal, smaller or larger than unity meaning respectively that sevoflurane binding is not conformational dependent, stabilizes the open structure or stabilizes the closed structure. Binding sites located nearby flexible protein regions for which the root-mean-square deviation (RMSD) between channel structures is larger than 4.0 Å are highlighted in red (cf. Computational Methods and S2 Fig for details) in vacuum can faithfully describe ligand interactions at protein sites. Specifically related to that assumption, we have considered the generated ensemble of docking solutions to estimate the location of binding sites δV j and the local distribution of the ligand ρ X (R|n j ). The generation of false positive hits is however a well documented drawback of docking algorithms as a result of limitations of the scoring function in describing ligand solvation energies and protein flexibility [25]. Given the same limitations of the scoring function, it is also not guaranteed that neither all binding hits, nor that ρ X (R|n j ) can be accurately known from docking. In this regard, although not considered here, it might be important to integrate docking results from different algorithms involving different scoring functions in order to characterize the bound ensemble. Still, thanks to the generality of the presented formulation, extension of the current investigation to sampling techniques other than docking, including all-atom flooding-MD simulations [9,11,12,14,16], might also be an important refinement in that direction.
Despite these sampling improvements that may eventually be obtained, the presented combination of extensive docking calculations against an ensemble of equilibrium receptor structures fit to handle protein flexibility, and FEP calculations based on fine force-fields to accurately estimate solvation energies are critical technical aspects of the applied methodology devised to minimize such drawbacks. Whereas docking was performed in vacuum, FEP calculations were carried out in presence of explicit all-atom lipids and water, therefore taking into consideration environmental effects in the estimated ligand binding free energies. In this regard, standard binding free-energies estimated from FEP (S1 and S2 Tables) are comparable or significantly smaller than the standard free energy for transferring the ligand from water to a pure lipid bilayer-supporting that sevoflurane is expected to partition preferentially into channel sites rather than the membrane. The partition coefficient (log K) and the related transfer free energy of sevoflurane between water and the lipid bilayer (POPC) amount respectively to 2.64±0.96 and -3.12±0.32 kcal.mol -1 as recently estimated by Tajkhorshid and coworkers [26]. Besides that, it is also important to note that the configuration space in FEP calculations overlap between channel structures at individual sites, i.e. sampling and binding affinities were evenly resolved between states (S2 Fig and Fig 4B)-meaning eventual biases or systematic errors were mitigated when comparing similar calculations between channel states according to the main goal here.
Under these technical considerations, we conclude that most of the identified binding sites are located nearby flexible protein regions for which the root-mean-square deviation between channel structures is larger than 4.0 Å. Then for the purpose of quantifying any direct ligand effect on channel energetics, the determined conformational dependence of sevoflurane binding to these gating-implicated protein regions appears robust and likely to impact function. Structural knowledge allied to solid electrophysiological data available for Kv1.2 make this channel an interesting model system for molecular-level studies of anesthetic action thereby justifying our choice. In detail, the atomistic structures account for most of the available experimental data characterizing closed and open conformations of the channel in the native membrane environment [20]. Previous findings support further that sevoflurane binds Kv1.2 to shift leftward its voltage-dependence and to increase its maximum conductance in a dosedependent manner [18]. Despite a similar pattern of interactions, we found here a clear conformational dependence for sevoflurane binding at multiple channel sites. The ligand binds preferentially the open-conductive structure to impact the C-O energetics in a dose-dependent manner as dictated by the classical equilibrium theory for chemical reactions embodied in eq (9). Front of the difficulty in conceiving and characterizing other, still more complex molecular processes that might impact channel energetics under applied anesthetics [5][6][7], the result is revealing by showing that in principle the isolated process of sevoflurane binding to Kv1.2 accounts for open-probability shifts as recorded in experiments. Within this scenario, the calculations reveal contrasting per-site contributions to the overall open probability of the channel. For instance, at 100mM concentration, binding of sevoflurane at the S4S5 linker and S4Pore interface significantly stabilizes the open structure of the channel overcoming the mild stabilization of the closed structure by ligand binding at the S6P-helix interface. By showing this non-trivial interplay of conformation-dependent modes of action involving distinct binding sites, the result is particularly insightful and should guide us to design novel site-specific mutagenesis and photolabeling experiments for further molecular characterization of anesthetic action.
Although not addressing the paucity of in vivo experimental evidences that a binding process to a specific molecular target as presented here is related to any clinically-relevant anesthetic outcome, our study adds support to the direct-site hypothesis by linking binding freeenergy and protein energetics. As such, our study treats and reveals a new layer of complexity in the anesthetic problem that brings us novel paradigms to think their molecular action and to design/interpret research accordingly. To the best of our knowledge, the main-text Figs 3 and 4 represent in the context of structural studies, a deeper and first revealed view on the intricate mode of interactions that might take place between general anesthetics and ion channels to impact function in general.

Theory
Anesthetic binding and channel energetics. Consider the voltage-gated channel embedded in a large membrane-aqueous volume that contains N ligand molecules under dilution. The protein is assumed to remain in a well-defined conformational state X, in which it presents s distinct binding sites for ligands. For simplicity, we consider that ligands dissolve uniformly across the membrane-aqueous region of the system from where they can partition into the protein sites. The lipid and aqueous phases thus provide with a bulk volume V occupied by ligands at constant density r and excess chemical potential m. We consider further that every site j = 1,. . .,s corresponds to a discrete volume δV j that can be populated by 0 ⩽ n j ⩽ n max j ligands. We denote by O � X ðn 1 ; . . . ; n s Þ the specific occupancy state featuring n j bound ligands at corresponding sites and by n = n 1 +. . .+n s the total number of bound ligands in this state.
Under these considerations, solution of ligand binding to multiple receptor sites relies fundamentally in determining the equilibrium constant K X (n 1 ,. . .,n s ) for the process O . . . ; 0 s Þ is the empty receptor state with all ligands occupying the bulk. As shown in previous work [21], at a fixed temperature β = (k B T) −1 , the binding constant K X (n 1 ,. . .,n s ) can be evaluated from MD-based free-energy perturbation (FEP) calculations in which m is the solvation free energy of the ligand in the bulk and W � X ðnÞ corresponds to the free-energy of n site-specific bound ligands relative to a gas phase state given that ligands i = 1,. . .,n are restrained with force constants k X (i) to occupy an effective site volume to be known in practice from free-energy calculations [27]. Here, the normalization condition appearing in the denominator of eq (2) runs over every occupancy state Oðn 0 1 ; . . . ; n 0 s Þ of the channel, ranging from O(0 1 ,. . .,0 s ) up to Oðn max 1 ; . . . ; n max s Þ. Note in eq (2) ρ X (n 1 ,. . .,n s ) depends on the number density or concentration of the ligand at the reservoir thus providing a useful equation for investigation of concentration effects.
To investigate any conformational dependence on ligand binding, we consider eq (2) in the context of conformational equilibrium of the channel over a range of TM voltages. Specifically, we assume the very same microscopic system submitted to a Nernst potential induced by nonsymmetrical electrolytes between membrane faces. The capacitive nature of the channel-membrane system ensures the Nernst potential accounts for a voltage difference V across the lipid bilayer. Accordingly, by denoting the entire set of instantaneous Cartesian coordinates of the channel as r P , the free energy of the protein F X (V) in a particular conformation X � X(r P ) e À bF X ðVÞ / R dr P d½X 0 ðr P Þ À X�e À b½Uðr P ÞþQðr P ÞV� ð3Þ can be written within an arbitrary constant, in terms of an effective potential energy of the protein U(r P ) + Q(r P )V when coupled to the external voltage V with charge Q(r P ) [28]. Note the integral runs over the entire configurational space accessible to r P , so long as it is compatible with conformation X, as indicated by the delta notation δ[X 0 (r P ) − X]. From eq (3), the open probability of the channel then reduces to in which, is the gating charge ΔQ = Q O − Q C resulting from differences in the effective protein charge in each conformational state and is the midpoint voltage in which ρ C (V) = ρ O (V) [28]. From eq (5), the equilibrium constant between protein states C and O then writes as KðVÞ ¼ e À bΔQ½V m À V� with Kð0Þ ¼ e À bV m ΔQ determining their equilibrium at 0 mV. In eqs (6) and (7), the voltageindependent free energies account for the microscopic potential energy of the channel and its solvation energy in each state whereas the corresponding voltage-dependent excess free energies are proportional to the applied voltage and associated protein charges. By combining eqs (2) and (5) through a generalized thermodynamic-cycle analysis dealing with all possible states of the ligand-free and ligand-bound receptor, binding effects on the channel energetics can be then explicitly expressed over a range of membrane voltages for the ensemble of occupancy states in each protein conformation. Eq (8) simplifies into the two-state Boltzmann equation now embodying the free-energy contributions arising from ligand binding (cf. Computational Methods for details).
In this contribution, we consider eqs (1), (5) and (9) to investigate the molecular binding of sevoflurane to open and closed structures of Kv1.2, and its functional impact on channel energetics.

Computational methods
A procedure was designed to solve the molecular binding of sevoflurane to the open-conductive (O) and resting-closed (C) structures of Kv1.2 for saturation conditions up to n max j ¼ 2, under equilibrium conditions. For both channel structures, the procedure consisted of (i) an extensive production of docking solutions for the ligand-receptor interaction, (ii) clustering of docking solutions into binding sites along the receptor structure and (iii) estimation of binding affinities using the free-energy perturbation (FEP) method. First completion of steps (i) through (iii) solved the ligand channel interaction for singly-occupied binding sites. Double occupancy of the receptor sites was investigated by inputting the first generated ensemble of docked structures into another round of (i) through (iii) calculations. In detail, step (i) was accomplished by docking sevoflurane as a flexible ligand molecule against an MD-generated ensemble of membrane-equilibrated structures of the protein receptor. Docking calculations included the transmembrane domain of the channel, free from the membrane surroundings.
Step (ii) provided the location of δV j volumes lodging docking solutions for the ligand along the channel structures. Each of these volumes was treated as binding site regions in FEP calculations. Here, FEP calculations started from equilibrium structures of the ligand-bound channel embedded in a explicit water-membrane environment. For the purpose of improving statistics, FEP estimates and the associated statistical errors were determined from at least two independent decoupling runs. Calculations were performed over~7.0 ns per site, per conformation and per replica to converge FEP estimates, in a total MD simulation time of~2.0 μs. Following this procedure, binding constants K X (n 1 ,. . .,n s ) for channel structures X � {C, O} were solved by inputting FEP estimates into eq (1), allowing for direct solution of state-dependent probability distributions ρ X (n 1 ,. . .,n s ) via eq (2). Here, binding constants K X (n 1 ,. . .,n s ) and the related standard binding free energies ΔG o X ðn 1 ; . . . ; n s Þ were solved for the condition of independent binding sites and relative to a homogeneous and diluted aqueous solution occupied by ligands at constant density r and excess chemical potential m. Ligand-free and ligand-bound open probability curves ρ X (V) were respectively computed from eqs (5) and (9) by taking into consideration the previously determined mid-point voltage and steepness of the open probability curve of Kv1.2 free of ligands, i.e. V m = −21.9 mV and ΔQ = 3.85e o as determined from best two-state Boltzmann fit of the measured conductance-voltage data of the channel [18]. Note that any ligand-induced shift in eq (9) is determined by the partition function ratio between open and closed structures and not by the choice of these reference parameters. Probabilities ρ X (n 1 ,. . .,n s ) and ρ X (V) were determined for sevoflurane concentrations in the range of 1mM-100mM (or in density units, 6.02x10 -7 Å -3 -6.02x10 -5 Å -3 ).
Membrane equilibrated channel structures. The Kv1.2 structure in the open-conductive (O) state was obtained from Treptow and Tarek [29]. The construct was previously acquired via molecular dynamics (MD) simulations of the published x-ray crystal structure [19]. The resting-closed (C) structure of Kv1.2 was obtained from Delemotte et al. [30]. Modeling details and validation can be found in the original papers. It is important to clarify that besides the structure by Delemotte et al., other resting-state models have been proposed for Kv1.2 [31][32][33][34]. Given that these structures were proven to provide with a consensus model for the resting state of the channel [35], we focus our investigation on the former model.
Structures C and O were embedded in the lipid bilayer for Molecular Dynamics (MD) relaxation and subsequent molecular docking of sevoflurane. Specifically, each structure featuring three K + ions (s4s2s0) at the selectivity filter was inserted in a fully hydrated and zwitterionic all atom palmitoyloleylphosphatidylcholine (POPC) phospholipid bilayer. After assembled, each macromolecular system was simulated over a MD simulation spanning~20 ns, at constant temperature (300 K) and pressure (1 atm), neutral pH, and with no applied TM electrostatic potential. The channel structures remained stable in their starting conformations throughout the simulations. Indeed, the RMSD profile of the channel in each simulation converges to a plateau value of approximately 4.0 Å, indicative of structural stability of the constructs [29,30]. Structures sampled in the steady phase of the trajectories were considered in subsequent docking and FEP calculations. Molecular Docking. We used AutoDock Vina [36] to dock sevoflurane against the MDgenerated ensemble of channel structures C and O. Each ensemble included 120 independent channel configurations at least. Docking solutions were resolved with an exhaustiveness parameter of 200, by searching a box volume of 100 x 100 x 100 Å 3 containing the transmembrane domain of the protein receptor. Sevoflurane was allowed to have flexible bonds for all calculations. Clustering of docking solutions was carried out following a maximum neighborhood approach.
Molecular dynamics. All MD simulations were carried out using the program NAMD 2.9 [37] under Periodic Boundary Conditions. Langevin dynamics and Langevin piston methods were applied to keep the temperature (300 K) and the pressure (1 atm) of the system fixed. The equations of motion were integrated using a multiple time-step algorithm [38]. Short-and long-range forces were calculated every 1 and 2 time-steps respectively, with a time step of 2.0 fs. Chemical bonds between hydrogen and heavy atoms were constrained to their equilibrium value. Long-range electrostatic forces were taken into account using the Particle Mesh Ewald (PME) approach [39]. The CHARMM36 force field [40] was applied and water molecules were described by the TIP3P model [41]. All the protein charged amino acids were simulated in their full-ionized state (pH = 7.0). All MD simulations, including FEP and voltage-driven simulations (see next), were performed on local HPC facility at LBTC.

Free-energy perturbation (FEP).
Eq (1) was simplified here for the condition of ligand interactions to multiple independent sites-a condition that appears to be fulfilled at the channel structures featuring sparse binding sites for sevoflurane. Within this scenario, binding constants for structures X � {C, O} were factorized as the product of independent equilibrium constants K X ðn 1 ; . . . ; n s Þ ¼ K X ðn 1 ; 0 2 ; . . . ; 0 s Þ � . . . � K X ð0 1 ; . . . ; 0 sÀ 1 ; n s Þ where, X ðn s ÞÀ n s m� denote respectively the binding constant of n j ligands to each of the j sites at structure X. Accordingly, the excess chemical potential m associated with coupling of the ligand from gas phase to bulk water and W � X ðn j Þ associated with coupling of n j ligands from gas phase to site j under restraints were quantified via FEP. Because computation of m does not depend on the choice of concentration, so long as the same thermodynamic state is used for the solution and gas phases, we estimated the excess potential by considering one sevoflurane molecule embedded into a water box of 60 x 60 x 60 Å 3 . W � X ðn j Þ was computed by taking into consideration the whole ligand-channel-membrane system.
All FEP calculations were performed in NAMD 2.9 [37] by considering the Charmm-based parameters for sevoflurane as devised by Barber et al. [42]. Starting from channel-membrane equilibrated systems containing bound sevoflurane as resolved from docking, forward transformation were carried out by varying the coupling parameter in steps of 0.01. Each transformation then involved a total of 100 windows, each spanning over 31800 steps of simulation. For the purpose of improving statistics, free-energy estimates and associated statistical errors were determined using the simple overlap sampling (SOS) formula [43] based on at least two independent FEP runs.
Specifically for ligand-protein calculations, the free-energy change W and k X 1 j q ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi trA j þ trB j p was determined [46] from the square root of the covariance matrices A j and B j associated respectively to C and O samples at site j. Specifically, A j and B j were computed as symmetric 3 × 3 covariance matrices for centroid positions R j of the ligand at site j X j ¼ hðR j À hR j iÞ:ðR j À hR j iÞ T i and their square roots were solved from the column major eigenvectors {R l , R 2 , R 3 } of the rotation matrix R and the associated eigenvalues {λ l , λ 2 , λ 3 }. Note that overlap is expectedly 1 for identical samplings and 0 for orthogonal configuration spaces. Absolute binding free energy and ensemble averages. An absolute binding free-energy  (2). Position-dependent probability densities. As demonstrated in reference [21], statedependent probabilities ρ X (n 1 ,. . .,n s ) for channel structures X � {C, O} can be mapped into the probability density ρ X (R) of any given ligand i to occupy position R in the system (regardless the position of the remaining N − 1 ligands). Given our original consideration that the reservoir is a homogeneous volume occupied by ligands with position-independent density r, the probability ρ X (R) simplifies to r X ðRÞ ¼ r j X ðRÞ; 8 R 2 dV j r; reservoir ( Binding of the general anesthetic sevoflurane to ion channels for every protein site j = 1,. . .,s. The determination of ρ X (R) thus reduces in practice to knowl-was determined in the form of the "electrical distance" where, F(V) is the local-electrostatic potential of the ion at the central cavity of the open channel. In practice, we solve δ � from two independent 2ns-long simulations at voltages V = 0mV and V = 600mV. For both runs, F(V) was estimated from the electrostatic potential map of the system and subsequently considered to solve δ � for δV = 600mV.
To investigate the conduction properties of Kv1.2 with sevoflurane bound to the main pore, the open channel structure was simulated under depolarized-membrane conditions using a charge-imbalance protocol [50].
Partition function decomposition. In the limit of s independent sites, binding constants can be factorized as the product of independent equilibrium constants then ensuring the associated partition function to be factorized in terms by noting that the voltage-dependent equilibrium constant between channel conformations is given by  S3 Movie. Ion conduction simulation in the presence of sevoflurane bound to channel central cavity. Molecular dynamics simulation with 600mV depolarizing potential was carried out to enhance sampling of conduction events. Kv1.2 channel, sevoflurane and K + ions can be seen in white, blue and yellow, respectively. Water molecules nearby potassium ions and sevoflurane are also shown. Note that the presence of sevoflurane does not hinder ion conduction. (MPG) S1 Table. FEP calculations and equilibrium binding constants for singly-and doubly-occupied sites j of the closed channel structure # . (PDF) S2 Table. FEP calculations and equilibrium binding constants for singly-and doubly-occupied sites j of the open channel structure # . (PDF)