Methodological problems in pressure profile calculations for lipid bilayers

Abstract

From molecular dynamics simulations of a dipalmitoyl-phosphatidyl-choline (DPPC) lipid bilayer in the liquid crystalline phase, pressure profiles through the bilayer are calculated by different methods. These profiles allow us to address two central and unresolved problems in pressure profile calculations: The first problem is that the pressure profile is not uniquely defined since the expression for the local pressure involves an arbitrary choice of an integration contour. We have investigated two different choices leading to the Irving-Kirkwood (IK) and Harasima (H) expressions for the local pressure tensor. For these choices we find that the pressure profile is almost independent of the contour used, which indicates that the local pressure is well defined for a DPPC bilayer in the liquid crystalline phase. This may not be the case for other systems and we therefore suggest that both the IK and H profiles are calculated in order to test the uniqueness of the profile. The second problem is how to include electrostatic interactions in pressure profile calculations when the simulations are conducted without truncating the electrostatic potential, i.e., using the Ewald summation technique. Based on the H expression for the local pressure, we present a method for calculating the contribution to the lateral components of the local pressure tensor from electrostatic interactions evaluated by the Ewald summation technique. Pressure profiles calculated with an electrostatic potential truncation (cutoff) from simulations conducted with Ewald summation are shown to depend on the cutoff in a subtle manner which is attributed to the existence of long-ranged charge ordering in the system. However, the pressure profiles calculated with relatively long cutoffs are qualitatively similar to the Ewald profile for the DPPC bilayer studied here.

Pressure Profile Calculation with Mesh Ewald Methods

“…Note that, because the contribution in reciprocal space is not pairwise additive, it is not possible to use an Irving-Kirkwood formulation either for the Ewald sum, or for its mesh-based approximations. 11 We now proceed to extending the reciprocal space contribution to mesh-based algorithms, by noting that Eq. (9) is, for µ = ν, the discrete inverse Fourier transform of the complex function S(m)f (m 2 )g µν (m), evaluated at the particles' positions, and scaled by the factor q i .…”

“…20 ). Sonne and collaborators 11 showed that the pressure profile calculated using the Irving-Kirkwood approach and a cutoff radius of 2 nm is qualitatively similar to the one calculated using the Harasima path with long-range electrostatic contributions, although with an indetermination of roughly 100 to 200 bar, which makes it only a qualitative test.…”

“…The Ewald summation 14 can, in principle, be used in such systems; the contribution of its correction term to the pressure profile was derived by Sonne and coworkers for the Harasima path. 11 They also showed that potentials that are not strictly pairwise additive, such as the correction term of the Ewald summation, cannot be used in combination with the Irving-Kirkwood profile. However, because of its poor scaling properties, the use of the Ewald summation is prohibitive for systems with a large number of charges.…”

“…In other words, a certain amount of the total lateral pressure is associated with each particle, making the calculation of the lateral pressure profile straightforward and computationally efficient. 11,12 Another serious technical problem in calculating pressure profiles in charged systems is how to take into ac- One molecule identified to be in the vapour phase by a clustersearch algorithm can be also seen above the first layer of the surface count the long range correction term of the electrostatic interaction. It is now well-established that neglecting this correction to the energy and forces can lead to important systematic errors.…”

A direct two-dimensional pressure formulation in molecular dynamics

“…Many simulated systems have inhomogeneity only in two dimensions (2D) such as defect nucleation in bulk and two-dimensional crystals, bio molecular assemblies such as lipid bilayers and membrane proteins, as well as thin film evaporation and heat transfer, and thus only require 2D pressure distribution. However, current literature on local pressure estimation is based on 3D [3] or 1D [4][5][6] pressure estimation. The 2D pressure distribution is obtained by averaging over the 3D pressure data, and is extremely computationally expensive [7] as it involves a 3D convolution.…”