Molecular Dynamics Simulations of Phospholipid Bilayers Under Deformation: A Comparison Between GROMACS & LAMMPS

INTRODUCTION

While macroscale anatomic lesion or organ dysfunction due to injury can often be identified in experimental or clinical environments, nanoscale cellular changes are more difficult to detect and can require in silico approaches. Among different in silico methods, multiple molecular dynamics (MD) studies have been conducted. By computationally modeling molecules and their dynamic evolution in time, MD elucidates cell membrane injury biomechanics and damage mechanisms under mechanical deformation, particularly with regards to stress-strain, mechanoporation and rupture response of the phospholipid bilayers - the fundamental biological structure of cell membranes. Several MD simulation software packages have been developed in both engineering and biomedical fields. However, depending on the MD simulator chosen, the examination of nanoscale deformation mechanisms of cellular structures could render drastically different results. Also, comparisons between these different MD simulators is typically an intricate task, requiring that all configurations be converted appropriately with a variety of reasonable parameter choices. Accordingly, the current study aims to perform and compare the results between two common MD programs (GROMACS and LAMMPS) for simulating a representative cellular membrane structure under deformation.

METHODS

A bilayer membrane, which consisted of 72 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC) phospholipids with a thick layer of approximately 9070 TIP3P water molecules, was equilibrated for ten nanoseconds using the CHARMM36 all-atom lipid force field. The equilibrated structure was then deformed under different strain states (equibiaxial, 2:1 non-equibiaxial, 4:1 non-equibiaxial, strip biaxial and uniaxial tension). Constant deformation velocities were applied in x and y dimensions, while the z dimension was allowed to adjust freely under a pressure of 1 atm and a temperature of 310 K. The results for the von Mises stress-strain, pore nucleation and growth, and damage behavior are compared between the respective GROMACS and LAMMPS simulations.

RESULTS AND CONCLUSIONS

In general, GROMACS and LAMMPS produced similar deformation behaviors, including for pore formation, damage evolution, and strain state detrimentality order. The general detrimental order was found as equibiaxial, 2:1 non-equibiaxial, 4:1 non-equibiaxial, strip biaxial, and lastly, uniaxial, of which the equibiaxial and 2:1 non-equibiaxial strain states caused the most damage and the uniaxial strain state led to the least damage. However, GROMACS produced lower stress values and higher damage values than LAMMPS. Multiple different setting options between GROMACS and LAMMPS, including barostat/thermostat algorithm variations, differing parameters for the water models, numerical precision and options for the implementation of long-range interactions, have been considered as potential factors leading to the observed differences. Overall, this study will aid in the cross-check of parameter settings and simulation results in future MD research, particularly on the mechanical deformation of phospholipid bilayers and other mechano-physiological damage mechanisms of cellular systems. Thereby, in future efforts, GROMACS and LAMMPS, as well as other MD software packages, could be exploited synchronously with better comparability and reproducibility.