Compatibility of GROMOS-Derived Atomic Parameters for Lipopolysaccharide Membranes with the SPC/E Water Model and Alternative Long-Range Electrostatic Treatments Using Single Nonbonded Cutoff and Atom-Based Charge Schemes

Abstract

Recent developments of GROMACS v.2016 ceased to support methodological approaches used in the development and validation of the GROMOS force field. We investigated the performance of a previously developed extension of the GROMOS force field for lipopolysaccharides to reproduce the structural dynamics of bacterial outer membrane (OM) using a single cutoff for nonbonded interactions and atom-based charge truncation. We further compared this setup for use with reaction field (RF) or particle mesh Ewald (PME) approximations in the presence of simple point charge (SPC) and extended simple point charge (SPC/E) water models. We find that the OM structural dynamics is well conserved in all simulated conditions, reproducing the available experimental data within the measurement uncertainty. The SPC/E model induces a small increase in OM fluidity, and when combined with the RF correction, shows a decrease in water orientation at the membrane surface. The present simulations support the compatibility of the GROMOS-derived lipopolysaccharide (LPS) parameters for use with single cutoff and atom-based charge schemes, either with RF or PME approximations, in GROMACS v.2016. Both SPC and SPC/E water models are suitable for use, but usage of SPC/E combined with the reaction field correction needs to be further investigated.

Keywords:
atom-based charge truncation; reaction field approximation; particle mesh Ewald; single cutoff; Verlet algorithm; GROMACS v.2106; bacterial outer membrane

Introduction

The Gram-negative outer membrane (OM) is a critical barrier for protection against xenobiotic agents such as antibiotics and host innate immune molecules. For this reason, the OM drastically limits the intracellular access of antibiotics, and plays a critical role in antibiotic resistance by these pathogens. The bacterial OM consists of an inner monolayer containing phospholipids, primarily phosphatidyl ethanolamine, and an outer monolayer largely containing a single glycolipid species: lipopolysaccharide (LPS). This highly negative molecule contains three regions (Figure 1):^3-5 the innermost Lipid-A region is a glucosamine diphosphate with 4-7 fatty acids attached to it, which anchors the LPS to the phospholipid inner leaflet. The Lipid-A is appended to the rough core region containing 8-12 variable sugars with 3-8 phosphorylation sites. These phosphate substituents are essential to membrane stability because their negative charge allows neighboring LPS molecules to be cross-linked by divalent cations. Another structural signature of this region is the highly conserved occurrence of the octasaccharide 3-deoxy-D-manno-oct-2-ulosonic acid (KDO). The rough core is connected to the O-antigen region formed by a variable repetition of 3-5 sugar units, typically hexoses. The O-antigen determines the serotype specificity of LPS molecules whereas Lipid-A is the activator of immune response in mammalians. Environmental stimuli or genetic mutations trigger bacteria to produce LPS lacking the O-specific chain, the rough LPS. This phenotype has core oligosaccharides of varying length and sequences (chemotypes). The chemical modifications associated to bacterial membrane remodeling are often expressed by LPS conformational transitions with a complex dependence on temperature and cation concentration.^6-8 The cation nature is also relevant since divalent but not monovalent cations lead to LPS aggregation and increase bacterial resistance.^9-13

Figure 1
(a) Chemical structure of the pentaacylated LPS molecule from P. aeruginosa^1,2 and (b) molecular representation of the simulated bacterial outer membrane. The outer leaflet is composed of LPS and the inner leaflet is made of DPPE. Carbohydrate moieties are shown in magenta and acyl chains from Lipid-A and DPPE are shown in grey. Carbon, oxygen and phosphorus atoms are colored grey, red and tan. Ca²⁺ counterions are represented by van der Waals spheres in green. For clarity, water molecules were removed.

The experimental characterization of the structural dynamics of lipid bilayers at atomic or near-atomic resolution is not easily achievable even for systems composed of a single type of lipid.^14-17 This is chiefly due to their great fluidity and lack of mid- and long-range order. The chemical complexity and dynamical polymorphism of LPS molecules makes the structural characterization of outer membranes an even greater challenge. Molecular dynamics (MD) simulation remains one of the primary tools to investigate the properties of model biological membranes at the atomic or near-atomic level. It is also a robust method to probe the spatial and temporal rearrangements of organization of lipid bilayers during membrane fusion, fission, poration and so on. Yet, the reliability of MD simulations relies on the accuracy of the interatomic potentials, or force field, to reproduce lipid structure and phases under variable conditions of temperature, hydration, pH and compositional heterogeneity for instance. Therefore, the true aptness of a given force field to correctly reproduce experimental data can only be assessed through systematic comparison of simulated properties against experimental measurements. It is, however, fundamental to distinguish direct from indirect “experimental” data since the latter is inferred from primary data and relies on assumptions of a given model (for a review see reference 15). Examples of primary measurements for lipid bilayers are X-ray and neutron form factors from scattering experiments, ², H¹³C and ³¹P relaxation rates and quadrupolar splittings in nuclear magnetic resonance (NMR) spectroscopy, and fluorescence intensity in fluorescence spectroscopy whereas cross-sectional area per lipid, lipid bilayer thicknesses, NMR bond order parameter and lipid diffusion coefficient are examples of indirect measurments.¹⁵ Direct experimental measurements for OM are still challenging. However, recent progresses in the development of supported asymmetric bilayers made OM models amenable to characterization via structural techniques.^18-20

The first MD simulation of a bacterial outer membrane was performed for the rough LPS of Pseudomonas aeruginosa using the GLYCAM93²¹ and AMBER95²² atomic parameters to describe carbohydrate and lipid moieties, respectively.²³ The simulated bilayer consisted of 16 LPS units in the outer leaflet, 40 DPPE (1,2-dipalmitoyl-sn-glycero-3-phosphoethanolamine) units in the inner leaflet and a total of 104 Ca²⁺ to neutralize the system. Each LPS molecule was assigned a protonation state equivalent to a charge of -13 e per molecule, and the system was simulated for 1 ns. These parameters were subsequently tested via several simulations of the rough²⁴ and smooth (rough LPS + O-antigen)²⁵ LPS membranes from P. aeruginosa and applied to the first simulation of an outer-membrane protein embedded in a bacterial outer membrane.^26,27 However, the GLYCAM93 atomic parameters underwent several developments from the original parameter set used in the first bacterial outer membrane simulations.²⁸ The latest GLYCAM06 parameter set ceased the use of 1-4 nonbonded scaling factors and the assignment of partial charges to every atom for which quantum mechanics (QM) molecular electrostatic potentials were obtained by excluding aliphatic hydrogen atoms during the fitting of partial charges.²⁸ In addition, the latest GLYCAM06 in conjunction with the PARM94 parameter set for van der Waals tackled some critical limitations of the GLYCAM95 set of great relevance for LPS membranes such as poor solvation behavior.^22,28 Therefore, an extension of the GLYCAM06 force field was specifically developed for LPS membranes. The new set of parameters for lipopolysaccharide molecules expanded the GLYCAM06 repertoire of monosaccharides to include phosphorylated N- and O-acetylglucosamine, 3-deoxy-D-manno-oct-2-ulosonic acid, L-glycero-D-manno-heptose and its O-carbamoylated variant, and N-alanine-D-galactosamine. The new atomic parameters were initially validated through 1 μs MD simulations of the rough LPS membrane of P. aeruginosa PA01 made of 72 LPS and 180 DPPE molecules.²⁹ The new assignment of protonation states to the LPS molecules led to a charge of -8 e per molecule as expected in neutral pH. The new parameter set has been further expanded to include four new chemotypes: rmlC, galU, LPS Re and Lipid-A (Figure 1).^30,31 Furthermore, these new atomic parameters have been used to examine the effect of different cations on the stability of the LPS membranes,³² the binding mechanism of antimicrobial peptides and proteins to LPS membranes,^33-35 the effect of chemical remodeling on the structural dynamics and electrostatic properties of OM,^30,31 and to validate coarse-grained parameter sets for LPS from P. aeruginosa.^36,37 Concurrently, extensions of the GROMOS force field have been developed for the penta- and hexa-acylated forms of Lipid-A, LPS Re and Lipid-A modified with 4-amino-4-deoxy-L-arabinose from P. aeruginosa^38,39 and for the Rd1 LPS from E. coli,⁴⁰ followed by an extension of the CHARMM force field for the smooth LPS from E. coli.⁴¹ The GROMOS extension for chemotypes of P. aeruginosa relied on the parameter set 53A6_GLYC⁴² to treat LPS carbohydrate moieties and on the parameter set 53A6 for the acyl chains^43,44 as adapted and validated for Lipid-A.³⁸ The GROMOS 53A6_GLYC parameter set contains dihedral potential corrections for hexopyranoses with all the remaining atomic parameters being the same as for parameter set 54A7.⁴⁵ Additionally, atomic charges for new chemical groups and new torsional potentials (e.g. for the P-glycosidic linkage and all C-C-C-O dihedrals) in the LPS chemotypes were calculated to maintain the compatibility with previous versions of GROMOS force fields.^46-48 Atomic parameters for van der Waals interaction terms were retrieved from GROMOS 45A4/53A6 functional form for carbohydrates.⁴⁶ Currently, complete sets of atomic parameters for most common carbohydrates found in LPS from different species are available within the major GLYCAM/AMBER, CHARMM and GROMOS force fields. The CHARMM-GUI web portal hosts CHARMM parameters for various bacteria species, including 37 Lipid-A types, 52 core oligosaccharide types, and 304 O-antigen polysaccharide types.⁴⁹

Despite important advances in experimental front, molecular-level information on LPS structure remains very limited and insufficient to validate the large number of LPS chemotype models.^10,50,51 Furthermore, lipid simulations are very sensitive to small variations in methodological schemes used to alleviate cutoff errors such as those involved in the treatment of long-range nonbonded interactions.^15,52 This is so because weak interactions, critical for the structure and dynamics of lipid aggregates, add up to significant contributions in lipid simulations. Force fields have been usually parameterized and validated for use with specific approximations, and for this reason their accuracy and overall quality depends critically on the strictness with which these approximations have been implemented and tested in available codes. The GROMOS-derived parameters for LPS^38,39 have been ported and validated using the GROMACS software,^53,54 which has recently undergone major developments with important implications for the reproducibility and transferability of the GROMOS force field within the framework of its validated settings (e.g. charge group, reaction-field and twin-cutoff schemes).^55,56 The latest version of GROMACS (v.2016)⁵⁷ replaced the leapfrog algorithm by a variation of the Verlet-I/r-RESPA as the multiple-step alternative method to separately integrate fast- from slow-evolving interactions. The advantage of the new Verlet method⁵⁷ is its time-reversibility and symplecticity, but at the cost of great dependency of the update frequency on the nature of the system.⁵⁵ On the other hand, the implementation of the new Verlet method in GROMACS precludes the use of the charge group scheme associated with GROMOS force field parameterization. An additional setback is the availability of GPU acceleration in GROMACS exclusively with the Verlet scheme.⁵⁷

We have performed a series of MD simulations of our reference OM model described with an extension of the GROMOS force field for LPS,^32,33 which includes dihedral potential corrections for hexopyranoses,^38,39,42 using two approximations to treat long-range electrostatics interactions (generalized reaction field and particle mesh Ewald) in the presence of simple point charge (SPC)⁵⁸ and extended simple point charge (SPC/E)⁵⁹ water models. The benchmark system is composed of rough LPS units in the outer leaflet and DPPE in the inner leaflet simulated using the conditions consistently adopted in the parameterization of the GROMOS force field:^{15,43,52,60-62} leapfrog integration with a time step of 2 fs, charge group scheme and the reaction field approximation for the treatment of long-range electrostatic interactions. Nonbonded interactions, however, were evaluated using a single cutoff scheme (instead of the twin cutoff used with the GROMOS force field) with a cutoff distance of 1.4 nm updated every 5 steps together with the atom pairlist. We report thereafter our findings and conclusions.

Experimental

MD simulations were performed for a bacterial OM model comprised of 72 pentaacylated rough LPS molecules in the outer leaflet, 180 DPPE (1,2-dipalmitoyl-sn-glycero-3-phosphoethanolamine) molecules and 288 Ca²⁺ counterions to neutralize a charge of -8 e per LPS unit. The GROMOS parameter set 53A6_GLYC^42,63 was used for the LPS carbohydrate moieties whereas the acyl chain parameters were taken from the GROMOS 53A6 force field,^43,44 adapted and validated for Lipid-A.³⁸ The GROMOS 53A6_GLYC parameter set contains dihedral potential corrections for hexopyranoses with all the remaining atomic parameters being the same as for parameter set 54A7.⁴⁵ Atomic charges for new chemical groups and new torsional potentials (e.g. for the P-glycosidic linkage and all C-C-C-O dihedrals) in the LPS chemotypes were calculated to maintain the compatibility with previous versions of GROMOS force fields.^46,47 Atomic charges were calculated at the HF/6-31G* level of theory, followed by a restrained electrostatic potential (RESP) fitting.⁶⁴ Atomic parameters for van der Waals interaction terms were retrieved from GROMOS 45A4/53A6 functional form for carbohydrates.⁴⁶ The bacterial OM was simulated in two different water models, the SPC⁵⁸ and SPC/E,⁵⁹ and subjected to two distinct long-range electrostatic approaches, the generalized reaction field (RF)⁶⁵ and particle mesh Ewald (PME).⁶⁶ The initial configurations of the simulated systems had the same number of solute atoms, but differed in the number of solvent molecules (i.e., SPC = 30699, SPC/E = 30507) in order to maintain the initial box dimensions equal in all systems. The topologies and atomic parameters used for the bacterial OM are presented in the Supplementary Information (SI, Table S1: atomic charges, dihedral torsional potentials, potentials for bond stretching, bond-angle bending and improper dihedral deformation). The simulated systems are presented in Table 1.

Table 1
Simulated systems and averaged structural properties calculated from the rough LPS bilayer simulations using different water models, long-range electrostatics (LRE)^a and charge cutoff (CC)^b schemes. Values were averaged over the full simulation length

System	Water	LRE	AL / nm2	DHH / nm	Volume / nm3	SCD	D / (× 10-5 cm2 s-1)	Coordination numberc
								Ion-OW	OW-OW	ION-PO4	OW-PO4
1	SPC	RF	1.53 ± 0.01	3.14 ± 0.01	2.39 ± 0.02	0.20 ± 0.04	2.3 ± 0.1	5.4 ± 0.3	2.732 ± 0.002	2.5 ± 0.3	0.57 ± 0.02
2	SPC	PME	1.54 ± 0.01	3.14 ± 0.02	2.40 ± 0.03	0.21 ± 0.04	2.4 ± 0.1	5.3 ± 0.2	2.944 ± 0.003	2.5 ± 0.2	0.56 ± 0.01
3	SPC/E	RF	1.51 ± 0.01	3.14 ± 0.01	2.35 ± 0.02	0.17 ± 0.04	1.5 ± 0.2	5.5 ± 0.2	2.728 ± 0.002	2.4 ± 0.2	0.57 ± 0.01
4	SPC/E	PME	1.50 ± 0.01	3.19 ± 0.02	2.37 ± 0.03	0.18 ± 0.04	1.5 ± 0.2	5.4 ± 0.2	2.943 ± 0.003	2.4 ± 0.2	0.56 ± 0.01

• ^a
  Long-range electrostatics treatment: RF: reaction field; PME: particle mesh Ewald;
• ^b
  charge cutoff scheme: charge- or atom-based truncation. A_L: area per lipid; D_HH: bilayer thickness; S_CD: deuterium order parameters for the acyl chains; D: diffusion constant;
• ^c
  coordination numbers are listed for atom pairs where ion corresponds to the Ca²⁺ metal, OW is the oxygen atom of the respective water molecule and PO₄ is the phosphate group. SPC: simple point charge; SPC/E: extended simple point charge.

MD simulations were performed in the NPT ensemble with a time step of 2 fs. Bond lengths within the solute and the geometry of water molecules were constrained using LINCS and SETTLE algorithms, respectively.^67,68 Initial velocities were taken from a Maxwell-Boltzmann distribution at 310 K and 1 bar. A separated coupling was used to control the temperature of solute and solvent with counterions to a velocity-rescaling thermostat with a relaxation time of 0.4 ps.⁶⁹ The pressure was maintained at 1 bar in the lateral and normal directions with respect to the bilayer by weakly coupling to a semi-isotropic pressure bath with a relaxation time of 0.1 ps and isothermal compressibility of 4.5 × 10^-5 bar^-1.⁷⁰ Long-range electrostatics interactions were treated via two approaches. The RF correction was applied to the electrostatic interactions beyond a cutoff of 1.4 nm in conjunction with a permittivity dielectric constant of 66.^71,72 The PME correction was used beyond a cutoff of 1.2 nm with a fourth order interpolation of charges on a 0.16 nm Fourier spacing. The nonbonded pair lists were updated every 5 steps using a single cutoff scheme for atom pairs beyond a cutoff of 1.4 and 1.2 nm in the RF and PME simulations, respectively.⁵⁶ Further, the charge group and atom-based cutoff schemes were used with the RF and PME simulations, respectively.⁵⁶ Trajectories were recorded at every 2 ps. The systems were previously equilibrated for 100 ns, when the area per lipid molecule converged, and subsequent data production was performed for additional 100 ns. All simulations were performed using GROMACS 4.6.7 software⁵³ running on both CPU and GPU platforms. Property analysis were performed with the GROMACS 4.6.7⁵³ and in-house software.

Results and Discussion

Assessment of MD-derived structural and dynamical properties for the bacterial OM model indicates that differences in response to varying methodological treatments are not significant (Table 1). Time averaged quantities calculated from the MD simulations are within a maximum error percentage of ca. 4% with respect to the corresponding quantity average for all four systems. However, these small differences can persist for some properties throughout the simulated times and reveal subtle trends as discussed thereafter.

Experimental area per lipid (A_L) values of 1.42, 1.56 and 1.82 nm² have been reported for the hexaacylated rough LPS in the liquid crystalline phase from X-ray diffaction.^10,73,74 The average A_L calculated from the MD simulations are within the spectrum of experimental values. MD simulations of OM with the SPC water shows slightly larger A_L values compared to the SPC/E model (Figures 2a-2b). Neutron scattering measurements⁵¹ and previous MD simulations^29,30 have shown that LPS bilayers are appreciably more hydrated than phospholipid bilayers, with water molecules penetrating deeper in the hydrophobic region. For this reason, the difference in A_L for the OM in the two water models is expected to result from increased hydration by the SPC water. The SPC and SPC/E models have identical geometry parameters, but differ in the atomic charges assigned to oxygen and hydrogen atoms which has a significant impact on the physical-chemical properties of the models (e.g. molecular dipole, diffusion and density).^59,75,76 Density profiles show a lower water density for the SPC model compared to the SPC/E (Figure 3), consistent with the larger dipole and more favorable potential energy E_pot of water molecules in the SPC/E model (Table 2). The average number of SPC water molecules inside the LPS leaflet is slightly higher than the SPC/E: 876 (2.9%) for SPC-RF, 907 (3.0%) for SPC-PME compared to 748 (2.5%) for SPC/E-RF and 815 (2.7%) for SPC/E-PME. In addition, the higher diffusion constant of the SPC model, and consequently higher water motional dynamics, may lead to some degree of swelling of the carbohydrate region, with increase of A_L (Table 2, Figure 4).³⁰ Furthermore, the average carbon-deuterium order parameter S_CD are higher for the SPC pair of simulations compared to the SPC/E (Table 1). The more ordered, less fluid OM in the SPC simulations has an increased acyl chain packing, which leads to higher A_L values.

Figure 2
(a) Time-dependent and (b) time-averaged area per lipid A_L; (c) time-dependent and (d) time-averaged bilayer thickness D_HH. Quantities were calculated from MD simulations: SPC-RF in black; SPC-PME in red; SPC/E-RF in green; SPC/E-PME in blue.

Figure 3
Mass density profiles for selected chemical groups. Simulated systems are (A) SPC-RF; (B) SPC-PME; (C) SPC/E-RF and (D) SPC/E-PME. Water molecules are represented in blue (a), LPS molecules in red (b), DPPE in green (c), phosphate groups in orange (d) and Ca²⁺ ions in purple (e).

Figure 4
(a) Average diffusion coefficients of water molecules as function of the time and (b) time average for the whole simulation. Simulated systems are SPC-RF in black, SPC-PME in red, SPC/E-RF in green and SPC/E-PME in blue.

Table 2
Force field parameters and selected properties of the simple point charge (SPC)⁵⁸ and extended simple point charge (SPC/E),⁵⁹ water models at 300 K and identical simulation conditions⁷⁵

Property	SPC	SPC/E
qO / e	-0.82	-0.8476
qH / e	0.41	0.4238
ΘHOH / degree	109.47	109.47
rO-H / nm	0.100	0.100
ε / (kJ mol-1)	0.65017	0.65017
σ / nm	0.31656	0.31656
Epot / (kJ mol-1)	-41.3 ± 0.16	-46.4 ± 0.18
Dipole / Debye	2.27	2.39
D / (× 10-5 cm2 s-1)	4.5 ± 0.08	2.8 ± 0.02

• qO: oxygen atomic charge; qH: hydrogen atomic charge; Θ_HOH: HOH angle; r_O-H: OH distance; ε: depth of the potential energy minimum; σ: separation of potencial energy minimum; E_pot: total potential energy of water molecules represented by the two different models. See reference 75 for more details; D: diffusion constant.

The OM average thickness is well conserved among the simulated systems (Figures 2c-2d). The system SPC-PME presents a negligible increase of bilayer thickness (D_HH) in ca. 0.05 nm. For comparison, the standard length of the covalent bond between carbon atoms is 0.154 nm. The S_CD calculated from the MD simulations shows a higher orientational dynamics of carbon-carbon bonds with respect to the bilayer normal in the SPC simulations (Figure 5a). The calculated S_CD are inversely related to the A_L. All S_CD profiles reproduce the characteristic pattern of saturated hydrocarbon chains with low order at the beginning and end of the acyl chains and increase in the middle (Figures 5b-5e). Attenuated total reflectance measurements coupled to Fourier transform infrared (ATR-FTIR) spectroscopy was previously used to estimate the orientational behavior of LPS molecules from the peak position x_s of the symmetric stretching vibration of the methylene group around 2850 cm^-1.^73,77-79 Accordingly with this technique, Lipid-A from P. aeruginosa exhibits S_CD of 0.28 in the liquid-crystalline phase.⁸⁰ The calculated S_CD indicates higher fluidity of the simulated bilayers with respect to the experimental estimate, but nonetheless the calculated S_CD values are representative of a bilayer in the liquid-crystalline La phase. However, S_CD has been shown to converge slowly, ca. 300 ns or longer, in previous MD simulations of LPS bilayers,²⁹ and longer simulations will be required before conclusive comparisons can be made.

Figure 5
Deuterium order parameters S_CD for LPS acyl chains. (a) Time averages over the last 50 ns of simulation for SPC-RF, SPC-PME, SPC/E-RF and SPC/E-PME. Acyl chain order parameter for (b) SPC-RF; (c) SPC-PME; (d) SPC/E-RF and (e) SPC/E-PME. Sn1 chain is represented in black, Sn2 in red and Sn3 in green.

The space- and time-averaged surface curvature angle (S_c) was also calculated for the simulated systems (Figure 6). Curvature angle distributions are within an interval of 0-35°, consistent with fully planar bilayers. The small variation in the curvature angle distribution from the SPC/E-RF simulation is too small to induce noticeable changes in membrane curvature; it only causes a small increase in the percentage of LPS molecules with curvature angles within 15-30° (Figure 6a). Likewise, the average curvature order parameter S_c is lower for SPC/E-PME, though it differs in only 1% from the largest S_c of SPC-PME. A perfectly flat surface has S_c = 1. Average water molecule orientation (S_water) with respect to the normal axis was also calculated for the simulated OM (Figure 7). The solvent has a random orientation in the bulk regions outside the OM, i.e., 0 to 2.5 nm and 9.5 to 12 nm, and becomes more ordered at the OM interface (Figure 7). The transition to a ordered regime occurs at the position of the phosphate groups from DPPE and Lipid-A phosphate group (Figure 7). Maximum S_water takes place inside the OM, in the region between 4 to 6.5 nm, which includes the hydrophobic region and the LPS inner-core region (Figure 7). S_water decreases steadily in the outer-core region until it reaches zero at the LPS leaflet surface. This profile is representative of all four systems. However, the time evolution of the S_water profiles differs as function of the water model and long-range electrostatics treatments. The initial and final average configurations from simulations SPC-RF and SCP-PME show similar S_water profiles, which are also representative of the intermediate configurations (Figures 7a-7b). These simulations also show a small increase in S_water with time. In contrast, SPC/E-RF and SPC/E-PME show greater differences with respect to the S_water profiles (Figures 7c-7d). This is most noticeable for S_water values at the phospholipid-water interface and inside the OM in the initial and final configurations. Whereas S_water decreases in a time-dependent manner for SPC/E-RF, it increases appreciably for SPC/E-PME (Figures 7c-7d). The apparent faster convergence of S_water for OM simulations with the SPC model may be due to the higher diffusion coefficient of this model (Table 2).^59,75 To this end, water diffusion coefficients D were calculated from the MD simulations (Figure 4). Unsurprisingly, the SPC water presents faster diffusion than the SPC/E. However, these diffusion coefficients were not responsive to the different treatments of long-range electrostatics interactions. Therefore, it appears that water diffusivity may explain only partially the S_water profiles for the SPC/E water. Since the only difference between the two water models are their charges and the resulting dipole, we assume that the dipole interactions with the bilayer will be the quantity influenced by the different long-range electrostatics treatment. This is currently under investigation.

Figure 6
(a) Surface curvature angle S_c distributions and (b) average curvature order parameters calculated from MD simulations. Systems are SPC-RF in black, SPC-PME in red, SPC/E-RF in green and SPC/E-PME in blue.

Figure 7
Orientational order S_water of water molecules at the outer membrane surface. Initial and final configurations are shown in black and red lines, respectively. Intermediate configurations sampled at every 5000 steps are shown in brown dashed lines. Simulated systems are (a) SPC-RF; (b) SPC-PME; (c) SPC/E-RF and (d) SPC/E-PME.

There has been a large body of work dedicated to the comparison of RF versus PME approximations for simulations of condensed-phase systems. The emerging consensus appears to be that PME is best suited for the treatment of periodic systems and for simulations of systems in large computational boxes with a high ion strength, which is effective in the screening of electrostatic interactions, and reduces the interactions between periodic solute copies. Nevertheless, in simulation conditions where periodic solute copies show unphysical interactions, the RF correction allows for a user-defined tuning of the lengthscale of allowed interactions.⁸¹ Previously, a systematic comparison of MD simulations of phosphatidylcholine bilayers using an extension of the GROMOS 53A6⁸² and the SPC model have not found any significant difference between structural and dynamical properties derived from simulations using the RF or PME truncation.⁵² This is consistent with our present findings for OM simulations, despite the difference in charge between the LPS (-8 e) and 1,2-dipalmitoyl-sn-giycero-3-phosphocholine (DPPC, zwitterionic) molecules. It is probable that charge-related artifacts in LPS-containing aggregates were minimized in previous and in the present simulations by the use of a large cutoff of 1.4 nm with the RF correction. However, the atomic parameters have lower performance when combined with the SPC/E water and the RF correction. It should be noticed that the length of the present simulations prevents a conclusive explanation of the effect of RF on the orientational dynamics of SPC/E at the OM surface. Nonetheless, our findings indicate the GROMOS-derived atomic parameters for LPS is compatible for use with the latest GROMACS v.2016 using PME to treat long-range electrostatic interactions (on GPU or CPU processing), nonbonded interactions using a single cutoff and atom-based charge scheme. A valuable bonus is that the present simulations on GPU cores are more than seven times faster than on CPU cores (Table S2, SI section).

Conclusions

Force fields for biomolecular simulations have often evolved in tandem steps. New parameter sets are developed against selected physical properties of a given number of small compounds and then characterized and validated for systems of interest. The choice of water model is crucial for the accurate reproduction of most structural and dynamical properties of biomolecules. Suitable water models for use as solvent in biomolecular simulation should be computationally efficient, reproduce the properties of bulk water satisfactory and be compatible with the force field employed to describe the solute interactions. The GROMOS force field has been parameterized for use with the SCP model which is the foundation on which the SPC/E model was built. The two models can accurately describe a wide range of water properties. However, the SPC/E model shows a dielectric permittivity and a diffusion coefficient close to the experimentally measured values which may be more suitable to reproduce the hydration dynamics of LPS membranes.

In praxis, force fields are developed and validated for use with specific approximations so that their accuracy and ability to reproduce experimental measurements depend on the reliability with which these approximations have been implemented and tested in available codes. The GROMOS-derived parameters for LPS^38,39 have been ported and validated using the GROMACS software with the GROMOS validation settings (e.g. charge group, reaction-field and twin-cutoff schemes). The latest GROMACS v.2016 no longer supports the methodological approaches used for GROMOS parameterization and validation. Therefore, we compared MD simulations of our reference outer membrane (OM) model, described with an extension of the GROMOS force field for LPS,^32,33 which includes dihedral potential corrections for hexopyranoses,^38,39,42 using the reaction-field (RF) and particle-mesh Ewald (PME) approximations to treat long-range electrostatics interactions in presence of the SPC⁵⁸ and SPC/E⁵⁹ water models. These simulations were performed with GROMACS v.2016 and GROMOS validation settings, and analyzed with respect to area per lipid (A_L), bilayer thickness (D_HH), density profiles, carbon-deuterium order parameters (S_CD), surface curvature angle (S_c), water orientational order (S_water) and diffusivity. The structural dynamics of the OM model was well conserved in the four simulation conditions, and were consistent with the available experimental data. However, the OM in SPC/E water exhibited small increase in fluidity as shown by higher A_L and lower S_CD values. The present simulations indicate the GROMOS-derived atomic parameters for LPS^32,33 is compatible for use with GROMACS v.2016 using PME to treat long-range electrostatic interactions (on GPU or CPU processing), nonbonded interactions using a single cutoff and atom-based charge scheme.

Supplementary Information

The GROMOS-compatible atomic parameters (atom type, bond, angle, and torsion parameters) and topologies for the pentaacylated rough LPS used in the OM model are available at http://jbcs.sbq.org.br as PDF file. Atomic parameters, topologies and atomic coordinates for the simulated OM models compatible with the GROMACS suite of programs are available for download at dqfnet.ufpe.br/biomat. Benchmarks for the simulations run on GPU and CPU cores are also available as supplementary information.

Acknowledgments

This work was supported by the Brazilian funding agencies FACEPE (APQ-0732-1.06/14), BioMol/CAPES (BioComp 23038.004630/2014-35), CNPq (INCT-FCx 465259/2014-6). Computational resources were partially provided by the High Performance Computing Center North (HPC2N) and the Santos Dumont Supercomputer Center at the Brazilian National Laboratory of Scientific Computing (LNCC). T. A. S. acknowledges CNPq for a productivity fellowship.

References

• ¹
  Sadovskaya, I.; Brisson, J. R.; Lam, J. S.; Richards, J. C.; Altman, E.; Eur. J. Biochem. 1998, 255, 673.
• ²
  Sadovskaya, I.; Brisson, J. R.; Thibault, P.; Richards, J. C.; Lam, J. S.; Altman, E.; Eur. J. Biochem. 2000, 267, 1640.
• ³
  Wilkinson, S. G.; Prog. Lipid Res. 1996, 35, 283.
• ⁴
  Brandenburg, K.; Wiese, A.; Curr. Top. Med. Chem. 2004, 4, 1127.
• ⁵
  Erridge, C.; Bennett-Guerrero, E.; Poxton, I. R.; Microbes Infect. 2002, 4, 837.
• ⁶
  Brandenburg, K.; Koch, M. H.; Seydel, U.; J. Struct. Biol. 1990, 105, 11.
• ⁷
  Seydel, U.; Koch, M. H. J.; Brandenburg, K.; J. Struct. Biol. 1993, 110, 232.
• ⁸
  Coughlin, R. T.; Haug, A.; McGroarty, E. J.; Biochemistry 1983, 22, 2007.
• ⁹
  Brandenburg, K.; Seydel, U.; Biochim. Biophys. Acta 1984, 775, 225.
• ¹⁰
  Snyder, S.; Kim, D.; McIntosh, T. J.; Biochemistry 1999, 38, 10758.
• ¹¹
  Garidel, P.; Rappolt, M.; Schromm, A. B.; Howe, J.; Lohner, K.; Andra, J.; Koch, M. H. J.; Brandenburg, K.; Biochim. Biophys. Acta 2005, 1715, 122.
• ¹²
  Tong, J.; McIntosh, T. J.; Biophys. J. 2004, 86, 3759.
• ¹³
  Walsh, A. G.; Matewish, M. J.; Burrows, L. L.; Monteiro, M. A.; Perry, M. B.; Lam, J. S.; Mol. Microbiol. 2000, 35, 718.
• ¹⁴
  Nagle, J. F.; Tristram-Nagle, S.; Biochim. Biophys. Acta 2000, 1469, 159.
• ¹⁵
  Poger, D.; Caron, B.; Mark, A. E.; Biochim. Biophys. Acta 2016, 1858, 1556.
• ¹⁶
  Fogarty, J. C.; Arjunwadkar, M.; Pandit, S. A.; Pan, J.; Biochim. Biophys. Acta, Biomembr. 2015, 1848, 662.
• ¹⁷
  Marrink, S. J.; Corradi, V.; Souza, P. C. T.; Ingólfsson, H. I.; Tieleman, D. P.; Sansom, M. S. P.; Chem. Rev. 2019, 119, 6184.
• ¹⁸
  Hsia, C.-Y.; Chen, L.; Singh, R. R.; DeLisa, M. P.; Daniel, S.; Sci. Rep. 2016, 6, 32715.
• ¹⁹
  Kaufmann, S.; Ilg, K.; Mashaghi, A.; Textor, M.; Priem, B.; Aebi, M.; Reimhult, E.; Langmuir 2012, 28, 12199.
• ²⁰
  Clifton, L. A.; Skoda, M. W.; Daulton, E. L.; Hughes, A. V.; Le Brun, A. P.; Lakey, J. H.; Holt, S. A.; J. R. Soc., Interface 2013, 10, 20130810.
• ²¹
  Woods, R. J.; Dwek, R. A.; Edge, C. J.; Fraser-Reid, B.; J. Phys. Chem. 1995, 99, 3832.
• ²²
  Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A.; J. Am. Chem. Soc. 1995, 117, 5179.
• ²³
  Lins, R. D.; Straatsma, T. P.; Biophys. J. 2001, 81, 1037.
• ²⁴
  Lins, R. D.; Vorpagel, E. R.; Guglielmi, M.; Straatsma, T. P.; Biomacromolecules 2008, 9, 29.
• ²⁵
  Soares, T. A.; Straatsma, T. P.; Lins, R. D.; J. Braz. Chem. Soc. 2008, 19, 312.
• ²⁶
  Soares, T. A.; Straatsma, T. P.; AIP Conf. Proc. 2007, 963, 1375.
• ²⁷
  Straatsma, T. P.; Soares, T. A.; Proteins: Struct., Funct., Genet. 2009, 74, 475.
• ²⁸
  Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; Gonzalez-Outeirino, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J.; J. Comput. Chem. 2008, 29, 622.
• ²⁹
  Kirschner, K. N.; Lins, R. D.; Maass, A.; Soares, T. A.; J. Chem. Theory Comput. 2012, 8, 4719.
• ³⁰
  Dias, R. P.; da Hora, G. C. A.; Ramstedt, M.; Soares, T. A.; J. Chem. Theory Comput. 2014, 10, 2488.
• ³¹
  Dias, R. P.; Li, L.; Soares, T. A.; Alexov, E.; J. Comput. Chem. 2014, 35, 1418.
• ³²
  Nascimento, A.; Pontes, F. J.; Lins, R. D.; Soares, T. A.; Chem. Commun. 2014, 50, 231.
• ³³
  Ravi, H. K.; Stach, M.; Soares, T. A.; Darbre, T.; Reymond, J. L.; Cascella, M.; Chem. Commun. 2013, 49, 8821.
• ³⁴
  Sinha, S.; Zheng, L.; Mu, Y.; Ng, W. J.; Bhattacharjya, S.; Sci. Rep. 2017, 7, 17795.
• ³⁵
  Reinhardt, A.; Wehle, M.; Geissner, A.; Crouch, E. C.; Kang, Y.; Yang, Y.; Anish, C.; Santer, M.; Seeberger, P. H.; J. Struct. Biol. 2016, 195, 387.
• ³⁶
  Ma, H.; Irudayanathan, F. J.; Jiang, W.; Nangia, S.; J. Phys. Chem. B 2015, 119, 14668.
• ³⁷
  Van Oosten, B.; Harroun, T. A.; J. Mol. Graphics Modell. 2016, 63, 125.
• ³⁸
  Pontes, F. J. S.; Rusu, V. H.; Soares, T. A.; Lins, R. D.; J. Chem. Theory Comput. 2012, 8, 3830.
• ³⁹
  Santos, D. E. S.; Pol-Fachin, L.; Lins, R. D.; Soares, T. A.; J. Chem. Inf. Model. 2017, 57, 2181.
• ⁴⁰
  Piggot, T. J.; Holdbrook, D. A.; Khalid, S.; J. Phys. Chem. B 2011, 115, 13381.
• ⁴¹
  Wu, E. L.; Engstrom, O.; Jo, S.; Stuhlsatz, D.; Yeom, M. S.; Klauda, J. B.; Widmalm, G.; Im, W.; Biophys. J. 2013, 105, 1444.
• ⁴²
  Pol-Fachin, L.; Rusu, V. H.; Verli, H.; Lins, R. D.; J. Chem. Theory Comput. 2012, 8, 4681.
• ⁴³
  Oostenbrink, C.; Soares, T. A.; van der Vegt, N. F.; van Gunsteren, W. F.; Eur. Biophys. J. 2005, 34, 273.
• ⁴⁴
  Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F.; J. Comput. Chem. 2004, 25, 1656.
• ⁴⁵
  Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; van Gunsteren, W. F.; Eur. Biophys. J. 2011, 40, 843.
• ⁴⁶
  Lins, R. D.; Hunenberger, P. H.; J. Comput. Chem. 2005, 26, 1400.
• ⁴⁷
  Soares, T. A.; Hunenberger, P. H.; Kastenholz, M. A.; Krautler, V.; Lenz, T.; Lins, R. D.; Oostenbrink, C.; van Gunsteren, W. F.; J. Comput. Chem. 2005, 26, 725.
• ⁴⁸
  Chandrasekhar, I.; Kastenholz, M.; Lins, R. D.; Oostenbrink, C.; Schuler, L. D.; Tieleman, D. P.; van Gunsteren, W. F.; Eur. Biophys. J. 2003, 32, 67.
• ⁴⁹
  Lee, J.; Patel, D. S.; Stahle, J.; Park, S. J.; Kern, N. R.; Kim, S. H.; Lee, J.; Cheng, X.; Valvano, M. A.; Holst, O.; Knirel, Y. A.; Qi, Y.; Jo, S.; Klauda, J. B.; Widmalm, G.; Im, W.; J. Chem. Theory Comput. 2018, 15, 775.
• ⁵⁰
  Abraham, T.; Schooling, S. R.; Nieh, M. P.; Kucerka, N.; Beveridge, T. J.; Katsaras, J.; J. Phys. Chem. B 2007, 111, 2477.
• ⁵¹
  Kucerka, N.; Papp-Szabo, E.; Nieh, M. P.; Harroun, T. A.; Schooling, S. R.; Pencer, J.; Nicholson, E. A.; Beveridge, T. J.; Katsaras, J.; J. Phys. Chem. B 2008, 112, 8057.
• ⁵²
  Poger, D.; Mark, A. E.; J. Chem. Theory Comput. 2012, 8, 4807.
• ⁵³
  Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E.; J. Chem. Theory Comput. 2008, 4, 435.
• ⁵⁴
  van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C.; J. Comput. Chem. 2005, 26, 1701.
• ⁵⁵
  Reißer, S.; Poger, D.; Stroet, M.; Mark, A. E.; J. Chem. Theory Comput. 2017, 13, 2367.
• ⁵⁶
  Silva, T. F. D.; Vila-Viçosa, D.; Reis, P. B. P. S.; Victor, B. L.; Diem, M.; Oostenbrink, C.; Machuqueiro, M.; J. Chem. Theory Comput. 2018, 14, 5823.
• ⁵⁷
  Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E.; Bioinformatics 2013, 29, 845.
• ⁵⁸
  Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., ed.; D. Reidel: Dordrecht, 1981, p. 331.
• ⁵⁹
  Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P.; J. Phys. Chem. 1987, 91, 6269.
• ⁶⁰
  Soares, T. A.; Daura, X.; Oostenbrink, C.; Smith, L. J.; van Gunsteren, W. F.; J. Biomol. NMR 2004, 30, 407.
• ⁶¹
  Reif, M. M.; Winger, M.; Oostenbrink, C.; J. Chem. Theory Comput. 2013, 9, 1247.
• ⁶²
  Poger, D.; Mark, A. E.; J. Chem. Theory Comput. 2010, 6, 325.
• ⁶³
  Pol-Fachin, L.; Verli, H.; Lins, R. D.; J. Comput. Chem. 2014, 35, 2087.
• ⁶⁴
  Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A.; J. Phys. Chem. 1993, 97, 10269.
• ⁶⁵
  Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F.; J. Chem. Phys. 1995, 102, 5451.
• ⁶⁶
  Darden, T.; York, D.; Pedersen, L.; J. Chem. Phys. 1993, 98, 10089.
• ⁶⁷
  Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M.; J. Comput. Chem. 1997, 18, 1463.
• ⁶⁸
  Miyamoto, S.; Kollman, P. A.; J. Comput. Chem. 1992, 13, 952.
• ⁶⁹
  Bussi, G.; Donadio, D.; Parrinello, M.; J. Chem. Phys. 2007, 126, 014101.
• ⁷⁰
  Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R.; J. Chem. Phys. 1984, 81, 3684.
• ⁷¹
  Essex, J. W.; Mol. Simul. 1998, 20, 159.
• ⁷²
  Glattli, A.; Daura, X.; van Gunsteren, W. F.; J. Chem. Phys. 2002, 116, 9811.
• ⁷³
  Brandenburg, K.; Funari, S. S.; Koch, M. H.; Seydel, U.; J. Struct. Biol. 1999, 128, 175.
• ⁷⁴
  Le Brun, A. P.; Clifton, L. A.; Halbert, C. E.; Lin, B.; Meron, M.; Holden, P. J.; Lakey, J. H.; Holt, S. A.; Biomacromolecules 2013, 14, 2014.
• ⁷⁵
  van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C.; J. Chem. Phys. 1998, 108, 10220.
• ⁷⁶
  Zielkiewicz, J.; J. Chem. Phys. 2005, 123, 104501.
• ⁷⁷
  Brandenburg, K.; Biophys. J. 1993, 64, 1215.
• ⁷⁸
  Brandenburg, K.; Seydel, U.; Eur. Biophys. J. 1988, 16, 83.
• ⁷⁹
  Brandenburg, K.; Seydel, U.; Eur. J. Biochem. 1990, 191, 229.
• ⁸⁰
  Brandenburg, K.; Andra, J.; Muller, M.; Koch, M. H.; Garidel, P.; Carbohydr. Res. 2003, 338, 2477.
• ⁸¹
  Reif, M. M.; Oostenbrink, C.; Theor. Chem. Acc. 2015, 134, 2.
• ⁸²
  Poger, D.; Van Gunsteren, W. F.; Mark, A. E.; J. Comput. Chem. 2010, 31, 1117.

Publication Dates

• Publication in this collection
  21 Oct 2019
• Date of issue
  Oct 2019

History

• Received
  5 Feb 2019
• Accepted
  29 May 2019