Microsecond molecular dynamics simulations of lipid mixing.

Molecular dynamics (MD) simulations of membranes are often hindered by the slow lateral diffusion of lipids and the limited time scale of MD. In order to study the dynamics of mixing and characterize the lateral distribution of lipids in converged mixtures, we report microsecond-long all-atom MD simulations performed on the special-purpose machine Anton. Two types of mixed bilayers, POPE:POPG (3:1) and POPC:cholesterol (2:1), as well as a pure POPC bilayer, were each simulated for up to 2 μs. These simulations show that POPE:POPG and POPC:cholesterol are each fully miscible at the simulated conditions, with the final states of the mixed bilayers similar to a random mixture. By simulating three POPE:POPG bilayers at different NaCl concentrations (0, 0.15, and 1 M), we also examined the effect of salt concentration on lipid mixing. While an increase in NaCl concentration is shown to affect the area per lipid, tail order, and lipid lateral diffusion, the final states of mixing remain unaltered, which is explained by the largely uniform increase in Na(+) ions around POPE and POPG. Direct measurement of water permeation reveals that the POPE:POPG bilayer with 1 M NaCl has reduced water permeability compared with those at zero or low salt concentration. Our calculations provide a benchmark to estimate the convergence time scale of all-atom MD simulations of lipid mixing. Additionally, equilibrated structures of POPE:POPG and POPC:cholesterol, which are frequently used to mimic bacterial and mammalian membranes, respectively, can be used as starting points of simulations involving these membranes.

Introduction

Lipid bilayers provide the framework of cellular compartmentalization and are essential for the function of many proteins embedded within them. The structural and dynamic features of bilayers have been investigated extensively using molecular dynamics (MD) simulations,[1−7] and considerable development in lipid force fields has been reported in the past few years.[8−13] Nonetheless, all-atom simulations still face the challenge posed by the limited time scale of MD and the slow lateral diffusion of lipids. Because of such slow diffusion, studying properties of complex membranes can be prohibitively expensive using atomistic MD simulations. This is exemplified by the mixing process of two or more types of lipids,[14−20] which generates mixed bilayers that better resemble cellular membranes than bilayers composed of a single lipid species. For this reason, lipid mixing is often studied using a number of other methods. For example, coarse-grained models[21−25] significantly reduce the computational cost and provide an effective approach in studying processes such as vesicle fusion.[26−28] Monte Carlo (MC) or hybrid MC–MD simulations have also been applied to investigate the mixing of different lipid species.[29,30]

In this study, we report all-atom MD simulations of lipid mixing performed on the special-purpose machine Anton.[31] Two mixtures, POPE:POPG (3:1) and POPC:cholesterol (2:1), are simulated for a total of 12 μs. The goal of the simulations is threefold: First, we aim to investigate the dynamics of lipid mixing and characterize the lateral distribution of lipids in converged mixtures. To this end, two sets of simulations are performed for each bilayer: in the first set, the lipids are mixed randomly; in the second set, the minor lipid species (POPG or cholesterol) are placed in a circle at the center of the bilayer. These two simulations enable us to minimize the effect of initial conditions and investigate the convergence of mixing. Our results indicate that POPE:POPG and POPC:cholesterol are each fully miscible at the simulated conditions, with the final states of the mixed bilayers very similar to a random mixture. Second, by simulating three POPE:POPG bilayers at different NaCl concentrations (0, 0.15, and 1 M), we also investigated the effect of salt concentration on lipid mixing. Although the area per lipid, tail order parameter, and lipid lateral diffusion are clearly affected by the increase in NaCl concentration, the final state of the mixtures, i.e., the ratio of unlike to like neighbors and the cluster size distribution, remains largely unaltered. These results are explained by the essentially uniform increase of Na+ ions around POPE and POPG. Finally, the microsecond-long Anton simulations reported here also enable us to measure water permeation directly from MD trajectories. Different water permeability is observed for bilayers with different lipid compositions and salt concentrations. Below we discuss these results in details and provide the equilibrated structures of all bilayers in the Supporting Information.

Methods

Simulation Protocols

System Preparation

A pure POPC bilayer is generated using the membrane plugin of VMD[32] with 85 molecules in each leaflet. The membrane normal is placed along the z-axis, and a 15 Å water (TIP3P) layer is added to each side of the bilayer. The final system consists of 170 POPC molecules and 6631 water molecules, corresponding to a water:lipid ratio of 39:1. This system also contains 19 Na+ and 19 Cl– ions, which yields a salt concentration of ∼0.15 M. To prepare the mixed POPC:cholesterol (PC:CHL) system, the CHARMM membrane builder[33] is used to generate an initial structure with 70 POPC and 35 cholesterol in each leaflet. The system contains 6682 water molecules, 19 Na+, and 19 Cl– ions (salt concentration of ∼0.15 M). As POPC and cholesterol are mixed randomly in this system, we will refer to it as PC:CHL-0.15r. A second PC:CHL system is then generated by rearranging all cholesterol molecules in a circle at the center of the bilayer. We will refer to this system as PC:CHL-0.15c.

To prepare the POPE:POPG (PE:PG) mixture, a pure POPE bilayer with 85 molecules in each leaflet is first constructed using the VMD membrane plugin. In each leaflet, 21 POPE molecules are selected randomly and mutated into POPG through atom deletion or renaming. Approximately equal numbers of L- and D-conformations were generated for the chiral carbon in POPG. The final system, which has a mixing ratio of 3:1 (POPE:POPG), contains 170 lipids and 6040 water molecules. As each POPG molecule carries a −1 charge, 42 Na+ ions are added to neutralize the system. Since no additional salt is added, we will refer to this system as PE:PG-0r. It should be emphasized that this system is not ion-free: Na+ ions are added to neutralize the charge of the system. We construct two other systems by adding NaCl to PE:PG-0 at a concentration of 0.15 M (17 Na+, 17 Cl–) and 1 M (113 Na+, 113 Cl–), which we will refer to as PE:PG-0.15r and PE:PG-1.0r, respectively. Structures of the PE:PG-0r, PE:PG-0.15r, and PE:PG-1.0r systems are shown in Figure S1. Finally, we also construct a PE:PG mixture with all POPG molecules arranged in a circle at the center of the bilayer. The salt concentration in this system is 0.15 M, and we will refer to it as PE:PG-0.15c. A complete list of bilayer systems studied in this paper is shown in Table 1.

Table 1

Lipid Bilayers Studied in This Work

system	no. lipids	mixing ratio	initial configa	NaCl (M)	t (μs)
PE:PG-0r	170	3:1	random	0	2.0
PE:PG-0.15r	170	3:1	random	0.15	2.1
PE:PG-0.15c	170	3:1	centered	0.15	2.1
PE:PG-1.0r	170	3:1	random	1.0	2.0
PC:CHL-0.15r	210	2:1	random	0.15	2.2
PC:CHL-0.15c	210	2:1	centered	0.15	2.2
POPC-0.15	170			0.15	1.0

Initial configurations of a mixture were either random (r), where lipid components are distributed randomly, or, centered (c), where the minor lipid species are placed in the center of the bilayer.

System Equilibration

As the POPC and PE:PG bilayers are constructed using the VMD membrane plugin, the majority of the lipid tails are in the trans conformation initially. Therefore, after a 5000-step minimization, a 1 ns constant temperature and constant pressure (NPT) simulation is performed to “melt” the tails, during which the phosphorus atoms are constrained in the z-dimension and the x:y ratio is kept constant. Another 5–10 ns simulation is then performed with the constraints on phosphorus atoms released. The PC:CHL bilayers constructed using the CHARMM membrane builder have “melted” lipid tails. Therefore, after minimization, a ∼10 ns NPT simulation is performed to equilibrate each PC:CHL system. The final snapshot of each equilibration is used as the starting structure for subsequent Anton simulations.

All minimization and equilibration are performed with the 2.8 release of NAMD[34] using the CHARMM36 force field.[8,35] A time step of 2 fs is used, with short-range forces calculated every step and long-range electrostatics calculated every two steps. Bonds involving hydrogen atoms were constrained using RATTLE,[36] and water geometries were maintained using SETTLE.[37] The cutoff for short-range nonbonded interactions is set to 12 Å, with a switching distance of 10 Å. Assuming periodic boundary conditions, the particle mesh Ewald (PME) method[38] with a grid density of at least 1/Å3 is employed for computation of long-range electrostatic forces.

Anton Simulations

All the Anton simulations are performed under the semi-isotropic NPT conditions (1 atm and 303.15 K). The Berendesn thermostat and barostat (Ber_NPT) are employed, with default compressibility (κ = 4.5 × 10–5), pressure relaxation time (τ = 2.0), and temperature relaxation time (τ = 1.0). Viparr 1.5.1 was used in the preparation of the simulation systems with the CHARMM36 force field. A time step of 2 fs is used throughout all simulations, with the bonded forces updated every step, and the nonbonded near and far forces updated every 1 and 3 steps, respectively. Trajectories are saved every 120 ps. All the calculations are performed using anton software version 2.6.4. The benchmarks for PE:PG and PC:CHL simulations are ∼5.6 and ∼5.1 μs/day, respectively.

Analysis

Radial Pair Distribution Function g(r)

The 2D radial pair distribution function g(r) is defined aswhere r is the (in-plane) distance between centers of mass of two particles, p(r) is the average number of pairs found at a distance between r and r + dr, A is the area of the system in the xy-plane, and N is the number of unique pairs: N = n·(n – 1) for n identical particles, and N = n1·n2 for n1 type 1 particles and n2 type 2 particles. The 3D g(r) computation routine in VMD[32] is used to perform the calculation, with only the x- and y-components of lipid center-of-mass considered. A scaling factor of 2r/Z is used to convert the 3D result to 2D, where Z stands for the simulation box length along the z-dimension. A resolution of 0.5 Å is adopted. Results for the upper and lower leaflets are calculated separately and then averaged to produce the final g(r).

To quantify the convergence of lipid mixing, we obtain g(r) for every frame in a simulation, i.e., g(r,t), and then determine its autocorrelation function C(r,t)The autocorrelation time τ(r) is then obtained:where Tzero(r) stands for the time that C(r,t) first reaches zero. The statistical inefficiency G(r) = 1 + 2τ(r) is used in the corresponding error analysis. The above calculation follows the protocol used by Chodera et al.[39] to treat noisy simulation data.

Clustering Analysis

Given its advantage in handling clusters with arbitrary shape, the density based spatial clustering of applications with noise (DBSCAN) algorithm[40] is used here to characterize the lateral distribution of lipids. Two parameters are used in our DBSCAN calculation: the parameter ϵ is the maximum distance for a point to be considered “density reachable” from another point, and the parameter MinNrs is the minimum number of neighbors required to form a cluster. A point that is not within any cluster is labeled as noise. Similar to the g(r) calculation, only the x- and y-components of lipid center of mass are considered in the clustering analysis. A sample DBSCAN script used in our analysis is included in the Supporting Information.

Electrostatic Potential Map

The PMEPOT plugin of VMD is used to calculate the electrostatic potential of the bilayers. As described by Aksimentiev and Schulten,[41] the plugin allows the approximation of a point charge by a spherical Gaussian:where the Ewald factor β represents the width of the Gaussian. The electrostatic potential is obtained by solving the Poisson equationThe above calculation is performed for every frame in the last 500 ns of selected simulations and then averaged to yield the 3D electrostatic potential. A grid spacing of 1.0 Å is used, along with a Ewald factor of 0.25. The resulting 3D electrostatic potential is averaged in both x- and y-axes in order to generate a 1D profile along the membrane normal.

Other Analysis

The area per lipid, the deuterium order parameter (SCD), and the lateral mean-square displacement (MSD) of lipids are calculated following protocols described in our previous work.[42] The lateral diffusion coefficient (D) is obtained from MSD data with a time interval between 20 and 200 ns. Water permeation events are calculated using the program VMD.[32] It is worth mentioning that due to the large spacing (0.12 ns) of saved frames in Anton trajectories, a water molecule can move from the vicinity of the upper lipid monolayer to the lower monolayer in just one frame. Therefore, care must be taken to distinguish whether a given water has passed through the bilayer or crossed the periodic boundary.

Results and Discussion

Mixing Dynamics on the Microsecond Time Scale

As described in the Methods section, we chose two different initial configurations for our mixing simulations: in the first configuration, lipids are placed randomly; in the second configuration, the minor lipid species is placed at the center of the bilayer. These two “opposite” configurations allow us to minimize the effect of initial conditions and examine the convergence of mixing. Figure 1 and Figure S2 are snapshots from these simulations, which indicate qualitatively that they have reached a similar final state. Quantitatively, the time evolution of g(r) for the two sets of simulations is shown in Figure 2 and Figure S3. These figures clearly demonstrate that simulations initiated from random configurations have essentially unaltered g(r) throughout the microsecond run, whereas simulations initiated from centered configurations have much larger changes in g(r). Nevertheless, g(r) of all lipid pairs converged in a few hundreds of nanoseconds in the latter simulations. Furthermore, when averaged over the last 500 ns of trajectories, the final g(r) is found to be very similar in simulations initiated from both configurations (Figure S4).

Figure 1

Snapshots of PE:PG-0.15c (a) and PC:CHL-0.15c (b) simulations. The minor lipid species (POPG: blue; CHL: green) were placed in the center of the bilayer at the beginning of both simulations. The simulation box and parts of its eight neighboring periodic images are shown.

Figure 2

Time evolution of the radial pair distribution functions g(r) in PE:PG-0.15c (top) and PC:CHL-0.15c (bottom) simulations. Calculated g(r) is averaged in 100 ns blocks and colored by simulation time, with blue and red indicating the beginning and the end of a simulation, respectively.

In order to further examine the dynamics of mixing, we calculated C(r,t), the autocorrelation function of g(r), and determined the time it takes for C(r,t) to reach zero (T_zero(r)). We then use T_zero(r) to quantify the convergence of mixing. The T_zero(r) obtained here can also be used as a benchmark to estimate the time scale of convergence for future all-atom simulations of lipid mixing. As expected from the time evolution of g(r), convergence is the fastest (T_zero(r) < 60 ns) in simulations initiated from random configurations, while simulations initiated from centered configurations require much longer time to converge (T_zero(r) < 600 ns). As shown in Figure 3, compared with PE:PG simulations, mixing is slightly slower in PC:CHL simulations, reflecting the slower diffusion of lipids in the latter systems—average lipid lateral diffusion coefficient is 3.0 × 10–8 and 2.6 × 10–8 cm2/s in PE:PG-0.15r and PC:CHL-0.15r, respectively.

Figure 3

Convergence of lipid mixing. (a, d) T_zero, the time that the autocorrelation function of g(r) first reaches zero. (b, e) C(t), the autocorrelation function of g(r) at r = 12.25 Å. The blue and red curves are running averages with a window size of 12 ns. (c, f) Running averages of C(t) in PE:PG-0.15c (c) and PC:CHL-0.15c (f) at four representative locations. The same window size is used as in (b, e).

It is worth mentioning that relatively large fluctuations are observed in g(r) data of PG:PG and CHL:CHL pairs (Figure S3). These fluctuations reflect the low sample numbers of the two minor lipid species. For instance, in the PE:PG mixture, the number of PG:PG pairs is only 21·20/(64·63) = 10.4% of that of PE:PE pairs, which significantly increases the statistical noise in g(r) data of the former. Similarly, the number of CHL:CHL pairs is only 24.6% of PC:PC pairs in the PC:CHL mixture. The g(r) data between major and minor lipid species also reflects such statistical noise, albeit to a lesser extent.

Lateral Distribution of Lipids at Long Distances

The rapid convergence of simulations initiated from a random configuration already hints that the final state of the bilayers is similar to a random mixture. In order to quantify the “randomness” of the equilibrated bilayers, we calculated the ratio of unlike neighbors to like neighbors (UL(r)) around a given lipid species a:where r is the radius of a circle centered at a given lipid, N and N are the total number of lipid species a and b, while g(r) and g(r) represent the radial pair distribution function between a and b, and a and a, respectively. N – 1, instead of N, is used in the denominator to be consistent with the exclusion of a lipid from its own neighbor list.

As pointed out previously[43] and clearly shown in Figure S3, a random mixture does not necessarily correspond to a featureless g(r). Indeed, at short distances, g(r) and UL(r) are heavily influenced by molecular size and shape, while their patterns at long distances reflect the miscibility of the components. The average UL(r), a metric similar to the mean excess neighbor used by Coppock and Kindt[43] in their study of DSPC and DMPC mixtures, allows us to quantify the similarity of our bilayers to a random mixture: if the mixing is indeed random, UL(r) at long distances should reflect the mixing ratio of lipid components, that is, UL(r) should be equal to the expected value, N/(N – 1). In the PE:PG mixtures, for instance, the expected ULPE(r) is 21/(64 – 1) = 0.33. A larger-than-expected UL(r) indicates that a lipid is more likely to find an unlike neighbor than a like neighbor, suggesting that there is an effective “attraction” between different lipid species, or, equivalently, a “repulsion” between lipids from the same species.

As shown in Figure 4, ULPE(r) from both PE:PG-0.15r and PE:PG-0.15c simulations is very close to the expected value at r > 10 Å, suggesting that the overall structure of the PE:PG mixture is indeed random. In the PC:CHL simulations, ULPC(r) tends to the expected value at approximately r > 28 Å, indicating that the PC:CHL bilayer also resembles a random mixture. However, one should note that UL(r) is always equal to the expected value when the entire bilayer is considered. Therefore, the size of a simulation box may affect the apparent convergence of UL(r). Both the PE:PG and PC:CHL bilayers studied here are approximately 70 Å by 70 Å in the xy-plane, which is comparable to similar mixtures reported in the literature.[16,17,44] However, the range of convergence in the latter system is somewhat close to the size of the bilayer. Thus, whether the PC:CHL system will deviate from a random mixture in larger simulation boxes remains to be determined.

Figure 4

Average ratio of unlike neighbors to like neighbors (UL(r)) around a given lipid species: (a) UL(r) around POPE in the PE:PG mixtures; (b) UL(r) around POPC in the PC:CHL mixtures. The expected values of UL(r) based on mixing ratios are shown as dashed lines.

Lateral Distribution of Lipids at Short Distances

While both PE:PG and PC:CHL bilayers resemble a random mixture, lateral distributions of lipids at short distances show clear patterns, as already revealed by the g(r) plots. To further examine these patterns, we performed clustering analysis for the minor lipid species (PG and CHL) using the DBSCAN algorithm. Two parameters are used in the calculation: ϵ, which determines whether a lipid is “density-reachable” from another lipid, and MinNrs, which controls the minimum size of a cluster. Here, we set MinNrs to 1, so that the smallest cluster is a dimer, and lipids classified as “noises” can be considered as clusters of size 1 (monomers). With a fixed MinNrs, the parameter ϵ determines the size of a cluster: if ϵ is too small, a large number of noise points will be produced (Figure S5a), whereas if ϵ is too big, small clusters may be lumped together even when they are separated by other lipids (Figure S5b). On the basis of a scan of ϵ (5 to 11 Å), we set it to 9 and 7 Å for the PE:PG and PC:CHL simulations, respectively. The choice of a smaller ϵ for the latter system reflects its higher lipid density.

As shown in Figure 5, on average, 43% out of the 21 POPG molecules exist as monomers in each leaflet of the PE:PG bilayer. Clusters of size 2–4 are frequently observed, which involve approximately 49% of the POPG population. The formation of these clusters is assisted by the headgroup coordination of Na+ ions and water, which counterbalances the electrostatic repulsion between the negatively charged lipids. In the PC:CHL mixtures, approximately 40% out of the 35 cholesterol molecules in a leaflet exist as monomers, while clusters of size 2–5 correspond to 54% of the cholesterol population. Figure S6 shows the histograms of the time it takes for two lipids within the same cluster to become separated. While generally occurring within tens of nanoseconds, this process may take up to ∼200 and ∼300 ns for POPG and CHL, respectively.

Figure 5

Clustering analysis of PE:PG and PC:CHL bilayers. The analysis was performed for the minor lipid species (PG and CHL) using the last 500 ns of each simulation. (a, d) Histograms of cluster size in the PE:PG-0.15r and PC:CHL-0.15r simulations as well as in two randomly constructed data sets. (b, e) Representative clusters identified in the PE:PG-0.15r (b) and PC:CHL-0.15r (e) simulations. Only the center-of-mass of a lipid molecule is shown. The major lipid species (PE and PC) are shown as small, white spheres, and the minor lipid species (PG and CHL) are shown as large, colored spheres. Each color represents a different cluster, with gray representing a cluster of size 1. (c, f) Clusters formed by two POPG (c) and four CHL (f). The POPG and CHL molecules are colored in blue and green, respectively.

Similar to using UL(r) as an evaluation of the “randomness” of our mixtures at long distances, we further analyzed the above cluster size distributions to evaluate the randomness of the mixtures at short distances. To this end, we performed DBSCAN analysis on two sets of random data. In the first set (rand1), 21 or 35 lipids from each monolayer of PE:PG-0.15c or PC:CHL-0.15c were randomly selected and designated as POPG or CHL, respectively. In the second set (rand2), we randomly placed 21 or 35 particles on a plane with the same size as the PE:PG-0.15c or PC:CHL-0.15c bilayers, respectively. No particles were placed within 3.5 Å of each other, in order to mimic the volume exclusion of lipids (see Figure 2 and Figure S3). For both random sets, a total of 50 000 samples were generated—for rand1, this was achieved by generating 1000 random selections for frames taken every 10 ns in the last 500 ns of the above two simulations. DBSCAN analysis was then performed on these two sets using the same settings as described above. Interestingly, as shown in Figure 5, the results of both sets are remarkably similar to the lipid mixtures simulated here. Along with the UL(r) data, these results indicate that mixing of our PE:PG and PC:CHL bilayers is random at both long and short distances. The slightly bigger difference between rand1 and PC:CHL-0.15c (Figure 5b) reflects the larger size difference of POPC and CHL, in comparison to POPE versus POPG. More interestingly, the comparison of both mixtures with the corresponding rand2 sets shows that only limited geometric constraints, such as the size of the bilayer plane and the minimum distance between particles, may be required to produce the distributions observed in our mixtures. It will be of interest to repeat the above comparison for other lipid mixtures and/or using other clustering algorithms in future studies.

The interaction between cholesterol and phospholipids has been the subject of extensive simulation studies in recent years.[23,30,45−51] The mismatch between the small headgroup and large nonpolar body of a cholesterol has been suggested to drive its preferential association with phospholipid molecules.[52] Our simulations show that the hydrophobic body of a cholesterol is indeed covered by neighboring PC molecules, with its hydroxyl group preferentially residing in the plane of glycerol backbones of PC (Figure 5f). Such a conformation is in clear contrast to that of PE:PG mixture (Figure 5c), where the headgroups of both lipids approximately reside on the same plane (Figure S7). The relatively fast “demixing” observed in the PC:CHL-0.15c simulation also reflects the preference of cholesterol to phospholipids. However, as clearly shown in Figure 5, small cholesterol clusters of size 2–5 can still form, as neighboring PC molecules provide sufficient coverage for the hydrophobic bodies of cholesterol.

Effect of Salt Concentration on Lipid Mixing

Salt concentration has been shown to affect the structural and dynamic properties of a bilayer, as well as the activity of antimicrobial peptides, which achieve their function by disrupting the integrity of cellular membranes.[53−56] Our PE:PG simulations under 0, 0.15, and 1.0 M NaCl revealed an area per lipid of 58.3, 57.7, and 56.0 Å, respectively. Such a decrease in area per lipid is accompanied by an increase in tail order parameters (Figure S8), and the lipid lateral diffusion coefficient drops from 3.4 × 10–8 cm2/s in PE:PG-0r to 3.1 × 10–8 cm2/s in PE:PG-0.15r and 2.7 × 10–8 cm2/s in PE:PG-1.0r.

Through calculating the electrostatic potential, we further analyzed the difference between the three PE:PG bilayers under increased NaCl concentrations. As shown in Figure S9, while the total electrostatic potential profiles are similar in all three systems, contributions from individual components differ significantly. In particular, the contribution from water dipoles decreases as the salt concentration increases, while the contribution of both Na+ and Cl– increases dramatically. These results are in line with a previous study on POPC bilayers by Böckmann et al.[57] However, it is worth noting that the total electrostatic potential is not completely flat in bulk water in the PE:PG-0r and PE:PG-1.0r systems (Figure S9a). This result may be explained qualitatively by the ion concentrations in these systems: in PE:PG-0r, the limited number of Na+ ions is insufficient to fully screen the negative charges of POPG, resulting in the “residual” electrostatic potential in bulk water; in PE:PG-1.0r, although sufficient Na+ ions are added, the exceedingly high ionic concentration requires a larger simulation box to level the electrostatic potential.

Despite their difference in structural and dynamic properties, the three systems described above demonstrate similar mixing behaviors: their final states all resemble a random mixture (Figure S10), and no significant difference is observed in the lateral distribution of lipids at short distances, as revealed by the nearly identical histograms of cluster size (Figure S11). This somewhat unexpected result can be explained by ion distributions in the PE:PG mixtures. Through integrating the radial pair distribution function between Na+ and the ester and phosphodiester oxygens of POPE and POPG, we obtain the average number of Na+ ions around the above two types of oxygens. The typical distance between Na+ and an oxygen is ∼3.3 Å. Integration performed up to this distance shows that in the PE:PG-0r system there is an average of 0.04 and 0.05 Na+ ion around an ester oxygen (atoms O22 and O32) from POPE and POPG, respectively. In the PE:PG-0.15r system, these numbers increase slightly to 0.05 and 0.06; when the salt concentration further increases to 1.0 M, these numbers become 0.09 and 0.10, respectively. The above results and the integration curves shown in Figure 6 clearly indicate that as the salt concentration increases, the number of Na+ ions around POPE and POPG increases rather uniformly; i.e., the increased salt concentration does not appear to distinguish the two lipid components. A similar trend is observed for Na+ ion distribution around the phosphodiester oxygens (atoms O13 and O14). As a result, while higher salt concentration produces the structural and dynamic changes described earlier, it makes no significant difference in the lateral distribution of the two lipid components. We should add that at the time the simulations were performed the NBFIX terms in the current CHARMM36 force field[58] were not implemented on Anton yet. To evaluate their effect on the results presented above, we extended in NAMD the PE:PG-0r, PE:PG-0.15r, and PE:PG-1.0r Anton runs by 100 ns. As Na+ binding has an average lifetime less than 0.2 ns,[58] these simulations are considered sufficient to relax Na+ around the lipids. As expected, the adoption of NBFIX resulted in reduced binding of sodium to both POPE and POPG (Figure S12). Nonetheless, similar to the Anton runs, these simulations revealed an essentially uniform increase in the number of ions around POPE and POPG as the salt concentration increases.

Figure 6

(a) Number density of ions in each 1 Å slab along the membrane normal. The number density profiles of ester oxygens and phosphodiester oxygens peak at ±17 and ±21 Å, respectively (not shown). (b, c) Average number of Na+ ion around an ester oxygen (b) and a phosphdiester oxygen (c) in POPE and POPG. The average number of Na+ is obtained as the integral (I(r)) of the radial pair distribution function between Na+ and the oxygen atoms.

Direct Measurement of Water Permeation

In addition to investigating lipid mixing, the microsecond-long Anton simulations presented here also enable us to directly measure water permeation across a membrane. As mentioned in the Methods section, due to the large spacing (0.12 ns) of saved frames in Anton trajectories, a water molecule may cross the periodic boundary and move from the vicinity of the upper lipid monolayer to the lower monolayer in just one frame (Figure S13). In order to avoid miscounting permeation events, we conducted a parameter scan and then set the central region of a bilayer as −9 Å ≤ z ≤ 9 Å. A water molecule must enter this region at least once in order to be counted as a permeation event. We note that while this treatment helps to remove false permeation counts, it could miss a true permeation event when a water crosses the central region in less than 0.12 ns. Nonetheless, while the latter can result in an underestimation of absolute water permeability, it is unlikely to affect the comparison of various bilayers significantly.

As there is no osmotic pressure across the bilayer in our simulations, water passes through the membrane from both directions with equal probabilities. The total number of permeation events, regardless of direction, reflects the intrinsic water permeability of a bilayer. As shown in Figure 7, on average, 1 permeation event is recorded every 43 ns in the PE:PG-0.15r and PE:PG-0.15c simulations. In the pure POPC bilayer, water permeation is recorded at an average rate of 1 per 17 ns. As expected, introducing cholesterol into the POPC bilayer significantly reduces its water permeability: the PC:CHL-0.15r simulation records 1 water permeation event every 114 ns. As the number of POPC molecules is similar in the two systems, the reduced water permeability largely results from the increased order of lipids in the presence of cholesterol.

Figure 7

Direct measurement of water permeation. (a, b) Number of water permeation events in PE:PG (a) and PC:CHL or POPC (b) bilayers. (c) Permeation time through the central region (−9 Å ≤ z ≤ 9 Å) of bilayers. The histogram is obtained from all simulations combined.

Compared with PC:CHL-0.15r, the PC:CHL-0.15c system shows a much higher water permeation rate (1 per 70 ns). This result is explained by the initially centered configuration of cholesterol molecules—with such a configuration, water can pass through the POPC part of the bilayer at a high rate. After the mixing is completed, its water permeability becomes comparable to the PC:CHL-0.15r system (Figure 7). It should also be noted that although it takes tens of nanoseconds for a water permeation event to occur, the permeation process itself is much faster. As shown in Figure 7c, water permeation through the central region (−9 Å < z < 9 Å) of a bilayer typically takes less than 1 ns, regardless of its lipid composition. This number reflects a high free energy barrier of water molecules in this region,[59−61] which prevents them from residing inside the bilayer for an extended period of time.

Notably, salt concentration clearly affects the water permeability of a bilayer—while only small difference is observed between PE:PG-0r and PE:PG-0.15r simulations, the PE:PG-1.0r simulation shows a significant decrease in the number of permeation events (Figure 7). This result can be explained by the structural changes caused by the increased salt concentration; i.e., the higher tail order and smaller area per lipid make the lipids more tightly packed, thereby reducing its water permeability. Additionally, the reduced water permeation may also arise from the increased free energy barrier of crossing the bilayer, since high salt concentration lowers the chemical potential of water in the bulk phase.

Conclusions

In this study, we performed microsecond-long all-atom MD simulations to investigate the dynamics of lipid mixing and characterize the lateral distribution of lipids in converged mixtures. Two mixtures, POPE:POPG (3:1) and POPC:cholesterol (2:1), as well as a pure POPC bilayer were each simulated for up to ∼2 μs. Our results indicate that POPE:POPG and POPC:cholesterol are each fully miscible at the simulated conditions, with the final states of the mixed bilayers similar to a random mixture. Through clustering analysis and g(r) calculation, we further examined the lateral distributions of lipids at both short and long distances. It is shown that small clusters of 2–4 POPG and 2–5 cholesterol form frequently throughout the simulations, which are also observed in two randomly constructed data sets. By simulating three POPE:POPG bilayers at different NaCl concentrations (0, 0.15, and 1 M), we also investigated the effect of salt concentration on lipid mixing. While the area per lipid, tail order, and lipid lateral diffusion are clearly affected by the increase in NaCl concentration, the lateral distribution of lipids remains unaltered. These results are explained by the largely uniform increase in the number of Na+ ions around POPE and POPG. Since POPE:POPG and POPC:cholesterol mixture are frequently used to represent bacterial and mammalian membranes, respectively, equilibrated structures from the simulations presented here can be used as starting points for future simulations involving these bilayers.