Atomistic Molecular Dynamics Simulations of Propofol and Fentanyl in Phosphatidylcholine Lipid Bilayers.

Atomistic molecular dynamics (MD) and steered MD simulations in combination with umbrella sampling methodology were utilized to study the general anesthetic propofol and the opioid analgesic fentanyl and their interaction with lipid bilayers, which is not yet fully understood. These molecules were inserted into two different fully hydrated phospholipid bilayers, namely, dioleoylphosphatidylcholine (DOPC) and dipalmitoylphosphatidylcholine (DPPC), to investigate the effects that these drugs have on the bilayer. We determined the role of the lipid chain length and saturation on the behavior of the two drugs. Pure, fully hydrated DOPC and DPPC bilayers were also simulated, and the results were in excellent agreement with the experimental values. Various structural and mechanical properties of each system, such as the area per lipid, area compressibility modulus, order parameter, lateral lipid diffusion, hydrogen bonds, and radial distribution functions, have been calculated to assess how the drug molecules affect the different bilayers. From the calculated results, we show that fentanyl and propofol generally follow similar trends in each bilayer but adopt different favorable positions close to the headgroup/chain interface at the carbonyl groups. Propofol was shown to selectively form hydrogen bonds at the carbonyl carbon in each bilayer, whereas fentanyl interacts with water molecules at the headgroup interface. From the calculated free-energy profiles, we determined that both molecules show a preference for the low-density, low-order acyl chain region of the bilayers and both significantly preferred the DOPC bilayer with propofol and fentanyl having energy minima at -6.66 and -43.07 kcal mol-1, respectively. This study suggests that different chain lengths and levels of saturation directly affect the properties of these two important molecules, which are seen to work together to control anesthesia in surgical applications.

Introduction

An important stage of the mechanism of action of any drug molecule is the diffusion through cellular membranes. The favored positions of these drug molecules within the target membranes also affect their transport and metabolism.[1] The most important membrane encountered by drug molecules is the plasma membrane into which these pharmaceuticals penetrate and with which they interact to reach the target cells.[2] The major components within these plasma membranes are the phosphatidylcholine (PC) lipids,[3] which are hence most commonly used in the modeling of drug interactions with cell membranes. Knowledge of the interactions between drug molecules and lipids is clearly essential to improve their mechanism of action and efficiency.

In this study, we have studied two drug molecules embedded in two different PC bilayers, namely, dipalmitoylphosphatidylcholine (DPPC) and dioleoylphosphatidylcholine (DOPC), which are displayed in Figure c,d, respectively. Both lipids consist of a zwitterionic headgroup region with a positively charged choline group and a negatively charged phosphate group. The difference between the lipids lies in the different lengths and saturation of the linear chains, (i.e., 16, fully saturated for DPPC and 18, unsaturated for DOPC).

Figure 1

Molecular structures of (a) propofol, (b) fentanyl, (c) dipalmitoylphosphatidylcholine, and (d) dioleoylphosphatidylcholine. The first and second carbons of the glycerol fragment are denoted as the sn-1 and sn-2 positions, respectively.

Since its introduction, propofol has been used routinely in general surgery as a general anesthetic to induce and maintain anesthesia. Propofol is also commonly used alongside opioid analgesics such as fentanyl. The use of these opioid analgesics has been shown to potentiate propofol,[4] thereby reducing the amount of drug required to induce anesthesia, which in turn reduces the unwanted side effects of propofol anesthesia, e.g., a large drop in blood pressure[5] and apnea.[6] Propofol is proposed to interact and bind at the γ-aminobutyric acid (GABA) receptors,[7] consisting of ligand-gated ion channels (GABAA) and G protein-coupled receptors (GABAB), which both control the flow of ions into the intracellular domain. There have also been studies that suggest that propofol acts at G protein-coupled receptors by inhibiting the function of M1 muscarinic acetylcholine receptors.[8] These studies suggest possible binding sites of propofol, but a complete mechanism for the action of anesthesia still remains unknown. Due to its properties, fentanyl is used in general surgery as the analgesic component to general anesthesia and also as the lone component.[9] Various studies have shown that when used in combination with propofol, it provides a safer and more satisfactory anesthesia method.[10−12] Fentanyl is known to act at various opioid sites throughout the body but predominately at the μ-receptor, although it also binds to δ-type and κ opioid receptors.[13] Although the binding sites for opioid drugs are relatively well known, it still remains unclear why fentanyl is able to potentiate propofol during anesthesia. Our previous work showed that fentanyl also modulates the Gloebacter violaceus ion channel, which is a known site for general anesthetic action.[14]

Propofol, as shown in Figure a, consists of a benzene ring core with two isopropyl groups at the 2 and 6 positions and one hydroxyl group in between. The hydroxyl group is reasonably well guarded by the steric bulk of the isopropyl groups, so hydrogen bonding will only be available in certain orientations. The structure of fentanyl, shown in Figure b, differs significantly from propofol. Fentanyl has an almost linear, flexible backbone consisting of two aromatic and one aliphatic six-membered rings and an acetyl group bonded to the terminal nitrogen.

To provide atomistic insight into the interaction of these drug molecules with the cellular membrane, molecular dynamics (MD) simulations have been used to investigate the thermodynamics and structural changes that these drug molecules induce on the lipid bilayer. Despite the enormous growth of the field of biomolecular simulation in recent years,[15] to the best of our knowledge, no studies have been reported, which directly compare, in a single study, these two drug molecules interacting with phospholipid membranes. Therefore, the primary aim of this paper is to investigate how these drug molecules behave within the bilayer systems by characterizing how the anesthetics alter the structural properties of the membranes themselves. In addition, we have also calculated the free-energy profile of diffusion of each drug molecule through each bilayer system, using the biased MD umbrella sampling method.

Results and Discussion

Area per Lipid

The area per lipid is considered to be one of the most important properties to describe the behavior of the bilayer and whether it is in the correct, biologically relevant Lα phase. This area can be calculated easily from our MD trajectories. The calculation involves dividing the xy cross-sectional area of the orthorhombic periodic cell (i.e., the lateral area of the bilayer) by the number of lipids.[16] The areas per lipid calculated for each pure bilayer simulation are shown in Table . The results obtained were within around 3–4% of the experimental values from the literature, which suggests that our simulated pure bilayers are in the correct phase and have equilibrated sufficiently to be used for the drug molecule simulations. The results for bilayers containing propofol are also shown in Table . Here, we see an increase in the area per lipid when propofol is simulated in the center of each bilayer (0.25 Å in DOPC and 0.27 Å in DPPC) due to the drug molecule interacting in the upper chain region around the carbonyl groups at the headgroup–chain interface. When propofol was simulated in the water phase, we observed natural diffusion into both types of lipid bilayers, leading to a larger increase in area due to disruption in the headgroup regions. This expansion of AL can be rationalized in terms of the lone-pair repulsions as propofol diffuses through the headgroup region. A similar pattern is observed when fentanyl is in the center of the bilayer as we observe further increases in the area, especially for DPPC. Our simulations of DPPC indicate that fentanyl does not lie as close to the headgroup region as propofol does but rather lies parallel and slightly below the carbonyl groups. We only observed spontaneous diffusion into the DOPC system for fentanyl, which caused a small 0.1 Å increase in the area per lipid. This observation can be explained by fentanyl diffusing in a linear conformation, which causes less disruption, whereas propofol is much more dynamic in its diffusion process, adopting more conformations. The calculated volume per lipid followed the same trend as the area per lipid as they are directly proportional.

Table 1

Averaged Structural Properties of Pure and Drug-Containing Lipid Bilayer Systems, Where Area per Lipid = (AL), Volume per Lipid = (VL), Isothermal Area Compressibility Modulus = (KA), Bilayer Thickness = (DHH), and Lateral Lipid Diffusion Coefficient = (D)a

lipid	drug	AL (Å2)	VL (Å3)	KA (m Nm–1)	DHH (Å)	Dxy (10–8 cm2 s–1)
DOPC experimental	none	67.4, 72.5[17]	1303[18]	265, 318[18]	35.3, 37.1[17]	11.5, 17[19]
DOPC calculated	none	70.1 ± 0.2	1274.2 ± 1.1	285.2 ± 19.7	37 ± 0.2	6.30
DPPC experimental	none	63.1, 64.3[20]	1232[18]	231[21]	38, 38.3[17]	12.5, 15.2[22]
DPPC calculated	none	62.9 ± 0.3	1175 ± 0.2	231.9 ± 22.7	37.8 ± 0.2	8.8
DOPC	propofol center	70.3 ± 0.3	1275 ± 1.5	419.5 ± 15.6	36.25 ± 1.1	12.75
	propofol water	70.94 ± 0.2	1279.7 ± 1.0	333 ± 30.1	36.75 ± 0.5	11.25
	fentanyl center	70.33 ± 0.2	1261.7 ± 1.3	227.6 ± 16.0	36.5 ± 1.0	7.9
	fentanyl water	70.2 ± 0.4	1262.3 ± 1.3	274 ± 20.5	37.25 ± 0.2	6.4
DPPC	propofol center	63.17 ± 0.4	1176.8 ± 1.5	319.5 ± 22.7	37 ± 0.7	17.5
	propofol water	63.50 ± 0.5	1177.8 ± 0.9	243.1 ± 16.0	37.25 ± 0.9	14
	fentanyl center	63.58 ± 0.4	1179 ± 1.1	236.2 ± 25.0	37.5 ± 0.8	8.3
	fentanyl water	62.86 ± 0.6	1179 ± 1.0	331.2 ± 13.0	37.25 ± 0.70	6.4

Water and center refer to the starting positions of the drug molecule.

Isothermal Area Compressibility Modulus

The isothermal area compressibility modulus (KA) is the stress required to induce an isotropic expansion in volume, which can be calculated from MD simulations aswhere kB is Boltzmann’s constant, T is the simulation temperature, ⟨AL⟩ is the average area per lipid, σA2 is the variance in the area per lipid over the simulation, and nlipid is the number of lipids in the simulation box. KA is a standard descriptor of the membrane phase, with large values being characteristic of the Lβ gel phase, as the chains would be fully extended and their van der Waals interactions would strengthen and lead to tighter, more ordered lipid packing.

Our results for the KA of the pure bilayers are in excellent agreement with the experimentally obtained values. When propofol starts in a position in the center of either membrane and diffuses toward the carbonyl region, we found an increase in KA, which induces a closer packing of the lipid chains. The conformational changes of the membrane chains increase the aforementioned interactions between them, causing the bilayer to stiffen (see Figures and 13). The largest increase in KA is observed for DOPC, owing to the presence of the double bond in the chain, in which the drug molecule inhibits the fluidity by binding in the upper chain region. During the diffusion process from the water phase into the membrane, we observe a small increase in KA due to the fast diffusion through the headgroup region. This observation suggests that the process of stiffening is not instantaneous, and the direction of diffusion plays a key role. In the DOPC bilayer, we calculate a lower KA for fentanyl compared to that for propofol, which indicates that fentanyl is causing more fluctuation within the headgroup region of the bilayer (see Figure ). For the natural diffusion of fentanyl into the DOPC bilayer, we observe that KA decreases by a small amount (11.2 m Nm–1) compared to that of the pure membrane, which shows that the diffusion does not allow the chains to pack particularly tightly together. For fentanyl in the DPPC bilayer, we calculate a small 4.3 m Nm–1 increase compared to that of the pure membrane, even though the drug molecule positions itself parallel to the headgroups at the interface, as it does in DOPC. As such, the difference in the chain structure plays a significant role in the way that the structural properties of lipid bilayers are affected by drug molecules. The most unexpected result found for fentanyl is the simulation with the drug molecule in the water phase where we saw no diffusion, even though the resulting KA shows an increase of 99.3 m Nm–1 over the pure bilayer. This result is interesting because fentanyl has no direct contact with the lipid chains and spends most of its time during the simulation at the water–headgroup interface. During the course of the simulation, we observe several partial diffusions, where fentanyl enters the headgroup region for 1–4 ns before returning into the water phase. These partial diffusions suggest that fentanyl induces a conformational change in the headgroups, which allows the chains to pack closer together and stiffen the bilayer. To confirm whether there was a change in the conformation of the headgroup environment when fentanyl partially diffuses into the headgroup region of the DPPC bilayer, the angles between the P–N vector and the normal of the bilayer were calculated using the MEMBPLUGIN tool[23] for VMD.[24] For the pure DPPC bilayer, an angle of 69.0 ± 2.4° was calculated, which is almost identical in the fentanyl-containing system. To study the specific conformation of the headgroups, which fentanyl causes disruptions when it partially diffuses, headgroup molecules were selected, which were 15 Å from the fentanyl molecule. The P–N vectors were calculated as 73.0 ± 2.0°, which shows that fentanyl causes a conformational change with the headgroup region. Results for propofol are not shown as there was no change in the P–N vector due to the fast diffusion through the headgroup region into the membrane interior.

Figure 4

SCD order parameters for (a) DOPC and (b) DPPC compared to experiment.[30−32]

Figure 13

SCD order parameters for DOPC (top) and DPPC (bottom) containing clinical concentrations of propofol and fentanyl.

Figure 2

Starting positions for fentanyl (A) and propofol (B) in the bilayers. End positions after the production simulation of fentanyl (C) and propofol (D); DOPCs shown as positions were similar in both systems. Headgroups are shown as yellow sticks, headgroup phosphate is shown as an orange sphere, and lipid chains are shown as gray sticks.

Lipid Lateral Diffusion

The lateral diffusion of a lipid bilayer is an important dynamical property, as it can affect many membrane parameters. For example, one study reported a connection between the lateral lipid diffusion and the viscosities in different parts of the membrane.[25] In this work, we have calculated the lateral diffusion coefficient of the lipids without (for reference) and with the drug molecules, using the mean-squared displacement (MSD) of the membrane, in the xy direction, over 20 ns window lengths and averaged over time origins separated by 200 ps (Figure ). This property is related to the Einstein equation in two dimensionswhere ⟨|Δr(t)|2⟩ is the MSD in the XY plane in time t, and D represents the lateral diffusion coefficient. The diffusion coefficient can be obtained from the gradient of the linear portion of the MSD plot. Table lists the calculated lateral diffusion coefficients for pure DOPC and DPPC, which are comparable to those calculated in the Lipid14 paper.[26] Our production runs were performed in the NPT ensemble using Langevin dynamics to control the temperature, which randomizes particle velocities and can therefore affect dynamic membrane properties such as lateral diffusion. Slightly more accurate results can be obtained using the microcanonical (NVE) ensemble, which was applied in the original lipid14 publication.[26] However, as our results for the pure bilayers were within the experimentally determined values, we decided to continue the drug molecule simulations in the NPT ensemble. With the addition of propofol within the bilayer, we observe a slight increase in the lateral diffusion coefficient for both lipid systems (Table ). This finding is unexpected since the drug molecules can penetrate the lipid bilayer where they occupy a portion of the free volume within the membrane and reduce the lateral diffusion.[27] This hypothesis is indeed observed when we have a spontaneous diffusion of propofol from the water phase as it occupies the free volume space within the headgroup region, whereas propofol starting within the bilayer does not penetrate fully into the headgroup space. For fentanyl, we observe a significant decrease in the lateral diffusion constant in both systems with respect to propofol. This difference in behavior can be explained by the larger size of the fentanyl molecule compared to that of propofol, which therefore occupies more free volume within the bilayer, hence limiting the fluidity of the membrane.

Figure 3

Time-averaged mean-squared displacement of the lipid center of mass (COM) versus simulation time.

Lipid Tail Order Parameters

The hydrocarbon tails that make up the hydrophobic interior of the bilayer are highly dynamic, which accounts for the overall fluidity of the membrane.[28] The mobility of these chains at individual carbon positions can be evaluated by measuring the order parameter. A value of SCD = 1 implies a position parallel to the bilayer normal, and a value of SCD = 0 implies total random motion. Experimentally, this property is determined by using deuterium NMR quadrupole splitting, which entails substituting the hydrogen atoms at each carbon position with deuterium and measuring their dynamic variations by 2H NMR. A detailed description of the procedure can be found here.[29] This property can also be calculated from our simulations, where the carbon–deuterium order parameter is determined by the tensor S in the equationwhere θ represents the angle formed between the Z direction and the bilayer normal. Simulated results for the pure bilayers are shown in Figure , with comparison to literature values. Our results obtained for the pure bilayer simulations are in good agreement with the trends observed experimentally. Figure shows the results obtained for the bilayers containing the drug molecules. For the DOPC sn-1 chain, we observe very little change when the drug molecules are present compared to the pure bilayer. A very small increase is observed between carbons 3 and 6, with respect to the pure bilayer, which suggests a slightly restricted motion of the chains at this position due to the presence of the two drug molecules. A small decrease in the order parameter for fentanyl between carbons 10 and 11 suggests more mobility at these positions. This behavior is also seen in the sn-2 chain for both anesthetic molecules, although note the lower values in carbons 12–14, which could account for the reduced area compressibility calculated for fentanyl due to the chains becoming more flexible and hence requiring less force to induce an isotropic expansion.

Figure 5

SCD order parameters for (a) DOPC sn-1 chain, (b) DOPC sn-2 chain, (c) DPPC sn-1 chain, and (d) DPPC sn-2 chain. Order parameters for drug-containing systems calculated when the drug was present in the hydrophobic phase of the bilayer.

In the DPPC system, we see clearer differences between the drug-loaded bilayers compared to that between the pure bilayer reference. In the sn-1 chain, we observe that propofol increases the order parameter from carbons 5 to 12, reducing the chain mobility and hence requiring more force to induce an isotropic expansion within the bilayer, in agreement with KA. The introduction of fentanyl causes only a negligible change in the order parameter of the pure DPPC bilayer. The sn-2 chain shows a clear difference in the 6–13 region compared to the pure membrane. Propofol causes a higher-order parameter for these carbons, as seen in the sn-2 chain, whereas fentanyl shows lower values at this point, which would suggest a more flexible chain and hence a lower KA. Table shows that there is no correlation between the lipid chain mobility and KA, which indicates that the flexibility in the headgroup region could be the link between them. We should note here that a similar study conducted on the local anesthetic articaine in a DMPC bilayer found higher-order parameter values for the neutral form of the molecule, similar to propofol in the DPPC bilayer.[33] The differences observed in that study are larger, which is most likely due to the higher concentration of drug molecules used in the simulations.

Hydrogen Bonds

Hydrogen bonds were observed between the drug molecules and water, and propofol was also observed to form hydrogen bonds with the headgroups of both lipid systems. The conditions used to determine if a hydrogen bond was formed were set by a bond distance of 2.5 Å between the donor hydrogen and the acceptor atom with a maximum angle of 30° between the donor hydrogen and acceptor hydrogen vectors.[34]Figure a shows that the pattern of hydrogen bonding to water is similar in both bilayer systems. Hydrogen bonds were mostly formed between the hydrogen of the propofol hydroxyl group and the oxygen of the water molecules. A few times during the simulations, the conformation allowed the hydrogen from another water molecule to coordinate with the propofol hydroxyl oxygen. The average number of hydrogen bonds formed per frame was similar for DOPC (0.120) and DPPC (0.119). Figure b shows the hydrogen-bonding pattern for fentanyl in both bilayers. Hydrogen bonds are formed faster within the DPPC bilayer as this is the most fluid membrane (Table , D), which allows the drug molecule to move faster toward the interface with water and hence form hydrogen bonds. Fentanyl has two extra possible hydrogen-bonding sites, i.e., carbonyl oxygen and both nitrogens, compared to propofol that has only one (hydroxyl oxygen) site, which explains the different number of hydrogen bonds formed. Fentanyl forms bonds between both the carbonyl oxygen and the neighboring nitrogen to water hydrogen. No bonds are formed to the piperidine nitrogen due to steric hindrance and the conformation adopted by fentanyl. Again, we calculate very similar hydrogen bonds in each bilayer, i.e., 0.758 and 0.715 for DOPC and DPPC, respectively. From Figure c, we can see that propofol forms hydrogen bonds to the lipid headgroups in both DOPC and DPPC, which fentanyl is unable to do. The carbonyl oxygen of the sn-1 chain is favored in both lipids over the sn-2, due to the conformation that the headgroup adopts throughout the simulations, which makes it better accessible for bonding. Figure d shows a snapshot from the simulation of DOPC, where we noted a number of water molecules entering the hydrophobic phase and forming “anchoring” hydrogen bonds to fentanyl. This behavior was observed in both bilayers but to a greater extent in the DOPC membrane. This phenomenon was not seen with propofol, which could account for its higher mobility within the bilayer compared to that of fentanyl. Table shows the average number of hydrogen bonds per frame between lipid headgroups and water molecules, which is used as a means to assess the effect of the drug molecules on the ability of the headgroups to form hydrogen bonds to the surrounding media. In the DPPC bilayer containing either of the two drug molecules, we observe very little change in the number of hydrogen bonds between the headgroups and water compared to the pure bilayer. The results for the DOPC system show a clear increase in the number of hydrogen bonds formed in the presence of fentanyl. This observation can be explained by the disruptions caused by fentanyl in the headgroup region (Figure ), which can accommodate more water molecules at the interface, thereby increasing the potential for hydrogen bonding. These results highlight the different behavior observed in the two different PC lipid bilayers despite their close structural similarity.

Figure 6

Hydrogen-bond plots of (a) propofol-to-water, (b) fentanyl-to-water, and (c) propofol-to-lipid headgroups; and (d) snapshot of fentanyl binding to water within the hydrophobic phase. Hydrogen-bond distance is shown in panel in Angstroms (d).

Table 2

Number of Hydrogen Bonds Formed between Lipid Headgroups and Water

lipid	drug	average bond per frame
DOPC	none	464.880
	fentanyl	481.300
	propofol	467.752
DPPC	none	396.409
	fentanyl	397.453
	propofol	399.372

Figure 7

Multiple water molecules moving through the DOPC lipid headgroup into the hydrophobic phase of the bilayer. Lipid headgroup is shown as yellow sticks, lipid tails are shown as gray sticks, phosphorous atom of the phosphate group is shown as an orange sphere, water molecules within 5 Å of fentanyl are shown as blue spheres, and fentanyl is shown as sticks.

Radial Distribution Function

To gain further insight into the distribution of water and the drug molecules around the bilayer, we have carried out radial distribution function (RDF) analyses. Figure shows the RDF plots for the oxygen atom of water around the phosphorous of the headgroup. From the insets of each graph, we can see that when fentanyl is added to the membrane there is a higher density of water molecules around the phosphate group. This difference is slightly larger in the DOPC system, which agrees with the observations that fentanyl causes more water molecules to penetrate into the headgroup region. RDFs were also computed for the atoms in each drug molecule, which exhibited hydrogen bonding to water and, in the case of propofol, the lipid headgroup. Figure a displays that water is more concentrated around the carbonyl oxygen of fentanyl in DOPC, with an average hydrogen-bond distance to water of 1.75 Å. The peak for the oxygen-to-oxygen density is a further 1 Å away, which indicates the length of the O–H bond in the water molecule. Figure b shows the same trend for DPPC as for DOPC, although the densities for N to H and N to O are closer, with the O-to-O density slightly higher by 0.05. Figure c shows the dominance of hydrogen bonding between propofol and the headgroup region. The conformation of the headgroup in DOPC increases the sn-1 lipid chain in such a way that it is better accessible to form hydrogen bonds to propofol. We found a higher density for sn-1 and sn-2 in DPPC due to the differences in chain mobility between the two systems; see Figure d. A sharp shoulder peak is observed at 1.95 Å in the propofol DPPC system, which is due to hydrogen bonds forming at slightly different distances.

Figure 8

Radial distribution function plots of (a) water oxygen to phosphorous of DOPC and (b) DPPC lipid bilayers. The insets show the increase in density caused by the addition of the drug molecules.

Figure 9

Radial distribution function plots of (a) fentanyl hydrogen-bonding atoms to water hydrogen and oxygen in DOPC; (b) fentanyl hydrogen-bonding atoms to water hydrogen and oxygen in DPPC; (c) propofol hydroxyl hydrogen to water hydrogen and oxygen, and sn-1 and sn-2 carbonyl oxygen in DOPC; and (d) propofol hydroxyl hydrogen to water hydrogen and oxygen, and sn-1 and sn-2 carbonyl oxygen in DPPC.

Clinical Concentration Simulations

By conducting these simulations at clinical concentrations, we can observe how these drug molecules perturb the membrane structure in a way that we cannot observe in the single-molecule simulations. Initially, the molecules were added to the water phase outside of the membrane, and over the course of the simulations, they diffuse into the membrane interior. Figure shows the final positions of the drugs within the membrane, which were very similar in DOPC and DPPC, so only one is shown. The average positions of the drug molecules over the course of the simulations can be calculated with decomposed electron density profiles computed using CPPTRAJ.[35] These profiles (Figure ) were constructed for each drug molecule in each bilayer using an average over all replicates. Both drug molecules are located predominantly under the lipid headgroups in the membrane interior. This observation is consistent with other molecules that possess anesthetic properties, such as alcohols,[36] benzocaine,[37] and halothane.[38] The two drug molecules prefer slightly different depths within the bilayers, with propofol having a density maximum at approximately 10.3 Å from the bilayer center in both DOPC and DPPC, while fentanyl has a density maximum at approximately 9.1 Å in DPPC and approximately 10.3 Å in DOPC. The density for fentanyl in both bilayers is higher than that for propofol in the center (0 Å), which indicates that the fentanyl molecules are able to cross the bilayer on the time scale of the simulation, as the free energy is lower in the center for fentanyl (Figure ).

Figure 10

End states of propofol (left) and fentanyl (right) in the lipid bilayers. Phosphorous atoms are shown as orange spheres, lipid headgroups are shown as yellow sticks, lipid tails are shown as gray sticks, and drug molecules are shown as VdW spheres.

Figure 11

Total and decomposed electron density plots for propofol in DOPC (A) and DPPC (B), and fentanyl in DOPC (C) and DPPC (D). Contributions from water, choline, phosphate (PO4), glycerol (Gly), carbonyl (COO), methylene (CH2), unsaturated CH = CH, and terminal methyls (CH3).

Figure 16

Potential of mean force profile for fentanyl in DOPC and DPPC bilayers.

To gain further insight into the hydration state and coordination of the drug molecules and the probability distribution of the water molecules around the drugs and the lipid headgroups, RDFs of the water oxygen around the phosphate of the lipid headgroup have been calculated for both drug molecules in both lipid bilayers, as shown in Figure . The RDF plots for the water oxygen around the phosphate group show similar behavior for both lipid systems, with the first minimum indicating the hydrogen bond between water and the phosphate oxygen. The RDFs for the systems containing the drug molecules are considerably higher, which can be understood from the density plots, which show that both drug molecules are mostly located in the upper part of the lipid tails and the ester group area. The presence of the drug molecules creates a larger area per lipid in all systems (Table ), which allows more space in the headgroup region that can be filled with additional water molecules, thereby increasing the hydration at the headgroup/extracellular region.

Figure 12

Radial distribution functions of water oxygen around the phosphate of the DOPC and DPPC lipid molecules in the pure system (black), propofol system (red), and fentanyl system (blue).

Table 3

Averaged Structural Properties of Pure and Clinical Concentration Systems, Where Area per Lipid = (AL), Volume per Lipid = (VL), Isothermal Area Compressibility Modulus = (KA), Bilayer Thickness = (DHH), and Lateral Lipid Diffusion Coefficient = (D)

system	AL (Å2)	VL (Å3)	KA (m Nm–1)	DHH (Å)	Dxy (10–8 cm2 s–1)
DOPC pure	70.1 ± 0.2	1274.2 ± 1.1	285.2 ± 19.7	37 ± 0.2	6.3
DPPC pure	62.9 ± 0.3	1175 ± 0.2	231.9 ± 22.7	37.8 ± 0.2	8.8
DOPC propofol	76.8 ± 0.4	1363.8 ± 2.3	666.7 ± 31.0	37.4 ± 0.3	5.8
DOPC fentanyl	71.0 ± 0.3	1283.4 ± 1.6	394.7 ± 15.9	37.25 ± 0.2	6.4
DPPC propofol	70.3 ± 0.4	1280 ± 0.9	312.1 ± 26.0	37.7 ± 0.4	6
DPPC fentanyl	66.0 ± 0.2	1242.2 ± 1.4	256.2 ± 14.0	37.6 ± 0.3	6.4

To assess how these drug molecules at clinical concentration alter the membrane ordering and dynamics, we have calculated headgroup tilt angles, SCD order parameters (Figure ), and lipid lateral diffusions (Table ). Order parameter plots show little or no changes when single drug molecules were incorporated into the bilayers (Figure ), but when clinical concentrations are used, we observe more significant changes in the membrane dynamics. The tilt angle of a phospholipid headgroup is an important property because of the dipole moment associated with the zwitterionic headgroup, which is involved in long-range electrostatics, which can affect many of the bilayer properties.[39] Drug molecules, which are located in or close to the headgroup region, may cause disruption of this angle, which was observed for lidocaine[40] and articaine.[33]Table shows the average calculated angles between the P–N headgroup vector and the bilayer normal for two pure reference systems and the systems containing clinical concentrations of the drug molecules. Results are obtained using the MEMBPLUGIN tool[23] for VMD. For these calculations, all lipids were included as opposed to lipids located close to each drug molecule due to the higher number of drug molecules in the systems. The data presented in Table show small influences from the drug molecules, more so for propofol due to the increased number of molecules in the system compared to fentanyl. Results for single-molecule simulations were identical to those for the pure systems, so they are not shown here. More significant results were obtained for the previously mentioned lidocaine and articaine as these molecules are charged, which caused a decrease in the angle by around 20° due to the charged molecules located within the lipid headgroup area, so the positive charge on these molecules causes a repulsion of the choline groups outwards toward the water phase. Both propofol and fentanyl in our simulations are neutral, and they therefore reside in the upper chain region near the ester group and slightly deeper toward the membrane center for fentanyl (Figure ), so there is no charge–charge interactions to significantly alter the P–N vector. The increased hydration of the headgroups (Figure ) can therefore be explained by the drug molecules causing separation in the headgroup region, which is seen as an increase in the area per lipid, shown in Table , which allows water molecules to penetrate deeper into the headgroup region where hydrogen bonds are formed to the drug molecules.

Table 4

Average DOPC and DPPC Angles between the P–N Vector and the Normal of the Bilayers for Pure and Drug-Containing Systems

system	P–N vector angle (deg)	standard error (±)
pure DOPC	68.88	2.30
propofol DOPC	70.60	2.31
fentanyl DOPC	70.14	2.30
pure DPPC	69.00	2.4
propofol DPPC	71.1	2.0
fentanyl DPPC	70.02	2.4

The calculated order parameters (Figure ) show that these drugs at clinical concentrations have a significant impact on the dynamics of the lipid chains. For DOPC-containing propofol, we calculate an increase in KA, which suggests a rigid/stiff bilayer, an increase in thickness, and a decrease in lipid lateral diffusion. The SCD plot shows a higher-order parameter in the upper chain region for propofol, more apparent in the sn-1 chain, resulting from the propofol occupying this space, which reduces the density of lipid tails here, causing them to straighten. Similar results are seen for fentanyl, in which KA and DHH also increase but to a smaller extent compared to those for propofol. The diffusion calculated for fentanyl in DOPC is almost identical to that of the pure system, which suggests that the number of molecules present in the system is a crucial factor in lipid lateral diffusion. The SCD order parameters for carbons 9 and 10 in both DOPC systems show a significant decrease; these carbons form the double bond in the DOPC lipid chain, which is a region where the drug molecules do not spend any notable time. We suggest that this decrease in order is due to the positioning of the drugs near the ester groups, which disrupts the conventional packing of the lipid chains. For the DPPC systems, we see similar trends in KA, DHH, and D compared to the DOPC systems. The upper region of the carbon chains in which propofol resides has an increased order parameter, which suggests stiffening of that region, as seen in the KA value and lower diffusion coefficients. For fentanyl, we see lower-order parameters for carbon 6 onwards in both the sn-1 and sn-2 chains, suggesting higher-chain mobility in this region, but this is not what we see in the diffusion coefficients, which show a decrease in membrane fluidity.

Potentials of Mean Force

The potential of mean force (PMF) is an important part of the overall membrane permeability, which can be obtained fromwhere R is the resistivity, P represents the permeability, β is the thermodynamic β (β = 1/kBT), W(z) is the PMF, D(z) is the local diffusivity coefficient, and z is the variable, which describes the position of the solute along the transmembrane axis. The PMFs were calculated by pulling the drug molecule from the center of the membrane into the water phase, which has been shown by Filipe et al.[41] to give faster convergence of the PMF compared to pulling from the water phase into the membrane. The PMFs were then symmetrized, i.e., the profiles were adjusted to be identical on either side of the bilayer center and both start and end at 0.[42] Our chosen method of computing the PMFs for our systems was the Umbrella sampling.[43] This approach involves applying a harmonic potential between the center of mass (COM) of the drug molecule and the bilayer, with a harmonic force constant of 2.5 kcal mol–1 rad2. Figure shows the PMF profiles for propofol in DPPC and DOPC. We can see that there is a large difference of 5.54 kcal mol–1 between the free energy at z = 0 in both bilayers. The negative PMF values for both bilayers indicate that z = 0 is a favorable position for propofol, which is not surprising as this molecule possesses a lone ionizable hydroxyl group with a pKa of 11,[44] while the remaining structural components are highly lipophilic. This high lipophilicity gives a poor water miscibility (150 μg L–1) and a log P value of 4.16.[45] As such, we would expect that propofol would be very stable within the hydrophobic, low-density, and low-order acyl chain region. A small dip in the minimum is observed at z = ±1.0 Å in the DOPC system, where the PMF drops to −6.66 kcal mol–1. A similar but more pronounced pattern is seen in the DPPC system, where the PMF drops to a minimum of −1.47 kcal mol–1 at 2 Å further away than in DOPC. This difference of 5.19 kcal mol–1 between the minima of the two systems occurs because propofol is closer to the high-density acyl chain region. We found that the higher order in DPPC is caused by the shorter chain length and full saturation, which allows closer packing of the chains. Propofol prefers to reside close to carbon 13 in DPPC and carbon 17 in DOPC at the PMF minima, which explains their energy difference. As propofol traverses toward the water phase, we observe a steep increase in free energy until it enters the water phase at 18 Å where the free energy is set to 0. For DPPC, we do observe a slight increase in free energy within the headgroup region to a maximum of 0.06 kcal mol–1. This is not unusual due to the highly structured nature of the headgroup region, characterized by strong specific interactions that restrict their motion, which can hinder the diffusion of hydrophobic solutes.[46] To investigate this behavior further, we have plotted in Figure the z-dependent diffusion coefficient profile. From the diffusion profile, we can see the decrease in diffusion within the DPPC system at, and close to, the headgroup region compared to that of DOPC, which accounts for the increase in energy observed compared to that obtained for DOPC. Figure shows a small barrier at 10 Å in the DPPC bilayer, which is most likely due to the excess hydration of the drug molecule, which is overcome as the drug penetrates deeper into the membrane interior; we can see this reflected in Figure , where propofol has diffusion coefficient slightly higher than in DOPC. These results show that the difference in chain length and saturation can significantly alter the drug diffusion, suggesting that propofol will be more selective toward the longer and unsaturated part of the cellular membrane.

Figure 14

Potential of mean force profile for propofol in DOPC and DPPC bilayers (left), and a histogram showing the suitable overlapping of windows (right).

Figure 15

Propofol z-dependent diffusion profile in DOPC and DPPC.

Figure depicts a trend similar to the lower PMF in the DOPC system for fentanyl compared to that for propofol. The minimum energy occurs when fentanyl is at z = 0.0 Å with a PMF of −43.07 kcal mol–1. In the DPPC system, the minimum occurs at z = ±2.0 Å but at 35 kcal mol–1 higher than in DOPC. This result is again expected because fentanyl is a large hydrophobic and lipophilic molecule. We can rationalize the difference in the PMF due to fentanyl’s close proximity to the high-density and less mobile acyl chain region in DPPC. In both lipids, we observe a linear increase in the PMF as fentanyl traverses from the membrane center into the water phase, where the PMF is defined as 0. Unlike propofol, we do not observe the increase in PMF within the headgroup region for fentanyl in the DPPC system, owing to the linear conformation adopted by fentanyl when it diffuses through the headgroups, as observed in our nonbiased simulations. Figure shows the z-dependent diffusion profile for fentanyl in both bilayers. Our umbrella sampling simulations indicate that fentanyl has a significantly lower diffusion and mobility over the range around 20 Å, leading to a higher PMF. However, the free energy is much lower in DOPC, because fentanyl adopts a more favorable conformation within this less ordered environment. As fentanyl moves into the water phase, we see a large increase in diffusion, which almost reaches the same values as it does in the DPPC bilayer. The diffusion profile of the molecules shows the greatest variance between propofol and fentanyl in this study. We should note here that the center of the bilayers are theoretically the most stable region for both drug molecules, as indicated by the PMF plots, but we have to remember that these simulations are carried out with harmonic restraints in each window. In the unrestrained simulations for single drug molecules and clinical concentrations, the drugs do not reside at the bilayer center during the simulations (Figure ) but are mostly located at the upper chain region/headgroup interface due to the highly dynamic nature of the lipid chains.

Figure 17

Fentanyl z-dependent diffusion profile in DOPC and DPPC.

Conclusions

In this study, we have utilized fully atomistic molecular dynamics simulations to model the physical and mechanical properties of pure and drug-containing lipid bilayers. Propofol and fentanyl show similar trends in the way they alter the general properties of both DOPC and DPPC. We determined that propofol prefers to form hydrogen bonds with the carbonyl oxygens of the lipid headgroups, especially with the sn-1 chain in both model membranes. Fentanyl prefers to orient itself parallel to the headgroups at the interface, where it forms hydrogen bonds with water, which are made available as a result of the disruption caused in the headgroup region due to the presence and positioning of fentanyl. Hydrogen-bonding analysis showed an increase in water molecules within the headgroups, most noticeably in the DOPC membrane. The calculated radial distribution functions also showed a higher density of water molecules around the phosphate group when fentanyl was present and to a lesser extent for propofol. Biased MD simulations in combination with umbrella sampling methodology were used to obtain the PMF and z-dependent diffusion profiles for both drug molecules in each lipid bilayer system. Our PMF and z-dependent diffusion profiles highlight the effects of chain length and level of saturation within the bilayer and how this could affect the selectivity for different parts of the cellular membrane by each drug. Our simulations involving the clinical concentrations of propofol and fentanyl show in detail that these drugs can cause significant perturbations to the membrane structure. At these concentrations, both drugs were shown to cause increased hydration of the lipid headgroups, stiffening of the acyl chains, and hence a decrease in the membrane fluidity. The resulting structural defects from our simulations could provide the basis for investigations into the indirect modulation of membrane protein function by both of these drugs, which are the main components of total intravenous anesthesia.

Methodology

Structure Generation

The initial structures of both lipid bilayers were constructed using the CHARMM-GUI Membrane Builder[47] at the experimentally observed hydration levels[48,49] of 32.8 and 30.1 water molecules per lipid for DOPC and DPPC, respectively. Sixty-four lipid molecules were placed in each leaflet to form bilayers consisting of 128 lipids. The PDB files were converted into the lipid14 forcefield format, with the charmmlipid2amber.py script implemented into the Amber simulation package.[50] Structures for propofol and fentanyl were generated using Avogadro[51] and geometry optimized using the Amber-GAFF forcefield[52] with the steepest decent method. Atom types, bond types, and atomic partial charges were assigned using the AM1-BCC semiempirical quantum mechanical method[53] implemented into the antechamber program in Amber.[54] In total, 14 systems were simulated, including two pure bilayer systems to generate equilibrated structures for use with the drug molecules and to reference our computational setup via comparison with experimental results. Four DOPC and DPPC systems were simulated, where each single molecule was placed individually inside either the bilayer hydrophobic phase or the water phase. Four additional simulations were carried out at clinical concentrations, where the drug molecules were added randomly to the water phase. The concentrations used for propofol and fentanyl were 7.1[55] and 1.0 μM,[56] which corresponds to approximately 36 and 4 molecules, respectively.

Simulation Details

All MD simulations in this study were carried out using the Amber16 simulation package,[54] with the lipid14,[26] ff14SB,[57] and GAFF[52] forcefields. The TIP3P water model[58] was used to describe the water molecules. More in-depth details of the general simulation procedure used can be found in the original lipid14 paper.[26] All production simulations were carried out in the isothermal–isobaric (NPT) ensemble, which maintains a constant number of particles, pressure, and temperature. The temperatures were maintained at the lipid-relevant temperature of 303 K for DOPC and 323 K for DPPC to be consistent with the experiment, using the Langevin thermostat[59] with a collision frequency of γ = 1.0 ps–1. The pressure was maintained at 1 atm using the anisotropic Berendsen method,[60] with a pressure relaxation time of 1.0 ps. Three-dimensional periodic boundary conditions with the usual minimum image convention were used. The SHAKE algorithm[61] was used 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 Å. The pure bilayer simulations were run for a total of 125 ns, with the first 25 ns discarded from the final analysis. Simulations containing the drug molecules were run from the equilibrated pure bilayer structure for a total of 100 ns each after a short 5 ns equilibration in which the drug molecules were harmonically restrained to allow equilibration of the solvent or bilayer. The biased MD umbrella sampling simulations were run using the GPU version of Amber16.[62,63] The drug molecule was pulled from the center of the bilayer into the water phase along the z-axis by steered molecular dynamics (SMD), for a total of 27 Å with a pulling rate of 1 Å per ns. The force constant for pulling was set at 1.1 kcal mol–1 rad2, and the force constant for restraint was set at 2.5 kcal mol–1 rad2 for the window simulations. From the SMD simulations, 27 windows were extracted (1 Å apart) along the z-axis from the center of the bilayer to the water phase. Each window was further simulated for 50 ns. The potential of mean force (PMF) profile was constructed using the weighted histogram analysis method (WHAM).[64] The analysis in this study was performed using the CPPTRAJ[35] program implemented into the AmberTools16[54] package, with visualization through tools supplied in VMD.[24] For simulations containing drug molecules, the 5 ns equilibration phase was not included in the analysis. In the single drug molecule simulations, lipids in the leaflet, where the drug molecule interacts, were used for the analysis of lipid properties; in the clinical concentration simulations, all lipids were used.