All simulations were performed
in the isothermal–isobaric ensemble, NPT,
at a pressure of 1 atm. The pressure was held constant by using the
Parrinello–Rahman barostat77 with
a coupling constant of 10.0 ps with an isothermal compressibility
of 4.5 × 10–5 bar–1. For
the bulk liquids an isotropic pressure coupling was used and for the
bilayer simulations a semi-isotropic pressure coupling scheme was
used. The temperature was kept constant by the Nosé–Hoover
thermostat78 ,79 (link) with a coupling constant of 0.5 ps. The
lipid bilayer and water were coupled separately to the thermostat.
Long-range electrostatic interactions were treated by a particle-mesh
Ewald scheme80 ,81 with a real-space cutoff at 1.4
nm with a Fourier spacing of 0.10 nm and a fourth-order interpolation
to the Ewald mesh. Single-atom charge groups were used. van der Waals
interactions were truncated at 1.5 nm and treated with a switch function
from 1.4 nm. Long-range corrections for the potential and pressure
were added.51 The inclusion of long-range
corrections should eliminate the LJ cutoff dependency in the simulations.
Due to the fact that lipid bilayers are inhomogeneous systems the
method introduced by Lagüe et al.82 to add long-range corrections could be applied instead. Periodic
boundary conditions were imposed in every dimension. A time step of
2 fs was used with a Leap-Frog integrator. The LINCS algorithm83 was used to freeze all covalent bonds in the
lipid, and the analytical SETTLE84 method
was used to hold the bonds and angle in water constant. The TIP3P
water model85 was the water model of choice.
The choice of water model can be explained by the fact that TIP3P
is the default water model in major FFs such as AMBER and CHARMM and
since one of the aims of the work presented here was to create a lipid
FF compatible with AMBER this was a natural choice. Further, earlier
work of Högberg et al.31 (link) has shown
that there is flexibility in the choice of water model for AA simulations
of lipid bilayers. Atomic coordinates were saved every 1 ps and the
neighbor list was updated every 10th step.
Bulk liquids were
simulated with a simulation box consisting of 128 molecules for the
larger alkanes and 256 for the smaller alkanes (hexane and heptane)
at a temperature of 298.15 K. The lipid bilayer systems were prepared
using the CHARMM-GUI86 (link),87 (link) with 128 lipids in total, 64
in each leaflet. In order to achieve proper hydration, 30 TIP3P water
molecules were added per lipid. Three different lipid types were simulated,
DLPC (12:0/12:0), DMPC (14:0/14:0), and DPPC (16:0/16:0). These system
were investigated under a range of temperatures; see Table 1 for an overview of all simulations performed. All
lipid bilayer systems were equilibrated for 40 ns before production
runs were initiated which lasted for 300–500 ns. All MD simulations
were performed with the Gromacs88 software
package (versions 4.5.3 and 4.5.4). All analysis were made with the
analysis tools that come with the MDynaMix software package.89 System snapshots were rendered and analyzed
with VMD.90 Neutron scattering form factors
were computed with the SIMtoEXP software.91 (link)The calculations of free energies of solvation in
water and cyclohexane
were performed by using thermodynamic integration over 35 λ
values in the range between 0 and 1. A soft core potential (SCP) was
used to avoid singularities when the solute is almost decoupled from
the solvent. The α-parameters used for the SCP and the simulation
workflow were set following the methodology described by Sapay and
Tieleman.92 (link) The amino acid analogues were
solvated with 512 and 1536 molecules of cyclohexane and water, respectively.