Regular paper

Condensed phase properties of n-pentadecane emerging from application of biomolecular force fields*

Maciej Bratek, Anna Wójcik-Augustyn, Adrian Kania, Jan Majta and Krzysztof Murzyn✉

Department of Computational Biophysics and Bioinformatics; Faculty of Biochemistry, Biophysics and Biotechnology; Jagiellonian University; Kraków, Poland

For over 20 years, the OPLS-All Atom (OPLS-AA) force field has been efficiently used in molecular modelling studies of proteins, carbohydrates and nucleic acids. OPLS-AA is successfully applied in computer modelling of many organic compounds, including decane and shorter alkanes, but it fails when employed for longer linear alkanes, whose chemical structure corresponds to hydrocarbon tails in phospholipids constituting cellular membranes. There have been several attempts to address this problem. In this work, we compare the ability to reproduce various condensed phase properties by six distinct sets of force field parameters which can be assigned to phospholipid hydrocarbon chains. In this comparison, we include three alternative sets of the OPLS-AA force field, as well as the commonly used CHARMM C36, Slipids, and Berger lipids’ parameters.

Key words: force field, n-pentadecane, condensed phase properties, lipid bilayers, phospholipids, molecular dynamics simulation

Received: 24 April, 2020; revised: 17 June, 2020; accepted: 25 June, 2020; available on-line: 17 September, 2020

✉e-mail: krzysztof.murzyn@uj.edu.pl

*Presented at the XLVII Winter School of the Faculty of Biochemistry, Biophysics and Biotechnology of the Jagiellonian University “Molecules, Pathways, and Games”, February 8–12, 2020, Zakopane, Poland.

Acknowledgements of Financial Support: This research project was supported by grants No 2011/01/B/NZ1/00081 (KM) and 2015/17/D/ST4/00555 (AWA) from the National Science Center (Poland) and in part by the PL-Grid Infrastructure. The authors have declared no conflicts of interest for this article.

Abbreviations: MD, molecular dynamics; FF, force field; PC, phosphatidylcholine; OPLS-AA, Optimized Potentials for Liquid Simulations All-Atom; Slipids, Stokholm lipids’ force field; BergerLips, Berger lipids’ force field; PE, phosphatidylethanolamine; SM, sphingomyelin; PS, phosphatidylserine; PG, phosphatidylglycerol; DPPC, dipalmitoylphosphatidylcholine; POPE, 1-palmitoyl-2-oleoylphosphatidylethanolamine; POPC, 1-palmitoyl-2-oleoylphosphatidylcholine; DPPE, dipalmitoylphosphatidylethanolamine

INTRODUCTION

Long aliphatic chains are important parts of phospholipids, glycolipids and sphingolipids, which constitute biological membranes. Correct reproduction of physicochemical properties of the hydrophobic core of a lipid membrane is critical for obtaining a valid computer model of the whole membrane, as evidenced by several works (Piggot et al., 2012; Siu et al., 2012; Kahn & Bruice, 2002). Over recent years, several biomolecular force fields (FF) have been developed and applied in molecular modelling of membrane lipids with varying success. One of the first commonly used phospholipid specific united-atom parameter sets was BergerLips (Berger et al., 1997). This set employed the bonded FF parameters from GROMOS87 (van Gunsteren & Berendsen, 1987), a set of newly fitted dihedral parameters for lipid aliphatic chains, and an eclectic choice of nonbonded parameters following combination rules from the OPLS united atom FF (Jorgensen et al., 1984). More recent approaches extend all-atom FF and include the CHARMM C36 (Klauda et al., 2010), Stockholm lipids (Slipids) (Jämbeck & Lyubartsev, 2012; Jämbeck & Lyubartsev, 2013), and OPLS-AA (Jorgensen et al., 1996; Murzyn et al., 2013).

The CHARMM C36 parameter set is the answer to problems that afflicted the previous version (C27) of the CHARMM FF, i.e. systematic overestimation of phospholipid bilayer melting points and underestimation of a surface area per lipid, and thus errors in the predicted deuterium order parameters (Klauda et al., 2010). Refitting of selected torsional and nonbonded parameters rendered CHARMM FF usable with typical phospholipids. One of the goals of the CHARMM FF development is to retain compatibility between the CHARMM and NAMD software packages (Lyubartsev et al., 2015). However, implementation of the CHARMM C36 parameters in the NAMD package gave unsatisfactory results (Klauda et al., 2010). Moreover, MD simulations with CHARMM C27 and C36 assigned parameters performed for lipid bilayers built of phosphatidylcholines with polyunsaturated fatty acid chains led to poor reproduction of various structural properties of lipid membranes when carried out in the isothermal-isobaric ensemble, such as the deuterium order parameters, X-ray form factors and surface area per lipid (Klauda et al., 2010). This study indicated that refinement of parameters for hydrocarbon chains is still an ongoing task (Klauda et al., 2010; Lyubartsev et al., 2015). Indeed, further improvements of CHARMM FF were made in the C36p set (Klauda et al., 2012; Lyubartsev et al., 2015).

Slipids parameters are a modern addition to the lipid force field family. Initially, they were developed for saturated phospholipids (Jämbeck & Lyubartsev, 2012) and then extended to include unsaturated phosphatidylcholine (PC) and phosphatidylethanolamine (PE) (Jämbeck & Lyubartsev, 2012), as well as sphingomyelin (SM), phosphatidylglycerol (PG), phosphatidylserine (PS), and cholesterol (Jämbeck & Lyubartsev, 2013; Ermilova & Lyubartsev, 2016). Inclusion of extensively validated parameters for a broad range of lipid types was a big advantage of Slipids, as many of the mainstream FF had mainly targeted PC and PE. It should be also noted that despite being focused solely on lipids, the Slipids parameters can be used in simulations of multicomponent biomolecular systems (Jämbeck & Lyubartsev, 2012), with nonlipid molecules being assigned the very popular and reliable AMBER all-atom FF parameters (Maier et al., 2015).

OPLS-AA is a generic FF, originally developed for small organic liquids (Jorgensen et al., 1996), and extended to nucleic acids (Pranata et al., 1991), carbohydrates (Damm et al., 1997; Kony et al., 2002), small heterocyclic compounds (McDonald & Jorgensen, 1998), amines (Rizzo & Jorgensen, 1999), and proteins (Kaminski et al., 2001). Although assignment of the generic OPLS-AA parameters to membrane phospholipids is possible, it results in a severe overestimation of the membrane phase transition temperatures and unsatisfactory reproduction of many other physicochemical properties of lipid bilayers (Siu et al., 2012; Kahn & Bruice, 2002). The OPLS-AA parameters have been occasionally refined. Studies on glycolipids, methyl acetate, dimethyl phosphate and n-pentadecane enabled improvements in condensed phase properties, such as densities, the enthalpy of vaporization and the Gibbs free energy of hydration (Murzyn et al., 2013). The refined OPLS-AA parameters were successfully used in MD simulations performed for a model of E. coli lipid A (Murzyn & Pasenkiewicz-Gierula, 2015) and dipalmitoylphosphatidylcholine (DPPC) (Murzyn et al., unpublished results). Another set of OPLS-AA parameters, which has been used in several MD studies on lipid bilayers (Pasenkiewicz-Gierula et al., 2012, Plesnar et al., 2012, Plesnar et al. 2013), was obtained by alteration of partial charges on the C and H atoms in order to reproduce the C-H bond dipole moment reported by Vereshchagin and Vul’fson (Vereshchagin & Vul’fson, 1967).

Throughout the history of development of FF parameters for lipids, n-alkanes have been traditionally used as model molecules, representing the aliphatic chains of membrane lipids (Murzyn et al., 2013). For instance, n-pentadecane (C15H32) may be used as a model of the palmitoyl chain, present in some of the most popular bio-membrane components: 1-palmitoyl-2-oleoylphosphatidylethanolamine (POPE), 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) and DPPC. It is generally acknowledged that good reproduction of bulk alkane properties is indicative of correct parametrization of the hydrophobic region of membrane lipids (Ryckaert & Bellemans 1975; Ryckaert & Bellemans 1978; Tuteja et al., 2007; Siu et al., 2012). Fortunately, a decent amount of experimental physical chemistry data is available for liquid hydrocarbons (Ben-Naim and Marcus, 1984; Haynes et al., 2011). As evidenced by the history of development of the OPLS-AA, the length of alkanes taken for validation is significant. The OPLS-AA parameters were originally validated only with hexane and shorter alkanes, resulting in poor reproduction of the long alkanes’ melting points, and consequently, the transition temperature of phospholipid bilayers.

Apart from the most obvious applications of biomolecular FF for membrane lipids, i.e. simulations of membranes at physiological temperatures, many interesting studies have taken place. Leekumjorn and Sum (Leekumjorn & Sum, 2007) have used MD simulations in combination with BergerLips FF to study the gel (L’β) to liquid crystal (Lα) phase transitions of DPPC and dipalmitoylphosphatidylethanolamine (DPPE) membranes, and characterized key bilayer properties in this context. Stepniewski and co-workers (Stepniewski et al., 2010) have conducted OPLS-AA simulations of phosphatidylcholines in the gel (L’β) and liquid crystal (Lα) phases and presented impact of the phase state on hydration and permeation by ions of the bilayer. This highlights the need for lipid-oriented force fields to be also applicable at non-physiological temperatures. Jämbeck and Lyubartsev (Jämbeck & Lyubartsev, 2012) have demonstrated the ability of Slipids to simulate either the gel or the liquid crystal phase of DPPC .

The aim of this study was to compare the quality of reproduction of several structural, thermodynamic and kinetic properties of n-pentadecane liquid phase, as modelled by popular biomolecular force fields employed in simulations of lipid bilayers. The molecule of n-pentadecane is of special interest in the area of lipid membrane simulations, being the equivalent of a typical saturated phospholipid chain. Validation of liquid phase alkane simulations primarily serves to detect and investigate any deficiencies in the quality of parametrization of long aliphatic chains in these force fields, whether employed in simulations of lipid membranes or any other systems.

This paper is organized as follows: in the Materials and Methods section we describe in detail the FFs investigated in our work, and then we outline crucial parameters used in the performed simulations. This is followed by a brief description of computational procedures involved in the comparison of FFs. In the Results and Discussion section, we analyze the outcome of simulations of n-pentadecane using chosen parametrizations to evaluate the quality of reproduction of critical thermodynamic, conformational and kinetic properties of the bulk phase for each of them. We try to put these results in a wider scope of lipid parameterization and draw comparisons between our work and other significant contributions to the field. In the Conclusion section, we summarize the main points established in this work.

MATERIALS AND METHODS

FF parameters. In our calculations of physicochemical properties of n-pentadecane we used three FF derived from OPLS-AA, labelled as follows: (1) OPLS-RAW: generic OPLS-AA FF parameters (Jorgensen et al., 1996), (2) OPLS-QQ: a set built on top of OPLS-RAW with –0.24e on carbon atoms of hydrocarbon chains (Pasenkiewicz-Gierula et al., 2012), and (3) OPLS-MP: an updated set of OPLS-AA FF parameters (Murzyn et al., 2013).

The general strategy for adapting parameters of n-pentadecane in three other FFs, which are included in our comparative study, relied on using topologies published for the palmitoyl chain of a phospholipid, e.g. POPE. In this way we conserve the link between simulation of n-pentadecane and validation of a FF dedicated for lipids. For Slipids, we downloaded both, the FF parameters and the recommended simulation protocol from the authors’ website. In the case of CHARMM C36, we downloaded the GROMACS adaptation of FF from the GROMACS website (section USER CONTRIBUTIONS). The BergerLips parameters were assigned to n-pentadecane according to Justin Lemkul’s tutorial on simulation of KALP­15 protein in a DPPC bilayer (Lemkul, 2012).

In calculations of Gibbs free energy of hydration of n-pentadecane, we used water models which are recommended for a given FF recognizing the fact that such models are integral parts of appropriate FFs (Pastor & MacKerell, 2011). Thus, we used the rigid TIP3P water model (Jorgensen, 1981) for all OPLS-AA FF sets and Slipids (Ermilova & Lyubartsev, 2016), the CHARMM-specific TIPS3P water model (Neria et al., 1996) for the CHARMM C36 set, as suggested by Piggot (Piggot et al., 2012), and the Simple Point Charge (SPC) water model (Berendsen et al., 1981) in case of the BergerLips set (van Buuren et al., 1993).

Simulation parameters. All the simulations were performed with the GROMACS package, version 4.5 (van Der Spoel et al., 2005; Hess et al., 2008). Simulations of the isothermal-isobaric (NPT) ensemble were 100 ns long, while the canonical ensemble (NVT) simulations were 300 ns long to ensure proper sampling needed for calculation of transport properties (Murzyn et al., 2013). The first 50 ns and 200 ns, respectively, were treated as an equilibration period and not considered during the trajectory analysis.

The temperature was maintained at either 277K or 298K with the Nose-Hoover algorithm (Nose, 1984; Hoover, 1985), with a time constant of 0.5 ps. The pressure was maintained at 1 bar using the Parrinello-Rahman algorithm (Parrinello & Rahman, 1981), with a time constant of 5 ps. In the condensed phase simulations, long-range electrostatics were treated with the Particle Mesh Ewald (PME) algorithm. Table 1 lists simulation parameters used for all considered FFs, as reported in the reference papers. In a recent work, Piggot and co-workers have highlighted several discrepancies within published protocols and parameters of many popular FFs for lipids, including CHARMM C36 and BergerLips (Piggot et al., 2012). In our work, we follow the Piggot’s recommendation for GROMACS MD simulations, as far as these two FFs are considered.

Condensed phase properties. The starting structure of 162 molecules of n-pentadecane for 277K and 298K condensed phase simulations was taken from an MD simulation at a temperature of 285K (Fig. 2), containing structural features reminiscent of both, the solid and liquid phases, to facilitate transition into either of them, depending on the MD simulation conditions and the choice of FF parameters.

The densities were calculated by averaging box volumes in the NPT MD simulations.

The enthalpy of vaporization (ΔHvap) was calculated according to the following equation:

where Egaspot and Eliqpot are potential energies in the gas and liquid state, respectively, R is the gas constant, and T is the absolute temperature. Additional gas phase simulations were performed with settings corresponding to those for condensed phase simulations; however, simple cutoff electrostatics were used instead of PME, appropriate for a single molecule in a vacuum. Because of the lack of periodic boundary condition option, no dispersion correction was used with BergerLips and Slipids, and the switch functions were replaced by unified cutoffs (1.4 nm for Berger, 1.2 nm for CHARMM C36). Block averages’ error estimates were calculated for both density and potential energies, using the g_analyze program from the GROMACS package (Hess, 2002).

The Gibbs free energies of hydration (ΔGhyd) were calculated with the Bennett Acceptance Ratio (BAR) method, implemented in GROMACS’s g_bar program (Bennett, 1976; Van Der Spoel et al., 2005; Hess et al., 2008). For each system, 40 lambda values, determining intermediate steps, were used, evenly spaced out by 0.025. The Coulomb and van der Waals parameters were perturbed simultaneously and a soft-core function was used for van der Waals interaction to prevent energy singularities at low distances, using parameters p = 1; σ = 0.3; α = 1. For each lambda value, a 1 ns simulation was performed, using the Stochastic Dynamics integrator. The remaining MD run parameters were identical to the regular condensed phase simulations described above. For our ΔGhyd calculations, a box of 1000 water molecules and one n-pentadecane molecule was used. Several additional simulations were performed to ensure that simultaneous nonbonded interactions’ decoupling could be used instead of a two-stage calculation. For details on the technical aspects of BAR calculations, please refer to other papers (Bennett, 1976; Christ et al., 2010; Pohorille et al., 2010; Villa & Mark, 2002).

Analysis of conformational states of n-pentadecane in a condensed phase was performed with an in-house Python program used to process output data from the GROMACS’s g_angle program. All the C-C-C-C dihedral angles of n-pentadecane were labeled as either gauche (+), gauche (–), trans or unclassified (with dihedral angles contained in the range of: 60 ± 30 degrees, –60 ± 30 degrees, 0 ± 30 degrees or beyond the aforementioned ranges, respectively). The content of gauche (G), end-gauche (GT-), double-gauche (-GG-) and kink (-GTG-) conformations was obtained by averaging over all appropriate dihedral angles of n-pentadecane molecules, with G designating either of the two gauche conformations. The obtained time series were averaged and their error estimated by the block averages method, using the g_analyze program from the GROMACS package. Figure 1 presents examples of the analyzed conformational states of n-pentadecane.

Kinetic properties. Calculations of kinetic properties (shear viscosity, η and self-diffusion coefficient, D) of n-pentadecane required MD simulations in the NVT ensemble. For each of the following FF parameter sets: OPLS-QQ, Slipids, BergerLips and CHARMM C36, we carried out three 300 ns MD simulations and used the trailing 100 ns parts of trajectories in calculations of the kinetic properties. We omitted the OPLS-AA FF set in reported calculations since at the temperature of 298.15K, which was kept in our MD simulations, n-pentadecane entered a solid state – not expected at standard conditions. Results for the OPLS-MP set are taken from our previous studies (Murzyn et al, 2013).

In calculations of shear viscosity, η, we have used the Green-Kubo formula (eqn. 2), as implemented in one of our Python in-house programs:

where Pαβ are the unique off-diagonal terms of the pressure tensor (Pxy, Pxz, Pyz, Pxx– Pyy, and Pyy– Pzz) extracted from trajectories by the GROMACS’s g_energy program program, V is the volume of the simulation box, kB is the Boltzmann constant, T is the absolute temperature, and 〈...〉 denotes averaging over consequtive time frames t0.

The self-diffusion coefficients were calculated using GROMACS’s g_msd program by fitting a linear function to the Einstein equation:

with ri being the instantaneous position of the center of mass of a molecule, df the number of translational degrees of freedom, and 〈...〉 denoting the average over all the molecules in the ensemble. The self-diffusion coefficients reported here were corrected for periodic boundary condition artifacts:

where η is a value of independently determined shear viscosity, L is the cubic box length, and ζ = 2.837298 is a constant resulting from the system’s symmetry (Yeh & Hummer, 2004). For more details regarding calculations of shear viscosity and self-diffusion coefficients, see our previous paper (Murzyn et al., 2013).

RESULTS AND DISCUSSION

We compared the performance of the selected FFs parameters in their ability to reproduce both, the condensed-phase and kinetic properties of n-pentadecane, using a series of MD simulations. Among condensed-phase properties of n-pentadecane we considered its melting point, density, enthalpy of vaporization, and Gibbs free energy of hydration. We also reported calculated values of shear viscosity and self-diffusion coefficients that being kinetic properties allowed us to evaluate how accurately various FFs described dynamics of liquid n-pentadecane. Finally, we summarize results of a detailed conformational analysis on n-pentadecane which allowed us to quantitatively describe the preferred conformations of n-pentadecane both, in a liquid and solid state.

Melting point

In Fig. 2 we show snapshots of the condensed-phase systems at the end of our NPT simulations. Presence of a large fraction of extended (all-trans) conformation was considered to be a hallmark of the solid phase. Although creation of uniform long-range crystalline ordering in the system is hindered by the system’s size and periodic boundaries, nucleation of highly ordered domains is an indicator of n-pentadecane tendency to solidify.

The experimentally determined melting point (Tm) for n-pentadecane is 283.1K (Lide, 2004). n-Pentadecane simulated with the OPLS-RAW parameters adopts a solid phase both at 277K and 298K. Thus, Tm is significantly overestimated. n-Pentadecane simulated with the OPLS-MP parameters is in a solid phase at 277K and in a liquid phase at 298K. This result is an expected outcome of the parametrization procedure: the torsional parameters in OPLS-MP were iteratively modified in order to reproduce the melting point for long alkanes (Murzyn et al., 2013). The assignment of CHARMM C36 parameters to n-pentadecane also leads to obtaining a correct melting point; however, in case of OPLS-QQ, Slipids and Berger FFs, a liquid state of n-pentadecane is observed both at 298K and 277K. In a series of additional short MD simulations with systematically lowered temperature, we estimated values of Tm to be approx. 266 ± 4K for both, the OPLS-QQ and Slipids sets, and 263 ± 4K for the BergerLips set. We included the determined values of Tm together with other properties of n-pentadecane simulated with various FF parameters’ sets in Table 2.

In Fig. 3, we show time evolution of the fraction of fully extended conformations of n-pentadecane simulated at 277K with OPLS-MP and CHARMM C36 sets. As can be seen, the equilibrium value of this property is reached within 20 ns (OPLS-MP) and 10 ns (CHARMM C36). Thus we can conclude that the trailing 50 ns part of MD trajectories can be safely used in further analyses of n-pentadecane in a thermally equilibrated solid state.

Density

Following criteria established in the large scale benchmark of OPLS-AA FF (Caleman et al., 2012), we consider values of density to be properly reproduced if they deviate from the experimental reference values by no more than 2%.

In the case of the OPLS-RAW set, the determined density of n-pentadecane (848.55 kg/m3) is markedly overestimated at 298K with respect to the experimental value (765.00 kg/m3, Tofts et al., 2000). This is expected since, as we have shown previously, the assignment of these FF parameters results in transition of n-pentadecane into a solid phase. OPLS-QQ, OPLS-MP, CHARMM C36, and Slipids sets give correct densities at 298K. Calculations of densities carried out with the Slipids set are the most accurate relatively to experimental data (+1.1%), followed by OPLS-QQ (–2.3%) and OPLS-MP (–2.6%). The density obtained with the CHARMM C36 parameters’ set is only 0.7% lower than the experimental value, but this value is burdened with a large error (45 kg/m3). It is noteworthy that only OPLS-RAW, OPLS-MP and CHARMM C36 sets correctly predict a solid phase for n-pentadecane at 277K. As evidenced in Fig. 2, simulations carried out with OPLS-QQ and Slipids retain a liquid state as far as 6K below the melting point. For the OPLS-QQ set, the density at 277K does not exceed the experimental density at the melting point (775.26 kg/m3, Eslami, 2000). Even though density calculated at 277K for the Slipids set is slightly higher, it does not exhibit a sharp transition observed for the OPLS-MP set. These facts further evidence that the OPLS-QQ and Slipids sets are unable to reproduce the solid state formation by n-pentadecane at 277K. Figure 4 shows the time evolution of densities at 277K and at 298K in simulations employing the studied FF parameters’ sets.

Enthalpy of Vaporization

Following the criteria established in the Caleman’s benchmark (Caleman et al., 2012), we consider values of the enthalpy of vaporization to be properly reproduced if they deviate from experimental reference values by no more than 11%. ΔHvap for the OPLS-RAW set is overestimated by 10.1 kcal/mol with respect to the experimental value (18.35 kcal/mol, Tofts et al., 2000, Table 2). Assignment of the OPLS-QQ and OPLS-MP sets results in ΔHvap of 18.44 and 18.45 kcal/mol, respectively. These values are in the best agreement with the reference value (+0.05%) among results obtained for all FF considered in this work. ΔHvap for other considered FF parameters’ sets are lower than the experimental value, while the assignment of the Slipids and CHARMM C36 sets results in ΔHvap being very close to the reference (–1.4% and –2.6%, respectively), while in the case of BergerLips we obtained a relatively large difference of 14.1%.

Gibbs Free Energy of Hydration

The calculated ΔGhyd are gathered in Table 2. The OPLS-RAW set yields ΔGhyd which differs by merely 0.31 kcal/mol from the reference of 4.13 kcal/mol (Ben-Naim & Marcus, 1984). The results obtained for the OPLS-MP set are in an even better agreement, but in this case we did not expect noticeable discrepancy since a sole change in parameters for a bonded term (i.e. dihedral potentials) seems to have a negligible effect on the interactions between water and a single molecule with no polar groups. On the other hand, doubling the partial charges on carbon atoms and the appropriate readjustment of hydrogen partial charges in the OPLS-QQ set resulted in an evidently lower ΔGhyd, which is presumably an effect of increased polarity of the n-pentadecane molecule. This result remains in the acceptable range, as it is only 0.9 kcal/mol lower than the relevant theoretical prediction. However, according to our calculations, ΔGhyd for Slipids is distinctly overestimated (87% over the theoretical reference value). It is worth noting that the ΔGhyd values for model compounds were not reported in the original Slipids papers (Jämbeck & Lyubartsev, 2012; Jämbeck & Lyubartsev, 2013). In the latter article (Jämbeck & Lyubartsev, 2013), the Slipids developers report detailed results for six nonpolar organic compounds, both in the polar and nonpolar solutions. Unfortunately, none of their model compounds could be used as a reference for long aliphatic chains in phospholipid molecules. Finally, the calculated values of ΔGhyd for the CHARMM C36 and BergerLips sets are similar to each other and overestimate ΔGhyd by approximately 38% and 27%, respectively.

Shear Viscosity and Diffusion

The reproduction of shear viscosity in the all-atom force fields is generally satisfactory, while the results obtained for the BergerLips set are notably different (1.10 mPa • s, Table 2) from the experimentally determined value (2.54 mPa • s, Tofts et al., 2000). There is a visible correspondence between the results obtained for the OPLS-MP and OPLS-QQ, which is not surprising since both of these FF sets share the same origin. On the other hand, results obtained with the CHARMM C36 and Slipids sets reveal a bigger difference, despite the fact that the substantial part of Slipids parameters was taken from CHARMM C36. In Fig. 5, we present the convergence of the integrals of Green-Kubo formula for different FF sets emerging from three independent MD simulations carried out in each case.

For the OPLS-MP set, reproduction of shear viscosity is the best among all compared FFs (+5.9% when compared to the reference). The shear viscosity obtained for other FFs, i.e. the OPLS-QQ (+13.8%), Slipids (+9.1%), and CHARMM C36 (–8.3%) sets are within an acceptable range.

The self-diffusion coefficients are in all cases somewhat underestimated. The largest discrepancies in reproduction of both considered kinetic properties of n-pentadecane are observed for the BergerLips set: this united-atom FF overestimates the diffusibility of n-pentadecane and underestimates the shear viscosity roughly by a factor of two. The data of MD simulation of bulk n-pentadecane in BergerLips are relatively scarce, as opposed to simulations of lipid bilayers, but it is worth noting that according to Piggot and co-workers (Piggot et al., 2012), the lateral diffusion coefficient of DPPC bilayer with BergerLips FF parameters is overestimated by a similar factor of two, which confirms a problem with reproduction of this kinetic property.

Chain Conformations

When reporting Tm values for n-pentadecane modelled by different FFs, we pointed out that the total number of gauche conformations in the molecule is incorrectly reproduced by the OPLS-RAW set at 298K (0.19 gauche conformations of C-C-C-C dihedrals per molecule) and it reveals no difference when determined at 277K (i.e. 0.15 gauche). While incorrect for the liquid phase, such low value is expected for the solid state in these simulations (cf. Fig. 2). For all other investigated FF sets, the total number of gauche conformations per molecule at 298K is close to the theoretical estimate by Holler and Callis (Holler et al., 1989). However, only for the OPLS-MP and CHARMM C36 sets the gauche content at 277K is low enough (1.16 and 1.06, respectively) to trigger formation of the solid phase. A closer look at Fig. 2 suggests that a mixture of a highly ordered solid-like phase and a liquid phase is present in the captured state of simulations. If this is the case, the number of gauche conformations at 277K would be an average of both, the liquid and solid fractions of n-pentadecane, with the fully ordered phase having gauche content close to zero and the liquid disordered fraction having the average gauche content of approximately four per molecule. In the simulations employing the OPLS-QQ, Slipids, and BergerLips sets, the total gauche content is virtually identical at both temperatures, consistently showing that no phase transition is taking place (see Table 3).

For all considered parameter sets that reproduce a liquid state of n-pentadecane at 298K, the content of gauche, end-gauche, double-gauche and kink conformations agree well with both, the relevant theoretical predictions and experimental results. However, due to the sparsity of experimental results available for 298K, quantitative comparisons with the calculated numbers of different conformations are not possible. There is a visible difference between the content of extended, gauche, end-gauche, double-gauche and kink conformations at 277K and 298K only for CHARMM C36 and OPLS-MP (Table 3), which is a consequence of different phase states at considered temperatures. It has been proven that different solid phases of long (C19 and longer) n-alkanes contain specific, varying amounts of the nonplanar defects described above (Marnocelli, 1982), so the approach described herein could be extended to investigate reproduction of this phase behavior in appropriate MD simulations, provided that the employed FF is suitable for long hydrocarbons.

CONCLUSION

Following the concept of graphical comparison of abilities of different FF parameters’ sets to reproduce the reference values, which was introduced by Botan and co-workers (Botan et al., 2015), we show in Fig. 6 performance of the OPLS-RAW, OPLS-MP, OPLS-QQ, BergerLips, Slipids, and CHARMM C36 sets. For the purpose of quantitative comparison of n-pentadecane chain conformations, we considered the subset of conformations from Table 3, which permits direct comparison with results obtained experimentally (i.e. the sum of DG, EG and kink conformations).

Summarizing all the results obtained in this work for n-pentadecane, we can conclude that assignment of the OPLS-MP parameters resulted in a very good agreement between all calculated condensed-phase properties (i.e. density, enthalpy of vaporization, conformation populations, Gibbs free energy of hydration and shear viscosity) and relevant experimental data (see Table 2). The results obtained with the CHARMM C36 set agree well with many of the considered properties. Our study indicates a weak point of employing the Slipids FF parameters, as they are unable to reproduce the melting point temperature of n-pentadecane and overestimate ΔGhyd approximately by a factor of two. In many aspects, the OPLS-QQ set increases reliability of modelling a liquid state of n-pentadecane as compared to the original OPLS-AA parameters, but nevertheless it shows some deficiencies, especially in terms of providing correct predictions of Tm . Moreover, assignment of the OPLS-QQ FF results in ΔGhyd and shear viscosity values characterized by a substantially larger error with respect to the expected value than the OPLS-MP set. The accuracy of simulations performed with BergerLips FF was the least satisfactory among all tested parameter sets, with all of the calculated properties drastically differing from the experimental data available for n-pentadecane.

REFERENCES

Ben-Naim A, Marcus Y (1984) Solvation thermodynamics of nonionic solutes. J. Chem. Phys. 81: 2016–2027. https://doi.org/10.1063/1.447824

Bennett CH (1976) Efficient estimation of free energy differences from Monte Carlo data. J. Comput. Phys. 22: 245–268. https://doi.org/10.1016/0021-9991(76)90078-4

Berendsen HJC, Postma JPM, Gunsteren WF, Hermans J (1981) Interaction models for water in relation to protein hydration. Intermolecular Forces 331–342. https://doi.org/10.1007/978-94-015-7658-1_21

Berger O, Edholm O, Jahnig F (1997) Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 72: 2002–2013. https://doi.org/10.1016/S0006-3495(97)78845-3

Botan A, Favela-Rosales F, Fuchs PFJ, Javanainen M, Kanduc M, Kulig W, Lamberg A, Loison C, Lyubartsev A, Miettinen MS, Monticelli L, Maättä J, Ollila OHS, Retegan M, Róg T, Santuz H, Tynkkynen J (2015) Toward atomistic resolution structure of phosphatidylcholine headgroup and glycerol backbone at different ambient conditions. J. Phys. Chem. B 119: 15075–15088. https://doi.org/10.1021/acs.jpcb.5b04878

Caleman C, Maaren P, Hong M, Hub JS, Costa LT, Spoel D (2012) Force field benchmark of organic liquids: density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant. J. Chem. Theory. Comput. 8: 61–74. https://doi.org/10.1021/ct200731v

Christ CD, Mark AE, Gunsteren WF (2010) Basic ingredients of free energy calculations: a review. J. Comput. Chem. 31: 1569–1528.
https://doi.org/10.1002/jcc.21450

Damm W, Frontera A, Tirado-Rives J, Jorgensen WL (1997) OPLS all-atom force field for carbohydrates. J. Phys. Chem. 18: 1955–1970. https://doi.org/10.1002/(SICI)1096-987X(199712)18:16<1955::AID-JCC1>3.0.CO;2-L

Dickson CJ, Madej BD, Skjevik ÅA, Betz RM, Teigen K, Gould IR, Walker RC (2014). Lipid14: the amber lipid force field. J. Chem. Theory Comput. 10: 865–879. https://doi.org/10.1021/ct4010307

Eslami H (2000) Equation of state for long-chain n-alkanes. Fluid Phase Equilib. 169: 19–30. https://doi.org/10.1016/S0378-3812(99)00331-3

Ermilova I, Lyubartsev AP (2016) Extension of the slipids force field to polyunsaturated lipids. J. Phys. Chem. B. 120: 12826–12842.
https://doi.org/10.1021/acs.jpcb.6b05422

Flory PJ (1969) Statistical mechanics of chain molecules (Interscience Publishers). https://doi.org/10.1002/bip.1969.360080514

Haynes WM (2011) CRC Handbook of Chemistry and Physics, 91st edn, Taylor and Francis Group, LLC, Boca Raton FL

Hess B, Kutzner C, Spoel D, Lindahl E (2008) Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 4: 435–447. https://doi.org/10.1021/ct700301q

Hess B (2002) Determining the shear viscosity of model liquids from molecular dynamics simulations. J. Phys. Chem. B. 116: 209–217.
https://doi.org/10.1063/1.1421362

Holler F, Callis JB (1989) Conformation of the hydrocarbon chains of sodium dodecyl sulfate molecules in micelles: an FTIR study. J. Phys. Chem. 93: 2053–2058. https://doi.org/10.1021/j100342a068

Hoover WG (1985) Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31: 1695–1697. https://doi.org/10.1103/PhysRevA.31.1695

Jämbeck JPM, Lyubartsev AP (2012) Derivation and systematic validation of a refined all-atom force field for phosphatidylcholine lipids. J. Phys. Chem. B. 116: 3164–3179. https://doi.org/10.1021/jp212503e

Jämbeck JPM, Lyubartsev AP (2012) An extension and further validation of an all-atomistic force field for biological membranes. J. Chem. Theory Comput. 8: 2938–2948. https://doi.org/10.1021/ct300342n

Jämbeck JPM, Lyubartsev AP (2013) Another piece of the membrane puzzle: extending slipids further. J. Chem. Theory Comput. 9: 774–784. https://doi.org/10.1021/ct300777p

Jämbeck JPM, Lyubartsev AP (2013) Implicit inclusion of atomic polarization in modeling of partitioning between water and lipid bilayers. Phys. Chem. Chem. Phys. 15: 4677–4686. https://doi.org/10.1039/C3CP44472D

Jorgensen W, Madura J, Swenson C (1984) Optimized intermolecular potential functions for liquid hydrocarbons. J. Am. Chem. Soc. 106: 6638-6646. https://doi.org/10.1021/ja00334a030

Jorgensen WL, Maxwell DS, Tirado-Rives J (1996) Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118: 11225–11236. https://doi.org/10.1021/ja9621760

Jorgensen WL (1981) Quantum and statistical mechanical studies of liquids. 10. Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. J. Am. Chem. Soc. 103: 335–340. https://doi.org/10.1021/ja00392a016

Kahn K, Bruice TC (2002) Parameterization of OPLS-AA force field for the conformational analysis of macrocyclic polyketides. Comput. Chem. 23: 977–996. https://doi.org/10.1002/jcc.10051

Kaminski GA, Friesner RA, Tirado-Rives J, Jorgensen WL (2001) Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. J. Phys. Chem. B. 105: 6474–6487. https://doi.org/10.1021/jp003919d

Klauda JB, Venable RM, Freites JA, O’Connor JW, Tobias DJ, Mondragon-Ramirez C, Vorobyov I, MacKerell AD, Pastor RW (2010) Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J. Phys. Chem. B. 114: 7830–7843.
https://doi.org/ 10.1021/jp101759q

Klauda JB, Monje V, Kim T, Im W (2012) Improving the CHARMM force field for polyunsaturated fatty acid chains. J. Phys. Chem. B. 116: 9424–9431. https://doi.org/10.1021/jp304056p

Kony D, Damm W, Stoll S, Van Gunsteren WF (2002) An improved OPLS–AA force field for carbohydrates. J. Phys. Chem. 23: 1416–1429. https://doi.org/10.1002/jcc.10139

Kirschner KN, Lins RD, Maass A, Soares TA (2012) A glycam-based force field for simulations of lipopolysaccharide membranes: parametrization and validation. J. Chem. Theory Comput. 8: 4719–4731. https://doi.org/10.1021/ct300534j

Leekumjorn S, Sum AK (2007) Molecular studies of the gel to liquid-crystalline phase transition for fully hydrated DPPC and DPPE bilayers. Biochim. Biophys. Acta 1768: 354–365. https://doi.org/10.1016/j.bbamem.2006.11.003

Lemkul J (2012) Gromas tutorial, kalp15 in dppc . URL http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/membrane_protein/index.html

Lide D (2004) CRC handbook of chemistry and physics : a ready-reference book of chemical and physical data (CRC Press, Boca Raton, Fla). https://doi.org/10.1021/ja041017a

Lim JB, Rogaski B, Klauda JB (2012) Update of the Cholesterol Force Field Parameters in CHARMM. J. Phys. Chem. B. 116: 203–210.
https://doi.org/10.1021/jp207925m

Maciejewski A, Pasenkiewicz-Gierula M, Cramariuc O, Vattulainen I, Rog T (2014) Refined OPLS all-atom force field for saturated phosphatidylcholine bilayers at full hydration. J. Phys. Chem. B. 118: 4571–4581. https://doi.org/10.1021/jp5016627

Maier JA, Martinez C, Kasavajhala K, Wickstrom L, Hauser KE, Simmerling C (2015) ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11: 3696–3713. https://doi.org/10.1021/acs.jctc.5b00255

Maroncelli M, Qi SP, Strauss HL, Snyder RG (1982) Nonplanar conformers and the phase behavior of solid n-alkanes. J. Am. Chem. Soc. 104: 6237–6247. https://doi.org/10.1021/ja00387a013

Marsh D (1991) General features of phospholipid phase transitions. Chem. Phys. Lipids 57: 109–120. https://doi.org/10.1016/0009-3084(91)90071-I

McDonald NA, Jorgensen WL (1998). Development of an all-atom force field for heterocycles. properties of liquid pyrrole, furan, diazoles, and oxazoles. J. Phys. Chem. B. 102: 8049–8059. https://doi.org/10.1021/jp981200o

Murzyn K, Bratek M, Pasenkiewicz-Gierula M (2013) Refined OPLS all-atom force field parameters for n-pentadecane, methyl acetate, and dimethyl phosphate. J. Phys. Chem. B. 117: 16388–16396.
https://doi.org/10.1021/jp408162d

Murzyn K, Pasenkiewicz-Gierula M (2015) Structural properties of the water/membrane interface of a bilayer built of the E. coli lipid A. J. Phys. Chem. B. 119: 5846–5856. https://doi.org/10.1021/jp5119629

Neria E, Fischer S, Karplus M (1996) Simulation of activation free energies in molecular systems. J. Phys. Chem. 105: 1902–1921. https://doi/org/10.1063/1.472061

Nose S (1984) A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 81: 511–519. https://doi.org/10.1063/1.447334

Parrinello M, Rahman A (1981) Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 52: 7182. https://doi.org/10.1063/1.328693

Pasenkiewicz-Gierula M, Baczyński K, Murzyn K, Markiewicz M (2012) Orientation of lutein in a lipid bilayer – revisited. Acta Biochim. Pol. 59: 115–118. https://doi.org/10.18388/abp.2012_2184

Pastor RW, MacKerell AD Jr (2011) Development of the CHARMM force field for lipids. J. Phys. Chem. Lett. 2: 1526–1532. https://doi.org/10.1021/jz200167

Piggot TJ, Pineiro A, Khalid S (2012) Molecular dynamics simulations of phosphatidylcholine membranes: a comparative force field study. J. Chem. Theory Comput. 8: 4593–4609. https://doi.org/10.1021/ct3003157

Plesnar E, Subczynski WK, Pasenkiewicz-Gierula M (2012) Saturation with cholesterol increases vertical order and smoothes the surface of the phosphatidylcholine bilayer: A molecular simulation study. Biochimica et Biophysica Acta (BBA) – Biomembranes 1818: 520–529. https://doi.org/10.1016/j.bbamem.2011.10.023

Plesnar E, Subczynski WK, Pasenkiewicz-Gierula M (2013) Comparative computer simulation study of cholesterol in hydrated unary and binary lipid bilayers and in an anhydrous crystal. J. Phys. Chem. B. 117: 8758–8769. https://doi.org/10.1021/jp402839r

Pohorille A, Jarzynski C, Chipot C (2010) Good practices in free-energy calculations. J. Phys. Chem. B. 114: 10235–10253. https://doi.org/10.1021/jp102971x

Pranata J, Wierschke SG, Jorgensen WL (1991) OPLS potential functions for nucleotide bases. Relative association constants of hydrogen-bonded base pairs in chloroform. J. Am. Chem. Soc. 113: 2810–2819. https://doi.org/10.1021/ja00008a002

Rizzo RC, Jorgensen WL (1999) OPLS all-atom model for amines: resolution of the amine hydration problem. J. Am. Chem. Soc. 121: 4827–4836. https://doi.org/10.1021/ja984106u

Ryckaert JP, Bellemans A (1975) Molecular dynamics of liquid n-butane near its boiling point. Chem. Phys. Lett. 30: 123–125. https://doi.org/10.1016/0009-2614(75)85513-8

Ryckaert JP, Bellemans A (1978) Molecular dynamics of liquid alkanes. Faraday Discuss. 66: 95–106. https://doi.org/10.1039/DC99191FX001

Senak L, Davies MA, Mendelsohn R (1991) A quantitative IR study of hydrocarbon chain conformation in alkanes and phospholipids: CH2 wagging modes in disordered bilayer and HII phases. J. Phys. Chem. 95: 2565–2571. https://doi.org/10.1021/j100159a084

Siu SWI, Pluhackova K, Böckmann RA (2012) Optimization of the OPLS-AA force field for long hydrocarbons. J. Chem. Theory Comput. 8: 1459–1470. https://doi.org/10.1021/ct200908r

Skjevik ÅA, Madej BD, Walker RC, Teigen K (2012). LIPID11: a modular framework for lipid simulations using amber. J. Phys. Chem. B. 116: 11124–11136. https://doi.org/10.1021/jp3059992

Skjevik ÅA, Madej BD, Dickson CJ, Teigen K, Walker RC, Gould IR (2015) All-atom lipid bilayer self-assembly with the AMBER and CHARMM lipid force fields. Chem. Commun. 51: 4402–4405.
https://doi.org/10.1039/c4cc09584g

Stepniewski M, Bunker A, Pasenkiewicz-Gierula M, Karttunen M, Rog T (2010) Effects of the lipid bilayer phase state on the water membrane interface. J. Phys. Chem. B. 114: 11784–11792. https://doi.org/10.1021/jp104739a

Tofts PS, Lloyd D, Clark CA, Barker GJ, Parker GJM, McConville P, Baldock C, Pope JM (2000) Test liquids for quantitative MRI measurements of self-diffusion coefficient in vivo. Magn. Reson. Med. 43: 368–374. https://doi.org/10.1002/(sici)1522-2594(200003)43:3%3C368::aid-mrm8%3E3.0.co;2-b

Tuteja A, Choi W, Ma M, Mabry JM, Mazzella SA, Rutledge GC, McKinley GH, Cohen RE (2007) Designing superoleophobic surfaces. Science 318: 1618–1622. https://doi.org/10.1126/science.1148326

Van Buuren AR, Marrink SJ, Berendsen HJC (1993) A molecular dynamics study of the decane/water interface. J. Phys. Chem. 97: 9206–9212. https://doi.org/10.1021/j100138a023

Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJC (2005) GROMACS: fast, flexible, and free. J. Comput. Chem. 26: 1701–1718. https://doi/org/10.1002/jcc.20291

Van Gunsteren WF, Berendsen HJC (1987) Groningen Molecular Simulation (GROMOS) Library Manual, Biomos, Groningen, The Netherlands, pp 1–221.

Vereshchagin AN, Vul’fson SG (1967) Bulletin of the Academy of Sciences of the USSR, Division of Chemical Science 16: 1186–1189.

Villa A, Mark AE (2002) Calculation of the free energy of solvation for neutral analogs of amino acid side chains. J. Comput. Chem. 23: 548–553. https://doi.org/10.1002/jcc.10052

Yeh IC, Hummer G (2004) System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with Periodic Boundary Conditions. J. Phys. Chem. B. 108: 15873–15879. https://doi.org/10.1021/jp0477147