Lipid Head Group Parameterization for GROMOS 54A8: A Consistent Approach with Protein Force Field Description

Membranes are a crucial component of both bacterial and mammalian cells, being involved in signaling, transport, and compartmentalization. This versatility requires a variety of lipid species to tailor the membrane’s behavior as needed, increasing the complexity of the system. Molecular dynamics simulations have been successfully applied to study model membranes and their interactions with proteins, elucidating some crucial mechanisms at the atomistic detail and thus complementing experimental techniques. An accurate description of the functional interplay of the diverse membrane components crucially depends on the selected parameters that define the adopted force field. A coherent parameterization for lipids and proteins is therefore needed. In this work, we propose and validate new lipid head group parameters for the GROMOS 54A8 force field, making use of recently published parametrizations for key chemical moieties present in lipids. We make use additionally of a new canonical set of partial charges for lipids, chosen to be consistent with the parameterization of soluble molecules such as proteins. We test the derived parameters on five phosphocholine model bilayers, composed of lipid patches four times larger than the ones used in previous studies, and run 500 ns long simulations of each system. Reproduction of experimental data like area per lipid and deuterium order parameters is good and comparable with previous parameterizations, as well as the description of liquid crystal to gel-phase transition. On the other hand, the orientational behavior of the head groups is more realistic for this new parameter set, and this can be crucial in the description of interactions with other polar molecules. For that reason, we tested the interaction of the antimicrobial peptide lactoferricin with two model membranes showing that the new parameters lead to a weaker peptide–membrane binding and give a more realistic outcome in comparing binding to antimicrobial versus mammal membranes.


INTRODUCTION
Cellular membranes are key promoters and regulators of many biological processes due to their crucial role in segregating the external world from the organism. Small molecule transport, drug permeation, intracellular signaling, and antibody response are all regulated by the cell membrane or by membrane-related components. 1−8 To fully comprehend and ultimately influence the bespoke processes, it is paramount to understand membranes and their constituting lipids in atomistic detail. However, due to the complexity of those systems, researchers have resorted to the use of simplified model membranes, which can be synthesized and characterized in vitro. This enables the individual contributions of the components involved to be disentangled. Indeed, for the cellular membrane to be able to perform different functions, its composition is necessarily complex. Lipids are one of the main components and can be present in up to hundreds of different species. 9 In addition, many transmembrane proteins tessellate the cell surface, promoting signaling pathways and influencing the membrane's structural and mechanical properties. 10,11 Phospholipid bilayers and micelles have been investigated, in particular, as these lipids represent the main components of the eukaryotic and the inner bacterial membranes. Both have been modeled selecting specific phospholipids to emulate the appropriate surface charge or to reproduce the human cell membrane fluidity by introducing, for example, cholesterol. 12,13 As these simplified membranes retain the core characteristics of their different biological templates, 14 they can be used to test the membrane interaction with proteins, peptides, antimicrobial molecules, or drugs.
Experiments can provide global properties of membranes and, despite the great accuracy of techniques like NMR and Xray scattering in measuring the average position of atoms in rigid structures, they face challenges when characterizing the biologically relevant fluid phase, as opposed to the gel one that emerges at lower temperatures. 15−18 Alongside experimental characterization, molecular dynamic (MD) simulations have played a central role in the investigation of the behavior of lipids, due to the atomistic spatial resolution they provide. Therefore, MD simulations complement our understanding of membranes' behavior and are also important for the study of lipid systems in combination with proteins, providing detailed insights into the mechanisms of their interactions. In the past, MD simulations have been successfully employed to reproduce typical phenomena in membranes, such as lipids' flip-flop, 19,20 vesicle formation, 21,22 aggregation into bilayers, 23−25 and stress-induced 26−29 and peptide-induced pore formations. 30−32 Moreover, the implementation of more realistic models of bacterial membranes, by including a more diverse set of components into the simulated systems, has been pursued 33,34 to test specific interactions with antimicrobial peptides and understand their selectivity. 35,36 The reliability of such simulations depends on the accurate parameterizations of lipids and proteins, which need to be validated against experimental data. Moreover, the two descriptions must be consistently integrated into the force fields used, i.e., be derived with the same parameterization procedure. Different approaches to the problem are possible, which resulted in the development of multiple force fields suitable for simulations of biomolecules: for example, the CHARMM 37−39 and AMBER 40 force fields are parameterized from quantum mechanics calculations, while GROMOS96 41 is calibrated to match global properties like the hydration free energy of chemical moieties. All of them have been constantly updated to meet the new experimental values available and more faithfully reproduce the different species involved.
However, it is a very difficult task to parameterize the constituents of a complex system so that all parameters are consistent with the rest of the force field and reproduce both the single-molecule observables and the collective behavior. In the present work, we consider the parameterization of phospholipids in the context of the GROMOS96 force field, 41 addressing some of the inconsistencies in the lipid head group parameters commonly used so far, particularly in consideration that these contribute to the description of recognition processes at the interface.
In the past, lipid simulations using the GROMOS96 force field suffered from difficulties involved in transferring the preexisting parameters, calibrated mainly for peptides in an aqueous environment, to the amphiphilic environment of the lipid assembly. This resulted in the failure to reproduce the membranes' behavior properly 42−44 and therefore a series of modifications were adopted, particularly in the choice of lipidspecific Lennard-Jones interactions 44−46 and partial charges. 42 In the light of recent reparameterizations of a set of choline moieties 47 and of phosphate-containing species, 48 we undertake the task of updating the parameters used for lipids, in particular, phosphocholines, as they contain both these chemical moieties. Within this work, we show that it is possible to integrate the recently computed partial charges within simulations while maintaining good agreement with the available experimental data. We also test the transferability of the new phosphate charges onto lipids without a choline head group, namely, phosphoethanolamine (POPE) and phosphatidylglycerol (POPG).
Most importantly, the new description of phosphocholine is consistent with the GROMOS96 parameterization philosophy, based on the decomposition of large molecules into smaller compounds and subsequently fitting their parameters to experimental hydration free energies. Together with adjustments to specific van der Waals potentials, we believe that the parameters presented here will contribute to improving the accuracy of the description of membrane−solvent and membrane−protein interactions. To this aim, we compared the available parameters with the one proposed in this work, simulating the interaction of an antimicrobial peptide with two model membranes, highlighting the differences in the mechanisms observed, and comparing them with the available experimental evidence.

METHODS
2.1. Background to Lipid Force Fields. The most recent iteration of the lipids' parameters commonly used in simulations with the GROMOS force field is the one by Poger and Mark. 44 They employed partial charges derived quantum-mechanically by Chiu et al., 42 combined with a modified repulsion between the choline methyl groups and the OM oxygen atoms in the phosphate with respect to the standard choline−OM one.
The original set of Chiu charges 42 was derived from ab initio Hartree−Fock self-consistent field calculations 52 and Mulliken population analysis. 53 Slight modifications were applied to make each individual charge group sum up to an integer value, following the GROMOS96 philosophy. Despite the resulting charges that differ substantially from the ones used for the same chemical groups in different chemical contexts, the GROMOS community employed this set as it gave results in closer agreement with the available experimental data.
The refinement of van der Waals parameters for aliphatic alkanes, together with the bond, bond angle, and torsional parameters for the ester groups, 49,50 prompted the reparameterization of the lipids' head group description: in particular, Chandrasekhar et al. 45 recomputed the head group torsional parameters from ab initio quantum-mechanical torsional profiles of each of the fragments composing the head group ( Figure 1).
The modifications above were included in the 53A6 version of the GROMOS force field. However, an additional change was necessary to match the experimental results. Therefore, Poger and Mark 44 introduced a change in the CH3 choline and OM repulsion. The C12 parameter (related to the Pauli repulsion) between the newly introduced atom-type CH3p (to

Journal of Chemical Theory and Computation
Article represent the united methyl atoms in the choline head group of lipids) and the oxygen-type OM, present in the phosphate group, was increased by a factor of 3.5. This modification was optimized and tested against experimental values, increases the spacing between individual lipids, and thus leads to the appropriate area per lipid (ApL). 51,54 The new atom-type CH3p has all of the characteristics of CH3, except for the bespoke parameter, i.e., the Lennard-Jones interactions involving OM. These Poger−Chiu parameters have been successful in reproducing membrane behavior and were used in many MD applications. 55−57 Later, Reif et al. 47 enhanced the methyl−methyl repulsion for both CH3 and CH3p in the 54A8 parameter set, which allowed for a decrease in the large repulsion value between CH3p and OM previously introduced 44 while still reproducing experimental values. The 54A8 parameter set contains two additional, nonlipid-specific, modifications important for this work: the choline Lennard-Jones parameters and partial charges, and the phosphate partial charges. The C12 Lennard-Jones repulsion term for the NL nitrogen atom type (present in the choline moiety) was increased to successfully prevent oversolvation. 47,54 To the same end, the +1 e total charge was evenly distributed over all five atoms, which resulted in a better approximation of the experimentally obtained hydration free energy in comparison to the 54A7 parameter set. Similarly, Margreitter et al. 48 calibrated the partial charges of four phosphate species and enhanced the reproduction of experimental data. The relevant phosphatecontaining species for this work is dimethyl-phosphate, a compound not directly present in force field versions prior to 54A8.
Another approach to lipid parameterization was proposed by Kukol, 46 namely, the use of the already available CH0 atom type for the ester carbons in place of the standard C atom type, in conjunction with the Chiu charges. This atom type, designed to describe a bare sp 3 carbon bound to four heavy atoms, has a repulsion energy term 10−40 times larger than a bare carbon bound to other atom types, enforcing a greater spacing between lipid molecules and thereby increasing the ApL. As this modification is also applicable in the absence of a choline head group and does not require the introduction of another atom type, this method can be used to parameterize POPE and POPG.
2.2. Parameterization Strategy. In an effort to enhance the consistency of the force field, we integrated the new partial charges for the choline and phosphate moieties [Reif− Margreitter (RM) charge set] into the lipid building blocks of GROMOS 54A8 so that the entire phosphocholine head group now follows the common GROMOS-like modeling approach (Figure 2). Only the partial charges of the ester groups remain as described in the Chiu set, a deviation from the canonical parameterization strategy necessary to match the experimental area per lipid values: Chandrasekhar et al. showed in ref 59 that the replacement of the ester charges with the standard ones for the ester moiety (parameterized to reproduce the experimental free energies of hydration of a series of alkane esters 60 ) resulted in a much smaller area per lipid, not compatible with the experimental values.  All are run for 500 ns and systems consisting of 512 lipid molecules (256 per layer), using a particle mesh Ewald (PME) long-range electrostatic scheme. b Charge set: Chiu from Ref 42, Reif−Margreitter (RM) as illustrated in the present work. c As a reference, the standard C12 parameter in 54A7/54A8 for CH3−OM is 4.44 × 10 −6 kJ mol −1 nm 12 . d The CH3−CH3 C12 parameters are 2.66 × 10 −5 for each parameter set.

Article
The introduction of the new head group charges required a refinement of the CH3p−OM Lennard-Jones repulsion, as the 54A8 value was set considering the original Chiu charges. Ideally, one would always try to keep the force field terms as much transferable as possible. Nevertheless, the complexity and anisotropic nature of some biological environments can be difficult to parametrize with single chemical groups, as the same chemical group can behave differently according to the context it is inserted in. Lipid systems are one of such examples, and to maintain the correct physical behavior of the system, we used specific C12 parameters for the CH3p−OM repulsion in phosphocholine lipid atoms. This allows for a more balanced description of the physicochemical properties of the lipid bilayer and a better match with the available experimental observables.
To prove the transferability of the new phosphate charges to other lipid species, which do not contain a choline head group (and thus an enhanced repulsion, which has an impact on the ApL), test simulations of a phosphoethanolamine (POPE) and a phosphoglycerol (POPG, Table 2) bilayer have been performed. These lipids have amine and glycerol head groups, respectively. The parameterization of both takes advantage of the Kukol approach 46 employing a CH0 atom for the ester moieties to enhance the repulsion between lipids. For POPE and POPG, simulations were run with the standard parameters from ref 33 (denoted as Piggot−Chiu in the present work) or with the updated phosphate partial charges (Supporting Information (SI) Table 2).
The evolution of simulation techniques seen in the recent years suggested two other changes in the simulation setup: first, the original set of parameters was designed to be used with a twin-range cutoff scheme and a reaction field long-range electrostatic contribution, 62 but the twin-range cutoff is no longer supported in the latest versions of the GROMACS software used for the present work. 63 Additionally, the PME algorithm 64 for long-range electrostatic treatment is currently the predominant method used for protein dynamics. In the context of unifying the two fields of protein and lipid simulations, we therefore opted for a PME long-range treatment, running a control simulation (on the DPPC bilayer) with a reaction field scheme to assess the impact of such a change (SI Table 1).
The other change we adopted in comparison to the earlier work was a larger system size. Due to computational limitations, the original parameterization was performed on a 128-lipid bilayer, 44 but recent advances allow for larger systems to be simulated and we therefore used membranes four times as large (512 lipids). This larger size allows to track larger undulations of the membrane, as the effect of periodic boundary conditions (PBCs) is less restrictive. Again, a control simulation on a 128 DPPC membrane has been run to test the relevance and the effect of this change (SI Table 1).
Finally, the improvements reached with the adoption of the new parameters are monitored through the comparison with experimental values, but it is useful to have benchmarks derived from other simulation experiments. For that purpose, we compare some key properties with the values obtained by the all-atom CHARMM36 force field. 37,65 Despite a thorough comparison is beyond the present work, it is relevant to observe whether the changes introduced by the new parameters are going in the direction of the outcomes proposed by other descriptions.
2.3. Simulation Systems. Seven pure lipid bilayers have been simulated, five of which contain phosphocholines, one phosphoethanolamine (POPE), and one phosphoglycerol (POPG), as described in Table 2. Every bilayer is formed by 512 lipids (256 per leaflet), generated by replicating an equilibrated 128-lipid system from the literature two times in the x and y directions (see Table 2).
Water molecules were added to reach a minimum distance of 7.5 nm between periodic copies of the membrane along the z-direction, with a ratio of 85−120 H 2 O per lipid. This distance is larger than the one used in the previous parameterization publications because we observed an enhanced undulatory behavior for larger membranes and therefore a higher distance is necessary to avoid interactions between periodic replicas in the z-direction.
2.4. Simulation Parameters. All simulations were run using the GROMACS software version 2016.3, 63,72,73 under periodic boundary conditions in a rectangular box. The temperature was maintained by coupling the membrane and the solvent independently to an external bath using the Berendsen thermostat 74 with a coupling time τ T of 0.1 ps, at the reference temperatures indicated in Table 2, which are above the gel−liquid phase transition temperature for each lipid. The pressure was kept at 1 bar with a semi-isotropic coupling using a Berendsen barostat, 74 applying isothermal 16:0/18:1c9 indicates that tail 1 has 16 carbons with no unsaturated bonds and tail 2 has 18 carbons with one unsaturated bond between carbons 9 and 10ester carbon counts as number 1.

Journal of Chemical Theory and Computation
Article compressibility of 4.5 × 10 −5 bar −1 and a coupling constant τ P of 0.5 ps. Covalent bond lengths of the lipids were constrained using the LINCS algorithm. 75 The geometry of the simple point charge water molecules was constrained using SETTLE. 76 A 2 fs time step was used, with a Verlet integration scheme. The PME 64 long-range treatment was applied to the electrostatic interactions beyond a 1.4 nm cutoff, and the reaction field scheme 62 control simulation was run with the same cutoff radius. A plain cutoff was used for van der Waals interactions, with a cutoff radius of 1.4 nm.
Each system was initially energy-minimized and then simulated at 50 K for 10 ps. Subsequently, the temperature was increased gradually over 500 ps until the final simulation temperature. The system was then simulated for 500 ns. The equilibration of the systems was monitored by examining the time evolution of the potential energy and the area per lipid: 200 ns is found to be sufficient to reach equilibration for all of the bilayers (SI Figure 2) so that the analysis has been performed over the last 300 ns of the production run, with frames stored every 100 ps. An overview of the simulations performed is given in Table 1 and SI Tables 1 and 2. 2.5. Analysis. To calibrate the lipid parameters, we used the observables listed below, as common practice in standard parameterization procedures. 61,77 2.5.1. Area per Lipid. For systems where the membrane is aligned to the xy plane, the area per lipid (ApL) can be computed from the product of the lateral dimensions of the simulation box divided by the number of lipids in one leaflet. As shown in SI Figure 2 for DPPC, after 100 ns of simulation, the ApL oscillates around a value with fluctuations of the same magnitude, indicating equilibration. To allow further time for local rearrangements, we restrict our analyses to the last 300 ns of the simulations.
The equilibration protocol was verified on the DPPC bilayer, repeating the computation of the ApL on two nonoverlapping time windows, specifically between 200 and 350 ns and between 350 and 500 ns. For all of the parameter sets, the two windows gave compatible values of the ApL, confirming the convergence of the simulations (SI Figure 3).
The above procedure is valid if the membrane is flat or has minor undulations only. To test this and verify that deviations from planarity are not influencing the results, the ApL was recomputed for DPPC taking into account membrane undulations according to the procedure outlined in ref 78. The differences with the values computed from the simulation box dimensions were between 0.20 and 0.46%, which is lower than the error derived from the standard deviation across the simulation for any of the area per lipid computed.
As such, our computations are of value in rating the results against experimental outcomes and/or to compare parameter sets, as a local measure would not significantly improve the comparison.
where k B is the Boltzmann constant, ⟨T⟩ and ⟨ApL⟩ are the ensemble averages of the temperature and the area per lipid, respectively, n L is the number of lipids in one leaflet, and σ ApL 2 is the variance of ApL.

Bilayer Thickness.
From the electron density profiles, the bilayer thickness can be evaluated in several ways and compared to the values from X-ray scattering experiments: the hydrophobic thickness (D HH ) is measured as the distance between the phosphorus peaks in the two layers, as these atoms have the highest electron density, while the Luzzati thickness (D B ) 61 is defined as where b z is the z-dimension of the simulation box and ρ W (z) is the volume fraction of water (vs other components) along z and normalized to 1 in the bulk water region where n W (z) is the time-averaged number of water molecules in a bin of width dz, V W is the specific volume of the water model used (taken from ref 79), and dV is the time-averaged volume of a slice. 2.5.4. Dipole Potential. The dipole potential along the zdirection (perpendicular to the membrane plane) can be computed from the charge density along z (ρ(z)) via a double integration 82 Several choices are possible for the two integration constants, 83 and for the present work, they are selected to set the dipole potential to zero in the middle of the bulk water region, at both sides of the membrane. 2.5.5. Deuterium Order Parameter of Lipid Chains. The deuterium order parameters S CD of the acyl chains for each lipid bilayer were calculated and compared between the different sets studied. S CD evaluates the average order of the lipid tails by measuring the orientation with respect to the bilayer normal of a carbon−hydrogen bond in a given position along the chain for each lipid in the bilayer. Their spread is evaluated according to the ensemble average As the GROMOS force field employs a united-atom representation, the tetrahedral positions of the hydrogens are constructed based on the neighboring carbons' positions. 58,80,81 2.5.6. Hydration of Head Groups. To estimate and compare the hydration of lipid molecules, we computed the distribution of the distances between the oxygen of water and the nearest lipid atom. For each individual chemical group, the distance between the water oxygen and the nearest atom within that group was calculated. A quantitative measure for hydration was obtained by integrating the distribution up to the first peak or second one (for phosphate and glycerol).
2.5.7. Orientation of Head Groups. We computed the orientation of the lipids' head groups as the distribution of the angle between the P−N vector (joining the phosphorus atom and the choline nitrogen) and the outward normal to the membrane. The orientation of the sn-1 and -2 carbonyl dipoles with respect to the bilayer normal has also been calculated.

Journal of Chemical Theory and Computation
Article 2.5.8. Lateral Diffusion. For each simulation, we extracted the trajectory of the phosphorus atom of every lipid in the top and bottom leaflets separately, removing the collective motion of the leaflet. These trajectories were used to compute the mean-square displacement (MSD) for each lipid as a function of time, discarding the first 200 ns of production. This figure was averaged over all of the lipids in the leaflet and, for a given interval of time, on all of the possible time windows of that length-fitting within the simulation time analyzed. The diffusion coefficient D was obtained from a linear fit of the average MSD profile, following the Einstein equation 84 The fit was performed discarding the first 50 ns of the profile, where the behavior is not linear, and the last 100 ns, where the poorer statistics leads to more noisy data. Coefficients obtained for the two leaflets were averaged to give the value reported. 2.5.9. Tilt Modulus. We computed the tilt modulus following the theoretical framework explained in ref 85. According to this, the angle θ a lipid forms with the local normal to the membrane follows the distribution where C is a normalization constant, k B is the Boltzmann constant, T is the temperature, and κ t l is the tilt modulus. This can be extracted from a fit of the distribution or, for computational reasons, from a fit of ln(P(θ)(sin θ) −1 ). The direction in which a lipid points is taken as the vector joining the center of mass of the terminal atoms of the tails and the center of mass of selected atoms in the head group. Specifically, the last three carbons of each tail are taken as the reference for the first group, and the phosphorus and the carbon from which the two tails divert for the second. The computation was performed using a dedicated python module 85 available on the openStructure platform. 86 2.6. Phase Transition. The set of new parameters performing best according to the previous observables was tested for sensitivity to temperature variations. A DPPC bilayer was chosen as the reference system and simulated at two additional temperatures: 303 and 333 K (SI Table 1), the first of which is below the experimentally determined liquid to gelphase transition temperature. 66−69 As DPPC has also been used to perform the other control simulations, we opted for this model membrane for consistency reasons.
Besides the standard analysis described before, the local area per lipid was computed using a Dirichlet tessellation 87 of the lipid tail positions projected onto the horizontal plane parallel to the membrane (one leaflet at the time). The tessellation divides the plane into polygons, each enclosing one tail position. Every polygon comprises the locations on the plane, which are closer to the position of the head enclosed by that polygon than to any other head.
Moreover, to quantify whether and how many lipids undergo face transition during the simulations, the regular packing of each of their chains was quantified by the hexagonal order parameter S 6 , as previously reported in the literature. 88 Specifically, a chain was represented by its position on the xy plane (parallel to the membrane surface), computed as the average x and y positions of its carbon atoms. For each chain j, the set of neighboring chains was defined as the ones within a 0.65 nm radius from j. Then, S 6 is defined as (8) with θ jk being the angle between the vector connecting j and k and the x axis (i is the imaginary unit). A chain is in the gel phase if it has an hexagonal order parameter larger than 0.72. 88

RESULTS AND DISCUSSION
In general, the parameters described in this work are shown to reproduce the available experimental target values well while, at the same time, are likely to allow for a better description of lipid−protein interactions, since the head groups are updated to the recent GROMOS force field. 3.1. Area per Lipid and Isothermal Area Compressibility Module. We report in Table 3 and Figure 3 the values of ApL for the simulation run. From such computations, it emerges that the increase of the CH3p−OM repulsion has a nonlinear effect on the area per lipid, as reported in ref 54. On the contrary, the comparison between simulation ID1 and the control ID0 for DPPC, which differ only in their partial charges, shows an almost identical ApL value (SI Figure 4). This suggests that the charge redistribution in the head group affects the global structure of the bilayer and the ApL less dramatically than the adopted value for the Lennard-Jones repulsion.
The comparison with the control simulation using a smaller membrane shows that larger systems allow for the evaluation of the ApL with a smaller error, as local fluctuations are averaged over a larger number of lipids. The standard deviation computed for the 128 lipids system is compatible with those reported in both the original 61 and a more recent publication, 89 in which the same system size was used.
The ApL from the simulation with a reaction field treatment for the long-range electrostatic term does not differ significantly from the one obtained with a PME treatment, being only slightly higher, which is in consistence with what was found in ref 90. Finally, the values found using the Chiu charge set and the 54A7 force field (ID1 in Table 1) are systematically lower than those obtained in the original publications, 44,61 despite employing the same charge set and force field, while a better agreement is shown with those obtained more recently by Reif et al. for DPPC. 54 We attribute this to the different versions of GROMACS used, as the integration algorithm has recently been updated, affecting the calculated properties. Moreover, the double-cutoff scheme is no longer supported, preventing a faithful reproduction of the simulation setup used in ref 61. The variability caused by these changes has been extensively investigated by Reißer et al. 89 and reflects the observed discrepancy between the present and previous results.
From the considerations above, we suggest parameter set RM/54A8_v1 as the one that best reproduces all of the tested lipids at once. For DOPC and POPC bilayers, however, parameter set RM/54A8_v2 performs slightly better: it must be noticed that these two species present unsaturated bonds along the tails, whose influence might not be fully represented by any of the parameter sets. Indeed, it has been suggested that only a polarizable force field would be able to correctly capture the dynamics of the hydrophobic region of the membranes, 91

Journal of Chemical Theory and Computation
Article taking in proper account the difference between saturated and unsaturated bonds.
For POPE and POPG, we resorted to the modification proposed by Kukol,46 i.e., the use of the CH0 atom type for the ester carbons (see Section 2). For both lipids, a good agreement with experimental ApL values could be achieved using the new partial charge parameters (Figure 3).
Along the same lines, when comparing the results with the ones obtained with the CHARMM36 force field in its original publication, 65 we find DOPC, presenting an unsaturated bond in each tail, to be the most diverging. In particular, CHARMM36 better captures the spacing between the lipids, enhanced due to the presence of the double bond, and we suspect that this is due to its all-atom description.
Results of the isothermal area compressibility calculations confirm the finding of refs 61 and 65 that K A values obtained from simulation are about 1.5−3 times larger than those measured experimentally (SI Table 3). This holds for all parameter sets tested. Set RM/54A8_v1 performs better than

Article
Chiu/54A7 for all of the lipids tested but DLPC, for which the results are equivalent. The overestimation of the compressibility is likely due to the underestimation of ApL fluctuations during dynamic simulations. The K A value computed for the small, 128 lipids, DPPC bilayer patch is smaller than the one computed for the 512 lipid ones (342 and 499 mN m −1 , respectively), as the small patch exhibits higher fluctuations of the ApL (see Section 3.1). It is thus evident that the size of the system plays a pivotal role in obtaining correct fluctuations and global properties.
3.2. Electron and Charge Density Profile. Across simulations with different parameters, the electron density qualitatively maintains the same profile for each phosphocholine lipid. In Figure 4, the density for parameter set 54A8_v1 and all of the lipids is displayed (SI Figures 10−13 show the same plot for the other parameter sets), while in SI Figures 5− 9, panel (b), the total and the phosphate group electron densities are shown for the five parameter sets tested, for one lipid at the time. The peak broadness shows a direct relationship with the packing density of the bilayer: larger ApL values correspond to a shallower profile of the density, due to fluctuations of the membrane along the z axis and to deeper penetration of water molecules into the bilayer.
The bilayer thickness was evaluated from the electron density profiles, as explained in Section 2. Our parameters are overall in better agreement with the Luzzati estimate of the thickness rather than the hydrophobic one, but altogether, these measurements (phosphate and Luzzati thickness) do not strongly discriminate between sets. In Table 4, the values for the Chiu/54A7 and RM/54A8_v1 sets are shown (see SI Table 4 for the complete results).
Further comparison of the dipole potential profiles, obtained from the charge density, shows how the RM/54A8_v1 charge set gives results closer to the ones obtained in all-atom simulations 82,92 (see SI Section 1 for a complete discussion).
3.3. Order Parameter of the Acyl Chains. For all of the lipids and parameter sets, S CD is lower than 0.25, which indicates that the tails are generally disordered and the membrane has not transitioned to a gel-like state, 93 even for the simulation with the lowest ApL. Comparing these different sets, simulations denoted by ID from 1 to 5 show a consistently decreasing S CD , in line with the increased area per lipid and decreased bilayer thickness. This indicates that when the lipid molecules are constrained in space, their tails tend to be stretched and ordered. The presence of unsaturated bonds in the DOPC and POPC lipids is captured, by all parameter sets, as a decrease in S CD at the positions related to those bonds. The main difference due to the introduction of the new charges is in the decreased order observed for the first and second carbon bonds of the sn-1 tail, which show S CD values smaller than the ones for the third carbon bond, while with the Chiu charges, a constant increase is observed with decreasing carbon index for tail sn-1.

Journal of Chemical Theory and Computation
Article Overall, the RM/54A8_v1 set is within the range of experimental values 94,95 ( Figure 5); in particular, it captures the low order of the first sn-2 carbon atom (numbered 2) well, while the Chiu/54A7 set presents closer values in the central region of the tails. However, it must be noticed that variability is found within the experimental data (see the different experimental values reported in Figure 5). Therefore, without aiming at a perfect fit to such a small pool of experimental data, we consider set RM/54A8_v1 as sufficiently accurate in representing the experimental findings, in particular, in better reproducing the regions in the vicinity of the head group, while the description of the hydrophobic core remains less accurate and subject to improvement.
3.4. Hydration of Head Groups and Glycerol/Carbonyl Moieties. The hydration of functional groups of lipids is a key characteristic for both their dynamics and potential interactions with other molecules, such as proteins. From the distribution of distances between water oxygens and the nearest atom of various lipid groups, it emerges that the new partial charges modify the hydration profile of the lipid head group ( Figure 6 shows the comparison between parameter sets for DPPC and SI Figures 17−20 for the other lipid bilayers).
The choline major peak at 0.38 nm and the phosphate one at 0.30 nm are higher and sharper when employing the RM charges rather than the Chiu ones, reflecting an increased average hydration of these two moieties. Additionally, for the simulations run with the RM charges, the choline profile does not display the first, low intensity, peak obtained with the Chiu set at 0.28 nm: indeed, the charges of choline and the modification of the C12 Lennard-Jones repulsion for the NL atom type introduced in parameter set 54A8 were optimized to successfully prevent oversolvation, repelling water from its core. 47,54 The profiles of the other components are partially influenced, as well. For the RM charges, the second peak for glycerol increases its value and the two ester peaks have more similar values between them ( Figure 6, panel (c), and SI Figures 17−20), which is consistent with deeper water penetration.
To quantify the observed differences, the hydration profiles were integrated up to the first peak or the second one in the case of phosphate and glycerol (SI Table 5). The results show that the average number of water molecules around the choline group is higher for the RM/54A8_v1 set than for the Chiu/ 54A7 one by one water molecule. This seems to contrast with the increased hydrophobicity of the newly parameterized choline moiety; however, this might partially be explained by the changed orientation of the head groups (see Section 3.5) and by the new parameterization of phosphate, 48 which accounted for the hydrogen bond potential of the most solvent-accessible atoms, leading to a better solvation of the head group in comparison to the Chiu/54A7 and Chiu/54A8 sets.
The integration up to the second peak of the distribution of distances between the water oxygens and any lipid head group atom gives values between 12 and 17 water molecules per lipid, which is in agreement with the experimental range of 10− 20. 96−99 Again, parameter set RM/54A8_v1 results in more hydrated head groups (about one water molecule more for each lipid) with respect to Chiu/54A7. Notably, the average number of water molecules increases, as expected, for the simulations resulting in a larger ApL (RM/54A8_v2 and RM/ 54A8_v3). The trends above are confirmed by solventaccessible surface area values, which are higher for the choline head groups described by the RM charge set with respect to the Chiu one, while the values are closer between parameter sets for the phosphate and glycerol moieties, which are more deeply buried (SI Figure 21).
The increased hydration might be of relevance when simulating interactions with peptides and proteins. Moreover, as shown in a recent comparison between different lipid force fields, 65 the Chiu/54A7 parameter set results in a slightly less hydrated head group with respect to the CHARMM36 37 and Lipid14 101 force fields; therefore, the new set of parameters achieves values closer to them.
3.5. Orientation of the Head Groups and Carbonyl Moieties. The orientation of the head groups, defined by the angle of the P−N vector with the outward bilayer normal, is similar for all of the lipids within the same parameter set (see SI Figure 22, top row). This indicates that the nature of the tails does not strongly affect the behavior of the head group, which is to be expected. Comparing different sets for DPPC ( Figure 7 and SI Figure 23), it emerges that with the Chiu charges, the distribution of P−N angles is bimodal, with preferred values around 60 and 90°, while the new charges restrict the motion to the 60°configuration. Recent experimental data support a value around 60°(see refs 100 and 102), as opposed to 90°as reported previously. 103 It is noteworthy that this property was not part of the calibration process, i.e., the agreement with the experimental

Journal of Chemical Theory and Computation
Article observables in ref 102 is most likely due to a more accurate description of the solvation of the choline and phosphate moieties. Simulations performed by Botan et al. 104 confirm that smaller angles with respect to the membrane normal are caused by a higher level of head group hydration, which is in line with conclusions from the previous section. This difference in the predominant configuration of the lipids' head group will most probably influence the interaction with proteins or peptides approaching the interfacial region, providing a different binding recognition landscape.
The orientation of the sn-1 and sn-2 carbonyl dipoles with respect to the bilayer normal is again similar across different lipids (SI Figure 22, middle and bottom rows). The introduction of the RM charges has a small effect on these dipoles, as a result of the spatial rearrangement of the nearby head group. The most probable value for CO1 is shifted from 110 to 120°(Chiu vs RM charges), while the one for CO2 from 135 to 150°.
3.6. Lipid Lateral Diffusion. To correctly reproduce the membrane and its functions, its dynamical characteristics are as important as its structural ones. To address this, lipid lateral diffusion can be measured and compared against experimental data. Lateral diffusion is influenced by the area per lipid, with a tighter packing preventing larger displacements but is not solely determined by it.
Lateral diffusion coefficients (D) measured from simulations are shown in Figure 8. As anticipated, the set with largest ApL (54A8_v3) presents the highest values; however, parameter set 54A8_v1 gives significantly higher diffusion coefficients than those obtained with the Chiu/54A7 set, despite the values of ApL being similar. SI Figure 24 depicts a comparison of the diffusion coefficient of DPPC between ID0 and ID1, which differ only in the partial charges of the head groups. It confirms that the RM charges (ID0) allow for more mobility of the lipids with respect to Chiu ones (ID1), independent from all other modifications to the force field.
Regarding the simulation conditions, the use of a reaction field scheme increases the mobility by 34%, whereas the size of the patch decreases it by a small but significant amount (19%; see SI Figure 24). It is known that periodic boundary conditions affect the evaluation of lipid diffusion; 105,106 therefore, the larger the system simulated, the more accurate the reproduction of the experimental values. However, the change in the D due to the electrostatic treatment and the patch size, taken in absolute terms (i.e., a difference of about 0.4 and 0.2 μm 2 s −1 , respectively), are small in comparison with the effect due to the adoption of the new parameters (between 2 and 6 μm 2 s −1 ).
The comparison with experimental values is challenging due to the fact that different experimental techniques report values, which are an order of magnitude apart. Poger et al. gave an overview of this variability for DPPC bilayers in Table 2  The plot shows that the tilt modulus κ t l varies between 3 and 5 × 10 −20 J nm −2 . In simulations resulting in larger ApL (e.g., parameter set 54A8_v3), the lipids are in a less dense environment and can better accommodate changes in their orientations resulting in a lower tilt modulus (the tilt modulus gives the energy necessary for tilting the lipids per unit area).
The comparison with the experimental values is very good for DMPC and DPPC, while it is poorer for DLPC and very poor for DOPC and POPC, which harbor unsaturated bonds in the tails. Results from the CHARMM36 simulations show more variability between different phosphocholines, but a similar if not lower agreement with the experiment (for example, for DOPC). In general, comparing the results together, we think that we achieved a sufficiently qualitative agreement with the previous computational literature.
The discrepancy with experiments (both for our results and for the ones from ref 111) is likely due to computational limitations: as the tilt is retrieved from an ensemble

Journal of Chemical Theory and Computation
Article distribution, larger and longer simulations are more likely to give a better result. Moreover, as briefly mentioned at various points in the manuscript, artifacts arising in the hydrophobic regions of lipids (such as the suboptimal modeling of unsaturated tails) will probably only be resolved using a polarizable description.
3.8. Phase Transition Behavior. The previous analysis points to parameter set 54A8_v1 as the one that best reproduces the experimental properties for each of the lipids simulated. Therefore, we test this set further to assess its ability to reproduce the change in lipid behavior under different temperature conditions. As mentioned in Section 2.6, we use as test system the DPPC bilayer patch.
The comparison between simulations at 323 and 333 K shows that the global area per lipid increases with temperature, consistently with what was expected. Parameter set 54A8_v1 captures the increase in lipid spacing, with a slight underestimation of ApL at 333 K with respect to the experiments (for 54A8_v1 and experiments, respectively: 0.624(6) vs 0.631(13) nm 2 at 323 K and 0.634(5) vs 0.650(13) nm 2 at 333 K). The experimental ApL is measured at a different temperature with the same experimental setup. 16 The simulation at 303 K shows the formation of a patch of ordered lipids, suggesting that the parameters can reproduce different phases: two nucleation sites for the gelification process are observed in both the upper leaflet and lower leaflet, in nonmatching positions, and the gel front extends over time.
To classify the phase a lipid belongs to, we computed the hexagonal order parameter S 6 for each chain from the last frame of the simulation (time point 400 ns). Figure 10 shows the position of the chains on the xy plane, for each leaflet separately, color-coded by S 6 : the regions where S 6 is larger correspond to a densely packed area with a quasi-hexagonal lattice. In particular, the center of the ordered patches has S 6 values larger than 0.72 (last two colors of the scale), i.e., it can be classified as a gel. Overall, 20% of chains have undergone this transition within the time simulated.
Averaging over all of the lipids, S 6 at 303 K is 0.45. As a comparison, we computed this average quantity on the last frame of the simulations performed at 323 and 333 K, finding 0.28 and 0.27, respectively (with only six and three chains above the gel threshold of 0.72).
A hexagonal order can be obtained when the tails are well ordered and parallel to each other, standing in a vertical straight conformation (SI Figure 25 shows a detail of a wellordered gel patch). We thus compute the S CD order parameter of the acyl chains averaging separately over the lipids for which at least one chain has an S 6 value larger than the 0.72 threshold (168 lipids overall) and for the others. The last 100 ns of the simulation time was used. These values are compared to the average S CD from the simulations at 323 and 333 K. Figure 11 shows highly ordered tails for the membrane simulated below the transition temperature, for both the gel and nongel lipids (classified according to the S 6 threshold). This suggests that the full patch is undergoing a phase transition, but the completion of the process is not seen due to the short simulation time scale. As a comparison, tails at 323 and 333 K are much more disordered, with a slight decrease in order with increasing temperature.
Finally, we computed the local area per lipid chain from a Dirichlet tessellation of the same set of points used to calculate the hexagonal order parameter. The average values over the gel-phase tails (multiplied by 2) give an ApL of 0.438 ± 0.038 nm 2 , while the remaining of the chains have widely spread values, correlated to their S 6 parameter, giving an average of 0.57 ± 0.18 nm 2 . The values found for the gel patches (at 303 K) are close to the experimental outcomes by Nagle et al. of 0.473 nm 2 at 293 K 112 and 0.479(2) at 297 K. 113 The value computed from simulations is smaller likely because it is computed only over the tails perfectly packed in a hexagonal lattice.
Altogether, these results prove that the newly developed parameters can successfully reproduce the gel phase when a lipid patch is simulated below the phase transition temperature.
3.9. Transferability to POPE and POPG. As mentioned above, the areas per lipid values of POPE and POPG simulations, where the phosphate partial charges have been replaced with the RM values and, in the case of POPE, the amine partial charges have been updated according to the 54A8 force field, are in good agreement with the available experimental data.
For POPE, a slightly enhanced hydration is obtained from the update of the phosphate charges (from 5.7 to 6.4 water Figure 10. Hexagonal order parameter S 6 for lipid acyl chains computed on the last frame of a DPPC bilayer simulation using 54A8_v1 parameters. Each point corresponds to the average position of the carbon atoms of the respective chain on the xy plane. Thus, every lipid is represented by two points in these plots. Solid black lines denote the boundaries of the simulation box and chains of the periodic images (used for the computation of S 6 boundaries) are shown faded out. Colors from red to blue denote an increasing S 6 value: the last two indicate gel-phase lipids. Figure 11. Order parameter for the acyl chain sn-2 for a DPPC bilayer simulated at a different temperature. The average is performed including both leaflets. For the simulation at 303 K, the lipids were split in two groups according to their hexagonal order parameter S 6 and the acyl chain order parameter S CD computed for each of them.

Journal of Chemical Theory and Computation
Article molecules per lipid with experimental values between 4 and 7; 114,115 see SI Table 7) with similar results in terms of thickness D B (3.89 nm for Chiu and 3.92 for RM set, with an experimental value of 4.13 nm; 116 SI Table 6 and SI Figures 26  and 27). Overall, these results confirm the transferability of the new phosphate charges to different types of phospholipids.

INTERACTION WITH PROTEINS
The adoption of the updated parameters enhances the consistency with the GROMOS parameters for protein simulations. To test how this affects the simulations of peptides interacting with a membrane, we performed additional simulations of a small antimicrobial peptide on the surface of two different model membranes.
The peptide selected is bovine lactoferricin (PDB code 1LFC). It has a length of 25 amino acids and adopts a βhairpin conformation in solution, with many aromatic hydrophobic residues on one side and charged amino acids distributed all over. 117 This peptide is antimicrobial and therefore found to preferentially bind bacterial membranes NA denotes no binding observed in the time simulated (100 ns). For POPC/OIII simulated with parameter sets Chiu/54A8, the binding time is in parentheses as LFC approaches the membrane but maintains a 0.4 nm distance (±0.02), which is higher than the threshold chosen to define a binding event.

Journal of Chemical Theory and Computation
Article versus mammal ones, as shown by NMR experiments on LFC subsequences (namely, LFC 4−9 118 and LFC 4−14 119 ). One can model the bacterial membrane by a mixture of zwitterionic and anionic lipids. The latter are characteristic of the cell wall of both Gram-positive and -negative bacteria. 120,121 In this study, we selected the mixture DLPC/DLPG with a 3:1 ratio that has been used to elucidate the antimicrobial activity of lactoferricin-derived peptides. 122 As for the mammal membrane description, we used POPC as it has been often used in molecular dynamics simulations with this purpose. 123−125 Despite being rather simple, these or similar model membranes have often been used in experiments to test, among others, the effects of antimicrobial peptides upon binding. 14 Molecular simulations can shed light on the differences in the binding process of LFC to antimicrobial and mammal model membranes. Our parameterization should then reflect a sensible difference in the binding behavior for these two cases.
The exact binding mechanism of LFC to a membrane is not fully understood, and a number of experimental papers have hypothesized binding modes for the interactions of this peptide with model membranes. A mutation study in LFC 1−15 suggests that Trp residues anchor the peptide to the membrane as the antimicrobial activity of the peptide was retained only when Trp was mutated in equally hydrophobic amino acids. 126 However, the role of Trp seems to be different in other antimicrobial peptides, where they reside at the lipid−water interface and form hydrogen bonds with the moieties nearby. 127 As experiments on the full-length peptide (25 amino acids) have not yet been reported and the full sequence contains additional charged and hydrophobic residues, it remains unclear whether this additional region would change the aforementioned binding mechanism. We therefore decided to use molecular dynamics simulations to elucidate molecular determinants in discriminating the binding of the peptide to mammal and bacterial membranes.
In order to not be biased by the initial configuration adopted in the simulation, we performed multiple simulations with different initial orientations of the peptide relative to the membrane. The hairpin main axis was aligned to the membrane plane and the peptide rotated around this axis in steps of 90°, leading to four different starting orientations named OI, OII, OIII, and OIV (SI Figure 29). This allows different segments of the sequence (and thus amino acids with different chemical characteristics) to face the membrane in the initial positioning. The initial minimum distance between the peptide and the lipids was set to 2 nm. The simulation length was 100 ns each, sufficient to see the binding process in all of the control cases.
The simulations have been performed for the proposed RM/ 54A8_v1 parameter set and the Chiu/54A8 one (control cases) to compare with the most recent set available in GROMOS and highlight the difference of the newly parameterized lipid head groups.
To quantify the outcome of the simulations, we monitored the time at which the peptide binds (always irreversibly) to the membrane as the time at which the minimum distance between the peptide and the membrane is below 0.3 nm ( Table 5). The cutoff was chosen, analyzing the configurations after LFC bound to the membrane, which resulted generally to stabilize around a minimum distance of 0.25 nm. The minimum distance was computed every 100 ps, and a running average was applied with a 10 frame window. Additionally, the insertion depth of each amino acid in the membrane has been calculated as the difference between the z position of the lowest atom of the amino acid and the average of the maximum z coordinate of the five lipids closest to it. Table 5 shows the different binding times for LFC against a mixed DLPC/DLPG or pure POPC membrane patch. For the mixed, anionic membrane, the new parameters favor a slower and weaker binding process. Indeed, with parameter set Chiu/ 54A8, the peptide is quickly sequestrated by the lipids due to the opposite charge interaction. This favors an unspecific binding, dependent on the sequence facing the membrane in the initial configuration ( Figure 12 and SI Movie 1). In Figure  13, the average insertion in the membrane after the binding is plotted for each amino acid: Chiu/54A8 favors a deep insertion of differently charged residues for different runs. The RM/54A8_v1 simulations produce a less inserted configuration of LFC and a more consistent protrusion of the hinge region out of the membrane, i.e., the central stretch of amino acids between Met and Leu ( Figure 13), as three out of the four simulations (all but OIII) show this behavior.
The angular orientation of the peptide around its axis has been computed as the angle formed with the z axis by the backbone carbon and nitrogen bonded via a hydrogen bond (amino acids 7 with 19 and 9 with 17), confirming that the new set of parameters allows for more freedom in the reorientation of the initial configurations, while the previous one tends to keep them close to the original configuration (SI Figure 30). Additionally, the new set of parameters seems to favor the reorientation of the peptide as to face the Trp residues toward the membrane surface in three out of the four simulations, in contrast with the results from the previous parameterization. This preference for the interfacial region is a known mechanism in the membrane binding of Trp-and Argrich peptides. 127 When simulating a pure POPC membrane (here considered as a mammal membrane model), the resulting binding poses obtained with the Chiu/54A8 and different initial conditions are consistent among each other; in particular, the three amino acids Lys12−Leu13−Gly14 located at the hinge of the hairpin promote the insertion, while the terminal region stays exposed in solution. Therefore, parameter set Chiu/54A8 discriminates between the two membranes as it suggests a weaker binding to the mammal one. However, with the parameter set RM/ Figure 13. Average insertion depth of each amino acid after binding of the LFC peptide as per Table 5. The zero value is the top of the membrane plane so that a negative depth means insertion into the membrane. Some reference amino acids are displayed at the bottom.

Journal of Chemical Theory and Computation
Article 54A8_v1, three out of the four simulations result in no binding at all, in agreement with experimental findings. 118 The remaining simulation (OI) shows a quick binding event, promoted by the terminal regions as observed for three out of four simulations with the DLPC/DLPG 3:1 membrane.
The results above highlight that the new parameters show a membrane-binding process less dependent on the initial conditions, allowing for a dynamical rearrangement of the protein at the membrane interface. This comes at the expenses of a longer sampling time needed to observe binding events for most of the configurations chosen. Future work will focus on systematic comparisons of available peptide−membrane simulations with other parameterizations and on longer simulated times.
The difference between the behavior on a model bacterial or mammal membrane is more pronounced for the new parameters, and this is consistent with the selective antimicrobial action of the peptide and its low hemolytic activity. 118,119 Overall, we think that these new parameters show promising characteristics for the simulation of membrane−peptide interactions within the GROMOS force field, particularly for the study of interfacial absorption.

CONCLUSIONS
In this work, we present a reparameterization of a range of phospholipids in the context of the GROMOS force field, taking advantage of recent optimizations reported for key chemical groups in these molecules. The effect of the newly adopted head group partial charges has been tested extensively to ensure that they match experimentally observable characteristics of lipid bilayers. In parallel, we tested the effect of the van der Waals repulsion between the choline methyl groups and the phosphate oxygens, as it was modified by Poger et al. to reproduce the experimental area per lipid values while using the partial charges derived by Chiu et al. A summary of the updated parameters and simulation conditions is available in SI Table 10.
The work proves that the new charges are suitable to describe all of the phosphocholine bilayers tested, matching the experimental values as successfully as the previous parameter set. The major advantage of the Reif−Margreitter set lies in the partial charges of the head group, which are derived by applying the GROMOS parameterization philosophy rather than quantum mechanics calculations, thereby providing a description, which is more consistent with the parameters adopted for other biomolecules such as proteins within this force field. By using the updated partial charges for the choline (more hydrophobic) and phosphate (more hydrophilic) groups, the parameters also show a better reproduction of the average head group orientation, which was recently reassessed by experiments. The value of the Lennard-Jones repulsion term found to best reproduce the experimental values is the one in set 54A8_v1, which is set to a value in between that of the 54A7 and 54A8 parameter sets.
In the Reif−Margreitter parameter set, only the partial charges of the ester groups remain as described in the Chiu charge set. Preliminary work has been started to test the influence of the ester charges in combination with the new ones for the head group but, in accordance to what was previously found by Chandrasekhar et al. (ref 59), the replacement of the ester charges with the standard ones for the ester moiety resulted in values of area per lipid too low with respect to the experimental findings. As mentioned previously, it is possible that this discrepancy with the rest of the force field can only be avoided by adopting a polarizable force field. 91 However, in the absence of further sophisticated changes to the force field parameterization, we are confident that the proposed parameters are a major step forward in the description of lipid head groups and they should enable improved modeling of the interaction of lipids with water and other soluble molecules.
The new phosphate partial charges have been proved to transfer well to other phospholipids not presenting a choline head. For those lipids, the Kukol modification, which takes advantage of a different atom type for the ester carbon, is adopted to obtain the correct area per lipid.
Finally, the performance in reproducing some specific peptide−membrane interactions was tested. In this respect, the new parameter set shows significant differences with respect to the latest Chiu/54A8 set: it better discriminates the binding of an antimicrobial sequence on a bacterial versus a mammal membrane. Additionally, it favors a weaker and more dynamic binding, which is less biased from the initial conditions of the simulations.
In conclusion, we believe that the new Reif−Margreitter charge set together with the GROMOS 54A8_v1 parameter set is a major improvement on the previous iteration of the GROMOS lipid force field and should be particularly suited for protein−membrane systems, such as studies including small antimicrobial peptides, which rely on an accurate peptide− membrane recognition.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00509.
Dipole potential computation Additional tables: control simulations; non phosphocholine simulations; compressibility, thickness, hydration profile for all parameter sets; summary of simulation parameters Additional plots: control simulations results; comparison of properties across parameter sets; electron density, order parameter, hydration profile and head group orientation for parameter sets not presented in the main text; phosphocholine moieties SASA; snapshot of DPPC patch partially in gel phase; non phosphocholine simulations results; scheme and measure of 1LFC orientation with respect to the membrane (PDF) Simulation of LFC on a bacterial membrane for both parameter sets Chiu/54A8 (left half) and RM/54A8_v1 (right half), starting from orientation OII (Movie 1) (MP4) POPG are available at http://fraternalilab.kcl.ac.uk/ wordpress/biomembrane-simulations. Note that these parameter files include the post-translationally modified residues specified for the GROMOS 54A8 force field earlier. 128

■ ACKNOWLEDGMENTS
The