despite well over 150 
 years of clinical use , a detailed molecular 
 mechanism of anesthesia remains unresolved .
the once prevailing hypothesis 
 suggested that anesthetics worked through a nonspecific , membrane 
 disruption mechanism , known as the meyer 
recent studies have shown that most , if not all , anesthetics affect a family of 
 ion channels in the nervous system known as the cys - loop family of 
 pentameric ligand - gated ion channels ( plgic ) .
while other potential 
 targets of anesthetics have been also found , e.g. , the hcn family 
 of channels , the two - pore k leak channels , and voltage - gated 
 cation channels , cys - loop receptors remain the 
 best characterized and most studied of membrane channels in the context 
 of anesthetic action .
these receptors act in response to the release 
 of neurotransmitters from the presynaptic terminal and are targets 
 for multiple pharmacological agents , including alcohol , barbiturates , 
 and benzodiazepines , in addition to anesthetics .  the structural changes associated with channel modulation , 
 however ,
are still poorly understood because of the challenge of resolving 
 high - resolution structures of these receptors .
recently , crystal structures 
 of glic ( g. violaceus ligand - gated ion channel ) and elic ( e. chrysanthemi ligand - gated ion channel ) , which are bacterial homologues of nervous system plgics , have been 
 solved in both the presence and absence of anesthetics and provide 
 a tremendous opportunity to resolve the molecular mechanism behind 
 plgic function and modulation .
several hurdles remain to understanding 
 anesthetic action , however , 
 including determining the energetics and kinetics of anesthetic binding 
 to plgics , elucidating how anesthetic binding disturbs the normal 
 dynamics of these channels , and uncovering the anesthetic binding 
 site ( or sites ) in ligand - gated ion channels .
multiple methods have 
 been used to try to elucidate the latter , including crystallography , mutagenesis , photolabeling , and molecular dynamics ( md ) ; however , multiple binding 
 sites have been reported and there is no consensus as to which site 
 modulates ion channel function .
while static structures provide insight 
 into possible binding sites of anesthetics , the modulation of channel 
 function is an inherently dynamic process that requires a dynamic 
 description of the anesthetic / plgic complex .
md simulation is a technique 
 that can probe the effects of bound anesthetics on ion channel dynamics with a high degree of spatial and temporal resolutions .
moreover , 
 the solution of atomic - resolution ( 2.93.1  ) glic structures 
 with bound anesthetics positions md simulations to be an integral 
 tool in understanding how these molecules affect ion channel dynamics 
 in atomic detail .
modern anesthetics comprise two main classes based on their delivery 
 method : inhaled ( volatile ) anesthetics and intravenous anesthetics . 

most inhaled anesthetics are heavily halogenated ( e.g. , halothane 
 and desflurane ) and can be difficult to model because of the unique 
 chemical properties of large halogens . compared to volatile anesthetics , 
 intravenous anesthetics are generally larger ( e.g. , propofol and lidocaine ) 
 and can comprise complex multi - ring structures ( e.g. , ketamine and 
 etomidate ) that require extensive computational resources to model 
 accurately , making parametrization difficult and time - consuming . because 
 of the availability of parameters for myriad chemical groups in standard 
 force fields , one common solution is to transfer 
 parameters from similar molecules in order to create an amalgamated 
 parameter set for each individual molecule . while this is generally 
 acceptable for bonds and angles , atomic charges and torsions are more 
 complex and not as portable
; moreover , standardized parameters for 
 halogenated ethers are nonexistent and therefore not amenable to this 
 technique .
models of halothane have been published previously 
 but were not designed for compatibility with biomolecular force fields , 
 such as charmm .
additionally , other models of modern anesthetics have 
 been published . however , with the exception of the isoflurane model developed by 
 hnin et al . , for which adequate 
 comparison to experimental data has been performed , the other models 
 lack extensive and rigorous testing against experimentally determined 
 physicochemical properties . in order to study the effect of anesthetics 
 on ion channels ,
a set of standardized parameters compatible with 
 biomolecular force fields that are validated against experimentally 
 verifiable physicochemical properties is necessary . in this 
 study
, we present a novel set of parameters for the ubiquitous 
 modern anesthetics desflurane , isoflurane , sevoflurane , and propofol 
 ( figure 1 ) , which is compatible with the widely 
 used charmm force field for biomolecules .
these four species were chosen because they represent four of the 
 most commonly used clinical general anesthetics that these bacterial 
 channels are also sensitive to .
the parameters 
 generated for each anesthetic were compared against five different 
 experimentally measured properties , namely , dipole moment , density , 
 heat of vaporization , free energy of solvation in water , and free 
 energy of solvation in oil , to assess the accuracy of the parameter 
 set developed . in all four cases ,
the atomic models reproduce experimental 
 physicochemical properties , signifying the robustness of the model 
 and reflecting their ability to reproduce the molecular behavior of 
 anesthetics in multiple environments .
moreover , we calculated 
 the free energy of partitioning of each 
 anesthetic into a 1-palmitoyl-2-oleyl - sn - glycero-3-phosphatidylcholine 
 ( popc ) bilayer to investigate how these molecules interact with the 
 membrane . in all cases , the anesthetics preferred the glycerol backbone 
 region of the membrane over either the membrane core or aqueous solution 
 by 45 kcal / mol .
we measured how the partitioned anesthetics 
 affected membrane structure by computing the atomic density profiles , 
 area per lipid , deuterium order parameters , dipole potential , and 
 lateral stress profile .
no difference was observed in these measurements 
 between a popc - only membrane and a popc membrane with partitioned 
 anesthetics .
the results from the physicochemical studies presented 
 suggest that , in physiological context , anesthetics partition to the 
 membrane , without disturbing the structure of the membrane itself , 
 prior to interacting with plgics .
the recent 
 development of the force field toolkit ( fftk ) by our group enables rapid and robust parametrization 
 of small molecules for application to biomolecular simulations employing 
 charmm - compatible force fields .
in addition to parameter optimization , 
 fftk provides several tools for assessing parameter quality throughout 
 the course of parameter development .
these tools , as well as simulation - based 
 calculations , were utilized extensively to develop optimal parameters 
 for a set of anesthetics as assessed by comparison to experimental 
 data . for a more detailed description of the parametrization process , 

two levels of theory were used to properly 
 model the different characteristics of the heavily halogenated inhaled 
 anesthetics and the large intravenous anesthetics .
ab initio calculations on desflurane , isoflurane , and sevoflurane 
 ( inhaled anesthetics ) utilized the hybrid b3lyp density - functional 
 theory with the 6 - 311+g(2d , p ) basis set for geometry 
 optimization and to calculate bonded , electrostatic , 
 and torsional interactions .
b3lyp is known to calculate accurate geometries 
 and hydrogen - bonded complexes for small molecules containing first 
 row elements and has been used to parametrize 
 halogenated molecules previously . for propofol , second - order mller 
plesset 
 theory ( mp2 ) and the 6 - 31g(d ) basis set 
 were used for geometry optimization and to calculate bonded , electrostatic , 
 and torsional terms .
parametrization of the anesthetics follows 
 the fftk workflow , which is based on 
 charmm general force field parametrization principles , so as to be consistent with the charmm force 
 field .
fftk does not currently support the optimization of lennard - jones 
 ( lj ) parameters de novo .
therefore , all lj parameters 
 for the models presented were adopted from chemically similar atoms 
 in the charmm force field . with the exception of isoflurane , in which 
 the well potential term (  ) was modified ,
once the geometry was optimized , partial charge assignments 
 were made based on minimum interaction energies and distances between 
 water and all accessible atoms ( e.g , an sp carbon would 
 not be accessible to water ) derived from quantum mechanical ( qm ) calculations . 

the qm dipole moment vector was used as an additional constraint in 
 optimizing charge distribution , with equal weighting applied to the 
 qm dipole and qm water interaction distance terms . in order to reproduce 
 condensed phase properties , parameters for md simulations should overestimate 
 the gas - phase dipole moment of the molecule by as much as 20% and an overestimation of the qm dipole by 1020% 

if the dipole moment was outside of this range , another iteration 
 of charge optimization was performed followed by another dipole moment 
 measurement .
the dipole moment was calculated according 
 to1where the bars denote the 
 vector length , n denotes the number of atoms in the 
 molecule , qi denotes 
 the charge of the ith atom , and ri denotes the 
 position vector of the ith atom .
scans of each torsion 
 were made at 3 steps , and the dihedral parameters were fit to 
 the resulting potential energy surface ( pes ) . in order to validate the quality of the parameters generated using 
 this scheme ,
md simulations were performed to calculate bulk thermodynamic 
 properties of each anesthetic , which were subsequently compared to 
 experimental data ( described in further detail below ) . following the optimization of bond , 
 angle , and dihedral terms against the qm data ,
the bulk properties 
 of density and enthalpy of vaporization were calculated for each anesthetic 
 and compared against the corresponding experimentally determined values 
 as an initial assessment of the accuracy of the model .
a temperature 
 of 298  k was utilized throughout the simulations because the 
 anesthetics are liquid at this temperature , and the experimental data 
 were measured at this temperature .
each simulated system comprised 
 216 copies of the anesthetic positioned on a 6  
 6  6 cubic grid each with a random initial orientation 
 about its center of mass .
the system was first minimized for 
 10  000 steps , then equilibrated for 100  ps , and finally 
 simulated in an npt ensemble for 2.5 ns .
the enthalpy of vaporization ( hvap ) was measured according to the following equation:2where 
 uliq is the potential 
 energy of the system in the liquid state , p is the 
 pressure of the system , vliq 
 is the volume of the liquid system , nmol is the number of molecules in the system , 
 and ugas is the potential 
 energy of the gaseous system .
the brackets  ...  denote 
 a time average over the length of the production simulation .
as indicated 
 in eq 2 , the pvliq term is negligible and is ignored in this calculation . assuming the intramolecular 
 energy is the same in gas and condensed phases ( i.e. , ugasintra = uliqintra ) , the intermolecular interactions 
 energies ( i.e. , electrostatic and van der waals interactions ) can 
 be computed from a single bulk phase simulation .
this is done according to3where 
 uliqinter 
 is the average intermolecular energy of the liquid , uliqsys is the average potential 
 energy for the entire system ( including inter- and intramolecular 
 interactions ) , and uliqintrai is the average intramolecular energy of the ith molecule .
the brackets  ...  denote a time average 
 over the length of the production simulation . for the sake of completeness
, 
 we also calculated the heat of vaporization as detailed in the charmm 
 general force field , using one 100 ns 
 gas phase calculation to calculate the intramolecular energy of the 
 molecule .
that the parameter set was not able to reproduce experimental 
 values for density and hvap within 
 the determined error ranges ( for example , in the case of isoflurane ) , 
 then the process of optimization was restarted beginning with charge 
 distribution optimization .
only when the parameter set was able to 
 reproduce experimental bulk properties was the parameter set tested 
 against the more computationally intensive and rigorous free energies 
 of solvation .
as discussed , the 
 calculated dipole 
 is expected to be 1020% higher than experimental gas phase 
 dipole moment to replicate condensed phase properties . because dipole 
 is calculated using only the equilibrium geometry ,
calculated values are presented 
 as the mean  standard deviation , where standard deviation was 
 calculated across five replicates .
values in parentheses are the heat 
 of vaporization calculated utilizing the method employed in the charmm 
 general force field , in which a single 
 but long gas phase trajectory for a copy of the molecular species 
 is presented for comparison .
because of the high boiling point 
 of propofol ( 529 k ) , the hvap varies 
 widely ; however , the calculated value is well within the experimental 
 range .
taking advantage of readily 
 available experimental data on the free energy of solvation in both 
 water ( gsolvh2o ) and olive oil ( gsolvoil ) for the anesthetics studied , we computed these values for each 
 of the anesthetics as the final criterion for determining the quality 
 of the parameters .
as olive oil is difficult to model because of its 
 unknown and variable composition , dodecane has been used as a substitute 
 in these types of calculations . in each 
 calculation , the anesthetic was placed at the origin of a pre - equilibrated 
 solvent box measuring 30  on each side containing either water 
 or dodecane .
alchemical free energy perturbation ( fep ) ( successfully 
 implemented for anesthetics previously ) was used 
 to calculate both gsolvh2o and gsolvoil for 
 each anesthetic by incrementally decoupling the molecule from the 
 surrounding solvent environment followed by incrementally recoupling 
 to the surrounding solvent environment .
each transformation proceeded 
 over 25 windows run both forward ( solvation ) and backward ( desolvation ) , 
 a total of 50 windows , to minimize the hysteresis .
each alchemical 
 window was equilibrated for 1 ps and simulated for 0.5 ns of data 
 collection , requiring a total of 25 ns for each decouple / recouple 
 cycle . to avoid the  end - point catastrophe
 , the zacharias soft - core potential was used in addition to shifting the lennard - jones 
 potentials by 5.5   and slowly decoupling the electrostatics 
 and van der waals interactions of the anesthetic from solution .
the 
 transformation was repeated 10 times for each molecule and the free 
 energy change calculated using the bennett acceptance ratio ( bar ) .
thus , the free energy calculations required 
 a total of 250 ns of fep simulation per anesthetic in each environment 
 ( a total of 2  s of simulation for all fep transformations ) . umbrella sampling ( us ) was used to 
 calculate the free energy profiles for partitioning of anesthetics 
 into a lipid bilayer based on similar protocols used for small molecules .
a total of 38 umbrellas with a 1  distance between biasing 
 potentials were defined along the membrane normal ranging from z = 0 ( center of membrane ) to z = 38 ( bulk 
 aqueous solution ) .
the center of mass of each anesthetic was restrained 
 to the center of each umbrella by a harmonic potential using a force 
 constant of 7.17 kcal mol  .
each umbrella was minimized for 5000 steps and simulated for 10 
 ns , recording the position of the molecule every 0.2 ps .
the free 
 energy profile was reconstructed with wham , utilizing the last 9.5 ns of the trajectory for the analysis .
therefore , 
 calculation of the profile for inserting a single anesthetic into 
 the membrane required a total of 380 ns of simulation ( totaling 1.52 
 s for all four anesthetics ) . in order to increase computational 
 efficiency and
gain better statistics for each anesthetic , the simulation 
 system contained two copies of each anesthetic offset in the z - direction by 38  to prevent unwanted interactions 
 between the two molecules .
once the pmf was calculated , the profile 
 was used to calculate the popc / water partition coefficient using the 
 following equation:4where w0 is the 
 width of the membrane taken to be the point along the membrane normal 
 ( z - axis ) where the choline density decays to zero 
 ( z = 25  ) . in order to better 
 understand how anesthetics affect membranes and membrane structure , 
 we simulated the interaction of a high concentration of anesthetic 
 with a biological membrane .
the popc membrane , measuring 96  
 96  in the xy - plane and constructed 
 using the membrane builder plug - in of vmd , was solvated and ionized 
 to 150 mm nacl using the solvate and autoionize plug - ins of vmd , respectively .
this gave the system initial dimensions of 96  96   120   with 216 lipids ( 108 lipids per leaflet ) .
the resulting 
 system was 
 minimized for 5000 steps and equilibrated for 3 ns in an npt ensemble ( p = 1.0 atm , t = 298  k ) , 
 giving final dimensions of 95  86  95  , and was used as the starting point for the membrane - only ( i.e. , 
 anesthetic and no protein ) simulations . by use of the equilibrated 
 system ,
anesthetic molecules , including desflurane , isoflurane , sevoflurane , 
 and propofol , were placed randomly in the aqueous phase at an anesthetic / lipid 
 ratio of 1:5 and each of the resulting systems was simulated for an 
 additional 100 ns .
all simulations 
 were performed 
 using namd2 together with the set of 
 force fields parameters developed herein for each anesthetic ( supporting information ) and the charmm36 parameters for the tip3p water model and popc .
topology and parameters for dodecane were generated by 
 analogy from hexane in the charmm36 parameter set .
the target pressure 
 was set to 1.0 atm and the temperature of the system was set to 298 
 k ( 25 c ) in order to compare directly to experimental physicochemical 
 properties .
hoover 
 langevin piston method , and constant temperature 
 was maintained using a langevin damping coefficient (  ) of 1 
 ps .
nonbonded interactions were cut off after 
 12  with a smoothing function applied after 10  .
long - range 
 electrostatic interactions were treated using the particle mesh ewald 
 ( pme ) method with a grid density greater 
 than 1  .
a time step of 2 fs was used with 
 bonded and nonbonded forces calculated at every time step , while pme 
 calculations were performed at every other time step .
the recent availability 
 of multiple atomic - resolution structures of glic , elic , and glucl ( a eukaryotic cys - loop receptor 
 from c. elegans ) has 
 enabled the detailed , atomic - level study of anesthetic - induced changes 
 in plgics by md simulations .
a set of widely available parameters 
 exhaustively tested against experimental physicochemical properties 
 is lacking , however , for most modern anesthetics .
therefore , we have 
 rigorously developed a set of parameters that are compatible with 
 the widely used charmm force field for four ubiquitous anesthetics : 
 desflurane , isoflurane , sevoflurane , and propofol .
these four anesthetic 
 species were chosen because they are known to inhibit the bacterial 
 channel glic and , therefore , can be widely 
 utilized in studying anesthetic action in addition to the fact that 
 these four molecules represent the most widely used general anesthetics . 

while desflurane and sevoflurane have not previously been parametrized 
 for the charmm force field , models for isoflurane and propofol have 
 recently been published , along with a desflurane 
 model for gromos .
a summary of recent 
 md simulations employing anesthestics together with the force fields 
 utilized is compiled in table s1 of supporting 
 information .
all four anesthetics were parametrized de novo 
 in this study ( see supporting information for charmm format parameter and topology files for all four anesthetics ) .
the lennard - jones ( lj ) terms for all atoms were taken by analogy 
 from the charmm force field and , in all cases except isoflurane , were 
 left unchanged in the models presented here .
the equilibrium geometry 
 was first calculated for each anesthetic using ab initio calculations 
 ( b3lyp/6 - 311+g(2d , p ) for halogenated ethers and mp2/6 - 31g(d ) for propofol ) . 

once the optimized geometry was established , the minimum energy distance 
 between a water molecule and every nonequivalent , accessible atom 
 in the molecule , was calculated and used to assign a partial charge 
 to each atom in the molecule ( see methods for 
 a description of accessible atoms ) . in order to assess the reliability 
 of the current models
, we have calculated several physicochemical 
 properties ( dipole moment , density , hvap , gsolvh2o , and gsolvoil ) for 
 each anesthetic model and have compared them to experimentally determined 
 values ( tables 1and 2 ) . in all cases ,
the parametrized model reproduced experimentally 
 measured properties within acceptable ranges as described in detail 
 below .
standard error in free energies 
 calculated using bar . standard error 
 in log  k values calculated using monte carlo 
 bootstrap error analysis .
although no data for free energy 
 of solvation in an aqueous solution exist for propofol , by use of 
 the aqueous solubility of 160 mg / l   and a vapor pressure of 0.01 mmhg , the 
 free energy of solvation could be predicted . in order to judge whether 
 the equilibrium geometry and charge distribution 
 reflected a chemically relevant state , the magnitude of the dipole 
 moment
the experimental dipole moment is generally determined in 
 the gas phase ; therefore , to reproduce the properties of the anesthetics 
 in the condensed phase ( i.e. , the relevant phase in biological simulations ) , 
 the models should overestimate the experimental dipole moment by up 
 to 20% . the results in table 1 show that the dipole moments of all anesthetics 
 overestimate the experimental gas phase dipole by 1520% , with 
 propofol being overestimated the most ( 20.0% ) and sevoflurane overestimating 
 the experimental dipole moment the least ( 16.7% ) .
once the charge 
 distribution and equilibrium geometry were validated 
 against the experimental dipole moment , the parameters for bonds and 
 angles were assigned based on a transformation of the ab initio hessian 
 matrix .
following assignment of bond and angle terms , torsion scans 
 were conducted for each molecule with a step size of 3 radiating 
 180 in each direction from the equilibrium geometry .
these scans 
 generated a potential energy surface for the molecule which was fit 
 by the dihedral parameters . with a full set of parameters
, we then 
 calculated the density and hvap ( eq 2 ) for each anesthetic from bulk phase 
 md simulations and compared them to experimental values ( table 1 ) .
the parameters were deemed of sufficient accuracy 
 to calculate gsolvoil if the density was within 5% error 
 of the experimental value and heat of vaporization was within a kt ( 0.59 kcal / mol at simulation temperature of 298 k ) of 
 the experimental value , standards set forth in parametrizing the charmm 
 general force field . the density , in 
 most cases , was consistently below 5% error .
the density is exquisitely 
 sensitive to lj terms , and the accuracy of the density for these molecules 
 justifies our transfer of these parameters from the charmm general 
 force field .
desflurane and isoflurane 
 had the least accurate hvap ; however , 
 the values were only off by 0.42 kcal / mol , less than a kt different from the experimental value . in the case of propofol , 

the density was 7.9% lower than experimentally determined , but the 
 hvap was well within the range 
 of experimental values , so the parameters were deemed acceptable . 

comparison to a previous propofol model , which included cmap terms to adjust the torsional term between the 
 phenol hydroxyl and isopropyl groups , shows that the parameters developed 
 here perform slightly better in terms of dipole moment , density , and 
 heat of vaporization ( table s2 ) .
it is 
 unclear , however , how these two models will perform with respect to 
 calculating protein 
all anesthetics , except 
 isoflurane , were able to be parametrized 
 in one iteration .
however , the first set of parameters for isoflurane 
 produced a density of 1.522 g / ml ( 1.7% error ) and the hvap was 8.79 kcal / mol , which differs from the 
 experimental value by 1.18 kcal / mol .
because larger halogens such 
 as the chlorine atoms in isoflurane are quite polarizable , and they 
 are more strongly affected by the molecular contexts , one expects 
 a higher degree of variability in their interaction with the environment . 

this in turn results in reduced transferability for parameters for 
 large halogens between different compounds .
a parameter set developed 
 for the electronic distribution of these atoms is more sensitive to 
 the chemical environment , and the parameters for these atoms are therefore 
 even less transferable than c , n , and o. the heat of vaporization 
 was too large in the first iteration of isoflurane parametrization , 
 but the dipole moment was accurate , suggesting that electrostatic 
 interactions are reasonable .
accordingly , and in order to achieve 
 a better heat of vaporization , we opted to modify the lj potential 
 well term (  ) for the chlorine atom (  = 0.343 
  0.225 ) . by use of the modified lj terms , a complete 
 second round of parametrization
was performed ( including geometry 
 optimization , charge distribution , and assignment of bonded terms ) , 
 yielding an isoflurane model that reproduced both density and hvap measures to within the acceptance criteria 
 discussed above ( i.e. , within 5% error for density and a kt for hvap ) ( table 1 ) .
while being able to reproduce bulk 
 properties is essential to an accurate parameter set , the bulk phase 
 is not necessarily representative of the biological environment and 
 further quantification of parameter assessment is needed . in biological 
 environments , anesthetics are likely to be found in aqueous solution , 
 partitioning into membranes or interacting with protein lumen
. therefore , 
 to assess how well each anesthetic behaves in these different environments , 
 the solvation free energies of each model in both aqueous and nonpolar 
 environments were calculated and compared to experimental measurements 
 ( figure 2 ) . computed free energy change as a function 
 of the coupling parameter , 
  , in both water ( blue ) and dodecane ( orange ) .
the latter mimics 
 an olive oil environment for desflurane ( top left ) , isoflurane ( top 
 right ) , sevoflurane ( bottom left ) , and propofol ( bottom right ) .
the 
 curve shown for each represents the average of 10 decouple / recouple 
 cycles , with the average and error calculated using the bennett acceptance 
 ratio .
the free energy of solvation of 
 each molecule was calculated for 
 water ( gsolvh2o ) and dodecane , the latter of 
 which acts as a mimic for oily environments ( gsolvoil ) .
this 
 was done using alchemical free energy perturbation ( fep ) , an expedient , albeit nonphysical , method to calculate the free energy 
 between two chemically distinct states . in this method ,
the molecule 
 is initially placed in either a box of water or dodecane and the anesthetic 
 is decoupled from its environment by slowly and incrementally scaling 
 the electrostatic and van der waals interactions between the molecule 
 and the environment to zero . with each increment , the free energy 
 cost to transform between these intermediate , nonphysical states can 
 be accurately estimated .
because of the immense computational cost 
 of calculating free energies ( a total of 2 s was needed for 
 fep data presented here ) , validation of the parameter set against 
 experimentally measured gsolvh2o and gsolvoil was the last step in determining if the parameters assigned to the 
 model were accurate .
the parameters were accepted if the absolute 
 difference in free energy was less than 1.5 kcal / mol and the sign 
 of the calculated free energy and the experimental free energy were 
 the same .
table 2 shows the change in 
 free energy 
 associated with the solvation of each anesthetic in water or dodecane . 
 in each case , the model parameters that were judged acceptable using 
 bulk properties produced accurate solvation free energies , some even 
 within a kt of experimental values .
the inhaled anesthetics 
 ( desflurane , isoflurane , sevoflurane ) are completely hydrophobic , 
 having gsolvh2o > 0 while gsolvoil < 
 0 .
in contrast , our data suggest that solvation of propofol is favored 
 in both environments but that propofol is more soluble in an oily 
 environment compared to an aqueous environment . in the case of propofol , 

solvation free energy data for water and oil are not currently available 
 for comparison .
however , the free energy of solvation in water can 
 be predicted from aqueous solubility and vapor pressure . by utilizing an aqueous solubility of 160 mg /
l   and a vapor pressure of 0.01 mmhg , the predicted gsolvh2o is 4.39 kcal / mol , which is extremely close to the value 
 of 3.93  0.30 kcal / mol computed for our model .
reproduction 
 of gsolv values in both water and 
 dodecane demonstrates that the parameter set is quite robust and of 
 such quality to be used in a variety of simulation protocols ( e.g. , 
 protein  anesthetic interaction , free energy calculations , binding 
 constant determination ,  flooding 
although 
 there is now a general consensus that anesthetics exert their effects 
 on plgics , there is currently a debate as to where the anesthetics 
 bind and how they exert their effects .
overton hypothesis suggests that 
 the membrane has some role in the binding and/or action of anesthetics . 
 in order to better understand how anesthetics interact with biological 
 membranes , we simulated the partitioning of anesthetics into a popc 
 membrane and measured the effect on membrane structure .
while the 
 partitioning of anesthetics between phases has been studied computationally 
 for various simplified systems , such as water / hexane and water / vacuum interfaces , their 
 behavior in physiologically relevant systems , such as that presented 
 here , has not been previously explored to the best of our knowledge . using umbrella sampling , we first calculated the free energy of 
 partitioning of each anesthetic into a popc membrane . our calculations 

demonstrate that there is a significant energetic minimum of approximately 
 45 kcal / mol relative to the bulk aqueous solution ( figure 3 ) where the fatty acyl tails meet the glycerol backbone . 
 at this depth in the membrane ,
the hydrocarbons composing the fatty 
 acyl tails provide a hydrophobic environment , while the ester moieties 
 and phosphate groups provide some polar nature that can interact with 
 either the ether on the halogenated ether species or the hydroxyl 
 group of propofol .
this energetic well suggests that anesthetics preferentially 
 partition into the membrane prior to binding plgics .
while the anesthetics 
 show a distinct preference for interfacial regions , there seems to 
 be almost no preference for the membrane core over bulk aqueous solution 
 or vice versa , with the largest difference in free energy being a 
 1 kcal / mol stabilization of sevoflurane in the membrane core compared 
 to aqueous solution . because detailed free energy profiles are not 
 measurable experimentally , the free energy profiles were transformed 
 ( figure 3 ) to popc / water partition coefficients , 
 in the form of log  k , and compared with experimentally 
 measured octanol / water partition coefficients .
as can be seen in table 2 , we are able to capture the partitioning of these 
 molecules into a nonpolar phase very well . some discrepancies arise , 
 including the relative trend of the magnitude of the partition coefficient 
 among the four anesthetics , and are potentially due to the comparison 
 of octanol / water partition coefficients ( experimental ) to popc / water 
 coefficients ( calculated ) and the differing properties of octanol 
 and popc .
the color of the atom groups in 
 this image corresponds to the color of the curves in the density profiles . 

( b ) density profile of the simulated systems used to demarcate the 
 regions of the membrane for analysis of anesthetic 
, total popc density is shown as the black dashed 
 line and water density is shown as the light blue line .
popc density 
 was further subdivided into tail ( gray ) , glycerol ( red ) , phosphate 
 ( gold ) , and choline ( blue ) density .
the colors of the curves correspond 
 to the color of the atoms shown in ( a ) .
( c ) pmf for inserting desflurane 
 ( blue ) , isoflurane ( green ) , sevoflurane ( orange ) , and propofol ( red ) 
 into the membrane .
all anesthetics show a distinct energy minimum 
 at the interfacial region , while there is negligible energy difference 
 between the bulk aqueous environment and midpoint of the lipid bilayer . 

( d ) representative snapshot of ( starting on left and moving right ) 
 desflurane , isoflurane , sevoflurane , and propofol in the umbrella 
 sampling simulations showing the low energy conformation at the amphipathic 
 boundary of the membrane .
while it is known that general anesthetics affect the conduction 
 of plgics , it is still unclear if anesthetics have drastic effects 
 on membrane structure and whether these effects might indirectly lead 
 to ion channel modulation . to better understand how anesthetics interact 
 with the membrane
, we simulated a high concentration ( initial aqueous 
 concentration of 100 mm or anesthetic / lipid ratio of 1:5 ) 
 of each anesthetic for 100 ns in a popc membrane system .
as the anesthetics 
 partition into the membrane rather rapidly , with 6090% of 
 anesthetics partitioned within 50 ns ( figure s1 ) , the length of the simulations was sufficient to observe partitioning 
 of a majority of anesthetic and to explore the effect of these drugs 
 on a popc membrane . to measure the effect of such a high concentration 
 of the anesthetic on membrane structure , we measured several physical 
 properties of the popc membrane ( atomic density profiles , area per 
 lipid , order parameters , dipole potential , and lateral stress profiles ) 
 after the anesthetics had partitioned into the membrane .
as can be 
 seen from figure 4 , the anesthetic has no measurable 
 effect on distribution of atomic groups in the membrane , as all atomic 
 groups have no shift in peak density or any difference in distribution 
 width upon partitioning of anesthetics .
moreover , the high concentration 
 of anesthetic does not seem to affect the local structure of the headgroups 
 or the palmitoyl / oleyl tails , with no changes in the order parameters 
 detected for any anesthetic simulated ( figure 4b , c ) . the area per lipid ( al )
was measured 
 for the membrane in both the absence and presence of anesthetics . 

there was no detectable change , as the al for the pure popc membrane was 65.05  3.15  versus an average of 65.27  0.42  across 
 all four anesthetic simulations ( desflurane , 65.47  3.19  ; isoflurane , 65.30  3.18  ; sevoflurane , 
 64.67  3.13  ; propofol , 65.62  3.19 
  ) .
( a ) atomic density profiles for membrane components in 
 the ( starting 
 left and moving right ) popc , desflurane , isoflurane , sevoflurane , 
 and propofol systems . the atomic densities of the tails ( black trace ) , 
 glycerol ( red trace ) , phosphate ( gold trace ) , and choline ( blue trace ) 
 moieties are shown as an average over the last 50 ns of the trajectory 
 ( accounting for insertion of 6090% of anesthetics )
. the atomic 
 density of water is also shown ( cyan trace ) and measures 1.0 g / ml , 
 which is the experimentally determined density of water at 298 k. 
 order parameters for both the headgroup ( b ) and lipid tails ( c ) are 
 shown for popc ( black ) , desflurane ( blue ) , isoflurane ( green ) , sevoflurane 
 ( orange ) , and propofol ( red ) .
additionally , the order parameters for 
 palmitoyl ( dashed line ) and oleyl ( solid line ) tails are shown separately . with those three measures of membrane 
 structure showing no change , 
 we made two additional measurements to see if the effect of anesthetics 
 on the membrane is more subtle . because anesthetics partition to the 
 membrane  water interface ( figure 5a 
d ) ) , 
 it is possible that anesthetics might perturb the organization of 
 water molecules at this layer , thereby disturbing membrane electrostatics 
 and structure . therefore , to account for this , we measured the dipole 
 potential of the popc membrane , a net positive potential at the center 
 of the membrane due to the ordering of water molecules at the membrane 
 interface , across all five simulations ( figure 5 ) .
these plots show that , to within error , there is no significant 
 change in the dipole potential of the membrane , suggesting that there 
 is no gross disturbance of water structure at the membrane surface . 

it is important to note that dipole calculations from md simulations 
 are known to overestimate the dipole potential , due to the lack of polarizability in classical force fields . 

comparison between the profiles calculated , however , remains valid , 
 as each simulation suffers from this deficiency .
moreover , since these 
 small molecules are partitioning into the membrane , it is conceivable 
 that partitioned anesthetics affect the interlipid dynamics in the 
 membrane , the effects of which can be measured in lateral stress profiles . 
 as can be seen in figure 5
, there is no statistically 
 significant difference in the lateral stress profile between popc 
 only and popc / anesthetic membranes for all cases . while these measures 
 do not constitute an exhaustive study of membrane behavior
, they do 
 provide significant evidence that anesthetics do not perturb membrane 
 structure , agreeing with previous experiments that have shown no effect 
 of anesthetics on membrane structure at clinically relevant concentrations .
noble gas molecules , such as xenon which has 
 been used as an anesthetic , have been shown to disturb membrane structure 
 in previous simulations .
these species , 
 however , are vastly different from the general anesthetics studied 
 here in that xenon is a hydrophobic atom , and much smaller in size 
 compared to the anesthetic molecules studied here , that readily partitions 
 to the membrane center , not the amphipathic membrane interface .
therefore , 
 it should be noted that the anesthetics studied here did not show 
 any disturbance of membrane structure , but anesthetics with completely 
 distinct chemical structures might have an effect on gross membrane 
 structure .
atomic density profiles 
 for desflurane ( a ) , isoflurane ( b ) , sevoflurane 
 ( c ) , and propofol ( d ) . the density of each anesthetic is further subdivided 
 into its constituents to show detail on the distribution of these 
 molecules in the membrane .
plot of the dipole potential ( e ) and the 
 lateral stress profile ( f ) of the popc membrane in the popc only ( black ) , 
 desflurane ( blue ) , isoflurane ( green ) , sevoflurane ( orange ) , and propofol 
 ( red ) .
the dipole potential and lateral stress profiles are symmetric 
 about the center of the membrane .
the error bars denote the standard deviation 
 from averaging the profiles across the two halves .
md simulations offer 
 an unprecedented opportunity to probe the 
 dynamics of the effect of anesthetics on ion channels with high spatiotemporal 
 resolution , but a lack of tested parameters has significantly hindered 
 progress in this arena . here , we have presented atomic - detailed models 
 for use in md simulations of four ubiquitous modern general anesthetics : 
 desflurane , isoflurane , sevoflurane , and propofol .
the models have 
 been shown to reproduce multiple experimentally determined physicochemical 
 properties that demonstrate the robustness of these models and their 
 ability to accurately describe the behavior of the anesthetics in 
 different environments , including aqueous and membranous environments . 
 thus , these parameters are suitable for use in multiple simulation 
 protocols such as probing protein  anesthetic interaction , free 
 energy and binding constant calculations , and identification of drug 
 binding sites , among others . using the parameters generated
these studies show that there is a 
 significant energy minimum ( 45 kcal / mol ) at the interface 
 between the tail region of the lipids and their glycerol moiety , supporting 
 previous studies that showed the same phenomenon in hexane / water and 
 vacuum / water systems .
in contrast to the meyer  overton hypothesis , however , there 
 is negligible difference in energy between the membrane core and aqueous 
 solution ( less than 1 kcal / mol ) .
it is interesting to note that the 
 glic anesthetic binding site lies at this interfacial level and these 
 data suggest an entrance mechanism whereby the anesthetic penetrates 
 the membrane prior to binding the protein .
additionally , flooding 
 the membrane with anesthetics does not perturb any of the membrane 
 structural measurements , including atomic density , order parameters , 
 area per lipid , dipole potential , and lateral stress profile of the 
 popc membrane .
while this suggests that anesthetics do not affect 
 membrane structure , more careful studies of membrane physical properties 
 in the presence of anesthetics need to be conducted .
this study provides 
 parameters for a set of anesthetics that have been cocrystallized 
 with plgics , allowing for rapid and expanded 
 study of anesthetic modulation of pentameric ligand - gated ion channels .
 while small molecules have been used 
 to induce anesthesia in a 
 clinical setting for well over a century , a detailed understanding 
 of the molecular mechanism remains elusive . in this study , we utilize 
 ab initio calculations to develop a novel set of charmm - compatible 
 parameters for the ubiquitous modern anesthetics desflurane , isoflurane , 
 sevoflurane , and propofol for use in molecular dynamics ( md ) simulations . 
 
 the parameters generated were rigorously tested against known experimental 
 physicochemical properties including dipole moment , density , enthalpy 
 of vaporization , and free energy of solvation . in all cases , 
 the anesthetic 
 parameters were able to reproduce experimental measurements , signifying 
 the robustness and accuracy of the atomistic models developed . 
 the 
 models were then used to study the interaction of anesthetics with 
 the membrane . 
 calculation of the potential of mean force for inserting 
 the molecules into a popc bilayer revealed a distinct energetic minimum 
 of 45 kcal / mol relative to aqueous solution at the level of 
 the glycerol backbone in the membrane . 
 the location of this minimum 
 within the membrane suggests that anesthetics partition to the membrane 
 prior to binding their ion channel targets , giving context to the 
 meyer  overton correlation . moreover , md simulations of these 
 drugs in the membrane give rise to computed membrane structural parameters , 
 including atomic distribution , deuterium order parameters , dipole 
 potential , and lateral stress profile , that indicate partitioning 
 of anesthetics into the membrane at the concentration range studied 
 here , which does not appear to perturb the structural integrity of 
 the lipid bilayer . these results 
 signify that an indirect , membrane - mediated 
 mechanism of channel modulation is unlikely . 
Introduction
Methods
Results
and Discussion
Conclusions