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 , and considerable development in lipid force fields has been reported 
 in the past few years .
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 , 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 significantly reduce the computational cost and provide an effective 
 approach in studying processes such as vesicle fusion .
md simulations have also been 
 applied to investigate the mixing of different lipid species . in this study , we report all - atom md simulations of lipid 
 mixing performed on the special - purpose machine anton
. 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 .
a pure popc bilayer is generated 
 using the membrane plugin of vmd 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 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 ,
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 .
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 . 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 510 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 using the 
 charmm36 force field .
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 , and water geometries were 
 maintained using settle .
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 with 
 a grid density of at least 1/ is employed for computation 
 of long - range electrostatic forces .
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 ) , 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 .
the 2d radial pair distribution function g(r ) is defined as1where 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 = n1n2 for n1 type 1 particles and n2 type 2 particles .
the 3d g(r ) 
 computation routine in vmd 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 
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)2the autocorrelation 
 time (r ) is then obtained:3where 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 . to treat noisy simulation data .
given its advantage in handling clusters with arbitrary 
 shape , the density based spatial clustering of applications with noise 
 ( dbscan ) algorithm 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 .
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
the pmepot plugin of vmd is used to calculate 
 the electrostatic potential of the bilayers . as described by aksimentiev 
 and
schulten , the plugin allows the approximation 
 of a point charge by a spherical gaussian:4where the ewald factor  represents the width of the gaussian . 

the electrostatic potential is obtained by solving the poisson equation5the 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 .
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 .
the lateral diffusion coefficient 
 ( d ) is obtained from msd data with a time interval 
 between 20 and 200 ns .
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 .
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 ) .
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 .
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 ) ) .
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 )
< 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 and 2.6  10 cm / s in pe : pg-0.15r 
 and pc : chl-0.15r , respectively .
( 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 .
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 
 2120/(6463 ) = 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 .
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:6where r is the radius of a circle centered at a 
 given lipid , na and nb are the total number of lipid 
 species a and b , while gab(r ) and gaa(r ) represent the 
 radial pair distribution function between a and b , and a and a , respectively .
na  1 , instead of na , is used in the denominator 
 to be consistent with the exclusion of a lipid from its own neighbor 
 list . as pointed out previously 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 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 , ula(r ) at long 
 distances should reflect the mixing ratio of lipid components , that 
 is , ula(r ) should be equal to the expected value , nb/(na  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
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 .
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 .
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 . 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 
 24 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 25 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 .
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 .
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 ) .
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 .
the mismatch between the small headgroup and large nonpolar body 
 of a cholesterol has been suggested to drive its preferential association 
 with phospholipid molecules .
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 ) .
observed in the pc : chl-0.15c simulation also reflects the preference 
 of cholesterol to phospholipids .
as clearly shown in figure 5 , small cholesterol clusters of size 25 
 can still form , as neighboring pc molecules provide sufficient coverage 
 for the hydrophobic bodies of cholesterol .
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 .
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 cm / s in pe : pg-0r to 
 3.1  10 cm / s in pe : pg-0.15r 
 and 2.7  10 cm / 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 bckmann 
 et al .
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 .
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 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 , 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 .
( 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 .
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 .
( 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 .
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 , 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 .
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 24 popg and 25 
 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 )
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 .
 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
Methods
Results and Discussion
Conclusions