Parameters for Martini sterols and hopanoids based on a virtual-site description

Sterols play an essential role in modulating bilayer structure and dynamics. Coarse-grained molecular dynamics parameters for cholesterol and related molecules are available for the Martini force field and have been successfully used in multiple lipid bilayer studies. In this work, we focus on the use of virtual sites as a means of increasing the stability of cholesterol and cholesterol-like structures. We improve and extend the Martini parameterization of sterols in four di ff erent ways: 1—the cholesterol parameters were adapted to make use of virtual interaction sites, which markedly improves numerical stability; 2—cholesterol parameters were also modified to address reported shortcomings in reproducing correct lipid phase behavior in mixed membranes; 3—parameters for ergosterol were created and adapted from cholesterols; and 4—parameters for the hopanoid class of bacterial polycyclic molecules were created, namely, for hopane, diploptene, bacteriohopanetetrol, and for their polycyclic base structure. C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http: // dx.doi.org / 10.1063 / 1.4937783]


I. INTRODUCTION
Cholesterol is a vitally important molecule for animals. It serves as a precursor to a wide range of hormones and is an important modulator of lipid membrane fluidity. Cholesterol has also been implicated in the raft hypothesis as one of the enriched constituents of lipid rafts, which have putative roles in protein localization and signaling. 1 Other related molecules play just as relevant roles in other kingdoms, namely, ergosterol in fungi and hopanoids in bacteria (Fig. 1). The latter, found throughout the geological record, have even been dubbed world's "most abundant natural product." 2 In both fungi and bacteria, the roles of these molecules are also tightly connected to the regulation of membrane fluidity. 3 The observation that they might help their hosts' membranes withstand high concentrations of metabolic by-products-such as ethanol-has made fungal sterols and bacterial hopanoids promising molecules for biotechnology applications, where it is hoped they can be used to increase industrial bioproduction yields. 4 Molecular dynamics (MD) simulations can provide unparalleled structural detail and have been extensively used to understand the role of sterols and hopanoids at a molecular level. [5][6][7] The use of coarse-grained (CG) methods, such as the Martini model, 8 allows a significant extension of the time and scale ranges of MD studies 9 compared to fine-grain (FG) approaches. The relevance of sterols has been recognized in the development of the Martini model, which includes parameters for cholesterol 8 and related molecules cholate and cholesteryl oleate. 10 The Martini cholesterol model has been successfully used in several studies and shown to reproduce lipid phase a) Electronic mail: s.j.marrink@rug.nl segregation in mixed-component bilayer systems. 11 These observations show that the Martini force field has enough detail to capture the cholesterol-phospholipid interactions responsible for macroscopic bilayer domain formation.
In spite of its successful applications, the Martini cholesterol is not without limitations. It is modeled as a stiff mesh of bonds and linear constraints ( Fig. 2(a)-Current CG model), which limits the numerical stability of the model. Furthermore, it has been reported that the Martini cholesterol fails to properly preserve fluidity of liquid-ordered (lo) domains. 12,13 In this work, we reparameterize the Martini cholesterol model to remedy these deficiencies. We stabilize it by reformulating the spatially constrained cholesterol topology using virtual interaction sites. We further improve the cholesterol model by adjusting its packing properties to allow a proper lo cholesterol phase behavior.
Besides cholesterol or closely related molecules, the Martini force field lacks models for other sterols or the hopanoid class. Based on the new cholesterol topology, we extend the Martini model representation to include new parameters for ergosterol and several hopanoids. Finally, we validate these new models against data on their behavior in membranes. Namely, we compare the spatial distribution of these molecules and their ordering effect on lipid tails.

Use of linear constraints
The Martini model is a coarse-grain approach that represents a system with fewer d.o.f. than counterparts with atomistic detail. Typically, four non-hydrogen atoms are FIG. 1. Chemical structures of the sterol and hopanoid molecules parameterized in this work. The ring hydrogens were included to clarify chirality but are not present in the united-atom models of these molecules.
center-of-mass-mapped to a single Martini super-atom or bead. There are, however, subtler ways to remove d.o.f. from a system. Linear constraints 14 are one such way. Each of these constraints removes a single d.o.f. from the system by fixing the distance between a particle pair. In FG simulations, bonds involving hydrogen atoms, with high-frequencies and lowamplitudes, are typical targets for constraining in biomolecular simulations.
The Martini model also makes use of linear constraints. Bonds, in Martini, typically span several Ångstrom with amplitudes over 2 Å, but can be much shorter and stiffer. In this context bonds with oscillation amplitudes smaller than 0.5 Å are good candidates for constraining, since the reproduction of such narrowly distributed distances will require bonds with relatively high spring constants, or equivalently steep potentials. Such interparticle distance distributions usually arise when mappings finer than 4-to-1 are used: the mapped centers-of-mass lie closer together and their mapped distance is affected by the rotations of fewer underlying atomistic bonds. This is the case of the Martini cholesterol ring moiety, where several beads are bound very close together at quite rigid interparticle distances.
Linear constraints have their own stability limits. Solvers become unstable when dealing with the tridimensional bond topologies 15 sometimes found in CG applicationsnamely, of particular interest to this work, in cholesterol ( Fig. 2(a)-current CG model). Because of this, the current Martini cholesterol topology is only partly constrained even though its narrow bead movement amplitude certainly warrants a fully constrained network. Unconstrained pairs are then held together by stiff bonds-with force constants above 10 4 kJ/(mol nm 2 ). This combination of a partial constraint network and stiff bonds stabilizes Martini cholesterol only up to simulation time steps of 20-30 fs-and even then large systems can spuriously become unstable, most noticeably if simulated over long timescales. Simulations at a 40 fs step, at which Martini lipids are stable, 16 fail within the nanosecond range (see the supplementary movie 1 and its description in Section S1, 17 for a visualization of a simulation failure of Martini cholesterol in a bilayer, at a 40 fs time step; notice how cholesterol's particles move at a much higher-and unneeded-frequency relative to that of the lipids).

Virtual interaction sites
One of the reasons why a complex mesh is required to keep Martini cholesterol together is that its beads have well-defined tridimensional positions relative to each other. Because bonds and linear constraints operate at a unidimensional level several are needed in order to effect the tridimensionally constrained cholesterol topology.
An alternative to constraints and stiff bonds is to absolutely specify the position of some of the molecule's particles as a function of other particles' positions. This is the so-called virtual interaction site approach, 15 and here we explore its applicability to the Martini sterols and the structurally related hopanoids. During the MD calculation of potentials and forces, a particle that is a virtual interaction site is considered normally: bonded and nonbonded forces acting on the particle are calculated from the derivative of the underlying potentials with respect to its position. However, the virtual particle itself will not be accelerated by those forces. Instead, the forces are propagated to the non-virtual particles that construct the virtual site. After force distribution, non-virtual particles are accelerated and displaced normally. The virtual site position is then recalculated from the new positions of the constructing particles.
Different implementations of virtual sites exist, some even specifically targeted at CG systems. 18 In this text, we will refer to the implementation in the GROMACS 4.6 simulation package, 19 which extends the work laid out by Stillinger and Rahman 20 and Feenstra et al. 15 This implementation is mainly intended to spatially constrain hydrogens and aromatic carbons, 15 as well as displaced charge points, 20 in atomistic systems. Nevertheless, the GROMACS interface is flexible enough that the application of the virtual-site approach to the Martini models is straightforward.
GROMACS provides different virtual site constructs to address different constraining needs. In this work, we make use of two such constructs (Fig. 2(c)). In the first, three nonvirtual particles define two vectors that are linearly combined to construct a coplanar virtual site. From constructing particle positions j, k, and l, the virtual site position i is given by where a and b are the two scalar coefficients of the combined vectors and the only needed parameters. The second construct also combines the vectors from three base particles to calculate i but allows it to be out-of-plane through addition of the cross product of − → j k with − → jl, In this construct, parameter c is the coefficient of the added cross product vector (the in-plane construct is effectively a particular case where c = 0). We refer the reader to Refs. 21 Bead colors indicate polar type, from black, most apolar, to blue, the hydroxyl-containing bead. The labeling of the ergosterol model highlights the differences to the cholesterol model; these also include a shorter and straighter C1-C2 tail bond. In virtual-site models, each virtual site is labeled as "V." Dashed lines connect each virtual site to its constructing frame. Green bonds indicate the constraints over which a hinge dihedral angle potential is applied (see Section S7 for details on the applied potential 17 ). The dashed regions on the united-atom model indicate the two atom groups to which planes were fitted to estimate the flexibility of this cholesterol hinge. In (c), particles i 1 and i 2 are constructed in-plane and out-of-plane from the positions of particles j, k, and l, according to Eqs. (1) and (2) (in this example for both particles a = −0.55 and b = 0.7; for particle i 2 c = 1). and 22 for further details on virtual site construction and implementation, namely, the propagation of forces from the virtual sites to the constructing particles.
In the virtual-site models used in this work, forces acting on virtual particles arise essentially from the intermolecular nonbonded Lennard-Jones interactions dictated by the virtual site's bead type (Coulombic forces will not act on the virtual particles since none of the Martini sterol or hopanoid models carry charges). In one case ( Fig. 2(a)-single-frame CG model, bead C1), a virtual particle is also involved in the calculation of intramolecular bonded forces.

B. Simulation and analysis details
The GROMACS 4.6 package 19 was used for all simulations (both parameterization and validation runs), except for those requiring restricted angle bending potentials, 23 for which GROMACS 5.1 was used.

General CG simulation conditions
CG data were generated with the Martini model and simulated at either 20 or 40 fs integration steps, depending on whether the less stable current cholesterol model was used or the newly parameterized models. Standard Martini Lennard-Jones and Coulomb potentials, shifted to 0 at 1.2 nm, were used 8 and neighbor list updates done every 10 steps with a 1.4 nm radius. For target CG runs or during parameterization of bonded parameters, temperature was coupled to 300 K using the Berendsen thermostat 24 with a 1 ps coupling time. Pressure was isotropically coupled to 1 bar using the Berendsen barostat 24 with a 3 ps coupling time. Bilayers were simulated with semi-isotropic pressure control, coupled independently in the x y and z directions. The standard Martini water model was used. 8 All bilayers were generated with the insane method, 25 with lateral randomization of lipid/sterol positions.
Hexadecane-water partition free-energies were calculated from the individual solvation free-energies into each solvent. In CG simulations to estimate solvation freeenergies, a time step of 20 fs was used; temperature and pressure were coupled with the v-rescale thermostat 26 and the Parrinello-Rahman barostat. 27 Single cholesterol/hopanoid molecules were decoupled from solvent boxes of either 5155 Martini waters or 320 hexadecane molecules. Decoupling was performed in 12 steps of 4 ns each, by scaling down solute-solvent Lennard-Jones interactions (there are no Coulombic interactions in these Martini systems). For bacteriohopanetetrol restricted bending potentials were switched off as the GROMACS 5.1 implementation does not yet allow their use together with free-energy calculations (even without these potentials the simulated systems were stable enough to withstand the short decoupling runs).
In the Martini model, first bonded neighbor beads are excluded from nonbonded interactions. In addition, all ring beads in the sterols/hopanoids are excluded from intramolecular nonbonded interactions with one another, as is already the case with the current cholesterol model. This means, in particular, that virtual sites will not interact in any way with the respective constructing particles, besides the aforementioned distribution of forces. The sterol/hopanoid ring beads can only establish intramolecular nonbonded interactions with tail beads. Further specific nonbonded exclusions will be detailed in the respective parameterization sections.
Martini ring moieties are represented in the Martini model by S-type beads, designed to promote correct packing Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 129.125.14.248 On: Tue, 01 Mar of these structures. These beads interact with one another with a shorter equilibrium Lennard-Jones distance (0.43 nm instead of 0.47 nm) and a potential well-depth scaled by 0.75. Interactions with regular Martini beads use unmodified potentials. 8

General FG simulation conditions
FG reference systems were simulated with the unitedatom GROMOS 54a7 force field 28 together with the single point charge water model. Systems were run at a 2 fs time step. Water bonds and angles were constrained using the SETTLE algorithm 29 and the sterols' O-H bonds using the LINCS algorithm. 14 Lennard-Jones interactions were evaluated with a twin-range scheme where interactions up to 0.9 nm were evaluated every step and between 0.9 and 1.4 nm every 10 steps (the same frequency as the neighbor list update). Coulombic interactions were evaluated up to 1.4 nm with the reactionfield correction; the reaction-field dielectric constant was 65.0 for systems in water and 2.14 for systems in hexadecane.
The v-rescale thermostat 26 was used to couple temperature to 300 K, with a 1 ps coupling time. Pressure was coupled to 1 bar using the Berendsen barostat, 24 with a coupling time of 3 ps. For free-energy estimates, the Parrinello-Rahman barostat 27 was used instead.

Target data
Target CG data for cholesterol virtual-site parameterization were obtained from a 3.2 µs Martini simulation of a single cholesterol molecule in a dipalmitoylphosphatidylcholine (DPPC) bilayer (127 lipids), coupled to a 323 K bath (above the gel-to-fluid transition temperature of DPPC).
FG target data were used when adjusting cholesterol parameters to create the ergosterol topology and when tuning the cholesterol's virtual-site frame hinge stiffness. Single molecules of either sterol were simulated in solvent boxes of water or hexadecane for at least 130 ns. Sterol and hexadecane parameters are the ones available at the Automated Topology Builder repository 30 (care was taken to choose manually submitted and validated parameters not automatically generated ones; the hexadecane parameters were used without any partial charges on the united carbons).
Just as for the FG cholesterol target data, the four hopanoids (the polycyclic moiety, hopane, diploptene, and bacteriohopanetetrol) were simulated in water or hexadecane for at least 60 ns. The used GROMOS 54a7 bacteriohopanetetrol and diploptene parameters are the ones described in Ref. 6. Minor changes were made to the diploptene parameters, within the GROMOS standard bonded potentials and atom types to obtain parameters for hopane and the truncated (chain-less) polycycle moiety.
As for CG, FG hexadecane-water partition free-energies were calculated from the individual solvation free-energies. The free-energies of cholesterol and hopanoid solvation were estimated at the FG level by gradually decoupling single molecules from solvent boxes of either 2495 single-pointcharge (SPC) waters 31 or 320 hexadecane molecules. No O-H bonds were constrained in these systems; in the case where they were present (cholesterol and bacteriohopanetetrol), the simulation time step was reduced to 1 fs. Each decoupling step was simulated for at least 4 ns. Only in the cholesterol-water, bacteriohopanetetrol-water, and diploptene-water systems are there solute-solvent charge interactions. In all other systems, either the solvent is chargeless (hexadecane) or the solutes are. In the latter cases a shorter (12 steps) decoupling process scaling only the Lennard-Jones interactions was used, versus a 24-step decoupling for the cholesterol/bacteriohopanetetrol/diploptene-water systems, where Coulombic and Lennard-Jones interactions were switched off sequentially.

Parameterization iterations
The tuning of CG bonded parameters (the stiffness of the hinges, and the bonds, angles, and dihedrals of hopanoid and ergosterol chains) was done by trial and error of different potential combinations. This repeated testing was performed in water with runs of at least 40 ns. Parameter tuning in other solvents was deemed unnecessary as only in one well-defined case were there any observed structural differences between runs in water and in hexadecane (namely, a high degree of intramolecular hydrogen bonding in the polyol tail of bacteriohopanetetrol was observed in hexadecane, where it caused a significant straightening of the chain. This effect was still predominant in water, but a small contribution of a bent tail configuration could also be observed).
For the tuning of bead type assignments, solvation freeenergies were calculated in water and hexadecane with the CG decoupling scheme described in Sec. II B 1.

Validation
The depth distribution and flip-flop rate of the virtualsite cholesterol model was compared to that of the current model in cholesterol-doped bilayers at a 3:1 lipid:sterol ratio. Palmitoyl-oleyl-phosphatidylcholine (POPC) was chosen as the host lipid due to an expected higher flip-flop rate (and hence better statistics) than in a fully saturated bilayer. Systems were simulated with a total of 720 lipids and for at least 1 µs.
Phase behaviour was studied on ternary DPPC:DLiPC: cholesterol mixed membranes (DLiPC: dilinoleoylphosphatidylcholine) with 42:28:30 molar ratio and a total of 2000 lipids. 30 fs time steps were used, consistently with previous Martini simulations of phase-separating membranes. 11 Spatial distributions of cholesterol neighbors in the bilayer plane were obtained from either DPPC or dioleoylphosphatidylcholine (DOPC) bilayers with 20% cholesterol, with a total 450 lipids. Simulations were run at 338 K; these conditions mimicked the ones used by Martinez-Seara et al. 32  representation of distearoylphosphatidylcholine). Only cholesterols within 1.0 nm of lipid phosphate groups were counted either as distribution references or as neighborsthis selection excludes horizontally tilted cholesterols in the bilayer core from analysis.

Analysis
Data analysis, such as the generation of bond, angle, and dihedral-angle distributions, was carried with tools provided in the GROMACS 4.6 suite. More complex procedures, such as moment-of-inertia (m.o.i.) fitting or flip-flop counting, were implemented in Python with extensive use of the MDAnalysis 34 and NumPy packages. All analyses discarded 1 ns of initial trajectory time as equilibration. Visualization and image generation were carried out using the VMD 1.9.2 package 35 and the Tachyon ray-tracer. 36 A constrained optimization by linear approximation algorithm, as implemented in the SciPy package, 37 was used for the optimization of the distribution of masses. The process minimized the sum of the squares of the differences between test and target m.o.i. tensor components.
Root-mean-square-displacements (RMSD) were calculated for each structure pair in an all-to-all cross-comparison between reference simulations and runs with the new models. The structures in a reference trajectory were also compared allto-all to themselves (to prevent correlation bias only structures at least 200 ps from each other were compared). RMSD values were calculated from the ring structures only (sterol and hopanoid chains were excluded), overlapped by rotationally and translationally fitting the frame beads.
Partition free-energies and associated error estimates were obtained using the multistate Bennett acceptance ratio (MBAR) method. 38 Flip-flop rates were calculated using a double threshold approach, whereby a flip was only counted when the polar bead of a cholesterol successively crossed imaginary horizontal planes lying at −0.7 nm and +0.7 nm from the bilayer center along the z-axis. This eliminates much of the noise that arises from a single threshold approach where a particle that lingers at the threshold position, flickering above and below it, can cause an artificial increase in measured flip-flop rates.
Order parameters (S) of CG bonds were calculated from their angle with the bilayer normal (θ) averaged over lipid molecules and trajectory frames, 39 A 0.95 confidence interval for flip-flop rates and bond order parameters was estimated from a bootstrap approach with 2000 resamplings, as implemented in the SciPyassociated SciKits package.
Distributions of structural parameters were always plotted as probability density functions (PDFs).

III. RESULTS AND DISCUSSION
A. Cholesterol reparameterization

Strategy
As stated, our goal in reparameterizing cholesterol is twofold: to stabilize its topology through the use of virtual sites and to improve its lo phase behavior. The reparameterization of cholesterol's 7-bead sterol moiety (which roughly corresponds to pregn-5-ene-3-ol) will also be used to directly update the existing parameters for cholate and cholesteryl oleate. The alkyl chain of cholesterol-its eighth bead-will be left unmodified, as it is connected by a standard low-frequency Martini bond.
Our objective is also to maintain the characteristics of the current Martini cholesterol model as much as possible, with the exception of those that influence lo phase fluidity. This is to ensure that the new model retains the correct behavior that the current model already displays in many applications. 11,40,41 No effort will be done to better approximate the current model to atomistic structural behavior (with one exception, see Sec. III A 3). Table I shows that the self-compared RMSD value of the current model's sterol moiety is very low (this is the average RMSD of an all-to-all comparison of all configurations in a trajectory with themselves) reflecting its rigid nature and foreboding a successful spatial constrainment approach. In the 7-particle layout of the Martini cholesterol ring topology, five beads are essentially coplanar while the Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 129.125.14.248 On: Tue, 01 Mar remaining two protrude to one side (beads R1 and R5, which correspond to two off-ring methyl groups; see Fig. 2(a)). In a first approach at converting this layout to virtualsite constructs, beads R3, R4, and ROH were used as the constructing particles ( Fig. 2(a)-single-frame CG model).

Adaptation to virtual sites
The three beads define a rigid triangular frame and were kept in place relative to each other using linear constraints. The positions of the remaining 4 beads were built from this frame using either the coplanar (beads R2 and C1) or the out-ofplane (beads R1 and R5) virtual-site constructs. Parameters for virtual site construction were obtained from a trajectory of the original Martini cholesterol in a fluid DPPC bilayer. The positions of cholesterol's ring particles were averaged after rotationally and translationally fitting the entire molecule to the R3-R4-ROH frame, and their placements relative to bead R3 decomposed as a function of vectors R3-R4 and R3-ROH, according to Eqs. (1) and (2). As with the current model, nonbonded intramolecular interactions between the ring particles were switched off. Section S5 details all virtual site construction parameters. 17 The resulting virtual-site cholesterol was indeed stable with a 40 fs integration step, and representing the ring moiety with a rigid set of beads proved to be a good approximation (see in Table I the quite low RMSD values between the configurations sampled by this model and those sampled with the current model). However, spurious instabilities still occurred at the multimicrosecond scale. We attributed these to cholesterol's rings now behaving like an absolutely rigid body. When torque is generated at one end of the moiety, the ensuing rigid-body rotation can cause a very large displacement on the opposite end, possibly resulting in unphysical overlap with other system components.

A shock-absorbing topology
To prevent the problems associated with having a large rigid body the virtual-site approach was decomposed in two. Instead of a single three-bead constructing frame, two were used, with beads R2-ROH-R3 and R3-R2-C1, respectively ( Fig. 2(a)-dual-frame CG model). Bead R1 was constructed from the first frame (out-of-plane) and beads R4 and R5 from the second frame (in-plane and out-of-plane, respectively). The two frames share an edge connecting beads R2 and R3 and were allowed freedom to hinge on it, but with a restraining potential that keeps them mostly coplanar. This hinge effectively works as a shock absorber that prevents rigid-body displacements from spanning the length of the molecule. Virtual site parameters were obtained as for the single-frame case, but the trajectory was separately fit to each frame for the determination of the respective parameters.
To parameterize the stiffness of the frame hinge, we opted to reproduce the atomistic flexibility of the ring moiety. Flexibility was inferred from combined FG simulations of cholesterol in water and hexadecane; a plane was fit to each set of atoms that map to each frame's bead group (protruding methyl groups were left out; see Fig. 2(a)-united-atom FG model). Dihedral potentials of different stiffness were tested in the CG model until dihedral angle distributions matching the FG ones were obtained (see Section S7 for the comparison of the dihedral angle distributions 17 ). The combination of a flexible hinge with the two-frame virtual-site approach eliminated the spurious instabilities observed with the singleframe model, over the same timescale.

Mass distribution
Virtual sites are not accelerated in a simulation and therefore any mass assigned to them becomes an irrelevant parameter. Since forces acting on virtual sites are propagated to the constructing beads, it is only the masses of those particles that need be considered. To keep the total effective mass of the ring moiety unchanged, we redistributed masses over the virtual-site-constructing beads. Mass distribution was optimized to best match the m.o.i. of the current model-in line with the least-change approach taken in this reparameterization of cholesterol. The time-averaged cholesterol structure in a DPPC bilayer, with the masses of the current model, was taken as the reference for the target m.o.i. The final masses for the virtual-site constructing beads are listed in Section S5. 17

Liquid-ordered phase fluidity
Use of the Martini cholesterol with phospholipids of different saturation level (DPPC and DLiPC) has been shown to reproduce expected phase separation. 11 However, attempts at reproducing this behavior have found that the DPPC-and cholesterol-enriched lo phase that forms is overly ordered: 12,13 instead of having a fluid nature, lipid tails and cholesterol pack into a hexagonal lattice characteristic of the so-called solid-order phase. In fact, instead of breaking the tail lattice, 42 the Martini cholesterol seems to stabilize it. 12,13 In these phases, favorable cholesterol-cholesterol stacking interactions also cause the Martini model to form unrealistically long and rigid cholesterol aggregates. A CG cholesterol model has been proposed where interactions with POPC are modulated by adjusting the offplane position of a bead representing the protruding ring methyls. 43 Here, we adapt that approach to address the overcondensation problem by adjusting the off-plane distance of beads R1 and R5. The virtual site construction makes this quite straightforward, as all that is needed is to scale the coefficient c in Eq. (2). Extra protrusion was tested in Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 129.125.14.248 On: Tue, 01 Mar increments of 10%, and at +30% the desired behavior was obtained. The adjustment was very subtle, considering that the off-plane protrusion of beads R1 and R5 in the new cholesterol model is a mere 1.10 Å and 1.41 Å, respectively, from 0.84 Å and 1.08 Å in the current model.
The dual-frame virtual-site cholesterol model with 30% extra off-plane bead protrusion is the final optimized cholesterol model on which validation will focus. The same model has also been adapted for the Dry Martini implicit solvent model 44 but with a 60% increase in protrusion. This is in line with the increased lipid-lipid interaction distances that the dry model requires.

Validation-Structural and dynamic behavior
Upon conversion to virtual sites, the visited conformational space of the new model is still very close to that sampled by the current model, as the cross RMSD comparisons in Table I show. This agreement, attained despite having less than half the number of d.o.f. of the current model, stresses the suitability of the virtual-site approximation.
Cholesterol flip-flop rates in 3:1 POPC:sterol bilayers were compared between the current and the virtual-site models, as a measure of how well the new model retains dynamic properties in relevant systems. The obtained rates of 4.7 ± 0.3 and 4.2 ± 0.2 flip-flops/µs/molecule, for the current and new models, respectively ( Fig. 3(a)), are in good agreement with each other, and with reported rates calculated from transmembranar free-energy profiles of cholesterol. 45 The distribution of membrane localization is unaffected by the conversion to virtual sites ( Fig. 3(b)).

Validation-Phase fluidity
Reproduction of lipid phase behavior was tested with ternary DPPC:DLiPC:cholesterol systems, either starting from a configuration with a randomized lateral lipid distribution or from a phase-separated system. A configuration with an overly ordered lo domain was chosen as the demixed starting point. Fig. 4 (Multimedia view) illustrates how the extra protrusion fluidizes the DPPC-enriched phase and decreases hexagonal packing and long-range order. The behavior of Martini cholesterol is now much more in agreement with that observed at full atomistic detail. 42 Indeed, it has already been successfully used in a large-scale study of mammalian plasma membrane behavior. 33

Validation-Neighbor distribution
The in-plane spatial distribution of neighbors around cholesterol in a bilayer was compared to FG data reported by Martinez-Seara et al. 32 Fig. 5 presents how DPPC acyl chains and cholesterol rings distribute around reference cholesterol molecules. A hallmark of the distributions obtained for atomistic models is a radially multilobed distribution of cholesterol's closest neighbors. 32 It can be seen from distribution symmetry. Not only does the increased asymmetry bring CG results closer to FG ones but also all three effects are likely to be mechanisms underlying the fluidization of over-ordered phases.

Validation-Tilt angle and lipid order
Finally, tilt-angle distributions of the current cholesterol model and the virtual-site model with extra protrusion were compared (Fig. 6(c)). The new model adopts a slightly more tilted position in POPC bilayers. This parameter has been shown to correlate with the ordering ability of cholesterol. 46 Indeed, in agreement with such observations, the current, straighter, cholesterol model orders POPC more than the new, more tilted, model (Fig. 6). This difference in cholesterol tilt and consequent weaker ordering effect is another factor that explains how the new model promotes the correct behavior of lo phases.

Strategy
Ergosterol is chemically very similar to cholesterol. Our aim when parameterizing its Martini model was to capture that Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 129.125.14.248 On: Tue, 01 Mar 2016 13:51:52 similarity to cholesterol, and, if possible, retain key differences that might yield different behavior in specific applications. To this end, we kept the sterol moiety virtual-site parameters from cholesterol and adjusted only (a) the bead type assignment of R2 and C1, to capture their double-bond nature, and (b) the bonded parameters of the alkyl bond, expected to change due to the introduction of a double bond and an extra methyl group in ergosterol (Fig. 2(b)).

Topology
Bead type reassignment was done according to the polarity scale of Martini, which prescribes a C4 type bead for conjugated double bonds (bead R2), and a C3 type for unconjugated ones (bead C1; see Fig. 2(b)). The difference in polarity between bead types C3 and C4 involves, among other changes, a weaker interaction of bead type C4 with bead type C1 (used for saturated aliphatic carbons). S-type variants of the C4 and C3 bead types were used, as in the current cholesterol model, for improved self interaction distances. Bond parameter adjustment was based on mapped atomistic data taken in polar and apolar media (SPC water and hexadecane) of cholesterol and ergosterol. The cholesterolmapping strategy used by Wassenaar et al. 47 was employed to obtain CG representations from atomistic data. A map for ergosterol, with the same number of target beads, was adapted with minor alterations from the one used for cholesterol.
The main alkyl bond parameter differing between ergosterol and cholesterol was its length, shorter by 0.3 Å in the CG-mapped ergosterol (see Section S6 for the bonded distributions and a more detailed description of the used mapping 17 ).
The tail angle with the rings along beads R5-C1-C2 also differed somewhat between cholesterol and ergosterol, being closer to collinear. This angle has no biasing potential in the Martini cholesterol model; a weak potential was applied in the ergosterol model to promote this difference in tail orientation (see Fig. S5 17 ).

Validation
Experimentally, cholesterol and ergosterol behave the same in many ways. Both sterols are known to increase lipid tail order to similar magnitudes. 48 Atomistic studies have revealed subtle differences whereby ergosterol promotes slightly higher lipid ordering. 5,[49][50][51] This difference has been attributed to a straighter tilt angle of the ergosterol rings and tail when in a bilayer, 49 owing to a higher stiffness of the alkyl tail.
We compared the behavior of the CG models of ergosterol and cholesterol regarding the ordering of POPC acyl chains at 30% sterol concentrations. Results show that both Martini sterols increase POPC tail order (Fig. 6). Ergosterol does order POPC tails more than the new cholesterol model, in a small but significant effect. The order of the unsaturated tail of POPC is more affected by ergosterol, likely as a consequence of the stronger interaction between the unsaturated beads of ergosterol and POPC. Tilt-angle distributions of the sterol tail (C1-C2 vector) show a difference compatible with that reported in FG studies, 49 which suggests that this parameter has the most influence on lipid ordering (Fig. 6). Concomitantly with the increased POPC tail order, and in agreement with the trend in FG results, 49 the ergosterolcontaining system also has a slightly smaller area compared to the cholesterol-containing system (a small but significant difference of 0.4%). The new Martini ergosterol model can therefore successfully capture-at least in a semi-quantitative way-differences relative to the effects of cholesterol (see Section S4 for a collection of properties of POPC bilayers at different cholesterol/ergosterol concentrations 17 ).

Strategy
Parameterization of the hopanoid class was carried out independently of existing CG models. The same virtual site, hinge, and mass-distribution approach was taken as in Secs. III A and III B, but the target positions were directly obtained from atomistic data mapped to eight CG centers ( Fig. 7 and Section S8 17 detail the mapping used for the hopanoids). Mass distribution took into account the full atomistic mass of the moiety and, to improve stability, optimization was constrained to let no bead weigh less than 72 amu. Since the hopanoid class is characterized by its pentacyclic moiety, we chose to first parameterize it independently of substituents. Learning from the cholesterol experience, we put emphasis on retaining protruding methyl features. The ring parameters were then reused as a base on which to substitute propanyl, propenyl, and tetrahydroxyoctyl groups when generating the parameters for hopane, diploptene, and bacteriohopanetetrol, respectively. As with cholesterol, the virtual-site approach proved to be an appropriate structural representation of the constrained configuration of the hopanoid ring moiety (see Table I).
Bonded parameters for the ring substituents were obtained by trial and error until satisfactory overlap with atomistic distributions was obtained (see Section S9 for all the relevant distribution comparisons 17 ). Bead type assignment followed the general Martini rules, with some exceptions detailed in Sec. III C 3. SC2 type beads (for ring structures with no double bonds) were used for the entire polycyclic moiety (see Fig. 7 for the detailed Martini bead assignment of each hopanoid). Polar-apolar partition free-energies were obtained by calculating individual solvation energies of each parameterized hopanoid into either water or hexadecane. The quality of the bead type assignment was validated by comparing the partition free-energies to reference atomistic data simulated under the same conditions.

Chain attachment
In the hopanoid topologies, the chain is bound to the polycyclic moiety via the R8-C1 bond. However, the mapping causes bead C1 to lie quite close to the ring beads and to be within repulsion range of bead R7 (for Martini, this occurs below approximately 0.5 nm; see Fig. 7). To prevent clashes, intramolecular nonbonded interactions between beads R7 and C1 were specifically excluded. Bead C1 is further held in place through an angle potential (over beads R7-R8-C1) and a dihedral potential (restraining planes R6-R7-R8 and R7-R8-C1).
In all hopanoids, the short and narrowly distributed connection between beads C1 and R8 would require a quite stiff bonded potential (in excess of 3 × 10 4 kJ/(mol nm 2 )) to be exactly reproduced. Typically, such a short and stiff connection would be replaced by a linear constraint, but we found that appending one to the hinged constraint frames resulted in unstable runs. The same happened if a stiff potential was used instead of a constraint. Finally, constructing this bead through the closest virtual-site frame also made the system unstable, presumably by exceeding the stable rigid-body size. The LINCS constraint solver 14 used by GROMACS 19 seems unable to solve this constraint/bond topology (bacteriohopanetetrol seemed particularly susceptible, perhaps because bead C1 is further bound to bead C2). This might indicate that the hinged frame approach is limited to two virtual-site frames, as any added constraints impair stability. A larger polycyclic molecule requiring more frames may need other methods to couple them together.
In our approach, the compromise to achieve a stable topology was to use lower bond and angle force constants, which still reproduce the bonded distribution well within typical Martini tolerances (see Section S9 17 ). A similarly broader distribution was also settled for in the case of the rather stiff R7-R8-C1 angle.

Bacteriohopanetetrol
The bacteriohopanetetrol chain was constructed with harmonic bond and angle potentials over beads C1, C2, and C3. Placement of the bead C1 relative to the rings followed the same considerations as for hopane and diploptene. Accurate reproduction of the chain configurations also required a weak torsion potential to be placed over the R7-R8-C1-C2 beads. This torsion has the potential of introducing numerical instabilities if angle R8-C1-C2 ever gets close to collinearity (a frequent problem of the softer potentials used when coarsegraining 23 ). The solution was to overlay a recently developed restricted bending potential 23 over angle R8-C1-C2 that ensures that sampled magnitudes remain in a numerically safe range.
Finally, the bacteriohopanetetrol polyol moieties were observed in atomistic reference simulations to frequently lock together through intramolecular hydrogen bonding ( Fig. 7(a)). This locking, almost permanent in hexadecane and very frequent in water, causes an obvious straightening of the chain and halves the intermolecular hydrogen-bonding potential of the involved groups. To properly capture this behavior focus was put on reproducing the straighter chain mode, and less polar Martini bead types were assigned to the polyol beads (P2, instead of the regularly used P4 for diols).

Validation-Partitioning free-energies
A first attempt at the parameterization of the hopanoid polycyclic moiety showed that its hexadecane-water partition free-energy largely exceeded the value obtained by atomistic simulations (a difference of over 20 kJ/mol; see Table II). Such a large difference diverges from the Martini principle of reproducing polar-apolar partition free-energies. A way around it would be to assign beads of a type more polar than SC2 to the polycyclic moiety. We felt this would be an artificial approach that would not guarantee the proper interaction behavior with Martini components other than water and hexadecane.
Our understanding of this problem is that, while Martini beads properly reproduce partition free-energies on their own, when tightly bound together-as is the case of the hopanoid ring structure-building-block partition additivity Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 129.125.14.248 On: Tue, 01 Mar breaks down. This is likely due to closely bonded beads requiring a significantly smaller solvent cavity than the sum of individual ones. The use of S-type Martini beads compounds this problem: this class of beads was designed with a scaled-down self-interaction distance and energy to allow for efficient packing of ring structures or otherwise finely mapped fragments. 8 S-beads, however, will still interact with regular beads with unscaled potentials (an SC1 bead will interact with a P4 one, for instance, just as a C1 would). From the water or hexadecane perspective (both represented by regular beads) the hopanoid ring structure is a dense cluster of unscaled interactions, with a consequent overall exaggerated solvent interaction energy. Similar discrepancies have been recently observed for Martini photosystem II cofactors, which also make extensive use of S-type particles. 52 Smaller ring compounds might suffer less from this problem but it has been recognized that a polarity adjustment is also required, for instance, for the phenyl group. 53 This mismatch prompted us to also check the partitioning behavior of cholesterol, and indeed a similar discrepancy of about 20 kJ/mol was found, independently of the use of virtual sites (Table II). Interestingly, the cholesterol model reproduces quite accurately the free-energy of insertion into a membrane from water, as estimated by Bennett et al. 45 However, even in that study, the same 20 kJ/mol energy difference to FG behavior is found when, instead of being at its preferred depth, cholesterol sits in the bilayer center (Figs. 1B and 1C in Ref. 45)-this being the membrane environment which the hexadecane solvent most closely mimics.
We opted to proceed with these hopanoid parameters in spite of the hexadecane-water partition mismatch, and to address the discrepancy in partition behavior in future work. The importance of this partition comparison, though a hallmark of the Martini model, might also be overestimated when it comes to this kind of membrane-inserted molecules. A proof of this is that the current model of cholesterol we now find to be too hydrophobic does perform very realistically in the many applications already mentioned. Note also that the relative hydrophobicities between the hopanoids, and in comparison to cholesterol, match the FG data very well (Table II).

Validation-Bilayer organization
Hopanoid organization in a POPC membrane was investigated and compared with published FG simulation results for diploptene and bacteriohopanetetrol. 6 There is a very good qualitative agreement of the preferred depth distributions of bacteriohopanetetrol and diploptene (Fig. 8), reached in time scales under 100 ns. Bacteriohopanetetrol prefers a vertical positioning in the bilayer, with the polyol tail sitting among the lipid headgroups, whereas diploptene prefers a more horizontal tilt and a positioning close to the bilayer core. Hopane displays essentially the same organization as diploptene. Two starting configurations were tested for hopane and diploptene, one with chains pointing towards the lipid headgroups and one where the chains point towards the bilayer core. The same preferred positioning was reached in both cases for both hopanoids.
Finally, the lipid ordering induced by bacteriohopanetetrol was compared to that of cholesterol (Figs. 6(a) and 6(b)). Bacteriohopanetetrol is able to order POPC tails to an extent similar to, but lower than that of cholesterol. This difference, most visible in the order of the saturated tail of POPC, is in agreement with the trend in published FG data 6 (see Fig. S3 for a collection of properties of POPC bilayers with different sterol/hopanoid concentrations 17 ).

IV. CONCLUSION
In this work, we present a successful extension of the use of virtual sites to coarse-grain parameterization. Even though the virtual-site approach was used with the Martini model, any particle-based model, no matter how fine or coarse, can employ it as long as a minimum number of spatially constrained particles are present.
The application of this method to the Martini cholesterol yielded a significant stability improvement. The collection of Martini sterols was extended by the addition of ergosterol, from the adaptation of a small number of bonded and nonbonded cholesterol parameters. Finally, the virtual-site approach could also be easily extended to the hopanoids, an entire new class of Martini molecules with promising applications.