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 X-ray 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 pre-existing 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 lipid-specific 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

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 ).

Figure 1

Evolution of the lipid parameters in the GROMOS force field. References—Chiu1995: ref (42); Chandrasekhar2001, 2002, 2003: refs (45, 49, 50); Poger2010: ref (44); Kukol2009: ref (46); Schmid2011: ref (51); Reif2012: ref (47); Margreitter2017: ref (48).

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 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.[ref55−ref57]

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 phosphate-containing 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 sp3 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.

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 ). 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 (56) 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[57]) resulted in a much smaller area per lipid, not compatible with the experimental values.

Figure 2

Partial charges for the phosphocholine head groups and the glycerol and ester moieties in the Chiu[42] scheme (left) and the one tested in the current work (right). Red font denotes values that have been changed between the two. Atoms belonging to the same charge groups are enclosed by the same dashed polygon.

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.

Aiming at this, and to disentangle the effect of charge parameterization versus the CH3p–OM repulsion, we tested three different values of such Lennard-Jones parameter with the new charges while control simulations were run using the Chiu partial charges and the GROMOS 54A7 or 54A8 parameter set for each lipid (Table ). The phosphocholines tested are 1,2-lauroyl-sn-glycero-3-phosphocholine (DLPC), 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC), 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC), and 2-oleoyl-1-palmitoyl-sn-glycero-3-phosphocholine (POPC), which have different tail lengths and numbers of unsaturated bonds as in previous works[44,58] (Table ).

Table 1

Table of Simulations for Phosphocholine Bilayersa

simulations of phosphocholine lipids
sim	chargesb	FF	CH3p–OM C12c (kJ mol–1 nm12)	CH3p–CH3p C12d (kJ mol–1 nm12)
1	Chiu	54A7	1.58 × 10–5	2.66 × 10–5
2	Chiu	54A8	6.93 × 10–6	6.48 × 10–5
3	RM	54A8_v1	1.10 × 10–5	6.48 × 10–5
4	RM	54A8_v2	1.58 × 10–5	6.48 × 10–5
5	RM	54A8_v3	4.50 × 10–5	6.48 × 10–5

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.

Charge set: Chiu from Ref (42), Reif–Margreitter (RM) as illustrated in the present work.

As a reference, the standard C12 parameter in 54A7/54A8 for CH3–OM is 4.44 × 10–6 kJ mol–1 nm12.

The CH3–CH3 C12 parameters are 2.66 × 10–5 for each parameter set.

Table 2

Details of the Systems Simulated: Lipid Name, Tail Composition, Initial ApL (and Reference from Which the Initial Coordinates Are Taken), Simulation Temperature, and Gel–Liquid Phase Experimental Transition Temperature

lipid bilayer systems
lipida	tailsb	ApL0 (nm2)	TMD (K)	TC (K)
DLPC	12:0/12:0	0.632[44]	303	276.4[63−66]
DMPC	14:0/14:0	0.616[58]	303	296.9[63−66]
DOPC	18:1c9/18:1c9	0.649[58]	303	255.7[63−66]
POPC	16:0/18:1c9	0.638[58]	303	270.5[63−66]
DPPC	16:0/16:0	0.631[58]	323	314.2[63−66]
POPE	16:0/18:1c9	0.568[33]	313	299.3[67]
POPG	16:0/18:1c9	0.602[33]	303	268.1[68]

DLPC: 1,2-lauroyl-sn-glycero-3-phosphocholine, DMPC: 1,2-dimyristoyl-sn-glycero-3-phosphocholine, DOPC: 1,2-dioleoyl-sn-glycero-3-phosphocholine, POPC: 2-oleoyl-1-palmitoyl-sn-glycero-3-phosphocholine, DPPC: 1,2-dipalmitoyl-sn-glycero-3-phosphocholine, POPE: 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine, POPG: 1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-(1′-rac-glycerol).

Example: 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.

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 ) 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,[59] but the twin-range cutoff is no longer supported in the latest versions of the GROMACS software used for the present work.[60] Additionally, the PME algorithm[61] 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,62] 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.

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 . 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 ).

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 H2O 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.

Simulation Parameters

All simulations were run using the GROMACS software version 2016.3,[60,69,70] 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[71] with a coupling time τT of 0.1 ps, at the reference temperatures indicated in Table , 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,[71] applying isothermal 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.[72] The geometry of the simple point charge water molecules was constrained using SETTLE.[73] A 2 fs time step was used, with a Verlet integration scheme. The PME[61] long-range treatment was applied to the electrostatic interactions beyond a 1.4 nm cutoff, and the reaction field scheme[59] 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 and SI Tables 1 and 2.

Analysis

To calibrate the lipid parameters, we used the observables listed below, as common practice in standard parameterization procedures.[58,74]

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 (75). 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.

Isothermal Area Compressibility Module

Following the protocol in ref (58), we computed the isothermal area compressibility module (KA) from the fluctuations of the ApL values according towhere kB is the Boltzmann constant, ⟨T⟩ and ⟨ApL⟩ are the ensemble averages of the temperature and the area per lipid, respectively, nL is the number of lipids in one leaflet, and σApL2 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 (DHH) 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 (DB)[58] is defined aswhere b 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 regionwhere nW(z) is the time-averaged number of water molecules in a bin of width dz, VW is the specific volume of the water model used (taken from ref (76)), and dV is the time-averaged volume of a slice.

Dipole Potential

The dipole potential along the z-direction (perpendicular to the membrane plane) can be computed from the charge density along z (ρ(z)) via a double integration[77]Several choices are possible for the two integration constants,[78] 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.

Deuterium Order Parameter of Lipid Chains

The deuterium order parameters SCD of the acyl chains for each lipid bilayer were calculated and compared between the different sets studied. SCD 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 averageAs the GROMOS force field employs a united-atom representation, the tetrahedral positions of the hydrogens are constructed based on the neighboring carbons’ positions.[55,ref77,ref78]

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).

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.

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[79] in two dimensionsThe 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.

Tilt Modulus

We computed the tilt modulus following the theoretical framework explained in ref (80). According to this, the angle θ a lipid forms with the local normal to the membrane follows the distributionwhere C is a normalization constant, kB is the Boltzmann constant, T is the temperature, and κtl 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[80] available on the openStructure platform.[81]

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 gel-phase transition temperature.[63−66] 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[82] 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 S6, as previously reported in the literature.[83] 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, S6 is defined aswith θ 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.[83]

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.

Area per Lipid and Isothermal Area Compressibility Module

We report in Table and Figure 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.

Table 3

Average Area per Lipid (in nm2) over the Last 300 ns of Simulations for Phosphocholine Bilayersa

ApL (nm2)
ID	charges/FF	DLPC	DMPC	DOPC	POPC	DPPC
1	Chiu/54A7	0.608(4)	0.591(4)	0.600(5)	0.604(4)	0.616(4)
2	Chiu/54A8	0.626(5)	0.612(4)	0.623(4)	0.623(4)	0.635(5)
3	RM/54A8_v1	0.631(4)	0.616(5)	0.625(6)	0.629(5)	0.638(5)
4	RM/54A8_v2	0.652(5)	0.643(6)	0.649(5)	0.650(5)	0.657(5)
5	RM/54A8_v3	0.690(6)	0.684(6)	0.690(6)	0.687(6)	0.693(5)
0	RM/54A7					0.603(5)
RF	Chiu/54A7					0.603(4)
small	Chiu/54A7					0.594(11)
experimentalb		0.608–0.632	0.589–0.660	0.674–0.725	0.643–0.683	0.570–0.717
CHARMM36[37]		0.644(4)	0.608(2)	0.690(3)	0.647(2)	0.629(3)

The number in parentheses is the standard deviation of the last digit. All simulations are run at 303 K, except for DPPC (run at 323 K). Analogous values for POPE and POPG are reported in SI Tables 6 and 8.

We report the maximum and minimum values of a review of experimental results given in Table 1 of ref (74). Only values referring to the temperature simulated are considered.

Figure 3

Area per lipid obtained for the five sets of parameters and seven lipid species. Error bars are the standard deviation over the 300 ns analyzed. Dashed lines indicate the range of experimental values from Table in refs (74) and (33). For the plot reporting POPE and POPG values, the black dashed lines refer to POPE and the blue one to POPG.

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[58] and a more recent publication,[84] 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 (85).

Finally, the values found using the Chiu charge set and the 54A7 force field (ID1 in Table ) are systematically lower than those obtained in the original publications,[44,58] 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 (58). The variability caused by these changes has been extensively investigated by Reißer et al.[84] 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,[86] 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 ). For both lipids, a good agreement with experimental ApL values could be achieved using the new partial charge parameters (Figure ).

Along the same lines, when comparing the results with the ones obtained with the CHARMM36 force field in its original publication,[62] 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 (58) and (62) that KA 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 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 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 ). It is thus evident that the size of the system plays a pivotal role in obtaining correct fluctuations and global properties.

Electron and Charge Density Profile

Across simulations with different parameters, the electron density qualitatively maintains the same profile for each phosphocholine lipid. In Figure , 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.

Figure 4

Electron density profiles of the hydrated DLPC, DMPC, DOPC, DPPC, and POPC bilayers (total) and of their individual components (Cho: choline, PO4: phosphate, gly + carb: glycerol and carbonyl groups, CH2: methylenes of the acyl chains, CH: CHdCH groups in the oleoyl chains, CH3: terminal methyls of the acyl chains) for simulation ID3 (54A8_v1 force field, Reif–Margreitter charges).

The bilayer thickness was evaluated from the electron density profiles, as explained in Section . 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 , the values for the Chiu/54A7 and RM/54A8_v1 sets are shown (see SI Table for the complete results).

Table 4

Bilayer Thickness for Phosphocholine Bilayers, Derived from the Electron Density Profiles (Example in Figure ) According to the Phosphate or Luzzati Methodsa

hydrophobic thickness DHH (nm)
ID	charges/FF	DLPC	DMPC	DOPC	POPC	DPPC
1	Chiu/54A7	2.83	3.59	3.05	3.30	4.30
3	RM/54A8_v1	2.72	3.48	2.89	3.22	4.06
experimentb		3.08	3.44–3.60	3.53–3.71	3.70	3.42–3.83

All simulations were run at 303 K, except for DPPC (323 K).

Values from ref (44) and Table 2 in ref (58).

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[77,87] (see SI Section 1 for a complete discussion).

Order Parameter of the Acyl Chains

For all of the lipids and parameter sets, SCD is lower than 0.25, which indicates that the tails are generally disordered and the membrane has not transitioned to a gel-like state,[88] even for the simulation with the lowest ApL. Figure and SI Figures 14–17 display the computed values for specific parameter sets, and SI Figures 5–9, panel (c), show a cross-parameter comparison for each lipid. Comparing these different sets, simulations denoted by ID from 1 to 5 show a consistently decreasing SCD, 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 SCD 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 SCD 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.

Figure 5

Deuterium order parameter SCD profiles of the sn-1 (solid curves) and sn-2 (dashed curves) fatty acyl chains of hydrated DLPC, DMPC, DOPC, DPPC, and POPC bilayers calculated from simulations ID1 (54A7 force field, Chiu charges) and ID3 (54A8_v1 force field, Reif–Margreitter charges). The SCD values are averaged over all of the lipid sn-1 and -2 acyl chains in the systems (proS hydrogen only). Experimental values Douliez1995 from ref (89) and Petrache2000 ones from ref (90).

Overall, the RM/54A8_v1 set is within the range of experimental values[89,90] (Figure ); 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 ). 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.

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 shows the comparison between parameter sets for DPPC and SI Figures 17–20 for the other lipid bilayers).

Figure 6

Distribution of the distance between the water oxygen and the nearest lipid head group atom for simulation DPPC. Cho: choline, PO4: phosphate, Gly: glycerol, CO1 and CO2: carbonyl groups at the sn-1 and sn-2 positions.

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 , 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 ) 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.[91−94] 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 solvent-accessible 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,[62] the Chiu/54A7 parameter set results in a slightly less hydrated head group with respect to the CHARMM36[37] and Lipid14[ref95] force fields; therefore, the new set of parameters achieves values closer to them.

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 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 (95) and (96)), as opposed to 90° as reported previously.[97]

Figure 7

Distribution of the P–N, CO1, and CO2 angles with respect to the outward normal to the bilayer.

It is noteworthy that this property was not part of the calibration process, i.e., the agreement with the experimental observables in ref (96) is most likely due to a more accurate description of the solvation of the choline and phosphate moieties. Simulations performed by Botan et al.[98] 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°.

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 . 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.

Figure 8

Lateral diffusion coefficient of DLPC, DMPC, DOPC, DPPC, and POPC bilayers for different parameter sets.

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;[99,100] 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 μm2 s–1, respectively), are small in comparison with the effect due to the adoption of the new parameters (between 2 and 6 μm2 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 of ref (74) and observed that values span from 0.5 to 50 μm2 s–1. In this view, the values obtained in the present work for DPPC are well within the range, regardless of the parameter set chosen. However, we report experimental values from ref (101) obtained through pulsed-field gradient nuclear magnetic resonance as a guide. Additionally, we report the values obtained with CHARMM36 in ref (62). The CHARMM36 benchmarks are present only for two of the phosphocholines analyzed in this work and show that the values obtained with this force field span a broader range. The consistently low values of D computed with the different GROMOS parameter sets in this work are in agreement with what was found in the literature.[102,103]

Tilt Modulus

We report in Figure the values of κtl obtained for each of the phosphocholines considered and each parameter set tested. For comparison, we plot the experimental values obtained by Nagle et al. in ref (104) and the results from simulations using the CHARMM36 force field.[105] Given that the data show quite a large spread in their values depending on the actual experimental setup used, for the comparison, we selected values, which were all obtained under the same conditions, for both the experiments and the computational results.

Figure 9

Tilt modulus κtl computed from the distribution of lipid tilt angles along the last 300 ns of the trajectories. The results are compared with the experimental values from ref (104) and the ones obtained (with the same procedure as employed here) from CHARMM36 simulations in ref (105).

The plot shows that the tilt modulus κtl 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 (105)) is likely due to computational limitations: as the tilt is retrieved from an ensemble 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.

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 , 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) nm2 at 323 K and 0.634(5) vs 0.650(13) nm2 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 S6 for each chain from the last frame of the simulation (time point 400 ns). Figure shows the position of the chains on the xy plane, for each leaflet separately, color-coded by S6: the regions where S6 is larger correspond to a densely packed area with a quasi-hexagonal lattice. In particular, the center of the ordered patches has S6 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.

Figure 10

Hexagonal order parameter S6 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 S6 boundaries) are shown faded out. Colors from red to blue denote an increasing S6 value: the last two indicate gel-phase lipids.

Averaging over all of the lipids, S6 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 well-ordered gel patch). We thus compute the SCD order parameter of the acyl chains averaging separately over the lipids for which at least one chain has an S6 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 SCD from the simulations at 323 and 333 K. Figure shows highly ordered tails for the membrane simulated below the transition temperature, for both the gel and nongel lipids (classified according to the S6 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.

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 S6 and the acyl chain order parameter SCD computed for each of them.

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 nm2, while the remaining of the chains have widely spread values, correlated to their S6 parameter, giving an average of 0.57 ± 0.18 nm2. The values found for the gel patches (at 303 K) are close to the experimental outcomes by Nagle et al. of 0.473 nm2 at 293 K[106] and 0.479(2) at 297 K.[107] 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.

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 molecules per lipid with experimental values between 4 and 7;[108,109] see SI Table 7) with similar results in terms of thickness DB (3.89 nm for Chiu and 3.92 for RM set, with an experimental value of 4.13 nm;[110] 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.[111] This peptide is antimicrobial and therefore found to preferentially bind bacterial membranes versus mammal ones, as shown by NMR experiments on LFC subsequences (namely, LFC4–9[112] and LFC4–14[113]).

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.[114,115] 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.[116] As for the mammal membrane description, we used POPC as it has been often used in molecular dynamics simulations with this purpose.[117−119] 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 LFC1–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.[120] 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.[121] 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 ). 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

Binding Time of Lactoferricin (LFC) Peptide to the Model Membranes in Examinationa

binding time (ns)
	DLPC/DLPG 3:1				POPC
	OI	OII	OIII	OIV	OI	OII	OIII	OIV
Chiu/54A8	0.5	3.8	6.3	2.8	13.9	21.9	(65.0)	13.9
RM/54A8_v1	2.8	75.2	62.0	9.6	1.0	NA	NA	NA

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.

Table 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 and SI Movie 1). In Figure , 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 ), as three out of the four simulations (all but OIII) show this behavior.

Figure 12

Final configurations of the simulations of LFC on a DLPC/DLPG 3:1 membrane, starting from four different initial orientations OI–OIV. OII, OIII, and OIV are obtained from OI with an anticlockwise rotation of, respectively, 90, 180, and 270° along the main axis (x axis in the top right panel). LFC is colored by residue type (blue charged, green hydrophilic, white hydrophobic); phosphorus atoms are shown in golden beads. The terminal and hinge regions of 1LFC are indicated (TER, hinge), together with GLN7 as a red dot to help the visualization. The insets shows a cartoon representation of the initial and final configurations, highlighting the positive patches as blue dots.

Figure 13

Average insertion depth of each amino acid after binding of the LFC peptide as per Table . 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.

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 Arg-rich peptides.[121]

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/54A8_v1, three out of the four simulations result in no binding at all, in agreement with experimental findings.[112] 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.[112,113] 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 (56)), 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.[86] 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.