phospholipids are amphiphilic molecules 
 that spontaneously form 
 bilayers in aqueous solution to minimize the interaction between their 
 hydrophobic tails and to maximize the interaction between their hydrophilic 
 heads . in recent years
, scientists have been working to exploit the 
 tunable properties of mixed lipid systems for a range of applications 
 including imaging , biosensors , and membrane trafficking .
in particular , mixed 
 lipid systems are being used as novel liposomal delivery carriers 
 for drugs , vaccines , genes , and antimicrobial agents . while scientists have 
 successfully demonstrated that liposomal formulations improve drug 
 efficacy , research efforts are still dedicated to 
 optimizing their structure and function .
for example , karve and co - workers 
 showed that liposomes composed of mixed lipid membranes could become 
 permeable to encapsulated drug molecules in response to a decrease 
 in ph .
mixed lipid systems are also being 
 used to create temperature - sensitive liposomes that can rapidly and 
 efficiently release drug molecules in response to heat . in order to fully realize the potential of these new technologies , 
 a better understanding of the structure and interactions between lipids 
 on a molecular level under various physiological conditions
this would enable researchers to design and develop liposomes with 
 optimal properties to release drug molecules in response to specific 
 stimuli .
an approach commonly used to gain insight into the behavior 
 of such systems on a molecular level is computer simulation . in a previous paper , we demonstrated how computer simulations of 
 lipid systems can complement experimental work by providing detailed 
 molecular - level information regarding their dynamics and self - assembly . in that work
, we introduced lime ( lipid intermediate resolution model ) , a new implicit solvent 
 lipid force field for use with discontinuous molecular dynamics ( dmd ) , 
 that enabled the rapid simulation of 1,2-dipalmitoyl - sn - glycero-3-phosphocholine ( dppc ) lipids in water .
the structural 
 properties of the dppc bilayer formed during our lime / dmd simulations , 
 including the area per lipid , bilayer thickness , bond order , and mass 
 density profiles , were in agreement with experimental observations 
 or atomistic simulation results . in this article , we take a 
 multiscale modeling approach to expand 
 lime to include parameters for ( 1 ) 1,2-dihenarachidoyl - sn - glycero-3-phosphocholine ( 21pc ) and ( 2 ) 1,2-distearoyl - sn - glycero-3-phospho - l - serine ( dsps ) at 310 k. we demonstrate 
 how discontinuous molecular dynamics ( dmd ) simulations with this expanded 
 lime force field can be used to model the behavior of lipid systems 
 containing mixtures of different molar ratios of dsps and 21pc over 
 long time scales at neutral ph and at a ph of 5.5 . throughout this 
 article
, we define a neutral ph as 7.07.4 and low ph as < 5.5 . 

we show that lime / dmd simulations can be used to study the phase separation 
 of mixed lipid systems into heterogeneous domains on a molecular level . 

our results are the first reported simulations of the ph - dependent 
 phase separation of mixed lipid bilayers . a number of different 
 approaches to modeling mixed lipid systems 
 with simulations have been described in the literature .
various levels 
 of detail are used in these models , which can be divided roughly into 
 two categories : high- and low - resolution models .
high - resolution , 
 or atomistic , models are based on a realistic representation of membrane 
 geometry and energetics and typically account for the motion of every 
 atom , including every solvent atom .
jiang and co - workers performed 
 gromacs atomistic simulations on bilayer 
 ribbons containing dimyristoylphosphatidylcholine ( dmpc ) , dihexanoylphosphatidylcholine 
 ( dhpc ) , and didecanoyl pc ( ddpc ) to investigate how each species partitions 
 between the flat and curved microenvironments of the bilayer and the 
 role dhpc plays in stabilizing the bilayer edge . in another study , hall and co - workers performed atomistic 
 molecular dynamics simulations to investigate the role of glycolipids 
 in lipid rafts ; systems of galactosylceramide ( galcer ) , cholesterol 
 ( chol ) , palmitoyloleoylphosphatidylcholine ( popc ) , and palmitoyl - sphingomyelin 
 ( psm ) were simulated using gromacs .
addition 
 of galcer increased the thickness of the raft membranes and induced 
 changes in the lateral diffusion of raft lipids , but it did not influence 
 the average area per lipid or the lipid conformation order . in another atomistic simulation , pandit and 
 co - workers studied mixed bilayers containing a 5:1 ratio of dppc and 
 dsps in nacl electrolyte solutions to learn how dppc and dsps interact 
 to form complexes .
gromacs was used to 
 perform all simulations based on force field parameters from the work 
 of berger .
these simulations showed that 
 nho and cho hydrogen bonding between 
 lipids serves as the basis of interlipid complexation and that dppc 
 alone is less likely to form interlipid complexes than in the presence 
 of bound ions or dsps .
low - resolution 
 models , which are based on a simplified representation 
 of molecular geometry and energetics , have also been used to study 
 the behavior of mixed lipid systems .
for example , illya and co - workers 
 developed a low - resolution model designed for use with dissipative 
 particle dynamics simulations to study the domain formation and the 
 material properties of a two - component amphiphilic bilayer membrane . in this solvent - explicit model , the amphiphile 
 geometry was chosen to represent the general shape of a typical phospholipid . 

the amphiphile h3(t6)2 was represented by 3 head spheres and 2 tails , 
 each composed of 6 spheres , and the amphiphile h3(t8)2 was represented 
 by 3 head spheres and 2 tails , each composed of 8 spheres .
the interaction parameters for the amphiphiles 
 were chosen to mimic the properties of lipid membranes , not calculated 
 from atomistic simulations of amphiphiles .
the authors calculated the area stretch modulus and bending rigidity 
 of the mixed amphiphilic bilayers and vesicles but could not compare 
 their results directly to any experimental observations on a specific 
 lipid mixture since the amphiphiles did not represent any specific 
 lipids .
another way to model mixed lipid systems in simulations 
 is to use 
 intermediate - resolution ( coarse - grained ) models , which represent the 
 geometry and energetics of lipids at a level of detail that is in 
 between that of high- and low - resolution models .
this allows for the 
 simulation of specific lipids while requiring shorter simulation times 
 than that for atomistic models .
this reduces the total 
 number of sites whose trajectories must be calculated and helps to 
 increase the speed of the simulation .
one example of an intermediate - resolution 
 explicit - solvent model is that of faller and marrink . in this model ,
distearoylphosphatidylcholine ( dspc ) and 
 dilauroylphosphatidylcholine ( dlpc ) are represented by 14 and 10 coarse - grained 
 sites , respectively .
lipid mixtures composed of dlpc and dspc were 
 simulated to study the phase separation that occurs between these 
 two lipids at phase transition temperatures such that the longer lipid 
 enters the gel phase and the shorter lipid remains in the fluid phase .
simulations were performed using gromacs with 
 the interaction parameters previously described by marrink and coworkers .
all simulations were started from a bilayer 
 composed of randomly assigned dspc and dlpc lipid .
results showed 
 that the higher the ratio of the number of longer lipids ( dspc ) to 
 the number of shorter lipids ( dlpc ) , the higher the gel 
fluid 
 transition temperature . in another intermediate - resolution explicit - solvent 
 model , risselada and
marrink performed simulations of the spontaneous 
 separation of a randomized mixture composed of dipalmitoyl - phosphatidylcholine 
 ( dic16-pc ) , dilinoleyl - phosphatidylcholine ( dic18:2-pc ) , and cholesterol into liquid - ordered and -disordered phases .
simulations of a planar membrane composed of 
 approximately 2000 lipid molecules and a liposome composed of approximately 
 3000 lipid molecules were performed . in this model , on average , 4 
 heavy atoms were represented by a single coarse - grained site , and 
 approximately 3 heavy atoms were used to represent the rings in cholesterol .
the simulations were performed with the gromacs 
 simulation package and the martini cg force field , version 2.0 .
intermediate - resolution implicit - solvent 
 models of mixed lipid 
 systems have also been described in the literature .
for example , lu 
 and voth used a multiscale coarse - graining approach to develop a solvent - free 
 coarse - grained model for a mixed bilayer composed of dioleoylphosphatidylcholine 
 ( dopc ) and dioleoylphosphatidyethanolamine ( dope ) . in this model , dopc and dope are each represented by 15 
 coarse - grained sites .
the results of the multiscale coarse - grained 
 model , including radial distribution functions , bilayer thickness , 
 and order parameters , show good agreement with atomistic results .
our motivation for studying 21pc / dsps 
 lipid systems is experimental 
 work by sofou and co - workers aimed at developing liposome - based drug 
 delivery devices for the treatment of cancer .
they examined the phase separation behavior of ph - triggered liposomes 
 composed of cholesterol and two different types of phospholipids : 
 21pc phospholipids , which have a neutral headgroup at both neutral 
 ( ph 7.47.0 ) and endosomal ph ( ph 5.55.0 ) , and dsps 
 phospholipids , which have a negatively charged headgroup at neutral 
 ph and a neutral headgroup at endosomal ph .
they showed that at neutral ph the vesicles were much less permeable 
 to encapsulated drug molecules than at endosomal ph . at endosomal 

ph , the cholesterol and 21pc and dsps lipids allowed much more of 
 the encapsulated drug molecules to leak out of the liposome than at 
 neutral ph .
this observation forms the basis of sofou and co - workers 
 proposed liposomal - based drug delivery scheme in which liposomes composed 
 of cholesterol , 21pc , and dsps are used to transport cancer drugs 
 throughout the body at neutral ph .
when the liposomes reach the endosome 
 in cancer cells , the decrease in ph triggers the formation of heterogeneous 
 domains and the encapsulated drug molecules leak out of them .
sofou s ph - tunable liposomes have additional features including 
 antibodies to target cancer cells . in this work
, we present 
 the results of simulations of 21pc / dsps 
 bilayers at neutral and low ph .
our goal was to show that the expanded 
 version of our lime / dmd force field could accurately model the phase 
 separation behavior of these bilayers at both neutral and low ph . 

we are currently running simulations of liposomes composed of 21pc / dsps , 
 that encapsulate drug molecules at neutral and low ph .
we hope to 
 soon report the results of these simulations , which will be the ultimate 
 test of the accuracy of our lime / dmd model for our mixed lipid system . in this article
, we describe the expansion of lime , the implicit - solvent 
 coarse - grained model previously used to simulate the spontaneous formation 
 of dppc bilayers at 325 k , to 21pc and dsps lipids at 310 k. we analyze the extent of phase separation that 
 occurs during simulations at both neutral and low ph in 21pc / dsps 
 bilayers containing 150 lipids and molar ratios of 1:1 or 9:1 .
the 
 lime geometrical and energetic parameters for 21pc and dsps at low 
 ph and 310 k are calculated using a multiscale modeling approach based 
 on evaluating radial distribution functions obtained from an atomistic 
 simulation of a fully hydrated bilayer composed of 64 dppc lipids 
 and 64 protonated dsps lipids performed with the gromacs simulation 
 package and the gromos96 53a6 force field .
the same method is used to extract parameters 
 at neutral ph ; the only difference is that the dsps lipids are unprotonated .
highlights of our results include the following : lime / dmd simulations 
 of equimolar bilayers composed of 21pc / dsps phase separate into heterogeneous 
 domains faster at a low ph value than at a neutral ph value regardless 
 of whether surface area restrictions are imposed .
in addition , without 
 any surface area restrictions , bilayers composed of a 9:1 molar ratio 
 of 21pc / dpps also phase separate faster at low ph than at neutral 
 ph .
this observation is supported by experimental data from the sofou 
 lab . according to sofou and co - workers ,
drug molecules encapsulated 
 within liposomes composed of 21pc and dsps escape from liposomes faster 
 at low ph than at neutral ph .
this is thought to 
 be a consequence of having faster heterogeneous domain formation at 
 low ph than at neutral ph .
in lime , 9 different coarse - grained types ( i  ix ) are used 
 to represent the groups of atoms that make up the 21pc and dsps molecules . 

21pc is composed of a polar headgroup that includes a choline , phosphate , 
 and two ester linkages and two nonpolar hydrophobic acyl tails .
dsps 
 is composed of a headgroup that includes carboxyl , amine , and phosphate 
 groups and two ester linkages and two nonpolar hydrophobic acyl tails . 
 at neutral ph ,
this negative charge causes electrostatic repulsion between the dsps 
 lipids , which helps to prevent their phase separation . at a ph of approximately 5.5 ,
when the dsps molecules lose their negative 
 charge , the electrostatic repulsion between these lipids stops , allowing 
 hydrogen bonds to form between their protonated amino groups and deprotonated 
 phosphate groups on the head groups . therefore , 
 at low ph , dsps lipids have a very strong attraction for each other . 

figure 1 shows the structure of a ( a ) 21pc 
 lipid and ( b ) dsps lipid at neutral ph . at low ph
structure of ( a ) 21pc and ( b ) dsps at neutral ph . at low 
 ph
figure 2 illustrates the 
 coarse - graining 
 of ( a ) a dsps molecule from 57 united atoms to 17 coarse - grained sites 
 and ( b ) a 21pc molecule from 60 united atoms to 18 coarse - grained 
 sites . in lime ,
each coarse - grained site is assigned a different coarse - grained 
 type , which is represented by a unique color .
this figure and all 
 other figures depicting lipid molecules throughout the article were 
 generated with visual molecular dynamics ( vmd ) .
type ia and ib ( black ) are used to represent the unprotonated 
 and protonated , respectively , carboxyl headgroup ( site 1 ) on dsps . 

the 
 phosphate groups on dsps and 21pc are represented by type iii ( yellow ) 
 and type ix ( yellow ) ; we represent these phosphate groups with different 
 types because the energy parameters calculated for these coarse - grained 
 sites were significantly different from each other and we therefore 
 felt it was important to treat them as unique coarse - grained types . 
 ester coarse - grained sites are assigned types iv ( red ) and v ( orange ) , 
 respectively .
two coarse - grained types for the ester are needed because 
 coarse - grained iii has one more carbon and hydrogen atom than coarse - grained 
 type iv .
coarse - grained type vi ( cyan ) is used to represent the tail 
 sites ( each contains 3 carbons and 6 hydrogens ) of dsps and 21pc . 

coarse - grained type vii ( gray ) is assigned to the terminal tail groups 
 of 21pc and dsps , representing 2 carbon and 5 hydrogen atoms .
finally , 
 coarse - grained type viii ( purple ) represents the choline entity on 
 21pc .
table 1 lists the atoms included in each 
 coarse - grained site , the type assigned to each coarse - grained site , 
 and the mass of each coarse - grained site .
united atom and coarse - grained 
 representations of ( a ) dsps and 
 ( b ) 21pc .
the color scheme is black ( negatively charged carboxyl group , 
 type ia ; protonated carboxyl group , type ib for dsps site 1 ) ; magenta 
 ( amine group , type ii for dsps site 2 ) ; yellow ( phosphate group , type 
 iii for dsps site 3 ) ; red ( ester group , type iv for dsps site 4 and 
 21pc site 3 ) ; orange ( ester group , type v for dsps site 11 and 21pc 
 site 11 ) ; cyan ( alkyl tail groups , type vi for dsps sites 59 
 and 1216 and for 21pc sites
49 and 1217 ) ; 
 gray ( terminal tail groups , type vii for dsps sites 10 and 17 and 
 for 21pc sites 10 and 18 ) ; purple ( choline entity , type viii for 21pc 
 site 1 ) ; and yellow ( phosphate group , type ix for 21pc site 2 ) .
the simulation method 
 used in this article is the discontinuous 
 molecular dynamics ( dmd ) algorithm , which is a very fast alternative 
 to traditional molecular dynamics simulation that is applicable to 
 systems of molecules interacting via discontinuous potentials , e.g. , 
 hard - sphere and square - well potentials .
for this reason , 
 all of the inter- and intramolecular interactions in our lipid model 
 are represented by a combination of hard - sphere and square - well potentials , 
 as opposed to the lennard 
a hard sphere 
 is an impenetrable solid sphere ; a square - well is a hard sphere surrounded 
 by an attractive well .
expressions for the hard - sphere ( hs ) and square - well 
 ( sw ) potentials between spheres i and j are , respectively12where r is the distance between 
 spheres , ij is the hard sphere 
 diameters , ij is the well 
 diameter , and ij is the well depth . 

these parameters define the strength of the interaction between the 
 coarse - grained spheres . unlike continuous potentials , such as the 
 lennard  jones potential , discontinuous potentials exert forces 
 only when particles collide .
this makes the simulation an efficient 
 event - scheduling algorithm , which allows sampling of longer time scales 
 and larger systems than that with traditional molecular dynamics . a dmd simulation proceeds in the following way .
molecules are placed 
 in an initial configuration consistent with excluded volume and angular 
 constraints .
particle trajectories 
 are followed by calculating the time between each collision and advancing 
 the simulation to the next event .
types of events include a collision 
 between two hard spheres , a bond event when the distance between two 
 bonded spheres reaches a minimum or maximum limit , and square - well 
 events when two spheres enter ( capture ) , unsuccessfully attempt to 
 escape ( bounce ) , or successfully leave ( dissociation ) a square well . in all lime / dmd simulations ,
the simulation temperature is 
 expressed 
 in terms of the reduced temperature t*3where kb is boltzmann s 
 constant , t is the temperature , and  * is the 
 reference interaction strength . for this 
 work ,
this reduced temperature 
 was chosen because it gave the most reasonable results for the bilayer 
 assembly ; at lower temperatures , the particles get kinetically trapped , 
 and at higher temperatures they have too much kinetic energy .
the coarse - grained parameters in lime were obtained using a multiscale 
 modeling technique .
multiscale modeling is a method in which several 
 atoms are grouped into a single coarse - grained site and the parameters 
 used to represent the behavior of this coarse - grained site are extracted 
 from atomistic simulations .
data used to calculate the lime parameters 
 was obtained by running united - atom explicit - solvent simulations at t = 310 k using the gromacs simulation package , version 4.5.4 , along with the gromos96 53a6 force field .
since the sofou lab conducted experiments at 310 k 
 to study how liposomes behave in the body , we parametrized our model 
 from gromacs simulations at 310 k. our gromacs simulations consisted 
 of bilayers composed of dppc and dsps , whereas our lime / dmd simulations 
 were of 21pc / dsps bilayers .
lime parameters for the 21pc head groups 
 were set equal to those of the dppc head groups since 21pc and dppc 
 have identical head groups .
similarly , the 21pc tail group sites were 
 parametrized from the alkyl tail groups of dppc and dsps . to 
 obtain the partial charges for dppc and dsps ( at neutral and 
 low ph ) for our gromacs simulations , we performed quantum mechanical 
 ( qm ) calculations using the gaussian09 software package . to facilitate the qm calculations , each lipid 
 molecule
( dppc , dsps neutral ph , and dsps low ph ) was split into two 
 parts .
a methyl group was used to cap each of the three head groups , 
 and a nh3po4 group was used to cap the common 
 tail group .
in addition , we calculated partial charges only for the 
 first 3 alkyl groups in each of the tails .
the two - stage restrained 
 electrostatic potential ( resp ) charge - fitting protocol was performed 
 to calculate the partial atomic charges .
the geometry 
 of each molecule was optimized at the hf/6 - 31 g * level , and the electrostatic 
 potential for each optimized molecule was calculated using the hf/6 - 31 g * 
 basis set .
the lime low - ph intermolecular ij , ij , 
 and ij between coarse - grained types 
 i  vi , viib , 
 and viii - x were obtained by analyzing results from a low - ph gromacs 
 simulation of a mixed bilayer composed of 64 dppc lipids and 64 dsps 
 lipids in explicit solvent .
the neutral - ph intermolecular values for 
 ij , ij , and ij between coarse - grained 
 types i  vi , viia , and viii - x were calculated from a neutral - ph 
 gromacs simulation composed of 64 dppc lipids and 64 nonprotonated 
 dsps lipids .
each gromacs simulation was started from an initially 
 homogeneous bilayer in a box with dimensions of 64  64  
 90  and was run at a temperature of 310 k and a 
 pressure of 1.0 bar for 10 ns . throughout both the neutral- and low - ph 
 simulations
the net charge of all gromacs 
 simulation systems was neutralized with the addition of nacl .
the 
 berendsen thermostat was used to keep 
 the temperature constant with a time constant of 0.1 ps . throughout 
 each gromacs simulation , the coordinates of each atom were written 
 to an output trajectory file every 1 ps .
these output files were used 
 to calculate the centers of mass for each of the coarse - grained sites .
the hard - sphere diameter ( hs ) for each pair of 
 interaction sites was determined by locating the smallest non - zero 
 separation between the two sites . the hs value for 
 a pair of coarse - grained types
was calculated by averaging all of 
 the hs values for pairs of interaction sites with 
 the desired coarse - grained types .
the depth of the square - well intermolecular 
 interactions in lime was calculated from the radial distribution functions 
 collected during the atomistic gromacs simulations using a one - step 
 boltzmann inversion approach .
( the lime intramolecular ij values were set to 0 , making all intramolecular 
 interactions hard - sphere collisions . )
the one - step boltzmann inversion 
 procedure involves the following steps : ( 1 ) the average radial distribution 
 function , g(r ) , between two nonbonded 
 coarse - grained sites separated by distance r is determined , 
 ( 2 ) the potential of mean force is calculated using5and ( 3 ) the minimum value of the potential 
 of mean force between the coarse - grained sites ,  , is chosen 
 to be the depth of the square - well potential .
mathematically ,  , 
 is expressed as6where g(r)max is the 
 maximum value of g(r ) in the radial 
 distribution function and t is the temperature of 
 the system .
if the  value between two 
 coarse - grained sites is greater than 0.01 ev , then the sites 
 are assumed to have a hard - sphere interaction (  = 0.0 ev ) . 

the  values for pairs of coarse - grained types were determined 
 by averaging all of the  values for pairs of interaction sites 
 with the same coarse - grained type .
we chose to use a different 
 procedure for calculating the square - well 
 diameter (  ) for each pair of interaction sites than 
 the approach described in our previous work . in this new procedure
, 
 the second virial coefficient associated with the potential of mean 
 force between coarse - grained sites ( which is calculated from atomistic 
 simulations ) is set equal to the second virial coefficient for the 
 site  site square - well potential .
the second virial coefficient , b2(t ) , for a system of molecules 
 interacting via a potential u(r ) 
 is defined to be7where k is the boltzmann 
 constant and t is the temperature .
if we take u(r ) to be the potential of mean force 
 calculated using eq 5 , then this becomes8if , instead , we take u(r ) to be a square - well potential as given in eq 2 , then we obtain the second virial coefficient for 
 a square - well interaction9where bo = 2/3 and  = 1/(kt ) .
we determined the value 
 of  for 
 each pair of coarse - grained sites by the following steps .
first , we 
 numerically integrated the right - hand side of eq 8 from 0 to 10  .
we chose to integrate only from 0 to 10  
 because the cutoff radius for our lennard  jones interactions 
 in our gromacs simulations was 10  .
the g(r ) in this step was the same as the g(r ) used in determining the hs and  
 values for that pair of coarse - grained types .
the square - well b2(t ) ( eq 9 ) was set equal to the value calculated from eq 8 .
the hs and  values previously calculated 
 for that pair of coarse - grained sites were used to solve for . 
 once  was calculated , it was multiplied by hs to determine the  for a pair of coarse - grained sites . 

the  value for a pair of coarse - grained types was calculated 
 by averaging all of the  values for coarse - grained sites 
 with the given coarse - grained types .
the minimum and maximum 
 bond lengths between two bonded coarse - grained 
 sites in lime were also determined from the gromacs simulation data 
 using the same method as was used in our previous paper . in lime , 

the stiffness of each lipid molecule is maintained by imposing pseudobonds , 
 which limit the fluctuation of coarse - grained sites to the angles 
 and torsional angles observed during the gromacs simulations .
bond 
 angles are maintained by imposing pseudobonds between all next - nearest 
 neighboring sites .
torsional angles are maintained with pseudobonds 
 between all next  next - nearest neighboring sites .
the minimum 
 and maximum pseudobond lengths were determined from the radial distribution 
 functions for the coarse - grained sites in the same manner as that 
 in our previous paper .
the bond and pseudobond values calculated from 
 the neutral - ph gromacs simulations were very similar to those calculated 
 from the low - ph gromacs simulations .
therefore , for dsps and 21pc 
 lipids , we choose to approximate the min and max values for interactions between coarse - grained types at 
 low ph with those calculated at neutral ph .
emblem , a c++ program 
 developed in the hall research lab , was used 
 to run all lime / dmd simulations .
emblem uses the discontinuous molecular 
 dynamics ( dmd ) algorithm to simulate the behavior of any type of molecule 
 or mixture of molecules .
work is in progress to 
 parallelize portions of emblem ; however , all of the simulations run 
 for this work were done in serial .
lime / dmd 
 simulations were run to study the behavior of bilayers 
 composed of 21pc / dsps lipids .
simulations were run on six systems , 
 referred to here as systems 16 , which differed in the ratio 
 of 21pc to dsps lipids , the ph , and simulation box length .
the six 
 systems are summarized in table 2 . for each 
 system , four replicate simulations , each composed of a total of 150 
 lipids , were run at t * = 0.20 .
each simulation was 
 started from a preformed bilayer in which the lipids of each species 
 were randomly placed into an area of 60  60  . to build each bilayer
, the following steps were performed : ( 1 ) 
 the phosphate group of a lipid molecule was placed at a random location 
 within the x  y plane at a 
 fixed distance from the bilayer normal ( for the top leaflet of the 
 bilayer ) , ( 2 ) a check was run to ensure that the lipid was not overlapping 
 with any other lipid , ( 3 ) steps 1 and 2 were repeated until no overlap 
 was detected , and ( 4 ) steps 13 were repeated for the bottom 
 leaflet of the bilayer . in simulations on systems 14 , the 
 bilayers , which spanned an area of 60  60  ,
were placed in the center of a box with dimensions of 120  
 120  120  .
we refer to the bilayers 
 in these systems ( 14 ) as free surface bilayers because the 
 surface area of these bilayers is free to expand or shrink by any 
 amount .
although interacting with its periodic image would not force 
 a bilayer to occupy a certain area , it might promote the bilayer to 
 span the simulation cell so that the lipids could maximize the number 
 of interactions that they have with each other .
in simulations 5 and 
 6 , the length of the simulation cell in the plane of the bilayer was 
 60  60  , which allowed the bilayers to interact 
 with their periodic images , thereby restricting the possible surface 
 area .
we chose to include surface area restrictions to ensure that 
 our simulation environment was a reasonably faithful mimic of the 
 experimental conditions in the sofou lab . in experiments ,
the osmotic 
 pressure inside liposomes composed of 21pc / dsps in a solution composed 
 of salt ions prevents the liposomes from shrinking to a large extent . 

these lipids would behave similarly to those in our simulations with 
 surface area restrictions , i.e. , with periodic boundary conditions .
simulations on systems 1 
 and 2 were run to investigate the behavior 
 of bilayers composed of equimolar ratios of 21pc : dsps , with free surfaces , 
 at both neutral and low ph .
figure 3 provides 
 snapshots ( aerial images ) of the bilayer configurations at 1 million 
 and 1 billion collisions for a single replicate in simulations on 
 systems 16 .
figure 3a , c shows that 
 after 1 million collisions little to no phase separation can be detected 
 between 21pc and dsps lipids in systems 1 and 2 .
figure 3b , d shows that after 1 billion collisions the 21pc and dsps 
 lipids in system 2 have separated to a greater extent than that in 
 system 1 .
simulations on systems 3 and 4 ( also with free surfaces ) 
 were run to investigate the heterogeneous domain formation of bilayers 
 composed of a 9:1 molar ratio of 21pc : dsps .
figure 3e , g shows that no detectable phase separation is achieved 
 in systems 3 and 4 after 1 million collisions .
however , after 1 billion 
 collisions , the snapshots provided in figure 3f , h show that system 4 has separated to a greater extent than system 
 3 .
snapshots ( aerial images ) of the bilayer formed at 1 million collisions 
 in simulations of systems 1 ( a ) , 2 ( c ) , 3 ( e ) , 4 ( g ) , 5 ( i ) , and 6 
 ( k ) and at 1 billion collisions in simulations of systems 1 ( b ) , 2 
 ( d ) , 3 ( f ) , 4 ( h ) , 5 ( j ) , and 6 ( l ) . the color scheme is lime ( 21pc 
 lipids ) and purple ( dsps lipids at neutral ph and low ph ) . figure 4 displays the normalized , 
 time - averaged 
 percentage of dsps lipid molecules within a 5  radius of each 
 dsps headgroup ( coarse - grained site 1 on each dsps lipid molecule ) . 

the higher the fraction of dsps lipid molecules within a 5  
 radius of each dsps headgroup , the greater the phase separation of 
 the bilayer .
figure 4 shows that the number 
 of dsps lipids with a dsps nearest neighbor is higher for systems 
 2 and 4 than it is for systems 1 and 3 .
this indicates that systems 
 2 and 4 separate to a greater extent by 1 billion collisions than 
 do systems 1 and 3 .
this is supported by the snapshots provided in 
 figure 3b , d , f , h .
normalized , time - averaged 
 percentage of dsps lipid molecules within 
 a 5  radius of each dsps headgroup ( coarse - grained site 1 on 
 each dsps lipid molecule ) .
the rate of phase separation was estimated by measuring the 
 change 
 in the average maximum height of the radial distribution function 
 over time between dsps head groups of coarse - grained types i ( a or 
 b ) and iii .
the larger the value of the maximum g(r ) for a given pair of coarse - grained types , the 
 closer , on average , those coarse - grained types would be to each other . 

the radial distribution function for coarse - grained types i ( a or 
 b ) and iii was chosen as our measure because this pair has the largest 
 interaction at low ph
. the maximum heights of each radial distribution 
 function were averaged over all replicates from each system .
table 3 provides the average maximum g(r ) at 3 different time points for systems 1 and 
 2 . at each time point , the g(r ) 
 was averaged over 1 million collisions . according 
 to table 3 , system 2 has a higher 
 average maximum g(r ) value than 
 that of system 1 at each time point .
this indicates that at each time 
 point system 2 had phase separated to a greater extent than that of 
 system 1 .
in addition , the change in the maximum g(r ) value between each time point was higher for 
 system 2 than that for system 1 .
this shows that system 2 separated 
 into heterogeneous domains at a faster rate than system 1 .
similar 
 to system 2 , systems 4 and 6 , which were also run at low ph , showed 
 an increase in the average maximum g(r ) with time .
however , systems 3 and 5 , which were run at neutral 
 ph , did not show a monotonic increase in g(r ) with time .
this shows that systems 4 and 6 are phase 
 separating to a greater extent and at a faster rate than are systems 
 3 and 5 . in order to learn
how bilayer surface area constraints 
 associated 
 with periodic boundary conditions impacts the phase separation in 
 21pc : dsps systems , we compare the results of simulations of systems 
 1 and 2 with those on systems 5 and 6 .
the goal was to determine if 
 allowing the bilayer to interact with its periodic image had an effect 
 on the rate at which the 21pc and dsps lipids separated into different 
 domains .
simulations on systems 1 , 2 , 5 , and 6 contained bilayers 
 with equimolar ratios of 21pc / dsps lipids , respectively .
each of these 
 simulations was started from an initial bilayer with an area of 60 
  60  .
the box lengths for simulations on systems 
 1 and 2 were 120  120  , whereas the box lengths 
 for simulations on systems 5 and 6 were 60  60  .
therefore , in simulations of systems 1 and 2 , the bilayers can not 
 interact with their periodic images , whereas in simulations of systems 
 5 and 6 , they can .
figure 3 provides snapshots 
 ( aerial images ) after 1 million and 1 billion collisions from one 
 of the replicates of the simulations that were run on systems 1 , 2 , 
 5 , and 6 . at 1 billion collisions , the area per lipid for systems 
 1 , 2 , 5 , and 6 was 46.2 , 46.2 , 48.0 , and 48.0  , 
 respectively
therefore , when the bilayers do not interact with their 
 periodic image , their surface area slightly decreases . in comparison 
 to the bilayer in system 1 ( no periodic boundary conditions , neutral 
 ph ) , the bilayer in system 5 ( periodic boundary conditions , neutral 
 ph ) appears to separate at approximately the same rate .
similarly , 
 the bilayer in system 2 ( no periodic boundary conditions , low ph ) 
 and the bilayer in system 6 ( periodic boundary conditions , low ph ) 
 appears to achieve about the same amount of phase separation . according 
 to figure 4 ,
systems 2 and 6 ( low ph ) have 
 the highest fraction of dsps lipid molecules within a 20  radius 
 of each dsps headgroup .
this indicates that systems 2 and 6 phase 
 separate to a greater extent than do systems 1 and 5 , which is supported 
 by the snapshots of each of these systems at 1 billion collisions 
 provided in figure 3 .
thus , we see that the 
 21pc / dsps bilayers at low ph separate faster than 21pc / dsps bilayers 
 at neutral ph regardless of whether the bilayers had a free surface 
 or periodic boundary conditions .
we 
 described the extension of lime , which was originally designed 
 for a single type of lipid ( dppc ) , at a temperature of 325 k , to be 
 applicable to dsps and 21pc lipids at a temperature of 310 k. lime 
 is an intermediate - resolution , implicit - solvent coarse - grained model 
 for phospholipid molecules , which is designed for use with dmd .
one 
 of the main advantages of combining lime with dmd is that it makes 
 it possible to sample much wider regions of conformation space and 
 longer time scales than that in traditional molecular dynamics . by 
 treating solvent implicitly , coarse - graining , and using a discontinuous 
 potential
, we can simulate the behavior of mixed lipid systems at 
 a much faster rate than that in traditional molecular dynamics simulations . 
 to our knowledge , our lime / dmd simulations are the first molecular - level 
 simulations that demonstrate ph - dependent phase separation in mixed 
 lipid bilayers .
the simulations that we describe in this work run 
 at a rate of approximately 0.017 cpu hours per 1 million collisions 
 ( intel xeon e5520 2.27 ghz ) .
the expanded version of lime was 
 used to simulate the behavior 
 of bilayers with and without surface area restrictions composed of 
 different molar ratios of 21pc and dsps at both neutral and low ph . 
 according to our simulations , 21pc / dsps systems separate faster at 
 low ph than at neutral ph .
these results are consistent with those 
 of sofou and co - workers , who found that the phase separation and the 
 formation of dsps - rich domains in lipid vesicles composed of 21pc , 
 dsps , 5 mol % cholesterol , and 1 mol % dspe - peg increases with a decrease 
 in ph .
the molecular mechanism driving 
 this phase separation is thought to be the following : ( 1 ) at neutral 
 ph ( 7.47.0 ) , the net negative charge on the dsps molecules 
 causes electrostatic repulsion between them , which hinders their ability 
 to phase separate , and ( 2 ) when the ph is lowered , the dsps head groups 
 become protonated and loose their negative charge , which allows these 
 molecules to hydrogen bond to each other and phase separate into heterogeneous 
 domains . the fact that our results are 
 consistent with those of sofou and co - workers
is not proof that lime / dmd 
 can be used to model the behavior of sofou s liposomes .
the 
 ultimate test will be to simulate the escape of drug molecules from 
 21pc / dsps liposomes at neutral and low ph using lime / dmd .
we found that bilayers interacting with their periodic images and 
 free surface bilayers phase separated to approximately the same extent 
 after 1 billion collisions .
bilayers interacting with their periodic 
 images are expected to behave more similarly to liposomes in experiments 
 because these bilayers will be promoted ( but not forced ) to maintain 
 a certain surface area and not to shrink or expand to a large extent . 
 in experiments
, liposomes composed of 21pc / dsps lipids in a solution 
 containing salt ions can not simply adjust to any desired size for 
 the following reason : as 21pc / dsps lipids begin to separate into heterogeneous 
 domains , the lipids try to pack closer together . as they do this , 
 the overall liposome size shrinks . during this process
, water molecules 
 can escape from the inner core of the liposome , whereas salt molecules 
 can not .
this results in an osmotic pressure that prevents the liposome 
 from shrinking to a large extent .
this scenario is best mimicked by 
 our simulations with bilayers interacting with their periodic images . to make better contact with the experiments on 21pc , dsps , and 
 cholesterol
, we are now calculating lime parameters for cholesterol 
 so that we can add it to our simulations to learn what role cholesterol 
 plays in the heterogeneous domain formation of mixed lipid systems . 

while our simulations were inspired by experiments done in the sofou 
 lab , we have not , as yet , included the specific molecular embellishments 
 associated with the sofou liposomes , including the targeting ligand , 
 fusion peptides , and cholesterol .
the multiscale modeling procedure that we used to 
 calculate lime parameters for dppc , 21pc , and dsps could , however , 
 be used to calculate parameters for these additional molecules . in 
 the future
one drawback of using dmd is that it is difficult to correlate 
 collision times with real time .
this is a result of the fact that 
 dmd simulations are event - driven rather than time - driven and that 
 solvent is modeled implicitly .
if information is available regarding 
 the experimental time required for a certain process or event to occur , 
 then this can be correlated to the simulation time required for the 
 same process or event to occur , allowing us to estimate a correlation 
 between simulation and real time .
additionally , treating solvent implicitly 
 prevents us from studying the diffusion and hydrodynamics associated 
 with the phase separation of the lipid membranes .
our model can not 
 be used to predict how solvent molecules will interact at the interface 
 of the bilayer with the 21pc and dsps head groups and affect the rate 
 of heterogeneous domain formation of these lipids .
another limitation 
 of our model is that we do not have a systematic way of correlating 
 temperature changes to real temperatures .
we know that by setting t * = 0.2 in our simulation the results seem to mimic the 
 behavior of 21pc / dsps systems at 310 k. multiscale modeling could 
 be used to determine the temperature dependence of the interaction 
 energies in our model , but this would require data from atomistic 
 simulations at various temperatures , which is beyond the scope of 
 this work . finally ,
although our lime / dmd simulations are currently 
 run using the nvt ensemble , we could perform npt simulations in the 
 future . to do this
dmd 
 simulations in which the volume change moves are made with monte carlo 
 and particle displacement moves are made with dmd .
this type of simulation 
 has been performed by previous members of the hall group in the past . despite these limitations
, we feel that our lime / dmd model does 
 a good job of modeling the behavior of liposomes observed experimentally 
 and that it could be a useful tool to researchers trying to study 
 the behavior of lipid systems on a molecular level .
 we 
 extend lime , an intermediate resolution , implicit solvent model 
 for phospholipids previously used in discontinuous molecular dynamics 
 simulations of 1,2-dipalmitoyl - sn - glycero-3-phosphocholine 
 ( dppc ) bilayer formation at 325 k , to the description of the geometry 
 and energetics of 1,2-distearoyl - sn - glycero-3-phospho - l - serine ( dsps ) and 1,2-dihenarachidoyl - sn - glycero-3-phosphocholine 
 ( 21pc ) and mixtures thereof at both neutral and low ph at 310 k. a 
 multiscale modeling approach is used to calculate the lime parameters 
 from atomistic simulation data on a mixed dppc / dsps system at different 
 ph values . in the model , 17 coarse - grained sites represent dsps and 
 18 coarse - grained sites represent 21pc . 
 each of these coarse - grained 
 sites is classified as 1 of 9 types . 
 lime / dmd simulations of equimolar 
 bilayers show the following : ( 1 ) 21pc / dsps bilayers with and without 
 surface area restrictions separate faster at low ph than at neutral 
 ph , ( 2 ) 21pc / dsps systems separate at approximately the same rate 
 regardless of whether they are subjected to surface area restrictions , 
 and ( 3 ) bilayers with a molar ratio of 9:1 ( 21pc : dsps ) phase separate 
 to form heterogeneous domains faster at low ph than at neutral ph . 
 
 our results are consistent with experimental findings of sofou and 
 co - workers ( bandekar et al . 
 mol . 
 pharmaceutics , 2013 , 10 , 152160 ; karve et al . 
 biomaterials , 2010 , 31 , 44094416 ) that more doxorubicin is released from 21pc / dsps liposomes at low 
 ph than at neutral ph , presumably because greater phase separation 
 is achieved at low ph than at neutral ph . 
 these are the first molecular - level 
 simulations of the phase separation in mixed lipid bilayers induced 
 by a change in ph . 
Introduction
Methods and Model
Results
Conclusions