A Molecular Dynamics Simulation Study - American Chemical Society

4876

J. Phys. Chem. B 2009, 113, 4876–4886

GM1 Ganglioside Embedded in a Hydrated DOPC Membrane: A Molecular Dynamics Simulation Study Pa´l Jedlovszky* Laboratory of Interfaces and Nanosize Systems, Institute of Chemistry, Eo¨tVo¨s Lorańd UniVersity, Pa´zmańy Pe´ter stny. 1/a, H-1117 Budapest, Hungary, and HAS Research Group of Technical Analytical Chemistry, Szt. Gelle´rt te´r 4, H-1111 Budapest, Hungary

Marcello Sega Department of Physics, UniVersity of Trento, Via SommariVe 14, I-38050 PoVo, Trento, Italy, and Frankfurt Institute for AdVanced Studies, J. W. Goethe UniVersity, Ruth-Moufang Str. 1, D-60438 Frankfurt, Germany

Renzo Vallauri Department of Physics, UniVersity of Trento, Via SommariVe 14, I-38050 PoVo, Trento, Italy ReceiVed: September 15, 2008; ReVised Manuscript ReceiVed: January 18, 2009

A long molecular dynamic simulation of a fully hydrated DOPC bilayer, containing one GM1 ganglioside molecule embedded in each of the two leaflets, has been performed. The location and conformation of the GM1 molecules as well as their effect on the properties of the membrane are investigated in detail. The simulation results reveal that the GM1 molecules are present in two equilibrium arrangements, differing in the orientation of one of their two headgroup branches. The existence of these two equilibrium arrangements of GM1 in the membrane is clearly demonstrated, although their relative population, and hence their free energy difference, cannot be inferred from the present results. A condensing effect on the membrane due to the presence of the GM1 molecules is observed, and the local changes in surface density are analyzed using Voronoi polygons. Although the DOPC molecules are packed more closely in proximity of the gangliosides, the analysis of the deuterium order parameter shows that the DOPC tails are less ordered when close to a GM1. 1. Introduction Phospholipid molecules are the main components of the membranes of living eukaryotic cells, constituting a matrix in which all the molecules that are involved in any communication between the living cell and its environment (e.g., signal transducing proteins, ion channels, etc.) are embedded. Thus, any interaction between living cells and the external world involves either the lipid matrix itself or other molecules that are surrounded by it. Among these molecules, gangliosides, i.e., sialic acid containing glycosphingolipids that predominantly occur in the outer leaflet of various vertebrate cells,1 are thought to be involved in various processes,1-3 such as cell growth regulation,4 cell differentiation,5 signal transduction,6 tissue development,7 and cell recognition.8-11 Thus, for instance, the binding of the cholera toxin to the headgroup of GM1 ganglioside is the initial step of the induction of cholera.8,9 The structural and dynamical properties of gangliosides in lipid membranes have been studied by a large variety of experimental methods, such as X-ray diffraction,12-16 nuclear magnetic resonance (NMR),17,18 dynamic light scattering,16 differential scanning calorimetry,13,19 fluorescence,20-23 electron paramagnetic resonance (EPR),24 and atomic force microscopy (AFM) measurements.25,26 On the other hand, the large number of existing experimental studies on ganglio* To whom correspondence should be addressed. E-mail: pali@ chem.elte.hu.

side-containing membranes are only accompanied by a handful of computer simulation investigations.27-30 In the detailed investigation of the structure of such membranes, computer simulation methods can, however, provide a particularly useful tool, since in a simulation the appropriately chosen model of the system of interest is seen at atomistic resolution, and hence these methods can provide such a detailed insight into the system that cannot be obtained by any kind of experiment. On the other hand, the relevance of the adopted model can only be proven a posteriori, by reproducing existing experimental data sufficiently well. In this respect, computer simulation methods and various experimental techniques well complement each other in the investigation of complex disordered systems, such as membranes. The simulation of a lipid bilayer is, however, a computationally far more demanding task than that of homogeneous, isotropic disordered systems. Indeed, the vast majority of the existing lipid simulation studies have been done in the past 10 years. Besides the early simulations of pure, one-component phospholipid membranes, built up by saturated phospholipids, such as dipalmitoylphosphatidylcholine (DPPC)31-37 or dimyristoylphosphatidylcholine (DMPC),38-44 unsaturated,45-48 branching,49-51 charged,52 or fluorinated53 lipid molecules, more recent studies targeted also membranes and micelles built up by two27-30,54-64 or three components,65,66 including cholesterol54-63,65 oritsderivatives,64 bileacids,65 sphingomyelin,66 organgliosides,27-30 and membranes containing various dissolved molecules of biological relevance (e.g., anesthetics,67-71 coenzymes,72 pep-

10.1021/jp808199p CCC: $40.75  2009 American Chemical Society Published on Web 03/10/2009

Simulation of GM1 in DOPC Membrane tides,73 or olygonucleotides74). Further, the free energy profile of various small, uncharged molecules of biological relevance (e.g., H2O, NH3, CO2, CO, O2, NO, etc.) has also been determined in various membranes by computer simulation methods several times.75-80 A few years ago, we performed molecular dynamics simulation of the fully hydrated membrane of GM3 ganglioside,81 analyzed its molecular level structure in detail, and compared its short-range structure with results of both small-angle81 and wide-angle X-ray scattering experiments.82 (The GM3 molecule is built up by a sphingosine and an arachic amide tail, and a headgroup consisting of a chain of three sugar rings, i.e., a glucose, a galactose, and a charged sialic acid residue.) We have also presented preliminary results of a simulation of a fully hydrated dioleylphosphatidylcholine (DOPC) bilayer containing one GM1 ganglioside molecule in each of the two layers.28 (GM1 differs from GM3 by containing two additional sugar residues, i.e., a galactosamine and another glucose, forming a branch of its headgroup.) Simulations of one single GM1 molecule embedded in one of the two layers of fully hydrated DMPC27 and DPPC29 membranes have been reported by Ray and Mukhopadhyay27 and by Patel and Balaji,29 respectively. Very recently, Patel and Balaji reported a set of molecular dynamics simulations of mixed DPPC/GM1 membranes of various compositions, including membranes containing GM1 at different concentrations in either one or both of the two layers.30 In the present study, we extended our previous simulation of the GM1-containing DOPC bilayer28 by almost an order of magnitude in time, and investigated the properties of a GM1 molecule in DOPC on a considerably larger time scale than what has been used in any previous study. (The schematic and spatial structure of the GM1 and DOPC molecules along with the numbering scheme of their atoms used throughout this paper is shown in Figure 1.) In particular, we focus our interest on the location of the various residues of the GM1 molecule along the membrane normal axis as well as on the orientation of the two branches of its headgroup relative to the axis perpendicular to the membrane. The position and conformation of gangliosides embedded in the membrane are of particular interest, because they are thought to influence the interaction with coreceptors83,84 as well as the mechanical85 and electrostatic properties of the membrane.17 We argue that the GM1 molecule has two preferred alignments in the membrane, characterized by different orientations of the neutral branch of its headgroup. Further, we also investigate the effect of the GM1 molecule on the order of the hydrocarbon tail and area per molecule on the nearby phospholipids, and address the apparent controversy between the results of various former studies at this point. The paper is organized as follows. In section 2, details of the simulations performed are given. In section 3, the obtained results are presented and discussed in detail, whereas in section 4 the main conclusions of this study are summarized. 2. Molecular Dynamics Simulation Molecular dynamics simulations of a fully hydrated membrane containing DOPC and GM1 molecules have been performed on the isothermal-isobaric (N,p,T) ensemble at 1 atm and 298 K. Each leaflet of the membrane consisted of 139 DOPC and 1 GM1 molecules, and the membrane was hydrated by 9991 water molecules. The net negative charge of GM1 was compensated by a Na+ ion in each layer. The DOPC molecules have been described by the force field introduced by Tieleman et al. for palmitoyl-oleyl phosphati-

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4877 dylcholine (POPC),86 substituting the palmitoyl tail by a second oleyl chain. For the GM1 molecules we used a modified version of the model we developed previously for the GM3 ganglioside.81 Similarly to the POPC model of Tieleman et al., this potential is also based on the GROMOS87 force field.87,88 The appropriate modifications introduced in this force field are described elsewhere.81,89 When compared to GM3, the oligosaccharide headgroup of GM1 contains two additional neutral sugar residues, i.e., a galactosamine and a glucose. The parameters describing these residues have thus been added to the force field-the galactose and glucose rings have been described in the same way as in the similar residues of the GM3 headgroup, and the parameters describing the acetamide group attached to the galactose ring have been taken from the GROMOS87 force field. The potential models used for GM1 and DOPC treat the CH, CH2, and CH3 groups as united atoms and includes Lennard-Jones and Coulomb interaction as well as bond angle bending, torsional, and Ryckaert-Bellmans terms, the explicit functional form of which is given elsewhere.90,91 The potential parameters describing the Na+ counterions have also been taken from the GROMOS87 library, whereas the water molecules have been modeled by the rigid, three-site SPC potential.92 The simulation was performed using the GROMACS molecular dynamics program package.91 A rectangular basic simulation box was used, and standard periodic boundary conditions were applied. The temperature and pressure of the system were kept constant by means of the weak coupling algorithms of Berendsen.93 The coupling parameters used in the thermostat and barostat were set to 0.1 and 1.0 ps, respectively. The lengths of the three box edges were allowed to fluctuate independent of each other, in order to decouple the equilibration of the diagonal pressure tensor components and let the internal stress tensor relax accordingly. The bond lengths of the DOPC and GM1 molecules were kept constant by means of the LINCS algorithm,94 while those of the water molecules were kept constant using the SETTLE method.95 All site-site interactions were truncated to zero beyond the cutoff distance of 9.0 Å. The long-range part of the electrostatic interaction was accounted for by employing the smooth particle mesh Ewald method,96 using a mesh spacing of 12.0 Å and spline order of 4. The initial configuration of the simulation was prepared using an equilibrated bilayer of 72 DOPC molecules.97 Four periodic images of this configuration were assembled to a larger basic box, containing 288 DOPC molecules. Since the original configuration was obtained by using a different force field, the system was briefly equilibrated in a short molecular dynamics run of constant temperature and volume. Then, five mutually neighboring DOPC molecules, together with the nearby waters, were removed from each of the two layers, and the two GM1 molecules, along with their Na+ counterions, were placed into the cavities created this way, with the axis of GM1 aligned to the membrane normal, and its headgroup sticking out of the membrane surface. The position of the two GM1 molecules on the XY plane was chosen so that their tails could not interact directly with each other. In order to remove the overlaps this procedure unavoidably resulted in, the length of the membrane normal edge of the basic box Z was modified accordingly. Further on, the energy of the system was minimized using the conjugate gradient method, and the internal degrees of freedom were further relaxed by performing a short molecular dynamics run at constant pressure, keeping the Z edge of the basic box unchanged. Finally, the Z edge length was increased to 85.0 Å, the system was filled with water molecules, its energy was minimized, and the internal degrees of freedom were further

4878 J. Phys. Chem. B, Vol. 113, No. 14, 2009

Jedlovszky et al.

Figure 1. Schematic and spatial structure of (a) the GM1 ganglioside molecule, and (b) the DOPC molecule. The numbering scheme of the atoms of these molecules used throughout this paper is also indicated.

relaxed in a short (0.6 ns long) simulation with a reduced (0.5 fs) integration time step. The configuration resulting from this procedure was considered relaxed enough to be used as the starting configuration of the simulation. In the entire course of the simulation, the integration time step of 2 fs has been used. The system was equilibrated by performing a 7 ns long molecular dynamics run. Then, in the 40 ns long production stage, 2000 sample configurations, separated from each other by 20 ps long trajectories, were saved for successive analyses, including estimation of statistical errors.

The evolution of the total potential energy of the system as well as the length of the Z (i.e., the membrane normal) axis and the XY (i.e., the membrane plane) cross-section area of the basic simulation box during the simulation are plotted in Figure 2. A snapshot of the starting configuration of the simulation as well as an equilibrium sample configuration taken out in the production stage is shown in Figure 3. The result of error analysis has been shown by displaying error bars in those plots where the difference between two curves is comparable with their oscillations, in order to assess the

Simulation of GM1 in DOPC Membrane

Figure 2. Time evolution of the potential energy of the system simulated (top), the length of the membrane normal edge Z of the basic box (middle), and the XY (i.e., membrane plane) cross-sectional area of the basic box (bottom) in the entire course of the simulation. The end of the equilibration stage is indicated by the dashed vertical line. The average equilibrium values of the cross-sectional area and Z edge length are 9320 ( 20 Å2 and 71.1 ( 0.3 Å, respectively.

physical significance of the data. The error bars represent the variance estimate of the average values, where the autocorrelation time of the observable has been taken into account by means of an automatic windowing procedure.98,99 3. Results and Discussion 3.1. Average Area per Molecule. The average crosssectional area of the basic simulation box resulted in 9320 ( 20 Å2, corresponding to the average area per molecule of 66.60 ( 0.15 Å2. However, such a naive calculation of the average area per molecule, i.e., simply dividing the average XY crosssectional area of the simulation box by the number of lipid molecules, is not appropriate in the case of more than one component, such as the present one. For instance, it might be expected that the large oligosaccharide headgroup of the GM1 molecule occupies a considerably larger area than the zwitterionic phosphatidylcholine group of DOPC. Moreover, it is also likely that the DOPC molecules located at the immediate vicinity of a GM1 might undergo conformational changes that lead them to occupy a different area than that of the DOPC molecules being far from gangliosides. Hence, a partitioning of the total surface area of the system between different components is needed, although there is no unique and obvious way to accomplish this task. Several different methods have indeed been proposed in the recent years for solving this problem.35,60,61 One possible way of doing this is to project a representative point (e.g., the center of mass) of the molecules to the XY plane of the basic box (i.e., the plane of the membrane) and calculate the average area of the Voronoi polygons (VP) corresponding to different components.35,60 For planar seeds, the Voronoi polygon of a given seed is the locus of the points that are closer to this seed than to any other one.100-102 In the present analysis, the seeds have been chosen as the projected molecular centers, and hence the area of a VP gives an estimate of the surface area that pertains to the corresponding molecule.

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4879 To distinguish between the DOPC molecules that are close to a GM1 and the ones that are far from GM1 we calculated the two-dimensional radial distribution function g2D(r) of the projections of the GM1 and DOPC centers of mass in the XY plane. The obtained function, shown in Figure 4, has its sharp first minimum at 7.75 Å, indicating the boundary of the first lateral coordination shell of DOPC molecules around a GM1. (It should be noted that, unlike in simple liquids, the second peak of the obtained g2D(r) function is as high as the first one. This particular feature reflects the fact that the system is not homogeneous and isotropic, and thus the directional preference of the interaction between the first neighbors gives rise to a peak lower than what one could expect in a simple, isotropic system.) The integration of the obtained g2D(r) function up to this minimum reveals that, on average, this shell consists of 2.2 DOPC molecules. This finding is in accordance with the recent experimental results, obtained on mixed DPPC/GM1 monolayers, evidencing condensed geometric complexes containing about 2.3 DPPC molecules per GM1.103 Therefore, in the following, we regard a DOPC molecule as being close to a GM1 (“near DOPC molecule”) if the lateral distance of their centers-of-mass is less than 7.75 Å and regard it as being far from GM1 (“far DOPC molecule”) otherwise. From this analysis, the mean VP area of the GM1 molecules resulted in 64.6 Å2. This value is somewhat (about 3%) smaller than the value of 66.6 Å2 obtained for the “far” DOPC molecules, and even slightly smaller than the mean VP area of the “near” DOPC’s of 65.6 Å2. This result might be somewhat surprising at first sight, if one considers that the GM1 headgroup is larger than that of DOPC. This difference in the size of the two headgroups is probably compensated by their different orientations relative to the membrane normal axis. Namely, at least some parts of the GM1 headgroup are expected to stick out of the plane of the membrane, while the zwitterionic headgroup of DOPC is likely to prefer to lie in the membrane plane. This point is discussed in detail in the following subsections. Nevertheless, the reduced occupancy in terms of surface area for the DOPC molecules that are located close to a GM1 is compatible with recent experimental findings on GM1/ phospholipid monolayers,104 where negative deviations from the ideal mixing behavior have been observed. 3.2. Density Profiles. The total electron and mass density profiles of the entire system as well as the mass density profile of the water, DOPC and GM1 molecules along the membrane normal axis Z are shown in Figure 5. In all the density profile plots the distances are reported with respect to the center-ofmass of the membrane, so that the middle of the bilayer corresponds to Z ) 0 Å. Since we are not aware of any experimental electron density data available for exactly this system, we compare the simulated electron density profile curve with those resulted from a recent X-ray diffraction measurements of the pure DOPC membrane.104 Due to the rather low density of GM1 molecules, however, the derived electron density profile is expected not to differ sensibly from the pure DOPC case. The comparison of the experimental and simulated electron density profiles is shown in Figure 5. The agreement between the simulated and experimental curves is excellent, giving us some degree of confidence in the results obtained in the present analysis, although care has to be taken in making strong conclusions on the model reliability by comparing electron density profiles.105 The mass density profile of the system shows four distinct regions along the membrane normal axis Z. In the aqueous phase, at |Z| g 24 Å the density of the system is equal to that


Jedlovszky et al.

Figure 3. Snapshot of (a) the starting configuration and (b) an equilibrium sample configuration of the system simulated. The atoms of the two GM1 molecules are shown by balls.

Figure 4. Two-dimensional radial distribution function of the projections of the centers of mass of the DOPC and GM1 molecules to the plane of the membrane XY. The dashed vertical line shows the limiting distance within which a DOPC molecule is regarded as being “close” to a GM1.

of bulk water. The high-density region of the hydrated headgroups extends to about |Z| ) 12 Å, where the density of the system drops below the bulk water density value. The hydro-

carbon phase of the membrane consists of a higher density region at 12 Å g |Z| g 5 Å, and a lower density region at |Z| e 5 Å. The former region contains orientationally ordered segments of the hydrocarbon chains, whereas the latter one is largely isotropic and is characterized by a large density of the chain terminal methyl groups. As it was shown by Marrink and Berendsen, these two regions of the hydrocarbon phase play markedly different roles in the permeability of the membrane.76,77 The peak-to-peak distance of the total mass density profile resulted in 36.0 Å. This value is considerably larger than the peak-to-peak distance obtained from the partial density profile of DOPC, i.e., 30.8 Å, indicating that a large amount of water is present in the high-density region of the lipid headgroups. Indeed, the density of water molecules is 0.4 g/cm3 at the position of the density peak of the entire system (Z ) ( 18.0 Å), and still 0.2 g/cm3 at the position of the density peak of the DOPC molecules (Z ) ( 15.4 Å). The interfacial width of the system can be estimated as the distance within which the water density drops from 90% to 10% of its bulk phase value along the membrane normal axis Z. This value resulted to be 11.5 Å from the present simulation. As is seen from Figure 5, all the profiles are almost perfectly symmetric to the middle of the membrane (i.e., Z ) 0 Å), with the exception of the GM1 density profile. Taking into account


J. Phys. Chem. B, Vol. 113, No. 14, 2009 4881

Figure 6. Time evolution of the position of the two GM1 O44 (solid lines), C74 (open circles), and O113 (filled circles) atoms along the membrane normal axis Z in the entire course of the simulation.

Figure 5. Top: mass (dashed line) and electron (solid line) density profile of the system simulated along the membrane normal axis Z. Experimental electron density results of Liu and Nagle,104 obtained for pure DOPC bilayer, are shown by full circles. The scale on the right refers to the electron density data. Bottom: mass density profile of the water (full circles), DOPC (dashed line), and GM1 (solid line) molecules in the system. The scale on the right refers to the GM1 density.

the fairly large simulation time, the asymmetric profile suggests that the two ganglioside molecules might be in two different metastable states, separated by a rather high free energy barrier. One feature of this density profile certainly needs particular attention, namely, the position of the GM1 density peak. In the leaflet characterized by negative Z values, this peak is located at Z ) -18 Å, whereas in the other leaflet at Z ) 12.5 Å. This indicates that the two GM1 molecules occupy different positions along the membrane normal axis: the oligosaccharide headgroup of the one at the negative side of the Z axis protrudes to the solvent phase (since the position of the GM1 density peak is located clearly farther from the membrane interior than that of the DOPC molecules). On the other hand, the headgroup of the other GM1 molecule is embedded into the headgroup region of the membrane, the density peak of this molecule being closer to the middle of the membrane than that of DOPC (see the bottom panel of Figure 5). This result is well illustrated by the equilibrium snapshot shown in Figure 3b. As is apparent, the headgroup of one of the GM1 molecules sticks out of the membrane protruding to the aqueous phase, while that of the other GM1 molecule remains within the headgroup region of the membrane. Similar behavior was observed by Patel and Balaji, who performed 11 different 11 ns long simulations of the hydrated bilayer containing 1 GM1 and 97 DPPC molecules.29 In three of these simulations, the entire GM1 molecule remained embedded in the membrane, whereas in the other eight simulations its headgroup protruded to the aqueous phase (see Figure 2 of ref 29). These authors could not explain the reason of this behavior and discarded the former three simulations as results of possible insufficient equilibration. While it is certainly true that longer simulations are needed to correctly sample the configurational space, it does not seem to be correct to simply discard the samples in the “embedded” configuration. Indeed, our simulation results, obtained for a slightly different system using markedly different starting configurations, clearly show that both states are quite stable on a 4 times longer time scale

of 40 ns. This suggests that the two conformations (“protruding” and “embedded”) might well be representative of two states with comparable Boltzmann weights at the equilibrium. This view is supported also by the fact that at the beginning of the present simulation both GM1 molecules have been placed into the membrane in such a way that their headgroups protruded into the aqueous phase (see the snapshot showing the initial configuration of the present run in Figure 3a), and one of the two GM1 molecules has, in fact, been soaked by the membrane during the simulation. On the contrary, in their study Patel and Balaji started the simulations from a configuration in which the GM1 molecule was embedded in the membrane, and thought that the total course of the simulation of 11 ns was not long enough for the GM1 headgroup to protrude to the water phase in three of their runs, therefore assuming the presence of only one stable state instead of the possibility of two statistically relevant metastable states.29 In order to further investigate this point, we have calculated the time evolution of the position of three GM1 atoms (Figure 6), namely O44 (i.e., the O atom that links the ceramide and oligosaccharide moieties, see Figure 1), O113 (i.e., the O atom that links the GalNAc4 and Glc5 residues, and C74 (on the NeuAc3 residue). As is seen, the O44 and C74 atoms explored a rather broad Z range during the total course of the simulation in both leaflets, and the |Z| ranges explored by the corresponding atoms in the two membrane layers largely overlap each other, while the O113 atoms span markedly different regions in the two leaflets. This indicates that the differences between the two sides of the GM1 density profile are related to different conformations of the headgroup. Therefore, the relevant degrees of freedom that characterize the two equilibrium arrangements in the phospholipid membrane have to be related to the orientation of the glycosidic headgroup with respect to the bilayer surface, and not only to the position of the center of mass along the Z direction. In order to investigate the structural consequences of this dual preference of the GM1 molecule, we have calculated the mass density profiles of the ceramide moiety and the oligosaccharide headgroup of the GM1 molecules separately, and also the separate mass density profiles of the five individual sugar residues of the headgroup. The obtained profiles are shown in Figure 7. It is seen that, while the density profile of the ceramide moiety is fairly symmetric, the two oligosaccharide headgroups give their density peaks at markedly different |Z| values. Thus, the asymmetry of the GM1 density profile is originated in the


Figure 7. Top: mass density profile of the GM1 ceramide group (i.e., apolar tails, solid line) and headgroup (dashed line) along the membrane normal axis Z. Middle: mass density profile of the Glc1 (solid line), Gal2 (dotted line), NeuAc3 (dashed line), GalNac4 (full circles with line), and Glc5 (open circles with line) residues of the GM1 headgroup along the membrane normal axis Z. Bottom: mass density profile of the N (full circles with line), P (open circles with line), C13 (dashed line), unsaturated C (i.e., C24, C25, C43, and C44, dotted line), and chain terminal C atoms (i.e., C51, C52, C53, and C54, solid line) of the DOPC molecules along the membrane normal axis Z.

asymmetric arrangement of their oligosaccharide headgroups. The resolution of the headgroup density peak to the contributions coming from the individual sugar residues can further refine this picture. Thus, the Glc1, Gal2, and NeuAc3 residues give their density peaks at about the same |Z| values in the two leaflets. (These peaks are located at the Z values of -10.7 and 11.0 Å (Glc1), -15 and 14.3 Å (Gal2), and -17.3 and 17.3 Å (NeuAc3), respectively.) On the other hand, there is a large difference between the |Z| position of the density peaks of the GalNAc4 and Glc5 residues in the two leaflets: the GalNAc4 peaks are at the Z values of -19.3 and 13.4 Å, whereas the density of the Glc5 residue peaks at -21.6 Å in one, and at 12.4 Å in the other leaflet. Considering the fact that the DOPC density peak is at about (15.4 Å, this result indicates that the charged Glc1-Gal2-NeuAc3 branch of the GM1 headgroup lies always inside the headgroup region of the membrane, and the neutral GalNAc4-Glc5 branch is the part of the headgroup that can protrude to the aqueous phase. Further, the observed dual preference of the GM1 position in the membrane reflects the dual preference of solely this GalNAc4-Glc5 headgroup branch, while the rest of the GM1 molecule has only one equilibrium arrangement inside the membrane. It should also be noted that the density peak of the GalNAc4 residue in the negative Z side of the membrane is bimodal: besides the main peak at Z ) -19.3 Å, it has another one at the Z value of -13.5 Å, i.e., almost exactly the same |Z| value where the peak of the GalNAc4 density is located in the other leaflet. Further, the Glc5 density peaks exhibit a shoulder in both leaflets at Z ) (16.4 Å. We want to stress that, since in this simulation the two GM1 molecules are actually sampling only a local minimum of their configurational space, the above result has to be interpreted only in a qualitative way. The length

Jedlovszky et al. of the simulation run, moreover, suggests that the free energy barrier separating the two states has to be rather high, and probably an ergodic sampling for the GM1 molecules is unattainable by means of brute force molecular dynamics. A quantitative estimation of the relative population of the two states could only be reached by performing extensive potentialof-mean-force calculations, employing possibly a suitable variant106,107 of the umbrella sampling method.108 Work in this direction is currently in progress. Finally, in order to characterize the average arrangement of the DOPC molecules along the membrane normal axis Z, we have calculated the mass density profile of their selected atoms. The obtained profiles of the N, P, C13 (i.e., the central carbon atom of the glycerol backbone), unsaturated (i.e., C24, C25, C43, and C44) and chain terminal (i.e., C52 and C54) carbon atoms of the DOPC molecules are shown in Figure 7. It is seen that the density peak of the N and P atoms are located at about the same position (i.e., at the |Z| value of 19.7 Å for the N, and 19.0 Å for the P atoms), suggesting that the dipole vector of the zwitterionic DOPC headgroup (pointing roughly from the P to the N atom) likely prefers to lie in the plane of the membrane. This point will be addressed in the following subsection. The density of the C13 atoms has a peak at |Z| ) 16.1 Å, whereas that of the unsaturated C atoms at |Z| ) 7.2 Å. As is seen from Figure 7, this |Z| value is close to the limit beyond which the GM1 headgroup atoms never penetrate. Finally, the chain terminal methyl groups are at high density within the |Z| value of 5 Å, in accordance with the conclusion drawn from the behavior of the total mass density profile of the system (Figure 5.). It should, however, be noted that the density peak of the chain terminal atoms extends to |Z| ) 15 Å, i.e., into the region of the headgroups, indicating that strongly bent or tilted arrangements of the lipid tails also occur with a small but nonvanishing probability. 3.3. Orientation of Residues Relative to the Membrane Normal. In order to obtain a deeper insight into the possible arrangements of the oligosaccharide headgroup of GM1 in the membrane, we have investigated the preferred alignments of its two branches relative to the membrane normal axis. For this purpose, we have calculated the cosine distributions of the angles γ formed by the membrane normal vector Z, pointing, by convention, to the aqueous phase, with the vectors pointing from the O44 atom (i.e., the oxygen that links the ceramide moiety to the oligosaccharide headgroup) to the C74 and to the C120 atom. These angles characterize the alignment of the Glc1-Gal2NeuAc3 and Glc1-Gal2-GalNAc4-Glc5 branches of the GM1 headgroup, respectively. The obtained distributions are shown in Figure 8. As is apparent, in the case of the Glc1-Gal2-NeuAc3 branch the distribution is unimodal, having its peak at cos γ ≈ -0.7. This value corresponds to the preferred alignment with the tilt angle of this branch of about 45° with respect to the membrane normal axis. As it has been demonstrated in the previous subsection, in spite of this tilted alignment, this branch of the GM1 headgroup is still embedded in the headgroup region of the membrane. A different picture is seen for the other (i.e., Glc1 f Glc5) branch of the GM1 headgroup. The cosine distribution is now clearly bimodal, having a peak at cos γ ) -1 and at cos γ ) 0, corresponding to the alignments in which this headgroup branch sticks straight to the aqueous phase, and in which it lies in the plane of the membrane, respectively. Such a result is also consistent with the obtained density profiles (Figure 7). The clear bimodal nature of this distribution points


Figure 8. Top: cosine distribution of the tilt angle γ formed by the Glc1-Gal2-NeuAc3 branch (i.e., the O44 f C74 vector, solid line), and Glc1-Gal2-GalNAc4-Glc5 branch (i.e., the O44 f C120 vector, dashed line) with the membrane normal vector Z pointing to the aqueous phase. The inset shows the cosine distribution of the angle between these two headgroup branches (i.e., the C74-O44-C120 angle). Bottom: cosine distribution of the tilt angle γ formed by the Glc1 (open circles), Gal2 (dash-dotted line), NeuAc3 (full circles), GalNac4 (dashed line), and Glc5 (solid line) residues of the GM1 headgroup (i.e., the O44 f O58, O58 f O96, O64 f C74, O96 f O113, and O113 f C120 vectors, respectively) with the membrane normal vector Z pointing to the aqueous phase. The GalNac4, NeuAc3, Gal2, and Glc1 data are shifted up by 0.06, 0.12, 0.18, and 0.24 units, respectively, for clarity. All the distributions have been obtained by averaging over the configurations of both GM1 molecules.

out that the angle γ, formed by the membrane normal vector Z and the vector describing the orientation of the Glc1-Gal2GalNAc4-Glc5 headgroup branch, is probably the most appropriate quantity to characterize the dual arrangement of the GM1 molecules in the membrane, and hence it could be the key parameter of any calculation that aims to estimate the free energy difference between the two preferred alignments. This bimodal character of the preferential alignment of the Glc1-Gal2-GalNAc4-Glc5 branch suggests that, due to steric reasons, the angle formed by the two headgroup branches cannot be too small, and hence the preferred tilted alignment of the charged Glc1-Gal2-NeuAc3 branch constrains the neutral Glc1Gal2-GalNAc4-Glc5 branch to certain orientations. In order to demonstrate this, we have calculated the cosine distribution of the angle γ0 formed by the two headgroup branches (i.e., by the O44 f C74 and O44 f C120 vectors). The obtained distribution, shown in the inset of Figure 8, is sharp and unimodal, having its peak at cos γ0 ≈ 0.75, indicating the strong preference of the two branches of forming an angle of about 45° with each other. In order to clarify the consequences of the observed preferred alignments of the two headgroup branches on their individual residues, we have also calculated the cosine distributions of the angles formed by the membrane normal vector Z with the vectors characterizing the alignment of the five sugar residues. We have chosen the O44 f O58, O58 f O96, O64 f C74,

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4883 O96 f O113, and O113 f C120 vectors to describe the orientation of the Glc1, Gal2, NeuAc3, GalNAc4, and Glc5 residues, respectively. The obtained distributions are also shown in Figure 8. Interestingly, the distributions obtained for the Glc1, Gal2, and NeuAc3 residues are all bimodal, having a peak at cos γ ) -1 (corresponding to the alignment sticking straight to the aqueous phase), and another, rather broad peak around the cos γ value of -0.3. This latter peak corresponds to a rather flat orientation, in which the sugar residue flips by only 15°-25° out of the plane of the membrane. Most probably, this bimodal character of the distributions of the Glc1 and Gal2 residues leads to the dual orientational preference of the entire Glc1-Gal2GalNAc4-Glc5 branch. Indeed, even the vector describing the alignment of these two groups together (i.e., the O44 f O96 vector) has a clear dual orientational preference (data not shown). Therefore, the bimodal shape of the P(cos γ) distribution of the charged NeuAc3 residue can be explained as the alignment of this residue has to be adjusted to that of the Glc1Gal2 unit in order to allow the tilted arrangement of the entire Glc1-Gal2-NeuAc3 branch. The distributions corresponding to the GalNAc4 and Glc5 residues show much less structure, being rather flat in a broad range of cos γ values, extending from -1 to about 0.5. Two small peaks can be distinguished, at the cos γ values of about -0.6 and around 0.15-0.25. The latter peak already corresponds to inward orientation of these residues. This finding indicates that the GalNAc4 and Glc5 headgroup residues are much more flexible than the other three, and hence the dual orientational preference of the Glc1-Gal2-GalNAc4-Glc5 branch is largely imposed by that of the Glc1 and Gal2 residues. In order to characterize the alignment of the DOPC headgroups and tails relative to the membrane normal, we also calculated the cosine distributions of the angles between Z and the vectors pointing from the P to the N atom (PN vector), and along the hydrocarbon tails (i.e., from C17 to C52 and from C36 to C54). In calculating these cosine distributions, we again distinguished between the “near” and “far” DOPC molecules. The obtained results are shown in Figure 9. For comparisons, the distributions characterizing the alignment of the GM1 sphingosine and arachic amide tails (defined through the C25 f C42 and C20 f C1 vectors, respectively) are also included in this figure. As is seen, the distribution of the PN vector is rather broad, having its peak at the cos γ value of 0, indicating that the PN vector, which approximately coincides with the dipole vector of the DOPC headgroup, clearly prefers to lie within the plane of the membrane. The distributions relative to the far and near molecules are very similar. Small differences, larger than the error bars, can be observed, although they do not change the qualitative picture. This finding is consistent with the conclusion drawn both from the average surface area occupied by the different molecules and from the density profiles of the P and N atoms (Figure 7). It is clear that, although the distributions are not symmetric, indicating that the N atoms are, on average, farther from the membrane interior than the P atoms, inward alignments of the PN vector (i.e., the ones corresponding to positive cosγ values) also appear with a substantial probability. The P(cos γ) distribution corresponding to the hydrocarbon tails of both molecules are very similar to each other, having their clear peak at cos γ ) 1 (corresponding to straight inward orientation), and vanishing at negative cos γ values, indicating the lack of even slightly outward oriented tails. The integration of these distributions shows that about 17% of the tails deviates less than 15° and 45% of them less than 30° from the membrane


Jedlovszky et al.

Figure 9. Top: cosine distribution of the tilt angle γ formed by the P f N vector of the DOPC molecules located near to a GM1 (dotted line) and far from GM1 (solid line) with the membrane normal vector Z pointing to the aqueous phase. Error bars are calculated as described in refs 98 and 99. Bottom: cosine distribution of the tilt angle γ of the tails (described by the C17 f C52 and C36 f C54 vectors, averaged over both tails) of the DOPC molecules that are located near to a GM1 (dotted line) and far from GM1 (solid line), as well as of the sphingosine and arachic amide tails of the GM1 molecules (described by the C25 f C42 and C20 f C1 vectors, full and open circles, respectively) relative to the membrane normal vector Z pointing to the aqueous phase.

normal axis, and this tilt angle is larger than 60° only in 14% of the cases. 3.4. Order of the Hydrocarbon Tails. The orientational order of the various segments of the hydrocarbon tails in lipid membranes is conventionally characterized through the deuterium order parameter profile of these tails. The deuterium order parameter SCD, accessible experimentally by nuclear magnetic resonance (NMR) measurements of selectively deuterated samples, is defined through the diagonal elements of the order parameter tensor as52,109

SCD )

2Sxx Syy + 3 3

(1)

where the ijth element of the order parameter tensor is

Sij )

〈3 cos ϑi cos ϑj - δij〉 2

(2)

Here δij is the Kronecker delta, the brackets 〈...〉 denote ensemble averaging, and ϑi and ϑj are the angles formed by the membrane normal axis Z with the ith and jth molecular axes, respectively. Since in the force field used in the present study the hydrogen atoms are not treated explicitly, the molecular coordinate frame fixed to the CH2 group considered is defined in the following way.52,109 Its z axis connects the two C atoms to which the CH2 group is chemically bound, axis y lies in the plane defined by these three C atoms and it is perpendicular to axis z, whereas axis x is perpendicular to the above two. Current interpretations of the existing results concerning the effect of GM1 molecules on the order of the phospholipid tails in mixed membranes are somewhat controversial. Thus, while electron paramagnetic resonance (EPR)24 and fluorescence

Figure 10. Deuterium order parameter profiles of the CH2 groups along the DOPC tails (averaged over the two tails of the DOPC molecule, top panel), GM1 sphingosine tail (middle panel), and GM1 arachic amide tail (bottom panel). The atoms are numbered from the middle of the molecules toward the end of the tails (see the text). Error bars are calculated as described in refs 98 and 99. Note that the unsaturated and chain terminal C atoms (i.e., CH and CH3 groups) are excluded from this analysis.

polarization22 studies indicated that GM1 increases the order of the membrane, studies of the membrane dynamics of intact cells using fluorescence probes led to opposite conclusions.23 Moreover, results of X-ray diffraction measurements did not show any effect of the GM1 molecules on the membrane order.12,13 Computer simulation results are not less conflicting in this respect: while in simulations of DPPC bilayers containing one single GM1 molecule Patel and Balaji found a decrease in the order parameters of the hydrocarbon chains in the neighborhood of GM1,29 in studying DPPC/GM1 mixed membranes of higher GM1 contents the same authors found that the DPPC chain order parameters increase with increasing GM1 content.30 In order to address this point, we have calculated the SCD profiles separately for the near and far DOPC molecules. In addition, the order parameters of the GM1 sphingosine and arachic amide tails have also been calculated. In these calculations, the C atoms of the chains have been numbered from the middle of the molecules toward the end of the tails. Thus, in the case of the DOPC tails the C17 and C36 atoms are the first, and the C51 and C53 atoms are the last (14th) ones. It should be noted that the unsaturated C atoms of DOPC have been omitted from this calculation. Similarly, in the GM1 sphingosine and arachic amide tails the order parameter of 12 and 18 carbon atoms have been calculated, their numbering starting from the C30 and C19 atoms, respectively. The obtained order parameter profiles are shown in Figure 10. As is seen, the profile of the DOPC molecules having a close GM1 neighbor goes clearly below that of the DOPC’s

Simulation of GM1 in DOPC Membrane that are far from GM1. The difference between the two profiles is the largest at the C atoms located farthest away from the end of the tails: in the case of the first (i.e., C17 and C36) atoms the SCD value of the near DOPC’s is 20% smaller than that of the far ones, but this difference gradually decreases approaching the end of the tails, and vanishes at the 12th (i.e., C30 and C49) carbon atoms. Interestingly, in spite of this disordering effect of the GM1 molecules on the nearby DOPC tails, the GM1 tails themselves are characterized by rather high order. Thus, the first two (i.e., C30 and C31) carbon atoms of the sphingosine tail are characterized by SCD values larger than 0.23. Moreover, this value is about 0.2 for the third to sixth (i.e., C17-C14) carbon atoms of the arachic amide tail. (In the case of DOPC the value of SCD never exceeds 0.19.) These atoms of GM1, characterized by a rather high orientational order, are located at the outer part of the hydrocarbon phase of the membrane, exactly where the disordering effect of the GM1 molecules on the nearby DOPC’s is the strongest. The present results thus suggest that the strong orientational order of the GM1 tails is either somewhat incompatible with that of the DOPC molecules, leading to a disordering effect, or that an orientational order different from the straight pointing to the middle of the membrane is induced by GM1 on the nearby segments of the neighboring DOPC tails through some specific interactions (e.g., hydrogen bonding between the ceramide OH or NH groups and the O atoms of the DOPC ester groups). Considering also the recent finding of Patel and Balaji,30 one should not exclude the possibility that, in spite of this disordering effect on the nearby lipid tails, at different concentrations the GM1 molecules could induce an overall ordering of the membrane, expressed in the higher order of the farther lipid tails. This conjecture could, however, be tested only in a systematic study of phospholipid/GM1 mixed membranes of various different compositions. Work in this direction is currently in progress. It should finally be noted that although the results of the present simulation, showing a decrease in surface area per molecule for DOPC molecules at the vicinity of a GM1, are in accordance with surface estimates from pressure-area curve measurements on GM1/phospholipid monolayers, the suggested picture103 of a decrease in area due to ordering effect of GM1 on the nearby hydrocarbon chains is not supported by the results of our present calculations of the tail order parameter. 4. Summary and Conclusions We performed long molecular dynamics simulations of one GM1 ganglioside molecule embedded in each of the two layers of a DOPC membrane, using a force field already tested in previous studies. The electron density profile of the simulated system turned out to be in an excellent agreement with the experimental data obtained by X-ray diffraction measurements.104 In order to shed some light to the effect of the presence of GM1 molecules on the properties of the membrane we focused our interest on the following aspects: (i) surface area condensing effect, (ii) position and orientation of the GM1 and DOPC headgroups, and (iii) orientational order of the lipid tails. The surface of the membrane was partitioned by constructing the Voronoi polygons obtained by projecting the center of mass of each solute molecule (i.e., GM1 and DOPC) onto the plane of the membrane. Such a construction allowed us to study the distribution of DOPC molecules around GM1. It turned out that there are, on average, 2.2 first lateral neighbors around a GM1 within the distance of about 7.75 Å, and that GM1 has a

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4885 condensing effect on the neighboring DOPC molecules, which thus occupy a smaller surface area than the farther DOPC molecules. Moreover, the area pertaining to GM1, irrespective of its bulky headgroup, is smaller than that of the DOPC molecules. This is in accordance with the finding that the DOPC headgroups lie preferentially within the plane of the membrane. The study of the GM1 conformation and residue density profiles revealed a peculiar organization pattern of these molecules in the two membrane layers. Thus, the GM1 molecules were found to have two preferred arrangements in the membrane, i.e., a “protruding” and an “embedded” one, differing in the orientation of the Glc1-Gal2-GalNAc4-Glc5 headgroup branch. In the “embedded” case, this headgroup branch is in contact with the membrane surface, whereas in the “protruding” case this branch is stuck out of the membrane to the aqueous phase. The observed dynamical behavior of selected GM1 atoms, which span a substantial range of distances, confirmed that equilibrium in the local free energy basins was attained. The relative population of these two states, however, cannot be reliably inferred from the results of the present simulations and would require the use of specific methods for the calculations of free energy profiles. Finally, we analyzed the orientational order of the hydrocarbon tails by examining the profile of the deuterium order parameters both for DOPC and GM1. Although the GM1 tails show a higher degree of ordering with respect to that of DOPC, the tails of the near DOPC molecules are less ordered than that of the far ones. This result somehow contradicts the usual picture that higher packing (i.e., lower area per molecule) is related to a higher order in the tail region. Acknowledgment. This work has been supported by the Hungarian OTKA Foundation under project Nos. 49673 and 75328, and partly by the bilateral grant of MTA and CNR. P.J. is an Eo¨tvo¨s fellow of the Balassi Institute, Hungary, which is gratefully acknowledged. The computations have been performed using the HPC facility of the Department of Physics, University of Trento. References and Notes (1) Zeller, C. B.; Marchase, R. B. Am. J. Physiol. 1992, 262, C1341. (2) Hakomori, S. I. Annu. ReV. Biochem. 1981, 50, 733. (3) Tettamanti, G.; Riboni, L. AdV. Lipid Res. 1993, 25, 235. (4) Yates, A. J.; Agudelo, J. D.; Sung, C. C. Lipids 1992, 27, 308. (5) Yu, R. K.; Bieberich, E.; Xia, T.; Zeng, G. J. Lipid Res. 2004, 45, 783. (6) Hakomori, S. I. Biochem. Soc. Trans. 1993, 21, 583. (7) Curatolo, W. Biochim. Biophys. Acta 1987, 906, 137. (8) van Heyningen, S. Science 1974, 183, 656. (9) van Heyningen, S. In Molecular Action of Toxins and Viruses; Cohen, P., van Heyningen, S., Eds.; Elsevier: New York 1982; pp 169190. (10) Lindberg, A. A.; Brown, J. E.; Strmberg, N.; Westling-Ryd, M.; Schultz, J. E.; Karlsson, K. A. J. Biol. Chem. 1987, 262, 1779. (11) Fukuta, S.; Magnani, J. L.; Twiddy, E. M.; Holmes, R. K.; Ginsburg, V. Infect. Immun. 1988, 56, 1748. (12) McIntosh, T. J.; Simon, S. A. Biochemistry 1994, 33, 10477. (13) Reed, R. A.; Shipley, G. G. Biophys. J. 1996, 70, 1363. (14) Majewski, J.; Kuhl, T. L.; Kjaer, K.; Smith, G. S. Biophys. J. 2001, 81, 2707. (15) Hayakawa, T.; Hirai, M. Eur. Biophys. J. 2002, 31, 62. (16) Hirai, M.; Iwase, H.; Hayakawa, T.; Koizumi, M.; Takahashi, H. Biophys. J. 2003, 85, 1600. (17) McDaniel, R. V.; McLaughlin, A.; Winiski, A. P.; Eisenberg, M.; McLaughlin, S. Biochemistry 1984, 23, 4618. (18) Arnulphi, C.; Levstein, P. R.; Ramia, M. E.; Martin, C. A.; Fidelio, G. D. J. Lipid Res. 1997, 38, 1412. (19) Bach, D.; Sela, B.; Miller, I. R. Chem. Phys. Lipids 1982, 31, 381. (20) Ollmann, M.; Schwartzmann, G.; Sandhoff, K.; Galla, H. J. Biochemistry 1987, 26, 5943.

4886 J. Phys. Chem. B, Vol. 113, No. 14, 2009 (21) Reed, R. A.; Mattai, J.; Shipley, G. G. Biochemistry 1987, 26, 824. (22) Hitzemann, R. J. Chem. Phys. Lipids 1987, 43, 25. (23) Ravichandra, B.; Joshi, P. G. Biophys. Chem. 1999, 76, 117. (24) Bertoli, E.; Masserini, M.; Sonnino, S.; Ghidoni, R.; Cestaro, B.; Tettamanti, G. Biochim. Biophys. Acta 1981, 647, 196. (25) Vie, V.; van Mau, N.; Lesniewska, E.; Goudonnet, J. P.; Heitz, F.; le Grimellec, C. Langmuir 1998, 14, 4574. (26) Yuan, C.; Johnston, L. J. Biophys. J. 2000, 79, 2768. (27) Roy, D.; Mukhopadhyay, C. J. Biomol. Struct. Dyn. 2002, 19, 1121. (28) Sega, M.; Jedlovszky, P.; Vallauri, R. J. Mol. Liq. 2006, 129, 86. (29) Patel, R. Y.; Balaji, P. V. Biochim. Biophys. Acta 2007, 1768, 1628. (30) Patel, R. Y.; Balaji, P. V. J. Phys. Chem. B 2008, 112, 3346. (31) Venable, R. M.; Zhang, Y.; Hardy, B. J.; Pastor, R. W. Science 1993, 262, 223. (32) Tu, K.; Tobias, D. J.; Klein, M. L. Biophys. J. 1995, 69, 2558. (33) Feller, S. E.; Zhang, Y.; Pastor, R. W. J. Chem. Phys. 1995, 103, 10267. (34) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871. (35) Shinoda, W.; Okazaki, S. J. Chem. Phys. 1998, 109, 1517. (36) Chiu, S. W.; Clark, M. M.; Jakobsson, E.; Subramaniam, S.; Scott, H. L. J. Comput. Chem. 1999, 20, 1153. (37) Åman, K.; Lindahl, E.; Edholm, O.; Håkansson, P.; Westlund, P. O. Biophys. J. 2003, 84, 102. (38) Pasenkiewitz-Gierula, M.; Takaoka, Y.; Miyagawa, H.; Kitamura, K.; Kusumi, A. Biophys. J. 1999, 76, 1228. (39) Jedlovszky, P.; Mezei, M. J. Chem. Phys. 1999, 111, 10770. (40) Pasenkiewitz-Gierula, M.; Takaoka, Y.; Miyagawa, H.; Kitamura, K.; Kusumi, A. J. Phys. Chem. A 1997, 101, 3677. (41) Zubrzycki, I. Z.; Xu, Y.; Madrid, M.; Tang, P. J. Chem. Phys. 2000, 112, 3437. (42) Jedlovszky, P.; Mezei, M. J. Phys. Chem. B 2001, 105, 3614. (43) Alinchenko, M. G.; Anikeenko, A. V.; Medvedev, N. N.; Voloshin, V. P.; Mezei, M.; Jedlovszky, P. J. Phys. Chem. B 2004, 108, 19056. (44) Jedlovszky, P.; Pa´rtay, L.; Mezei, M. J. Mol. Liq. 2007, 131132, 225. (45) Heller, H.; Schaefer, M.; Schulten, K. J. Phys. Chem. 1993, 97, 8343. (46) Hyvo¨nnen, M. T.; Rantala, T. T.; Ala-Korpela, M. Biophys. J. 1997, 73, 2907. (47) Saiz, L.; Klein, M. L. J. Chem. Phys. 2002, 116, 3052. (48) Rabinovich, A. L.; Balabaev, N. K.; Alinchenko, M. G.; Voloshin, V. P.; Medvedev, N. N.; Jedlovszky, P. J. Chem. Phys. 2005, 122, 084906. (49) Husslein, T. J.; News, D. M.; Pattnaik, P. C.; Zhong, Q.; Moore, P. B.; Klein, M. L. J. Chem. Phys. 1998, 109, 2826. (50) Shinoda, W.; Mikami, M.; Baba, T.; Hato, M. J. Phys. Chem. B 2003, 107, 14030. (51) Shinoda, W.; Mikami, M.; Baba, T.; Hato, M. J. Chem. Phys. 2004, 121, 9648. (52) Lo´pez-Cascales, J. J.; Garcıá de la Torre, J.; Marrink, S. J.; Berendsen, H. J. C. J. Chem. Phys. 1996, 104, 2713. (53) Smondyrev, A. M.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 9864. (54) Tu, K.; Klein, M. L.; Tobias, D. J. Biophys. J. 1998, 75, 2147. (55) Smondyrev, A. M.; Berkowitz, M. L. Biophys. J. 1999, 77, 2075. (56) Pasenkiewitz-Gierula, M.; Ro´g, T.; Kitamura, K.; Kusumi, A. Biophys. J. 2000, 78, 1376. (57) Chiu, S. W.; Jakobsson, E.; Scott, H. L. J. Chem. Phys. 2001, 114, 5435. (58) Hofsa¨ss, C.; Lindahl, E.; Edholm, O. Biophys. J. 2003, 84, 2192. (59) Jedlovszky, P.; Mezei, M. J. Phys. Chem. B 2003, 107, 5311. (60) Jedlovszky, P.; Medvedev, N. N.; Mezei, M. J. Phys. Chem. B 2004, 108, 465. (61) Falck, E.; Patra, M.; Karttunen, M.; Hyvo¨nnen, M. T.; Vattulainen, I. Biophys. J. 2004, 87, 1076. (62) Falck, E.; Patra, M.; Karttunen, M.; Hyvo¨nnen, M. T.; Vattulainen, I. J. Chem. Phys. 2004, 121, 12676. (63) Alinchenko, M. G.; Voloshin, V. P.; Medvedev, N. N.; Mezei, M.; Pa´rtay, L.; Jedlovszky, P. J. Phys. Chem. B 2005, 109, 16490. (64) Smondyrev, A. M.; Berkowitz, M. L. Biophys. J. 2000, 78, 1672. (65) Marrink, S. J.; Mark, A. E. Biochemistry 2002, 41, 5375. (66) Pandit, S. A.; Vasudevan, S.; Chiu, S. W.; Mashl, R. J.; Jakobsson, E.; Scott, H. L. Biophys. J. 2004, 87, 1092. (67) Bassolino-Klimas, D.; Alper, H. E.; Stouch, T. R. J. Am. Chem. Soc. 1995, 117, 4118.

Jedlovszky et al. (68) Lo´pez-Cascales, J. J.; Hernańdez-Cifre, J. G.; Garcıá de la Torre, J. J. Phys. Chem. B 1998, 102, 625. (69) Tu, K.; Tarek, M.; Klein, M. L.; Scharf, D. Biophys. J. 1998, 75, 2123. (70) Koubi, L.; Tarek, M.; Klein, M. L.; Scharf, D. Biophys. J. 2000, 78, 800. (71) Chau, P. L.; Hoang, P. N. M.; Picaud, S.; Jedlovszky, P. Chem. Phys. Lett. 2007, 438, 294. (72) So¨derha¨ll, J. A.; Laaksonen, A. J. Phys. Chem. B 2001, 105, 9308. (73) Tieleman, D. P.; Berendsen, H. J. C. Biophys. J. 1998, 74, 2786. (74) Duong, T. H.; Mehler, E. L.; Weinstein, H. J. Comput. Phys. 1999, 151, 358. (75) Bandyopadhyay, S.; Tarek, M.; Klein, M. L. J. Phys. Chem. B 1999, 103, 10075. (76) Marrink, S. J.; Berendsen, H. J. C. J. Phys. Chem. 1994, 98, 4155. (77) Marrink, S. J.; Berendsen, H. J. C. J. Phys. Chem. 1996, 100, 16729. (78) Jedlovszky, P.; Mezei, M. J. Am. Chem. Soc. 2000, 122, 5125. (79) Jedlovszky, P.; Mezei, M. J. Phys. Chem. B 2003, 107, 5322. (80) Shinoda, W.; Mikami, M.; Baba, T.; Hato, M. J. Phys. Chem. B 2004, 108, 9346. (81) Sega, M.; Vallauri, R.; Brocca, P.; Melchionna, S. J. Phys. Chem. B 2004, 108, 20322; Erratum: J. Phys. Chem. B 2005, 109, 6036. (82) Sega, M.; Vallauri, R.; Brocca, P.; Cantu`, L.; Melchionna, S. J. Phys. Chem. B 2007, 111, 10965. (83) Sonnino, S.; Cantu`, L.; Corti, M.; Acquotti, D.; Venerando, B. Chem. Phys. Lipids 1994, 71, 21. (84) Lemieux, R. U. In IUPAC Frontiers of Chemistry; Laidler, K. J., Ed.; Pergamon Press: Oxford, UK, 1982; pp 3-25. (85) van Gorkom, L. C. M.; Cheetham, J. J.; Epand, R. M. Chem. Phys. Lipids 1995, 76, 103. (86) Tieleman, D. P.; Robertson, K. M.; MacCallum, J. L.; Monticelli, L. Int. J. Quantum Chem. 2004, 100, 1071. (87) Hermans, J.; Berendsen, H. J. C.; van Guntseren, W. F.; Postma, J. P. M. Biopolymers 1984, 23, 1513. (88) van Guntseren, W. F.; Berendsen, H. J. C. Groningen Molecular Simulation (GROMOS) Library Manual; Biomos: Groningen, Germany, 1987. (89) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72, 2002. (90) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43. (91) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306. (92) Berendsen, H. J. C.; Postma, J. P. M.; van Guntseren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, C., Ed.; Reidel: Dordrecht, 1981; p 331. (93) Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (94) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463. (95) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952. (96) Essman, U.; Perera, I.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (97) Feller, S. E.; Yin, D.; Pastor, R. W.; MacKerrel, A. D. Biophys. J. 1997, 73, 2269. (98) Wolff, U. Comput. Phys. Commun. 2004, 156, 143. (99) Wolff, U. Comput. Phys. Commun. 2007, 176, 383. (100) Voronoi, G. F. J. Reine Angew. Math. 1908, 134, 198. (101) Medvedev, N. N. The Voronoi-Delaunay Method in the Structural InVestigation of Non-crystalline Systems, SB RAS: Novosibirsk, Russia, 2000 (in Russian). (102) Okabe, A.; Boots, B.; Sugihara, K.; Chiu, S. N. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams; John Wiley: Chichester, UK, 2000. (103) Frey, S. L.; Chi, E. Y.; Arratia, C.; Majewski, J.; Kjaer, K.; Lee, K. Y. C. Biophys. J. 2008, 94, 3047. (104) Liu, Y.; Nagle, J. F. Phys. ReV. E 2004, 69, 040901. (105) van Gunsteren, W. F.; Dolenc, J.; Mark, A. E. Curr. Opin. Struct. Biol. 2008, 18, 149. (106) Mezei, M. J. Comput. Phys. 1987, 68, 237. (107) Mezei, M. Mol. Simul. 1989, 3, 301. (108) Torrie, G.; Valleau, J. P. J. Comput. Phys. 1977, 23, 178. (109) Mezei, M.; Jedlovszky, P. In Methods in Molecular Biology, Vol. 400. Methods in Membrane Lipids; Dopico, A. M., Ed.; Humana Press: Totowa, NJ, 2007; pp 127-144.

JP808199P

A Molecular Dynamics Simulation Study - American Chemical Society

Recommend Documents