Common Force Field Thermodynamics of Cholesterol

Four different force fields are examined for dynamic characteristics using cholesterol as a case study. The extent to which various types of internal degrees of freedom become thermodynamically relevant is evaluated by means of principal component analysis. More complex degrees of freedom (angle bending, dihedral rotations) show a trend towards force field independence. Moreover, charge assignments for membrane-embedded compounds are revealed to be critical with significant impact on biological reasoning.


Introduction
Apart from a prominent role in cardiovascular disease [1,2], cholesterol is also known for its largely varying membrane concentration in different tissues. For example, red blood cells exhibit concentrations of 50 mol%, significantly higher than the average 22 mol% in regular vertebrate cell membranes [3]. Even higher levels are seen in brain and nerve tissue, the metabolic significance of which is yet to be determined [4,5]. Cholesterol is also an essential metabolic precursor for the biosynthesis of bile acids [6], steroid hormones [7], and vitamin D [8].
Owing to its importance, cholesterol has been the subject of intensive research over the last 70 years [9,10]. Its complicated cellular location has led to a recent intensification of auxiliary methods with many of them entirely based on computer simulation. Indeed, many instructive examples have been given and the underlying computational models have advanced to a reasonable level of quality [11][12][13][14][15][16][17].
A general requirement for all force fields is that different parameter sets applied to the same molecule-cholesterol, for example-should yield comparable results and should not reveal significant differences. From the many comparisons available, those focussing on dynamic aspects appear to be particularly interesting [26][27][28][29]. Here we want to add to this type of dynamic assessment and present data for the thermodynamics of cholesterol using four different force fields from the above list. Applying principal component analysis (PCA) [30] we identify the most "influential" bonds/angles/dihedrals in a particular parameter set and compare them to each other, thus emphasizing the dynamic 2 The Scientific World Journal character of cholesterol when described by common force fields.

AMBER(RESP).
A single cholesterol structure was model-built and optimized using Gaussian-09 [31] at the AM1 and B3LYP/3-21 G * level. The minimized geometry was considered at the HF/6-31 G * level and electrostatic potentials (ESPs) were computed for subsequent RESP assignment of atomic partial charges using ANTECHAMBER together with the GAFF force field (AMBER [18,32] version 11). A single copy of RESP/GAFF-parameterized cholesterol was loaded into XLEAP and output in appropriate AMBER formats (prmtop/inpcrd), then minimized (2500 steps, cutoff 20Å), and heated to 300 K target temperature (vacuum conditions) using Langevin dynamics, collision frequency = 1 ps −1 , 12Å cutoff, a time step of 1 fs, and no constraints on any bonds (i.e., SHAKE switched off). Identical conditions were applied during 5 ns of MD simulation where conformational snapshots were saved every 5000 steps to create a sample of 250 structures.

AMBER(bcc).
Cholesterol was model-built and optimized at the HF/6-31 G * * level using Gaussian-03 [31]. Atomic partial charges were assigned following the approach of AM1-bond charge correction (bcc) available in ANTECHAMBER of the AMBER package [18,32] (version 8). The optimized structure was minimized (20 steps, steepest descent) and heated to the target temperature of 300 K within 100 ps of equilibration MD. Production MD over 5 ns used a time step of 1 fs, SHAKE constraints on XH bonds, a Berendsen thermostat, vacuum conditions without periodicity, and AMBER version 10.

CHARMM.
Cholesterol parameters were employed as reported previously [33]. A system containing a single copy of cholesterol was set up and heated to 300 K based on straight dynamics and CHARMM36 all-hydrogen lipid topology/CHARMM27 all-hydrogen lipid parameters [19]. Production MD was extended over a period of 5 ns using a time step of 1 fs, no SHAKE constraints, and the default cutoff of 12Å specified in the cholesterol parameter file.

GROMACS.
Cholesterol parameters were used as reported previously [34]. A single molecule of cholesterol was put into a cubic box (58.219Å 3 ) and simulated at constant volume using periodic boundary conditions. All bonds were constrained [35], a time step of 2 fs was applied, total translation/rotation was periodically removed every 1000 steps, neighbour lists were updated every 5 steps, a cutoff of 10Å was used, and the system was maintained at 300 K target temperature by means of a Nose-Hoover thermostat.
2.5. PCA. All 250 snapshots collected by all the 4 different MD simulations were automatically converted from pdb format to xyz format and analyzed frame by frame for geometrical relationships (bonds/angles/dihedrals) with the help of TINKER [25] (module ANALYZE, option d) using the MM3 force field [36]. However, the MM3 force field was only a technical requirement for monitoring the actual dimensions of bonds/angles/dihedrals of a particular frame; hence these parameters were never really used or did never influence or alter any of the original geometries obtained by the 4 different force fields examined. PCA was then performed following the standard protocol [30]. A stand-alone solution was implemented in ANSI C supported by LAPACK [37].

ESPs (Membrane Internal).
Mean structures were computed via superposition of all the 250 snapshots of all the 4 force fields on arbitrarily chosen reference frames. Massweighted all-atom fitting including H-atoms was performed using TINKER [25] (module SUPERPOSE) based on dummy employment of the MM3 force field [36] (see note in the previous section). The "most representative" structure for each of the 4 force fields was then determined as that frame that showed minimum RMSD from the average structure. A continuum description of the cellular membrane containing cholesterol was employed following the experimental findings of Ashcroft, Subczynski, and coworkers [38,39].
Thus the OH-group and a small number of adjacent atoms on the cholesterol ring system (positions 2, 3, and 4 in Figure 1(a)) were assigned to the polar head group domain of the membrane which was represented by methanol [38]. The remaining part of cholesterol was assigned to the hydrophobic core domain of the membrane which was represented by cyclohexane [38]. Program POLCH [40] was used to compute ESPs following previous reports [15,16,41].

Results and Discussion
The standard nomenclature is adopted as schematically illustrated in Figure 1(a). Numbers are assigned to C-atoms only; hence missing H-atoms need to be considered whenever implicated in any of the identified bonds/angles/dihedrals. We started out to construct 4 different data sets of cholesterol conformations composed of 250 snapshots obtained from 5 ns of MD simulation based on AMBER(RESP) [18], AMBER(bcc) [32], CHARMM [19], and GROMACS [22]. The initial two descriptions differ with respect to the charge model applied. Either bonds (77 in total), angles (157), or dihedrals (259) were extracted from each of these 250 structures and written into separate data matrices which then became subject to PCA [30]. The top-ranked principal components, that is, specific linear combinations of the 77 bonds (or the other variables examined), will then identify those bonds (or the other variables examined), that experience the largest fluctuations, and hence are most relevant to the thermodynamics. Due to only marginal separation of eigenvalues corresponding to the top-ranked PCs, we took into account a subset, , of PCs capable of reestablishing 90% of the original data set. We then isolated the top 10% components as follows, The Scientific World Journal where refers to a particular bond/angle/dihedral, , part of the top 10% ranked coefficients, , of the PCs selected and weighted by their corresponding eigenvalue, . PCA results for bonds are shown in Figure 1(b) where different colours symbolize different force fields. Independence of parameter sets is seen when a particular bond is identified several times; for example, the bond between C 25 and C 26 was top-ranked both by CHARMM and GROMACS (hence its colour of half blue, half green). The model most prone to bond variations was AMBER(bcc) involving exclusively C-C bonds (see yellow substructure in Figure 1(b)). In contrast, only C-H bonds were ranked by AMBER(RESP) (see red numbers in Figure 1(b)). However, this just reflects the employment/avoidance of SHAKE [42] constraints (see computational methods). CHARMM did identify a small set of both types, while GROMACS revealed only a minor group of bonds at the terminal end of cholesterol, a consequence of the restraints on all bonds [35] imposed during MD. In general, most of the bonds forming the ring system were not implicated in any of the top lists, an indication that the tetracyclic ring system maintains a rather rigid structure.
Next we turned our attention to the PCA of angles, and corresponding results are summarized in Figure 1(c). Again, the ring system was characterized as rather rigid, and many thermodynamically relevant angles did involve H-C-H groups with even exclusive preference by AMBER(RESP) (see red numbers in Figure 1(c)). Several hot spots of angle variation were identified to be located at terminal CH 3groups (see, e.g., C 19 and C 27 ). The overall impression gained was that there is a trend towards force field independence with increasing complexity of the type of motion analyzed. For example, the number of multiply identified sites of topranked angles did increase considerably when compared to the number of top-ranked bonds indicated by more than one force field. Even more impressive in this respect was the PCA of dihedrals (Figure 1(d)). Virtually all sites were marked by all force fields with the exception of only isolated positions inside the ring system.
In an attempt to link up our PCA data with MM energies we next determined trends of kinetic and potential energies and relative contributions to the nonbonded energy of the three types of motion studied. Results are summarized in Figure 2. The first interesting finding was that all the 4 force fields delivered net potential energies of positive sign but to a different relative extent. The two AMBER descriptions were roughly comparable, while CHARMM tended to an equipartition of energy between kinetic and potential forms, and GROMACS put more emphasis on potential energy (see relative contributions of red and blue graphs in Figures  2(a)-2(d)). In addition, significant differences were seen with respect to the degrees of fluctuations affecting kinetic and potential energies. Here, CHARMM revealed the trajectory progressing most smoothly, while AMBER and GROMACS exhibited comparable levels of "thermal noise. " Focussing on the relative contribution of the three types of motion analyzed by PCA, an approximately equal role of angles and dihedrals was observed for GROMACS, AMBER(RESP), and AMBER(bcc), while CHARMM prioritized contributions due to angles (compare sizes of blue and green areas of the pie charts in Figures 2(a)-2(d)). This may reflect the extra Urey-Bradley term of the CHARMM force field that was included in the contribution of angles when forming averages (green sector in Figure 2(c)). Relative contributions resulting from bonds (red sectors in the pie charts of Figure 2) turned out to be comparable between CHARMM and AMBER(RESP) but were less significant in AMBER(bcc) and entirely absent in GROMACS. The latter effect is the consequence of putting restraints on all types of bonds (GROMACS default mode [35]), not just on C-H bonds as commonly done by the other force fields. Given the rather pronounced differences between AMBER(RESP) and AMBER(bcc) it appeared interesting to also examine different charge assignments and their corresponding electrostatic potentials (ESPs). Because we are interested in the effect of cholesterol on structure and function of cellular membranes, we applied a corresponding continuum description of such an environment [38,39]. In practice, this means that only the OH-group and a small number of adjacent atoms of the cholesterol ring system (positions 2, 3, and 4 in Figure 1(a)) become exposed to the domain of polar head groups of the cellular membrane and thus will be assigned a medium of higher dielectric constant (i.e., comparable to methanol [38]), while the rest of the cholesterol structure will face the hydrophobic core environment of the membrane interior (modelled by cyclohexane) [15,38]. We computed mean structures from the 250 explicit snapshots for all the 4 force fields tested and picked those snapshots that came closest to the average structure. Membrane specific ESPs were then computed based on the continuum description outlined above using an efficient computer model [40] that had recently been applied in a number of related studies of membrane associated biochemistry [15][16][17]41]. Results are summarized in Figure 3  develop significant areas of nonneutral ESPs except for a single spot of positive ESPs in the vicinity of the OH-group. This is the natural consequence of a very restricted charge assignment by GROMACS involving only the OH-group and the anchor C-atom of cholesterol (position 3 in Figure 1(a)). Such marked differences in membrane internal ESPs may gain significant importance in properly explaining basic modes of receptor activation and signal transduction [43].

Conclusions
In conclusion, the comparison of common force fields described here reveals a largely unifying picture of the structural dynamics of cholesterol and an increasing tendency of force field independence with more complex degrees of freedom such as angle bending and dihedral rotations ( Figure 1). The methodic focus on dynamic aspects highlights the usefulness of nonenergy based techniques like, for example, PCA [30]. Our results clearly demonstrate that such a thermodynamic similarity is far from being obvious when strictly taking into account only partial contributions of kinetic and potential energies (Figure 2). In addition, we point out that particular care must be taken of realistic charge assignments for membrane-embedded compounds ( Figure 3) since the effect on biomolecular interactions may be profound and consequences on biological reasoning may be severe [43].