Headgroup Mediated Water Insertion into the DPPC Bilayer: A

Mar 8, 2011 - sides of the bilayer in the extracellular and intracellular fluid. Lipid molecules contain polar headgroups and nonpolar hydro- carbon t...
0 downloads 0 Views 4MB Size
Download PDF
ARTICLE pubs.acs.org/JPCB

Headgroup Mediated Water Insertion into the DPPC Bilayer: A Molecular Dynamics Study Prithvi Raj Pandey and Sudip Roy* Physical Chemistry Division, National Chemical Laboratory, Pune-411008, India

bS Supporting Information ABSTRACT: Molecular dynamics simulation was performed on the 1,2-dipalmitoyl-sn-phosphocholine (DPPC) bilayerwater system using the GROMOS96 53a6 united atom force field. The transferability of force field was tested by reproducing the area per lipid within 3% accuracy from the experimental value. The simulation shows that water can penetrate much deeper inside the bilayer almost up to the starting point of the aliphatic chain. There is significant evidence from experiments that water goes deep in the DPPC bilayer, but it has not been reported from theoretical work. The mechanism of insertion of water deep inside the lipid bilayer is still not clear. In this report, for the first time, the mechanism of water insertion deep into the bilayer has been proposed. Water transport occurs by the headgroup and its first solvation shell. The trimethyl ammonium (NMe3) group (headgroup of DPPC) has two stable conformations at the bilayer-water interface, one outside the bilayer and another inside it. The NMe3 group has a large clustering of water around it and takes the water molecules inside the bilayer with it during its entry into the bilayer. The water molecules penetrate into the bilayer with the help of the NMe3 group present at the headgroup of DPPC and eventually form hydrogen bonds with carbonyl oxygen present deep inside the bilayer. Structural characteristics at the bilayerwater interface region are also reported.

1. INTRODUCTION S. J. Singer and G. L. Nicolson proposed in their famous fluid mosaic model1 of the biological membranes that it contains oriented globular proteins and lipids. Water is present near both sides of the bilayer in the extracellular and intracellular fluid. Lipid molecules contain polar headgroups and nonpolar hydrocarbon tails. They assemble to form a lipid bilayer. Being polar, the headgroups are hydrophilic in nature and point toward the extracellular fluid and the cytoplasm. The tails of the two layers face each other and form a hydrophobic region. Phospholipids are the major kind of lipid molecules present in the biological systems. Various kinds of phospholipids present in biological systems are phosphatidylserine, phosphatidylethanolamine, and phosphatidylcholine. They differ in the structure of their headgroups. 1,2-Dipalmitoyl-sn-phosphocholine (DPPC, Figure 1) contains a choline group at its headgroup. Since phospholipids are the building blocks of biological membranes, their structure and dynamics have been studied extensively.2-7 In these experimental studies, X-ray crystallography and NMR have largely been used and proved to be helpful techniques for elucidating the structure of various kinds of crystalline and amorphous phospholipids. This has been a subject of interest for almost two decades, but interesting questions still remain unsolved. How the headgroup of the lipid molecules in the bilayer interact with water is one among them? Or, in other words, how at the lipid-water r 2011 American Chemical Society

interface the hydrophilic headgroup of lipid molecules interact with adjacent water molecules? In this regard, NMR8-12 and neutron scattering13-15 and diffraction16 experiments have shown that water molecules permeate the bilayer interface with a steeply decreasing concentration as they proceed toward the inner hydrophobic region of the bilayer. Studies with other experimental techniques (e.g., FTIR) have also provided similar insight.44 For a quantitative understanding of the interactions of water molecules with polar headgroups at the lipid-water interface at the molecular level, molecular dynamics (MD) simulations have proved to be useful. MD simulations have also enforced the facts, just like experimental studies, that the majority of the water molecules residing in the polar region of the membrane are hydrogen (H-) bonded to the phosphate groups, with a smaller fraction binding to the carbonyl groups located deeper in the bilayer but with more extensive structural details.17-19,43,45 One of the initial reports that presented the study of the lipid-water interface on the dilauroylphosphatidylethanolamine lipid bilayer with MD simulations was by Damodaran et al.20 They calculated the velocity autocorrelation function for both lipid and water and also the orientational correlation function for water to understand the dynamics of the system and diffusive properties of the Received: September 21, 2010 Revised: February 13, 2011 Published: March 08, 2011 3155

dx.doi.org/10.1021/jp1090203 | J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B

ARTICLE

Figure 1. Stucture of the DPPC molecule with atom numbering as used in the text.

system, respectively. MD simulations have also shown that organization of water molecules at the lipid-water interface depends upon the type of lipid headgroup.21 The most critical component of any MD simulation is the force field used, which consists of a set of mathematical functions and parameters that describes the bonded and nonbonded interactions between the particles (atoms, united atoms) of the system under study. Recently, Kukol in his study22 has proposed a model for DPPC, which reproduces the area per lipid of the lipid bilayer within 3% error without the assumption of a constant surface area or the inclusion of surface pressure as provided by experimental study.26,27 In Kukol’s model, modification of the original model included in the GROMOS96 53a628 force field has been done with two changes in the topology. First, the partial charges on the lipid headgroup due to Chiu et al.33 were implemented, with a subdivision into four charge groups as suggested by Chandrasekhar et al.41 Second, ester-carbonyl carbon atom type was changed to “CHO” from “C”, resulting in an increase in the van der Waals radius for ester-carbonyl carbon atom to 0.664 nm as opposed to 0.336 nm before. These changes resulted in better values for area per lipid and also increased the penetration of water into the lipid bilayer headgroup region. Gierula et al. have studied the H-bonding dynamics between water and dimyristoylphosphatidylcholine (DMPC) using MD simulation.46 In their study, they have observed that, among all the DMPC oxygens, the largest ordering of water is around nonester phosphate oxygen but to a much lesser extent around estercarbonyl oxygens. The dynamics of water molecules at the fully hydrated 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC), 1-palmitoyl-2-oleoyl-phosphatidyl ethanolamine (POPE), and 1-palmitoyl-2-oleoyl-phosphatidylglycerol (POPG) bilayers and the effect of headgroups on this motion have been studied by Murzyn et al.31 They have shown that the headgroup plays an important role in the dynamics of water molecules at the interface. Inter- and intralipid interactions were studied by Zhao et al. for anionic palmitoyloleoylphosphatidylglycerol (POPG) lipid.32 A recent experimental study by Sovago et al.42 showed that, for lipid monolayers of lauric acid (LA), octadecyl trimethyl ammonium bromide (OTAB), 1,2-myristoyl-sn-glycerol-3-phosphoserine (DMPS), DPPC, 1,2-dipalmitoyl-sn-glycero-3-phosphoethanolamine (DPPE), and 1,2-dipamitoyl-3-trimethyl ammonium-propane (DPTAP), buried water molecules exist between the lipid headgroups and their alkyl tail. Also, they have shown that these water molecules have a weaker H-bond network than that of bulk water. As water penetrated much deeper in the bilayer has been proved by experiment and theory,22 in this work, we have studied for the first time to the best of our knowledge the mechanism of further insertion of water molecules, which cross through the headgroup region of the DPPC bilayer and remain in between the head and tail of lipid molecules. We have also studied the change in hydrogen bonding environment, using the modified topology proposed by Kukol, around different electronegative groups of DPPC molecules. This has been studied by clustering of water molecules and also the

Figure 2. DPPC-water interface.

H-bonding of the water molecules attached to non-ester carbonyl oxygens, phosphoryl oxygens, and NMe3 of the lipid molecules. We have shown that the insertion of these water molecules inside the bilayer occurs with the help of the NMe3 group present at the headgroup of DPPC. All the structural details obtained points to the mechanism as proposed in this work.

2. COMPUTATIONAL METHODS The initial structure of the DPPC bilayer and water system was taken from the last frame of the 40 ns run from Kukol’s work. Each DPPC molecule was taken as the united atom model with 50 atoms as discussed by Kukol. However, the all atomistic model of DPPC contains a total of 130 atoms. The choice of the united atom model was made, as it is the best known model at this stage for DPPC and is even better than any other atomistic models, and certainly provides faster and longer time scale simulation, and hence better statistics. The initial coordinates consisted of a DPPC bilayer having 128 DPPC molecules hydrated with 3655 water molecules (Figure 2). The previous model of DPPC based on GROMOS-53a6 has failed to reproduce the experimental data with high accuracy, which has already been shown in previous studies.55 The force field by Berger et al.23-25 has also been used for DPPC but without much success in the case of membrane protein. In this work, the DPPC bilayer was simulated for 40 ns using the GROMACS 4.0.729,30 package and the GROMOS53a628 force field, with the slightly modified topology proposed by Kukol,22 which happens to be the best parameters available for the united atom DPPC system and tested for proteins and bilayer assemblies. The SPC water model was considered for the simulation. The NPT ensemble and periodic boundary conditions were used. Temperature was kept constant at 325 K using a v-rescale thermostat for the whole simulation time. At this temperature, DPPC exists in liquid crystalline state, as its phase transition temperature is 314.5 K.34,35 Semi-isotropic pressure coupling was applied using a Berendsen barostat36 with separate coupling to the xy-plane and the z-direction (the bilayer normal). The pressure coupling time constant was 2.0 ps in order to maintain a constant pressure of 1.0 bar. For Lennard-Jones interaction, a cutoff at 1.4 nm was applied; electrostatic interactions were taken care of with the particle mesh Ewald (PME)37 method and a real space cutoff of 0.9 nm. Electrostatic interactions were treated with the PME method, which did not introduce artificial ordering like reaction field cutoff methods. 3156

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B

ARTICLE

Figure 3. Area per lipid.

The lipid molecule bonds were constrained with the LINCS algorithm. All of these parameters were in agreement with Kukol’s work. The MD simulation for 40 ns was carried out, and a trajectory was written every 0.5 ps. The trajectory was further analyzed to calculate radial distribution functions between different atoms of DPPC and water. Partial density profiles are reported as a two-dimensional plot where we have considered the interface along the z-axis. Dihedral distributions for a few specific dihedrals are reported to illustrate the mechanism of water insertion into the bilayer. All other nonconvetional analysis methods are described along with the results and discussions.

3. RESULTS AND DISCUSSION 3.1. Area per Lipid. Area per lipid is the most widely used property for characterization of lipid bilayers, as it can also be measured by experiments. Also, area per lipid is related to various other properties like lateral diffusion, membrane elasticity, etc. In the present work, the lipid bilayer was along the Z-axis. Thus, the area per lipid was calculated from the lengths of the box in the X and Y direction by the following equation ðbox length in XÞðbox length in Y Þ area per lipid ¼ no: of lipids in one lamella of bilayer ði:e:; 64Þ

ð1Þ The plot of area per lipid as a function of time for a total of 40 ns is shown in Figure 3. The trend of the area per lipid was similar to that reported in Kukol’s work. The average area per lipid for the 40 ns trajectory was calculated to be 0.64063 ( 0.02 nm2. This value is within a range of 3% of the experimental value of 0.64 nm2. This provides the validation of the force field and also shows the reproducibility of our simulation. We have also performed a self-assembly simulation starting from randomly oriented DPPC in water, using the same force field, and we have obtained the same result for area per lipid after selfassembly. 3.2. Partial Density Profile. We have plotted the partial density profile for DPPC, water, and various atoms present at the head and neck region of DPPC versus the length of

Figure 4. Partial density profile.

simulation box in the Z-direction in Figure 4. Interfaces between DPPC and water are clearly visible in the plot, one from 0.5 to 2.5 nm and the other from 4 to 6 nm approximately. A small dip in the plot of DPPC (at around 3.5 nm) shows the intersection between the hydrophobic regions of the two monolayers. The peak of nitrogen (N) appears nearest to the water molecules at the interface. This shows that N or the NMe3 group faces most of the bulk of water, which is quite obvious. For the phosphoryl group (PdO), both of the non-ester oxygens of PdO were treated equivalently. The partial density of PdO is highest among all the atoms in the head and neck region of DPPC. This is possibly because we have treated both of the phosphoryl oxygens as equivalent atoms. The difference between the positions of the peaks of N and PdO is negligible. However, with N being the outmost atom, we expect the peak of N in the partial density plot to be present more toward the bulk water at the bilayer-water interface. This issue has been clarified and 3157

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B

ARTICLE

Figure 5. RDFs for clustering of oxygens (water) with headgroup nitrogen, phosphorus, and carbonyl carbons.

discussed later. The position of the peak of O16 (atom number used as in Figure 1) appears more toward the inside of the bilayer than N and PdO. The position of the peak of O35 appears more toward the tail region of the bilayer than the others. O35 instead of being present at a difference of one carbon atom (C32) than O16 from the headgroup, the partial density plot shows that O35 is still present at the water-DPPC interface region. This shows that water penetrates deep into the DPPC bilayer. 3.3. Environment of the Headgroup (NMe3) and Its Neighboring Group (PdO). 3.3.1. Clustering of Water Molecules. Figure 5 shows the radial distribution functions (RDFs) of oxygen atoms of water molecules with carbon atoms (C15, C34) of carbonyl groups, phosphorus (P), and nitrogen (N) atoms present in DPPC. Thus, the RDFs represent the clustering of water molecules around these atoms in DPPC. Integration of the first RDF peak shows that clustering of water molecules is largest around N (as integration of its first peak has a value of 13.64), which is expected as N is located at the external part of the headgroup of DPPC. As a result, it sees most of the bulk of water. The heights of the first RDF peaks for P and C15 with oxygen of water molecules appear at similar values, whereas the integration of the first peaks, which is the measure of the average number of water molecules present in the first solvation shell, is 3.27 for P and 1.87 for C15. This shows that a considerable number of water molecules penetrated into the bilayer and there is a rapid decrease in water concentration from NMe3 to C15. This is also visible from the partial density profile, where we have observed that O35, despite being situated at the beginning of the tail of DPPC, is present in the interfacial region. The peak of C34 is smallest of all and solvated by only an average of 1.02 water molecules. This is because C34 is located deeper inside the bilayer where the hydrophobic region starts to form and is also at a difference of one carbon atom (C32) with respect to C15. Thus, clustering of water molecules around various atoms in the headgroup and neck region is in the following order as interpreted from the values of the number of water molecules present in the first solvation shell. N > P > C15 > C34 The size of the solvation shells around these atoms (as interpreted from the first minima of the peaks) are in the

Figure 6. (a) RDFs between carbonyl oxygens and non-ester oxygens of phosphorus with hydrogen of water molecules. (b) RDFs between nitrogen-O16, N-O (water), O16-H (water), and carbonyl C(C15)-O (water).

following order N > P  C15 > C34 3.3.2. Water Environment around O16 (Carbonyl Oxygen), PdO (non-ester Oxygen of Phosphorus), and O35 (Carbonyl Oxygen). Gierula et al. in their study suggested that the largest ordering of water in the DMPC bilayer is around phosphoryl non-ester oxygens (PdO). They observed less ordering around non-ester carbonyl oxygens. Also, Lopez et al.19 have obtained similar results with DMPC from a longer trajectory. We have plotted RDFs between hydrogen atoms of water with non-ester oxygens of PdO (O9 and O10), O35, and O16 (Figure 6a). The height of the first peak in the RDF plot is maximum for O16, which suggests ordering of a greater number of water molecules around O16 than O9 and O10 (non-ester oxygens of phosphorus atom). This is in direct contradiction to the results obtained by Gierula et al. The first peak of O9 and O10 of phosphorus is smaller than O16, and hence suggests a lesser ordering of water molecules around PdO than O16. As expected, the number of 3158

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B water molecules present in the first solvation shell of O35 is the least among others present at the interface. It appears to be obvious as O35 is located near the starting of the tail and is at a difference of one carbon atom (C32) with respect to O16 from the headgroup region. Also, the peak integrations indicate the same fact. Integration of first RDF peaks for O9 þ O10, O16, and O35 are 1.04, 1.49, and 0.61, respectively. Although we have observed and has also been proposed by Kukol22 that water penetrates considerably into the bilayerwater interface region, the higher number of water molecules near O16 than O9 þ O10 (non-ester oxygen of phosphorus) is quite interesting. A question arises of how much water molecules bypass the non-ester oxygens of phosphorus and take part in the solvation shell of O16 which is located deeper inside the bilayer than the non-ester oxygens of phosphorus? Since PdO sees more of the interfacial water as compared to O16, water is expected to be more around PdO as compared to O16, as observed by Gierula et al. for DMPC. Thus, the question is why do we get a higher first peak for O16 than PdO in the RDF plot? To address this issue and find the mechanism of water insertion deep inside the bilayer, we have compared the RDF plots between N and O16 (combined intermolecular and intramolecular (inter þ intra) RDF and only intermolecular (inter) RDF) with RDFs between N and water O (N-Owater), O16 and water H (O16-Hwater), and C15 and water O (C15Owater) in Figure 6b. We can readily observe that the minimum of the first RDF peak of N and water O appears at a similar distance as that of inter þ intra and only inter RDF plot between N and O16. It shows that the solvation shell of O16 merges with the clustering shell of the NMe3 group present in the same molecule or on a neighboring molecule. Similarly, the minimum of the first peak between C15 (carbonyl carbon) and water O appears within the first peak of N and O16. This shows that the size of the clustering shell of C15 is within the distance of closeness of N and O16 irrespective of whether they belong to the same DPPC molecule (intramolecular) or to neighboring DPPC molecules (intermolecular). Thus, it appears that the solvation shell of O16, which is attached to C15, merges with the first solvation shell of the NMe3 group. The peaks with respect to N appear at a higher distance than the O16 peak because N is surrounded by three CH3 (Me3) groups, and hence water molecules can arrange themselves just outside the -CH3 groups. Similar is the case with C15, where C is attached with O16 so water molecules arrange around O16, whereas O16 is naked and directly sees water without any intervening group. At a first intuition, it might look unusual that water molecules arrange around hydrophobic CH3 groups, but actually it is rather common. Nitrogen being positively charged (in DPPC) attracts the polar water molecules toward itself. The Coulombic attraction between N and water molecules outweighs the hydrophobic mismatch between -CH3 and water. Lum et al.39 and Chandler et al.40 have calculated the hydrophobic interactions between solute and water. Also, Godawat et al. have shown the density fluctuations at the interface of water and a range of hydrophilic and hydrophobic solutes.38 Integration of the peaks also shows that N has a larger first solvation shell than O16. Revisiting our question of why the first RDF peak is higher for O16 than PdO in Figure 6a, it can be suggested that, since the first solvation shell of C15 (and hence O16) merges with the first solvation shell of NMe3 (Figure 6b), also we have seen that clustering of water around N is the highest (Figure 5), it is possible that O16 shares the solvation environment of the NMe3

ARTICLE

group in certain intervals. This possibility can only be true if the NMe3 group periodically bends toward O16, and as a result of which O16 shares the solvating water molecules of NMe3. Dihedral bending of NMe3 is represented by the dihedral angle distribution plot of N-C5-C6-O7 sampled from the 40 ns trajectory in Figure 7a, where C5 and C6 are the connecting methyl groups between headgroup N and ester oxygen of the phosphoryl group. Two peaks in the plot suggest two conformations of NMe3, one of which corresponding to toward O16 and the other away from it. Krishnamurty et al.47 have shown using density functional theory (DFT) calculations that, for DMPC in the lowest energy conformation, the headgroup (choline group) is bent toward carbonyl in the neck region. Also, it has been proposed previously49-53 that the bent headgroup configuration should minimize the intramolecular electrostatic interaction in the polar region of phospholipids. In other words, the internal electrostatic energy resulting from repulsion between the phosphoryl group and carbonyl oxygen is minimized by attraction between choline and carbonyl groups. In Figure 6a, we have also observed that the first RDF peak of PdO (non-ester oxygen of phosphorus) with water H is in the region of the first RDF peak of O16 with water H. Hence, from this, we can also say that O16 (carbonyl oxygen near the phosphorus) shares the solvation shell with PdO. Leekumjorn et al. have shown by varying the percentage of DPPE in DPPC that the amine group in DPPE strongly interacts with phosphoryl and carbonyl groups through intra- or intermolecular H-bonding.54 Also, a surface sumfrequency generation spectroscopy study of buried water molecules inside phospholipid membranes by Sovago et al.42 suggested that water molecules interact with the zwitterionic lipid DPPC with its O-H group pointing toward bulk water and is positioned just below the choline group. This leads to the fact that the interfacial water may have a variety of hydrogen bonding environments. This fact is further confirmed by IR spectroscopic study of various kinds of lipids present in the biological system by Hubner et al.48 Further, in RDF plots between N and O16 (inter þ intra and inter) sampled for a 40 ns trajectory (Figure 6b), we have observed two peaks: one large peak at around 0.4 nm and another peak at around 0.85 nm. The origin of the predominant two peaks is due to the sampling of dihedrals (N-C5-C6-O7) at different conformations. There is a second dihedral C5-C6O7-P which may also have some effect on the bending of the NMe3 group toward the carbonyl oxygen (O16). The dihedral angle distributions for the above dihedral angles are plotted in the same figure (Figure 7a). We observe that the C5-C6-O7-P dihedral is sampling much more angular space than N-C5C6-O7. To understand the time span of each conformation produced by the combination of these two dihedrals (N-C5C6-O7 and C5-C6-O7-P), we have plotted the dihedral angles for a single DPPC molecule as a function of time in Figure 7b. In the same figure, we have also depicted the distance between N and O16 (carbonyl oxygen). From the figure, it is observed that indeed the dihedral related to the headgroup samples only two major angles and it is spending a maximum of 0.5 ns in each conformation, whereas the other dihedral (C5C6-O7-P) samples angles from -180 to 180° with much more rapid motion than the other. Therefore, the distance between N and O16 changes mainly due to the headgroup dihedral motion assisted by the much more frequent dihedral angle C5-C6O7-P motion. Taking these observations for one molecule into account, distributions of the C5-C6-O7-P dihedral were 3159

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B

ARTICLE

Figure 7. (a) Dihedral angle distribution N-C5-C6-O7 and C5-C6-O7-P. (b) Dihedral angles (N-C5-C6-O7 and C5-C6-O7-P) and distance between N and O16 as a function of time. (c) Conditional probability distribution of dihedral angle C5-C6-O7-P when N-C5-C6-O7 dihedral angle samples positive and negative angles. (d) Distribution of dihedral residence time for N-C5-C6-O7 and C5-C6-O7-P for both positive or negative directions.

plotted (Figure 7c) to understand the motion of the C5-C6O7-P dihedral with respect to the N-C5-C6-O7 dihedral under two conditions, first when the N-C5-C6-O7 dihedral lies between -180 and 0° (i.e., negative) and second when it lies between 0 and 180° (i.e., positive). These dihedral distributions show that the C5-C6-O7-P dihedral samples dihedral angles from -180 to 180° in both positive and negative regions of the N-C5-C6-O7 dihedral. It is also interesting to note that the C5-C6-O7-P dihedral distribution has a small peak around 90° when the N-C5-C6-O7 dihedral lies in the positive region and at -90° when the N-C5-C6-O7 dihedral lies in the negative region. This shows the selectivity of a conformation. In Figure 7b, we have observed that the N-C5-C6-O7 dihedral stays in the positive or negative region maximum for 0.5 ns for one molecule. Whereas in this time the C5-C6-O7P dihedral fluctuates very frequently and aquires a larger space. To make this fact statistically more acceptable, we have plotted the distribution of residence time of both the dihedrals (N-C5C6-O7 and C5-C6-O7-P) in the positive and negative space in Figure 7d. We recorded the time for which a dihedral stays in one direction for all DPPC molecules and counted how many dihedrals stay in a particular direction for a particular span of time. Finally, we have converted it into a distribution which is normalized by the number of frames and number of DPPC molecules. It is clearly visible from the plot that the N-C5C6-O7 dihedral distribution decays much more rapidly than the

C5-C6-O7-P dihedral. That is, the N-C5-C6-O7 dihedral stays in either the positive or negative direction for much longer time than the C5-C6-O7-P dihedral. As a result of these observations, we predict that the distances obtained for different dihedrals and dihedral angle conformations correspond to the closest distance of approach of N toward O16 of the same molecule. This means that the second peak of the RDF plot between N and O16 in Figure 6b has a contribution from both intramolecular and intermolecular distances between them. That is, at the headgroup region of the DPPC bilayer, N and O16 of same molecule and also of the neighboring molecules come close to each other. As a result, O16 not only sees the solvation sphere of N of the same DPPC molecule but also the solvation sphere of N of neighboring DPPC molecules. This fact also contributes to the higher RDF peak of O16 with water H than PdO with water H in Figure 6a. This results in a larger solvation shell for O16. A two-dimensional density map of DPPC and water (Figure 8) shows that the distribution of water is continuous at the interface region from the 40 ns simulation. This continuous environment of water indicates the mobile nature of the water molecules even deeper in the interface. Bending of the headgroup toward inside the bilayer is further proved by plotting the angle between the O7-N bond vector and the z-axis with time (Figure 9). When the angle is 90°, i.e., the vector is perpendicular to the z-axis, it is parallel to the interface. Also, an angle greater than 90° means the vector and hence the NMe3 group is inside the bilayer. In 3160

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B

ARTICLE

Figure 10. Distribution of the number of common water molecules between the solvation shells of N and O16.

Figure 8. Two-dimensional density map of the DPPC-water interface in the XZ-plane.

Figure 9. Angle between the O7-N vector and the z-axis of the box with time (90° angle is shown by a horizontal line).

Figure 9, we observe that the concerned angle is greater than 90° (shown by a straight line) for a considerable amount of time. This is also a very firm proof of bending of the NMe3 group inside the bilayer. In Figure 10, a distribution of the common water molecules between the solvation shell of NMe3 and O16 has been depicted. The distribution is made using a small part of the trajectory (2530 ns) and was normalized with respect to number of frames in that time. From each frame number of common water molecules between the two solvation shells were counted and divided by the total number of DPPC, i.e., 128 in the system. Then, the distribution of this number of common water molecules per DPPC was plotted. The plot shows one broad peak at 0.5, which proves the presence of 0.5 common water molecules between the solvation shell of NMe3 and O16 molecules. The plot acts as a

proof of our postulate for the mechanism of penetration of water molecules into the bilayer from bulk. In addition to the water insertion deep into the bilayer, we have also observed the bridging H-bonded water molecules between O16 and non-ester oxygen attached to the phosphorus atom. These water molecules simultaneously form H-bonds to both O16 and PdO. The average number of such H-bonds in each frame is 15.6. Thus, out of 128 DPPC molecules present in the system, around 16 DPPC molecules are involved in such H-bonding. To calculate this number, we counted the number of water molecules appearing within H-bonding distance (0.28 nm observed from Figure 6a) of O16 and PdO for each DPPC molecule in each frame. From this, we counted the number of water molecules common to both O16 and PdO and divided it by the number of frames. The number appears to be low with respect to the total number of DPPC molecules present in the system. It has been mentioned in the text that various experimental studies have shown that water has a different H-bonding environment around the headgroup of the lipid bilayer. Mean square displacement (MSD) of different water molecules, present initially in different regions, as a function of time has been checked. For this plot, we first recorded the water molecule present near the carbonyl group (CdO), PdO, and in bulk separately. Then, we tracked these water molecules in the Zdirection for 500 ps of the trajectory and plotted MSD in the Zdirection, XY-plane, and XYZ-directions separately, for those who stayed in their initial region with respect to the Z-axis (i.e., the axis along the interface) for more that 200 ps. Finally, we have reported a mean MSD plot out of these various MSDs plotted for the movement of water in the XY-plane and all XYZ-directions in the Supporting Infomation. In the two plots, namely, the average MSD of water in the XY-plane (SI1) and in all XYZ-directions (SI2), a clear distinction between the three regions is visible. This indicates different environments of water and its confinement near the headgroup region and also after entering the bilayer. Dynamics and entropy of such water molecules have been studied recently by Debnath et al.56 which supports our finding.

4. CONCLUSION Molecular dynamics simulation for 40 ns was performed using GROMOS96 53a6 force field parameters but with modified 3161

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B

ARTICLE

H-bonding water molecules points to the various H-bonding environments in the headgroup region of DPPC. This study also opens a scope to investigate the dynamics of such bridging H-bonds using the same force field as used in this work. Sharp peaks in the RDF plots of O16 with water H and also PdO with water H point to the fact that some ordered structure of water might be present in the continuous distribution of water at the interface. Further studies in this direction are required to study the possibility and nature of the ordering of the water molecules present at the interface. Figure 11. Cartoon depicting the overlap of the solvation shells of NMe3 and O16. Blue and red spheres represent the solvation shells of N and O16, respectively.

topology, as suggested by Kukol22 on the DPPC-water bilayer system. The simulations provided the area per lipid value for DPPC in 3% error range from experimental values. A better match to the experimental area per lipid leads to the penetration of water molecules deep into the bilayer, giving us the scope to unveil the mechanism of water insertion into the bilayer. We have observed from our study that the NMe3 group at the headgroup of DPPC moves inside and outside of the bilayer. Density plots show that water is distributed in a continuous fashion at the bilayer-water interface and water penetrates deep into the bilayer up to the start of the hydrophobic tail of DPPC molecules. Water keeps moving from the bilayer-water interface region to bulk water and vice versa. Due to this, it might be possible that, during the process of insertion of the headgroup, water molecules clustered around the NMe3 group enter into the deep interfacial region. Taking this as a postulate, we propose that, as a result of bending, the NMe3 group acts as a water carrier and pumps in and out water molecules in the interfacial region of the bilayer. The mechanism is shown in the cartoon (Figure 11) which is drawn from a snapshot taken from the simulation. The cartoon depicts the mechanism of how the headgroup bends inside the bilayer interface and overlaps with the solvation shell of carbonyl oxygen (O16). The NMe3 group bends inside the bilayer with water molecules solvating it. At the inside conformation of the NMe3 group, the solvation shells of NMe3 and O16 merge. There is a possibility of exchange of water molecules inside the bilayer-water interface with the solvation shell of headgroup atoms; i.e., water molecules clustering around the NMe3 group get hydrogen bonded to O16 or O9 and O10. This exchange is possible because, around the NMe3 group, water molecules are clustered outside the methyl groups, and hence in addition to the attraction with N due to its positive charge, there exists a hydrophobic mismatch with the methyl groups. However, when the water molecules approach near O16 or O9 or O10, they see only bare oxygen and are involved in H-bonding with it. This exchange is not necessarily between the solvation shells of NMe3 and O16 of the same molecule, but exchange between neighboring molecules is also possible. Complex motion of proteins or lipids is strongly affected by the amount of solvent interacting with them. Our mechanism may help to verify how much water is interacting during its entry into the bilayer mediated by headgroup. Hydration at the protein-lipid interface is required for enzymatic activity; therefore, the effect of water on the enzymatic activity might be dependent on the amount of water present at the interface, which can be verified with our mechanism. In addition to the mechanism proposed, observation of bridging

’ ASSOCIATED CONTENT

bS

Supporting Information. Average mean square displacement (MSD) plots for water molecules staying near CdO, PdO, and in bulk for more than 200 ps. Also table containing values of diffusion constants. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Phone: þ91 (020) 25902735. Fax: þ91 (020) 25902636.

’ ACKNOWLEDGMENT We gratefully acknowledge CSIR and NCL for financial support for the project. ’ REFERENCES (1) Singer, S. J.; Nicolson, G. L. Science 1972, 175, 720–731. (2) Hauser, H.; Pascher, I.; Pearson, R. H.; Sundell, S. Biochim. Biophys. Acta 1981, 650, 21–51. (3) Hauser, H.; Guyer, W.; Pascher, I.; Skrabal, P.; Sundell, S. Biochemistry 1980, 19, 366–373. (4) Hauser, H.; Pascher, I.; Sundell, S. Biochemistry 1988, 27, 9166–9174. (5) Hong, M.; Zimmermann Schmidt-Rohr, H. Biochemistry 1996, 35, 8335–8341. (6) Bruzik, K.; Harwood, S. J. Am. Chem. Soc. 1997, 119, 6629–6637. (7) Aussenac, F.; Laguerre, M.; Schmitter, J. M.; Dufourc, E. J. Langmuir 2003, 19, 10468–10479. (8) Yeagle, P. L.; Martin, R. B. Biochem. Biophys. Res. Commun. 1976, 69, 775. (9) Schmidt, C. F.; Barenholz, Y.; Huang, C.; Thompson, T. E. Biochemistry 1977, 16, 3948–3954. (10) Pope, J. M.; Comell, B. A. Chem. Phys. Lipids 1979, 24, 27. (11) Smaby, J. M.; Hermetter, A.; Schmid, P. C.; Paltauf, F.; Brockman, H. L. Biochemistry 1983, 22, 5808–5813. (12) Zhou, Z.; Sayer, B. G.; Hughes, D. W.; Stark, R. E.; Epand, R. M. Biophys. J. 1999, 76, 387. (13) Fitter, J.; Lechner, R. E.; Dencher, N. A. J. Phys. Chem. B 1999, 103, 8036–8050. (14) Sotomayor, C. P.; Aguilar, L. F.; Cuevas, F. J.; Helms, M. K.; Jameson, D. M. Biochemistry 2000, 39, 10928–10935. (15) Poolman, B.; Spitzer, J. J.; Wood, J. M. Biochim. Biophys. Acta 2004, 88, 1666. (16) Jacobs, R. E.; White, S. H. Biochemistry 1989, 28, 3421–3437. (17) Takaoka, Y.; M. Pasenkiewicz-Gierula, M.; Miyagawa, H.; Kitamura, K.; Y. Tamura, Y.; Kusumi, A. Biophys. J. 2000, 79, 3118. (18) Pasenkiewicz-Gierula; Takaoka, M. Y.; Miyagawa, H.; Kitamura, K.; Kusumi, A. Biophys. J. 1999, 76, 1228. (19) Lopez, C. F.; Nielsen, S. O.; Klein, M. L.; Moore, P. B. J. Phys. Chem. B 2004, 108, 6603–6610. 3162

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163

The Journal of Physical Chemistry B (20) Damodaran, K. V.; Kenneth, M. Merz, Jr.; Gaber, B. P. Biochemistry 1992, 31, 7656–7664. (21) Polyansky, A. A.; Volynsky, P. E.; Nolde, D. E.; Arseniev, A. S.; Efremov, R. G. J. Phys. Chem. B 2005, 109, 15052–15059. (22) Kukol, A. J. Chem. Theory Comput. 2009, 5, 615–626. (23) Beevers, A. J.; Kukol, A. J. Mol. Graphics Modell. 2006, 25, 226–233. (24) Liu, X.; Xu, Y.; Li, H.; Wang, X.; Jiang, H.; Barrantes, F. J. PLoS Comput. Biol. 2008, 4, e19. (25) Cuthbertson, J. M.; Bond, P. J.; Sansom, M. S. P. Biochemistry 2006, 45, 14298–14310. (26) Nagle, J. F.; Tristram-Nagle, S. Curr. Opin. Struct. Biol. 2000, 10, 474–480. (27) Kucerka, N.; Tristram-Nagle, S.; Nagle, J. F. Biophys. J. 2006, 90, L83–L85. (28) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. J. Comput. Chem. 2004, 25 (13), 1656–1676. (29) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G. J. Comput. Chem. 2005, 26, 1701–1718. (30) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435–447. (31) Murzyn, K.; Zhao, W.; Karttunen, M.; Kurdziel, M.; Rog, T. Biointerphases 2006, 1 (3), 98–105. (32) Zhao, W.; Rog, T.; Gurtovenko, A. A.; Vattulainen, I.; Karttunen, M. Biophys. J. 2007, 92 (4), 1114–1124. (33) Chiu, S. W.; Clark, M.; Subramaniam, B. S.; Scott, H. L.; Jacobson, E. Biophys. J. 1995, 69 (10), 1230–1245. (34) Bakouche, O.; Gerlier, D.; Letoffe, J. M.; Claudy, P. Biophys. J. 1986, 50 (1), 1–4. (35) Lopez, J. J.; Cascales; Otero, T. F.; Romero, A. J. F.; Camacho, L. Langmuir 2006, 22, 5818–5824. (36) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81 (8), 3684–3690. (37) Patra, M.; Karttunen, M.; Hyvonen, M. T.; Falck, E.; Lindqvist, P.; Vattulainen, I. Biophys. J. 2003, 84 (6), 3636–3645. (38) Godawat, R.; Jamadagni, S. N.; Garde, S. Proc. Natl. Acad. Sci. U. S.A. 2009, 106, 15119–15124. (39) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570–4577. (40) Chandler, D. Nature 2002, 417, 491. (41) Chandrasekhar, I.; Kastenholz, M.; Lins, R. D.; Oostenbrink, C.; Schuler, L. D.; Tieleman, D. P.; van Gunsteren, W. F. Eur. Biophys. J. 2003, 32 (1), 67–77. (42) Sovago, M.; Vartiainen, E.; Bonn, M. J. Chem Phys. 2009, 131, 161107. (43) Zhou, F.; Schulten, K. J. Phys. Chem. 1995, 99, 2194. (44) Volkov, V. V.; Palmer, D. J.; Righini, R. J. Phys Chem. B 2007, 111, 1377. (45) Damodaran, K. V.; Merz, K. M., Jr. Langmuir 1993, 9, 1179–1183. (46) Pasenkiewicz-Gierula, M.; Takaoka, Y.; Miyagawa, H.; Kitamura, K.; Kusumi, A. J. Phys. Chem. A 1997, 101, 3677–3691. (47) Krishnamurty, S.; Stefanov, M.; Mineva, T.; Bgu, S.; Devoisselle, J. M.; Goursot, A.; Zhu, R.; Salahub, D. R. J. Phys. Chem. B 2008, 112, 13433–13442. (48) Hubner, W.; Blume, A. Chem. Phys. Lipids 1998, 96, 99–123. (49) Hauser, H.; Pascher, I.; Pearson, R. H.; Sundell, S. Biochim. Biophys. Acta 1981, 650, 21–51. (50) Hauser, H.; Pascher, I.; Sundell, S. Biochemistry 1988, 27, 9166–9174. (51) HongK, M.; Zimmermann; Schmidt-Rohr, H. Biochemistry 1996, 35, 8335–8341. (52) Bruzik, K. S.; Harwood, J. S. J. Am. Chem. Soc. 1997, 119, 6629–6637. (53) Aussenac, F.; Laguerre, M.; Schmitter, Jean-Marie; Dufourc, E. J. Langmuir 2003, 19, 10468–10479. (54) Leekumjorn, S.; Sum, A. K. Biophys. J. 2006, 90 (11), 3951. (55) Chandrasekhar, I.; Bakowies, D.; Glattli, A.; Hunenberger, P.; Pereira, C.; Van Gunsteren, W. F. Mol. Simul. 2005, 31, 543–548.

ARTICLE

(56) Debnath., A.; Mukherjee, B.; Ayappa, K. G.; Maiti, P. K.; Lin, S.T. J. Chem. Phys. 2010, 133, 174704.

3163

dx.doi.org/10.1021/jp1090203 |J. Phys. Chem. B 2011, 115, 3155–3163