/
Coarse Grained Model for Semiquantitative Lipid Simula Coarse Grained Model for Semiquantitative Lipid Simula

Coarse Grained Model for Semiquantitative Lipid Simula - PDF document

faustina-dinatale
faustina-dinatale . @faustina-dinatale
Follow
429 views
Uploaded On 2015-04-17

Coarse Grained Model for Semiquantitative Lipid Simula - PPT Presentation

Marrink Alex H de Vries and Alan E Mark Department of Biophysical Chemistry Uni ersity of Groningen Nijenborgh 4 9747 AG Groningen The Netherlands Recei ed August 22 2003 In Final Form October 29 2003 This paper describes the parametrization of a ne ID: 51042

Marrink Alex

Share:

Link:

Embed:

Download Presentation from below link

Download Pdf The PPT/PDF document "Coarse Grained Model for Semiquantitativ..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.


Presentation Transcript

CoarseGrainedModelforSemiquantitativeLipidSimulationsSiewertJ.Marrink,*AlexH.deVries,andAlanE.MarkDepartmentofBiophysicalChemistry,UniersityofGroningen,Nijenborgh4,9747AGGroningen,TheNetherlandsed:August22,2003;InFinalForm:October29,2003Thispaperdescribestheparametrizationofanewcoarsegrained(CG)modelforlipidandsurfactantsystems.Reductionofthenumberofdegreesoffreedomtogetherwiththeuseofshortrangepotentialsmakesitcomputationallyveryefficient.Comparedtoatomisticmodelsagainof34ordersofmagnitudecanbeachieved.Micrometerlengthscalesormillisecondtimescalesarethereforewithinreach.Toencourageapplications,themodeliskeptverysimple.Onlyasmallnumberofcoarsegrainedatomtypesaredefined,whichinteractusingafewdiscretelevelsofinteraction.Despitethecomputationalspeedandthesimplisticnatureofthemodel,itprovestobebothversatileinitsapplicationsandaccurateinitspredictions.Weshowthatdensitiesofliquidalkanesfromdecaneuptoeicosanecanbereproducedtowithin5%,andthemutualsolubilitiesofalkanesinwaterandwaterinalkanescanbereproducedwithin0.5oftheexperimentalvalues.TheCGmodelfordipalmitoylphosphatidylcholine(DPPC)isshowntoaggregatespontaneouslyintoabilayer.Structuralpropertiessuchastheareaperheadgroupandthephosphatephosphatedistancematchtheexperimentallymeasuredquantitiesclosely.Thesameistrueforelasticpropertiessuchasthebendingmodulusandtheareacompressibility,anddynamicpropertiessuchasthelipidlateraldiffusioncoefficientandthewaterpermeationrate.Thedistributionoftheindividuallipidcomponentsalongthebilayernormalisverysimilartodistributionsobtainedfromatomisticsimulations.Phospholipidswithdifferentheadgroup(ethanolamine)ordifferenttaillengths(lauroyl,stearoyl)orunsaturatedtails(oleoyl)canalsobemodeledwiththeCGforcefield.Theexperimentalareaperheadgroupcanbereproducedformostlipidswithin0.02.Finally,theCGmodelisappliedtononbilayerphases.Dodecylphosphocholine(DPC)aggregatesintosmallmicellesthatarestructurallyverysimilartoonesmodeledatomistically,andDOPEformsaninvertedhexagonalphasewithstructuralparametersinagreementwithexperimentaldata.1.IntroductionSolutionsoflipidsandsurfactantsdisplayaveryrichphasebehavior,includingmicellar,lamellar,hexagonal,andcubicphases.Thestructuralcharacteristicsofthesephasesandthedetailsofthephasediagramsstronglydependonthephysico-chemicalnatureoftheconstituents.Understandingtherelationbetweenmolecularstructureandaggregationbehavioristhereforeofgreatimportance.Computersimulationshaveprovidedausefultooltoelucidatethestructureoflipidphases,especiallythemicellarandlamellarphases.Atomisticsimula-tionsrevealmaximumdetailbutarerestrictedtosmalllengthandtimescales.Althoughspontaneousaggregationprocessescanalsobemodeledatomistically(e.g.,seeref1),amorecomputationallyefficientsimulationmodelisrequiredtofullyinvestigatetherichphasebehavioroflipidandsurfactantCoarsegrained(CG)models,inwhichsmallgroupsofatomsarerepresentedbysingleinteractionsites,arebecomingincreasinglypopulartostudysystemsoflipidsandsurfactants(forrecentreviews,seerefs2and3).CGmodelscanbeeitheron-latticeoroff-lattice.Whereasthefirstisfaster,theoff-latticetypemodelsaremoreversatileandrealistic.TheCGmodelthatwepresentinthispaperisanoff-latticemodel.Thefirstoff-latticeCGlipidmodelwasdevelopedbySmitetal.Intheirmodeltwodifferentinteractionsitesaredistinguished:eitherwater-likeoroil-like.TheyinteractviaaLennard-Jonestypeinteractionpotentialthatispurelyrepulsiveforwaterinteractionsorshortrangeattractiveforlikeparticles.Stringsoftheseparticlescansuccessfullyrepresentlipidandsurfactantlikemolecules.SimulationsbasedontheSmitmodelhavesubsequentlybeenappliedtostudytheformationandstructureofmicelles(e.g.,seerefs5and6)andbilayers(e.g.,seerefs7and8).Omittingexplicitsolventinthesimulationfurtherincreasestheaccessiblelengthandtimescalesandallowsstochasticsimulationsoftheconcentrationdependentmicelleorvesicleformation,forinstance.Anotherstochasticapproach,whichhasitsorigininsimulationsofblockisthemethodofdissipativeparticledynamics(DPD).InDPDthelipidsaremodeledassoftbeads,represent-ingfluidelementsratherthanrealparticles.Thebeadsinteractpairwiseviaacombinationofrepulsive,dissipative,andrandomforces.AswiththeSmitmodel,springsdefinemoleculararchitecture.TheDPDtechniquehasfoundrecentapplicationsintheareaoflipidsandsurfactants,withsimulationsofbilayerelasticproperties,poreformation,vesicleformationandfusion,andtheconstructionofacompletephasediagramforasimpleABtypesurfactant.AlthoughthecurrentCGmodelsusedinMDandDPDarebecomingverypowerfulinunderstandingstructuralaspectsoflipidandsurfactantphasesandtherelativephasestabilities,themodelscitedabovearequalitativeratherthanquantitativeintheirpredictions.Illustrativeofthisfactistheuseofdimension-lessunitstomeasurethelength,time,andenergyscales.Wherealinkismadetorealisticsystemsthemappingontophysicalmeasuresisdoneinhindsight.Thequestionis,canaCGtype *Correspondingauthor.E-mail:marrink@chem.rug.nl.J.Phys.Chem.B10.1021/jp036508gCCC:$27.502004AmericanChemicalSocietyPublishedonWeb12/02/2003 modelbespecificallyparametrizedtomodelarealisticlipidsysteminadvanceandcanitthenbeusedtomakequantitativepredictions?ArecentattempttodosowasmadebyShelleyetTheyparametrizedaCGmodeltoexplicitlymodelDMPC.LJinteractionsforsolventsitesandtailgroupsareoptimizedtoreproducethebulkpropertiesofwaterandalkanes.AtomisticsimulationsofDMPCbilayersareusedtooptimizetheparametersforthelipidheadgroupandtheanglepotentialsofthelipidtails.Electrostaticinteractionsarealsotakenintoaccountexplicitlyinthismodel.MonteCarloandMDsimula-tionswiththismodelshowtheassemblyofDMPCintobilayerswithpropertiesinsemiquantitativeagreementwithatomisticAlthoughsuggestingthatmorequantitativeresultscouldbeobtainedwithCGmodels,theShelleymodelhasseverallimitingfeatures:itisoptimizedforaparticularlipidinaparticularphaseonly,itusescomplicatedinteractionfunctionsrequiringshortintegrationtimesteps,anditevaluateslong-rangeinteractions.ToimproveontheCGmodelscurrentlyavailable,fouraspectsaredeemedimportant:speed,accuracy,applicability,andversatility.IntheCGmodelpresentedinthispaperthesefouraspectsareoptimizedsimultaneously:(i)Speedisobtainedbyincludingonlyshort-rangeinteractionsandbytheuseofsmoothpotentialssuchthatlargeintegrationstepscanbeused.(ii)AccuracyismaximizedbymatchingCGresultstoatomisticsimulationsasmuchaspossible,foravarietyofcomponentsandphasessimultaneously.Inthisprocessstructuralanddynamicalaswellasthermodynamicaldataareused.(iii)Applicabilityisenhancedthroughthesimplicityoftheforcefield,theuseofstandardinteractionpotentials,andfewparameters.Furthermore,theparametersarephysicallymean-ingfulandusedinaconsistentmanner.(iv)Versatilityisimpliedastheforcefieldleavesenoughroomtoaccommodatestructuraldetailofmolecules,andbecausethereisnorestrictiontothephaseofthesystem.Therestofthispaperisorganizedasfollows.Anextensivedescriptionofthemodelispresentedinthenextsection.Applicationsofthemodelarethenshownforanumberofsystems:alkane/water,asaltsolution,lipidbilayers,andnonbilayerphases.Ashortdiscussionfollowsinwhichthemeritsandshortcomingsofthemodelareexamined.2.Model2.1.InteractionSites.Afour-to-onemappingisusedtorepresentthemoleculesinthesimplifiedmodel;i.e.,onaverage,fouratomsarerepresentedbyasingleinteractioncenter.Thisruleisnotstrict,assometimesitisappropriatetomapthree,five,ormoreatomsintooneinteractioncenter.Becauseoftheirsmallsizeandmass,hydrogenatomsarenotconsideredatall.Tokeepthemodelsimple,wecurrentlyconsiderfourmaintypesofinteractionsitesonly:polar(P),nonpolar(N),apolar(C),andcharged(Q).Polarsitesrepresentneutralgroupsofatomsthatwouldeasilydissolveintowater(e.g.,ethyleneglycol),apolarsitesrepresenthydrophobicmoieties(e.g.,butane),andnonpolargroupsareusedformixedgroupswhicharepartlypolar,partlyapolar(e.g.,propanol).Chargedsitesarereservedforionizedgroups(e.g.,ammonium).ForparticlesoftypeNandQfoursubtypes(0,d,a,andda)arefurtherdistinguished.Thesubtypesallowfine-tuningoftheinteractionsonthebasisofthechemicalnatureoftheatoms,whicharerepresentedbytheCGgroups.Subtype0appliestogroupsinwhichnohydrogenbondingcapabilitiesexist,dandatogroupsthatcouldactasahydrogenbonddonororacceptor,respectively,anddatogroupswithbothdonorandacceptoroptions.Realisticmassescanbeassignedtotheparticles,butforreasonsofcomputationalefficiency,thesamemassescanbeused.Intheapplicationsdescribedhereamassof72amu(correspondingtofourwatermolecules)isassignedtoeachsiteunlessotherwisestated.ThemappingofsomeofthemoleculesstudiedinthispaperisdepictedgraphicallyinFigure1.2.2.NonbondedInteractions.ThenonbondedinteractionsbetweeninteractionsitesaredescribedbytheLennard-Jones(LJ)potentialrepresentingtheeffectiveminimumdistanceofapproachbetweentwoparticlesandthestrengthoftheirinteraction.Thestrengthoftheinteraction,determinedbythevalueofissummarizedinTable1forallpossibleinteractionpairs.Notethatonlyfivelevelsofinteractionaredefined:attractive(I,5kJ/mol),semiattractive(II,4.2kJ/mol),intermediate3.4kJ/mol),semirepulsive(IV,2.6kJ/mol),andrepulsive(V,1.8kJ/mol).ThelevelIinteractionmodelsstrongpolarinteractionsasinbulkwater,levelIIImodelsthenonpolarinteractionsinaliphaticchains,andlevelVmodelsthehydrophobicrepulsionbetweenpolarandnonpolarphases. Figure1.Mappinginthecoarsegrainedmodelforwater,ions,butane,hexadecane,DPC(dodecylphosphocholine),andDPPC(dipalmi-toylphosphatidylcholine).AllatomsexceptnearestneighborsinteractthroughaLennard-Jonespotential.Fourmainatomtypesaredistin-guished:polar(P),nonpolar(N),apolar(C),andcharged(Q).Nearestneighborsareconnectedbyaweakharmonicspring,nextnearestneighborsinteractthroughaharmonicanglepotential.Chargedgroupsalsointeractthroughashort-rangeelectrostaticpotential.Seetextfor TABLE1:InteractionMatrix PNCQgroupsubtype0dada0dadaPIIVIIIIIIIIVIIIIN0IVIIIIIIIIIIIIIIIIIIIIIIIIIIIdIIIIIIIIIIIIIVIIIIIIIIIIaIIIIIIIIIIIIIVIIIIIIIIIIdaIIIIIIIIIIVIIIIIIIICVIIIIVIVVIIIVVVVQ0IIIIIIIIIIIIIVIIIIIIIIIIIdIIIIIIIIIIIVIIIIIIIIIaIIIIIIIIIIIVIIIIIIIIIdaIIIIIIIIIVIIIIILevelofinteractionI(attractive),II(semiattractive),III(intermedi-ate),IV(semirepulsive)orV(repulsive).Fourdifferentgroupsareconsidered:polar(P),nonpolar(N),apolar(C),andcharged(Q).BothgroupsNandQhavefoursubtypes:0fornohydrogenbondingcapabilitiespresent,dforgroupsactingashydrogenbonddonor,aforgroupsactingashydrogenbondacceptor,anddaforgroupswithbothdonorandacceptoroptions. r)12-(óij CGModelforSemiquantitativeLipidSimulationsJ.Phys.Chem.B,Vol.108,No.2,2004 NotethatthelevelVinteractionisactuallystillattractivebutislabeledrepulsiveasitresultsinthedemixingofwaterandoilphases.LevelsIIandIVareofintermediatestrength.Forallfiveinteractiontypesthesameeffectivesizeisassumed,0.47nm.InthesimulationstheLJinteractionpotentialiscutoffatadistance1.2nm,correspondingtoap-proximately2.5.Toreducethecutoffnoise,theLJpotentialissmoothlyshiftedtozerobetweenadistance0.9nm.WiththestandardGromacsshiftfunctionboththeenergyandtheforcevanishatthecutoffdistance.Analkane/watersystemwasusedasareferencesystemtocalibratethelengthandenergyscalesoftheLJparameters.Atrialanderrorprocedurewasusedtosimultaneouslyoptimizefourbasicparameters,namelythesizeoftheparticlesplustheinteractionenergiesbetweenpolarparticles(levelI),apolarparticles(levelIII),andbetweenpolarandapolarparticles(levelV).Ourgoalwastoreproducetheexperimentaldensitiesofpurewaterandalkanesystemsaroundroomtemperature,themutualsolubilityofoilandwater,andtherelativediffusionrates.Resultsarepresentedintheapplicationsection.InadditiontotheLJinteraction,chargedgroups(typeQ)interactviathenormalelectrostaticCoulombicpotentialwithrelativedielectricconstant20forexplicitscreening.Thispotentialisalsoshiftedsmoothlyfrom0.0nmtothesamecutoffdistanceasusedfortheLJinteractions,1.2nm.ThestandardshiftfunctionofGromacsisusedinwhichboththeenergyandforcevanishatthecutoffdistance.Shiftingoftheelectrostaticpotentialinthismannermimicstheeffectofadistance-dependentscreening.NotethatinteractionsitesoftypeQareintendedforgroupsbearingfullorclosetofullchargesonly.InteractionsarisingfromsmallpartialchargesarerepresentedbytheLJpotential.Thechargesinprinciplerepresentthetruechargeofagroup;however,inthecaseofsmallhydratedatoms(e.g.,ions)areducedchargeisusedtotakeintoaccounttheeffectofanimplicithydrationshell.2.3.BondedInteractions.Bondedinteractionsbetweenchemicallyconnectedsitesaredescribedbyaweakharmonic)withanequilibriumdistanceTheLJinteractionisexcludedbetweenbondedparticles.Bondedparticles,onaverage,aresomewhatclosertoeachotherthanneighboringnonbondedparticles(forwhichtheequilibriumdistanceis2).Theforceconstantoftheharmonicbondingpotentialis1250kJmol.Thisforceconstantallowsconsiderabledeviationsfromtheequilibriumbondlength15%)atthecostof.Torepresentchainstiffness,aweakharmonicpotential)ofthecosinetypeisusedfortheThebasicequilibriumbondangle,withaforceconstantof25kJmol.Suchasmallforceconstantallowsanangledeviationof30atthecostof.Notethatthisanglereproducesthepropertiesofaliphaticchains.Foranglesofnonaliphaticnaturecomparisontoatomisticmodelsisusedtooptimizetheparameters(seethevariousapplicationsbelow).LJinteractionsbetweensecondnearestneighborsarenotexcluded.2.4.SimulationParameters.Themodelisdesignedforuseclosetoroomorphysiologicaltemperatures.Thetemperaturedependenceofthemodelhasnotbeeninvestigatedindetail.Standardcouplingschemescanbeusedforbothtemperatureandpressure.Thelargestfeasibleintegrationtimestepformostsystemsisd50fs,butsometimesasmallertimestepisrequiredforstability(3040fs).Theneighborlistcanbeupdatedevery10stepsusinga1.2nmneighborlistcutoff.AlthoughthevariablesoftheCGsystem(densities,lengthscales,energies,temperature,pressure)keeptheirphysicalmeaning,thisisnotstrictlytrueforthetimescale.ThedynamicsarefasterbecausetheCGinteractionsaremuchsmoothercomparedtoatomisticinteractions.OnthebasisofcomparisonofdiffusionconstantsintheCGmodelandinatomisticmodels,theeffectivetimesampledusingCGis36timeslarger.NotethatthisfactoraffectsALLthedynamicspresentinthesystem.Therelativedynamicspresentwithinthesystemappearstobewellpreserved(withinafactorof2).WheninterpretingthesimulationresultswiththeCGmodel,onecantoafirstapproximationsimplyscalethetimeaxis.Thestandardconver-sionfactorweuseisafactorof4,whichistheeffectivespeedupfactorinthediffusionaldynamicsofCGwatercomparedtorealwater.Thesimulationtimesreportedintheremainingofthepaperareeffectivetimes,unlessotherwisestated.ThesimulationswereperformedwiththeGromacssimulationpackageversion3.0.Foratypicalsystemcontaining100000CGparticles,arateof10000timesteps(2nseffectivetime)perCPUhouronadualPentium1GHznodeisachieved.Theparametersandexampleinputfilesoftheapplicationsdescribedinthispaperareavailableathttp://md.chem.rug.nl/3.Applications3.1.Alkane/Water.ThepropertiesofbulkwatercanbereproducedbycoarsegrainingfourwatermoleculesintooneLJbeadoftypeP.At300Kand1baradensityof1gcmisobtainedforasystemcontaining400CGparticles(1600realwaters).Fromthevolumefluctuationsofthissystemtheisothermalcompressibilitycoefficientisfoundtobe6.0,similartotheexperimentalvalue).Theself-diffusionconstantoftheCGwatersitesat300Kis.TheCGwatersites,however,representthecenter-of-massoffourrealwatermolecules.Theaveragemeansquareddisplacementofthecenter-of-massoffourmoleculesis4timeslessthantheaveragemeansquareddisplacementoffourindependentlydiffusingmolecules.effectivediffusionconstantofindividualwatermoleculesasrepresentedbytheCGparticlesistherefore4timeslarger.Theself-diffusionrateofwaterasmodeledbytheCGparticlesistherefore2.Forpurewatertheexperimentaldiffusionconstant(at300K).TodeterminethefreezingpointoftheCGwater,asmallicecubewassimulatedsurroundedbyliquidwater.Atatemperature290Kthesolidphaseappearsinequilibriumwiththeliquidphase.Belowthistemperaturetheicecubeinducesfreezingofthewholesystem,andabovethistemperaturetheicecubemelts.ThemeltingtransitiontemperatureoftheCGwatermodelis5K,slightlyhigherthanforrealwater(273K).Spontaneousfreezingofliquidwaterisonlyobservedfortemperaturesbelow240K,however.WithinthetemperaturerangeofrealliquidwatertheCGwaterthereforealsoremains 4ð0rr(2)Vbond(R))1 2Kbond(R-Rbond)2(3)Vangle(õ))1 J.Phys.Chem.B,Vol.108,No.2,2004Marrinketal. inafluidstate.Belowatemperatureof290KtheCGwaterisinasupercooledstate.BulkpropertiesofliquidalkanescanbereproducedusingstringsofconnectedCparticles(eachCparticlerepresentingfourmethyl/methylenegroups)withstandardbondandanglepotentials(seeabove).Wehaveperformedbulkalkanesimula-tionswithoursimplifiedforcefieldforbutane,octane,dodecane,hexadecane,andeicosane(1through5Cparticles,respectively).Thesystemscontained100alkanemoleculesexceptforoctane(200)andbutane(400).Thesimulationswereperformedat300Kand1barusingaParrinelloRahmancouplingscheme.Thedensities,compressibilities,andself-diffusionratesobtainedfromthesimulationsarepresentedinTable2.Inthistablewealsopresentresultsfortheintermediatealkanessuchashexane,decane,etc.,whichdonothaveamultipleoffourmethyl/methylenegroups.ThepropertiesofthesemoleculescanbestbeapproximatedbyassumingmappinglessthanfouratomisticgroupsintoaCGsite;e.g.,decaneismodeledas3Cparticlesimplying3.3to1ratherthan4to1mapping.ThusasimulationofCCrepresentsdecaneaswellasdodecane.OnlythemassesassignedtoeachCGparticlearedifferent,resultingindifferenteffectivedensitiesanddiffusionrates.Especiallyforthelongeralkanes(dodecaneandup)theexperimentaldensitiesarewellreproduced(within5%)withtheCGmodel.Pushingthelimitsofacoarsegrainedmodel,fortheshorteralkanesthedifferenceincreasesto10%.TheexperimentalisothermalcompressibilitiesarealsoapproximatedwiththeCGmodeloverthewholerangeofchainlengths.Thediffusionconstantsforthecoarsegrainedalkanesystemsareobtainedfromtheslopeofthemeansquaredisplacementasafunctionoftime,usingtheeffectivetimescalesetbythediffusionrateofpurewater(seeprevioussection).Theself-diffusionofthealkanesappears3timesslowerthanexperimentallydetermined.EspeciallyforthelongeralkanestheagreementbetweenthemodelandrealalkanesissatisfactoryforaCGmodel.ThechainstiffnessoftheCGalkanesiscomparabletothechainstiffnessofalkanesmodeledatomistically.ForcomparisonweperformedatomisticsimulationsofalkanesystemswiththesamenumberofmoleculesandsimulationconditionsasfortheCGmodel.ForthesesimulationsweusedthestandardGromacsforcefield,whichisfullyatomisticexceptforthemethylandmethylenegroupsthataretreatedasunifiedatoms.IntheCGmodel,theaverageCCangleinabulkphaseofstringsofthreetofiveCsitesis140withastandarddeviationof(onlyveryweaklychainlengthdependent).Atomisticsimulationsofdodecane,hexadecaneandeicosanegive(136,takingthecenter-of-massoffouradjacentmethylenegroupsasthecorrespondingpositionoftheCGsite.Acisdoublebondcanbeeffectivelymodeledbychangingtheequilibriumbondanglefrom180to120.Thechainstiffnessisalsoincreasedslightly(from25to35kJmol).Inasimulationofstringsof5CparticlesinwhichthemiddleCangleismodeledasadoublebondtheaveragecentralangleis122.Forcomparison,inanatomisticsimulationofoctadecenethecorrespondingangleis131Apartfromacorrectrepresentationoftherelativedensitiesoftheoilandwaterphases,itisimportantininterfacesystemstoreproduceoil/watermutualsolubilities.Tomodelthealkane/waterinterface,randomlymixedsystemscontainingeither100hexadecane(C4)or400butane(C)moleculesand1600watermolecules(400Pparticles)wereprepared.Thesimulationswererunfor1s(CG,effectivesimulationtime)or5ns(atomistic),300Kandanisotropicpressureof1bar.Allsystemsquicklyphaseseparatetoformanaqueousslabandanoilslab.InFigure2wecomparetheresultingrelativedensitiesofhexadecaneandwaterobtainedfromsimulationswithsimplifiedandwithatomisticforcefields.Thedensityprofilesareverysimilar.IfthedensitywithinaCGsitewouldberepresentedmorerealisticallyassmearedoutoveritsentirevolume,thematchingisevenbetter.Theinterfacesforthebutane/watersystemslookverysimilar(notshown).Thesolubilityofwaterinoilcanbecomputedstraightfor-wardlyfromtheequilibriumdensityofwaterintheoilphase.Thefreeenergyofpartitioningofwaterbetweentheaqueousandtheoilphaseisgivenby),whereisthedensityofwaterintheaqueousphase.Withthesimplifiedmodel,simulationscanbeeasilyextendedintothemicrosecondrange,enoughtoobtainstatisticallyreliableresultson.Thepartitioningfreeenergyofwaterinhexadecaneforoursimpli-fiedmodelwascalculatedtobe2kJmol25kJmolLikewiseonecouldcomputethesolubilityofhexadecaneinwater.Thesolubilityoflongeralkanesinwaterissolow(50kJmolhowever,thatevenwiththesimplifiedmodelspontaneouspartitioningintotheaqeuousphaseisnotobserved.Ofcourseonecoulduseothersimulationmethodstocomputethefreeenergyofsuchaprocess,butitwasdeemedmoreusefultolookatthebutane/watersystem.Theexperimentalpartitioningfreeenergyforbutaneinwateris21kJmoltheobservedequilibriumdensityofbutaneinwaterinourCGsimulation,weestimate2kJmol,alsoingoodagreementwiththeexperimentalresult.BoththesolubilityofwaterinalkaneandofanalkaneinwatercanbeaccuratelyreproducedintheCGforcefield.3.2.SaltSolution.Pairdistributionfunctionsofsodiumchlorideatphysiologicallyrelevantconcentrations(0.11.0M)canbeapproximatedwithQparticlesbearingareducedcharge TABLE2:PropertiesofWaterandAlkaneswiththeCGModelComparedtoExperimentalValues systemCGmodelgcmwaterP0.99(0.99)6(4.5)2.0(2.3)butaneC0.68(0.58)28(17)1.9(hexaneCC0.58(0.66)14(17)0.7(4)octaneCC0.77(0.70)14(13)0.6(2)decaneCC0.67(0.73)12(11)0.35(1)dodecaneCC0.80(0.75)12(10)0.3(tetradecaneCC0.71(0.76)12(9)0.25(hexadecaneCC0.81(0.77)12(9)0.2(octadecaneCC0.74(0.78)11()0.2(0.3)eicosaneCC0.82(0.79)11()0.15(Propertiesat300K,unlessspecified.Experimentalpropertiesaregiveninparentheses.Experimentaldensitiesat293K.mentalisothermalcompressibilitiesfromref20.ThevaluesfromsimulationsarecomputedfromthevolumefluctuationsinanNPTDiffusionrateswereobtainedfromtheslopeofthemeansquareddisplacement(MSD)curveinthelongtimelimit.Experimentalvaluesextrapolatedfromtemperature-dependentdata. Figure2.Hexadecane/waterinterface,fromatomistic(dashed)andsimplified(solid)simulations.CGModelforSemiquantitativeLipidSimulationsJ.Phys.Chem.B,Vol.108,No.2,2004 0.7.Themagnitudeofthisreducedcharge(0.7)ismerelyaparameter,optimizedtomodelthepairdistributionfunctionsobtainedfromatomisticsimulations.Inphysicalterms,thereducedchargemimicstheimplicitscreeningofthefirsthydrationshell;i.e.,theCGparticlerepresentsahydratedion.Experimentally,thefirsthydrationshellofbothsodiumandchloridecontainsaroundsixwatermolecules.SixhydrationwatersarethereforeconsideredtobeimplicitintheCGions.Theimplicithydrationshellimpliesstronghydrogenbondingpossibilities.Theionsarethereforemodeledassubtypeda.Giventheshortrangenatureoftheelectrostaticinteractioninourmodel,thecoarsegrainedionsareexpectedtomodelrealisticionsinthelimitofmoderatetohighionicstrengthonly(above0.1M).Althoughoptimizedforsodiumchloride,theCGionscouldbeusedforothersaltsaswell.InFigure3cumulativepairdistributionfunctionsobtainedwiththeCGmodelarecomparedtothoseobtainedfromareferenceatomisticsimulation.Theatomisticsystemcontained16sodiumchloridepairsand2135watermolecules,aconcen-trationof0.4M.A10nssimulationwasperformedatKusingPMEtoaccountforthelong-rangeelectrostaticinteractions.TheCGsystemofthesamesolutionconsistsof32Qparticles,16ofwhichcarryapositivechargeof0.7and16anegativeoneof0.7.ThenumberofsolventparticlesintheCGsolutionis485.AssumingfourwatersperCGwaterparticleandsixhydrationwatersimplicitintheQparticletheconcentrationisalso0.4M.TheCGmodelwassimulatedfors.Figure3showsthatthecoordinationnumbersforbothionandionwaterpairsobtainedfromatomisticsimula-tionscanbecloselyreproducedbytheCGmodelfordistanceslargerthan0.6nm.ForsmallerdistancestheatomisticnatureoftheparticlesbecomesimportantwhichtheCGmodelisunabletomodel.Similaragreementisfoundforionconcentra-tionsintherange0.11.0M.Apparently,withinthisrange,theelectrostaticscreeningissufficientlyhighthatlong-rangeinteractioncanbeneglected.3.3.LipidBilayers:DPPC.3.3.1.CGModel.Thephos-pholipidDPPCismodeledusing12CGsites(seeFigure1).Thetwotails(representing15methyl/methylenegroupspertail)arebasedontheCGmodelforhexadecane.Theglycerolesterlinkageismodeledbytwononpolarparticles,typeN.ThesubtypeawaschosenbecauseofthehydrogenbondacceptorcapabilitiesofthecarbonyloxygensrepresentedbytheseCGparticles.ThezwitterionicPCheadgroupismodeledbyapositivelychargedQparticlerepresentingthecholinemoiety,andanegativelychargedQparticlerepresentingthephosphategroup.Thepossibilityforthephosphategrouptoactasahydrogenbondacceptorrequiresmodelingassubtypea,whereasthecholinegrouplackshydrogenbondingoptionsandisthereforemodeledassubtype0.Theanglepotentialsintheheadgroupregionwereoptimizedtoreproducethevaluesobtainedinatomisticsimulations.ForthetripletsGLYC-C1-C2andPO4-GLYC-C1thesameanglepotentialsareusedasforthelipidtails,i.e.,anequilibriumangleof180andaforceconstantof25kJmol.Ananglepotentialwithasmallerequilibriumangleof120wasusedtomodeltheglycerolbackbonePO4-GLYC-GLYC.Theforceconstantwasagain25kJmol3.3.2.SpontaneousAggregation.Placedrandomlyintosolu-tion,thecoarsegrainedDPPCmoleculesspontaneouslyformbilayers.Figure4showsanexampleofthisaggregationprocessforasystemcontaining1600DPPCmoleculesand60000CGwaters.Adefectfreebilayerisformedafter200nssimulation.Forsystemscontainingfewerlipidsthebilayerformationprocessisfasterandcomparabletotherateofformationobservedforatomisticmodels.Systemswith128lipids,forinstance,self-assemblewithintensofnanoseconds.Theinter-mediatestagesofbilayerformationarealsosimilartothestagesfoundusingatomisticmodels:arapidlocalclusteringoflipidsintoirregularlyshapedaggregatesonthesubnanosecondtimescalewhichcoalescetoformacontinuouslipidaggregate.Thisaggregatesubsequentlytransformsitselfintoabilayerconfig-uration,stillcontainingdefects.Waterporesappearmetastableontimescalesrangingbetween3and25ns.Incontrasttoatomisticallysimulatedbilayerassembly,withtheCGapproachanothertypeoflonglivingdefectisobserved,namelyalipidbridge(indicatedbythearrowinFigure4).Alipidbridgeisanelongatedstalkpassingthroughthewaterphaseconnectingtwobilayers,ormorepreciselythebilayertoitsperiodicimage.Itisthetopologicalequivalentofawaterporethatconnectsthewaterphasesonbothsidesofthemembranerunningthroughthemembrane.Forsmallsystemstheselipidbridgesarestable Figure3.Cumulativepairdistributionfunctionsfora0.4MNaClsolution.SolidandopensymbolsdenotedataobtainedwiththeCGmodelandatomisticmodel,respectively.Circlesshowsodiumtrianglessodiumchloride,squaressodiumsodium,anddiamondschloridedistributionfunctions.NotethatwithintheCGmodelsodiumandchlorideareindistinguishableexceptfortheoppositecharge.CGwatersarecountedasfourrealwatermolecules.SixhydrationwatersareconsideredtobeimplicitintheCGions. Figure4.SpontaneousaggregationofDPPClipidsintoabilayer.Lipidheadgroupsarecoloredyellow,lipidtailspurple,waterblue.Thesimulationboxisindicatedwithgraylines.Thearrowpointsatametastablelipidbridge(orelongatedstalk).J.Phys.Chem.B,Vol.108,No.2,2004Marrinketal. upto20ns,typicallyslightlylessstablethanthewaterpores.ForlargersystemssuchastheoneshowninFigure4,however,itcantakehundredsofnanosecondsbeforeruptureoccurs.3.3.3.StructuralProperties.Oncethesemetastabledefectshavedisappeared,thebilayerrelaxesintoitsequilibriumstructure.TheequilibriumareaperlipidforDPPCintheCGmodelat323Kis0.64nm.Thismatchestheexperimentalestimateof0.64nmFigure5comparestheelectrondensityprofilesofanequilibratedcoarsegrainedDPPCbilayertoprofilesobtainedforonesimulatedatomistically.Bothdensitydistributionswereobtainedfromasimulationofabilayerpatchcontaining256DPPClipidsatfullsaturation(46watersperlipid).TheatomisticsimulationwasperformedusingastandardatomisticlipidforcefieldandwithPMEforthelong-rangeelectrostaticinteractions.TocomparetotheCGdistributions,thecenter-of-massoftheatomisticgroupscorrespondingtoaCGparticlewereusedtocalculatethedensities,exceptinthecaseofwaterforwhichthisisnotpossible.BoththeCGandatomisticsimulationswereperformedatthesametemperature323K)withthelateralandnormalpressuresindependentlycoupledtoapressureof1bartomimicconditionsofzerosurfacetension.ThepressurecouplingschemewasParrinelloThepeaksintheelectrondensitydistributionscoincidetowithin0.1nm,foreachofthemembranecompo-nents.ThewaterpenetratesthebilayertoalargerextentintheatomisticsimulationsthanintheCGsimulations,butthisismainlyaresultofrepresentingfourrealwatermoleculesbythecenterofasingleCGsolventparticle.IfthedensityweresmearedoutequallyacrossthediameteroftheCGparticle(almost0.5nm),thedepthofwaterpenetrationbecomesverysimilarinbothcases.TheinterfaceintheCGmodelisslightlymorestructuredthanintheatomisticmodel.ThethicknessoftheCGbilayer,measuredfromthepeaksofthephosphatedistribution,measures4.00.1nm,closetotheexperimentallydeterminedbilayerthicknessof3.85nmforthelamellarphaseofDPPCintheliquid-crystallinephase.Theclosecor-respondencebetweenthestructureoftheCGbilayerandtheatomisticoneisfurtherillustratedinFigure6,whichshowsanalignmentofsnapshotsofbothsystems.Forfurthercomparison,inFigure7lipidorderparametersareshownforboththeCGandtheatomisticDPPCbilayer.Thesecond-rankorderparameter(3cos1)wascomputedforconsecutivebonds,withtheanglebetweenthedirectionofthebondandthebilayernormal.Perfectalignmentwiththebilayernormalisindicatedby1,perfectanti-alignmentwith0.5,andarandomorientationwith0.Fortheatomisticsimulationthecenter-of-massoffouradjacentmethylenegroupswastakenasthepositioncorre-spondingtothepositionofaCGsite.Theprofilesareverysimilar.Boththephosphatecholinebondandtheglycerollinkagehaveapredominantlyparallelorientationwithrespecttothesurfacenormal,whereastheotherbondspreferaligningalongwiththesurfacenormal.Towardtheendofthetails,theorderdecreases.Notethattheorderparametersobtainedherecannotbedirectlycomparedtoexperimentalbondorderparameters,astheinformationaboutatomisticbondsisnotpresentintheCGmodel.However,goodagreementbetweenexperimentalandatomisticallysimulatedorderparameterprofileshavebeenshowninthepast.3.3.4.ElasticProperties.Apartfromtheabilitytoreproduceimportantstructuralquantities,collectiveelasticpropertiessuchasthebilayerbendingmodulus,theareacompressibility,andthelinetensionareimportantparameterstojudgethequalityofthemodel.Thebendingmoduluscanbeobtainedfromtheundulationspectrumofalargepatchofbilayer,asdescribedin Figure5.Comparisonofelectrondensitydistributionsobtainedfromatomistic(solidlines)andcoarsegrained(dashed)simulations.Thewaterdensityisshownincyan,thedensityofNC3groupsinblue,groupsinorange,theglycerolbackboneinpurple,andtheterminaltailgroupsingreen.Thetotaltaildensityisshowninblack.Thebilayercenterislocatedat0nm. Figure6.SnapshotsofaDPPCbilayersimulatedwithanatomisticmodel(left)andwiththeCGmodel(right).Thetopimagesshowsideviewsofthebilayer,thebottomimagesshowtopviews(withthewateromitted).Thecholinegroupsarecoloredblue,phosphatidylgroupsyellow,theglycerolbackbonered,thetailsgray.Adarkershadeofgrayisusedtodistinguishtheterminalmethylgroup.Wateriscolored Figure7.orderparameterofconsecutivebondswithrespecttothebilayernormal.ResultsobtainedfortheCGbilayeraremarkedwithcircles,thosefortheatomisticbilayerwithsquares.Dataareaveragedoverbothtails.CGModelforSemiquantitativeLipidSimulationsJ.Phys.Chem.B,Vol.108,No.2,2004 ref28.Forthispurposethesmallbilayercontaining256lipidswascopiedfivetimesinbothlateraldirectionstogenerateapatchconsistingof6400lipidswithlateraldimension45nm.Thesystemwassimulatedfor250nsatthesameconditionsasdescribedabove.Ittook100nsfortheundulatorymodestofullydevelop.Figure8showstheundulatoryspectrumobtainedfromthefinal150nssimulationofthissystem.Theintensityoftheundulationswithamplitudebehaviorinthelongwavelengthregime,aspredictedbyacontinuum.Hereisthemembraneareadenotesthebendingmodulus.Fittingthisequationtothedataobtainedfromthesimulationresultsinabendingmodulusof4J,wheretheuncertaintiesaredeterminedfromdifferentfittingprocedures.Thisvalueisclosetothevalueof5JobtainedforanatomisticallysimulatedDPPCbilayerandalsotothevalueof(5.6JobtainedexperimentallyforDMPC.Theareacompressibilitymoduluscanbecalculatedfromthefluctuationsinthemembraneareaperlipid,wheredenotesthenumberoflipidspermonolayeranddenotestheequilibriumarea.Forthelargebilayerpatchcontaining6400lipidstheareacompressibilitymodulusisfoundtobe40mN/m.Forthesmallersystemconsistingof256lipidsthevalueishigher,30mN/m.ThedifferenceisduetothecontributionofundulatorymodesinthelargesystemwhichlowertheareaThevalueobtainedforthelargesystemcompareswellwiththeexperimentalvalueforDPPCwhichisreportedtobe20mN/m.Itisreasonabletoassumethattheexperimentalvaluealsocontainscontributionsfromundulatorymodes.Notethatatomisticsimulationsperformedonsmallbilayerpatchesusuallyalsoreportvaluesthatareafactortwolargerthantheexperimentalestimate(e.g.,seerefThelinetension,oredgeenergy,ofthelipidmembranecanbeobtainedfromthecriticaltensionatwhichporescanbestabilizedinsideamembrane.Accordingtoatheoreticalmodeltheenergyofaporeinsideamembraneisgivenbythe,wheredenotestheporeradius,thelinetension,andthesurfacetension.Atlowtensiontheporesareunstable.Theedgeenergydrivestheclosureofthepore.However,atacriticaltension*theedgeenergyisovercome.Porescanthenbestabilizedat.Forlargertensionsporegrowthisunlimited,leadingtomembranerupture.Thetension-dependentopeningandclosingofaporeinthemembranehasrecentlybeensimulatedforDPPCusinganatomisticmodel.Itwasfoundthatporesofradius0.9nmappearedstableatacriticaltensionof40mN/m(dependingondetailsoftheforcefield),correspondingtoalinetensionof.ThisvalueisofthesameorderofmagnitudeastheestimatesfromexperimentsusingelectroporationtechniquesfornaturallecithinsandforSOPC,andfrommechanicalruptureeventsDOPCandSOPC.Repeatingthesimulationsasdescribedinref33,usingtheCGforcefield,quantitativelysimilarresultsarefound.Largetensionsarerequiredtospontaneouslyinduceporeformation,immediatelyfollowedbymembranerupture.Atatensionof100mN/maporeappearswithinnanoseconds;atatensionof65mN/mittakesseveralmicroseconds.Atsmallertensions(50mN/m)themembranethins,butnoporeformationisobservedontimescalesupto50.Thethinningofthemembraneatlowtensions(510mN/m)iscomparabletothatobservedexperimentally(1PreexistingporesinthemembranecanbestabilizedatacriticaltensionofmN/m.Theporesaresimilarinshapetothosereportedfortheatomisticsimulations,i.e.,anhourglassshapedwaterchannellinedbylipidheadgroups.Figure9showsacloseupofsuchaporestabilizedinacoarsegrainedDPPCmembrane.Theradiusoftheporefluctuatesbetween1.5and2.0nm,fromwhichthelinetensionisestimatedtobe.Thisvalueforthelinetensionisslightlyhigher,butofthesameorderofmagnitude,astheexperimentalvaluescitedabove.3.3.5.DynamicalProperties.Asdiscussedabove,therateofspontaneousaggregationoflipidsintobilayersusingtheCGmodelissimilartotheratesfoundforatomisticsimulations.ThelateraldiffusionconstantofthelipidscanbeobtainedfromthelimitingslopeoftheMSDcurve.Microsecondsimulationsofthebilayercontaining256DPPCmoleculesallowstheobservationoftruelylongtimediffusivebehavior.AtK,thelateraldiffusionratesofDPPCequals(3,ofthesameorderofmagnitudeasexperimentallymeasured(valuesaretypicallyreportedtobearound1attemperaturescloseto323K;e.g.,seerefs38andAnotherinterestingdynamicquantityforwhichexperimentaldataexist,istherateofpermeationofwateracrossthebilayer.Duringmicrosecondsimulationsofa256DPPCbilayerpatchthespontaneouspermeationofwatermoleculesisobserved.Themagnitudeoftheunidirectionalfluxis1CGwatermoleculeper100ns,correspondingto4realwaters.Fromthisfluxthepermeationratecanbecomputedascms,where55.5moldmistheeffectivewaterconcentrationanddenotesthemembranearea( Figure8.Spectralintensityforundulatorymodesasafunctionof.Thesolidlinerepresentsafittothedataatlongwavelengths,obeying Figure9.CloseupofaporestabilizedinacoarsegrainedDPPCmembraneatasurfacetensionof25mN/m.J.Phys.Chem.B,Vol.108,No.2,2004Marrinketal. ).ThemagnitudeofthepermeationrateisofthesameorderofmagnitudeastheexperimentallymeasuredpermeabilitycoefficientsforpureDPPCvesicles.3.3.6.CrystallinePhase.Uponcoolingofthe256DPPCbilayerpatchbelowatemperatureof283K,formationofacrystallinephasecanbeobserved.InFigure10thisprocessisshowninaseriesofsnapshots.At283Kthebilayerremainsinafluidphasefortensofnanoseconds,beforethetransitiontakesplace.Thetransitionisobservedtobehighlycooperative,takingplacewithinafewnanoseconds.TheareaperheadgroupintheCGcrystallinephaseis0.47nm,withthetailspackedinahexagonallattice.Somepackingdefectsexist,however,ascanbeseenfromthesnapshots.Theabsenceoflipidlateraldiffusivemotionimpliesacrystallinestateratherthanagelphase.Thelipidheadgroupsremainmobile.Experimentally,DPPCformsacrystallinephase(L)belowtemperaturesof273Kwithanareaperheadgroupof0.46nm.IntheexperimenttheDPPCtailsaretiltedwithrespecttothebilayernormalwithatiltanglecloseto35.Althoughthesimulationboxisfullyadjustableallowinganypreferentialpacking,nospontaneoustiltisobserved.Uponheatingofthesystem,meltingisobservedattemperaturesabove300K.Whetheragelphasecanalsobestabilizedasanintermediatephaseisnotyetclear.3.4.LipidBilayers:OtherPhospholipids.otherthanDPPCcanalsobereadillymodeledusingtheCGapproach.Thetaillengthcanbemodifiedbyremovingoraddingadditionalapolarsites.Onthebasisoftheoptimalmappingforalkanes(seesection3.1),alauroyltailismodeledwith3Csites,bothmiristoylandpalmitoylwith4,andstearoylwith5.Cisdoublebondscanbemodeledinasimilarmannerasinalkenes(seesectiononalkane/water)bysettingtheequilibriumangleoftheCGtriplet,whichincludesthedoublebondto120andincreasingthechainstiffnessslightlyto35kJmol.Aphosphatidylethanolamine(PE)headgroupcanbemodeledbytwoparticlesjustasforPC.Tomimicthestronghydrogenbonddonorcapabilitiesandindirectacceptorcapabili-ties(viawatermeditatedhydrogenbonds)oftheethanolaminesitecomparedtothecholinesite,theparticlerepresentingtheethanolaminegroupisofthemorepolarsubtypeda.TheareaperheadgroupobtainedfromCGsimulationsforanumberofcommonphospholipidsisreportedinTable3.Thesimulationswereperformedonpatchesof256lipidswithahydrationlevelof40waters/lipid,closeto,orbeyond,theswellinglimitofmostphospholipids.Thesimulationscoverarangeofphysiologicallyrelevanttemperatures.Zerosurfacetensionconditionsapplytoallsystems.Giventheuncertaintyinexperimentalmeasurements,theCGmodelperformswell,atthesamelevelofaccuracyasatomisticforcefields.CrucialeffectssuchasthecondensationofthePEheadgroupwithrespecttothePCheadgroup(compareDPPC/DPPEat338K),theareaexpansionupontemperatureincrease(DPPCat300338K)oruponchainunsaturation(e.g.,DOPC/DPPCat300K),aswellasthenegligibleeffectofchainelongation(DLPC-DPPC-DSPCat338K)areallreproduced.3.5.NonlamellarPhases.Inthissectionweshowthatthecoarsegrainedmodelcanalsobeappliedtostudynonlamellarphases.Althoughoptimizedandtestedforbilayersystems,thereisnobiasinthemodelfavoringlamellarphases.Thissectiondescribestwosuccessfulapplicationsinwhichlipidsspontane-ouslyaggregateintoeitheraninvertedhexagonalphaseoramicellarphaseinagreementwiththeexperimentallydetermined3.5.1.InertedHexagonal.amine(DOPE)experimentallyformsaninvertedhexagonalphaseattemperaturesabove300K(dependingonhydrationconditions).Themaximumamountofwaterthatcanbetakenupbythehexagonalphaseisabout20watersperlipid.TostudythephasepreferenceoftheCGmodelforDOPE(seeprevioussection),asystemconsistingof1024DOPElipidsrandomlysolvatedinto4224solventsites(16realwaters/lipid)wasgenerated.Thepressurescalingwasperformedcompletelyanisotropically,allowingthedeformationoftheboxshapeandthedevelopmentofhexagonalsymmetry.Thesystemwassimulatedattwodifferenttemperatures:273and318K.Inagreementwiththeexperimentalbehavior,atthelowertem-peratureDOPEaggregatesintoalamellarphase,whereasatthehighertemperatureaninvertedhexagonalphaseformsspontaneously.Repeatedsimulationswithdifferentrandomstartingconditionsshowthesamebehavior.Theaggregationintotheinvertedhexagonalphasecanbedividedintofoursubstages:(i)rapidclusteringofwaterintoanetworkofconnectedinvertedmicelles,(ii)disappearanceofsomeandgrowingofotheraqeuousconnectionsleadingtotheformationoftwoinvertedcylindricalmicelles,(iii)connectionofthesemicelleswiththeirperiodicimageresultingintheformationoftwowaterchannels,and(iv)relaxationofthewaterchannelsintoaperfecthexagonallattice.Thefirststageiscompletewithinafewnanoseconds.Stagestwoandthree,whicharenotalwayswelldistinguished,togethertake100ns.Thefinalrelaxationisaccomplishedwithintensofnanoseconds.Aswiththeformationofbilayers(seethesectiononDPPC),thedynamicsofself-aggregationisverydependentonsystemsize.Spontane-ousaggregationrunsforsystems8timeslargerremaintrapped Figure10.Thermalphasetransitionfromliquid-crystallinetocrystalphase.ADPPCbilayerpatchinthefluidphase(left)wascooledto283K.After20nsthesystemunderwentacooperativephasetransition(middle)andrelaxedintoacrystalphase(right). TABLE3:AreaperLipidforCommonPhospholipidsintheCGModel simulatedarea,experimentalarea,DPPC0.47(283K)0.46(273K),0.48(293K)42,26DPPC0.59(300K)DPPC0.64(323K)0.64(323K)26DPPC0.66(338K)0.64(333K),0.67(338K)43,44DLPC0.60(300K)0.57(293K),0.63(303K)43,44DLPC0.66(338K)0.64(333K),0.71(338K)43,44DSPC0.66(338K)0.65(333K),0.66(338K)43,44DOPC0.67(300K)0.72(303K)26DPPE0.62(338K)0.61(342K)44DOPE0.61(273K)0.65(271K)45DLPE0.55(300K)0.51(308K)26POPE0.59(300K)0.57(303K)45Simulationresultsareaccuratetowithin0.01nm.Experimentalresultsaresubjecttosubstantialinterpretationuncertainties(typicallynotknownwithaprecisionofmorethan0.02nmGelphase.Extrapolatedestimates.CGModelforSemiquantitativeLipidSimulationsJ.Phys.Chem.B,Vol.108,No.2,2004 atsubstagetwo(withmanyinvertedcylindricalmicelles)atleastuptothemicrosecondtimescale.Asimulationofasystemconsistingofeightcopiesoftherelaxedinvertedhexagonalphase(morethan8000lipids)wasstable,however.Thefinalsnapshot(at0.5s)ofthissystem,witheightindependentwaterchannels,isshowninFigure11.Thewaterchannelsareclearlyarrangedinahexagonalpattern.Thelipidshaveadoptedaninvertedgeometrywiththeheadgroupsconfinedintothesmallareaaroundthechannels,leavingamplespaceforthetailsthatresisttightpackingduetothepresenceofadoublebondineachtail.Remarkably,thegeometryadaptedbythewaterchannelsisnotspherical,buthexagonal.Experimentallythedetailsofthepackingarenotknown,butbothasphericalandahexagonalarrangementhavebeenproposed.Thehexago-nalpackingaroundthechannelsalsoshowsupattheborderbetweenthetwomonolayerleaflets,reflectedbythedistributionoftheterminaltailgroup.Thehexagonalspacinginthesimulationsis5.7nm,whichisinagreementwiththeexperi-mentalspacing(estimatedtobe5.75.8nmat318Kwith16waters/lipid)determinedfromX-rayanalysisoftheinvertedhexagonalphaseofDOPE.Amoreextensiveanalysisofthestructureofsimulatedinvertedhexagonalphasesandamoreelaboratecomparisontotheavailableliteraturedatawillbepublishedseparately.3.5.2.Micelle.Removingtheglycerollinkageaswellasoneofthetails,andonesiteoftheremainingtailfromDPPC,thefiveremainingsitesmodeldodecylphosphocholine(DPC).ThemappingoftheDPClipidintheCGmodelisshowninFigure1.Experimentally,DPCformsmicellesatconcentrationsabovethecriticalmicelleconcentration(cmcAsystemcomposedof400DPClipidsrandomlysolvatedby125,000CGwaterswassimulatedforamicrosecond.Theconcentrationoflipidsis0.04M,wellabovethecmc.TheDPClipidsindeedspontaneouslyformsmallmicelles.Themicellarsizedistributiondoesnotconvergewithinamicrosecond,however.Encountersbetweenmicellesarediffusion-limitedandthereforerare.Theexchangerateoflipidsfrommicellesintosolutionisalsotooslowtoallowforarapidequilibration.ThefinaldistributionisshowninFigure12.AfewsmallerDPCclusters(20)stillexist,butmostlipidsareorganizedinsmallsphericalmicelles.Thesizeofthesemicellesrangesfrom40to70.ExperimentallythesizeofDPCmicellesisdeterminedas565(at0.02MAggregationofasmallersystemconsistingof54DPClipidsand5600CGwaters(0.12M)wasalsostudiedtocomparetherateofaggregationtoresultsobtainedwithatomisticsimulationsforasimilarsystem.TheinitialrateofclusteringofsinglelipidsintosmallclustersofafewlipidstakesplacewithinananosecondforboththeCGandtheatomisticsystem.Subse-quentformationoflargerclustersandmicelles,however,isclearlyfasterintheatomisticsimulations.Theformationtimeofthefinalmicelle(containingall54lipids)intheatomisticsimulationwas12ns,whereaswiththeCGmodelittakesbetween50and150nsinfourindependentattempts.Thereasonforthedifferentdynamicsislikelytheabsenseoflong-rangeattractionintheCGmodel.TheradialstructureofthefinalCGmicelleisverysimilartothatoftheatomisticone.InFigure13wecomparetheradialatomdensitydistributionoftheCGmicelletoanatomisticonesimulatedbyTielemanetal.,containing54DPCmolecules.Thedistributionofthevariousatomgroupsinsidethemicelleisverysimilar.Onlythedipindensityinthecenterofthemicelleappearsmorepronouncedintheatomisticsimulation.Themicelleradiusis2nm(tothephosphatepeak)inbothsystems.4.DiscussionIntheintroductionfourcornerstonesforasuccessfulcoarsegrainedmodelweredefined:speed,accuracy,applicability,andversatility.Thecoarsegrainedmodelpresentedinthispaperhasbeenoptimizedtomatchthesecriteriaaswellaspossible.(i)Speed:Aspeedupfactorbetween10and25isobtainedduetotheuseofa50fstimestepinsteadof25fsmaximumforatomisticsimulations.Thereducednumberofinteractionsitesallowssimulationsapproximatelybetween5(forforcefieldwithunifiedmethylenegroups)and10timesfaster(forforcefieldsthathaveexplicithydrogens).Furthermore,theshort-rangenatureoftheinteractionsspeedsupthecalculationsbetween3times(comparedtosimulationsusingtypicalcutoffsof1.4nm)and10times(whenEwaldsummationtypetechniquesareused Figure11.SnapshotofaninvertedhexagonalphaseofDOPEat318K.Thesystemcontainsmorethan8000lipidsandat16waters/lipid.Theethanolaminegroupsarecoloredorange,phosphatidylgroupsyellow,theglycerolbackbonered,thetailsgray.Alightershadeofgrayisusedtodistinguishtheterminalmethylgroup.Wateriscoloredcyan.Thesystemisviewedalongthedirectionofthewaterchannels.Notethehexagonalstructureofthechannelsthemselves. Figure12.DistributionofcoarsegrainedDPCmicellesafter1sofsimulation.Thecholinegroupsaredepictedblue,thephosphateyellow,andthetailsgray.Wateriscoloredcyan.Oneperiodicimageofthesimulationboxisalsodepicted.Theboxedgelengthis25nm.J.Phys.Chem.B,Vol.108,No.2,2004Marrinketal. totreatthelong-rangeinteractions).TheincreaseddynamicsintheCGmodelresultsinanothereffectivespeedupwithafactorof4.Together,thesefactorsresultinaspeedupof4ordersofmagnitudewiththeCGmodelcomparedtocurrentlyusedatomisticsimulationtechniques.AlsocomparedtothesemiquantitativeCGmodelofShelleyetal.thecurrentmodelismuchfasterandcomparabletothespeedobtainedwithDPDtechniques.(ii)Accuracy:TheCGmodelisshowntobeaccurateatleastatasemiquantitativelevelforstructural,elastic,dynamicaswellasthermodynamicpropertiesforarangeoflipidsystems.Structuralpropertiesoflipidsystemssuchastheareaperheadgroupinthelamellarphaseorthehexagonalspacingintheinvertedhexagonalphaseagreewellwiththeavailableexperimentaldata.Comparedtoresultsobtainedwithatomisticsimulations,atomdensitydistributionsareverysimilarinallcasesconsidered.ElasticpropertiescomputedforaDPPCbilayersuchasthebendingmodulus,thelinetensionortheareacompressibilityareofthesameorderastheexperimentalmeasurements.TheabsolutedynamicsintheCGmodelisfastercomparedtorealsystemsorsimulationswithatomisticforcefields.Therelativedynamics,however,appearswellpreserved.Withatimeconversionfactorof4avarietyofdynamicpropertiessuchasself-diffusioninbulkphasesandlipidlateraldiffusion,thepermeationrateofwateracrossamembrane,orthelipidself-aggregationratearemeaningfulatasemiquanti-tativelevel.Thermodynamicallyspeaking,theCGmodelhasencouragingpropertiestoo.Notonlyisthemutualsolubilityofwaterandalkanewellreproduced,butmoreimportantlylipidsaggregateintothecorrectphaseswhetherlamellar,micellarorhexagonal.ForDPPCafreezingtransitionwasevenobservedwhenthetemperaturewasloweredbelowthemainphasetransitiontemperature.(iii)Applicability:Theuseofphysicallymeaningfulparametersmakesthemodeleasytointerpret.Insteadofdealingwithreducedunitsthatneedconversionafterward,itisimmediatelyclearwhatthestateconditionsare.Thelimitednumberofparticletypesandinteractionlevelsprovideforasmallsetofbuildingblocksfromwhichrelatedmoleculesthatareexpectedtoperformatthesamelevelofaccuracyastheexamplesgiveninthispapercanbeeasilyconstructed.(iv)Versatility:Althoughoptimizedmainlyforlipidsystemsinthelamellarphase,theCGmodelhasnobuiltinrestrictionsastothephaseofthesystem.ThelipidsintheCGmodelareveryflexible,freetoadapt(almost)anyconformationatareasonableenergycost.Twoapplicationstononlamellarsystemswereusedtoillustratethisversatility.Recently,theCGmodelhasalsobeensuccessfullyappliedtostudythespontaneousformationofvesicles,andtothefusionmechanismofvesicles.Apartfromphospholipidsystems,theCGmodelcaninprinciplebeappliedtoothertypesofmoleculesaswell.Forinstance,currentlyaCGtopologyisbeingdevelopedforcholesterol.OfcoursetheCGmodelaspresentedinthispaperislimited.Therewillbemanyapplicationsforwhichitisnotwellsuited.FurtheroptimizationsarepossiblebutaswithanyCGmodelitsutilityisinherentinitssimplicity.Anyapplicationforwhichlong-rangeelectrostaticforcesareimportantshouldbeconsid-eredwithcare.Finechemicaldetailisinaccessibleinanycoarsegrainedapproach.Oneshouldbecarefulnottooverinterprettheresultsonaquantitativelevel.TheCGmodelisnotatooltoreplaceatomisticsimulations,butratherthetwoshouldbeusedsidebyside.WiththeCGmodelthelongtime-scaleorlength-scalepropertiesofthesystemofinterestcanbeexplored,whereaswithatomisticmodelsthedetailscanbestudied.ResultsformatomisticsimulationsshouldbeusedasmuchaspossibletojudgethequalityoftheCGforcefieldintheapplicationatReferencesandNotes(1)Marrink,S.J.;Lindahl,E.;Edholm,O.;Mark,A.E.J.Am.Chem.,8638.(2)Shelley,J.C.;Shelley,M.Curr.Opin.ColloidInterfaceSc.,101.(3)MuÈller,M.;Katsov,K.;Schick,M.J.Pol.Sci.B,1441.(4)Smit,B.;Hilbers,P.A.J.;Esselink,K.;Rupert,L.A.M.;VanOs,N.M.;Schlijper,A.G.,624.(5)Smit,B.;Esselink,K.;Hilbers,P.A.J.;VanOs,N.M.;Rupert,L.A.M.;Szleifer,I.,9.(6)Palmer,B.J.;Liu,J.,746.(7)Goetz,R.;Lipowsky,R.J.Chem.Phys.,7397.(8)DenOtter,W.K.;Briels,W.J.J.Chem.Phys.,4712.(9)VonGottberg,F.K.;Smith,K.A.;Hatton,T.A.J.Chem.Phys.,9850.(10)Noguchi,H.;Takasu,M.Phys.Re,1.(11)Groot,R.D.;Madden,T.J.;Tildesley,D.J.J.Chem.Phys.,9739.(12)Shillcock,J.C.;Lipowsky,R.J.Chem.Phys.,5048.(13)Ayton,G.;Voth,G.A.Biophys.J.,3357.(14)Venturoli,M.;Smit,B.Phys.Chem.Commun.,1.(15)Groot,R.D.;Rabone,K.L.Biophys.J.,725.(16)Yamamoto,S.;Maruyama,Y.;Hyodo,S.J.Chem.Phys.,5842.(17)Jury,S.;Bladon,P.;Cates,M.;Krishna,S.;Hagen,M.;Ruddock,N.;Warren,P.Phys.Chem.Chem.Phys.,2051.(18)Shelley,J.C.;Shelley,M.;Reeder,R.;Bandyopadhyay,S.;Klein,M.L.J.Phys.Chem.B,4464.(19)Lindahl,E.;Hess,B.;VanderSpoel,D.J.Mol.Mod.(20)Lide,D.R.CRCHandbookofChemistryandPhysics,72nded.;CRCPress:BocaRaton,FL,1992.(21)Douglass,D.C.;Mccall,D.W.J.Phys.Chem.,1102.(22)Parrinello,M.;Rahman,A.J.Appl.Phys.,7182.(23)Schatzberg,P.J.Phys.Chem.,776.(24)Ben-Naim,A.ationThermodynamics;PlenumPress:NewYork,1987.(25)Copestake,A.;Neilson,G.;Enderby,J.J.Phys.CSolidState,4211.(26)Nagle,J.F.;Tristram-Nagle,S.Biochim.Biophys.Acta(27)Lindahl,E.;Edholm,O.Biophys.J.,426.(28)Marrink,S.J.;Mark,A.E.J.Phys.Chem.B,6122.(29)Safran,S.A.StatisticalThermodynamicsofSurfaces,Interfaces,andMembranes;Addison-Wesley:Reading,MA,1994.(30)Rawicz,W.;Olbrich,K.C.;Mcintosh,T.;Needham,D.;Evans,Biophys.J.,328.(31)AneÂzo,C.;deVries,A.H.;HoÈltje,H.D.;Tieleman,D.P.;Marrink,S.J.J.Phys.Chem.B,9424.(32)Glaser,R.W.;Leikin,S.L.;Chernomordik,L.;Pastushenko,V.F.;Sokirko,A.I.Biochim.Biophys.Acta,275.(33)Leontiadou,H.;Mark,A.E.;Marrink,S.J.Biophys.J.,inpress. Figure13.RadialdistributionfunctionsofaDPCmicelle,fromatomisticsimulations(dashedlines)andCGsimulations(solidlines).CGModelforSemiquantitativeLipidSimulationsJ.Phys.Chem.B,Vol.108,No.2,2004 (34)Chernomordik,L.V.;Kozlov,M.M.;Melikyan,G.;Abidor,I.;Markin,V.;Chizmadzhev,Y.Biochim.Biophys.Acta,643.(35)Zhelev,D.V.;Needham,D.Biochim.Biophys.Acta(36)Evans,E.;Heinrich,V.Compt.Rend.,265.(37)Winterhalter,M.Curr.Opin.ColloidInterfaceSc.,250.(38)Sheats,J.R.;McConnell,H.M.Proc.Natl.Acad.Sci.U.S.A.,4661.(39)Kuo,A.L.;Wade,C.G.,2300.(40)Carruthers,A.;Melchior,D.L.,5797.(41)Andrasko,J.;Forsen,S.Biochem.Biophys.Res.Commun.,813.(42)Nagle,J.F.;Wiener,M.C.Biochim.Biophys.Acta,1.(43)BalgavyÂ,P.;DubnickovaÂ,M.;Kucerka,N.;Kiselev,M.A.;Yaradaikin,S.P.;UhrõÂkovaÂ,D.Biochim.Biophys.Acta,40.(44)Petrache,H.I.;Dodd,S.;Brown,M.Biophys.J.,3172.(45)Rand,R.P.;Parsegian,V.A.Biochim.Biophys.Acta(46)Rand,R.P.;Fuller,N.Biophys.J.,2127.(47)Hamm,M.;Kozlov,M.M.Eur.Phys.J.B,519.(48)Lauterwein,J.;BoÈsch,C.;Brown,L.R.;WuÈthrich,K.Biophys.Acta,244.(49)Marrink,S.J.;Tieleman,D.P.;Mark,A.E.J.Phys.Chem.B,12165.(50)Tieleman,D.P.;VanderSpoel,D.;Berendsen,H.J.C.J.Phys.Chem.B,6380.(51)Marrink,S.J.;Mark,A.E.J.Am.Chem.Soc.,inpress.(52)Marrink,S.J.;Mark,A.E.J.Am.Chem.Soc.,11144.J.Phys.Chem.B,Vol.108,No.2,2004Marrinketal.