In silico studies of the interactions between propofol and fentanyl using Gaussian accelerated molecular dynamics

Abstract

Fentanyl is a potent opioid analgesic, which for decades has been used routinely in surgical and therapeutic applications. In addition to its analgesic properties, fentanyl also possesses anesthetic properties, which are not well understood. Fentanyl is used in the general anesthesia process to induce and maintain anesthesia in combination with the general anesthetic propofol, which fentanyl is known to potentiate. As the atomic-level mechanism behind the potentiation of propofol is unclear, we have used classical molecular dynamics simulations to study the interactions of these drugs with the Gloeobacter violaceus ion channel (GLIC). This ion channel has been identified as a target for many anesthetic drugs. We identified multiple binding sites using flooding style and Gaussian accelerated molecular dynamics (GaMD) simulations, showing fentanyl acting as a stabiliser that holds propofol within binding sites. Our extensive GaMD simulations were also able to show the pathway by which propofol blocks the channel pore, which has previously been suggested as a mechanism for ion channel modulation. General anesthesia is a multi-drug process and this study provides the first insight into the interactions between two different drugs in the anesthesia process in a relevant biological environment.

Communicated by Ramaswamy H. Sarma

Keywords:

• Ion channel
• molecular dynamics simulation
• enhanced sampling
• anesthetic
• opioid

1. Introduction

The phenomenon of general anesthesia is a medically induced loss of consciousness which renders the patient immobile and non-responsive to painful stimulus, thereby allowing invasive surgical procedures to be performed under the safest possible conditions. The discovery of molecules which possess properties that induce anesthesia has caused significant advances within the field of medicine, but to date a lack of understanding remains of the molecular mechanisms which lead to general anesthesia. Early investigations into the mechanisms of anesthesia were carried out by Meyer (Meyer, 1899) and Overton (Overton, 1901) who proposed a lipid-mediated mechanism due to the correlation between the solubility of anesthetic molecules in oil and their potency. Consequently, studies have been performed on a large range of anesthetic molecules in different solvent environments, (Johansson & Zou, 1999; Taheri et al., 1991) which have followed the hypothesis laid down by Meyer and Overton, but several exceptions have been identified. (Cardoso et al., 1999; Fang et al., 1997; Koblin et al., 1995, 1994; North & Cafiso, 1997) Several studies have proposed that anesthetics cause alteration of specific membrane properties, such as the lateral pressure profile, which will affect membrane proteins whose function depends on a transition between conformational states that is accompanied by a depth-dependent change in cross sectional area. (Cantor, 1999; van den Brink-van der Laan et al., 2004) The specific locations adopted by anesthetics within the cellular membrane have also been investigated by Pohorille et al., who proposed that anesthetics reside at the lipid-water interface. (Pohorille et al., 1998) Although there is a wealth of published work examining the role of the lipid membrane, the main theory which has attracted the most attention in recent years is the direct binding of anesthetics to pentameric ligand-gated ion channels (pLGICs) or cys-loop receptors. pLGICs are of vital importance to life as they are responsible for converting synaptic chemical signals into electrical impulses within the peripheral and central nervous systems (Corringer et al., 2012) and as such they are targets for general anesthetics. (Campagna et al., 2003; Hemmings et al., 2005) Prokaryotic homologues of pLGICs have provided high-resolution crystal structures such as those from Gloeobacter violaceus (GLIC). (Bocquet et al., 2009; Hilf & Dutzler, 2009) Crystal structures of GLIC have also been identified that show binding of various general anesthetics, including ethanol, (Sauguet et al., 2013) ketamine, (Pan et al., 2012) desflurane, (Nury et al., 2011) and propofol, (Fourati et al., 2018; Nury et al., 2011) to the extracellular and transmembrane domains of the protein structure. These discoveries have provided the strongest evidence thus far that anesthetics exert their main function by modulating the activity of these ion channels.

The Gloeobacter violaceus ion channel is a pentameric ligand-gated ion channel which is composed of five identical subunits (S1 - S5) forming a homopentameric ion channel, where each subunit possesses an extracellular domain (ECD) and a transmembrane domain (TMD). The TMD of each subunit consists of four α-helices (M1 - M4) which span the entirety of the phospholipid membrane in which the GLIC is inserted. The M2 α-helices are orientated towards the center of the pore which forms the ion-conducting channel. This is illustrated in Figure 1. Propofol is an intravenous general anesthetic, which is widely used to induce and maintain general anesthesia, often in combination with the opioid analgesic fentanyl, (Bajwa et al., 2010; Bell et al., 1994; Todd et al., 1993) as it facilitates rapid and safe anesthesia with stable hemodynamics and relatively few negative side effects. The role of fentanyl is to provide the analgesic component to the general anesthesia process, although fentanyl does possess anesthetic properties itself and has been used as the main component in general anesthesia. (Sebel et al., 1981; Stanley & Webster, 1978) Fentanyl has also been shown to potentiate the anesthetic effect of propofol, (Vuyk, 1997; Lichtenbelt et al., 2004; Smith et al., 1994) but knowledge of the specific interactions between propofol and fentanyl at the molecular level is lacking. Understanding the interactions between these drugs at their site of action, however, could lead to further development of safer and more efficient anesthetic drugs which will reduce the risk of anesthesia-related complications and help uncover the mechanisms behind this phenomenon.

Figure 1. (A) Graphical representation of the transmembrane α helical domains showing the 5 symmetric subunits and the ion conducting M2 helices (red). (B) GLIC structure in the membrane environment, some lipids removed for clarity of TMD.

In this paper we utilise two enhanced sampling molecular dynamics (MD) simulation techniques, namely conformational flooding and Gaussian accelerated molecular dynamics (GaMD), to probe the interactions between fentanyl and propofol at the GLIC channel. GLIC was chosen for study as several structures have been published that show propofol binding, (Fourati et al., 2018; Nury et al., 2011) and we propose that fentanyl could also act in this channel in combination with propofol. To the best of our knowledge, no MD simulations have been performed on this combination of drug molecules with GLIC. From our simulations we were able to identify multiple sites in GLIC where both drug molecules appear to act and alter the protein mobility and flexibility.

2. Materials and methods

2.1. System preparation and MD simulations

All systems prepared for the MD simulations were based on the high-resolution open-state crystal structure (PDB: 4HFI) of GLIC. (Sauguet et al., 2013) The structure was embedded in a DOPC lipid bilayer and fully hydrated with TIP3P (Jorgensen et al., 1983) water using the CHARMM-GUI web server (Jo et al., 2008; Wu et al., 2014) and data from the orientation of proteins in membranes (OPM) database. (Lomize et al., 2012) DOPC was chosen as it was used previously by the publishers of the crystal structure. (Sauguet et al., 2013) Each system has the approximate dimensions of 128 × 130 x 180 Å where the z-axis is normal to the bilayer, and contains 1 GLIC, 307 DOPC lipids, and ∼54,000 water molecules. Sodium and chloride ions were added to give a 150 mM salt concentration. 50 propofol and 10 fentanyl were added to the aqueous solution for the flooding and GaMD simulations. Each system contained ∼230,000 atoms.

Protonation states of the ionisable residues of the protein were assigned based on the calculations performed by Bocquet et al. (Bocquet et al., 2009) Parameters and atomic point charges for each drug molecule were assigned using the antechamber program which is freely available with AmberTools program, and is compatible with the General Amber ForceField (GAFF) (Wang et al., 2004) for small organic molecules.

All MD simulations were performed using the CUDA enabled Amber16 (Case et al., 2016) molecular dynamics package. The ff14SB forcefield (Maier et al., 2015) was used for the protein, solvent and ions, the lipid14 forcefield (Dickson et al., 2014) was used for the lipids and the GAFF2 forcefield (Wang et al., 2004) was used for the drug molecules. The systems were energy minimised for 20,000 steps (10,000 steepest descent then 10,000 conjugate gradient). Heating was conducted in two stages, where stage 1 involved heating the system to 100 K in the NVT ensemble with the protein, lipids and drug molecules restrained with a force constant of 10.0 kcal/mol Ų. Stage 2 involved heating the system slowly to 303 K in the NPT ensemble, again with the protein, lipids and drugs restrained with a force constant of 10.0 kcal/mol Ų. Equilibration was performed initially with multiple 2 ns simulations where the lipid restraints were gradually reduced from 10 to 0 kcal/mol Ų over 22 ns. This procedure was then repeated while reducing restraints on the protein backbone totalling another 22 ns (44 ns in total). The unrestrained conformational sampling simulations were run for 1 μs and the GaMD simulations were run for 500 ns. 3 repeats for each method were carried out. One simulation was carried out on the pure GLIC structure as a control to determine if the system was constructed correctly, which was run for 500 ns. All simulations were carried out at a constant pressure of 1.0 atm, which was maintained using the anisotropic Berendsen method with a pressure relaxation time of 1.0 ps. The temperature was maintained at 303 K using the Langevin thermostat with a collision frequency of γ = 1.0 ps^{– 1}. Three-dimensional periodic boundary conditions with the usual minimum image convention were used. The SHAKE algorithm was employed to constrain covalent bonds to hydrogen, allowing the use of a 2 fs time step. Electrostatic interactions were treated with the PME method using a cutoff of 10 Å. Data analysis was carried out using VMD, (Humphrey et al., 1996) CPPTRAJ (Roe & Cheatham, 2013, 2018) and in house scripts. Images were rendered using UCSF Chimera (Pettersen et al., 2004) and VMD. (Humphrey et al., 1996)

2.2. Flooding simulations

Simulation studies on the GLIC system in the past have mostly focused on the anesthetic binding site from the published crystal structure, (Nury et al., 2011) which often involved docking the drug molecule within the observed site. However, this methodology does not take into account the pathway of the drug molecule, or any other binding sites that can be occupied. We have investigated the possibility of these alternative binding sites and pathways by implementing three independent 1 μs flooding simulations in which 50 propofol molecules and 10 fentanyl molecules were added to the solvent phase of the last frame of the pure GLIC control simulation, which was used as the starting structure; this methodology is similar to other studies conducted on different anesthetics. (Arcario et al., 2017; Brannigan et al., 2010) The drug molecules were added to the system at random initial positions using the gmx insert-molecules tool available with the GROMACS package. (Van Der Spoel et al., 2005) Propofol was used in excess in our simulations to mimic the higher concentration of propofol in the general anesthesia process. We note that the ratio of opioid to anesthetic is higher here than would be used in surgery, but the goal of using a higher ratio is to increase the probability of finding a site of opioid-anesthetic interaction during the simulations. The MD flooding approach provides a dynamic, flexible and physical dock which is a significant advantage over traditional docking procedures. (Brannigan et al., 2010; Raju et al., 2013)

2.3. Gaussian accelerated molecular dynamics

GaMD is an accelerated MD method which allows unconstrained sampling and accurate reconstruction of the free-energy profiles, thereby allowing long timescale events which would normally occur over hundreds of microseconds to milliseconds to be probed in more realistic simulation times. (Miao & McCammon, 2017; Palermo et al., 2017) A full mathematical description of the methodology can be found in the original publication (Miao et al., 2015) and only a brief description will be given here. The enhanced conformational sampling provided by GaMD is obtained by adding a harmonic boost potential which reduces the energy barriers in the system. When the total system potential V($\stackrel{\rightharpoonup }{r}$) is lower than the reference energy E, the new potential V*($\stackrel{\rightharpoonup }{r}$) of the system is calculated as: (1) $$V^{*}(\stackrel{\rightharpoonup }{r})=V(\stackrel{\rightharpoonup }{r})+\Delta V(\stackrel{\rightharpoonup }{r})$$(1) (2) $$\Delta V(\stackrel{\rightharpoonup }{r})=\{\begin{array}{cc}\frac{1}{2}k{(E-V(\stackrel{\rightharpoonup }{r}))}^2,\hfill & V(\stackrel{\rightharpoonup }{r})<E\hfill \\ 0,\hfill & V(\stackrel{\rightharpoonup }{r})\ge E,\hfill \end{array}$$(2) where k is the harmonic force constant. The two adjustable parameters E and k are automatically determined based on three principles. First, for any two arbitrary potential values V${}_A^{*}$($\stackrel{\rightharpoonup }{r}$) and V${}_B^{*}$($\stackrel{\rightharpoonup }{r}$) found on the initial energy surface, if V${}_A^{*}$($\stackrel{\rightharpoonup }{r}$) < V${}_B^{*}$($\stackrel{\rightharpoonup }{r}$), ΔV should be a monotonic function which allows the relative order of the biased potential to remain constant. Second, if V${}_A^{*}$($\stackrel{\rightharpoonup }{r}$) < V${}_B^{*}$($\stackrel{\rightharpoonup }{r}$), then the potential difference on the flattened energy surface should be smaller than that of the original. By a combination of the first two principles, we obtain: (3) $$V_{\mathit{max}}\le E\le V_{\mathit{min}}+1/k,$$(3) where V_min and V_max are the system minimum and maximum potentials. The standard deviation of ΔV has to have a narrow distribution, so accurate re-weighting by cumulant expansion to the second order can be performed: ${\sigma }_{\Delta V}$ = k(E - V_avg)${\sigma }_V$ $\le$ ${\sigma }_0,$ where V_avg and ${\sigma }_V$ are the average and standard deviation of the system potential energies. ${\sigma }_{\Delta V}$ is the standard deviation of ΔV and ${\sigma }_0$ is a user-specified upper limit for accurate re-weighting.

We have carried out extensive tests to obtain suitable parameters for our systems. The system threshold energy was set to E = V_max for all GaMD simulations. The dual-boost scheme was used, in which an acceleration potential is applied to the torsional terms and to the entire system potential. V_max, V_min, V_avg and ${\sigma }_V$ were obtained from an initial 8 ns NPT simulation without boost potential. Each GaMD simulation then underwent a 40 ns run in which the boost potential was updated every 1.6 ns, thus reaching equilibration values. Production runs were carried out for 500 ns for all three simulations, in which each was started with random particle velocities. The same number of drug molecules as in the flooding simulations were included in the GaMD simulations. Trajectory frames were saved every 1 ps.

3. Results and discussion

3.1. Flooding

Almost instantly during the production runs, both drug molecules begin to partition into the membrane (Figure 2), often forming small clusters in which multiple propofol molecules bind to one fentanyl molecule. To assess the structural stability of the GLIC structure within the membrane, the root mean squared deviation (RMSD) of the backbone of the protein was calculated to be less than 3.0 Å over the full length of the simulations and the average RMSD was 1.922 ± 0.209 Å, based on which the simulations were deemed to be stable as the protein remained in its native structure. Multiple binding regions were identified for propofol, where protein residues interacted with the anesthetic for several ns of the simulation. Of these potential interaction sites, only one was observed over multiple subunits for several hundred ns. This interaction site was located in the transmembrane domain (TMD) close to the experimentally identified site, see Figure 3. The other sites, where shorter interactions were observed, most likely represent non-specific adhesion of propofol to the protein surface rather than a functional binding site. The site where most stability is observed is identical to that observed in a study carried out by Arcario et al., (Arcario et al., 2017) in which they observed desflurane binding to GLIC by forming a hydrogen-bond to Tyr-254 during their simulations. The stability associated with propofol at this site (Figure 3B) is due to the formation of a hydrogen-bond between the hydrogen of the Tyr-254 hydroxyl group and the oxygen of the propofol hydroxyl group, which was calculated for the propofol interacting at subunit 3. The hydrogen-bond remains dominant over the full time in which the propofol molecule remains stable, which is approximately 400 ns. This stability was shown by calculating the distance between the propofol molecules observed to interact at this site in each subunit, and the center of mass (COM) of the experimental anesthetic binding site, in which propofol has been shown to bind. (Nury et al., 2011) The S3 propofol remains bound from ∼ 200 ns to ∼ 650 ns, where the hydrogen-bond to Tyr-254 contributes to propofol’s overall stability here. Three propofol molecules were observed to bind to the same site at the same time, and this asymmetric ligand binding has been shown to play a role in channel function. (Mowrey et al., 2013) Our distance calculations have shown that this binding site holds the propofol molecules at ∼ 10 Å from the deeper transmembrane site, indicating a possible transition pathway, although this was not observed to happen in any of our simulations. Throughout the flooding simulations, multiple bound anesthetics were observed at multiple subunits which helps to explain why large amounts of general anesthetic are required to induce the desired clinical effect. This observed behaviour strengthens the argument that these drugs indiscriminately bind to several targets. (Eckenhoff et al., 2002; Harrison, 2000; Herold & Hemmings, 2012)

Figure 2. (A) First and (B) final frame of 1 μs of production simulation. Blue represents fentanyl and red represents propofol.

Figure 3. (A) Distance between interacting propofol molecules and COM of experimental binding site. (B) Snapshot showing the binding site in which propofol (licorice representation) interacts with Tyr-254 (licorice representation) and Asn-307 (yellow). The anesthetic binding site residues are shown; Tyr-193 (green), Ile-202 (blue), Met-205 (red), Val-242 (orange), Thr-255 (pink) and Ile-259 (blue). M1-M4 are the transmembrane helices located in each subunit. Hydrogen-bond between propofol and Tyr-254 shown by green dashed line.

An interesting observation made during our flooding simulations was propofol occupying the exact same binding site in the extracellular domain that the general anesthetic ketamine has been shown to occupy. (Pan et al., 2012) This site is approximately 10 Å below the orthosteric agonist-binding site in the extracellular domain. Ketamine binding in this site was also shown to modulate the function of the GLIC channel. (Ion et al., 2017) Previously identified sites for propofol were located in the upper transmembrane domain within a subunit (Nury et al., 2011) which is very different from that identified for ketamine, which is located ∼ 10 Å below the orthosteric agonist-binding site. We observed propofol binding in this extracellular site within the first 100 ns of simulation and it remains in this site for the rest of the 1 μs simulation, totalling ∼ 900 ns of residence time within the site. Propofol remains fairly mobile within the site with a RMSD of 5.0 ± 1.4 Å, which is due to the relatively small molecular volume of propofol (∼ 191 ų) compared to the ketamine site (∼ 248 ų) The residues within 4 Å of propofol displayed a RMSD of 1.7 ± 0.5 Å, which shows that the structure in this flexible loop region is fairly resilient to propofol binding. Identification of multiple binding sites presents a compelling case for the allosteric action of molecules which possess anesthetic properties.

3.2. GaMD

Our GaMD simulations were able to capture complete binding of one propofol molecule into the channel pore just above the hydrophobic gate region. Binding was captured within ∼ 50 ns in all three independent GaMD simulations, whereas this binding mode was not observed in any of our 1 μs flooding simulations. We found that propofol diffuses from the bulk solvent to the β-3 loop in the upper extracellular domain of subunit one, before moving to the β-8 loop between subunits one and two in the extracellular domain, and finally into the center of the M2 pore above the Ile-239 gating residues, see Figure 4B. By aligning the C-terminal domain of the GLIC, the RMSD of the diffusing propofol molecule relative to its position in the pore in the last simulation frame reaches a minimum of 1.1 Å. It forms hydrophobic interactions with several residues around the pore opening, see Figures 4A and C.

Figure 4. GaMD simulations captured pore binding of propofol: (A) 2 D (N_contact, propofol RMSD) PMF calculated by reweighting the 500 ns GaMD simulation. (B) Pathway of propofol diffusion into the M2 helix pore in all three GaMD simulations. (C) Hydrogen-bonds formed between propofol and the GLIC, where bonds formed after 10000 steps are in the pore. Protein refers to short lived hydrogen-bonds formed on the protein exterior. Data combined from all three replicates. (D) Final frame snapshot showing propofol hydrogen-bonding to Thr-236 in S1.

The boost potential applied during the GaMD simulations follows a Gaussian distribution and its distribution anharmonicity γ equals to 1.94 × 10^{– 2}. The average and standard deviation of ΔV are 15.3 kcal/mol and 4.7 kcal/mol, respectively. Such narrow distribution will ensure accurate re-weighting for the free energy calculations using cumulant expansion to the 2^nd order. Using the RMSD of propofol relative to the final frame pose and the number of protein heavy atoms that are within 5 Å of propofol (N_contact), a 2 D potential of mean force (PMF) profile was calculated by re-weighting the three 500 ns GaMD simulations using the PyReweighting toolkit. (Miao et al., 2014) As the pore blocking pathway was observed in all three replicates, the data were combined to produce the PMF in Figure 4A. The re-weighted PMF allows us to identify four low-energy states: the unbound (”U”), intermediate 1 (”I”), intermediate 2 (”II”) and bound (”B”) states. The maximum energy was capped at 15 kcal/mol for clarity in identifying relevant states. The unbound state is located in a local energy well centred at ∼(10, 60 Å), the intermediate 1 state is centred at ∼(50, 40 Å), the 2nd intermediate state is centred at ∼(20, 20 Å) and the bound state is located at the global energy minimum ∼(60, 0 Å). We should note here that propofol binding within the M2 helix pore was observed in all three 500 ns GaMD simulations, so we can assume that calculated energies between individual states have a low degree of error due to limited sampling. Propofol forms hydrogen-bonds within the channel pore with residues Thr-236, Thr-1480, Thr-858 and Thr-1169 (Figure 4C). Hydrogen-bonds are also formed along the pathway to the pore with residues Ile-854, Arg-1345, Asn-1319 and Glu-1305.

Anesthetic transmembrane binding site

Most studies conducted on this system focus on the binding site which was observed in the original crystal structure, (Nury et al., 2011) for example using docking methods. To assess whether any propofol molecules in our GaMD simulations entered this site, we defined the binding site by the following residues: Tyr-193, Ile-198, Met-201, Val-238, Thr-251 and Ile-255. We used the center of mass (COM) of these binding site residues and computed the distance between this COM, and any propofol molecules which were identified as being close to the subunit by visual inspection in VMD. (Humphrey et al., 1996) During the calculations, a step size of 10 was used. Figure 5 shows the distances of the identified propofol molecules to the anesthetic transmembrane binding site. No molecules were identified for S2. Figure 6 shows multiple molecules competing at S1, although we notice that none of the molecules get close to the binding site. Stable regions are located between 10 and 20 Å from the experimentally observed site. We initially suspected that this distance was due to the propofol molecules interacting with each other and forming a small cluster, as was observed in our previous work. (Faulkner et al., 2020) Upon further investigation we observe this same distance for subunits S3 and S4 in which only one molecule has stable distances between 10 and 20 Å, which rules out the hypothesise that this is caused by intermolecular interaction. The data shown in Figure 5 were taken from one replicate in which the most interactions were observed. When we investigate the site at which the propofol molecules are stable, we observe a site formed by residues Phe-261, Phe-295 and Val-264 (Figure 6B). These residues block the passage of the propofol molecules into the upper region of the transmembrane domain. The Phe residues form π-π interactions with each other and the propofol molecule, which contributes to the stability observed here. If we are to take the crystal structure binding site as true, we must assume that longer time scales are required to observe binding, or an allosteric mechanism in which another binding site has to be occupied to allow passage of propofol through the Phe-261 Phe-295 site into the crystal binding site, although none of this was observed in 3 μs of cMD and 1.5 μs of GaMD. We should note here that no propofol molecules were observed to enter from the other side of the M4 helix (between M4 and M1). To investigate the dynamical consequences on GLIC induced by the drug molecules, we computed the distance between Asp-32 and Arg-192 which form a salt bridge. It has been suggested that salt bridge perturbation is important for channel functions; this was shown in a nAChR channel, (Xiu et al., 2005) and a GLIC where the volatile anesthetic halothane caused disruption of this salt bridge, which led to instability within the channel. (Chen et al., 2010) We did observe for multiple subunits, where drug molecules were present, that this salt bridge did not hold.

Figure 5. Calculated distances between the binding site COM and identified propofol molecules. Each line in the graphs represents a propofol molecule.

Figure 6. Interaction site of propofol at GLIC subunits: (A) Sites of propofol interaction in relation to the anesthetic binding site in which residues are coloured; Tyr (green), Ile (blue), Met (red), Val (orange), Thr (pink) and His (yellow). Propofol molecules are shown in licorice representation. (B) Blocking residues inhibiting the path to the anesthetic binding site, propofol (violet) shown in licorice representation.

To determine whether fentanyl and propofol share binding sites, we constructed density maps using the VOLMAP VMD plugin. (Humphrey et al., 1996) These density maps, shown in Figure 7, were averaged over the whole simulation trajectories, apart from the first 50 ns which were discarded for system equilibration. Identified regions of high density therefore reflect locations, where drug molecules resided for a prolonged period of time. Density maps were computed for both the cMD flooding simulations and the GaMD simulations, which allowed us to determine if the accelerated simulations were able to discover binding sites which are unattainable from conventional simulations. The main difference that we observed between the classical molecular dynamics flooding simulations and the GaMD simulations is the occupation of the ion channel pore by propofol, which was observed in all 3 GaMD simulations. The pore is occupied by one propofol molecule which resides above the hydrophobic gating Ile-239 residues that form the top of the hydrophobic gate in GLIC (Zhu & Hummer, 2010). The molecule is mobile within the pore with a RMSD of 6.06 ± 0.76 Å. This propofol molecule blocks the majority of the pore, which disrupts the flow of water molecules through the pore. As stated, the pore-blocking molecule is mobile within the pore, where it mostly forms hydrogen-bonds with the threonine residues located at the top of the M2 helices. Hydrogen-bonding calculations show that these residues form the majority of the hydrogen-bonding interactions within the pore (Figure 4C). Only one propofol molecule was found to block the pore in each GaMD simulation. The general anesthetic isoflurane has been shown to block the pore of GLIC and nAChR as a dimer, (Brannigan et al., 2010) but this behaviour was not observed for propofol in any of our simulations, most likely due to the larger size of propofol compared to isoflurane. LeBard et al. (LeBard et al., 2012) predicted the physical blocking of the pore in GLIC by propofol, which was shown to bind to the pore with high affinity, as calculated from free energy perturbation (FEP) calculations. They predicted that propofol would bind as a dimer in a similar fashion to isoflurane, but we did not observe this in around 5 μs of conventional and Gaussian accelerated MD simulations. A structure of propofol bound to the pore of GLIC has not been resolved experimentally, most likely due to the binding of detergents (dodecylmaltoside, DDM) within the pore. (Bocquet et al., 2009; LeBard et al., 2012)

Figure 7. Regions of prolonged occupation by propofol (red) and fentanyl (blue) in the GLIC. (A and B) represent the TM and EC domains of the GLIC, respectively, from the GaMD simulations. (C and D) represent the TM and EC domains of the GLIC, respectively, from the cMD flooding simulations.

Our GaMD simulations allowed us to identify the pathway by which propofol binds to the channel pore, Figure 8 shows the sites associated with the intermediate states calculated from the reweighted GaMD simulations. The first intermediate state of the propofol binding pathway is located in the β3 loop region at the top of the extracellular domain, where propofol forms short-lived hydrogen-bonds to the glutamic acid residues before moving to the second intermediate state. This second binding state is located between the β1-β2 loop and the β4-β5 loop in the inner region of the extracellular vestibule, where propofol remains for ∼36 ns before diffusing into the M2 helix pore to remain there for the rest of the simulation. Propofol molecules that were observed to bind near the transmembrane intra-subunit anesthetic binding site displayed interesting dynamics within the binding site during the GaMD simulations, where the molecule would” hop” from one intra-subunit site to the adjacent intra-subunit site. Propofol was observed to traverse through a binding tunnel between the M3 helix of one subunit and the M1 helix of the adjacent subunit into an inter-subunit site initially, then into the intra-subunit site (Figure 9). We also observed multiple occupancy in transmembrane sites where up to three propofol molecules occupied an intra-subunit site. This multiple occupancy was initiated by small clusters of propofol molecules forming within the lipid bilayer region and then diffusing towards the TMD where they then bind to the TMD helices. Fentanyl was not observed to bind within the channel pore; We have previously shown that fentanyl binds to an inter-subunit site in the extracellular domain, (Faulkner et al., 2019) when propofol is not present. In this study we observe sites where fentanyl shares a binding site with propofol, mostly shown in our flooding simulations in the extra-cellular domain (Figure 7D). In the extra-cellular domain we initially observe several propofol molecules interacting with a single fentanyl molecule, before this cluster moves towards the extra-cellular region of the ion channel, where several propofol molecules diffuse away from the fentanyl molecule and the remaining propofol molecule forms a very stable hydrogen-bond with Asp-1418 located on loop C (Figure 10C). This hydrogen-bond remains stable for the majority of the simulations, where the fentanyl molecule remains located beside the propofol molecule, which we believe acts as a stabiliser for the bound propofol which has an RMSD of 2.04 ± 0.4 Å. Table 1 shows the RMSD values for propofol when bound to Asp-1418, before and after the interaction with fentanyl. We can see that there is a stabilisation effect on propofol when fentanyl is present at the binding site. When we look at the flexibility of loop C, where propofol binds (Figure 10C), we see that the flexibility is almost identical to the pure structure when fentanyl is present. When only propofol binds here, we see that the flexibility of this region increases significantly. This shows that fentanyl not only stabilises the propofol molecule, but also a structurally important section of the GLIC extracellular region. We also analysed the flexibility of the cys-loop (Figure 10A) in the pure system and compared it to the systems in which propofol and fentanyl are bound. The cys-loop is thought to be a crucial component in the gating of pLGICs, (Cederholm et al., 2009; Nys et al., 2013; Schofield et al., 2003) and in our analysis we can clearly see the increased mobility of the cys-loop (Figure 10B) when the drug molecules are bound. This shows that the combination of these drug molecules could play a role in modulating the channel function. The backbone of the fentanyl molecule remains flexible throughout the simulations, but the N-phenyl-propanamide ring remains at an ∼90^° angle to the aromatic ring of propofol, forming a proposed π-π interaction to stabilise the propofol in this site. From the density maps (Figure 7) we can see that there is overlap of the two drugs in the transmembrane domain, as well as in the extracellular domain. Within the TM domain we observe similar interactions as those shown in Figure 10C, where fentanyl forms π-π interactions with propofol and holds it in a specific site. The transmembrane sites where we observe these interactions are located near the bottom of the M2 helices, close to the phospholipid headgroup region. The fentanyl molecules stabilise the propofol molecule in the lower M2 intra-subunit site before diffusing to the next subunit. This behaviour was observed in both the flooding style simulations and the GaMD simulations.

Figure 8. Graphical representations of the intermediate states of propofol (green sticks) calculated from reweighted GaMD simulations. Contact residues shown explicitly in stick representation.

Figure 9. (A) Graphical representation of the propofol subunit hopping pathway in the TMD. (B) RMSD plot showing the change in position of propofol from the initial site (1) to the intermediate, intersubunit site (I) and then to the adjacent intrasubunit site (2). Reference for the RMSD calculation was taken from propofol in site 1. Data shown from one trajectory.

Figure 10. (A) Structure of one subunit highlighting key areas of the protein. (B) Cys-loop RMSD plot for pure and drug-bound systems. (C) π-π interaction of the fentanyl N-phenyl-propanamide ring and the propofol aromatic ring and the hydrogen-bond between propofol and Asp-1418 on loop C.

Table 1. RMSD values for propofol binding to Asp-1418 with and without fentanyl, and loop C which Asp-1418 is part of.

Region/molecule	RMSD (Å)
Propofol + Fentanyl^a	2.04 ± 0.40
Propofol^a	2.84 ± 0.44
Loop C (Pure)^b	2.03 ± 0.17
Loop C (Propofol + Fentanyl)^a	1.97 ± 0.36
Loop C (Propofol)^a	3.41 ± 0.32

^a Data averaged over flooding and GaMD trajectories where interaction was observed.

^b Data obtained from pure GLIC trajectory.

4. Conclusions

In this study we have applied flooding style and Gaussian accelerated molecular dynamics simulations to study how the general anesthetic propofol and the opioid analgesic/anesthetic fentanyl interact with the GLIC channel. Our GaMD simulations reveal a detailed pathway for the pore-blocking mechanism by propofol, which has not been observed previously using conventional MD simulation methodologies. The GaMD simulations also showed that propofol is able to hop between transmembrane intra-subunit binding sites by diffusing through a binding tunnel formed by helices of the occupied and adjacent transmembrane subunits. Fentanyl was shown to play an important role in the stabilisation of propofol molecules in several binding sites. Initial interactions between fentanyl and propofol lead to the accumulation of a small number of propofol molecules with a fentanyl molecule. This small cluster either diffuses rapidly into the membrane region where most propofol molecules leave the fentanyl, with one remaining propofol binding to the transmembrane domain of GLIC where fentanyl acts as a stabiliser. The other propofol molecules, diffuse towards the extracellular domain of GLIC, binding close to the orthosteric agonist-binding site, where we again observe ligand stabilisation by fentanyl. The results presented here show the first evidence of the opioid fentanyl interacting with, and altering, propofol binding in GLIC. The structural results of propofol and fentanyl interacting at GLIC present a compelling case for the allosteric action of anesthetics and opioids at pLGICs. We consider that this study, along with our previous work, (Faulkner et al., 2019) calls for a comprehensive examination of the role of the opioid in the general anesthesia process.

Acknowledgements

This research was undertaken using the Supercomputing Facilities at Cardiff University operated by ARCCA on behalf of CSF and Supercomputing Wales. This equipment is part-funded by the Welsh European Funding Office (WEFO) as well as by Cardiff University. Information on the data underpinning the results presented here, including how to access them, can be found in the Cardiff University data catalogue at http://doi.org/10.17035/d.2020.0106797733

Disclosure statement

There are no conflicts to declare

Additional information

Funding

We acknowledge the Engineering and Physical Sciences Research Council (Grant No. EP/R512503/1) and AWE for funding.