Atomistic models of general anesthetics for use in in silico biological studies.

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 4-5 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

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,[1] known as the Meyer–Overton hypothesis (or Meyer–Overton correlation). Recent studies[2−4] 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,[5,6] 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.[2−4] 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)[7−12] and ELIC (E. chrysanthemi ligand-gated ion channel),[13,14] 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,[9,12] mutagenesis,[15,16] photolabeling,[17,18] and molecular dynamics (MD);[19,20] 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[19−25] with a high degree of spatial and temporal resolutions. Moreover, the solution of atomic-resolution (2.9–3.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.

Crucial to these simulations are well-parametrized models of anesthetics. 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,[26,27] 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[28−30] 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.[20,31,32] However, with the exception of the isoflurane model developed by Hénin et al.,[32] 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.[26,33−36] 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.[37] 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.

Figure 1

Chemical structure and formula of the four anesthetics parametrized in this study.

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 4–5 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.

Methods

The recent development of the force field toolkit (ffTK)[38] 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, see Mayne et al.[38]

Parametrization

Two levels of theory were used to properly model the different characteristics of the heavily halogenated inhaled anesthetics and the large intravenous anesthetics. The Gaussian09 program[39] was used for all ab initio calculations. Ab initio calculations on desflurane, isoflurane, and sevoflurane (inhaled anesthetics) utilized the hybrid B3LYP density-functional theory[40,41] 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[42,43] and has been used to parametrize halogenated molecules previously.[30,32,44] For propofol, second-order Møller–Plesset theory (MP2)[45] 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,[38] which is based on CHARMM general force field parametrization principles,[26] so as to be consistent with the CHARMM force field. ffTK does not currently support the optimization of Lennard-Jones (LJ) parameters de novo. These parameters, however, are usually the most transferable between chemical species. 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, all LJ parameters were unmodified. 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 sp3 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%[26] and an overestimation of the QM dipole by 10–20% was used as the acceptance criterion to continue to bonded terms. If the dipole moment was outside of this range, another iteration of charge optimization was performed followed by another dipole moment measurement. This process was repeated until the computed dipole moment was within acceptable ranges. The dipole moment was calculated according towhere the bars denote the vector length, N denotes the number of atoms in the molecule, q denotes the charge of the ith atom, and r denotes the position vector of the ith atom. The bonds and angles terms were calculated directly from the Hessian. Scans of each torsion were made at 3° steps, and the dihedral parameters were fit to the resulting potential energy surface (PES).[38,46] 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).

Bulk Properties

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. Initial grid spacing was based on the experimentally measured density of the molecule. 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. Each simulation was repeated five times to ensure proper sampling and averaging.

The enthalpy of vaporization (ΔHvap) was measured according to the following equation:where ⟨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 P⟨Vliq⟩ term is negligible[26] 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.[47] This is done according towhere ⟨Uliqinter⟩ is the average intermolecular energy of the liquid, ⟨Uliq⟩sys is the average potential energy for the entire system (including inter- and intramolecular interactions), and ⟨Uliqintra⟩ 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,[26] using one 100 ns gas phase calculation to calculate the intramolecular energy of the molecule. These values are included in Table 1. In the case 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.

Table 1

Bulk Properties of Parametrized Anesthetics

	dipole (D)		density (g/mL)		ΔHvap (kcal/mol)
	calca	exp	calcb	exp	calcc	exp
desflurane	3.39	2.87	1.48 ± 0.01	1.46	8.02 ± 0.03 (8.21)	7.60
isoflurane	2.91	2.47	1.48 ± 0.02	1.49	8.03 ± 0.01 (8.27)	7.61
sevoflurane	2.72	2.33	1.46 ± 0.01	1.52	7.87 ± 0.03 (7.59)	7.89
propofol	1.92	1.60	0.98 ± 0.01	1.03	16.27 ± 0.05 (16.05)	12.27–20.58d

As discussed, the calculated dipole is expected to be 10–20% higher than experimental gas phase dipole moment to replicate condensed phase properties. Because dipole is calculated using only the equilibrium geometry, there is no uncertainty attached to this measurement.

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,[26] 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.

Free Energy of Solvation

Taking advantage of readily available experimental data on the free energy of solvation in both water (ΔGsolvH) 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.[32] 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[20,32]) was used to calculate both ΔGsolvH 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”,[48] the Zacharias soft-core potential[49] 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).[50] 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

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.[51−53] 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–1 Å–2. 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,[54,55] 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:[56,57]where 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 Å).

Membrane/Anesthetic Simulations

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 Å2 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.[58] This gave the system initial dimensions of 96 × 96 × 120 Å3 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 Å3, 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.

Simulation Protocols

All simulations were performed using NAMD2[59] together with the set of force fields parameters developed herein for each anesthetic (Supporting Information) and the CHARMM36[26,35] parameters for the TIP3P water model[60] 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. Constant pressure was maintained using the Nosé–Hoover Langevin piston method,[61,62] and constant temperature was maintained using a Langevin damping coefficient (γ) of 1 ps–1. 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[63] with a grid density greater than 1 Å–3. 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.

Results and Discussion

Validation of the Parametrized Models

The recent availability of multiple atomic-resolution structures of GLIC,[7−12] ELIC,[13,14] and GluCl (a eukaryotic Cys-loop receptor from C. elegans)[64] 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[37] 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,[20,32] along with a desflurane model for GROMOS.[65] 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, ΔGsolvH, 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.

Table 2

Calculated Free Energies of Parametrized Anesthetics

	ΔGsolvH2O (kcal/mol)		ΔGsolvOil (kcal/mol)		log K
	calca	exp	calc	expb	calcc	exp
desflurane	0.24 ± 0.20	0.92	–2.45 ± 0.42	–1.74	2.48 ± 0.88	n.d.
isoflurane	1.03 ± 0.25	0.31	–2.77 ± 0.36	–2.72	3.24 ± 0.81	2.22
sevoflurane	1.47 ± 0.28	0.27	–3.12 ± 0.32	–2.28	2.64 ± 0.96	3.00
propofol	–3.93 ± 0.30	–4.39d	–7.28 ± 0.38	–6.32	2.14 ± 1.01	3.63

Calculated values are presented as the mean ± standard error. Standard error in free energies calculated using BAR.[50] Standard error in log K values calculated using Monte Carlo bootstrap error analysis.

Experimental values are determined for anesthetics in olive oil. Calculations were performed using dodecane as an approximation of olive oil.

Calculated values are presented as the mean ± standard error. 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 [68] and a vapor pressure of 0.01 mmHg,[69] the free energy of solvation could be predicted.[67]

In order to judge whether the equilibrium geometry and charge distribution reflected a chemically relevant state, the magnitude of the dipole moment was calculated for each molecule according to eq 1. 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%.[26] The results in Table 1 show that the dipole moments of all anesthetics overestimate the experimental gas phase dipole by 15–20%, 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.[26] 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.[26] 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,[20] 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–anesthetic interaction energies, which is beyond the scope of the present study.

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).

Robust Anesthetic Models Reproduce Solvation Free Energies in Different Environments

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).

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 (ΔGsolvH) and dodecane, the latter of which acts as a mimic for oily environments (ΔGsolvOil). This was done using alchemical free energy perturbation (FEP),[20,32,48,66] 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 ΔGsolvH 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 ΔGsolvH > 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.[67] By utilizing an aqueous solubility of 160 mg/L [68] and a vapor pressure of 0.01 mmHg,[69] the predicted ΔGsolvH 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” simulations).

Interaction of Anesthetics with the Membrane

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.[9,16,19,20,70,71] Moreover, the trend observed in the Meyer–Overton hypothesis[1,72] 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[31,73] and water/vacuum[32] 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 4–5 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.

Figure 3

(a) Snapshot of the POPC membrane used in the simulations. Water and ions are omitted for clarity. 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–membrane interactions. Here, 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 60–90% 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 Å2 versus an average of 65.27 ± 0.42 Å2 across all four anesthetic simulations (desflurane, 65.47 ± 3.19 Å2; isoflurane, 65.30 ± 3.18 Å2; sevoflurane, 64.67 ± 3.13 Å2; propofol, 65.62 ± 3.19 Å2).

Figure 4

(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 60–90% 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,[35] 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.[74] Noble gas molecules, such as xenon which has been used as an anesthetic, have been shown to disturb membrane structure in previous simulations.[75] 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.[76]

Figure 5

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 two halves of each profile were averaged and presented here. The error bars denote the standard deviation from averaging the profiles across the two halves.

Conclusions

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, we probed the interaction of these molecules with a POPC membrane. These studies show that there is a significant energy minimum (4–5 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.[31,32,77] 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,[9] allowing for rapid and expanded study of anesthetic modulation of pentameric ligand-gated ion channels.