Systematic implicit solvent coarse-graining of bilayer membranes: lipid and phase transferability of the force field

We study the lipid and phase transferability of our recently developed systematically coarse-grained solvent-free membrane model. The force field was explicitly parameterized to describe a fluid 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC) bilayer at 310 K with correct structure and area per lipid, while gaining at least three orders of magnitude in computational efficiency (see Wang and Deserno 2010 J. Phys. Chem. B 114 11207–20). Here, we show that exchanging CG tails, without any subsequent re-parameterization, creates reliable models of 1,2-dioleoylphosphatidylcholine (DOPC) and 1,2-dipalmitoylphosphatidylcholine (DPPC) lipids in terms of structure and area per lipid. Furthermore, all CG lipids undergo a liquid–gel transition upon cooling, with characteristics like those observed in experiments and all-atom simulations during phase transformation. These studies suggest a promising transferability of our force field parameters to different lipid species and thermodynamic state points, properties that are a prerequisite for even more complex systems, such as mixtures.

changes of lipid topology or temperature-without further re-parameterization-result in lipid membranes with satisfactory physical properties, thus making headway towards the ultimate goal of studying complex membranes using efficient CG force fields. We will argue that this aim indeed seems achievable.

Coarse-grained (CG) model and force field
Each 1,2-dipalmitoylphosphatidylcholine (DPPC), POPC and 1,2-dioleoylphosphatidylcholine (DOPC) lipid molecule is mapped onto a structure consisting of 15, 16 and 17 CG sites, respectively. This requires eight different bead types for characteristic chemical moieties (see figure 1), which we denote as follows: CH for the choline entity of the head group; PH for the phosphate group; GL for the glycerol backbone; E1 and E2 for the ester groups of the sn-1 and sn-2 tails of the lipid, respectively; AS, AD and AE for the hydrocarbon groups (CH 2 CH 2 CH 2 ) , (CH 2 CH 2 ) and (CH 2 CH 3 ), respectively. It is important to distinguish between the unsaturated alkyl groups AD and the saturated AS groups, because lipid chain orientational order and bilayer phase behavior are sensitive to the position and number of alkyl unsaturated bonds [7,40], as we will also see in the remainder of the paper.
The CG model for POPC has been parameterized via systematic structure-based coarsegraining, as explained in detail in [1]. Briefly, the bonded and non-bonded interactions, together with the effective cohesion mimicking the hydrophobic effect, have been derived by matching the structural properties of a POPC bilayer from experiments and AA bilayer simulations. In the present paper, we transfer the parameters of this CG force field to the topologically different lipids DOPC and DPPC without any additional modification.

All-atom (AA) simulations
Our AA simulations were performed with the parallel Molecular Dynamics (MD) program NAMD [41]. All systems were simulated using the fully atomistic CHARMM27 [42] parameters to describe the lipid interactions. The visualizations were performed with VMD [43].
Periodic boundary conditions were applied and a constant temperature was maintained using a Langevin thermostat with a damping coefficient of 0.5 ps −1 (T = 310 K for both POPC and DOPC membranes and T = 323 K for the DPPC membrane). Constant pressure (P zz = 1 atm) in the simulations was obtained using a Langevin-piston barostat with a piston period of 2 ps, a damping time of 2 ps and a fully anisotropic pressure coupling [44,45]. The long-range electrostatic interactions were computed every time step with the Particle Mesh Ewald (PME) algorithm [46], employing a real-space cutoff of 12 Å. The integration time step was 1 fs.
The smallest square membrane system consisted of 72 lipids fully hydrated with approximately 4000 water molecules and equilibrated for approximately 30 ns in the N P zz AT ensemble, with the pressure normal to the membrane P zz = 1 atm and the projected membrane area A = L x × L y corresponding to an area per lipid of 68.3 Å 2 for a POPC membrane, 72.4 Å 2 for a DOPC membrane and 64.3 Å 2 for a DPPC membrane, which are the saturated areas measured in experiment [10,47]. Some of the initial configurations of the small membranes were obtained from previous publications [48]- [50]. A larger membrane system with 288 lipids was then created from this base system by replicating it along both the x-and y-direction.

CG simulations
The CG MD simulations were performed in both the N V T and the N T ensembles, where is the lateral bilayer tension. The CG simulation units are ς = 1 Å (length), ε = k B T 4.28 × 10 −21 J 2.58 kJ mol −1 0.617 kcal mol −1 (energy) at T = 310 K, m = 1 Da (mass) and τ = ς √ m/ε 0.062 ps (bare time corresponding to instantaneous dynamics such as the kinetic energy). The reduced molecular friction in the CG system implies a different time mapping at longer time scales. Using the asymptotic lipid diffusion as the basis of mapping, one finds an approximate scaling of the order of τ 0.6-3.2 ps. We employed a CG integration time step of δt = 0.1τ .
The constant temperature in all the CG simulations was achieved via a Langevin thermostat [51] using a friction constant = 0.2τ −1 . Constant-area simulations were implemented using a quadratic (L x = L y ) cross-sectional area in the bilayer x y-plane. In the N T ensemble the tension was controlled via a modified Andersen barostat [52], allowing for box resizing only in xand y-direction (with a box fiction box = 4 × 10 −5 τ −1 and a box mass Q = 10 −5 -10 −4 m). Periodic boundary conditions were applied in all three directions. All CG simulations were performed using ESPResSo [53].

Lipid transferability
Compared to top-down CG bilayer models [54]- [64], our systematic solvent-free CG bilayer model not only possesses an improved computational efficiency, but also preserves much of the chemical specificity and structural accuracy [1]. As in other bottom-up CG models [28]- [39], this enables its application to a specific rather than a generic problem, though the ability to transfer parameters from one lipid or thermodynamic state point to another then becomes crucial. The reason is not just that otherwise the costly parameterization work would need to be repeated for every new system-even more crucially, when simulating mixtures the system might phase separate into regions of varying composition, thus requiring a force field that gives a fair description of all the coexisting phases simultaneously.  Figures 2 and 3 show a detailed comparison of RDFs among all CG types between AA and CG simulations of the DOPC bilayer at T = 310 K and the DPPC bilayer at T = 323 K, respectively. The higher temperature in the DPPC case was chosen because the membrane would otherwise be in the gel phase (see below), while with the present choice both DOPC and DPPC are in the fluid L α phase, consistent with experiment [10].

Radial distribution functions (RDFs)
The comparison between AA and CG RDFs in figures 2 and 3 shows that all pair distributions (36 for the DOPC bilayer with eight types of CG beads in figure 2 and 28 for the DPPC bilayer with seven types of CG beads in figure 3) obtained from AA simulations are quantitatively reproduced by the CG force field. The differences in RDFs between AA and CG simulations for DOPC and DPPC bilayers are quantitatively comparable with those for a POPC bilayer (see [1]), even in the worst-converging case with the largest difference in RDF (see figures 2 and 3 for CH-E1). This indicates that although the CG force field was derived to match structures of a fluid POPC bilayer at T = 310 K, a good structural transferability is achieved for the fluid DOPC and DPPC lipid bilayers. Figure 4 compares the RDFs between the DOPC bilayer at T = 310 K and the DPPC bilayer at T = 323 K. These pair distributions are based on the same set of interaction potentials, but they do not coincide due to the different macroscopic bilayer characteristics that emerge from the change in lipid (one additional double-bond bead 'AD' in each tail of DOPC) and the slightly different temperature. The biggest effect is probably the change in area per lipid, combined with the increased orientational order in the tail beads for DPPC.

Saturated area per lipid
The saturated area per lipid is a central quantity measuring the lateral organization of bilayers. We have previously shown that our CG model yields a saturated area per lipid of A = (69.8 ± 0.5) Å 2 at T = 310 K for a tensionless POPC bilayer, which is consistent with the experimental   value A = (68.3 ± 1.5) Å 2 at T = 303.15 K [10]. The area per lipid was calculated from the cross-sectional area of the simulation box (the x y-plane of the bilayer) divided by the number of lipids per leaflet (144 lipids in our simulations). Without modification of the CG force field, our model yields a saturated area per lipid of A = (79.7 ± 0.8) Å 2 at T = 310 K for a DOPC bilayer and A = (65.9 ± 0.6) Å 2 at T = 323 K for a DPPC bilayer (see figure 5). Table 1 illustrates how these values compare to experimental results (also taking into account the thermal area expansivity of bilayers) [6,10,47,65,66]. The thermal area expansivity of the DOPC bilayer membrane was obtained from [47]. To the best of our knowledge, experimental values for the thermal area expansivity of POPC are not available. In table 1, as a reasonable estimate [67], the thermal area expansivity of POPC is taken as the average of DMPC [65] and DOPC [47]. We see that the value for DPPC is almost exactly matched (2.5% error), even though both the lipid and the temperature have been changed. In the case of DOPC, the agreement is fair but not as good (8% error).
Given that these results are obtained without any re-parameterization of our force field, we consider them to be a rather encouraging illustration of the fact that transferability can  be obtained. This is by no means trivial: the saturated area per lipid emerges as a balance of tensile forces in the interfacial region between heads and tails of lipids and the pressure of the compacted tails and head groups, until the integral across the bilayer over the difference between tangential and normal stress vanishes. In our case, this implies a balance of bonded and nonbonded interactions, together with the phenomenological cohesion mimicking the hydrophobic effect. But even if we can, in principle, capture density changes in the tail region fairly well, the same is not necessarily true for the head group region: here the CG interactions have an implicit water contribution parameterized in, and since the hydration of the head groups changes with area per lipid, the water contribution to this interaction might change, an effect one cannot expect to be included in our potentials. That the DPPC and DOPC cases are nevertheless remarkably well captured is thus unexpected but gratifying.

P 2 order parameter
The orientational order parameter of each lipid bond, P 2 , is usually defined as a measure of the lipid alignment with the bilayer average normal [7]:  6. Comparison of the orientational order parameter P 2 (see equation (1)) of intra-molecular bonds between AA (red) and CG (blue) simulations for DOPC (a) and DPPC (b) bilayer membranes. The bond numbers listed on the horizontal axis are indexed in the schematic diagram. The DOPC bilayer with one unsaturated aliphatic AD group per tail shows a significantly lower bond order in its tail than the DPPC bilayer with two fully saturated hydrocarbon tails.
where θ is the angle between the unit vector along some particular CG bond and the average bilayer normal. For a completely aligned bond, S bond = 1; for a completely random (isotropic) bond, S bond = 0; and for a bond perfectly perpendicular to the bilayer normal, S bond = − 1 2 . In figure 6, we compare the P 2 order parameter between CG and AA simulations for DOPC (figure 6(a)) and DPPC ( figure 6(b)). The orientational order of lipids for the unsaturated DPPC @ 323 K DOPC @ 310 K POPC @ 310 K hydrocarbon tails is seen to drop compared to fully saturated chains, and this is captured in the CG model. Similar to the result of POPC (see [1]), the comparison of the P 2 order parameter shows good overall agreement for all lipid bonds between CG and AA simulations. Although the CG force field leads to a slight underestimation of the tail order, it reproduces important structural features corresponding to the lipid bond order from AA simulations. A comparison of the P 2 order parameter for all three types of bilayer membranes is shown in figure 7. At the same temperature, there is no pronounced difference in the P 2 order of the unsaturated sn-2 hydrocarbon tails between DOPC and POPC. However, comparing DPPC with POPC, we find that even though DPPC is at a higher temperature (323 K) than POPC (310 K), the saturated sn-1 hydrocarbon tail of POPC is less ordered than the equally saturated sn-1 tail of DPPC. Evidently, the unsaturated sn-2 tail in POPC and its impact on the overall bilayer state also lower the degree of order in the saturated sn-1 tail.

Line tension
The line tension of a bilayer is the excess free energy of an open edge. To avoid exposing the hydrophobic tail region to the aqueous solvent, lipids curve around the edge, and the balance between core exposure and lipid deformation away from the flat bilayer state increases the free energy [57], [68]- [70]. The resulting structure of the membrane edge has been studied both with experiments and AA simulations [70]- [75]. If a bilayer is put under lateral stress, the opening of a pore constitutes a balance between the strain energy in the bulk membrane and the line energy of the open pore, which implies that the line tension also affects the bilayer rupture tension, even though in an unusual area-dependent way [57,61], [68]- [70].
In our CG model, the hydrophobic effect of the aqueous environment is implicitly accounted for by a phenomenological cohesion between non-headgroup beads, which is optimized to reproduce the neat bilayer state, but not an open edge. Hence, both the line tension and, as a consequence, the overall stability of the bilayer could be adversely affected. To test whether the membrane is overstabilized by the cohesion and to evaluate the transferability of the cohesion in other lipids, we computed the line tension for POPC, DOPC and DPPC membranes.
The line tension can, in principle, be determined from the stress-strain relation of a bilayer all the way up to the regime where a pore opens. The balance of the line tension and the area tension gives rise to a characteristic curve that can be fitted by a simple (ground state) theory [61,69]. Here, we use a simpler method. By spanning a membrane across the x-direction of a periodic box that is not large enough to also span the y-direction, two stable linear open edges appear (see [1]). These edges will exert a force equal to twice the line tension along the x-direction. Note that since the width of the membrane can relax by choosing the distance between the two open edges appropriately, no additional force due to surface tension exists. Hence, the line tension γ is given by the simple formula where σ x x is the x x-component of the stress tensor, and L y and L z are the side lengths of the simulation box in the yand z-direction, respectively. The integration time step chosen was five times smaller than the usual 0.1 τ , since the forces are very small and subject to stronger time-step discretization artifacts than structural observables. This has also been observed previously when estimating the bending modulus from a measurement of the small axial forces in cylindrical membrane tethers [76]. After warm-up and an additional equilibration time of 2 × 10 5 τ , averages were obtained by simulating up to 10 7 τ and measuring the stress σ x x every 2τ (after which the stress auto-correlation has essentially decayed to zero). The measured values of γ are: at T = 310 K, γ = (22 ± 5) pN for POPC and γ = (18 ± 4) pN for DOPC bilayer membranes; at T = 323 K, γ = (19 ± 4) pN for DPPC, γ = (17 ± 4) pN for POPC and γ = (14 ± 4) pN for DOPC bilayer membranes. Typical experimental line tensions for a phosphatidylcholine (PC) bilayer membrane have been reported in the range of 6.5-30 pN [71]- [74]. AA simulations found (12 ± 9)-(35 ± 10) pN for a DMPC membrane [75]. While measurements of the line tension show a considerable spread, preventing a more accurate comparison, it is certainly true that our values are well within the range of both experiments and AA simulations of PC membranes, thus showing that we have not overstabilized the bilayer phase.
It appears that, if same-temperature pairs are compared, the line tension for POPC is larger than that for DOPC, while it is smaller than that for DPPC, suggesting that unsaturated tails provide more flexibility for the hydrocarbon chains to rearrange at the edge of the membrane and more efficiently reduce the excess free energy. Moreover, the line tension of POPC appears to drop with an increase of temperature, again in agreement with the expectation that a greater tail flexibility will lower the edge free energy.

Self-assembly
We followed the self-assembly of a system of 288 POPC lipids randomly dispersed in a box of size (99.17 Å) 3 . The box size was chosen since (99.17 Å) 2 /144 is the saturated area per lipid at 310 K. The initial configuration was generated by choosing both the position and orientation of every lipid randomly (of course without mutual overlap). Parallel tempering was implemented with exchanges between eight replicas at the temperatures 279, 294.5, 310, 325.5, The lipids self-assemble to a disordered liquid-crystalline bilayer at T = 310 K. At T = 279 K, an ordered-disordered mixed domain formed instead. Due to the smaller area per lipid, the bilayer now exhibits a pore. 341, 356.5, 372 and 387.5 K. This way, we not only avoid kinetic trapping, but also achieve self-assembled structures at various temperatures (even though not at the correct area per lipid at that temperature).
At 310 K or higher, the POPC lipids self-assembled into a disordered liquid-crystalline (L α ) bilayer, see figure 8(a). At the two colder temperatures 279 and 294.5 K, the lipids selfassembled into bilayers that showed a gel-fluid coexistence or even a pore, as seen for 279 K in figure 8(b). The pore defect is a straightforward consequence of the mismatch between the actual box size and the size that would correspond to the reduced area per lipid. This suggests that even though our CG force field has been derived to represent a POPC bilayer in the fluid phase, it might be transferable to the simulations of a bilayer in a gel phase.

Annealing
To further test the transformation of a fluid into a gel phase, we studied the temperatureinduced phase transformation of DOPC, POPC and DPPC bilayers, i.e. we perform temperature scans (both cooling and heating) over the region in which the actual transition is expected to occur and monitor the evolution of observables characterizing the bilayer phase. We point out that this is different from a study of the actual phase transition [26,27], because the phase transformation process includes kinetic aspects, e.g. the evolution of intermediates and growth of nuclei in a mother phase, etc. In consequence, hysteresis is usually observed during a phase transformation [77], while phase transition only refers to an equilibrium phase diagram, where an accurate phase transition temperature is usually obtained via an equilibrium free-energy computation [21,78].
We performed a series of annealing simulations (cooling and subsequent re-heating) for DOPC, POPC and DPPC membranes in the N ( = 0)T ensemble. For each type of lipid membrane, more than ten different initial configurations were chosen from a pre-equilibrated system at a high temperature (344 K for DOPC, 366 K for POPC and 350 K for DPPC). The system was first cooled from this high temperature to about 250 K, followed by a re-heating from the low to the high temperature. The rates used were in both cases 1.86 K per 20 000τ , i.e. (see figures 10(a) and (b)) and DPPC membranes (see figures 11(a) and (b)). We encountered a more pronounced hysteresis during heating simulations than has been observed in AA simulations [26,27] (see the curves upon heating in figures 9-11). This is, however, consistent with previous findings from explicit solvent CG simulations and was explained as a system size effect [77]. Briefly, a small bilayer size under periodic boundary condition tends to stabilize the gel phase, leading to a stronger hysteresis upon heating, while the effect is weaker on the phase transformation during cooling. Consider that heating a quite rigid gel phase can lead to an increase in box size only if the entire system melts in one go (as indeed observed in our simulation), whereas upon cooling a fluid phase, more than one nucleus might form and thus give rise to domain boundaries that ultimately need to anneal, but these would not substantially affect the overall area or the value of the order parameter.
With our CG model, the phase transformation happens over the range of 265-295 K for a DOPC membrane, 280-310 K for a POPC membrane and 300-330 K for a DPPC membrane. The DPPC interval contains the experimental values of the actual phase transition temperature as obtained in experiments (315 K, see [79]), whereas the experimental transition temperatures for DOPC and POPC are slightly below our bracketed intervals (about 253 K for DOPC, see [80]; and about 270 K for POPC, see [81]). Concerning the effect of an unsaturated bond on the melting temperature, the right trend is achieved with our CG model: from DPPC to POPC to DOPC, the phase transition temperature drops substantially.
The phase transformation of a lipid bilayer from the liquid-crystalline to a gel phase implies many structural changes. The evolution of the bilayer structure during cooling and heating is shown in figures 12-14. Note that the bilayer size we used, 288 lipids, is not large enough to form a ripple phase [25]. Rather, when cooled, our bilayer displays a coexistence of ordered and disordered domains.
For a DPPC bilayer, the mixed domains were shown at T = 308 and 250 K in cooling and 304 and 332 K in heating (see figure 12). We did not find a tilted lipid arrangement in the ordered domain in the range of our simulated temperatures. The ordered domain consists of lipids that are fully stretched with no overlap between lipid tails between two leaflets (untilted non-overlapping sections. At a lower temperature T = 267 K in cooling, a tilt arrangement was observed with little or no interdigitation between two leaflets. These structures of DPPC and POPC at the molecular level have all been previously observed in AA simulations [26,27]. For a DOPC bilayer, we observed that a disordered L α phase transforms to an ordered phase with all double-bond kinks aligned upon cooling (see figure 14(b)), at the temperature 250 K during cooling and at the temperature 270 K during heating).
Fitting the slopes in figure 9(a) from 295 to 340 K, we obtain an area change of (0.52 ± 0.05) Å 2 K −1 for a DOPC bilayer in the L α phase, which at 310 K and A = 73.8 Å 2 translates to a thermal expansivity of 0.0070 K −1 . Doing the same for the POPC membrane and fitting the slope between 310 and 340 K in figure 10(a) leads to an area change of (0.54 ± 0.05) Å 2 K −1 , which at 310 K and thus 69.8 Å 2 translates to 0.0077 K −1 . For the DPPC membrane, fitting the slope between 325 and 350 K in figure 11(a) leads to an area change of (0.62 ± 0.09) Å 2 K −1 , which at 323 K and thus 65.9 Å 2 translates to 0.0094 K −1 . These values are of the right order of magnitude, but about a factor of two larger than expected from experiments [47,65,66]. Note, however, that their relative ordering comes out correctly: with increasing saturation, the expansivity becomes larger.

Conclusions
Aiming to perform mesoscopic quantitative bilayer simulations, we previously presented a solvent-free CG lipid model with efforts to preserve both computational efficiency and local structural and chemical information. The force field of the CG model was systematically parameterized to represent the L α phase of a fully hydrated POPC bilayer at 310 K. Since a reliable performance of any CG model away from its point of parameterization, or the usage of the same interaction potentials for different molecular topologies, is never guaranteed, we investigated the lipid and phase transferability of this solvent-free CG lipid model. Concerning lipid transferability, we simulated bilayers consisting of either DOPC or DPPC lipids, which were generated from the original POPC lipid by merely switching the tails. A comparison of bilayer properties such as RDFs, P 2 order, saturated area per lipid, and line tensions between CG simulations and experiments, together with AA simulations, suggested a very good force field transferability.
Concerning phase transferability, self-assembly of POPC from parallel tempering CG simulations with eight different temperatures showed that the CG force field drives a random lipid dispersion into disordered liquid-crystalline bilayers at 310 K (or warmer), whereas lipids formed bilayers with gel-fluid coexistence at temperatures below about 294.5 K. We further performed a series of annealing simulations between 250 and 360 K for DOPC, POPC and DPPC bilayer membranes. Hysteresis curves of area per lipid and P 2 order of hydrocarbon bonds showed that the phase transformation temperatures lie in the range of 265-295 K for DOPC, 280-310 K for POPC and 300-330 K for DPPC bilayer membranes. Fitting the temperature behavior of the area per lipid during the cooling scan, we found values for the thermal area expansivity that are approximately a factor of two larger than typical experimental values, but in the expected relative ordering (expansivity drops with increasing unsaturation). The temperature-induced change in structure of DOPC, POPC and DPPC bilayers (288 lipids) was recorded during cooling and heating scans. We observed the same lipid structures of mixed order-disorder domains as in AA simulations at a similar bilayer size. Taken together, these results indicate that good force field transferability can be achieved in systematically structurebased CG lipid models.