Table Of ContentTHEJOURNALOFCHEMICALPHYSICS134,055109(2011)
On the proper calculation of electrostatic interactions in solid-supported
bilayer systems
In-ChulYeha) andAndersWallqvist
BiotechnologyHighPerformanceComputingSoftwareApplicationsInstitute,Telemedicine
andAdvancedTechnologyResearchCenter,U.S.ArmyMedicalResearchandMaterielCommand,
FortDetrick,Maryland21702,USA
(Received21September2010;accepted6January2010;publishedonline4February2011)
Modeling systems that are not inherently isotropic, e.g., extended bilayers, using molecular simu-
lation techniques poses a potential problem. Since these methods rely on a finite number of atoms
andmoleculestodescribethesystem,periodicboundaryconditionsareimplementedtoavoidedge
effects and capture long-range electrostatic interactions. Systems consisting of a solvated bilayer
adsorbed on a solid surface and exposed to an air/vacuum interface occur in many experimental
settingsandpresentsomeuniquechallengesinthisrespect.Here,weinvestigatedtheeffectsofim-
plementing different electrostatic boundary conditions on the structural and electrostatic properties
ofaquartz/water/vacuuminterfaceandasimilarquartz-supportedhydratedlipidbilayerexposedto
vacuum. Since these interfacial systems have a net polarization, implementing the standard Ewald
summationwiththeconductingboundaryconditionfortheelectrostaticlong-rangeinteractionsin-
troducedanartificialperiodicityintheout-of-planedimension.Inparticular,abnormalorientational
polarizations of water were observed with the conducting boundary condition. Implementing the
Ewald summation technique with the planar vacuum boundary condition and calculating electro-
static properties compatible with the implemented electrostatic boundary condition removed these
inconsistencies.Thisformulationisgenerallyapplicabletosimilarinterfacialsystemsinbulksolu-
tion.©2011AmericanInstituteofPhysics.[doi:10.1063/1.3548836]
I. INTRODUCTION The Ewald summation formula contains a surface
term that depends on the dielectric constant of the medium
As many chemical, biological, and physiological pro-
surrounding the periodic images.2,5,14,15 This surface term
cessessensitivelydependonlong-rangeelectrostaticforces,it
vanishes for the conducting boundary condition where the
isessentialthatthemoleculardescriptionsusedtoprobethese
dielectric constant of the surrounding medium is infinite.2,5
systems correctly account for these forces.1 Periodic bound-
However, for vacuum boundary conditions where the sur-
ary conditions (PBCs) are the natural choice in computer
rounding medium is characterized by a dielectric constant
simulations toavoid edge effects due tothefinite number of
of 1, the surface term is expressed as a function of the total
atoms and molecules located in a unit cell.2 Treatments of
dipole moment or polarization of the system.2,5 The surface
electrostatic interactions under PBCs can strongly influence
termhasaslightlydifferentfunctionalformdependingonthe
calculated properties and, hence, the physical and biological
orderormodeofsummation.14,15Itplaysanimportantrolein
interpretationofthemodeledsystems.Anumberofdifferent
modelingmolecularsystemsthatarenotinherentlyisotropic,
methods3,4 have been implemented for calculating electro-
e.g.,extendedbilayers,orpolarizedinterfacialsystems.16For
staticinteractionsinsystemsdescribedbyPBCs.Themethod
systems with planar interfaces, the planar vacuum boundary
of choice, the Ewald summation method,2,3,5,6 combined
condition where periodic images are surrounded by planar
withtheparticlemeshEwaldtechnique,7providesanefficient
vacuum boundaries would be a natural choice. For this
mean to calculate electrostatic interactions for molecular
boundary condition, the Ewald summation is carried out in
systems subject to PBCs and has been broadly implemented
the planar manner, i.e., summing in planar directions first
in biomolecular simulation studies.8 Solid-supported lipid
and then progressing in the out-of-the plane direction. The
bilayer films exposed to an air/vacuum interface have been
resulting surface term depends only on the total dipole mo-
studiedextensivelyexperimentally9andcomputationally.10,11
mentintheout-of-theplanedirection.14YehandBerkowitz16
However, the implementation of the Ewald sum in such
demonstrated that simulations with the surface term cor-
nonisotropic systems presents some unique problems. An
responding to the planar vacuum boundary condition in
improper consideration of electrostatic boundary conditions
conjunctionwithanextravacuumspacealongtheout-of-the
forsuchsystemscanresultinnonphysicalartifacts.12,13Here,
plane direction reproduce results obtained by the rigorous,
we present an implementation of the electrostatic boundary
but computationally expensive, two-dimensional Ewald
condition that includes a straightforward correction term to
summation technique17 almost exactly. This approximation
removetheseartifacts.
holds as long as the length of the extra vacuum space along
a)Author to whom correspondence should be addressed. Electronic mail: theout-of-theplanedirectionislongerthantheperiodicbox
icy@bioanalysis.org. lengths in planar directions.16 Subsequently, simulations of
0021-9606/2011/134(5)/055109/9/$30.00 134,055109-1 ©2011AmericanInstituteofPhysics
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it
does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
SEP 2010 2. REPORT TYPE 00-00-2010 to 00-00-2010
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
On the proper calculation of electrostatic interactions in solid-supported
5b. GRANT NUMBER
bilayer systems
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
U.S. Army Medical Research and Materiel Command,Biotechnology REPORT NUMBER
High Performance Computing Software Applications
Institute,Telemedicine and Advanced Technology Research Center,Fort
Detrick,MD,21702
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON
a. REPORT b. ABSTRACT c. THIS PAGE Same as 9
unclassified unclassified unclassified Report (SAR)
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
055109-2 I.YehandA.Wallqvist J.Chem.Phys.134,055109(2011)
manyinterfacialsystemshavesuccessfullyimplementedthis molecules.Theforcefieldforaquartz(011)crystalwastaken
surface term for the planar vacuum boundary condition.18,19 fromtheonedevelopedbyLopesetal.26Aprimitiveunitcell
However, the vast majority of molecular dynamics (MD) ofquartz(011)of∼15-Åthicknesswasreplicatedintwodi-
simulations are still being performed with the conducting mensions to construct a quartz (011) crystal with a surface
boundary conditionregardless whether thesimulatedsystem area of 36.5 × 37.0 Å2, as in the previous study by Lopes
hasanetpolarizationornot.Thus,itisimportanttohighlight et al.26 One side of the constructed crystal was covered
theconsequenceofusingconductingboundaryconditionsin with hydrophilic silanols (Si–OH), and silicon atoms on the
simulationsofsystemswithnonzeronetpolarization. other side were saturated with hydrogens (Si–H). Crystal
It is important when simulations are performed using a atoms, except for the O–H groups of silanols, were held
specificelectrostaticboundaryconditionthatthecalculations fixed to maintain the (011) crystal geometry. Simulations of
ofelectrostaticpropertiesareconsistentwiththeimplemented the quartz/water interface and the quartz-supported lipid bi-
boundary condition. Sachs et al.20 and Yeh and Hummer21 layer embedded in vacuum were performed with fixed sim-
introduced an accurate formula to calculate the electrostatic ulation box sizes using 3D PBCs. However, the box length
potential profile from simulations performed with the con- (L ) in the z-direction normal to the interface was chosen to
z
ducting boundary condition. This formula includes a pre- be large enough for the system to form interfaces with the
viously overlooked integration constant, which is necessary vacuum. This allowed the systems to adjust to an optimum
to enforce the three-dimensional (3D) PBC. Interestingly, as density without constant pressure simulations. The penetra-
will be shown below, this integration constant has the same tion of water molecules into the vacuum phase was negli-
functional form as the surface term for the planar vacuum gible under all systems and conditions studied in this work.
boundaryconditionbutwithanoppositesign.Asaresult,the Short-rangeinteractionsoutsidea10-Åcutoffweretruncated.
traditional formula for calculating the electrostatic potential Long-rangeelectrostaticinteractionswerecalculatedwiththe
without the integration constant was found tobe appropriate particle mesh Ewald method with and without the surface
for simulations performed with the planar vacuum boundary term for the planar vacuum boundary condition,16 referred
condition. to as EW3DC and EW3D, respectively. The EW3D method
Earlierconsiderationsoftheplanarvacuumboundarydid isequivalenttoimplementingtheEwaldsummationwiththe
notspecificallyaddressthesystemsconsideredhere,inpartic- conducting boundary condition. The contribution to the po-
ular,thosewithintrinsicallyisotropiccomponentsmixedwith tentialenergy(U )fromtheEW3DCcorrectiontermisgiven
c
a strongly polarized local component. These systems with bythefollowingsurfacetermfortheplanarvacuumboundary
nonzeronetpolarizationareencounteredevenintheabsence condition:
of charged surfaces or external electric fields, e.g., in net- 1 2π 1
neutral asymmetric lipid bilayers or quartz crystal microbal- U = M2 = M2, (1)
c 4πε V z 2ε V z
ancemeasurementsperformedonpolarizedcrystals.22Inthis 0 0
report,weexaminedthestructuralandelectrostaticproperties whereε ,V,andM arethevacuumpermittivity,thevolume
0 z
ofinterfacialsystemswithnonzeronetpolarizationcalculated of the unit simulation cell, and the z-component of the total
withdifferentelectrostaticboundaryconditions.Weobserved dipolemomentoftheunitsimulationcell.Thecontributionto
abnormal orientational polarizations of water in simulations theforce(F )actingonatomiwithanelectrostaticchargeof
z,i
ofaquartz/waterinterfaceandasimilarquartz-supportedhy- q fromthistermisasfollows:
i
drated dimyristoylphosphatidylcholine (DMPC) bilayer em-
∂U q
bsuemddmedatiinonvawcuituhmth,peecrofonrdmuecdtinwghibloeuunsdinagrythceosntdaintidoanr.dIEmwpalled- Fz,i =−∂zic =−ε0Vi Mz. (2)
menting the Ewald summation technique with the planar
We implemented the EW3DC correction term in NAMD by
vacuum boundary condition removed these abnormal polar-
adding forces due to the following effective electric field
izations of water. Here, we give the correct formulas for
(EEW3DC)inthez-direction16:
obtaining electrostatic properties consistent with different
boundary conditions. We also confirmed that the modeled M
EEW3DC =− z . (3)
quartz-supportedhydratedlipidbilayerreproducedtheessen- ε V
0
tial features of the lipid bilayer in bulk water. Finally, we
This electric field can also be related to the z-component of
showedthattheplanarvacuumboundaryconditioncanbeap-
polarization(P )asfollows:
pliedtobulkphasesimulationsofaquartz/waterinterfacein z
ordertoremoveartifactsintroducedbytheconductingbound-
P
arycondition. EEW3DC =− z. (4)
ε
0
BondsinvolvinghydrogenswereconstrainedwiththeSHAKE
II. METHODS algorithm.27 A time step of 2 fs was used for the time
integration. The temperature was maintained at 310 K
A. Simulations
using Langevin dynamics with a damping coefficient of
MDsimulationswereperformedwiththeNAMDMDsim- 10 ps−1.28 Coordinates were saved at every 1 ps for further
ulation program23 and CHARMM all atom force field24 for analysis.Errorsreportedinthisworkwerebasedonthe95%
DMPC.TheTIP3Pwatermodel25wasusedtodescribewater confidenceintervalestimatedas1.96timesthestandarderror.
055109-3 Electrostaticsinsolid-supportedbilayer J.Chem.Phys.134,055109(2011)
FIG. 2. A snapshot from the simulation of a hydrated dimyristoylphos-
phatidylcholine(DMPC)lipidbilayerdepositedonaquartzcrystalsurface.
TheblueboxrepresentstheunitsimulationcellwithLzof200Å.Thesim-
ulationcellisrepeatedperiodicallyinallthreedirections.
prismwiththe CHARMM-GUIMembraneBuilder.29 Thesur-
face area of the rectangular base of the simulation cell was
chosen to be same as the area of the quartz (011) surface
so that the area per lipid headgroup was ∼60 Å2, close to
the experimental value for DMPC.30 The prepared hydrated
DMPC bilayer was equilibrated under NPAT (constant nor-
mal pressure and constant lateral surface area of the bilayer
and constant temperature) condition with the CHARMM MD
simulationprogram.31 AfterabriefrununderthesameNPAT
FIG.1. Asnapshotfromthesimulationofaquartz/waterinterfaceembedded ensemble with NAMD, the lateral periodic boundaries were
invacuum.Theblueboxrepresentstheunitsimulationcellwithaboxlength graduallychangedtothoseofthequartz(011)surface.There-
inthez-direction(Lz)of100Å.Thesimulationcellisrepeatedperiodically sultinghydratedlipidbilayerwasbroughtclosertothequartz
inallthreedirections.Waterisincontactwiththehydrophilicsideofthe
(011)surface,andthethicknessofthewaterlayerthatinter-
quartz(011).Emptyspace(vacuum)isusedtoseparatewaterandtheother
faceswiththemembraneandvacuumwasdoubled.Figure2
sideofthequartzcrystalinthez-direction.
showsthesnapshotfromMDsimulationwiththissetup.Both
theEW3DandEW3DCmethodswereusedtocalculateelec-
Standarderrorswereestimatedfromananalysisoffiveequal trostaticinteractions.ThesameL of200Åwasusedinboth
z
timeblocksfromthetrajectoryofeachproductionrun. methods.Figure2showsthatthelengthofthevacuumregion
with L of 200 Å was clearly larger than lateral box lengths.
z
Each simulation lasted for 200 ns, of which the initial 70 ns
B. Quartz/waterinterface
wasdiscardedduetoequilibrationconsiderations.Toinvesti-
Aquartz/waterinterfacewaspreparedbyadding985wa- gatepossiblefinitesystem-sizeeffects,weperformedanaddi-
termoleculestothehydrophilicsideofthecrystal.Foursim- tional100-nssimulationwiththeEW3DCmethod,wherethe
ulationswithdifferentL valuesof100,150,200,and300Å systemsizewasquadrupledtohave176lipidsbyduplicating
z
were performed with the EW3D method. The smallest box the simulation system along each of two lateral dimensions.
length Lz of 100 Å was used for a simulation with the The value for Lz was increased to 250 Å to ensure that the
EW3DC method. A snapshot of the prepared quartz/water length of the empty space in z-direction is greater than the
systemforaL of100ÅisshowninFig.1.Awater/vacuum increasedlateralboxsizes.Similarly,simulationofaDMPC
z
interface was formed due to the extra empty space between bilayer in bulk water under the NPAT ensemble with the in-
water and the other sideof the crystal inthe simulation cell. creasedsystemsizewasalsoperformedfor100nstocompare
Figure1showsthatthelengthoftheextraemptyspaceinthe propertiesofthelipidbilayerinbulksolutionandnearapo-
simulation cell with the smallest box length is clearly larger larizedquartzinterface.Theinitial30nsofthesesimulations
thanlateralboxlengths,whichisrequiredforthesuccessful wasdiscardedasequilibration.
applicationoftheEW3DCmethodtoeffectivelyremovethe
periodicity in the z-direction.16 Each simulation lasted for
20ns,ofwhichtheinitial1nswastreatedasanequilibration D. Quartz/waterinterfaceinbulksolution
periodanddiscarded.
We prepared a quartz/water interface in bulk solution
withoutvacuumboundariesbyadding1400watermolecules
to both the hydrophilic and hydrophobic sides of the quartz
C. Quartz-supportedlipidbilayer
(011) crystal. In addition to the O–H groups on the hy-
A hydrated lipid bilayer composed of 44 DMPC drophilic side, hydrogen atoms on the hydrophobic side of
molecules (22 molecules on each leaflet) and 1400 water thecrystalincontactwithwaterwereallowedtomove.Sim-
moleculeswaspreparedinthesimulationcellofarectangular ulationswereperformedfor20nsunderthe NPAT ensemble.
055109-4 I.YehandA.Wallqvist J.Chem.Phys.134,055109(2011)
(cid:2)
Tthheecionnitdiaulc1tinngswboausnddiasrcyarcdoenddaistioenquwiliitbhratthioenE.WWe3Dusmedebthoothd E(z)= −zLz/2ρεq(z(cid:3))dz(cid:3) + εPz, (11)
0 0
andtheplanarvacuumboundaryconditionwiththeEW3DC
method to highlight the physical effects introduced by the (cid:3) (cid:4) (cid:5)
choiceofboundarycondition. 1 z P L
φ(z)=− (z−z(cid:3))ρ (z(cid:3))dz(cid:3)− z z+ z .
ε0 −Lz/2 q ε0 2
E. Electrostaticpropertyprofilesacrosstheinterface (12)
Distributions of the electrostatic potential and electric Similar expressions were recently obtained by Gurtovenko
field along the z-direction were calculated from the charge and Vattulainen.32 In their study, they examined the elec-
densitydistributionρq(z)obtainedfromthesimulations.Ex- trostatic potential calculated for asymmetric lipid bilayers
pressionsforcalculatingtheelectrostaticpropertiesweredif- with nonzero net polarization using the conducting bound-
ferentdependingontheelectrostaticboundaryconditionused ary condition.32 For this case, they suggested that calcu-
in the simulations. First, let us consider the case when the lations of the electrostatic potential should be done using
conducting boundary condition wasapplied asintheEW3D formulasappropriatefortheplanarvacuumboundarycondi-
method.StartingfromPoisson’sequation, tion.However,Eqs.(11)and(12)shouldbeusedtocalculate
d2φ(z) ρ (z) electrostaticpropertiesfromsimulationsperformedusingthe
=− q , (5) EW3Dmethodundertheconductingboundarycondition,re-
dz2 ε
0 gardlessofwhetherthesystemispolarizedornot,aswewill
whereφ(z)istheelectrostaticpotentialatpositionz.Anelec- show in Sec. III. However, when using the EW3DC(cid:6)method
tricfield[E(z)]wascalculatedbyintegratingPoisson’sequa- undertheplanarvacuumboundarycondition,the P ε term
z 0
tiononce,asfollows: inEq.(11)iscancelledoutateverytimestepbytheeffective
(cid:2)
E(z)=−dφ(z) = −zLz/2ρq(z(cid:3))dz(cid:3) +C , (6) peloescittreicsifigenl.dTEhEeWre3fDoCreo,fthtehefoslalmoweinmgaegxnpitruedsesiobnustwfoirththteheeloepc--
dz ε0 1 tricfieldandelectrostaticpotentialshouldbeusedinstead:
where C is an integration constant. Sachs et al.20 and Yeh (cid:2)
1 z ρ (z(cid:3))dz(cid:3)
and Hummer21 showed that this integration constant is not E(z)= −Lz/2 q , (13)
ε
necessarily zero but can be determined by enforcing PBCs. 0
Integrating the above equation once more, we obtained the
(cid:3)
followingexpressionfortheelectrostaticpotential:
1 z
(cid:3) φ(z)=− (z−z(cid:3))ρ (z(cid:3))dz(cid:3). (14)
1 z ε q
φ(z)=−ε (z−z(cid:3))ρq(z(cid:3))dz(cid:3) 0 −Lz/2
0 (cid:4)−Lz/2 (cid:5) Equations(13)and(14)arethusappropriatefortheEW3DC
−C z+ Lz +C . (7) methodusingtheplanarvacuumboundarycondition.
1 2
2
Ifwechoosez=−L /2asareferencepointandsetφ(−L /2)
z z III. RESULTSANDDISCUSSION
= 0, the integration constant C becomes 0. Imposing PBCs
2
bysettingφ(−L /2)=0=φ(−L /2+L )=φ(L /2)yielded A. Quartz/waterinterface
z z z z
thefollowingexpressionforC :
1 1.Structuralproperties
(cid:3) (cid:4) (cid:5)
C =− 1 Lz/2 Lz −z(cid:3) ρ (z(cid:3))dz(cid:3). (8) Figure 3(a) shows a comparison of mass density distri-
1 ε0Lz −Lz/2 2 q butions of water along the z-direction using different elec-
trostaticboundaryconditionsandboxlengths.Overall,these
Equation(8)wasexpandedtogivethefollowing:
distributionswereinsensitivetochangesintheboundarycon-
(cid:3) (cid:3)
C1 =−21ε0 −LLzz//22ρq(z(cid:3))dz(cid:3)+ ε01Lz −LLzz//22z(cid:3)ρq(z(cid:3))dz(cid:3). dtfhaicteieownaasttezorr=dbe0onxsÅitlyewnaagtst∼hcs1lo0usseÅedtoawitnhaeythbferuolsmkimwtuhaleatetqiroudnaersnt.zsI/iwntyaatolelfrc0ian.9stee9rs3-,
(9) g/cm3 at 310 K, demonstrating that a well-defined and sta-
blequartz/waterinterfacewasformedwiththeadditionofthe
The first term vanishes because the net charge of the system
water/vacuum interface. However, simulations using smaller
iszero.Then,
(cid:3) (cid:3) boxlengthswiththeEW3Dmethodresultedinslightlylower
C = 1 Lz/2 z(cid:3)ρ (z(cid:3))dz(cid:3) = 1 Lz/2 z(cid:3)ρ (z(cid:3))Adz(cid:3) water densities (see Table I). Small differences in water
1 ε0Lz −Lz/2 q ε0ALz −Lz/2 q density also existed near the water/vacuum interface. Lopes
et al.26 obtained similar water density distributions near the
M P
= z = z. (10) quartz/waterinterfacefromsimulationsofwaterenclosedbe-
ε V ε
0 0 (cid:6) tween quartz surfaces. The distributions are also similar to
BysubstitutingC = P ε intoEqs.(6)and(7),weobtained those obtained by Wander and Clark33 using different force
1 z 0
thefollowingexpressions: fieldparametersrepresentingthequartz/waterinterface.
055109-5 Electrostaticsinsolid-supportedbilayer J.Chem.Phys.134,055109(2011)
(a) withtheEW3DCmethodbutwaslargerthanzeroforsimula-
tions with the EW3D method. These abnormal orientational
1.5 polarizations of water with the EW3D method were more
3m) pronounced at the water/vacuum interface, where the water
g/c1.0 density was low. The results shown in Fig. 3(b) and Table I
y ( show that even though the difference between water ori-
sit
n entation distributions decreased with increasing L values
De 0.5 z
in the EW3D simulations, noticeable differences persisted
even for large L values. This indicates that there were
z
0.0
0 10 20 30 significant differences in electrostatic forces experienced by
z (Å) watermoleculesinsimulationswiththeEW3DandEW3DC
(b) protocols.
EW3D L = 100 Å
0.6 Z
EW3D L = 150 Å
Z
>Z 0.4 EW3D LZ = 200 Å 2.Electrostaticproperties
θdipole, 0.2 EEWW33DD CL ZL = =3 0100 0Å Å termTooleicnuvleesstiignasteimthuelautinounssuwalitohrtihenetEatWion3aDlmoredtehroindg,woef cwaal--
os Z
<c culated electrostatic properties across the quartz/water inter-
0.0
face.Figures4(a)and4(b),showprofilesoftheelectrostatic
potential and electric field calculated by applying Eqs. (13)
-0.2
0 10 20 30 and (14), compatible with the planar vacuum boundary con-
z (Å)
FIG.3. Structuralpropertiesofthequartz/waterinterface.(a)Distributions
(a)
ofmassdensityofwaterunderdifferentelectrostaticboundaryconditions. 20
(b)Averagevaluesofcosineoftheanglebetweenthewaterdipoleandthe EW3D LZ = 100 Å
z-axis((cid:4)cosθdipole,z(cid:5))acrosstheinterface. z), V1105 EEWW33DD LLZZ == 125000 ÅÅ
(φ 5 EW3D LZ = 300 Å
Theaveragecosineoftheanglebetweenthewaterdipole 0 EW3DC LZ = 100 Å
-20 -10 0 10 20 30 40
and the z-axis ((cid:4)cos θdipole,z(cid:5)) is a good indicator of orienta- z (Å)
(b)
tionalpolarizationsofwaterandissensitivetothechoiceof
0.1
electrostaticboundaryconditionsusedinthesimulations.13,34
Å 0.0
Figure3(b)shows(cid:4)cosθ (cid:5)ofwaterasafunctionofthe V/
distancefromthequartzsduiproflea,zce.Underisotropicconditions, Ez(), --00..21
as in bulk water without any external electric field, (cid:4)cos
-0.3
θ (cid:5) is expected to be zero. In previous simulations -20 -10 0 10 20 30 40
dipole,z z (Å)
of water in contact with solid surfaces using the EW3DC (c)
15
method,thewaterphasebecameisotropicroughly10Åaway
forfobmultkhleikseuwrfaacteer.1d8e,3n5sTityhe(1v0alÅue<ofz(cid:4)<co1s5θÅdip)owle,az(cid:5)sactlotsheetroegzieorno z), V10
(φ 5
0
-20 -10 0 10 20 30 40
TABLE I. Values of (cid:4)cos θdipole,z(cid:5) and the average mass density at the z (Å)
bulklikeregion(10Å<z<15Å)ofwaterinthequartz/waterinterface (d)
0.2
as a function of different electrostatic boundary conditions and box size
variationsinthez-direction(Lz). V/Å 0.1
z), 0.0
Electrostatic E(-0.1
boundary
-0.2
condition Lz(Å) (cid:4)cosθdipole,z(cid:5) Density(g/cm3) -20 -10 0 10 20 30 40
z (Å)
EW3D 100 0.0506(0.0001) 0.9909(0.0009)
EW3D 150 0.0301(0.0003) 0.9949(0.0008) FIG.4. Electrostaticpropertiesofthequartz/waterinterface.(a)and(b)Pro-
EW3D 200 0.0214(0.0002) 0.9945(0.0005) filesofelectrostaticpotential[φ(z)]andelectricfield[E(z)],respectively,cal-
culatedbyapplyingEqs.(13)and(14)suitablefortheplanarvacuumbound-
EW3D 300 0.0135(0.0002) 0.9967(0.0009)
ary condition to simulation data generated under the conducting periodic
EW3DC 100 0.0000(0.0003) 0.9972(0.0005)
boundarycondition(EW3D).(c)and(d)Thesameprofilesarecalculatedby
aEW3DandEW3DCrefertoEwaldsummationmethodswithoutandwiththecorrec- applyingEqs.(11)and(12),i.e.,applyingtheconsistentconductingperiodic
boundaryconditiontoboththesimulation(EW3D)andthecalculationof
tiontermfortheplanarvacuumboundarycondition,respectively.Errorsinparentheses
werebasedonthe95%confidenceintervalestimatedas1.96timesthestandarderror. electrostaticproperties.Allresultsfromsimulationsperformedandanalyzed
Standarderrorswereestimatedfromananalysisoffiveequaltimeblocksfromthe usingtheplanarvacuumboundaryconditionwiththeperiodicboundarycon-
trajectoryofeachproductionrun.(cid:4)cosθdipole,z(cid:5)valuesaretheaveragevaluesofcosine ditioncorrection(EW3DC)areshownassolidlinesandwereobtainedby
oftheanglebetweenthewaterdipoleandthez-axis. usingEqs.(13)and(14).
055109-6 I.YehandA.Wallqvist J.Chem.Phys.134,055109(2011)
dition. In this case, the electrostatic potential increased with (a)
increasing z in the bulklike water region, resulting in nega-
tive electric field values, which were not consistent with the 35.6
positive(cid:4)cosθdipole,z(cid:5)valuesshowninFig.3(b).Figures4(c) Å) 35.4
and4(d),showtheresultsobtainedbyapplyingformulasap- (P
P35.2
d
propriatetotheelectrostaticboundaryconditionsusedinthe
35.0
simulations, namely, Eqs. (11) and (12) for simulations us-
ingtheconductingboundarycondition(EW3D)andEqs.(13) 34.8
0 50 100 150 200
and (14) for the simulation using the planar vacuum bound- Time (ns)
ary condition (EW3DC). The electrostatic potential profile (b)
obtained with the EW3D method was nearly flat in the re- 1.5
gionoccupiedbythebulklikewater,resultinginnearzerobut 3m) Water
c Water
slightlypositiveelectricfieldvalues,consistentwiththepos- g/1.0 DMPC
itive(cid:4)cosθdipole,z(cid:5)valuesshowninFig.3(b).However,these sity (
relativelyweakelectrostaticfieldswerenotstrongenoughto en 0.5
D
alter the density distributions of water molecules embedded
0.0
on the hydrogen-bonded network of the liquid water phase. 0 20 40 60 80
In contrast, the average electric field in the region with the z (Å)
(c)
bulklikewaterdensitywiththeEW3DCmethodwascloseto 1.0
zero, also consistent with the (cid:4)cos θ (cid:5) values shown in EW3D
dipole,z
> 0.5 EW3DC
Fig.3(b)andTableI.Thisalsodemonstratesthatthesmallest Z
boxsizeof100Åwaslargeenoughtoobtaincorrectelectro- θdipole, 0.0
static properties when we use the proper correction term for s
o
theplanarvacuumboundarycondition,asimplementedinthe <c -0.5
EW3DCmethod.
-1.0
Gurtovenko and Vattulainen32 claimed that using 0 20 40 60 80
Eq. (12), derived for conducting boundary conditions, for z (Å)
simulation systems with nonzero net dipole moments is not
FIG. 5. Structural properties of the quartz-supported DMPC lipid bilayer.
appropriate even if the simulations are performed with the (a)dPP,theaveragedistancebetweenphosphateatomsinupperandlower
conductingboundarycondition.Instead,theysuggestedusing leaflets,asafunctionoftime.Circlesandsquaresrepresent10-nsblockav-
Eq. (14), which was derived for the planar vacuum bound- eragesforsimulationswiththeEW3DandEW3DCmethods,respectively.
Dashed and solid lines represent cumulative averages for the EW3D and
ary condition. However, our results in Figs. 3 and 4 clearly
EW3DCmethods,respectively.(b)Distributionsofmassdensitycontributed
show that self-consistent electrostatic properties can be ob- bytheDMPClipidbilayerandwatermoleculescalculatedfromsimulations
tained only if we use the formula corresponding to the elec- withtheEW3D(dashedlines)andEW3DC(solidlines)methods.(c)Com-
trostaticboundaryconditionusedinthesimulation. parisonof(cid:4)cosθdipole,z(cid:5)valuesobtainedfromsimulationswiththeEW3D
(dashedlines)andEW3DC(solidlines)methods.
B. HydratedDMPCbilayeronthequartzcrystal bulkphase,thewaterdipoleinitiallypointstowardthemem-
brane at the interfaces (8.4 Å < z < 19.3 Å or 50.2 Å < z
1.Structuralproperties < 61.9 Å). For water molecules penetrating further into the
Figure 5(a) shows the thickness of the DMPC bilayer membrane (19.3 Å < z < 50.2 Å), their dipole points away
represented by d , the average distance between phosphate fromthemembranecenter.Asobservedforwatersinthepure
PP
atomsinupperandlowerleaflets,asafunctionoftime.Both quartz/watersystem,watermoleculesinthesimulationswith
10-nsblockandcumulativeaveragesofd showthatthesim- theEW3Dmethodweremorepolarizedcomparedwiththose
PP
ulatedsystemsrelaxedwithin70ns,ourchoiceoftheequili- withtheEW3DCmethodeveninthebulklikeregion.
brationtime.Figure5(b)showsthecontributionstothemass
densitydistributionsbyDMPCandwateracrossthehydrated
2.Electrostaticproperties
DMPC bilayer deposited on the quartz crystal. Density dis-
tributionscalculatedfromsimulationsusingeithertheEW3D Figure6showselectrostaticpropertiescalculatedforthe
orEW3DCmethodswerequitesimilar.Densitydistributions hydrated DMPC bilayer supported on the quartz crystal. As
ofwaterbetweenthequartzsurfaceandtheDMPCbilayerat shown in Fig. 6(a), the electrostatic potential dropped no-
0Å< z<25Åwereinfluencedbyinteractionswithboththe ticeably across the membrane region and the vacuum phase
quartzandtheDMPCbilayer.Ontheotherhand,thethicker with the EW3D method, resulting in a larger electric field
water layer between the DMPC bilayer and vacuum phase [as shown in Fig. 6(b)] compared with those obtained with
contained a significant bulklike region between z = 65 and the EW3DC method. At a first glance, these potential drops
75Å. across the membrane region and vacuum region may seem
Figure 5(c) shows the values of (cid:4)cos θ (cid:5) as a func- physicallyunreasonable.However,thisisarealconsequence
dipole,z
tionofthedistancefromthequartzsurfaceanddemonstrates ofhavingPBCsundertheconductingboundarycondition.An
thatwhenwatermoleculesapproachthelipidbilayerfromthe increasedpotentialinonepartofthesystem(apolarizedcrys-
055109-7 Electrostaticsinsolid-supportedbilayer J.Chem.Phys.134,055109(2011)
15 (a) 3m) 2.0 (a)
z(), Vφ105 EEWW33DDC Density (g/c 0101....0055-3W0ate-2r0 -10DM0PC10 20Wat3e0r
z (Å)
0 (b)
0 50 100 150 200 > 1.0
Z
(b) z (Å) θdipole, 00..05
0.2 s -0.5
o
<c -1.0
Å 0.1 -30 -20 -10 0 10 20 30
z), V/ 0.0 (c) z (Å)
E(-0.1 0.8
V 0.6
-0.2 0 50 z (Å10)0 150 200 z(), φ 000...024
-30 -20 -10 0 10 20 30
FIG.6. Electrostaticpropertiesofthequartz-supportedDMPClipidbilayer. z (Å)
(a) and (b) Profiles of φ(z) and E(z) with the EW3D (dashed lines) and (d)
EW3DC(solidlines)methods. Å 0.1
V/
z), 0.0
E( -0.1
talinthiscase)hastobecompensatedforelsewheretosatisfy
-30 -20 -10 0 10 20 30
theimposedPBCs. z (Å)
Theaverageelectricfieldneartheregionoccupiedbythe (e)
bulklikewaterbetweenDMPCbilayerandthewater/vacuum 0.2
interface between z = 65 and 75 Å in simulations with the CD
S 0.1
EW3D method had a slight positive value of 0.95 (±0.15) -
0.0
mV/Å. With the polarization P in the z-direction estimated 2 4 6 8 10 12 14
z
fromthe(cid:4)cosθ (cid:5)valueinthebulklikeregionofwaterin Carbon atom number
dipole,z
Fig. 5(b), we can estimate the dielectric constant ε of water FIG.7. ComparisonofDMPClipidbilayersinsolution(dashedlines)and
bythefollowingrelation: withthequartz/waterinterface(solidlines)withalargersystemsizehaving
176lipids.Dotted-dashedlinesrepresenttheresultswiththeoriginalsys-
P temsizehaving44lipids.Distributionswiththequartz/waterinterfacewere
ε =1+ z . (15)
ε E shiftedwithrespecttothemembranecentertocoincidewiththefreesol-
0 vatedbilayer.(a)DistributionsofmassdensitycontributedbyDMPClipids
Theestimatedvalueofε was103(±17),incloseagreement andwatermolecules.(b)(cid:4)cosθdipole,z(cid:5)valuesacrosstheinterface.(c)φ(z)
values.(cid:7)(d)E(z)va(cid:8)lues.(e)Orderparameterofacylchainsn-1calculatedby
withthevalueof94obtainedfromaseriesofsimulationsof
TIP3P bulk water with different electric fields.36 This con- SCD= 3cos22θ−1 ,whereθ istheanglebetweentheC–Hvectorandthez-
axis.CirclesandsquaresrepresentresultswithalargersystemsizeforDMPC
firms that Eq. (11) produces an electric field consistent with
bilayersinsolutionandwiththequartz/waterinterface,respectively.Pluses
theobservedorientationalpolarizationofwaterinsimulations representtheresultswiththeoriginalsystemsize.
performedwiththeEW3Dmethod.
C. Comparisonofbilayersinsolutionandatthe lipidacylchainorderparameter,37 respectively,inbulksolu-
quartz/waterinterface
tion, in the presence of the quartz/water interface, and with
We compared structural and electrostatic properties of a smaller system size. Even though small differences ex-
the DMPC bilayer solvated in bulk water and near the isted, the properties across the bilayers were, overall, very
quartz/water interface. Simulations with a larger system size similar. These results demonstrated that a quartz-supported
having 176 lipids were used for a more realistic compari- lipid bilayer deposited on top of a thin water film does not
son.TheresultsfromtheEW3DCsimulationwiththeorigi- introduce any gross differences to a completely solvated bi-
nal system size shown in Figs. 5 and 6 were also compared. layer. Thus, the simulations indicate that the deposited bi-
Figure 7 shows the results of this comparison. Distributions layercouldbeanequivalentmodeltostudythepropertiesof
for the solid-supported bilayer were shifted with respect to alipidbilayerinsolution.However,thepropertiesoflipidbi-
the average membrane center to facilitate comparisons with layers could be affected if a thinner water layer between the
thebulksimulation.Figure7(a)showsthatthecontributions solidsupportandthelipidbilayerismodeledassuggestedby
to the mass density distribution from DMPC were similar experiments9 and simulations.11 Most importantly, with the
in bulk water and near the quartz/water interface. Mass den- rigorous treatment of electrostatic boundary conditions de-
sitydistributionsofwaterwerealsosimilareventhoughthey veloped in this study for the quartz-supported lipid bilayer,
fluctuated around the bulk density near the quartz surface. avarietyoflipid/cholesterol/peptide/proteininteractionswith
Figures7(b)–7(e)showthecomparisonsoftheprofilesofpo- membranescannowbeaccuratelysimulatedwhilepreserving
larization of water, electrostatic potential, electric field, and theinherentanisotropyofthesystem.
055109-8 I.YehandA.Wallqvist J.Chem.Phys.134,055109(2011)
D. Quartz/waterinterfaceinbulksolution condition in Fig. 8(b) was more pronounced compared
with when embedded in vacuum. In addition, the water
The planar vacuum boundary condition or the EW3DC
density distribution near the hydrophobic side (z < −7 Å)
method is commonly applied to simulations of interfacial
inFig.8(a)wasalsoaffectedbythechoiceoftheelectrostatic
systems embedded in vacuum such as the quartz/water
boundary condition due to the high degree of ordering of
interfaceorthesolid-supportedlipidbilayerdescribedabove.
water with the EW3D method. The distribution of water
However,the P/ε terminEq.(11)isapplicabletoanysim-
z 0 density near the hydrophobic side was also slightly different
ulations performed under conducting boundary conditions,
compared with that obtained from a previous simulation26
regardless whether the simulated system is in bulk solution
whereelectrostaticinteractionsweretruncated.Asexpected,
orembeddedinvacuum.Therefore,weinvestigatedwhether
theunusualorderingofwaterwascompletelyremovedwhen
artificialPBC-inducedorientationalorderingwaspresentfor
we used the planar vacuum boundary condition with the
a fully solvated quartz interface using conducting boundary
EW3DC method. Distributions of the electrostatic potential
conditions and how it can be removed by implementing the
andelectricfieldshowninFigs.8(c)and8(d),calculatedwith
EW3DC planar vacuum boundary condition. Figures 8(a)
formulascorrespondingtotheemployedboundarycondition,
and 8(b), show a comparison of distributions of mass den-
were consistent with the (cid:4)cos θ (cid:5) values shown in
sity and orientational polarization of water, respectively, dipole,z
Fig. 8(b). A similar relationship between (cid:4)cos θ (cid:5) and
across the quartz/water interface in bulk solution obtained dipole,z
the electric field was observed in the previous simulation of
from simulations with the EW3D and EW3DC methods.
TIP3P bulk water.36 Our results demonstrate that the planar
Density distributions of water near the hydrophilic side (z
vacuum boundary condition can be applied to more general
> 7 Å) of the quartz crystal for both methods were similar
interfacialsystemsinbulksolutiontoremovethenonphysical
to those shown in Fig. 3(a) for the quartz/water interface
orderingobservedundertheconductingboundarycondition.
embedded in vacuum. However, the unusual ordering of
waterwiththeEW3Dmethodusingtheconductingboundary
IV. SUMMARYANDCONCLUSIONS
The general system consisting of a solvated bilayer ad-
(a)
3m) 1.5 sorbedonasolidsurfaceandexposedtoanair/vacuuminter-
c Water Water faceoccursinmanyexperimentalsettingsandpresentssome
sity (g/ 10..05 Quartz uhnowiqutheecsotmrupcututartailoannadlcehleacllteronsgteasti.cInprtohpisersttiuedsyo,fwpeoleaxriazmediniend-
n
e
D 0.0 terfacialsystemsareaffectedbythechoiceoftheelectrostatic
-30 -20 -10 0 10 20 30
boundaryconditionandtheimportanceofanalyzingelectro-
z (Å)
static properties with the correct formulas corresponding to
(b)
1.0 the deployed boundary conditions. We showed that imple-
>
θZdipole, 00..05 mboeunntdinagrythceonstdaintidoanrdfoErwinatledrfsaucmiaml asytisotnemwsitwhitthheaconnedtupcotlianrg-
os -0.5 ization introduced an abnormal orientational polarization of
c
<
-1.0 water. Our analyses suggest that the polarized component in
-30 -20 -10 0 10 20 30
z (Å) the simulated system (e.g., a polarized quartz crystal or an
(c) asymmetric bilayer) affects the properties of other parts of
15.0
thesystembyacompensationmechanismwhenweapplythe
V 10.0 0.0 conducting boundary condition. Implementing the Ewald
z(), φ 5.0 --00..21 EEWW33DDC summationtechniquewiththeplanarvacuumboundarycon-
-30 -20 -10 dition and evaluating electrostatic properties with formulas
0.0
-30 -20 -10 0 10 20 30 corresponding to the electrostatic boundary condition does
z (Å) notintroduceanysuchcompensatoryeffect,andtheinfluence
(d) of the polarized component is more realistically accounted
0.04 for.Wealsoconfirmedthatthemodeledquartz-supportedhy-
Å 0.02
V/ drated lipid bilayer reproduced the essential features of the
Ez(), -00..0020 lipidbilayerinbulkwater.
-0.04
-30 -20 -10 0 10 20 30
z (Å) ACKNOWLEDGMENTS
We thank Dr. Ramanathan Nagarajan and Dr. Hyung-
FIG.8. Structuralandelectrostaticpropertiesofthequartz/waterinterface
in bulk solution without vacuum interfaces. Results with the EW3D and JuneWoofortheirhelpfulsuggestionsanddiscussions.This
EW3DCmethodsareshownbythedashedandsolidlines,respectively.Dis- project was funded in part by a competitive In-house Labo-
tributions were calculated with respect to the center of the quartz crystal. ratory Independent Research (ILIR) grant by the U.S. Army
(a)Distributionsofmassdensitycontributedbywatermolecules.(b)(cid:4)cos
Assistant Secretary of the Army for Acquisition, Logistics,
θdipole,z(cid:5)valuesacrosstheinterface.(c)φ(z)values.(d)E(z)values.Theinset
inCshowsamagnifiedviewofφ(z)profilesinthebulklikeregionofthe and Technology (ASAALT). Funding support for this work
waterphase. also came from the Department of Defense (DoD) High
055109-9 Electrostaticsinsolid-supportedbilayer J.Chem.Phys.134,055109(2011)
Performance Computing (HPC) Modernization Program Of- 14E.R.Smith,Proc.R.Soc.London,Ser.A375,475(1981).
fice,undertheHPCSoftwareApplicationsInstituteinitiative, 15H.D.Herce,A.E.Garcia,andT.Darden,J.Chem.Phys.126(12),124106
(2007); B. P. van Eijck and J. Kroon, J. Phys. Chem. B 101(6), 1096
the U.S. Army Medical Research and Materiel Command.
(1997).
ComputationaltimewasprovidedbytheU.S.ArmyResearch 16I.-C.YehandM.L.Berkowitz,J.Chem.Phys.111(7),3155(1999).
Laboratory and Navy DoD Supercomputing Resource Cen- 17D. E. Parry, Surf. Sci. 54, 195 (1976); D. M. Heyes, M. Barber, and
ters. The opinions or assertions contained herein are the pri- J.H.R.Clarke,J.Chem.Soc.,FaradayTrans.273,1485(1977);E.Spohr,
J.Chem.Phys.107,6342(1997).
vateviewsoftheauthorsandarenottobeconstruedasofficial
18J.JanecˇekandR.R.Netz,Langmuir23(16),8417(2007).
orasreflectingtheviewsoftheU.S.ArmyoroftheU.S.De- 19V.Buch,A.Milet, R.Vacha, P.Jungwirth, andJ. P.Devlin,Proc. Natl.
partmentofDefense.Thispaperhasbeenapprovedforpublic Acad.Sci.U.S.A.104(18),7342(2007);D.BostickandM.L.Berkowitz,
release. Biophys.J.85(1),97(2003).
20J.N.Sachs,P.S.Crozier,andT.B.Woolf,J.Chem.Phys.121(22),10847
(2004).
21I.-C.YehandG.Hummer,Proc.Natl.Acad.Sci.U.S.A.101(33),12177
1B.HonigandA.Nicholls,Science268(5214),1144(1995);N.Sinhaand (2004).
S.J.Smith-Gill,Curr.ProteinPept.Sci.3(6),601(2002);P.Koehl,Curr. 22K.A.Marx,Biomacromolecules4(5),1099(2003).
Opin.Struct.Biol.16(2),142(2006);K.D.Collins,G.W.Neilson,andJ. 23L. Kale, R. Skeel, M. Bhandarkar, R. Brunner, A. Gursoy, N. Krawetz,
E.Enderby,Biophys.Chem.128(2–3),95(2007). J.Phillips,A.Shinozaki,K.Varadarajan,andK.Schulten,J.Comput.Phys.
2M.P.AllenandD.J.Tildesley,ComputerSimulationofLiquids(Oxford 151(1),283(1999).
UniversityPress,NewYork,1987). 24A.D.MacKerell,D.Bashford,M.Bellott,R.L.Dunbrack,J.D.Evanseck,
3P.Ewald,Ann.Phys.64,253(1921). M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy,
4L. Greengard and V. Rokhlin, J. Comput. Phys. 73(2), 325 (1987); J. L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo,
Lekner,PhysicaA176,485(1991);X.WuandB.R.Brooks,J.Chem. D.T.Nguyen,B.Prodhom,W.E.Reiher,B.Roux,M.Schlenkrich,J.C.
Phys.122(4),044107(2005). Smith,R.Stote,J.Straub,M.Watanabe,J.Wiorkiewicz-Kuczera,D.Yin,
5S.W.deLeeuw,J.W.Perram,andE.R.Smith,Proc.R.Soc.London,Ser. andM.Karplus,J.Phys.Chem.B102(18),3586(1998).
A373,27(1980). 25W.L.Jorgensen,J.Chandrasekhar,J.D.Madura,R.W.Impey,andM.L.
6S.W.deLeeuw,J.W.Perram,andE.R.Smith,Proc.R.Soc.London,Ser. Klein,J.Chem.Phys.79,926(1983).
A388,177(1983). 26P.E.M.Lopes,V.Murashov,M.Tazi,E.Demchuk,andA.D.MacKerell,
7T.A.Darden,D.M.York,andL.G.Pedersen,J.Chem.Phys. 98,10089 J.Phys.Chem.B110(6),2782(2006).
(1993);U.Essmann,L.Perera,M.L.Berkowitz,T.Darden,H.Lee,andL. 27J.P.Ryckaert,G.Ciccotti,andH.J.C.Berendsen,J.Comput.Phys.23,
G.Pedersen,J.Chem.Phys.103,8577(1995). 327(1977).
8M.Karttunen,J.Rottler,I.Vattulainen,C.Sagui,andE.F.Scott,Current 28I.C.YehandA.Wallqvist,J.Phys.Chem.B113(36),12382(2009).
TopicsinMembranes(Academic,NewYork,2008),Vol.60,p.49. 29S. Jo, T. Kim, V. G. Iyer, and W. Im,J. Comput. Chem. 29(11), 1859
9C.E.Miller,J.Majewski,T.Gog,andT.L.Kuhl,Phys.Rev.Lett.94(23), (2008).
238104(2005);E.B.Watkins,C.E.Miller,D.J.Mulder,T.L.Kuhl,and 30H. I. Petrache, S. W. Dodd, and M. F. Brown, Biophys. J. 79(6), 3172
J. Majewski, Phys. Rev. Lett. 102(23), 238101 (2009); T. H. Anderson, (2000).
Y.J.Min,K.L.Weirich,H.B.Zeng,D.Fygenson,andJ.N.Israelachvili, 31B.R.Brooks,R.E.Bruccoleri,B.D.Olafson,D.J.States,S.Swaminathan,
Langmuir25(12),6997(2009). and M. Karplus, J. Comput. Chem. 4, 187 (1983); B. R. Brooks, C. L.
10D.R.Heine,A.R.Rammohan,andJ.Balakrishnan,Mol.Simul.33(4– Brooks,A.D.Mackerell,L.Nilsson,R.J.Petrella,B.Roux,Y.Won,G.
5), 391(2007); S. V. Bennun,A. N. Dickey, C. Y. Xing, and R. Faller, Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R.
FluidPhaseEquilib.261(1–2),18(2007);C.Y.XingandR.Faller,J.Phys. Dinner,M.Feig,S.Fischer,J.Gao,M.Hodoscek,W.Im,K.Kuczera,T.
Chem.B112(23),7086(2008);C.Y.Xing,O.H.S.Ollila,I.Vattulainen, Lazaridis,J.Ma,V.Ovchinnikov,E.Paci,R.W.Pastor,C.B.Post,J.Z.Pu,
andR.Faller,SoftMatter5(17),3258(2009);C.Y.XingandR.Faller, M.Schaefer,B.Tidor,R.M.Venable,H.L.Woodcock,X.Wu,W.Yang,
ibid.5(22),4526(2009);C.Y.XingandR.Faller,J.Chem.Phys.131(17), D.M.York,andM.Karplus,ibid.30(10),1545(2009).
175104(2009);M.I.Hoopes,M.Deserno,M.L.Longo,andR.Faller,J. 32A. A. Gurtovenko and I. Vattulainen, J. Chem. Phys. 130(21), 215107
Chem.Phys.129(17),175102(2008). (2009).
11M. Roark and S. E. Feller, Langmuir 24(21), 12469 (2008); S. V. 33M. C. F. Wander and A. E. Clark, J. Phys. Chem. C 112(50), 19986
Bennun, M. I. Hoopes, C. Y. Xing, and R. Faller, Chem. Phys. Lipids (2008).
159(2), 59 (2009); M. I. Hoopes, C. Xing, and R. Faller, Biomembrane 34S.H.LeeandP.J.Rossky,J.Chem.Phys.100(4),3334(1994).
Frontiers,editedbyR.Faller,M.L.Longo,S.H.Risbud,andT.Jue(Hu- 35R. S. Neves, A. J. Motheo, R. P. S. Fartaria, and F. Fernandes, J. Elec-
mana,NewYork,2009),p.101. troanal.Chem.612(2),179(2008);F.Goujon,C.Bonal,B.Limoges,and
12S.E.Feller,R.W.Pastor,A.Rojnuckarin,S.Bogusz,andB.R.Brooks,J. P.Malfreyt,Langmuir25(16),9164(2009).
Phys.Chem.100,17011(1996);I.-C.YehandM.L.Berkowitz,J.Chem. 36I.-C.YehandG.Hummer,Biophys.J.86(2),681(2004).
Phys.110(16),7935(1999). 37L.S.Vermeer,B.L.deGroot,V.Reat,A.Milon,andJ.Czaplicki,Eur.
13J.C.ShelleyandG.N.Patey,Mol.Phys.88,385(1996). Biophys.J.36(8),919(2007).
