Subscriber access provided by University of Oklahoma
Article
Phase Separation in a Lipid/Cholesterol System: Comparison of Coarse-Grained and United-Atom Simulations Davit Hakobyan, and Andreas Heuer J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp312245y • Publication Date (Web): 07 Mar 2013 Downloaded from http://pubs.acs.org on March 18, 2013
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
TITLE Phase Separation in a Lipid/Cholesterol System: Comparison of Coarse-grained and United-atom Simulations AUTHOR NAMES. Davit Hakobyan,1 and Andreas Heuer*,1
AUTHOR ADDRESS. 1Institute of Physical Chemistry, Corrensstr. 28/30, Muenster D-48149, Germany Tel.: +49 251 83-29177 Fax: +49 251 83-29159 E-mail:
[email protected] 1 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 40
ABSTRACT. Ternary mixtures of saturated and unsaturated phospholipids and cholesterol constitute a well-known model system to study raft formation in membranes. This phenomenon is, e.g., observed in cell membranes. Here, coarse-grained (CG) and microscopic united-atom (UA) simulations are performed to investigate the phase separation of a membrane bilayer for the ternary system of saturated 16:0 (DPPC) and unsaturated 18:2 (DUPC) phospholipids and cholesterol (CHOL). The results of a 9 microsecond UA simulation demonstrate the onset of phase separation and can be compared with properties of the corresponding CG system. While sharing the structural features of phase separation in the CG model, the onset of de-mixing for the UA model is nearly two orders of magnitude slower. This factor can be rationalized by the different short-time diffusion constants. Various system properties such as order parameters, lipid – CHOL and CHOL - CHOL interactions are analyzed and compared between the UA and the CG models. From the structural perspective the phase separation process appears to be rather similar for both models, which illustrates once more the power of using CG approaches.
KEYWORDS. raft, membrane, DPPC aggregation, dilinoleoyl-phosphocholine
2 ACS Paragon Plus Environment
Page 3 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
INTRODUCTION Biological membranes are complex entities containing self assembled mixtures of different lipids and proteins. Among all these lipids, cholesterol (CHOL) is playing quite a special role in the creation of membrane rafts.1,2 The issues related to the nature and organization of rafts in natural membranes are far from being clarified and agreed upon.3,4 Since biological membranes are complex mixtures that are very difficult to analyze, many investigations are performed on model membranes containing either pure components or wellcontrolled mixtures of either two or three components.5 Direct visualization of raft-like domains in model bilayer membranes has provided a tangible proof for the coexistence of liquid-ordered (lo) and liquid disordered (ld) phases.6-8 Of particular interest are synthetic membranes containing three components: saturated phospholipids, unsaturated phospholipids and cholesterol. In these model membranes one can observe raft like domains enhanced with cholesterol and saturated phospholipids. The domain sizes range from nanoscale to microscale. Computer simulations of mixtures containing cholesterol and phospholipids employ different models that describe the molecular interactions on a different level of detail. Enormous computational demands, especially in the field of membrane simulations, frequently require the use of coarse-grained techniques. With help of the well-established coarse-grained (CG) MARTINI potential9,10 the process of raft formation between DPPC (1,2-dipalmitoyl-sn-glycero-3-phospocholine), DUPC (1,2dilinoleoyl-sn-glycero-3-phosphocholine) and CHOL could be demonstrated.11 On the other hand, united-atom (UA) simulations of lipid systems reported so far could cover only a few hundred nanoseconds (ns) due to highly intensive calculation requirements. While the MARTINI CG model merges, normally, 4 heavy (non-hydrogen) atoms into one bead, the UA model uses an atomic level description except that the C-H covalently bonded atoms are represented as a single carbon atom of a special type. Pandit et al.12 reported results of 200 ns simulations for a relatively 3 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 40
small system consisting of sphingomyelin (SM), 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) and cholesterol where only local aggregation was observed. Whether this aggregation described the initial stage of phase separation or was only a full extent of lateral aggregation, remained open. Other investigations were limited to either pre-formed “raft” domains where the system was rather static13 or were based on randomly distributed palmitoyloleoyl-phosphatidylcholine (POPC), SM and cholesterol configuration. The properties of the lo phase such as order parameters and diffusion coefficients were calculated for a mixed system.14,15 To the best of our knowledge we report for the first time the onset of lipid phase separation for a UA ternary DPPC/DUPC/CHOL system, i.e. for a system with microscopic resolution. For this purpose 9 microsecond (µs) simulation of a UA ternary DPPC/DUPC/CHOL system have been conducted. The DUPC was used as an unsaturated lipid for our system since the DPPC/DUPC/CHOL system has been shown to phase separate for 30 % CHOL at 295K temperature for the MARTINI CG model.11 By direct experiments16 as well as by interpolation of experimental results17,18 the DPPC/DOPC system with 15 % CHOL was shown to phase separate below the miscibility transition temperature of 30-35 o
C. Compared with DOPC, DUPC has a higher disorder, and, therefore, the DPPC/DUPC with 15 %
CHOL composition under equal conditions is indeed expected to phase-separate. In the present CG system (CGs) and UA system (UAs) the small cholesterol mole fraction of 0.15 is chosen to minimize the influence of cholesterol on lipid diffusion and to be sufficient for lipid de-mixing. With large-scale simulations of the UAs together with the corresponding simulations of the CGs important information about the nature of phase separation can be obtained. More specifically, two aspects are of key interest: (i) To which degree is the nature of phase separation in the CGs backed up by the more detailed information from the UA model? Of particular interest are the relative time scales in both models. (ii) How does the microscopic properties change during the process of phase separation? Please note that the terms “phase separation” and “lipid de-mixing” are used interchangeably throughout the text since 4 ACS Paragon Plus Environment
Page 5 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
the former are usually guided by lipid decomposition. The 9 µs UA simulation is supplemented by a preformed raft-like UA system (UAraft-like) to compare the properties of phase separated CGs and UAs.
METHODS For the CG simulations the MARTINI potential 2.0 was used for lipids, cholesterol and water.9 The DUPC structures for the different models are shown in Figure 1. The CG DUPC structure differs from the CG DPPC structure in the second and third bead (marked with thickened circles) of each chain. The DUPC structure for the UA model was prepared using the model of 1-palmitoyl-2-linoleyl-snglycero-3-phosphatidylcholine19 (PLPC) by symmetrically changing the parameters of its 16:0 saturated chain according to the parameters of the 18:2 di-unsaturated chain. The extensively tested force field parameters and partial charges used for DPPC in the works of Tieleman et al.20 and Berger et al.21 were also used in the UA model here. The cholesterol model parameters and topology22 were taken from the GROMACS web site (http://www.gromacs.org/Downloads/User_contributions/ Molecule_topologies). The regular water bead9 and simple point charge23 (SPC) models were used for water in the CGs and the UAs, respectively. The concentrations 0.34:0.51:0.15 were used for the DPPC/DUPC/CHOL composition for both UA and CG bilayer systems. The bilayers were kept perpendicular towards the z direction during the simulations. The x/y side lengths of the CG bilayer and the UA bilayer were ~20 and ~16 nm, respectively. The UA/CG systems were composed of 324/510 DPPC, 486/810 DUPC lipids, and 142/238 cholesterols. GROMACS version 4.5.124,25 was used for all simulations. Electrostatic interactions were calculated through the particle mesh Ewald26 and the cutoff (with the shift truncation scheme) methods for the UAs and the CGs, respectively. The cutoff distance of 1.2 nm was used for van der Waals and electrostatic interactions.
5 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 40
The systems were simulated by using an NPT ensemble with a semiisotropic pressure tensor of 1 bar and time step of 2 femtoseconds (fs) and 20 fs for the UAs and the CGs, respectively. A temperature of 300 K for the UAs was chosen to keep the latter above the room temperature at which phase separation was observed16-18 and below the miscibility transition temperature of 30-35 oC.16,27 The CGs was simulated at 295 K at which phase separation has been shown elsewhere.11 A random initial distribution of lipids in the UAs and the CGs were obtained by a 20 ns simulation at a temperature of 450 K. Bond lengths were constrained by the linear constraint solver (LINCS) algorithm.28 The results for the CGs and the UAs encompass simulation data of 12 µs and 9 µs of effective (nonscaled) times. For the UAs this is at least one order of magnitude longer than usual bilayer simulations reported so far. Five independent runs of the initial CGs were simulated with different initial velocities to evaluate the impact of stochastic effects as well as to enhance the statistics. Moreover, a de-mixed UAraft-like system of composition of 392 DPPC, 582 DUPC and 164 CHOL has been simulated. The UAraft-like system was prepared from a pure DUPC bilayer by replacing (radially from the center of each leaflet) a certain number of DUPC by DPPC and CHOL molecules which resulted in two circular DPPC-CHOL raft-like regions (surrounded by DUPC) each positioned in different leaflets and facing each other. The UAraft-like system was simulated for 900 ns with the same parameters as for the UAs. Nearest Neighbors. For the UAs and the CGs we considered the positions of the three atoms of both lipid and cholesterol shown in Figure 2 for classifying a lipid or a cholesterol as a neighbor of another lipid or cholesterol. For the UA and the CG lipids the atoms C22, C41 and C13 (Figure 2A) and atoms C2A, C2B (D2A and D2B for DUPC lipid) and GL1 (Figure 2B) have been used. For cholesterol the atoms O6, C1, C13 and the atoms ROH, R1, R2 have been chosen for the UAs (Figure 2C) and the CGs (Figure 2D), respectively.
6 ACS Paragon Plus Environment
Page 7 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Two lipids are considered neighbors if any of the distances between the pairs C13 – C13 (GL1 – GL1 for the CG model), C22/C41 – C22/C41 (C2A/D2A/C2B/D2B – C2A/D2A/C2B/D2B for the CG model) is within the distance of the first minimum of the corresponding radial distribution function (RDF). Similarly, the nearest neighbor distances for the lipid - CHOL pair are checked against the RDF reference distances of the C13 – O6 (GL1 – ROH) and C22/C41 – C1/C13 (C2A/D2A/C2B/D2B – R1/R2) pairs. The RDF pairs for the CHOL – CHOL nearest neighborhood are the O6 – O6 (ROH – ROH), C1 – C1 (R1 – R1), C13 – C13 (R2 – R2) and C1 – C13 (R1 – R2). The distances which correspond to the first minimum of the RDF at 1 µs of effective times between different pairs for the UAs and the CGs are shown in Table 1. In order to compensate the system compression due to the condensing effect of CHOL29 as well as the temperature change from initial 450K to 300K the nearest neighbor distances for other times were corrected by scaling them with the square root ratio of the system box area (in X and Y directions) at that time to the reference box area for which the RDF distances (Table 1) are derived. In this way the artificial increase of number of nearest neighbors over time is avoided. Cholesterol faces and orientations. The well known cholesterol asymmetry, due to its two methyl groups, has different ordering effect and affinity to saturated and unsaturated lipids.30-32 The smooth face of cholesterol (CHOLα) has higher affinity and ordering effect on saturated lipids such as DPPC while unsaturated lipids favor interacting with the rough face of cholesterol (CHOLβ). This discrepancy is diminished in the MARTINI cholesterol model as a natural consequence of coarsegraining. Nevertheless, in the present analysis smooth and rough faces are assigned to CHOL for the CG model. The R1 bead of CG CHOL is considered to represent the rough face as it encompasses the methyl group with the C1 atom of the corresponding UA model while P1 represent a point of the smooth face and is the projection of R1 coordinates onto the plane defined by the beads R2, R3 and R4 (for more information please refer to the legend of Figure 3). The orientation of the cholesterol face to 7 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 40
a lipid is calculated by measuring the minimal distance from the lipid plane defined by lipid atoms C13. C22, C41 and C2A, C2B, GL1 for the UA model and the CG model, respectively (Figure 2, atoms D2A, D2B and GL1 for CG DUPC), and CHOL atoms C1 and C2 (the latter is bonded to the methyl C1 atom) for the UA and R1 bead and P1 hypothetical point for the CG models, respectively. The CHOL is assumed to be facing with the rough side to a lipid when the plane–C1 (plane–R1 for the CG) distance is found to be smaller than the plane–C2 (plane–P1) distance, and is assumed to be facing with the smooth face, otherwise. The relative orientation of two cholesterols is based on the angles defined as shown in Figure 3 and explained in the corresponding legend. The tilt angle of CHOL is calculated as the angle between the line connecting the C5 and C21 atoms of CHOL for the UAs and the R2 and R4 beads of the CGs and the bilayer normal (z axis). Order Parameter. The order parameter between A and B particles of the lipid chains is calculated as S AB = (3 cos 2 θ z − 1) / 2,
where θz is the angle between the AB bond and the z direction. The order parameter of a lipid is the average over all the consecutive bonds in the chain of a CG lipid and the average of all the imaginary bonds in UA lipid chains taking A and B atoms to be every second consecutive atoms in the chains. Diffusion Coefficients. The mean square displacement (MSD) and diffusion coefficients were calculated using the g_msd tool of Gromacs24,25 taking the average values of two leaflets as it has been reported that the relative diffusion of two leaflets may have considerable impact on overall MSD and consequently on diffusion coefficients.33 Simulation trajectories for all the systems were processed after removing the center of mass translation per leaflet. The diffusion calculations started from the first 15 and 500 ns for the CGs and the UAs, respectively. The “trestart” parameter was adjusted 8 ACS Paragon Plus Environment
Page 9 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
according to the time range (from several tens to several hundreds of nanoseconds) for which the diffusion was calculated, taking into account only the linear part of the mean square displacement over the time. Graphics. The VMD package34 was used for CHOL, lipids and bilayer presentation.
RESULTS AND DISCUSSION Temporal Evolution of the UAs and the CGs. Initial and final snapshots of the UAs and one of the five CG systems are shown in Figure 4. The initial configurations correspond to equilibrium structures for a temperature of 450K where for entropic reasons phase separation is fully suppressed. The lipids in the CGs completely separate into the green lo region (DPPC) and the blue ld region (DUPC) as already reported.11 Most interestingly, also the UAs after several microseconds shows the onset of lipid de-mixing. The intermediate CG state at 100 ns time is shown to demonstrate its similarity with the final configuration of the UAs. Just from visual inspection it is clear that the phase separation in the CGs evolves much faster than in the UAs. For example, already after 1 µs the CGs basically shows full phase separation whereas for the UAs even 9 µs are too short to reach this state. In order to compare the properties of the CGs and the UAs after the phase separation the UAraft-like system was also simulated the initial and final snapshots of which are shown on the right in Figure 4. In the final snapshot the raft region is somewhat distorted compared to the initial circular form, however, no significant dissociation of the raft is observed. The CHOL molecules which were initially randomly positioned in the central DPPC region reorganize around DPPC over the simulated 900 ns having some of the CHOL molecules being migrated into the DUPC region. This suggests that even after 900 ns the CHOL flux from the central region is not fully equilibrated in the UAraft-like system, although such properties as lipid order parameters and CHOL tilt angle have shown to be around their equilibration (saturation) values. 9 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 40
The process of the phase separation for the CG model by means of the change of number of nearest DPPC – DPPC neighbors are shown in Figure 5 for all the five CG runs separately. The decrease of two DPPC neighbors and increase of four DPPC neighbors may be observed. Importantly, it might be noted that all the five independent runs demonstrate rather concordant dynamics in time which suggests that the stochastic effects (which are expected to be present if the phase separation is driven via nucleation) may be neglected Relative Time Scales for Mobility and Phase Separation. The short and long time diffusions have been already addressed numerously in the literature.14,18,35-37 In order to differentiate the long-time lipid motion from the short-time cage-motion (also known as subdiffusion) the mean square displacement (MSD) was calculated for the CGs and the UAs as shown in Figure 6A. For the MSD calculation the DUPC was selected since diffusion of DUPC is less effected from the phase separation dynamics and, thus, the logarithmic scale of MSD at long times better approximates to the straight dashed lines with α = 1, where MSD = 4D·tα. These lines are shifted from each other by a factor of 41 which means that also the diffusions of the CGs and the UAs in the mixed state differ by a factor of 41. The diffusive regime is reached around ~15 ns and ~500 ns for the CGs and the UAs, respectively. Note that the onset of diffusive dynamics scales very similarly as the diffusion coefficient (500/15 vs. 41). It is worth noting that for short times (up to an order of nanosecond) the MSD for the UAs and the CGs basically agree. However, the forward-backward rattling is much less pronounced for the CGs as reflected by the less pronounced subdiffusive regime. This results from the strongly reduced ruggedness of the energy landscape of the CGs and, correspondingly, the lower backdragging tendency. In any event, this effect leads to the different diffusion coefficients. In the diffusive regime (t > 15 ns and t > 500 ns for the CGs and the UAs, respectively) a timedependent diffusion coefficient is determined with the help of the g_msd tool program for Gromacs (Figure 6B) by calculating the MSD for different time intervals which, of course, need to be 10 ACS Paragon Plus Environment
Page 11 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
sufficiently long to rule out any effects from the short-time subdiffusive behavior. Naturally, simulation of the diffusivity for the UAs requires simulations significantly beyond 500 ns to reach the linear regime of the MSD. Thus, comparison of the diffusivities of the CGs and the UAs requires large-scale UA simulations. The diffusion coefficients of lipids in two models in a mixed state differ by a factor of 42 as shown by the dotted lines. This factor is rather sensitive to the UA diffusion values which need to be estimated as precisely as possible. Being very close to the above estimated factor of 41 is a good consistency check. In what follows we will use a factor of 40 as the ratio of CG and UA diffusivities in the mixed state. As would be expected, for both systems the diffusion of DUPC lipids in the mixed state (homogeneous mixture) is nearly the same as of the DPPC lipids. The black arrows show the values of DUPC (3.5 ± 0.1 × 10-8 cm2s-1) and DPPC (0.5 ± 0.1 × 10-8 cm2s-1) lipids for the UAraft-like system as the limiting values for the fully separated UAs. It might be noted that the scaling factor of 40 between lipid diffusions in the mixed state differs from the factor for fully separated CGs and UAraft-like system. While for DPPC lipids the factor of 40 works well, for DUPC lipids rather a factor of 20 needs to be applied. This discrepancy might be related to a rather strong condensing behavior of CHOL in the present UA force field which, however, would require a separate study for clarification. For systems in thermodynamic equilibrium typical scaling factors between CG simulation and atomistic simulation as well as experiment are mentioned to be between 2 and 10.9-11 It is important to note that lateral diffusion of lipids in a MARTINI CG bilayer system was predominantly compared with experimental results where mainly the order of the diffusion was mentioned to agree with experiment.10,11,14 Comparison of lipid lateral diffusion of the CGs with diffusion of the UAs (atomistic) system is barely available nowadays and the present analysis might appear the first such attempt. Moreover, the diffusion of CG lipids was compared with experimental results mostly for non phase-separable system or for already phase-separated system.11,14 On the contrary, here the diffusion 11 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 40
of lipids in the mixed state is estimated i.e. at stronger non-equilibrium thermodynamic conditions which need not necessarily be scaled by the same factors. Last but not least not only the diffusion of CG lipids but also of UA lipids seems to require scaling as compared to experiment. For example, in the present CGs and UAs the diffusion coefficients of lipids in the mixed state (~60 × 10-8cm2s-1 and 1.45 × 10-8cm2s-1, respectively) differ from the value measured experimentally for a similar system (for DPPC/DOPC with 10 % CHOL a homogeneous mixture is observed with diffusion coefficient of 3.9 ± 0.2 × 10-8cm2s-1)17,18 by factors of ~20 and ~0.5, respectively. The de-mixing process is governed by two factors, as reflected by Fick’s law J ∝ -D × grad(µ), where J is the diffusion flux, D is the diffusion coefficient and µ is the chemical potential. The first factor reflects the relevant time scale for elementary dynamic processes, the second factor the thermodynamic driving force for phase separation. Actually, the Cahn-Hilliard theory of spinodal decomposition can be derived by complementing Fick’s law by an appropriate expression for the chemical potential.38 Thus after scaling the time by the respective short-time diffusion constant (in respect to the time when the global diffusion starts) one is sensitive on the thermodynamic driving force. For a better comparison of the CGs and the UAs we therefore multiply the time scale of the CG simulation by the factor of 40 introduced above. In the remaining of the paper and figures the factor of 40 is applied to the CG times. Detailed Analysis of Phase Separation Process. In the next step we display some characteristic features of the phase separation. First we show how the numbers of nearest DPPC-CHOL and DPPCDPPC neighbors increase with time (see Figure 7A and 7B). Its temporal evolution clearly indicates the process of aggregation of DPPC lipids and CHOL. For example, for the UAs the number of DPPC lipids with four DPPC neighbors is roughly doubled between 1 ns and 9 µs. This reflects the status of phase separation as already seen in Figure 4 at the end of the UA simulation.
12 ACS Paragon Plus Environment
Page 13 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Furthermore, in the phase separation regime the DPPC order parameter increases with the number of nearest DPPC neighbors (Figure 7C and 7D). For times t < 102 ns the order parameter hardly depends on the neighborhood. Thus, the build-up of local order is only possible in the phase separated regime. For t > 103 ns one observes a general increase of order. The increase is strongest if a DPPC lipid is surrounded by 4 other DPPC lipids. However, even for the DPPC lipid with only one DPPC neighbor, i.e. most likely in the ld phase, a small but significant increase is visible. The combination of Figure 7B and 7D reveal two reasons for the strong increase of the average DPPC order parameter (data not shown). First, there exists a general increase in Figure 7D, second, the onset of DPPC clustering (Figure 7B) gives more and more weight to the upper curves in Figure 7D for the overall average. The fluctuation of DPPC order parameter with one DPPC neighbor observed for the phase separated CGs (Figure 7D) is a result of low statistics since in the most of the five CG runs no DPPC with only one DPPC neighbor was encountered. The DPPC order parameter values derived from the UAraft-like system and shown with arrows present the final values of the fully separated UAs. Obviously, among the observables considered so far the order parameter reflects the onset of phase separation in the UAs most sensitively. Furthermore, its dependence on time and numbers of nearest DPPC neighbors behaves very similarly as compared to the CGs. As a natural consequence of coarsegraining the calculation of the order parameters between the UA model and the CG model is somewhat different. Thus, the absolute values of the order parameters between the models are not directly comparable. However, the order parameter differences between the final UAraft-like values and the values in the mixed UA state are rather consistent with the corresponding order parameter differences of the CGs. Please note that after temporal scaling of the CG data, as suggested by the initial diffusion constants, the onset of phase-separation at times around 0.5 µs is quite similar for both systems. Thus, at least semi-quantitatively the thermodynamic driving forces are comparable for both systems. This gives 13 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 14 of 40
additional trust to the relevance of the CGs to reflect the key properties of the underlying microscopic system. Furthermore, from simple extrapolation one may estimate that a full phase separation of the UAs on the presently considered length scales of 16 nm can be expected at sub-millisecond time scales. In contrast to the DPPC order-neighborhood dependence during the phase separation for both models, the order parameter of CG DUPC shows nearly no dependence on its neighborhood and only a weak dependence is observed for the UA DUPC (Figure 8). Interestingly, the order of CG DUPC lipids increases after the phase separation. For a possible explanation of such an increase an accompanying simulation of a pure CG DUPC bilayer at 295K was performed from which the average order parameter of DUPC was estimated to be 0.154 ± 0.004 which is higher than the average CG DUPC order parameter before 500-600 ns (before ~15 ns of CG effective times) and lower than the order parameter for later times (over ~200 ns of CG effective times). The lower order parameter value before 600 ns might be explained by the disorder of DUPC after the initial high temperature dynamics (which preceded the production run) since the first 600 ns of the scaled CG times predominantly belong to the subdiffusive time range. A higher DUPC order after the phase separation than the order of a pure DUPC bilayer might be explained by the presence of a more rigid raft region which is big enough to be comparable to the system size in the de-mixed DUPC environment in the former case. Acting as a large obstacle the overall order of the system including the DUPC lipids is increased, i.e. introducing a global rather than a local change to the order parameter. The order of the UA DUPC increases before the phase separation sets in which is, likely, a result of recovering from the previous high temperature randomization dynamics. In contrast to the CG DUPC the UA DUPC lipids demonstrate a weak dependence on their DUPC neighborhood. The average order parameter of UAraft-like DUPC is 0.23 which seems reasonable since the order parameter of unsaturated lipid is supposed to decrease with time as the system de-mixes.18 This different behavior of DUPC 14 ACS Paragon Plus Environment
Page 15 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
when comparing the CGs with the UAs might be explained by the stronger ability of UA DUPC to order when mixed with DPPC/CHOL. Regardless of the CG or the UA model the average order of DUPC should increase once an obstacle (such as a raft) is introduced into the DUPC environment. However, because, on the other hand, the presence of DPPC/CHOL increases the initial order of DUPC of the UAs (in the mixed state) at larger extent than of the CGs, the consequent phase separation leads to an overall decrease of the UA DUPC order (which is roughly a sum of a decrease of the order due to isolation of DPPC/CHOL and an increase due to formation of a rigid raft region) and to an increase of the order in the CG model as compared to the DUPC order in the mixed state. CHOL – Lipid Interactions. The interaction of CHOL with both types of lipids is shown in Figure 9A. More specifically, we show the number of CHOL molecules around either DPPC or DUPC for both systems. A clear DPPC – CHOL aggregation can be observed for both systems. In comparison to DUPC the preference of CHOL to interact with saturated lipids as compared to DOPC is less apparent. For example, in a 200 ns simulation of a UA SM/DOPC/CHOL system with all equal concentrations no preference between SM and DOPC was observed for the CHOL.12 In contrast, for the DPPC/DOPC/CHOL system a DPPC-CHOL preference has been observed via MCMD simulations, showing, additionally, a direct correlation with CHOL tilt angle.39 In the following we analyze the interaction with CHOL somewhat closer by taking into account that CHOL has a smooth α-face and a less smooth β-face. In a number of simulation projects on UA models it has been reported that the α-face tends to orient towards the saturated lipid chains while the β-side prefers unsaturated chains.30-32 Furthermore the effect of higher packing ability of POPC at presence of CHOL over DPPC and DOPC was attributed to a combination of the anisotropic chain structure of POPC and the anisotropic faces of CHOL40 which might implicitly indicate the matching of the smooth α-face with the saturated POPC chain and the β-face with the unsaturated chain.
15 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 16 of 40
For the UAs we find (see Figure 9B) that the preference of the DPPC lipids for the given time range to interact with the α-face of CHOL is indeed higher than with the β-face. The same level of preference is preserved for the UAraft-like system. In contrast, the CG CHOL model shows no real affinity of either α- or β-face to the DPPC or DUPC lipid (Figure 9C). Since the R1 bead of CG, presenting the rough side of CHOL, is positioned closer to the plane of CHOL than the methyl C1 atom in the UA model, it is expected that the face specificity of CG CHOL is diminished in comparison to the UA model. Taking into account that the CGs fully separates and that the UAs also shows continues DPPC-CHOL aggregation, one might suggest that the preference of DPPC to interact with the α-face of CHOL is not significant for the de-mixing of DPPC from DUPC. CHOL – CHOL Interactions. The direct CHOL – CHOL interaction (within the first coordination shell) in a mixture with distearoyl phosphatidylcholine (DSPC) or DOPC was shown to be lower than the distribution of CHOL – CHOL in the second coordination shell for concentrations of 20 % and above indicating that cholesterols avoid being located in nearest-neighbor positions.32 The DPPC – DPPC and CHOL – DPPC aggregation shown in Figure 7A and 7B for the UAs naturally implies some type of CHOL – CHOL aggregation. However, as seen in Figure 10 the amount of nearest CHOL – CHOL decreases as compared to the initial high-temperature configuration. At the end of the simulation on average 5% of all the CHOL molecules has nearest CHOL neighbors. This has to be contrasted with an average value of 90% for the second nearest-neighbor shell which, furthermore, remains nearly constant during the whole simulation (data not shown). Thus, in agreement with previous work32 we see that for the UAs CHOL molecules avoid direct interactions. The amount of 15 % of CHOL having nearest CHOL molecules in the CGs prior to de-mixing in comparison to the rapid decrease of the CHOL – CHOL neighbors of the UAs suggests that the direct interaction of CHOL molecules for the CG model is less restricted. For the latter system an increase of the amount of nearest CHOL – CHOL pairs is observed during de-mixing. This increase is, likely, a 16 ACS Paragon Plus Environment
Page 17 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
result of an apparent crowding of CHOL molecules in the DPPC environment i.e. in a smaller area during the phase separation. The ratio of the CHOL – CHOL concentration in the nearest neighborhood to the concentration in the second nearest neighborhood was 0.21 for the mixed system while for the phase separated system the ratio reached 0.3. This increase of the ratio supports the suggestion of crowding of CHOL molecules. As some preference of CHOL to interact with DPPC has already been shown for the UAs (Figure 7A and 9A) it may be expected that some increase of nearest CHOL – CHOL neighbors would also take place in the further course of the raft formation. This, however, could not been checked against the UAraft-like system since even after 900 ns of simulation the diffusion of CHOL into the DUPC environment from the central raft region is observed i.e. a continuous decrease of CHOL – CHOL neighbors (from 56 to 41 %) indicating that the nearest CHOL – CHOL neighborhood will stabilize in microscale times. CHOL Tilt. In Figure 11 the temporal evolution of CHOL tilt angles are also shown for both systems. The CGs shows a decrease of the average tilt angle of CHOL during the phase separation. The start of the decrease of the tilt angle (at 10-20 ns) might be related to the cooperation of CHOL with the DPPC lipids since aggregation of CHOL with DPPC (Figure 7A) seems to be the only process which if not precedes then at least starts at nearly the same time as the decrease of CHOL tilt angle. This agrees with the results of MCMD simulations of a UA DPPC/DOPC/CHOL system reported elsewhere39 where it was shown, the other way around, that the DPPC/DOPC ratio through the random lipid mutation increases at lower tilt angles of CHOL which itself was controlled by an external potential. Note that the starting time of the decrease of the CHOL tilt angle is significantly smaller than the onset of DPPC – DPPC aggregation (500 – 1000 ns; see Figure 7B). Thus, it might be regarded as a precursor dynamics to the final DPPC – DPPC aggregation. The tilt angle of CHOL in the UAs behaves differently. After the initial decrease it remains nearly constant. Obviously, here the CGs and the UAs behave somewhat differently. The average CHOL tilt 17 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 18 of 40
angle (12.5o) for the UAraft-like system is shown with an arrow and presents the limiting CHOL tilt angle at fully de-mixed UAs. Second Nearest Neighborhood of CHOL. As already discussed in the context of Figure 10, nearby CHOL molecules predominantly reside in the second coordination shell.32 Here we study the distribution and relative orientations of CHOL – CHOL pairs in the second nearest neighborhood as shown in Figure 12. The horizontal axes display the angle (α or β) between the atoms of the methyl group of one CHOL and the center of mass of another CHOL (the angles and the center of masses are defined as described in the legend of Figure 3). By definition the distributions of α and β angles are identical. The distributions look rather similar between the UAs and the CGs. The results in Figure 12 clearly show that the angles of 0o and 180o are strongly suppressed. It might be noted that for both systems the CHOL-CHOL orientations of highest probabilities in the second nearest neighborhood are among the peaks of the curves i.e. the most probable orientations might be constructed roughly by combinations of α and β angles taking the values of 30, 150, 210 and 330 degrees. Importantly, this effect is very similar for both the UAs and the CGs. Additionally, it turns out that the respective values of α and β are largely uncorrelated (correlation coefficients of 0.05 and 0.002 between the angles α and β for the UAs and the CGs, respectively). These distributions basically remain constant during the whole simulation time so that they reflect very local CHOL – CHOL properties independent of the occurrence of phase separation. CG and UA Performances in Observing Raft Formation. Finally, we would like to estimate the collective CPU speed up of observing the phase separation in the CGs in comparison to the UAs of the present DPPC/DUPC/CHOL composition. The analytical estimation of the speed up from using the CG approach would roughly scale to the product of ratio of 10 of the chosen elementary time steps between the CGs and the UAs and the factor O(NUA·ln(NUA)) / O(NCG) ~ 60 which comes from the ratio of using PME and cutoff methods for the electrostatics (likely, the strongest factor which affects 18 ACS Paragon Plus Environment
Page 19 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
the calculation times), where the NUA = 152969 and NCG = 30704 are the number of particles in the UAs and the CGs, respectively.26,41,42 Taking into account the additional scaling factor of 40 used for scaling the CG times, the speed up from using the CG model to observe the raft formation for the present composition becomes ~24000. In case of parallel simulations performed on equal number of CPU cores the times which are required to simulate 1 ns are equal to 0.018 and 3.4375 hours for the CGs and the UAs, respectively. Taking into account the time scaling factor of 40 for the CG model we obtain ~7640 for the real speed up factor. The factor of 7640 is in the range of 3-4 order of magnitude of speed up of the CG model as compared to atomistic simulation techniques evaluated elsewhere.10 Around 8300 and 1815000 CPU hours (230 and 10100 of real hours) were used to simulate 12 and 9 µs of effective times for the CGs and the UAs, respectively. The tremendous difference in performance between the present CG method and UA method makes the CG approach an extremely attractive and indispensable for simulations of such large systems as membrane bilayers.
CONCLUSIONS With a large-scale simulation of ternary lipid/cholesterol mixture a clear tendency of phase separation for UA model is shown. The UA results are based on 9 µs simulation which is at least one order of magnitude longer than simulations of membrane systems reported so far. This simulation length was essential in order to grasp the onset of phase-separation for the UA system. While sharing the same propagation line of phase separation, the times needed for lipid de-mixing in the UA model differs from the times predicted by the equivalent MARTINI CG model by nearly two orders of magnitude which, however, just reflects the ratio of diffusion coefficients of the lipids between the systems for short times (in mixed states). We conclude that the effective thermodynamic driving force have to be quite similar. 19 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 20 of 40
A clear preference of CHOL – DPPC over CHOL – DUPC interaction is shown for the systems. It is further shown that for the UAs DPPC prefers to interact more with the smooth α-face rather than with the rough β-face of the CHOL, while for the DUPC no such preference is observed. For the CGs no particular preference for CHOL face was observed for either DPPC or DUPC. Since both systems demonstrate either full phase separation or an onset of phase separation it might be suggested that the preference of CHOL face does not have notable impact on the process of raft formation. Furthermore, temporal re-distribution of the cholesterols that avoid localization in the first coordination shell is demonstrated for the UAs. The CGs shows an increase of the amount of CHOL – CHOL pairs in the nearest neighborhood during the phase separation which might be related to crowding of CHOL in a smaller area in the DPPC environment. Finally, the relative orientations of CHOL – CHOL pairs in the second nearest neighborhood for the UAs and the CGs are shown where a rather similar distribution behavior is observed for both systems. Despite some certain differences between the UAs and the CGs, from the structural perspective the phase separation process appears to be rather similar across the models for the given lipid and cholesterol compositions. Given the observed increase of three orders of magnitude in simulation for the CG model, it is highly promising that the properties of the CGs are strongly correlated with those of the microscopic UAs. This strongly suggests that the CG model indeed displays a realistic scenario for the raft formation. At the same time it is interesting to observe that for some structural features also differences are observed. This may, on the one hand, indicate some intrinsic limitations of the CG model and, on the other hand, suggest that these specific features are not of primary importance for the physical understanding of the raft formation.
AUTHOR INFORMATION Corresponding Author 20 ACS Paragon Plus Environment
Page 21 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
*E-mail:
[email protected] ACKNOWLEDGEMENT The authors thank H.-J. Galla and V. Gerke for helpful discussions. This work was supported by German Research foundation (DFG) via SFB 858. The resources for the High Performance Computing (PALMA) at the University of Muenster are gratefully acknowledged.
REFERENCES (1) Simons, K.; Ikonen, E. Nature 1997, 387, 569-572. (2) Lingwood, D.; Simons, K. Science 2010, 327, 46-50. (3) Hancock, J. F. Nat. Rev. Mol. Cell Biol. 2006, 7, 456-462. (4) Pike, L. J. J. Lipid Res. 2006, 47, 1597-1598. (5) Simons, K.; Vaz, W. L. C. Annu. Rev. Biophys. Biomol. Struct. 2004, 33, 269-295. (6) Bagatolli, L. A.; Gratton, E. Biophys. J. 2000, 79, 434-447. (7) Dietrich, C.; Bagatolli, L. A.; Volovyk, Z. N.; Thompson, N. L.; Levi, M.; Jacobson, K.; Gratton, E. Biophys. J. 2001, 80, 1417-1428. (8) Dietrich, C.; Volovyk, Z. N.; Levi, M.; Thompson, N. L.; Jacobson, K. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10642-10647. (9) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. J. Phys. Chem. B 2007, 111, 7812-7824. (10) Marrink, S. J.; de Vries, A. H.; Mark, A. E. J. Phys. Chem. B 2004, 108, 750-760. (11) Risselada, H. J.; Marrink, S. J. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 17367-17372. (12) Pandit, S. A.; Jakobsson, E.; Scott, H. L. Biophys. J. 2004, 87, 3312-3322. (13) Pandit, S. A.; Vasudevan, S.; Chiu, S. W.; Mashl, R. J.; Jakobsson, E.; Scott, H. L. Biophys. J. 2004, 87, 1092-1100. (14) Apajalahti, T.; Niemela, P.; Govindan, P. N.; Miettinen, M. S.; Salonen, E.; Marrink, S. J.; Vattulainen, I. Faraday Discuss. 2010, 144, 411-430. (15) Niemela, P. S.; Ollila, S.; Hyvonen, M. T.; Karttunen, M.; Vattulainen, I. PLoS Comp. Biol. 2007, 3, e34. (16) Mills, T. T.; Tristram-Nagle, S.; Heberle, F. A.; Morales, N. F.; Zhao, J.; Wu, J.; Toombes, G. E. S.; Nagle, J. F.; Feigenson, G. W. Biophys. J. 2008, 95, 682-690. (17) Kahya, N.; Scherfeld, D.; Bacia, K.; Schwille, P. J. Struct. Biol. 2004, 147, 77-89. (18) Scherfeld, D.; Kahya, N.; Schwille, P. Biophys. J. 2003, 85, 3758-3768. 21 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
(19)
Page 22 of 40
Bachar, M.; Brunelle, P.; Tieleman, D. P.; Rauk, A. J. Phys. Chem. B 2004, 108, 7170-
7179. (20) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871-4880. (21) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72, 2002-2013. (22) Holtje, M.; Forster, T.; Brandt, B.; Engels, T.; von Rybinski, W.; Holtje, H. D. Biochim. Biophys. Acta, Biomembr. 2001, 1511, 156-167. (23) Berendsen, H. J. C. P., J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction models for water in relation to protein hydration. In Intermolecular Forces; Pullman, B., Ed.; Reidel, Dordrecht, 1981; pp 331-342. (24) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435–447. (25) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701-1718. (26) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577-8593. (27) Veatch, S. L.; Keller, S. L. Biophys. J. 2003, 85, 3074-3083. (28) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463-1472. (29) Edholm, O.; Nagle, J. F. Biophys. J. 2005, 89, 1827-1832. (30) Rog, T.; Pasenkiewicz-Gierula, M.; Vattulainen, I.; Karttunen, M. Biochim. Biophys. Acta, Biomembr. 2009, 1788, 97-121. (31) Pandit, S. A.; Scott, H. L. Biochim. Biophys. Acta, Biomembr. 2009, 1788, 136-148. (32) Martinez-Seara, H.; Rog, T.; Karttunen, M.; Vattulainen, I.; Reigada, R. Plos One 2010, 5, e11162. (33) Wohlert, J.; Edholm, O. J. Chem. Phys. 2006, 125, 204703. (34) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graph. 1996, 14, 33-38, 27-38. (35) König, S.; Pfeiffer, W.; Bayerl, T.; Richter, D.; Sackmann, E. J. Phys. II 1992, 2, 15891615. (36) Tabony, J.; Perly, B. Biochim. Biophys. Acta 1991, 1063, 67-72. (37) Oradd, G.; Westerman, P. W.; Lindblom, G. Biophys. J. 2005, 89, 315-320. (38) Cahn, J. W.; Hilliard, J. E. J. Chem. Phys. 1958, 28, 258-267. (39) de Joannis, J.; Coppock, P. S.; Yin, F.; Mori, M.; Zamorano, A.; Kindt, J. T. J. Am. Chem. Soc. 2011, 133, 3625-3634. (40) Pandit, S. A.; Chiu, S. W.; Jakobsson, E.; Grama, A.; Scott, H. L. Langmuir 2008, 24, 6858-6865. (41) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089-10092. (42) Petersen, H. G. J. Chem. Phys. 1995, 103, 3668-3679.
22 ACS Paragon Plus Environment
Page 23 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Table 1. Pair wise maximum distances of nearest neighbors for the UA and CG models. lipid/CHOL UA
DPPC/DUPC
CHOL
atoms
C13
C22, C41
O6
C1
C13
C13
0.72
-
0.65
C22, C41
-
0.75
-
O6
0.65
-
0.65
C1
-
0.57
-
0.55
0.7
C13
-
0.8
-
0.7
0.8
atoms
GL1
C2A, C2B
D2A, D2B
ROH
R1 / R2
GL1
0.68
-
-
0.72
-
C2A, C2B
-
0.75
0.75
-
0.61 / 0.64
D2A, D2B
-
0.75
0.75
-
0.65 / 0.63
ROH
0.72
-
-
0.7
-
R1 / R2
-
0.61 / 0.64
0.65 / 0.63
-
0.77
-
DPPC/DUPC
CHOL
CG
DPPC/DUPC
0.57
0.8 -
CHOL The dashes indicate non-neighbor pairs. All the distances are in nm units.
23 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 24 of 40
Figure 1. DUPC structure in the CG and UA models. Each bead in the CG DUPC structure is labeled by the bead type and the bead name. The second and third beads in each chain of the CG DUPC structure are more apolar (type C4, marked with thickened circles) than the other beads of the chains and reflect the double-bond regions. The UA PLPC lipid is show in the middle. The UA DUPC lipid is prepared from the PLPC lipid by making the 16:0 saturated chain be identical to the 18:2 diunsaturated chain of PLPC. Figure 2. The atoms of the lipids and the cholesterols for the UA model and the CG model which are relevant for the nearest neighbor and CHOL tilt axis: (A) UA DPPC/DUPC lipid (B) CG DPPC/DUPC lipid. Atom names in brackets correspond to DUPC lipid. (C) CHOL for the UA model and (D) CHOL for the CG model. Figure 3. Cholesterol models with schematic illustration of the angles for an orientation of a pair of cholesterols. For the UA model the angle is formed by the atoms C1, C2 of one CHOL and the center of mass point defined by the atoms C1, C2, C3 and C9 of another CHOL (A). For the CG model the projection of the center of R1 bead onto the CHOL plane formed by the beads R2, R3 and R4 (shown as a triangle on two sides of which the centers of the R2, R3 and R4 beads reside) is used for the definition of the central point (P1) of the angle (B). The positions of the bead R1, projection point P1 of one CHOL and the center of mass point of the beads R1 and R2 of the other CHOL define the angle. In the schematic (C) two cholesterols (CHOL1 and CHOL2) are shown in grey (view from top). The methyl groups of CHOL (the short grey sides) are positioned asymmetrically with respect to the CHOL plane (the long grey sides). The C1 and C2 points correspond to the atoms of the UA CHOL model and the R1 and P1 shown in braces correspond to the R1 bead and the constructed P1 point of the CG model (B). The arcs show the angles α and β which define the relative orientations of the cholesterols. 24 ACS Paragon Plus Environment
Page 25 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
For simplicity, the angle formed with a neighbor CHOL is shown to be connected to the C2 (P1) point, although the real point is defined according to the center of mass of the atoms (or beads) defined above. Figure 4. Snapshots of one of the bilayer leaflets of the UAs and the CGs at initial, intermediate and final times. On the right side the initial and final snapshots of one of the leaflets of the UAraft-like system are shown. The DPPC, DUPC lipids and CHOL are shown in green, blue and red, respectively. Side lengths of the UA, UAraft-like and CG systems are ~16, ~18 and ~20 nm, respectively. Figure 5. Number of nearest DPPC neighbors of DPPC lipids for five independent CG runs. The number in the legend indicate the number of DPPC neighbors. Each CG run is colored separately and is differentiated by a corresponding superscript number. Figure 6. The MSD of DUPC and the diffusion coefficients of lipids of the CGs and the UAs. The CG data is averaged over five independent runs. The shown MSD for DUPC (A) is the average of MSD of two leaflets where the displacement of only GL1 beads (CG) and C13 atoms (UA) relative to the overall movement of the corresponding leaflet are taken into account. The dashed lines represent the limiting slopes where the ratio lg(MSD) / lg(t) is equal to 1 and help to identify the times where the global diffusions of the DUPC lipids (as well as for the DPPC) start. The diffusion coefficients of DPPC and DUPC lipids (B) are shown starting from 15 and 500 ns for the CGs and the UAs, respectively. The time intervals for diffusion calculation around each data point varied in the range of 50-200 ns and 900-1200 ns for the CGs and the UAs, respectively. Only the linear parts of MSD curves for those time intervals after 15 ns and 500 ns were taken into account for the diffusion. The upper and lower horizontal dotted lines crosses the first points of the diffusion coefficients of DUPC lipids. The vertical dotted line indicates the difference (factor of 42) of diffusion coefficients between the CGs and UAs for the mixed systems. The arrows represents the diffusion coefficients of 3.5 × 10-8 25 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 26 of 40
cm2s-1 and 0.5 × 10-8 cm2s-1 for of UAraft-like system after 900 ns simulation for the DUPC and DPPC lipids, respectively, thus showing the limiting values for the fully separated UAs. Figure 7. Numbers of nearest neighbors and order parameters of DPPC lipids. The CG data is averaged over five independent runs. The fraction of DPPC lipids with a certain number of CHOL (A) and DPPC (B) neighbors, respectively, as a function of time. The average order parameters of DPPC for the UAs (C) and the CGs (D), respectively. The DPPC order parameters are shown as a function of the number of nearest DPPC lipids which are shown in the legends. The vertical dotted lines separate the regions of negligible (left) and considerable (right) impact of local neighborhood. The vertical lines are shown at 500 ns time and originate primarily from the differentiation of DPPC order parameters with respect to the number of nearest DPPC neighbors (C and D). For the sake of clarity the same vertical line is shown in all the subsequent figures up to Figure 11 inclusive. The CG times are scaled by a factor of 40. Figure 8. The order parameters of DUPC lipids for the UAs and the CGs. The CG data is averaged over five independent runs. The DUPC order parameters are shown as a function of two and four nearest DUPC lipids which are labeled accordingly. The dashed line presents the average order parameter of DUPC for the UAs The vertical dotted line separate the regions of negligible (left) and considerable (right) impact of local neighborhood. The CG times are scaled by a factor of 40. Figure 9. The average number of nearest CHOL neighbors of DPPC and DUPC for the UAs and the CGs (A). The averaged ratio of CHOL around DPPC and DUPC facing a lipid with either smooth (α) or rough (β) faces for the UAs (B) and the CGs (C) . The CG data is averaged over five independent runs. The ratio is determined as a fraction of all the cholesterols with a given face over the total number of cholesterols among the nearest neighbors of lipids. The vertical dotted lines in all the
26 ACS Paragon Plus Environment
Page 27 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
subfigures separate the regions of negligible (left) and considerable (right) impact of local neighborhood. The CG times are scaled by a factor of 40. Figure 10. The amount of CHOL molecules which have one or more nearest CHOL neighbors. The CG data is averaged over five independent runs. The vertical dotted line separate the regions of negligible (left) and considerable (right) impact of local neighborhood. The CG times are scaled by a factor of 40. Figure 11. The average tilt angles of CHOL for the UAs and the CGs. The CG data is averaged over five independent runs. The vertical dotted line separate the regions of negligible (left) and considerable (right) impact of local neighborhood. The CG times are scaled by a factor of 40. Figure 12. The distribution of the second nearest CHOL – CHOL neighbors and their relative orientations. The CG data is averaged over five independent runs. The distributions are shown for four consecutive time intervals for both UAs (A) and CGs (B) where an angle step of 30o is used for the horizontal axis. The vertical axes in both subfigures show the probability of having a certain angle between any given CHOL and a CHOL from the second nearest neighborhood. The CG times shown in the legend are scaled by a factor of 40.
27 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 28 of 40
Figure 1
28 ACS Paragon Plus Environment
Page 29 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Figure 2
29 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 30 of 40
Figure 3
30 ACS Paragon Plus Environment
Page 31 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Figure 4
31 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 32 of 40
Figure 5
32 ACS Paragon Plus Environment
Page 33 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Figure 6
33 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 34 of 40
Figure 7
34 ACS Paragon Plus Environment
Page 35 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Figure 8
35 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 36 of 40
Figure 9
36 ACS Paragon Plus Environment
Page 37 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Figure 10
37 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 38 of 40
Figure 11
38 ACS Paragon Plus Environment
Page 39 of 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Figure 12
39 ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 40 of 40
Table of Contents (TOC)
40 ACS Paragon Plus Environment