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

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 ad-dress 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. no conflicts of interest for this article. Abbreviations : MD, molecular dynamics; FF, force field; PC, phosphatidylcholine; OPLS-AA, Optimized Potentials for Liquid Simu- lations All-Atom; Slipids, Stokholm lipids’ force field; BergerLips, Berger lipids’ force field; PE, phosphatidylethanolamine; SM, sphin- gomyelin; PS, phosphatidylserine; PG, phosphatidylglycerol; DPPC, dipalmitoylphosphatidylcholine; POPE, 1-palmitoyl-2-oleoylphos- phatidylethanolamine; POPC, 1-palmitoyl-2-oleoylphosphatidylcho-line; 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) , 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 Lyubartsev et al., 2015).
Slipids parameters are a modern addition to the lipid force field family. Initially, they were developed for saturated phospholipids  and then extended to include unsaturated phosphatidylcholine (PC) and phosphatidylethanolamine (PE) , as well as sphingomyelin (SM), phosphatidylglycerol (PG), phosphatidylserine (PS), and cholesterol 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 , 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 npentadecane 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 , 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, npentadecane (C 15 H 32 ) may be used as a model of the palmitoyl chain, present in some of the most popular biomembrane 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 ap-plicable at non-physiological temperatures.  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.
The general strategy for adapting parameters of npentadecane 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 CON-TRIBUTIONS). 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 Properties of n-pentadecane emerging from application of biomolecular force fields (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 coworkers 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 (ΔH vap ) was calculated according to the following equation: where E gas and E liq are potential energies in the gas and liquid state, respectively, R is the gas constant, and T is the absolute temperature. Additional gas phase simula-tions 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).
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 ΔG hyd 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 Table 1. Simulation parameters for FFs used in this study. OPLS-AA herein jointly designates the OPLS-RAW, OPLS-QQ, and OPLS-MP sets. PME grid spacing and all of the cutoffs are given in nm. E+P denotes dispersion correction applied to both, the energy and pressure. Footnotes: (a) A switch function was used between 1.4 and 1.5 nm to ensure a smooth decay to zero. (b) Parameters based on the optimal simulation protocols recommended by Piggot and coworkers (Piggot et al., 2012). (c) A switch function was used between 0.8 and 1.2 nm. (d) Parameters recommended by the CHARMM C36 authors (Klauda et al., 2010)  (1) pot pot 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 (P xy , P xz , P yz , P xx -P yy , and P yy -P zz ) extracted from trajectories by the GROMACS's g_energy program program, V is the volume of the simulation box, k B is the Boltzmann constant, T is the absolute temperature, and 〈...〉 denotes averaging over consequtive time frames t 0 .
The self-diffusion coefficients were calculated using GROMACS's g_msd program by fitting a linear function to the Einstein equation: with r i being the instantaneous position of the center of mass of a molecule, d f 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 condensedphase 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. (2) (4) Properties of n-pentadecane emerging from application of biomolecular force fields

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 (T m ) 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, T m 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 T m 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 T m 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 simu- lated 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/m 3 ) is markedly overestimated at 298K with respect to the experimental value (765.00 kg/m 3 , 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/m 3 ). 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/m 3 , 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%. ΔH vap 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 ΔH vap 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. ΔH vap for other considered FF parameters' sets are lower than the experimental value, while the assignment of the Slipids and CHARMM C36 sets results in ΔH vap 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 ΔG hyd are gathered in Table 2. The OPLS-RAW set yields ΔG hyd 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 ΔG hyd , 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, ΔG hyd for Slipids is distinctly overestimated (87% over the theoretical reference value). It is worth noting that the ΔG hyd values for model compounds were not reported in the original Slipids pa- We considered a conformation of n-pentadecane molecule to be a fully extended one if all but its terminal C-C-C-C dihedral angles are in the trans conformation, i.e. 180o ± 30o. Time profiles for OPLS-MP set is shown in black and for CHARMM C36 in gray.  . In the latter article , 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 ΔG hyd for the CHARMM C36 and BergerLips sets are similar to each other and overestimate ΔG hyd 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.
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 unitedatom 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 T m 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 Table 2. Calculated physicochemical properties of n-pentadecane simulated with various FF parameters' sets. Densities (d) are given in kg/m 3 , enthalpies of vaporization (ΔH vap ) and Gibbs free energy of hydration (ΔG hyd ) in kcal/mol, T m (melting point) in K, self-diffusion coefficients (D) in cm 2 /μs and shear viscosity (η) in mPa • s. Unless otherwise noted, quantities were calculated at 298K. Errors of d and ΔH vap (in parentheses) were estimated by the GROMACS's block averages method, implemented in the g_analyze program. Errors of ΔG hyd were calculated by the GROMACS's g_bar program. Errors of kinetic properties are Standard Errors of the Mean (SEMs, in parentheses) and were calculated as arithmetic averages from results of three independent MD simulation runs. Further notes: (a) kinetic properties for OPLS-RAW were not calculated due to the forming of solid-like phase at 298K, (b) reported previously (Murzyn et al., 2013), (c) experimental values (Lide, 2004), (d) a value extrapolated from the linear dependence of enthalpy of hydration on the number of carbon atoms (Ben-Naim & Marcus, 1984), (e) experimental values (Tofts et al., 2000).  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 solidlike 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 (C 19 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 ΔG hyd 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 T m . Moreover, assignment of the OPLS-QQ Table 3. Statistics of extended (E), gauche (G), end-gauche (EG), double-gauche (DG) and kink (kink) conformations per one molecule in condensed phase simulations. In case of columns dedicated to the Extended (E) conformation, i.e. E(298K) and E(277K) -the values represent the percentage fraction of molecules that exhibits the extended conformation, averaged over simulation time. For the remaining conformations, i.e. G, EG, DG, kink -the values correspond to the average content of the conformation per molecule, and thus they are averaged over all of the C-C-C-C torsional angles of the n-pentadecane molecule, all the molecules in the ensemble, and the production simulation time. The error estimates are in parentheses. If omitted, they are less or equal to 0.01. Additional notes: (a) An estimate for n-pentadecane based on the RIS (Rotational Isomeric State) model (Flory, 1969) -36% gauche, by Holler and Callis (Holler et al., 1989), (b) 0-10% gauche in the gel phase phospholipids, estimated by Marsh (Marsh, 1991), (c) Theoretical values for n-pentadecane, based on the RIS model by Senak and co-workers (Senak et al., 1991), (d) Experimental values for n-tridecane (C 13 H 28 ) in 298K, by Holler and Callis (Holler et al., 1989)  FF results in ΔG hyd 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. Figure 6. Performance of the considered FF parameters' sets.
The ability of a given FF parameters' set to reproduce the enthalpy of vaporization (ΔH), density (d), Gibbs free energy of hydration (ΔG), shear viscosity (η), diffusion self-coefficient (D), number of chain conformation defects (C), and melting point (T m ) of n-pentadecane is color-coded with GREEN, ORANGE, and RED corresponding to, excellent, good, and acceptable agreement, respectively, between the calculated and the reference value. BLACK denotes failure in qualitative reproduction of an appropriate reference value. Different threshold values were defined for each of the compared properties of n-pentadecane. For d, GREEN, ORANGE, RED, and BLACK correspond to the relative error of less than 2%, between 2% and 5%, between 5% and 10%, and 10% or more, respectively. For ΔH, ΔG, D, η, and these threshold values were defined as: <10%, 10-20%, 20-40%, and ≥40%. When evaluating reproduction of T m , we used the absolute errors of 5K, 5-10K, 10-20K, 20K or more. The Σ rows/columns contain values of a scoring function, which correspond to the overall performance of a given FF parameters' set (the lower the better) or the relative difficulty in reproducing a given physicochemical parameter (the higher the more difficult). The scores were calculated by assigning 0, 1, 3, and 6 pts, to GREEN, ORANGE, RED, and BLACK categories, respectively.