Surface Tension Parameterization in Molecular Dynamics Simulations of a Phospholipid-bilayer Membrane: Calibration and Effects

Abstract

Molecular dynamics simulations of interfacial systems such as phospholipid bilayers require experimental estimates of either the surface area per lipid or the surface tension (y). To investigate the mode of selection of yvalues, we have carried out the molecular dynamics (MD) simulations described here for a neat DMPC bilayer comprising 90 lipids using NP N γT and NP N AT ensembles. Influence of y values coupled with the effect of force truncation for both NP N γT and NP N AT ensembles was analyzed in multiple simulations. The membrane structural properties were studied comparatively using different y values in NP N γT ensembles with a spherical cutoff of 18 A and particle mesh Ewald method (PME) that takes long-range electrostatic interactions into account. In parallel, we tested the sensitivity of calculated y value in simulations using NP N AT ensembles with four different spherical cutoffs and PME. Our studies suggest that when PME is used to calculate the nonbonded forces, the value of surface tension for the system of this size should be in the range of 25 to 30 dyn/cm (per interface). On the other hand, for simulations with a cutoff of 18 A, shrinkage of box size of about 7-10% is observed for y values ranging from 0 to 85 dyn/cm. The surface tension value has to be raised to 625 dyn/cm to keep the surface area/lipid close to the experimental estimate for neat DMPC bilayers. Molecular dynamics simulations were also carried out with the opioid peptide dynorphin A(1-13) embedded within the DMPC bilayers. Several y values were tested with two electrostatic schemes. Our results demonstrated that for y values 0, 45, and 65 dyn/cm, due to the shrinkage of box size, the α-helical structure in the N-terminal region of the peptide is distorted with an 18 A cutoff. When the y value is increased to 625 dyn/cm with an 18 A cutoff or when PME is used with y = 25 dyn/cm, the peptide's α-helical structure is maintained and its orientation is similar to that obtained from earlier simulations 3 5 using NVE ensembles. The aromatic side chains in the peptide (Tyr-1 and Phe-4) interacted similarly in all the simulations except when γ=0 is used with a cutoff of 18 A. All constant surface tension simulations were compared with constant area ensembles in which the surface area/lipid remains fixed. Our results suggest that neglect of long-ranged electrostatic interactions and/or incorrect surface tension value results in serious artifacts in membrane simulations and suggest that caution is required when constant surface tension simulations are carried out on neat bilayers. In simulations of membranes with an embedded solute, the choice of an appropriate y value is more complex, as it will affect both the structure and the interactions of the solute within the bilayer.