Interaction of neurotransmitters with a phospholipid bilayer: A molecular dynamics study

Interaction of neurotransmitters with a phospholipid bilayer: A molecular dynamics study

Chemistry and Physics of Lipids 184 (2014) 7–17 Contents lists available at ScienceDirect Chemistry and Physics of Lipids journal homepage: www.else...

2MB Sizes 9 Downloads 75 Views

Chemistry and Physics of Lipids 184 (2014) 7–17

Contents lists available at ScienceDirect

Chemistry and Physics of Lipids journal homepage: www.elsevier.com/locate/chemphyslip

Interaction of neurotransmitters with a phospholipid bilayer: A molecular dynamics study Günther H. Peters a, *, Mikkel Werge a , Maria Northved Elf-Lind a , Jesper J. Madsen a , Gustavo F. Velardez a , Peter Westh b, ** a b

Department of Chemistry, Technical University of Denmark, Kongens Lyngby 2800, Denmark NSM, Research Unit for Functional Biomaterials, Roskilde University, Roskilde 4000, Denmark

A R T I C L E I N F O

A B S T R A C T

Article history: Received 5 July 2013 Received in revised form 14 August 2014 Accepted 22 August 2014 Available online 23 August 2014

We have performed a series of molecular dynamics simulations to study the interactions between the neurotransmitters (NTs) g-aminobutyrate (GABA), glycine (GLY), acetylcholine (ACH) and glutamate (GLU) as well as the amidated/acetylated g-aminobutyrate (GABAneu) and the osmolyte molecule glycerol (GOL) with a dipalmitoylphosphatidylcholine (DPPC) bilayer. In agreement with previously published experimental data, we found the lowest membrane affinity for the charged molecules and a moderate affinity for zwitterionic and polar molecules. The affinity can be ranked as follows: ACH– GLU << GABA < GLY << GABAneu << GOL. The latter three penetrated the bilayer at most with the deepest location being close to the glycerol backbone of the phospholipids. Even at that position, these solutes were noticeably hydrated and carried 30–80% of the bulk water along. The mobility of hydration water at the solute is also affected by the penetration into the bilayer. Two time scales of exchanging water molecules could be determined. In the bulk phase, the hydration layer contains 20% slow exchanging water molecules which increases 2–3 times as the solutes entered the bilayer. Our results indicate that there is no intermediate exchange of slow moving water molecules from the solutes to the lipid atoms and vice versa. Instead, the exchange relies on the reservoir of unbounded (“free”) water molecules in the interfacial bilayer region. Results from the equilibrium simulations are in good agreement with the results from umbrella sampling simulations, which were conducted for the four naturally occurring NTs. Free energy profiles for ACH and GLU show a minimum of 2–3 kJ/mol close to the bilayer interface, while for GABA and GLY, a minimum of respectively 2 kJ/mol and 5 kJ/mol is observed when these NTs are located in the vicinity of the lipid glycerol backbone. The most important interaction of NTs with the bilayer is the charged amino group of NTs with the lipid phosphate group. ã 2014 Elsevier Ireland Ltd. All rights reserved.

Keywords: Simulations Lipid membrane Neurotransmitters Umbrella sampling Water dynamics

1. Introduction Many small molecules have intracellular targets, and hence they must pass across one or more phospholipid bilayer membranes to reach the intracellular targets and elicit a response to their pharmacological action. Therefore, small molecule–lipid interactions are inevitable, and a wide range of small molecules have been demonstrated to interact with lipid membranes (Berquand et al., 2005; Barcelo et al., 2004; Hidalgo et al., 2004; Preetha et al., 2007). The nature of these interactions has been studied using different biophysical techniques (Klacsová et al., 2011; Wanderlingh et al., 2010; Peetla et al., 2009) (and reference

* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (G.H. Peters), [email protected] (P. Westh). http://dx.doi.org/10.1016/j.chemphyslip.2014.08.003 0009-3084/ ã 2014 Elsevier Ireland Ltd. All rights reserved.

therein). However, an understanding on a molecular level of how binding to the lipid evokes the biological response remains limited since it is difficult to probe these interactions on a single-molecule level (Peters, 2004; Karplus, 2012; Srinivas and Klein, 2004). Here, computational methods such as atomic-level molecular dynamics simulations have been used to elucidate the structure and dynamics of lipid bilayer membranes and to probe partitioning and permeation of small molecules (Stouch, 1997; Xiang and Anderson, 2002; Boggara and Krishnamoorti, 2010; Tieleman, 2006; Bemporad et al., 2004). A wide range of solutes ranging from small polar molecules such as amino acids, anesthetics to organic solvents and drugs have been studied for their interactions with different types of lipid bilayer membranes (Klacsová et al., 2011; Rodgers et al., 2010; Mojumdar and Lyubartsev, 2010; Marrink and Berendsen, 1994; Shinoda et al., 2004; Tu et al., 1998; Henin et al., 2010; Högberg and Lyubartsev, 2008; Pitman et al., 2004; Pandit et al., 2008; Bennett et al., 2009; Pedersen et al., 2007; Peters et al.,

8

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

2009; Dickey and Faller, 2007; Cerezo et al., 2011; MacCallum et al., 2007; Johansson and Lindahl, 2008; MacCallum et al., 2008; Mukhopadhyay et al., 2004; Norman and Nymeyer, 2006; Bemporad et al., 2005; Li et al., 2008; Ulander and Haymet, 2003; Boggara and Krishnamoorti, 2010). In particular, the mechanism by which anesthetics work has drawn much attention, and two apparently incompatible theories have emerged. One suggests that anesthetic action arises from direct anesthetic– protein interactions (Franks and Lieb, 1984; Franks and Lieb, 1985; Slatter et al., 1993; LaBella et al., 1998; Franks and Lieb, 1982) including ligand-gated ion channels (Weng et al., 2010; Nury et al., 2011), whereas the other one suggests that the lipids of the neuronal membranes are the prime site of anesthetic action. The indirect, lipid-mediated mechanism has been discussed since the anesthetic potency of a chemical species correlates with its octanol–water partition coefficient known as the Meyer–Overton rule (Meyer, 1899; Overton, 1901). This led to the proposition that anesthetics affect the postsynaptic membrane by modulating its physical properties through interaction with its lipid component (Cantor, 1997; Milutinovic et al., 2007). The lipid bilayer perturbation as the primary event is then transmitted to a membrane protein (Richards et al., 1978). The lipid-mediated mechanism has recently also being linked to the action of neurotransmitters (NTs) (Cantor, 2003; Sonner and Cantor, 2013). Sonner and co-workers provided evidence that ligand-gated ion channels function can be modulated by coreleased NTs in a similar fashion as observed for anesthetics (Milutinovic et al., 2007). The authors showed that non-native NTs can affect receptor function by modulating the mechanical properties of the bilayer. These mechanical changes affect consequently the conformational equilibrium of ligand-gated ion channel receptors and thereby their response to the native agonist. Lipid–NT interactions have also been hypothesized to influence neural transmission if the membrane accumulates a reservoir for neurotransmitters and thereby facilitates binding of these molecules to the target protein (Scheidt and Huster, 2008; Hemmings et al., 2005; Vautrin et al., 2000; Vautrin and Barker, 2003; Jodko-Piorecka and Litwinienko, 2013). The idea of a lipid-mediated mechanism is further supported by a recent work that has suggested that different types of neurotransmitters have affinity for lipid bilayer membranes. Aromatic transmitters such as serotonin and dopamine appears to have high

affinity (Peters et al., 2013; Jodko-Piorecka and Litwinienko, 2013; Orlowski et al., 2012), while smaller more hydrophilic NTs such as glycine, glutamate and GABA may show moderate affinity which depends strongly on the composition of the lipid membrane (Wang et al., 2011). These observed affinities appear to be a necessary requirement for an indirect role of bilayer–NT interactions in neural transmission. Analysis of this putative effect clearly requires a deeper insight into e.g., driving forces and structural characteristics of NT–lipid interactions. The affinity of the aromatic NTs was suggested to rely on contact between the (cationic) primary amine of 5-HT and the lipid phosphate group. This provided a strong affinity of 5-HT for the membrane interfaces in spite of the fact that serotonin is hydrophilic with an oil–water partitioning coefficient well below unity. In the current work, we investigate the origin of the weaker membrane affinity found for the non-aromatic NTs with a phosphatidylcholine bilayer using equilibrium molecular dynamics (MD) simulations and umbrella sampling simulations. Specifically, we have calculated the energy cost (potential of mean force (PMF)) for partitioning of NTs into the bilayer, identified the most favorable NT–lipid contacts and analyzed positional distribution and hydration of interfacially located NTs. The NTs studied were the zwitterionic neurotransmitters g-aminobutyrate and glycine, as well as the charged NTs acetylcholine and glutamate. For comparison and to address the effect of charges on the absorption properties to the bilayer, we also included amidated/ acetylated g-aminobutyrate and the osmolyte molecule glycerol in our study. 2. Methods 2.1. Equilibrium simulations MD simulations were performed for systems composed of dipalmitoylphosphatidylcholine (DPPC)/water/solutes. The solutes were: glycerol (GOL); the zwitterionic neurotransmitters: g-aminobutyrate (GABA) and glycine (GLY); the charged neurotransmitters: acetylcholine (ACH) and glutamate (GLU); and the amidated/acetylated g-aminobutyrate (GABAneu). Their structures are displayed in Fig. 1. The abbreviations given in parentheses are introduced for convenience and will be used throughout the text when referring to a certain solute molecule. The bilayer consisted of 72 DPPC molecules (36 per leaflet) and was fully hydrated with

Fig. 1. Structure of the different solutes: g-aminobutyrate (GABA), glycine (GLY), acetylcholine (ACH), glutamate (GLU) and acetylated/amidated g-aminobutyrate (GABAneu) and glycerol (GOL) as well as 1,2-dihexadecanoyl-sn-glycero-3-phosphocholine (DPPC). For clarity, the hydrogen atoms in the DPPC molecule are not shown. The phospholipid atoms N (choline group, lipidN), P (phosphate group, lipidP), C1 (carbonyl carbon of glycerol backbone, lipidC1) and C16 (carbon in the methyl group (tail), lipidC16) are chosen for the calculation of the probability distributions.

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

3400 water molecules (47 water molecules per lipid molecule). The initial DPPC/water configuration was taken from a previous simulation study (Sonne et al., 2007). The center of mass (COM) of the bilayer was located at the origin of the coordinate system (x = y = z = 0 Å) with z being the normal direction to the bilayer plane. For each system, three to five independent simulations were performed to assess the statistical uncertainties in the calculated quantities. The different starting configurations were obtained by placing one NT molecule at different distances (z) away from the interfacial bilayer region. Placing the NT molecule was carried out using the Visual Molecular Dynamics (VMD) 1.9.1 software package (Humphrey et al., 1996). Water molecules within 2.4 Å of the NT molecule were deleted using VMD. To neutralize the charged systems DPPC/ACH and DPPC/GLU, one water molecule in each system was randomly picked and replaced with chloride and sodium ion, respectively. Simulations were performed using the NAMD 2.5 software package (Phillips et al., 2005) with a modified CHARMM27 force field for DPPC (Sonne et al., 2007) and the TIP3P water model (Jorgensen et al., 1983). We previously re-parameterized the force field for the DPPC head group region based on the experimental data from Nagle et al., (1996). The velocity Verlet algorithm with a time step of 0.001 pico second (ps) was used to solve Newton’s equations of motion. Simulations were carried out in the NPT ensemble at constant pressure (P = 1 atm) and temperature (T = 323 K; which is above the experimental gel-to-fluid transition temperature of 314 K). Temperature and pressure were controlled respectively by the Langevin thermostat (damping coefficient: 5 ps1) and Nosé–Hoover Langevin barostat. Anisotropic pressure regulation was applied in all simulations, and piston damping coefficient, piston period and piston decay for the Langevin barostat (Feller et al., 1995) were set respectively to 5 ps1, 0.200 ps and 0.500 ps, as previously applied (Peters et al., 2009; Sonne et al., 2007). The long-ranged electrostatic forces were calculated using the Particle Mesh Ewald method (Essmann et al., 1995). The grid spacing was approximately 1 Å, and a fourth order spline was used for the interpolation. Electrostatic forces were updated every fourth femto second. Van der Waals interactions were cut off at 12 Å in combination with a switching function starting at 10 Å. Periodic boundary conditions were applied in x,y and z-directions. The potential energy in all systems was initially minimized using 2000 steps of the conjugated gradient method. Depending on the neurotransmitter, the length of the individual simulations varied between 24–106 nanoseconds (ns) amounting to a total simulation length of 97–307 ns (Table 1). The time evolution of the area/ phospholipid molecule was used to determine the length of the equilibration phase. Depending on the neurotransmitter, this resulted in omitting the initial 5–20 ns of a particular trajectory for analysis. Analyses were performed using VMD (Humphrey et al., 1996).

9

2.2. Water dynamics From the equilibrium MD simulations, we assessed the dynamics of the water molecules by computing their mean residence times in the first hydration shell of solutes and selected groups in the phospholipids. The definition of residence time used follows that of Makarov et al., (2000) and Schröder et al. (2006) in which one considers the binary function xi,a(t) where i denotes a particular solvent molecule, a refers to solutes or lipid atoms and t denotes the time. The function is one when the solvent molecule is within a 3.5 Å cutoff distance from the selected entity and zero otherwise. The calculation of the autocorrelation function ra(t) for the sequence xi,a(t) was based on sliding windows of 800 ps. The overall autocorrelation function r(t) for adsorbed solvent molecules was calculated as the average ra(t) over the particular solute or selected lipid group, and a bi-exponential model (Makarov et al., 2000) was fitted to the function r(t). The model, which reproduced the sampled r(t) well, is defined by

rbiexp ðtÞ ¼ a1 expðt=t 1 Þ þ a2 expðt=t 2 Þ

(1)

where a1, a2, t1 and t 2 are adjustable parameters. t1 and t 2 are the residence times of slowly and rapidly exchanged water molecules, respectively. The relative populations of slowly and rapidly exchanged water molecules, Pslow and Pfast, can be estimated from the relative magnitudes of a1 and a2 (Makarov et al., 2000). 2.3. Umbrella sampling simulations Umbrella sampling simulations were performed with the NAMD 2.7 software package (Phillips et al., 2005) using tcl script for adding biasing forces (NAMD’s TclForces functionality) at the constraining z-position during the simulations. To conduct umbrella sampling, a series of configurations was generated by placing NTs at different z-position (i.e., along the bilayer normal); 44 Å < z < 44 Å; Dz = 1 Å. When possible, initial configurations were taken from the equilibrium simulations. For z positions of solute molecule in the bilayer, which were not reached during equilibrium simulations, a pulling force was applied to the center of mass of the solute molecule to translate a solute molecule closest to the desired (umbrella sampling) z position to that z position. We have used the protocol previously successfully applied by MacCallum et al., (2008) by placing two neurotransmitters in the simulation cell with maximum distance in zdirection from each other and using semi-isotropic pressure regulation. The force constant for the harmonic umbrella potential was 8.374 kJ/mol. Force field parameters and simulation conditions were the same as used in the equilibrium molecular dynamics simulations described above. The system consisted of two NTs, a bilayer of 72 phospholipid molecules and 5200 water molecules. Umbrella sampling simulations were carried out for 40 ns at each z-position, and the potential of mean force profiles along the bilayer normal (PMF(z)) were constructed using the weighted

Table 1 Details about the different simulations: number of simulations for each system and length of simulations. ACH: Acetylcholine, GLU: glutamate, GLY: glycine, GABA: g-aminobutyrate, GABAneu: acetylated/amidated g-aminobutyrate, GOL: glycerol, DPPC: 1,2-dihexadecanoyl-sn-glycero-3-phosphocholine.

DPPC/ACH DPPC/GLU DPPC/GLY DPPC/GABA DPPC/GABAneu DPPC/GOL DPPC a

Number of simulations

Simulation time (ns)

Total simulation time (ns)

3 3 5 5 3 5 1

48, 25, 24 106, 102, 99 69, 36, 33, 30, 27 50, 26, 24, 24 56, 39, 34 72, 37, 35, 30, 28 108

97 307 195 124 129 202 108a

Started from a previous study (Sonne et al., 2007).

10

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

histogram analysis method (WHAM) (Grossfield, 2012). In all simulations, the umbrella histograms of configurations overlap well with the adjacent windows (data not shown). Convergence of the simulations was tested by calculating the profiles in sliding windows of 0.5 ns. Errorbars were estimated using the blocking method as described by Flyvbjerg and Petersen (1989). The time series were divided in different block sizes, PMFs were calculated for each block, and the variance were estimated. For all systems, the variance reached a constant value at 25 block transformations. 3. Results and discussion 3.1. Area per lipid molecule The simulation of a pure DPPC membrane was continued from a previous simulation (Sonne et al., 2007), and the simulation was carried for an additional 110 ns. The time evolution of the area per phospholipid molecule for the pure DPPC bilayer is shown in Fig. S1 (Supplementary material). The average area per lipid molecule is 68  2 Å2, which is within 5% of the experimental value of 64 Å2/ lipid molecule (Nagle and Tristram-Nagle, 2000; Kucerka et al., 2006). Addition of solutes resulted generally in a slight increase in the fluctuation of the area per phospholipid molecule. Averaged over all DPPC/small organic molecule simulations, the area was found to be 69  3 Å2/lipid molecule, which is within the statistical uncertainties of the area determined for the pure DPPC bilayer. 3.2. Distribution along bilayer normal We monitored the time evolution of selected phospholipid atoms, water molecules and solutes. We selected four atoms in the lipid molecule which are shown in the image for the DPPC molecule in Fig. 1. The lipid atoms are (i) the nitrogen atom in the choline group (lipidN); (ii) the phosphorous atom in the phosphate group (lipidP); (iii) the carbonyl carbon atoms in the ester groups (lipidC1); and (iv) the carbon atoms in the methyl groups located at the end of the aliphatic chains (lipidC16). These selected lipid atoms encompass the molecule’s extent in the z-direction. For the water molecules and the solutes, we selected the oxygen atom and the center of mass (COM), respectively. From the time evolution of the coordinates, we determined the positional probability distribution of water molecules, solutes and selected lipid–molecule atoms in the z-direction, perpendicular to the membrane surface. Selected distributions are displayed in Figs. 2 and 3. The water profile shows an almost constant number density outside the membrane in the z-range 25 to 34 Å. The mass density of water estimated from the GABAneu simulations is 0.997  0.025 g/cm3, which is in close agreement with the experimental value of 0.98807 g/cm3 at 323 K (Weast et al., 1985). The penetration of the water molecules into the membrane was quite substantial and extended down to the region of the carbonyl group in the ester groups of the lipids. From 10 Å down to 5 Å, the profile of the water molecules follows the profile of the carbonyl groups very closely, which may indicate that the water molecules in that region are tightly bound to the lipid carbonyl groups. This is supported by our analysis of the dynamics of water molecules that will be discussed in more detail below. A general feature of the positional probability distribution functions for all the selected lipid atoms is that they are quite broad with a fullwidth-at-half-maximum 7–8 Å. This feature implies that the membrane surface is quite rough, a result which is consistent with other simulations including different membrane systems (Nagle and Tristram-Nagle, 2000; Edholm and Nagle, 2005; Bhide and Berkowitz, 2005; Broemstrup and Reuter, 2010; Peters et al., 2009; Pedersen et al., 2010; Debnath et al., 2010; Murzyn et al., 2001). The position of the maximum in the distribution functions for the

Fig. 2. Positional probability distributions for selected groups along the bilayer normal, z, extracted from the different trajectories after the area/phospholipid molecule reached a constant value. Profiles are shown for the systems: (a) DPPC/ GABA and (b) GABAneu. Profiles are shown for solutes (center of mass, COM), lipidN, lipidP, lipidC1, lipidC16 and water (oxygen atom, w0). See Fig. 1 for details. Profiles were smoothened with a running average of five. An interval of Dz = 0.5 Å was applied for calculating the profiles. Profiles are shown for the systems: (a) DPPC/ GABA and (b) DPPC/GABAneu.

nitrogen and phosphorous atoms are only a few Angstroms apart. This implies that the dipole vector between them is almost parallel to the membrane surface, which is in good agreement with the reported NMR results (Seelig et al., 1977). To quantify membrane affinity, we defined the membrane:bulk ratio, R, of solute R Zin PðzÞdz R ¼ R Z0 (2) out PðzÞdz 0 where P(z) is the positional distribution function. Solute molecules were considered to be in the membrane when zCOM, solute < 22 Å (i.e., zin = 22 Å in the integration limits in Eq. (2) and otherwise in the water phase (zout = 40 Å). This resulted in the following R values: 0.05  0.02% (ACH), 0.04  0.02% (GLU), 0.10  0.03% (GABA), 0.15  0.05% (GLY), 0.27  0.09% (GABAneu) and 0.63  0.14% (GOL). The standard deviations of the mean were based on the replica simulations. This picture of a low or moderate affinity of NTs for a pure PC membrane is in line with a recent experimental study (Wang et al., 2011). The charged neurotransmitters, ACH and GLU, had very low R values (i.e., very low affinity for zwitterionic PC bilayers). The affinity for the zwitterionic neurotransmitters, GABA and GLY, increased where

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

Fig. 3. Positional probability distributions for selected groups along the bilayer normal, z, extracted from the different trajectories after the area/phospholipid molecule reached a constant value. Profiles are shown for the systems: (a)0 DPPC/ ACH and (b) DPPC/GOL. See Fig. 2 caption for details.

GLY showed a higher affinity towards the membrane. This balance between the affinity of GABA and GLY, which has also been found experimentally (Wang et al., 2011), is interesting because GABA has more hydrophobic surface area than GLY (Fig. 1). This suggests that hydrophobicity is less important for membrane partitioning of NTs than polar interactions and molecular size. GLY is smaller than GABA and may therefore more easily enter the interfacial bilayer region. Reducing the charge on the neurotransmitter had a significant effect on the absorption of NT. Neutralization of the charges on GABA by amidation and acetylation (GABAneu) resulted in a significant increase in the affinity and led to a deeper penetration into the interfacial bilayer region (Fig. 2). Glycerol had the highest affinity between the solutes considered here. Membrane–glycerol systems have been extensively studied because glycerol is a common osmolyte; i.e., a solute that accumulates in organisms to regulate osmotic pressures and mitigate deleterious effects on membranes and other structures during exposure to drought and cold (Yancey et al., 1982). Perhaps the most fundamental observation is that the PC bilayers were moderately permeable to GOL (Mitragotri et al., 1999), and hence the solute must diffuse into the bilayer albeit to a limited extent. This was confirmed in equilibrium studies showing that the bilayer–water partitioning coefficient for GOL was detectable but far below unity (Katz and Diamond, 1974). Other thermodynamic measurements have suggested a weak preferential hydration of

11

bilayer vesicles (Westh, 2003); that is a slightly reduced glycerol concentration at the membrane interface compared to the bulk composition. In contrast to this, X-ray reflectivity has suggested accumulation of glycerol at the interface of lipid monolayers (Pocivavsek et al., 2011). One interesting aspect of this discussion is the significant co-partitioning of water found for GOL, GLY and GABAneu (see Table 2 and the discussion below). These solutes drag along 5–12 water molecules even when located quite deeply in the membrane (zCOM, solute < 15 Å). This implies that the (molal) concentration of solute in the bulk is in fact decreased as a result of partitioning (because more water than solute is removed from the bulk phase upon partitioning). This behavior could be the origin of the observed discrepancies between thermodynamic and “structural” methods (such as partitioning measured by radioactive tracers). Thus, based on the current results an equilibrium approach would report preferential hydration of the membrane (Schellman, 1993; Timasheff, 1998) because the bulk concentration and hence, the chemical potential of GOL, GLY and GABAneu increases upon addition of lipid (more water than solutes is removed from the bulk). However, a structural approach would detect an increased presence of these solutes in the interfacial region and hence report accumulation. Such co-partitioning of water and small polar solutes has been proposed earlier (Patra et al., 2006) and may be an integral part of these molecules’ effect on lipid membranes. Further penetration of these NTs was not observed and was probably restricted due to the hydrophobicity of the lipid chain region and the flexibility of the tail of the aliphatic chains. There was a small probability that the lipid tails (methyl groups, C16; Figs. 2 and 3) reached the glycerol backbone region (C1; Figs. 2 and 3). This entropic contribution may restrict how far these molecules can penetrate into the bilayer. This is consistent with our umbrella sampling simulations (discussed below), where we found a relatively high energy barrier for solutes to cross the middle region (at z  0) of the bilayer. The entering of solutes had only small effects on the conformation of the molecules and only a 15% decrease in the extended conformation was observed when they entered the bilayer (data provided in Supplementary material). This is expected due to their relatively small sizes. 3.3. Solute–phospholipid interactions To identify key interactions between the solutes and the bilayer, we monitored the distances between selected atoms of the solutes and phospholipids. For the lipid, we selected the phosphorous atom (lipidP) in the phosphate groups and the nitrogen atom (lipidN) in the choline groups. For the solutes, we selected the amino group, carboxyl group or polar groups depending on the structure of the molecule. The analysis was only performed for molecules that have entered the bilayer; i.e., zCOM, solute < 22 Å. The

Table 2 Number of water molecules within 3.5 Å of acetylcholine (ACH), glutamate (GLU), glycine (GLY), g-aminobutyrate (GABA), acetylated/amidated g-aminobutyrate (GABAneu), and glycerol (GOL) determined at different distance from the center of the bilayer, which is located at x = y = z = 0. The surface normal of the bilayer lays parallel to the z-axis. For the analysis, the center of mass of the solutes was chosen to define their z-position, and three regions were chosen: 0 < zCOM, solute < 15 Å, 15 < zCOM, solute < 22 Å and zCOM, solute > 22 Å. Subscript on the number of water molecules refers to the standard deviation. Molecule

0 < zCOM < 15 Å

15 < zCOM < 22 Å

zCOM > 22 Å

ACH GABA GABAneu GLU GLY GOL

– – 104 – 122 52

243 162 136 183 122 63

263 192 234 212 152 163

12

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

probability distributions for the DPPC/GLY system are shown in Fig. 4 for distances between lipidP/lipidN and the amino group (GLYN; Fig. 1)/carboxyl group (GLYC; Fig. 1). There was a high probability for finding strong interactions between the lipidP and the positively charged amino group in the neurotransmitter GLYN at a minimum distance of 3.5 Å. Other interactions are less pronounced than lipidP–GLYN. Interactions between the positively charged lipid choline group and the negatively charged carboxyl group in GLY were seen, but the distribution is much broader and the maximum in the probability distribution is 4 times lower than observed for the lipidP–GLYN distribution. Broad distributions are also observed for lipidP–GLYC and lipidN–GLYN. This is expected since the analysis was restricted to GLY located in the interfacial bilayer region. Probabilities and standard deviations calculated for the solutes are listed in Table S1 (Supplementary material). Besides for ACH and GABAneu, similar results were observed for the other neurotransmitters; i.e., relatively high probabilities are seen for lipidP–NTN interactions at 3.5 Å (Supplementary material). For ACH that showed low affinity for the PC bilayer, and GABAneu that penetrated the head group region of the bilayer, the probabilities of the different interactions are similar with distances between 5– 8 Å. In general, the standard deviations are relatively high, which reflects the dynamical behavior of the solutes in the interfacial bilayer plane. The general picture that emerges is that once the NTs enter the interfacial bilayer region, the molecules interact strongest with the phosphate groups. However, as discussed before, NTs showed significantly different affinity towards the PC bilayer. Therefore, lipidP–NTN interactions are important, but not strong enough to anchor the NTs at the site of the phosphate groups. The most pronounced interactions between GOL and the lipid molecules were seen for lipidP with the hydroxyl groups in GOL. Within the uncertainties, there was no difference between the hydroxyl groups in GOL. Again, this is an indication for the dynamical behavior of the solute in the interfacial bilayer region. Throughout the course of the equilibrium MD simulations, we have not observed any 360 degrees rotation of NTs located in the interfacial bilayer plane. Probability distributions of the tilt angle of NTs calculated from the angle between the long-axis of NT and

Fig. 4. Probability distribution of distances between selected GLY atoms (C and N; Fig. 1) and selected phospholipid atoms (phosphorus atom (P) in the phosphate groups and nitrogen atom (N) in the choline groups; Fig. 1). Data were extracted from the different replica simulations after the area/phospholipid molecule reached a constant value. Error bars were estimated from the replica simulations (Table 1). An interval of Dd = 0.5 Å was applied for calculating the distribution function.

the bilayer normal displayed maxima around 30–60 . That means that NTs located in the bilayer are tilted with the amino group pointing towards the bilayer interface. These configurations were used as the starting configurations in the umbrella sampling simulations as discussed below. 3.4. Hydration and water dynamics We recorded the hydration of solutes as a function of the COM position of the molecules along the bilayer normal (i.e., along the zaxis). Any water molecule within 3.5 Å of the solutes was considered to be in the first hydration shell. The cutoff corresponds to the first minimum in the oxygen-oxygen radial distribution function for bulk water (Jorgensen et al., 1983). The analysis was done for three regions: 0 < zCOM, solute < 15 Å, 15 Å < zCOM, solute < 22 Å and zCOM, solute > 22 Å (Table 2). A general trend can be seen. As shown in Fig. 5, the hydration of the molecules was reduced as the solutes entered the bilayer. Besides for GABAneu and GOL, the number of water molecules in the hydration shell reduced by 10– 20% as the molecules entered the interfacial bilayer plane (15 Å < zCOM, solute < 22 Å). For GABAneu and GOL, the decrease in the number of water molecules in the first hydration shell amounted to 40–50%. Only GABAneu, GLY and GOL reached a z-position in the bilayer below 15 Å. Even at that location in the bilayer, these molecules were noticeable hydrated (Table 2 and Fig. 5) and carried 40% (GABAneu), 80% (GLY) and 30% (GOL) of the bulk water along. The largest contribution to the hydration stemmed from the polar groups (data not shown). To gain further insight into the behavior of water molecules, we studied the water dynamics. Water dynamics at hydrophobic and hydrophilic surfaces have been a subject of experimental and computer simulation studies for decades, and these studies suggested that translational motion of water molecules near the surface is restricted (Gruenbaum and Skinner, 2011; König et al., 1994; Jensen et al., 2003; Jensen et al., 2004; Foglia et al., 2010; Murzyn et al., 2006). Similarly, simulation studies on water

Fig. 5. Number of water molecules (No_H2O) within 3.5 Å of acetylcholine (ACH), g-aminobutyrate (GABA), acetylated/amidated g-aminobutyrate (GABAneu), glutamate (GLU), glycine (GLY) and glycerol (GOL). Three intervals were chosen along the surface normal of the membrane (parallel to the z-axis) to define the position of the molecules. The position was based on the center of mass of the molecule, and the intervals are as follows: 0 < zCOM < 15 Å, 15 Å < zCOM < 22 Å and zCOM > 22 Å. The center of the bilayer is at x = y = z = 0 Å. The number above the bars at 15 < zCOM < 22 Å is the ratio No_H2O (zCOM > 22 Å)/No_H2O (zCOM < 22 Å) and reflects the reduction in the number of water molecules within 3.5 Å of solute molecules as the solute molecules enter the interfacial bilayer region.

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

13

Table 3 Residence times (t 1 and t 2) and relative populations (Pslow and Prapid) of slowly and rapidly exchanged water molecules calculated for the different neurotransmitter, glycerol and selected groups in the head group region of the phospholipids. PC–CO: lipid ester group, PC–N: lipid choline group, PC–P: lipid phosphate group, ACH: acetylcholine, GLU: glutamate, GLY: glycine, GABA: g-aminobutyrate, GABAneu: acetylated/amidated g-aminobutyrate, GOL: glycerol. For comparison, residence time and relative population data are included for a pure water simulation. Subscripts indicate standard deviation, estimated from the replica simulations.

PC–CO PC–N PC–P ACH GLU GLY GABA GABAneu GOL Pure water

t 1 bulk (ps)

t2 bulk (ps)

Pslow bulk

Pfast bulk

t1 membr. (ps)

t2 membr. (ps)

Pslow membr.

Pfast membr.

%-in membr.

– – – 11517 12518 11532 10427 12221 14175 600200

– – – 122 161 141 144 122 135 405

– – – 244 203 206 264 252 175 141

– – – 764 803 806 747 754 836 861

71631 67040 52232 502437 18648 19048 21634 35434 32067 –

394 212 381 1814 104 1713 2215 361 2710 –

721 481 483 416 478 425 495 495 6013 –

281 521 523 596 538 585 515 515 4013 –

– – – 33 41 189 67 2417 6729 –

dynamics near POPC bilayer surface showed that orientational order of water molecules extends up to 6 Å from the interfacial membrane surface into the aqueous phase (Bhattacharyya and Bagchi, 2000; Bellissent-Funel, 1998; Marrink et al., 1996; Rog et al., 2002). Inside this region, anisotropic lateral diffusion and relatively short rotational correlation time of neighboring water were observed (Murzyn et al., 2001; Åman et al., 2003). Similar results were reported for negatively charged phosphatidylserine and ganglioside bilayers which clearly demonstrate the importance of heterogeneous surface on the dynamics of water molecules (“boundary effects”) (Bhide and Berkowitz, 2005; Sega et al., 2005). In the present study, the dynamics of water molecules was determined by computing the residence time of water molecules in the first hydration shell of solutes and selected phospholipid groups: choline group (selected atoms: N, C), phosphate group (P, O) and carbonyl group of the ester groups in the two aliphatic chains (C, O). The two residence times observed for water molecules in the interfacial membrane plane are different by an order of magnitude. The slow and fast dynamics occurred respectively in the range 500–700 ps and 20–40 ps (Table 3). Similar residence times have been observed in POPC bilayers (Rog et al., 2002; Bhide and Berkowitz, 2005; Gruenbaum and Skinner, 2011). These values should be compared to the results obtained for pure water simulations, where t 1 = 0.6  0.2 ns, t 2 = 40  5 ps and Pslow = 14  1% (Table 3) (Wedberg et al., 2012). Although no trend could be observed in the residence times determined for water molecules residing at phospholipid atoms, the percentage of slow vs. fast moving water molecules varied. The percentage of slow exchanging water molecules increased as water molecules entered deeper into the membrane, which indicates relatively strong binding of water molecules to lipid atoms, and hence an increase in ordering of water molecules in this region, which also been reported for other bilayer systems (Debnath et al., 2010; Murzyn et al., 2001). The effect was strongest for water molecules that were hydrogen bonded to the carbonyl oxygen atoms (PC– CO), where the ratio between slow and fast moving water molecules, Pslow/Pfast, is 2.6. For water molecules located around the phosphate (PC–P) and choline groups (PC–N), Pslow/Pfast = 0.9. In pure water, Pslow/Pfast = 0.2. Hence, the percentage of slow moving water molecules follows: carbonyl oxygen atoms > phosphate and choline groups >> pure water. The observed decrease in the mobility of water molecules in the interfacial compares well with recent experimental data reported by Osborne and co-workers. The authors applied multidimensional infrared spectroscopy and could show that the dynamics of water molecules in the interfacial bilayer is nearly three-fold decreased relative to bulk water molecules (Osborne et al., 2013). Although it is not directly comparable, it is intriguing that we observe for the

pure bilayer that the percentage of slow exchange water molecules increases 3 times; Pslow,bulk = 14% to Pslow,bilayerPC–P = 48% (Table 3). For the solutes, we performed the analysis for zCOM, solute > 22 Å (bulk phase) and zCOM, solute < 22 Å (bilayer region). For z > 22 Å, we found t 1 120 ps, t 2  10 ps and Pslow  20%. As the solute molecules enter the bilayer, Pslow increases by 40–60% (Table 3). In both regions, the residence times for fast exchanging water molecules were similar. However, t1 referring to the slow dynamics of water molecules in the first hydration shell of solutes increases 2–3 times in the interfacial bilayer region. When compared to the lipid atoms, t1 for PC–CO, PC–P and PC–N was 2–3 times lower than t 1 observed for the solutes. Since the time scales are different, it appears that there was no intermediate exchange of slow moving water molecules from the solutes to the lipid atoms and vice versa. The slow exchange of water molecules only happened between unbounded (“free”) water molecules in the interfacial bilayer region and water molecules in the hydration shell of solutes. Hence, even though the hydrated solutes enter the interfacial bilayer plane, water molecules are still mobile, which aids dehydration of the solutes. 3.5. Umbrella sampling simulations Umbrella sampling simulations were conducted for the four naturally occurring NTs. A series of configurations was generated by placing NTs at different z-position (i.e., along the bilayer normal); 44 Å < z < 44 Å; Dz = 1 Å. Each window was simulated for 40 ns. The dynamics (tumbling) of NTs around the z-position was monitored by calculating the probability distribution of the

Fig. 6. Probability function (P(u)) of the orientation of GLY at different z position given in the inset.

14

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

Fig. 7. PMF profiles for (a) DPPC/ACHO, (b) DPPC/GABA, (c) DPPC/GLU and (d) DPPC/GLY. Errorbars show the standard deviations as estimated by blocking method described by Flyvbjerg and Petersen (1989).

orientation of the molecules (P(u)) with respect to the bilayer normal (z-axis). An example is of P(u) is shown in Fig. 6 for GLY. Probability distribution for the other NTs is provided in Fig. S2 (Supplementary material). Free tumbling is indicated by a horizontal line. As indicated Fig. 6, free tumbling is observed in the bulk water phase and close to the bilayer interface. As NTs enter the membrane, sampling at u < 300 and u > 1500 is lower than in the region 300 u < 1500. This is expected, since it is energetically favorable that NTs align with the alkyl chain of the bilayer lipids. The umbrella sampling simulation results are presented in Fig. 7, and minima in the PMF profiles are seen at 29 Å (3.4 kJ/mol) for ACHO, 26 Å (2.2 kJ/mol) for GLU, 22 Å (1.6 kJ/mol) for GABA and 21 Å (5.0 kJ/mol) for GLY. ACHO and to a certain extent GLU are favorably located at the choline group, whereas GABA and GLY are favorably located at the glycerol backbone/phosphate moiety. Within 3–4 Å, the position of the minima observed in the PMF profiles agrees with the maxima observed in the positional probability distribution of the equilibrium dynamics simulations (Figs. S3, Supplementary material). We also calculated the free energy from the equilibrium simulations using the process of binding unbinding events. A binding event occurred when the center of mass of NTs where below z < 30 Å. The calculated energies are 0.3  0.9 kJ/mol for ACHO, 0.2  1 kJ/mol for GLU, 0.2  0.6 kJ/mol for GABA and 1.3  0.9 kJ/mol for GLY. The statistical uncertainties are significant and the only pattern that

emerges is that GLY has a significant lower free energy than the remaining NTs. 4. Conclusion We have used MD simulations to study interactions of neurotransmitters (GLY, GABA, ACH and GLU) and osmolytes (GOL, GLY and to a lesser extent GABA) with a DPPC bilayer. The osmolytes or protectants are involved in membrane stabilization during common environmental stresses such as desiccation or freezing. The interaction of NTs and GOL with DPPC varied and reflects the nature of the solutes. Very low affinity was seen for the charged molecules, and an increased penetration was observed for the zwitterionic and polar molecules. Positional probability distribution and PMF profiles indicate that GABAneu, GLY and GOL penetrate deeper into the bilayer than the other NTs and were located in the vicinity of the lipid glycerol backbone. No deeper penetration was observed, which is probably due to the hydrophobic environment of the aliphatic chains and the entropic contribution from the methyl groups (tails) of the aliphatic chains. This was also reflected in relatively high energy barrier for crossing the middle of the bilayer region. The general picture that emerges for NTs is that electrostatic interactions between the positively charged amino group in NTs and the negatively charged lipid phosphate group play a major role.

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

Although these interactions were important, they were not strong enough to anchor NTs at the site of the phosphate groups. Furthermore, our findings suggest that the solutes were relatively mobile in the interfacial bilayer region. With respect to the solute– lipid interactions, the positively charged lipid choline group played only a minor role even for solutes with an anionic moiety. Analysis of the water dynamics in the bulk phase, in the interfacial bilayer region, and in the first hydration shell of solutes identified rapidly and slowly exchanging water molecules, Pfast and Pslow. The percentage of slow exchanging water molecules in the pure bilayer increases as water molecules penetrate deeper into the bilayer and follows: carbonyl oxygen atoms > phosphate and choline groups >> pure water. In other words, the water molecules are most tightly bound to the carbonyl group. Significant amounts of water co-partitioned with the solutes into the membrane interface and even relatively deeply located solutes carried 5–12 water molecules. Two residence times could also be observed for water molecules in the first hydration shell of the solutes. In the bulk phase, Pslow  20%, which increased to 40– 60% once the solutes entered the interfacial bilayer region. We found that the slow exchange of water molecules relies on the reservoir of unbounded (“free”) water molecules in the interfacial bilayer region. The results from the umbrella sampling simulations support the conclusion drawn from the equilibrium MD simulations. The position of the minima observed in the PMF profiles is in reasonable agreement with the maxima observed in the density profiles. PMF profiles for ACH and GLU show a minimum of 2– 3 kJ/mol close to the bilayer interface. GLY and GABA penetrate deeper into the bilayer than the former, and a minimum of 2– 6 kJ/mol is observed at a position close to the lipid glycerol backbone. Conflict of interest The authors declare that there are no conflicts of interest. Transparency document The Transparency document associated with this article can be found in the online version. Acknowledgments This work was supported by grants from the Lundbeck Foundation, the Danish Research Agency (Grant 272-06-0505), and the Danish National Research Foundation through the establishment of MEMPHYS, Center for Biomembrane Physics. Simulations were performed at the Danish Center for Scientific Computing at the University of Southern Denmark and the Technical University of Denmark. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.chemphyslip.2014.08.003. References Åman, K., Lindahl, E., Edholm, O., Håkansson, P., Westlund, P.-O., 2003. Structure and dynamics of interfacial water in an La phase lipid bilayer from molecular dynamics simulations. Biophys. J. 84, 102–115. Barcelo, F., Prades, J., Funari, S.S., Frau, J., Alemany, R., Escriba, P.V., 2004. The hypotensive drug 2-hydroxyoleic acid modifies the structural properties of model membranes. Mol. Membr. Biol. 21, 261–268. Bellissent-Funel, M.-C., 1998. Structure and dynamics of water near hydrophilic surfaces. J. Mol. Liq. 78, 19–28.

15

Bemporad, D., Luttmann, C., Essex, J.W., 2005. Behaviour of small solutes and large drugs in a lipid bilayer from computer simulations. Biochim. Biophys. Acta 1718, 1–21. Bemporad, D., Luttmann, C., Essex, J.W., 2004. Computer simulation of small molecule permeation across a lipid bilayer dependence on bilayer properties and solute volume size, and cross-sectional area. Biophys. J. 87, 1–13. Bennett, W.F., MacCallum, J.L., Tieleman, D.P., 2009. Thermodynamic analysis of the effect of cholesterol on dipalmitoylphosphatidylcholine lipid membranes. J. Am. Chem. Soc. 131, 1972–1978. Berquand, A., Fa, N., Dufrene, Y.F., Mingeot-Leclercq, M.-P., 2005. Interaction of the macrolide antibiotic azithromycin with lipid bilayers effect on membrane organization, fluidity, and permeability. Pharm. Res. 22, 465–475. Bhattacharyya, K., Bagchi, B., 2000. Slow dynamics of constrained water in complex geometries. J. Phys. Chem. A 104, 10603–10613. Bhide, S.Y., Berkowitz, M.L., 2005. Structure and dynamics of water at the interface with phospholipid bilayers. J. Chem. Phys. 123, 224702-1–224107-16. Boggara, M.B., Krishnamoorti, R., 2010. Partitioning of nonsteroidal antiinflammatory drugs in lipid membranes: a molecular dynamics simulation study. Biophys. J. 98, 586–595. Broemstrup, T., Reuter, N., 2010. Molecular dynamics simulations of mixed acidic/ zwitterionic phospholipid bilayers. Biophys. J. 99, 825–833. Cantor, R.S., 1997. Lateral pressures in cell membranes: a mechanism for modulation of protein function. J. Phys. Chem. B 101, 1723–1725. Cantor, R.S., 2003. Receptor Desensitization by Neurotransmitters in Membranes: Are Neurotransmitters the Endogenous Anesthetics? Biochemistry 42, 11891– 11897. Cerezo, J., Zuniga, J., Bastida, A., Requena, A., Ceron-Carrasco, J.P., 2011. Atomistic Molecular Dynamics Simulations of the Interactions of Oleic and 2-hydroxyoleic acids with phosphatidylcholine bilayers. J. Phys. Chem. B 115, 11727–11738. Debnath, A., Mukherjee, B., Ayappa, K.G., Maiti, P.K., Lin, S.-T., 2010. Entropy and dynamics of water in hydration layers of a bilayer. J. Chem. Phys. 133, 174704-1– 174704-14. Dickey, A.N., Faller, R., 2007. How alcohol chain-length and concentration modulate hydrogen bond formation in a lipid bilayer. Biophys. J. 92, 2366–2376. Edholm, O., Nagle, J.F., 2005. Areas of molecules in membranes consisting of mixtures. Biophys. J. 89, 1827–1832. Essmann, U., Perera, L., Berkowitz, M.L., Darden, T., Lee, H., Pedersen, L.G., 1995. A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593. Feller, S.E., Zhang, Y., Pastor, R.W., Brooks, B.R., 1995. Constant pressure molecular dynamics simulation: the Langevin piston method. J. Chem. Phys.103, 4613–4621. Flyvbjerg, H., Petersen, H.G., 1989. Error estimates on averages of correlated data. J. Chem. Phys. 91, 461–466. Foglia, F., Lawrence, M.J., Lorenz, C.D., MacLain, S.E., 2010. On the hydration of the phosphocholine headgroup in aqueous solution. J. Chem. Phys. 133, 145103-1– 145103-10. Franks, N.P., Lieb, W.R., 1982. Molecular mechanisms of general anaesthesia. Nature 300, 487–493. Franks, N.P., Lieb, W.R., 1984. Do general anaesthetics act by competitive binding to specific receptors? Nature 310, 599–601. Franks, N.P., Lieb, W.R., 1985. Mapping of general anaesthetic target sites provides a molecular basis for cutoff effects. Nature 316, 349–351. Grossfield, A., 2012. WHAM: the weighted histogram analysis method. version Version 2.0.6. Gruenbaum, S.M., Skinner, J.L., 2011. Vibrational spectroscopy of water in hydrated lipid multi-bilayers. I. Infrared spectra and ultrafast pump-probe observables. J. Chem. Phys. 135, 75101-1–075101-13. Hemmings Jr., H.C., Akabas, M.H., Goldstein, P.A., Trudell, J.R., Orser, B.A., Harrison, N.L., 2005. Emerging molecular mechanisms of general anesthetic action. Trends Pharmacol. Sci. 26, 503–510. Henin, J., Brannigan, G., Dailey, W.P., Eckenhoff, R., Klein, M.L., 2010. An atomistic model for simulations of the general anesthetic isoflurane. J. Phys. Chem. B 114, 604–612. Hidalgo, A.A., Caetano, W., Tabak, M., Oliveira, O.N., 2004. Interaction of two phenothiazine derivatives with phospholipid monolayers. Biophys. Chem. 109, 85–104. Högberg, C.-J., Lyubartsev, A.P., 2008. Effect of local anesthetic lidocaine on electrostatic properties of a lipid bilayer. Biophys. J. 94, 525–531. Humphrey, W., Dalke, A., Schulten, K., 1996. VMD – Visual molecular dynamics. J. Mol. Graphics 14, 33–38. Jensen, M.Ø., Mouritsen, O.G., Peters, G.H., 2004. The hydrophobic effect: molecular dynamics simulations of water confined between extended hydrophobic and hydroiphilic surfaces. J. Chem. Phys. 120, 9729–9744. Jensen, T.R., Jensen Ø, M., Reitzel, N., Balashev, K., Peters, G.H., Kjaer, K., Bjørnholm, T., 2003. Water in contact with extended hydrophobic surfaces: direct evidence of weak dewetting. Phys. Rev. Lett. 90, 86101-1–086101-4. Jodko-Piorecka, K., Litwinienko, G., 2013. First experimental evidence of dopamine interactions with negatively charged model biomembranes. ACS Chem. Neurosci. 4, 1114–1122. Johansson, A.C., Lindahl, E., 2008. Position-resolved free energy of solvation for amino acids in lipid membranes from molecular dynamics simulations. Proteins 70, 1332–1344. Jorgensen, W.L., Chandrasekhar, J., Medura, J.D., Impey, R.W., Klein, M.L., 1983. Comparison of simple potential models for simulating liquid water. J. Chem. Phys. 79, 926–935. Karplus, M., 2012. Molecular dynamics simulations of biomolecules. Acc. Chem. Res. 35, 321–323.

16

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17

Katz, Y., Diamond, J.M., 1974. Thermodynamic constants for nonelectrolyte partition between dimyristoyl lecithin and water. J. Membr. Biol. 17, 101–120. Klacsová, M., Bulacu, M., Kucerka, N., Uhrikova, D., Teixeira, J., Marrink, S.J., Balgavy, P., 2011. The effect of aliphatic alcohols on fluid bilayers in unilamellar DOPC vesicles—a small-angle neutron scattering and molecular dynamics study. Biochim. Biophys. Acta 1808, 2136–2146. König, S., Sackmann, E., Richter, D., Zorn, R., Carlile, C., Bayerl, T.M., 1994. Molecular dynamics of water in oriented DPPC multilayers studied by quasielastic neutron scattering and deuterium-nuclear magnetic resonance relaxation. J. Chem. Phys. 100, 3307–3316. Kucerka, N., Tristram-Nagle, S., Nagle, J.F., 2006. Closer look at structure of fully hydrated fluid phase DPPC bilayers. Biophys. J. 90, L83. LaBella, F.S., Stein, D., Queen, G., 1998. Occupation of the cytochrome P450 substrate pocket by diverse compounds at general anesthesia concentrations. Eur. J. Pharm. 358, 177–185. Li, C., Yi, M., Hu, J., Zhou, H.X., Cross, T.A., 2008. Solid-state NMR and MD simulations of the antiviral drug amantadine solubilized in DMPC bilayers. Biophys. J. 94, 1295–1302. MacCallum, J.L., Bennett, W.F., Tieleman, D.P., 2007. Partitioning of amino acid side chains into lipid bilayers: results from computer simulations and comparison to experiment. J. Gen. Physiol. 129, 371–377. MacCallum, J.L., Bennett, W.F., Tieleman, D.P., 2008. Distribution of amino acids in a lipid bilayer from computer simulations. Biophys. J. 94, 3393–3404. Makarov, V.A., Andrews, B.K., Smith, P.E., Pettitt, B.M., 2000. Residence times of water molecules in the hydration sites of myoglobin. Biophys. J. 79, 2966–2974. Marrink, S.J., Berendsen, H.J.C., 1994. Simulation of water transport through a lipid membrane. J. Phys. Chem. 98, 4155–4168. Marrink, S.J., Tieleman, D.P., van Buuren, A.R., Berendsen, H.J.C., 1996. Membranes and water: an interesting relationship. Faraday Discuss. 103, 191–200. Meyer, H., 1899. Zur Theorie der Alkoholnarkose. Arch. Exp. Pathol. Pharmacol. 42, 109–118. Milutinovic, P.S., Yang, L., Cantor 2, R.S., Sonner, S.J.M., 2007. Anesthetic-like modulation of a g-aminobutyric acid type A, strychnine-sensitive glycine and Nmethyl-D-aspartate receptors by coreleased neurotransmitters. Anesth. Analg. 105, 386–392. Mitragotri, Johnson, M.E., Blankschtein, D., Langer, R., 1999. An analysis of the size selectivity of solute partitioning diffusion and permeation across lipid bilayers. Biophys. J. 77, 1268–1283. Mojumdar, E.H., Lyubartsev, A.P., 2010. Molecular dynamics simulations of local anesthetic articaine in a lipid bilayer. Biophys. Chem. 153, 27–35. Mukhopadhyay, P., Vogel, H.J., Tieleman, D.P., 2004. Distribution of pentachlorophenol in phospholipid bilayers: a molecular dynamics study. Biophys. J. 86, 337–345. Murzyn, K., Rog, T., Jezierski, G., Takaoka, Y., Pasenkiewicz-Gierula, M., 2001. Effects of phospholipid unsaturation on the membrane/water interface: a molecular simulation study. Biophys. J. 81, 170–183. Murzyn, K., Zhao, W., Karttunen, M., Kurdziel, M., Rog, T., 2006. Dynamics of water at membrane surfaces: effect of headgroup structure. Biointerphases 1, 98–105. Nagle, J.F., Tristram-Nagle, S., 2000. Structure of lipid bilayers. Biochim. Biochim. Acta 1469, 159–195. Nagle, J.F., Zhang, S., Tristram-Nagle, S., Sun, W., Petrache, H.I., Suter, R.M., 1996. Xray structure determination of fully hydrated Ladipalmitoylphosphatidylcholine. Biophys. J. 70, 1419–1431. Norman, K.E., Nymeyer, H., 2006. Indole localization in lipid membranes revealed by molecular simulation. Biophys. J. 91, 2046–2054. Nury, H., Van Renterghem, C., Weng, Y., Tran, A., Baaden, M., Dufresne, V., Changeux, J.-P., Sonner, J.M., Delarue, M., Corringer, P.J., 2011. X-ray structures of general anaesthetics bound to a pentameric ligand-gated ion channel. Nature 469, 428– 431. Orlowski, A., Grzybek, M., Bunker, A., Pasenkiewicz-Gierula, M., Vattulainen, I., Männistö, P.T., Róg, T., 2012. Strong preferences of dopamine and L-dopa towards lipid head group: importance of lipid composition and implication for neurotransmitter metabolism. J. Neurochem. 122, 681–690. Osborne, D.G., Dunbar, J.A., Lapping, J.G., White, A.M., Kubarych, K.J., 2013. Sitespecific measurements of lipid membrane interfacial water dynamics with multidimensional infrared spectroscopy. J. Phys. Chem. B 117, 15407–15414. Overton, C.E., 1901. Studien über die Narkose Zugleich ein Beitrag zur Allgemeinen Pharmakologie. Gustav Fischer Verlag, Jena, Germany. Pandit, S.A., Chiu, S.W., Jakobsson, E., Grama, A., Scott, H.L., 2008. Cholesterol packing around lipids with saturated and unsaturated chains: a simulation study. Biophys. J 24, 6858–6865. Patra, M., Salonen, E., Terama, E., Vattulainen, I., Faller, R., Lee, B.W., Holopainen, J., Karttunen, M., 2006. Under the influence of alcohol: the effect of ethanol and methanol on lipid bilayers. Biophys. J. 90, 1121–1135. Pedersen, U.R., Peters, G.H., Schrøder, T.B., Dyre, J.C., 2010. Correlated volume-energy fluctuations of phospholipid membranes: a simulation study. J. Phys. Chem. B 114, 2114–2130. Pedersen, U.R., Peters, G.H., Westh, P., 2007. Molecular packing in 1-hexanol-DMPC bilayers studied by molecular dynamics simulation. Biophys. Chem. 125, 104– 111. Peetla, C., Stine, A., Labhasetwar, V., 2009. Biophysical interactions with model lipid membranes: applications in drug discovery and drug delivery. Mol. Pharm. 6, 1264–1276. Peters, G.H., 2004. Computer simulations: a tool for investigating the function of complex biological macromolecules. In: Svendsen, A. (Ed.), Enzyme Function-

ality: Design, Engineering, and Screening. Marcel Dekker Inc., New York, pp. 97– 147. Peters, G.H., Hansen, F.Y., Moeller, M.S., Westh, P., 2009. The effects of fatty acid inclusion in a DMPC bilayer membrane. J. Phys. Chem. B 113, 92–102. Peters, G.H., Wang, C., Cruys-Bagger, N., Velardez, G.F., Madsen, J.J., Westh, P., 2013. Binding of serotonin to lipid membranes. J. Am. Chem. Soc. 135, 2164–2171. Phillips, J.C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., Chipot, C., Skeel, R.D., Kale, L., Schulten, K., 2005. Scalable molecular dynamics with NAMD. J. Comp. Chem. 26, 1781–1802. Pitman, M.C., Suits, F., MacKerell Jr., A.D., Feller, S.E., 2004. Molecular-level organization of saturated and polysaturated fatty acids in a phosphatidylcholine bilayer containing cholesterol. Biochemistry 43, 15328. Pocivavsek, L., Gavrilov, K., Cao, K.D., Chi, E.Y., Li, D., Lin, B., Meron, M., Majewski, J., Lee, K.Y.C., 2011. Glycerol-induced membrane stiffening: the role of viscous fluid adlayers. Biophys. J. 101, 118–127. Preetha, A., Huilgol, N., Banerjee, R., 2007. Effect of fluidizing agents on paclitaxel penetration in cervical cancerous monolayer membranes. J. Membr. Biol. 219, 83–91. Richards, C.D., Martin, K., Gregory, S., Keightley, C.A., Hesketh, T.R., Smith, G.A., Warren, G.B., Metcalfe, J.C., 1978. Degenerate perturbations of protein structure as the mechanism of anesthetic action. Nature 276, 775–779. Rodgers, J.M., Webb, M., Smit, B., 2010. Alcohol solubility in a lipid bilayer: efficient grand-canonical simulation of an interfacially active molecule. J. Chem. Phys. 132, 064107-1–064107-10. Rog, T., Murzyn, K., Pasenkiewicz-Gierula, M., 2002. The dynamics of water at the phospholipid bilayer surface: a molecular dynamics simulation study. Chem. Phys. Lett. 352, 323–327. Scheidt, H.A., Huster, D., 2008. The interaction of small molecules with phospholipid membranes studied by 1H NOESY NMR under magic-angle spinning. Acta Pharmacol. Sin. 29, 35–49. Schellman, J.A., 1993. The relation between the free energy of interaction and binding. Biophys. Chem. 45, 273–279. Schröder, C., Rudas, T., Boresch, S., Steinhauser, O., 2006. Simulation studies of the protein–water interface: I. Properties at the molecular resolution. J. Chem. Phys. 124, 234907-1–234907-18. Seelig, J., Gally, G.U., Wohlgemuth, R., 1977. Orientation and flexibility of the choline head group in phosphatidylcholine bilayers. Biochim. Biophys. Acta 467, 109– 119. Sega, M., Vallauri, R., Melchionna, S., 2005. Diffusion of water in confined geometry: the case of a multilamellar bilayer. Phys. Rev. E 72, 041201-1–041201-4. Shinoda, W., Mikami, M., Baba, T., Hato, M., 2004. Molecular dynamics study on the effects of chain branching on the physical properties of lipid bilayers: 2. Permeability. J. Phys. Chem. B 108, 9346–9356. Slatter, S.J., Cox, K.J.A., Lombardi, J.V., Ho, C., Keily, M.B., Rubin, E., Stubbs, C.D., 1993. Inhibition of protein kinase C by alcohols and anaesthetics. Nature 364, 82–84. Sonne, J., Jensen, M.Ø., Hansen, F.Y., Hemmingsen, L., Peters, G.H., 2007. Reparameterization of all-atom dipalmitoylphosphatidylcholine lipid parameters enables simulation of fluid bilayers at zero tension. Biophys. J. 92, 4157– 4167. Sonner, J.M., Cantor, R.S., 2013. Molecular mechanisms of drug action: an emerging view. Annu. Rev. Biophys. 42, 143–167. Srinivas, P.R., Klein, M.L., 2004. Computational approaches to nanobiotechnology: probing the interaction of synthetic molecules with phospholipid. Nanotechnology 15, 1289–1295. Stouch, T.R., 1997. Solute transport and partitioning in lipid bilayers: molecular dynamics simulations. Prog. Colloid Polym. Sci. 103, 116–120. Tieleman, D.P., 2006. Computer simulations of transport through membranes: passive diffusion, pores, channels and transporters. Clin. Exp. Pharmacol. Physiol. 33, 893–903. Timasheff, S., 1998. Control of protein stability and reactions by weakly interacting cosolvents: the simplicity of the complicated. Adv. Protein Chem. 51, 355–432. Tu, K., Tarek, M., Klein, M.L., Scharf, D., 1998. Effects of anesthetics on the structure of a phospholipid bilayer: molecular dynamics investigation of halothane in the hydrated liquid crystal phase of dipalmitoylphosphatidylcholine. Biophys. J. 75, 2123–2134. Ulander, J., Haymet, A.D., 2003. Permeation across hydrated DPPC lipid bilayers: simulation of the titrable amphiphilic drug valproic acid. Biophys. J. 85, 3475– 3484. Vautrin, J., Barker, J.L., 2003. Presynaptic quantal plasticity: Katz’s original hypothesis revisited. Synapse 47, 184–199. Vautrin, J., Maric, D., Sukhareva, M., Schaffner, A.E., Barker, J.L., 2000. Surfaceaccessible GABA supports tonic and quantal synaptic transmission. Synapse 37, 38–55. Wanderlingh, U., D’Angelo, G., Conti Nibali, V., Crupi, C., Rifici, S., Corsaro, C., Sabatino, G., 2010. Interaction of alcohol with phospholipid membrane: NMR and XRD investigations on DPPC-hexanol system. Spectroscopy 24, 375–380. Wang, C., Ye, F., Velardez, G.F., Peters, G.H., Westh, P., 2011. Affinity of Four Polar Neurotransmitters for Lipid Bilayer Membranes. J. Phys. Chem. B 115, 196–203. Weast, R.C., Astle, M.J., Beyer, W.H. (Eds.), 1985. CRC Handbook of Chemistry and Physics. CRC Press Inc., Boca Raton, Florida, USA. Wedberg, R., Abildskov, J., Peters, G.H., 2012. Protein dynamics in organic media at varying water activity studied by molecular dynamics simulation. J. Phys. Chem. B 116, 2575–2585. Weng, Y., Yang, L., Corringer, P.J., Sonner, J.M., 2010. Anesthetic sensitivity of the Gloeobacter violaceus proton-gated ion channel. Anesth. Analg. 110, 59–63.

G.H. Peters et al. / Chemistry and Physics of Lipids 184 (2014) 7–17 Westh, P., 2003. Unilamellar DMPC vesicles in aqueous glycerol: preferential interactions and thermochemistry. Biophys. J. 84, 341–349. Xiang, T.-X., Anderson, B.D., 2002. A computer simulation of functional group contributions to free energy in water and a DPPC lipid bilayer. Biophys. J. 82, 2052–2066.

17

Yancey, P.H., Clark, M.E., Hand, S.C., Bowlus, R.D., Somero, G.N., 1982. Living with water stress: evolution of osmolyte systems. Science 217, 1214–1222.