Binding of General Anesthetics 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 misses. We focus on clinically used sevoflurane interaction to anesthetic-sensitive Kv1.2 mammalian channel to resolve if sevoflurane binds the 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 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 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 reassuring as demonstrates that the process of multiple 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, voltagegated and nongated ion channels are bestknown players Franks, 2008;Franks and Honoré, 2004). Primary exemplars are the Cysloop nicotinic acetylcholine and aminobutyric acid class A receptors, the voltagegated sodium and potassium channels, and the tandem pore potassium channels. γ An extensive series of electrophysiological studies corroborate the hypothesis by demonstrating inhibition to potentiation effects 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 clinicallyrelevant anesthetic outcome. For instance, the knockout of the nongated tandem pore potassium channel trek1 produces an animal model (Trek1/) resistant to anesthesia by inhalational anesthetics (Heurteaux et al., 2004). 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 contributing for anesthetic action (Cantor, 1997;FinolUrdaneta et al., 2010;Roth et al., 2008), crystallography and manifold studies involving molecular dynamics support that anesthetics bind ion channels at clinical concentrations (Arcario et al., 2017;Barber et al., 2011Barber et al., , 2014aBrannigan et al., 2010;Jayakar et al., 2013;Kinde et al., 2016;LeBard et al., 2012;Nury et al., 2011;Raju et al., 2013). Binding interactions have been evidenced in anesthetic containing systems of mammalian voltage and ligandgated channels, and bacterial channel analogs as well. Specifically, partitioning of anesthetics in the membrane core allows accessibilityto and bindingto multiple transmembrane (TM) protein sites featuring single or multiple occupancy states a process that might depend further on chemotypes, channeltypes and conformations. Although some progress has been made in one or more aspects of the current view, a firstprinciple demonstration that anesthetics bind ion channels to affect protein equilibrium and function as recorded in experiments still misses. Here, we focus our efforts on the haloether sevoflurane and its molecular interaction to Kv1.2, a mammalian voltagegated potassium channel. Experimental work supports that sevoflurane potentiates the channel in a dosedependent manner (Annika F. Liang et al., 2015). Effects on current tracings include a leftward shift in the conductancevoltage relationship of the channel and an increased maximum conductance. Among all other aspects that might impact channelanesthetic interactions in general, we are specifically interested in determining if sevoflurane binds the wellcharacterized openconductive (O) and restingclosed (C) structures of Kv1.2 (Long et al., 2005;Stock et al., 2013) in a conformationdependent manner to impact protein equilibrium. Very recently, we went through an innovative structurebased study (Stock et al., 2017) of concentrationdependent binding of small ligands to multiple saturable sites in proteins to show that sevoflurane binds the openpore structure of Kv1.2 at the S4S5 linker and the S6Phelix interface a result largely supported by independent photolabeling experiments (Bu et al., 2017;Woll et al., 2017). Here, we aim at extending these previous calculations to investigate sevoflurane interactions with the entire TMdomain of the channel and more importantly, to resolve any conformational dependence for its binding process to channel structures. Accordingly, in the following sections, we first provide the theoretical framework to study binding of sevoflurane to a fixed conformation of the channel under equilibrium conditions. A statedependent strategy is put forward to describe anesthetic binding in terms of occupancy states of the channel that embodies multiple saturable sites and concentration effects. The strategy is then generalized to account for ligand effects on the CO equilibrium, allowing for reconstruction of voltagedependent 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 openconductive state shifting leftward the open probability of the channel at diluted ligand concentrations.

Theory
Anesthetic Binding and Channel Energetics. Consider the voltagegated channel embedded in a large membraneaqueous volume that contains N ligand molecules under dilution. The protein is assumed to remain in a welldefined conformational state X providing with 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 ρ and excess chemical potential μ . 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 j max 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 Κ X (n 1, ... , n s ) for the process O X * (0 1 , ... ,0 s )+n L ⇔O X * (n 1 , ... ,n s ) where, O X * (0 1 , ... ,0 s ) is the empty receptor state with all ligands occupying the bulk. As shown in previous work (Stock et al., 2017), Κ X (n 1 ,... ,n s ) can be evaluated from MDbased free energy perturbation (FEP) calculations in which μ is the solvation free energy of the ligand in the bulk and W X * (n) corresponds to the freeenergy of n sitespecific bound ligands relative to a gas phase state given that ligands i=1,... ,n are restrained to occupy an effective site volume [1] is solved for the thermodynamic limit N ≫n and 1 n 1 !...n s ! corrects the binding constant for equivalent configurations of n j indistinguishable ligands within the site volumes δV j . Within this formulation, knowledge of Κ X (n 1 ,...,n s ) ensures the probability of any occupancy state ρ X (n 1 ,... , n s )=ρ (n 1 +...+n s ) Κ X (n 1 ,... , n s ) to be known in practice from freeenergy calculations (Chipot and Pohorille, 2007). Note in eq. [2], ρ X (n 1 ,...,n s ) depends on the number density or concentration of the ligand at the reservoir thus providing us with 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 consider the very same microscopic system submitted to a Nernst potential induced by nonsymmetrical electrolytes between membrane faces. The capacitive nature of the channelmembrane system ensures the Nernst potential to account for a voltage difference V across the lipid bilayer. Accordingly, by denoting as r P the entire set of Cartesian coordinates of the channel, the free energy of the protein F X (V ) in the particular conformation X ≡ X (r P ) 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 ) (Roux, 2008). From eq. [3], the open probability of the channel then reduces to for the case of a voltagegated channel with two conformational states X≡{ C ,O } connected by the reaction process C ⇔ V O . In terms of chemical freeenergies of the receptor F C (V =0) and F O (V =0) and the corresponding excess freeenergies Δ F C (V ) and Δ F O (V ) , eq. [4] simplifies into the familiar twostate Boltzmann equation 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 ) (Roux, 2008). From eq.
[5], the equilibrium constant between protein states C and O then writes as with Κ(0)=e −βV m Δ Q determining their equilibrium at 0 mV. In eq. [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 voltagedependent excess free energies are proportional to the applied voltage and associated protein charges.
By combining eq. [2 and 5] through a generalized thermodynamiccycle analysis dealing with all possible states of the ligandfree and ligandbound 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 of the protein conformations. Eq.
[8] simplifies into the twostate Boltzmann equation embodying now the freeenergy contributions arising from ligand binding. Note that eq.
[8] is achieved by rewriting the state probability densities in terms of the reference state ρ C (0 1 ,... , 0 s ,V ) .
In the following, we consider eq. [1, 5 and 9] to investigate the molecular binding of sevoflurane to open and closed structures of Kv1.2, and its functional impact on the channel energetics.

Results and Discussion
Binding of Anesthetics to Multiple Channel Sites. We applied largescale 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, S6Phelix interface and at the extracellular face, next to the selectivity filter. Further docking sites were resolved within the voltagesensor, at the S4Pore interface and at the centralcavity of the channel. 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 is up to 2×3 21 occupancy states of the channel structures that might contribute for sevoflurane binding and functional effects. To evaluate this quantitatively, we performed an extensive series of FEP calculations to estimate the persite binding affinity for one and twobound ligands against the channel structures (Fig. S1, TableS1 and S2). Binding constants for the individual sites are heterogeneous and take place under 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 S6Phelix 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 Κ X (n 1 ,...,n s ) for every occupancy state of the channel were then reconstructed from the persite affinities to determine state probabilities via eq. [2]. Here, estimates of Κ X (n 1 ,...,n s ) were determined for the condition of independent binding sites as minimum sitetosite distances of ~15 Å demonstrated their nonoverlap distributions in each of the channel structures. At low 1mM concentration, ρ X (n 1 ,... ,n s ) are largely dominated by the probability of the empty state ρ X (0 1 ,... ,0 s ) implying only a small fraction of bound states with nonnegligible occurrences (Fig. S2). 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 ) at the individual sites ( Fig. 2). At higher concentrations, there is a clear shift of ρ X (n 1 ,...,n s ) towards states of the channel that enhances significantly 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 S6Phelix interface. In contrast, ρ X (n j ) for sites within the voltagesensor, at the central cavity and nearby the extracellular face of the channel remains negligible over all concentrations. For completeness, note in TableS1 that equilibrium constants for doublyoccupied sites are comparable to or even higher than estimates for onebound 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 positiondependent density ρ X j (R) of sevoflurane in each binding site j (cf. eq. 18 And 20). As shown in Fig. 3 and supplementary Movies S1 and S2, the density of sevoflurane makes sense of the results by showing the concentration dependent population of bound ligands. Projection of ρ X j (R) along the transmembrane direction z of the system, ρ X j (z ) , stresses further the results (cf. eq. 20). Note from ρ X j (R) that sevoflurane binds channel structures in a concentration dependent manner, binding preferentially the S4S5 linker and the interfaces S4Pore and S6Phelix over a range of concentrations.
So far, our calculations demonstrate that sevoflurane binds Kv1.2 structures over a range of concentrations, preferentially at the linker S4S5 and at the segment interfaces S4Pore and S6Phelix. From a physicalchemical point of view, spots at these channel regions are primarily dehydrated lipidaccessible amphiphilic pockets providing with favorable interaction sites for the polar lipophilic sevoflurane molecule (Fig. S3). It is worth mentioning that these findings recapitulate recent photolabeling experiments demonstrating that photoactive analogs of sevoflurane do interact at the S4S5 linker and at the S6Phelix interface of the open conductive Kv1.2 channel (Bu et al., 2017;Woll et al., 2017). In detail, Leu317 and Thr384 were found to be protected from photoactive analogs, with the former being more protected though. As shown in Fig. S4, atomic distances of bound sevoflurane to these aminoacid side chains are found here to be respectively 7. 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 also recapitulate the stronger protection of Leu317 in the sense that relative to sites at S6Phelix, the affinity of sevoflurane is found to be higher at the S4S5 linker given its stable occupancy by one or two ligands. The stable occupancy of the linker by one or two ligands as computed here, is consistent with recent floodingMD simulations of the homologous sodium channel NaChBac (Annika F Barber et al., 2014a) and more importantly, with previous Ala/Valscanning mutagenesis showing a significant impact of S4S5 mutations on the effect of general anesthetics on family members of K + channels (Barber et al., 2011). In special, a single residue (Gly329) at a critical pivot point between the S4S5 linker and the S5 segment underlines potentiation of Kv1.2 by sevoflurane (Liang et al., 2015). Sevoflurane is close to that amino acid when bound at the S4S5 linker. In contrast to the aforesaid spots, sites within the voltagesensor, at the main pore and nearby the extracellular face of the Kv1.2 structures are primarily hydrated lipidinaccessible amphiphilic pockets that weaken sevoflurane interaction as reflected in the state and spacedependent densities shown in Fig. 2 and 3. The binding probabilities at these sites thus support that the nonnegligible 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 (Annika F. Liang et al., 2015), we read the negligible or absent density of the ligand at the centralcavity of the channel as a selfconsistent 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 the ionic flux conditions in which potentiation takes place. Accordingly, as shown in Fig. S5, the binding affinity of a potassium ion at the central cavity overcomes that of sevoflurane due its binding and excess freeenergies under applied voltages. Once bound, the ion destabilizes sevoflurane interactions and the molecule is not expected to bind the channel cavity at low concentration. As also shown in Fig. S5 supplementary Movie S3, 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 (Annika F Barber et al., 2014a) due major structural differences between Na + and K + channels. Specifically, the pore of potassium channels lacks lipidaccessible openfenestrations 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 ρ X j (R) reveals for most sites, an obvious differential affinity of sevoflurane across Kv1.2 structures ( Fig. 2 and 3). The overall consequence for sevoflurane binding is then clear: the average number of bound ligands to the openconductive channel exceeds systematically that number for the restingclosed channel over the entire concentration range. There is therefore a remarkable conformational dependence for the anesthetic interaction, with sevoflurane binding preferentially the openconductive structure.
Implications for Kv1.2 energetics were then investigated by quantifying modifications of the open probability ρ O (V ) of the channel 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 eq. [5 and 9] shows that sevoflurane shifts leftward the open probability of Kv1.2 in a concentrationdependent manner voltage shifts amount from 1.0 mV to 30.0 mV with concentration increase of the ligand in solution. 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 nontrivial interplay of conformationdependent modes of action involving distinct sites (cf. eq. [23 and 24]). In detail, binding of sevoflurane at the low affinity sites within the voltagesensor, central cavity and next the extracellular face of the channel are mostly conformation independent and do not impact open probability (ratio ≈ 1). On the other hand, conformationdependent 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 S6Phelix and reflected in rightward shifts of ρ O (V ) (ratio > 1). The overall conformationdependent binding process is therefore encoded differentially across distinct channel regions. [18], this 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 from eq. [19] by coarsegraining over state probabilities in Fig. S2 whereas, ρ X ( R|n j ) was calculated from the centroid distributions of docking solutions shown in Fig. 1B and 1C. (B) Projection of ρ X j ( R) along the transmembrane direction z of the system, ρ X j ( z) . Projections were determined as prescribed in eq. [20].
Potentiation of Kv1.2 by sevoflurane has been attributed to stabilization of the openconductive state of the channel (Liang et al., 2015). Given the critical role of S4 and S4S5 linker on the gating mechanism of the channel (Long et al., 2005), it is likely that sevoflurane interactions with these segments as found here are at the origins of the experimentally measured voltagedependent component of anesthetic action. While restricted to sevoflurane interactions with the restingclosed and openconductive structures, the presented twostate binding model only embodies left or rightward shifts in the open probability of the channel and therefore, it cannot clarify any molecular process accounting for maximum conductance increase as experimentally recorded and shown in Fig.  4A. As supported by a recent kinetic modeling study (Annika F. , generalization of eq. [9] to include a third non conducting open state yet structurally unknown is needed to account for such conductance effects. We then speculate that binding of sevoflurane at the S4Pore and S6Phelix interfaces might interfere allosterically with the pore domain operation thus affecting channel's maximum conductance. A working hypothesis also raised in the context of anesthetic action on bacterial sodium channels (Barber et al., 2014a;Raju et al., 2013), assume indeed that nonconducting states of the selectivity filter are implicated.
Corroboration of a such assumption from a molecular perspective is however not trivial and will necessarily involve further structural studies to demonstrate how ligand binding might impact nonconducting open states of the channel to affect maximum conductance.

Concluding Remarks
Here, we carried out extensive structurebased calculations to study conformationdependent binding of sevoflurane to multiple saturable sites of Kv1.2 structures X ≡{C , O} 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 j max =2 . Our study relied on the assumption that molecular docking 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 (Deng et al., 2015). In this regard, the combination of extensive docking calculations against an ensemble of equilibrium receptor structures to handle protein flexibility and FEP calculations based on fine forcefields to accurately estimate solvation energies are critical technical aspects of the applied methodology to minimize such drawbacks (Deng et al., 2015). Given the same limitations of the scoring function, it is also not guaranteed that 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 allatom floodingMD simulations (Arcario et al., 2017;Barber et al., 2014a;Brannigan et al., 2010;LeBard et al., 2012;Raju et al., 2013), might also be an important refinement in that direction (manuscript in preparation). Despite these sampling improvements that may eventually be obtained, it is worth mentioning that the configuration space in FEP calculations overlap between channel structures at individual sites, meaning that sampling and binding affinities were evenly resolved between states ( Fig.  S1 and 4B). Besides that, most of the identified binding sites are located nearby flexible protein regions for which the rootmean 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 binding sevoflurane at these gatingimplicated 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 molecularlevel studies of anesthetic action thereby justifying our choice. In detail, the atomistic structures complain most of the available experimental data characterizing closed and open conformations of the channel in the native membrane environment (Stock et al., 2013). Previous findings support further that sevoflurane binds Kv1.2 to shift leftward its voltagedependence and to increase its maximum conductance in a dosedependent manner (Liang et al., 2015). 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 openconductive structure to impact the CO energetics in a dosedependent 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 (Cantor, 1997;FinolUrdaneta et al., 2010;Roth et al., 2008), the result is reassuring by showing that in principle the isolated process of sevoflurane binding to Kv1.2 accounts for openprobability shifts as recorded in experiments. Within this scenario, the calculations reveal unexpectedly, contrasting persite 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 S6Phelix interface. By showing this nontrivial interplay of conformationdependent modes of action involving distinct binding sites, the result is particularly insightful and should guide us to design novel sitespecific 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 clinicallyrelevant anesthetic outcome, our study adds support to the directsite 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 maintext Fig. 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.

Computational Methods
A procedure was designed to solve the molecular binding of sevoflurane to the openconductive (O) and restingclosed (C) structures of Kv1.2 for saturation conditions up to n j max =2 . For both channel structures, the procedure consisted of (i) an extensive production of docking solutions for the ligandreceptor interaction, (ii) clustering of docking solutions into binding sites along the receptor structure and (iii) estimation of binding affinities using the freeenergy perturbation (FEP) method. First completion of steps (i) through (iii) solved the ligand channel interaction for singlyoccupied binding sites. Double occupancy of the receptor sites was investigated by inputing 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 MDgenerated ensemble of membraneequilibrated 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 were treated as binding site regions in step (iii) calculations. FEP calculations were carried out by taking into consideration the whole ligandchannelmembrane system.
Following this procedure, binding constants Κ 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 statedependent probability distributions via eq. [2]. Here, affinity constants were solved for the condition of independent binding sites. Ligandfree and ligandbound open probability curves were respectively computed from eq. [5 and 9] by taking into consideration previously determined experimental values of V m =−15.1 mV and Δ Q=3.9 e o for Kv1.2 (Liang et al., 2015). Estimates were determined for sevoflurane concentrations in the range of 1mM 100mM (or in density units, 6.02x10 7 Å 3 -6.02x10 5 Å 3 ). A detailed description of the calculations is provided bellow.
Membrane Equilibrated Channel Structures. The Kv1.2 structure in the openconductive (O) state was obtained from Treptow and Tarek (Treptow and Tarek, 2006). The construct was previously acquired via molecular dynamics (MD) simulations of the published xray crystal structure (Long et al., 2005). The restingclosed (C) structure of Kv1.2 was obtained from Delemotte et al. (Delemotte et al., 2011). Modeling details and validation can be found in the original papers. 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 an 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. The root meansquare deviation (rmsd) values for the channel structures range from 1.5 to 3.5 Å, which agrees with the structural drift quantified in previous simulation studies (Delemotte et al., 2011;Treptow and Tarek, 2006). Molecular Docking. We used AutoDock Vina (Trott and Olson, 2010) 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 (Phillips et al., 2005) 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 timestep algorithm (Izaguirre et al., 1999). Short and longrange forces were calculated every 1 and 2 timesteps respectively, with a time step of 2.0 fs. Chemical bonds between hydrogen and heavy atoms were constrained to their equilibrium value. Longrange electrostatic forces were taken into account using the Particle Mesh Ewald (PME) approach (Darden et al., 1993). The CHARMM36 force field (Huang and MacKerell, 2013) was applied and water molecules were described by the TIP3P model (Jorgensen et al., 1983). All the protein charged amino acids were simulated in their fullionized state (pH=7.0). All MD simulations, including FEP and voltagedriven simulations (see next), were performed on local HPC facility at LBTC amounting to a total run time of ~2.0 μs. FreeEnergy 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 Κ X (n 1 , ... ,n s )=Κ X (n 1 ,0 2 ... ,0 s )×...×Κ X (0 1 ,... , 0 s−1 , n s ) where, ...
denote respectively the binding constant of n j ligands to each of the j sites at structure X . Accordingly, the excess chemical potential μ 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 μ does not depend upon 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 considering the whole ligandchannelmembrane system. All FEP calculations were performed in NAMD 2.9 (Phillips et al., 2005) by considering the Charmmbased parameters for sevoflurane as devised by Barber et al. (Barber et al., 2014b). Starting from channelmembrane 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, freeenergy estimates and associated statistical errors were determined using the simple overlap sampling (SOS) formula (Lu et al., 2004) based on at least two independent FEP runs. Specifically for ligandprotein calculations, the freeenergy change W X * (1 j ) for singlyoccupied sites j was computed as a FEP process that involves ligand coupling to a vacant site. Differently, for doublyoccupied sites, W X * (2 j ) was computed as a twostep FEP process involving ligand coupling to a vacant site W X * (1 j ) followed by binding of a second ligand at the preoccupied site is a state function, the stepwise approach is equivalent to a singlestep process involving simultaneous coupling of two ligands to the protein site that is, W X * (2 j )=W X * (1 j )+W X * (2 j |1 j ) . The colvars module  in NAMD 2.9 was used to apply the harmonic restraint potentials when computing these quantities. The value of W X * (n j ) depends on the parameters of the restraint potential adopted in the FEP calculation ie., the reference positions of the ligands in the bound state { R X * (1 j ),... ,R X * (n j )} and the magnitude of force constants { k X (1 j ),..., k X (n j )} . By minimizing the contribution of the restraint potential to the binding freeenergy W X * (n j ) , Roux and coworkers (Roux et al., 1996) devised optimum choices for the parameters and in which, ⟨ R X (1 j )⟩ ,... ,⟨ R X (n j )⟩ and ⟨[δ R X (1 j )] 2 ⟩ ,... , ⟨[ δ R X (n j )] 2 ⟩ are respectively the equilibrium average positions for each of the n j bound ligands at site j and their corresponding meansquare fluctuations when interacting to structure X . Here, these parameters were estimated from the space of docking solutions and the resulting force constants, in the range of 0.03 to 1.35 kcal/mol/Å 2 , were considered for computations of the bound state.
The equilibrium binding constant (eq. [10 and 11]) and following results are derived in the limit of a homogeneous diluted reservoir occupied by ligands at constant density ρ and excess chemical potential μ . Given that, we treated the system reservoir as a homogeneous aqueous solution despite its intrinsic inhomogeneity provided by the solvated lipid bilayer. An excess chemical potential of 0.1 kcal.mol 1 was estimated here as the reservoir potential for sevoflurane in bulk water.
Sampling Overlap. Here, a persite measure of sampling overlap o( A j , B j ) between FEP configurations in structures C and O was determined (Hess, 2002) 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 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 freeenergy Δ G X o (n 1 ,... ,n s ) (Gilson et al., 1997) associated with state O X * (n 1 ' ,..., n s ' ) can be defined as where it is understood that this refers to the free energy of binding n ligands to the protein structure X≡{ C ⟨ A X ⟩ (n 1 ' , ..., n s ' ) ρ X (n 1 ' , ...,n s ' ) as the ensemble average of any thermodynamic property of the system A X (n 1 ' ,... ,n s ' ) for state O X * (n 1 ' ,... , n s ' ) can be known from eq. [16].
PositionDependent Probability Densities. As demonstrated in reference (Stock et al., 2017), 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 positionindependent density ρ , the probability ρ X (R) simplifies to for every protein site j=1,... ,s . The determination of ρ X (R) thus reduces in practice to knowledge of the persite density ρ X j (R) where, ρ X (R|n j ) is the local density at site j when occupied exactly by n j molecules and ρ X (n j ) is the probability for this occupancy state. In eq.
[18], ρ X (R|n j ) describes the local equilibrium density of the ligand, conditional to a specific number of bound molecules that satisfies ∫ δV j d Rρ X (R|n j )=n j . In contrast, ρ X (n j )= ∑ n 1 ' ,... , n s ' δ n j ' ,n j ρ X (n 1 ' ,...,n s ' ) denotes the marginal probability of site j to be occupied by n j ligands regardless the occupancy of the other sites. Eq.
[18] establishes a formal relation between spacedependent and statedependent densities of the system. At a fine level, this relation involves the set of equilibrium constants Κ X (n 1 ,...,n s ) satisfying ρ X (n j ) . From eq.
[18], spatial projections of ρ X (R) along the transmembrane z direction of the system can be achieved as where, A (z)=Δ x Δ y is the total area of the membraneaqueous region along the Cartesian x and y directions.
CoarseGraining Over States. Consider any macrostate O X * (n) of the system mapping an ensemble of accessible states O X * (n 1 , ... ,n s ) in which n ligands bind the receptor regardless their specific distributions over the binding sites. Because O X * (n) is degenerate, the probability density of the macrostate ρ X (n)= ∑ n 1 ' ,... ,n s ' δ n ' , n ρ X (n 1 ' ,... ,n s ' ) can be determined by coarsegraining over the receptor states O X * (n 1 , ... ,n s ) featuring exactly n=n 1 +...+n s bound ligands. Here, the Kronecker delta function δ n ' ,n ensures summation over states accessible to O X * (n) only. ]. An inwater excess potential of 69.52 kcal.mol 1 was estimated for potassium.

Binding of Potassium and
Specifically for K + , a total binding freeenergy was obtained by summing up its absolute binding free energy with its charge ( q ) excess free energy ( q ϕV ) under an applied external voltage V (Dong et al., 2013;Souza et al., 2014). The voltage coupling ϕ was determined in the form of the "electrical distance" where, Φ(V ) is the localelectrostatic potential of the ion at the central cavity of the open channel. In practice, we applied the charge imbalance protocol (see next) to solve δ ϵ from two independent 2nslong simulations at voltages V =0 mV and V =600 mV . For both runs, Φ(V ) was estimated from the electrostatic potential map of the system and subsequently applied into eq.
To investigate the conduction properties of Kv1.2 with bound sevoflurane at the main pore, the open channel structure was simulated under depolarizedmembrane conditions using a chargeimbalance protocol (Delemotte et al., 2008).
Partition Function Decomposition. In the limit of s independent sites, binding constants can be factorized as the product of independent equilibrium constants eq.
[23] is useful to estimate the persite contributions impacting the opening probability of the channel as defined in eq. [9]. For any given site j , ratio values