Ganglioside-Lipid and Ganglioside-Protein Interactions Revealed by Coarse-Grained and Atomistic Molecular Dynamics Simulations.

Gangliosides are glycolipids in which an oligosaccharide headgroup containing one or more sialic acids is connected to a ceramide. Gangliosides reside in the outer leaflet of the plasma membrane and play a crucial role in various physiological processes such as cell signal transduction and neuronal differentiation by modulating structures and functions of membrane proteins. Because the detailed behavior of gangliosides and protein-ganglioside interactions are poorly known, we investigated the interactions between the gangliosides GM1 and GM3 and the proteins aquaporin (AQP1) and WALP23 using equilibrium molecular dynamics simulations and potential of mean force calculations at both coarse-grained (CG) and atomistic levels. In atomistic simulations, on the basis of the GROMOS force field, ganglioside aggregation appears to be a result of the balance between hydrogen bond interactions and steric hindrance of the headgroups. GM3 clusters are slightly larger and more ordered than GM1 clusters due to the smaller headgroup of GM3. The different structures of GM1 and GM3 clusters from atomistic simulations are not observed at the CG level based on the Martini model, implying a difference in driving forces for ganglioside interactions in atomistic and CG simulations. For protein-ganglioside interactions, in the atomistic simulations, GM1 lipids bind to specific sites on the AQP1 surface, whereas they are depleted from WALP23. In the CG simulations, the ganglioside binding sites on the AQP1 surface are similar, but ganglioside aggregation and protein-ganglioside interactions are more prevalent than in the atomistic simulations. Using the polarizable Martini water model, results were closer to the atomistic simulations. Although experimental data for validation is lacking, we proposed modified Martini parameters for gangliosides to more closely mimic the sizes and structures of ganglioside clusters observed at the atomistic level.

Introduction

Glycolipids are amphiphilic molecules composed of a pair of hydrophobic alkyl tails, anchoring the lipid to the bilayer, and a hydrophilic oligosaccharide headgroup. The headgroup resides on the bilayer surface and interacts with aqueous and plasma membrane constituents from the same cell or neighboring cells.[1−5] Various species of glycolipids exist with different oligosaccharide headgroups, lipid backbones, and alkyl tails. In glycoglycerolipids, the headgroups are attached to glycerol, and in glycosphingolipids (GSLs), they are attached to ceramide. Gangliosides are GSLs with one or more sialic acids, hence one or more negative charges, in their oligosaccharide headgroups. In GSLs, the headgroups are connected to ceramide backbones via the 1-hydroxyl moiety of one of the sugars.[6] Gangliosides are involved in many physiological processes,[3,7,8] such as cell adhesion,[9] pathogen recognition and viral infection,[10,11] signal transduction,[12] cell proliferation, and neuronal protection.[8,13] The diversity of their biological functions is related to the diversity of their oligosaccharide headgroups.[6] Two common GSLs are GM1 and GM3 (Figure ), both with one sialic acid but different headgroup sizes. GM1 is enriched in the gray matter of the central nervous system and plays a key role in neurotrophic and neuroprotective processes[8,14] as well as neurodegenerative diseases.[15,16] GM3 is thought to regulate cell growth, adhesion, and motility.[17] Abnormal lipid microdomains involving GM3 are associated with cancer development[17] and insulin resistance in type 2 diabetes.[18−20]

Figure 1

Structures of GM1, GM3, and POPC.

In bilayer environments, gangliosides segregate into microdomains (rafts) enriched in sphingomyelin, cholesterol, and gangliosides,[21] and glycosynapses, which are formed by specific proteins and gangliosides but without cholesterol.[22−24] These microdomains are highly dynamic, both spatially and temporally. Gangliosides can influence the lateral distributions of lipids in membranes.[25] They can modulate the morphology of bilayers by inducing membrane curvature and are involved in membrane fusion and viral budding.[26−29] Gangliosides specifically interact with and regulate the functions of various proteins, e.g., the epidermal growth factor receptor (EGFR)[30] and insulin receptor.[31] They also influence the lateral sorting of transmembrane lipid raft proteins.[32]

An important structural feature of gangliosides is their hydrogen bonding ability through their headgroups, which results in extensive ganglioside and protein-ganglioside interactions. However, the negative charges, steric hindrance, and hydration shell around the headgroups modulate these interactions. Other characteristics that contribute to ganglioside and protein-ganglioside interactions include the higher transition temperature due to long saturated tails, possible hydrogen bonding between the amide groups of the ceramide tails, and interdigitation between leaflets due to their long tails.[24] Differences in ganglioside structures may result in significant differences between their interactions. GM1 and GM3 differ by two sugar units (GalNAc and Gal2, shown in Figure ), but GM1 and GM3 clusters are separated from each other in membranes.[33,34]

Ganglioside and protein-ganglioside interactions have been investigated in model lipid mixtures both experimentally and with molecular dynamics (MD) simulations. The experiments have focused on detecting lipid domains containing gangliosides, as reviewed in refs (27 and 35), for example. The primary advantage of MD simulations is that they can resolve detailed lipid interactions.[36−38] Manna et al.[39] recently reviewed glycolipid MD simulations, highlighting the structural characterizations that affect lipid organizations in various environments. Hall et al. analyzed interactions in lipid mixtures of galactosylceramide (GalCer), a glycosphingolipid, and other lipid raft components.[40,41] Sega et al. explored properties of GM3 bilayers by MD simulation.[42] Patel et al.[43,44] and Jedlovszky et al.[45] investigated the conformations and distributions of GM1 in 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) bilayers. Mori et al. compared GM1 clusters in 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayers and lipid mixtures containing sphingomyelin and cholesterol.[46]

A notable part of the simulation studies has been conducted at the coarse-grained (CG) level using the Martini force field.[47−49] In Martini, four heavy atoms are generally combined to form one CG bead. The force field is parametrized based on free energy landscapes of molecular interactions, particularly the partitioning free energies of molecules between oil- and water-like phases. Martini enables simulations of larger systems and/or longer time scales due to the reduced level of detail while still maintaining a reasonable degree of chemical accuracy. For instance, the plasma membrane simulations by Ingólfsson et al. found dynamic ganglioside clusters,[50] whereas another simulation of complex lipid mixtures revealed correlation between bilayer curvature and GM3 clusters.[51] Shorthouse et al.[52] reviewed a coarse-grained model of GM3 and its application in protein-ganglioside interactions. Prasanna et al.[53] investigated binding of GM1 to the serotonin 1A receptor. Jong et al. simulated sorting of proteins to liquid ordered (Lo) domains induced by gangliosides and lipid anchors.[54] Basu et al.[55] and Kociurzynski et al.[56] simulated phase separation and transitions of lipid mixtures containing gangliosides.

Despite extensive studies of ganglioside-containing lipid mixtures, the strength of ganglioside and protein-ganglioside interactions is still poorly known (e.g, the size and stability of ganglioside clusters). In this study, we conducted both equilibrium MD simulations and free energy calculations of gangliosides, or gangliosides and proteins, in POPC bilayers both at the atomistic and CG levels to explore their interactions. We simulated GM1 and GM3 (Figure ) as examples of gangliosides with different headgroup sizes and WALP23 and aquaporin (AQP1) as examples of membrane proteins. Our simulations provide detailed insight on the strength of ganglioside and protein-ganglioside interactions in a lipid bilayer environment and reveal how different types of ganglioside headgroups affect ganglioside aggregation.

Methods

System Setup

Both equilibrium simulations and PMF calculations were performed to investigate ganglioside and protein-ganglioside interactions. All simulations are summarized in Table .

Table 1

Summary of Simulations

	system composition	simulation length (μs)	simulation levela	water modelb	electrostatic interactionsc
Equilibrium Simulations
GM1-POPC	upper leaflet: 24 GM1, 120 POPC	6.5	CG	W	shift
	lower leaflet: 144 POPC	6.5	CG	PW	shift
		6.5	CG	PW	PME
		3	AA	SPC	RF
GM3-POPC	upper leaflet: 24 GM3, 120 POPC	6.5	CG	W	shift
	lower leaflet: 144 POPC	6.5	CG	PW	shift
		6.5	CG	PW	PME
		3	AA	SPC	RF
GM1-GM3-POPC	upper leaflet: 12 GM1, 12 GM3, 120 POPC	2	AA	SPC	RF
	lower leaflet: 144 POPC
WALP23-GM1- POPC	protein: 1 WALP23	15	CG	W	shift
	upper leaflet: 23 GM1, 116 POPC	15	CG	PW	shift
	lower leaflet: 139 POPC	15	CG	PW	PME
		2	AA	SPC	RF
AQP1-GM1-POPC	protein: 1 AQP1	15	CG	W	shift
	upper leaflet: 41 GM1, 207 POPC	15	CG	PW	shift
	lower leaflet: 248 POPC	15	CG	PW	PME
		2	AA	SPC	RF

CG: coarse-grained simulation; AA: atomistic simulation.

W: standard Martini water model; PW: polarizable Martini water model; SPC: simple point charge atomistic water model.

PME: particle mesh Ewald; RF: reaction field.

For equilibrium simulations, GM-POPC binary mixtures containing 17% GM1 or GM3 in the upper leaflet were simulated to investigate ganglioside interactions. We also simulated a GM1-GM3-POPC ternary mixture with 8.5% GM1 and 8.5% GM3 in the upper leaflet to illustrate the interactions between GM1 and GM3. Simulations of WALP23 (sequence: GWW(LA)8LWWA) and aquaporin 1 (AQP1) were conducted in binary mixtures with 17% GM1 in the upper leaflet to investigate protein-ganglioside interactions. Equilibrium simulations were conducted at both atomistic and CG levels. Free energy calculations were only conducted at the CG level. PMFs of the interaction between two GM1 lipids, two GM3 lipids, WALP23 and GM3, and AQP1 and GM3 in POPC bilayers were calculated as a function of distances in the x–y plane. For both the equilibrium simulations and PMF calculations, three parallel simulations were conducted at the CG level with different water models and parameters for electrostatic interactions (see below for detailed parameters).

The CG systems were built using the insane script.[57] The CG model of GM lipid tail has the bead order AM1-AM2-T1A-C2A-C3A-C1B-C2B-C3B-C4B. The CG POPC lipid has the bead order NC3-PO4-GL1-GL2-C1A-C2A-D3A-C4A-C1B-C2B-C3B-C4B. WALP23 was built as an ideal helix and then equilibrated for 100 ns in a POPC bilayer before being used in the ganglioside simulations, whereas the AQP1 simulations were started from the crystal structure (PDB ID: 1J4N).[58] Atomistic simulation systems were backmapped from equilibrated snapshots of the corresponding CG simulations.[59] GM(d20:1/20:0) lipids were used in the atomistic simulations (i.e., the tail was constituted by a long-chain base of 20 carbons and an eicosyl fatty acid residue, as shown in Figure ).

CG Simulation Details

CG simulations were performed with the Martini force field[49] and the ganglioside lipid extension.[48,50] Simulations were done with the GROMACS 4.6.5 package[60,61] using an integration time step of 20 fs. Then, 0.15 M NaCl was added to mimic the physiological ionic strength. Temperature was maintained at 310 K using the v-rescale algorithm[62] with a relaxation time of 2.0 ps. Semi-isotropic pressure coupling was applied with a reference pressure of 1 bar using the Parrinello–Rahman algorithm[63,64] with a relaxation time of 12.0 ps; van der Waals interactions were turned off between 0.9 and 1.2 nm using the shift function in GROMACS. Three parallel simulations were conducted for each system with different Martini water models and parameters for electrostatic interactions. In the first set, called CG W, the standard Martini water[49,65] was employed. The shift function in GROMACS was used to turn off the electrostatic interactions between 0 and 1.2 nm, and a dielectric constant of 15 was used. In the second set, called CG PW, polarizable Martini water was used with the same cutoff distance and a dielectric constant of 2.5.[66] In the third and last set, called CG PW PME, polarizable Martini water and a dielectric constant of 2.5 were also used but additional long-range electrostatic interactions were included using the PME method.[67,68] ElNeDyn approach was used to restrain the protein structures.[69] All simulations and conditions are listed in Table .

Atomistic Simulation Details

Atomistic simulations were conducted with the GROMACS 4.0.7 package[60,61] using an integration time step of 2 fs. Parameters for POPC and the lipid part of gangliosides are based on the GROMOS 53A6 united-atom force field,[70] whereas parameters of ganglioside headgroups are based on the GROMOS hexopyranose force field.[48,71,72] The systems were solvated in SPC water, and both counterions and 0.15 M NaCl were added. Temperature was maintained at 330 K by coupling the bilayer/protein and solution separately via the Berendsen thermostat algorithm[73] with a relaxation time of 0.1 ps. A temperature of 330 K was used in atomistic simulations because we observed gel phase for gangliosides at 310 K. A 330 K simulation of the GM3-POPC mixture at the CG level was conducted as a control (see Discussion for details). Semi-isotropic pressure coupling with a relaxation time of 0.5 ps was applied to maintain the pressure at 1 bar.[73] Nonbonded interactions within 0.9 nm were calculated every step based on a pair list that was updated every five steps. Interactions between 0.9 and 1.4 nm were calculated every five steps.[74] Long-range electrostatic interactions were included using the reaction-field algorithm with a relative dielectric permittivity of 54.[75] In most of the simulations, the proteins and lipids were free to move, but in the AQP1 simulations, position restraints of 10 kJ mol–1 nm–2 were applied on the protein backbone to maintain a backbone RMSD of ∼0.22 nm.

PMF Calculations

Umbrella sampling[76] was employed to calculate the free energies of ganglioside and protein-ganglioside interactions at the CG level in POPC bilayers (Table ). Starting structures for each window, spaced 1.0 Å, were taken from an initial 0.5 μs simulation in which one ganglioside molecule was pulled away from a second ganglioside or the protein in the plane of the bilayer. A force constant of 1000 kJ mol–1 nm–2 was used to control the distances between the centers of mass of gangliosides (or distances between ganglioside and protein) in the x–y plane. Twenty-three to thirty-six windows were used. MD simulations for each window were performed for 1.5–5.5 μs depending on the simulation system (see Table ). The first 0.5 μs of each simulation was taken as equilibration, and the remaining parts were used to generate the PMF profiles based on the weighted histogram analysis method (WHAM).[77] For the convergence to be verified, the trajectories were divided into blocks, and PMF profiles based on each block were calculated and compared. For GM3-AQP1 PMF, extensive sampling of GM3 around AQP1 is computationally too costly due to the large size of the protein. Therefore, PLUMED[78] was used to control the orientation between GM3 and AQP1, so that GM3 approached AQP1 in one specific direction (see Figure S1) and a one-dimensional PMF profile as a function of the distance between them was obtained. GM3 was used to calculate the protein-ganglioside interactions instead of GM1 because the GM3 lipid has a smaller headgroup and the convergence of the PMF profile is faster. Given the simulation lengths required for the GM PMFs, accurate atomistic PMFs are currently not feasible with this protocol.

Modifications to the Martini Ganglioside Force Field

Because the aggregation of gangliosides in the CG simulations was significantly stronger than in the atomistic simulations (see section of the results for details), we reconsidered the CG Martini GM1 and GM3 parameters. The structure and size of ganglioside clusters, and the radial distribution functions of gangliosides in the POPC bilayer at the atomistic level, were used as criteria. The atomistic to CG mapping (described in López et al.[48]) was kept unchanged, the bonded and nonbonded interactions were modified slightly (see section of the results and the Supporting Information for detailed changes). For the bonded interactions, atomistic trajectories were converted to the CG level every 0.5 ns, and the distributions of the bonds, angles, and dihedrals were analyzed and used as references. Several new angles and dihedrals were introduced to better control ganglioside conformations at the CG level. For nonbonded interactions, different types of CG beads were assigned to the headgroups to reduce their interactions. Equilibrium simulations of GM1-POPC and GM3-POPC binary mixtures were conducted using the modified CG force field and compared with the previous parameters.

Data Analysis

The initial part of each simulation was excluded from the analyses. For the CG simulations of 6.5, 10, and 15 μs in length, the last 4, 7, and 12 μs were used for analysis, respectively, whereas for the atomistic simulations of 2 and 3 μs in length, the last 1 and 1.5 μs were used, respectively. The analysis of ganglioside cluster sizes was performed using the g_aggregate tool.[79] Gangliosides were considered to be in the same cluster if the distance between their centers of mass in the x–y plane was within 1.5 nm. Radial distribution functions (RDFs) were calculated based on the distances between the centers of mass of lipids/ganglioside headgroups/tails and proteins in the x–y plane. RDF profiles were only calculated for the lipids in the upper leaflet. Hydrogen bonds were defined in the following way: the distance between the two heavy atoms is smaller than 3.5 Å, and the angle between donor, hydrogen atom, and acceptor is smaller than 150°. Order parameters of the lipid tails at the CG level were calculated based on the angle between the bonds of the lipid tails and the bilayer normal (approximated as the box z axis) using the g_ordercg tool.[80]

Results

First, the ganglioside-ganglioside interactions in simulations of mixed ganglioside-POPC bilayers using the original Martini force field and the GROMOS force field are compared. Sizes and structures of ganglioside clusters in the CG W, CG PW, CG PW PME, and atomistic simulations are described in detail. Next, proteins are included and the protein-ganglioside interactions at the CG and atomistic levels are examined. The distributions of GM lipids around WALP23 and AQP1 are presented. Lastly, modifications of the Martini ganglioside force field are conducted to mimic the ganglioside clustering in the atomistic POPC bilayer.

Interactions between Gangliosides

First, we analyze the ganglioside-POPC mixtures (Table ). The bilayer properties (Table S1) and the density profiles of lipids (Figure S2) are calculated for the atomistic and three CG simulations. The polarizable Martini (PW) water and long-range electrostatic interactions only slightly affect bilayer properties. We focus here on ganglioside-ganglioside interactions and the formation of ganglioside clusters.

Ganglioside Distribution in POPC Bilayers

The distribution of GM1 and GM3 is explored in POPC lipid simulations with 17% GM lipids in the upper leaflet (Figure ). Under all of the tested conditions, the gangliosides cluster, but the degree of aggregation varies, as shown by the snapshots at the end of each simulation (Figure A). The averaged ganglioside cluster size and the number of clusters are calculated as a function of simulation time (Figure B and Table ). GM1 aggregates most in the CG W simulations. Using PW decreases the cluster size significantly, whereas also including long-range electrostatic interactions (PME) has almost no additional effect. The average GM1 cluster sizes in the three CG simulations are 18, 5, and 4 lipids for the W, PW, and PW PME systems, respectively (Table ). In the atomistic simulation, a large GM1 cluster exists at the beginning, as the initial conformation is backmapped from an equilibrated CG W snapshot. However, these large clusters quickly disassociate. After equilibration, more than 10 clusters are found over most of the simulation time with an average cluster size of only 2 lipids, indicating weaker GM1 interactions at the atomistic level (Figure B and Table ). Simulations of GM3 reveal similar results, i.e., including PW weakens ganglioside interactions dramatically, and GM3 interactions at the atomistic level are the weakest among all of the four simulations. The GM3 cluster sizes in the three CG simulations (W, PW, and PW PME) and the atomistic simulation are 6, 3, 3, and 3 lipids, respectively (Table and Figure B). GM3 clusters are somewhat smaller than GM1 at the CG level, whereas the opposite is observed at the atomistic level.

Figure 2

Ganglioside interactions in the GM-POPC binary mixtures. (A) Snapshot of the upper leaflet at the end of the simulations. GM lipids and POPC are shown in red and gray, respectively. (B) Averaged sizes and numbers of ganglioside clusters as a function of simulation time. (C) PMF calculations of ganglioside interactions at the CG level. Abbreviations: CG, coarse-grained simulations; W, standard Martini water; PW, polarizable Martini water; PME, PME algorithm for long-range electrostatic interactions; AA, atomistic simulations.

Table 2

Average Size and Number of Ganglioside Clusters in POPC Bilayersa

lipid type		CG W	CG PW	CG PW PME	AA	CG W New FF	CG PW New FF
GM1	size	18 ± 7	5 ± 2	4 ± 1	2 ± 1	3 ± 1	2 ± 1
	number	2 ± 1	6 ± 2	6 ± 2	11 ± 2	10 ± 2	11 ± 2
GM3	size	6 ± 2	3 ± 1	3 ± 1	3 ± 1	4 ± 1	3 ± 1
	number	5 ± 1	8 ± 2	8 ± 2	10 ± 2	7 ± 2	10 ± 2

Abbreviations: CG W, CG simulations with standard Martini water; CG PW, CG simulations with polarizable water; CG PW PME, CG simulations with polarizable water in which PME was used for long-range electrostatic interactions; AA, atomistic simulations; New FF, simulations with reoptimized Martini ganglioside force field. Values are shown as average ± standard deviation.

The strength of ganglioside interactions revealed by three equilibrium simulations at the CG level is consistent with the PMF calculations: the association energy of GM1-GM1 is significantly larger than that of GM3-GM3 (in CG W simulations, the values are approximately −12 and −5.5 kJ/mol, respectively; Figure C). The calculated energies with the PW are roughly half of the values obtained with the standard Martini water (approximately −6.5 and −2.5 kJ/mol for GM1-GM1 and GM3-GM3, respectively).

Structures of Ganglioside Clusters

In addition to the cluster analyses, we calculated radial distribution function (RDF) profiles (Figure A and B; examples of convergence are shown in Figure S3) to investigate the structures of ganglioside clusters. RDFs for all CG simulations show peaks at ∼0.8 nm, but the magnitudes are decreased significantly by PW for both GM1 and GM3, which is consistent with the cluster sizes. However, in the atomistic simulations, the peak of the GM1-GM1 RDF is shifted to ∼1.6 nm and is much lower. For GM3, there is a strong peak at ∼0.6 nm and a second weak peak around 2.0 nm. Although the first peak of the GM3-GM3 RDF shows a magnitude comparable with that of CG PW, its width is much narrower, implying clusters with tighter packed structures.

Figure 3

Structures of ganglioside clusters. (A) RDFs between the centers of mass of GM lipids and RDFs between the centers of mass of GM and POPC. (B) RDFs between the GM headgroups and between the GM tails. (C) GM1 and (D) GM3 clusters in the CG PW PME simulations. (E) GM1 and (F) GM3 clusters in the atomistic simulations. Gangliosides and POPC are shown in yellow and cyan, respectively. The ganglioside headgroups at the CG level are shown in pink. The sialic acid groups of gangliosides are shown in green. Hydrogen bonds are shown with dotted lines. Abbreviations: CG, coarse-grained simulations; W, standard Martini water; PW, polarizable Martini water; PME, PME algorithm for long-range electrostatic interactions; AA, atomistic simulations.

The RDFs of the headgroups and tails are plotted separately to further illustrate the ganglioside interactions (Figure B). Peaks are found for both headgroups and tails in all of the CG simulations. The headgroup peaks of GM1 are shifted to the right compared with GM3, which is consistent with the larger size of the GM1 headgroup. In the atomistic simulations of GM3, peaks are found for the RDFs of both the headgroups and tails at ∼0.6 nm, although the peak for tails is much weaker. However, for GM1, only the headgroup RDF shows a weak peak at ∼1.1 nm. Thus, at the atomistic level, both headgroups and tails interact with each other in GM3 clusters but only the headgroups aggregate in GM1 clusters. Minor differences in ganglioside cluster sizes and significant differences in the shapes of RDFs indicate different structures of ganglioside clusters.

The GM-POPC RDFs show peaks around 0.9 nm in the atomistic simulations, but the peaks are much weaker in the CG simulations (Figure A). In the atomistic simulations, these peaks are wider for GM1 than GM3 due to the larger spaces between GM1 lipid tails that may allow accommodation of POPC lipids.

Typical GM1 and GM3 clusters at both the CG and atomistic levels are shown in Figure C–F. Clusters in atomistic simulations consisting of two or three (Figure E and F) gangliosides are shown together with associated POPC molecules. In the GM1 clusters, Neu5Ac, and Gal2, GalNAc extend in different directions and form hydrogen bonds, whereas Gal, Glc, and ceramide groups are further away from each other (Figure E), consistent with the fact that peaks are only found for RDFs of headgroups but not for the RDFs of tails (Figure B). However, in the GM3 clusters, the headgroups are almost planar and tightly packed parallel to each other (Figure F). In this way, the ceramide groups are able to approach each other close enough to form hydrogen bonds (Figure F). The POPC lipid tails are inserted into the ganglioside clusters, whereas the phosphate groups form hydrogen bonds with the Glc and ceramide of the gangliosides, as shown in Figure E and F. The clusters in the CG simulations do not show distinctly different structures between GM1 and GM3 (Figure C and D), and both the headgroups and tails have extensive contacts with each other, which are consistent with the peaks in the RDFs in Figure B.

Hydrogen Bond Interactions Involving Gangliosides

Because the ganglioside headgroups can act as both hydrogen bond donors and acceptors, we analyzed intermolecular hydrogen bond interactions in the atomistic simulations; the results are shown in Table . Neu5Ac (the sialic acid) contributes most to the hydrogen bond interactions between gangliosides for both GM1 and GM3, which are consistent with Mori et al.’s simulations.[46] In all of the atomistic simulations, hydrogen bonds involving Neu5Ac account for at least 55% of the total number of hydrogen bonds between the gangliosides. For GM3, these hydrogen bonds are mainly formed between Neu5Ac of different molecules, whereas for GM1, Gal2 also forms hydrogen bonds with Neu5Ac in addition to the hydrogen bonds between Neu5Ac groups. Although GM1 has a larger headgroup and more hydrogen bond donors and acceptors, the averaged intermolecular hydrogen bond numbers are smaller than for GM3.

Table 3

Hydrogen Bond Interactions between Gangliosides in Atomistic Simulations

hydrogen bonds between gangliosides
simulation system	GM	Neu5Ac	Neu5Ac	Neu5Ac	Neu5Ac	Gal2	Gal2	Gal1
	GMa	Neu5Ac	Gal2	GalNAc	Gal1	Gal2	GalNAc	Glc
GM1-POPC	1.97 ± 0.73	27.5%	26.4%	7.1%	6.1%	8.1%	5.1%	<2%
GM3-POPC	2.73 ± 1.62	65.5%			30.2%			3.4%

Average numbers of hydrogen bonds between one molecule in the first group and all of the molecules in the second group and breakdown of the H-bonds in percentages due to different fragments. If two groups are the same, the hydrogen bond number between one molecule and all of the other molecules in the group is calculated. The nomenclature for the sugars is shown in Figure .

Phos: Phosphate group.

Cer: Ceramide.

Hydrogen bonds are also found between the linkage parts of gangliosides, in particular for GM3, which packed tightly to form clusters. The average number of hydrogen bonds between amide groups in the GM3 simulations is ∼1.6, which means there are one or two pairs of GM3 lipids forming this kind of hydrogen bond in each simulation frame (for GM1, this is only 0.4). These are notably common interactions when considering the ratio of ganglioside to POPC lipids in this system.

Extensive hydrogen bond interactions are also found between gangliosides and POPC lipids. Each ganglioside forms ∼2.5 hydrogen bonds with POPC lipids, and most of them are hydrogen bonds between the phosphate groups of the POPC and the Glc and ceramide backbones of gangliosides (Table ). Hydrogen bonds between phosphate and Glc, and hydrogen bonds between phosphate and ceramide, each account for 28–35% of the total number of hydrogen bonds. The results are similar for GM1 and GM3.

Interactions between GM1 and GM3

For the interactions between GM1 and GM3 to be investigated, an atomistic MD simulation of a GM1-GM3-POPC ternary mixture was performed. Both GM1 and GM3 clusters are observed (Figure ), as in the simulations of binary mixtures, as well as small mixed clusters of GM1 and GM3. Peaks are found for the GM3-GM3 RDFs for entire lipids, headgroups, and tails (Figure B–D) with comparable positions and magnitudes as in the binary mixtures (Figure A and B). For GM1-GM1, no clear peak is observed for the entire lipid RDF, but a weak peak for the headgroup RDF is found at ∼1.2 nm (Figure C) as in the binary mixtures (Figure B). These results do not reveal notable differences in structure and size of ganglioside clusters in binary and ternary mixtures for either GM1 or GM3. No strong interactions between GM1 and GM3 were found. None of the GM1-GM3 RDFs for lipids, headgroups, and tails show defined peaks (Figure B–D), suggesting limited GM1-GM3 clustering. GM1/GM3-POPC RDFs are shown in Figure E.

Figure 4

Ganglioside distributions in the atomistic GM1-GM3-POPC ternary mixture. (A) Snapshots of the simulation system. POPC lipids are shown in gray, and GM1 and GM3 are in red and blue, respectively. RDFs between the centers of mass of (B) GM lipids, (C) GM headgroups, (D) GM tails, and (E) GM and POPC lipids are shown.

We calculated the number of contacts between the gangliosides, averaged every 0.5 μs (Figure S4 and Table S2). The initial structure of the ternary mixture was created based on the last snapshot of the GM1-POPC simulation by changing some of the GM1 headgroups to GM3 headgroups; during the 2 μs simulation, the number of contacts between GM1 and GM3 decreased, whereas the GM3-GM3 contacts increased, and the number of GM1-GM1 contacts fluctuated. Our results are consistent with experiments conducted by Fujita et al.,[33,34] which revealed separated GM1 and GM3 clusters, supporting the reliability of our atomistic simulations.

Protein-Ganglioside Interactions

Interactions between WALP23 and GM1

2D density maps (Figure A) and RDF profiles (Figure C) revealed that WALP23-GM1 interactions at the CG level are stronger than at the atomistic level. Using PW reduces the interactions dramatically, similarly to what we found in the lipid mixtures. Using PW PME makes the WALP23-GM1 interactions stronger than in the PW simulation without PME. In the simulation with PW (without PME), the RDF does not show peaks and is comparable to the corresponding RDF of the atomistic simulation (Figure C). In these two simulations (PW and AA), GM1 is depleted from WALP23. The GM1-protein RDF of the PW PME simulation shows a peak at ∼1.4 nm (Figure C). The 2D density maps are consistent with RDF profiles, as shown by the positions of high GM1 densities (Figure A). RDFs of the headgroups and tails in Figure D reveal that GM1 tail-WALP23 interactions are weaker than the GM1 headgroup-WALP23 interactions in all of the CG and atomistic simulations. The strength of GM-WALP23 interactions suggested by free energy calculations are consistent with the equilibrium simulations. The association energies of GM3-WALP23 are −2.5, −0.7, and −0.7 kJ/mol in the W, PW, and PW PME CG simulations, respectively (Figure B).

Figure 5

Ganglioside-protein interactions in the GM1-POPC-protein simulations. (A) 2D density maps of GM1 lipids in the bilayer plane. The relative densities are shown. Structures of WALP23 and AQP1 are included to show their orientations during the density calculations. (B) PMF calculations of protein-ganglioside interactions at the CG level. (C) RDFs of GM1 and POPC lipids relative to the mass centers of the proteins. (D) RDFs of GM1 headgroups and tails relative to the mass centers of the proteins.

In the atomistic simulations, POPC is enriched around WALP23 compared with GM1 (see RDFs of GM1 and POPC lipids in Figure C). Figure A shows the snapshot at the end of the atomistic simulation. Most of the GM1 lipids are away from WALP23 except for one or two GM1 lipids (Figure A). The tails of that GM1 lipid do not come close to the peptide, but its headgroup tilts and interacts with the positively charged N-terminus via Neu5Ac (the sialic acid ,which has one negative charge) and GalNAc (Figure C).

Figure 6

Interactions between GM1 lipids and WALP23 or AQP1 in atomistic simulations of GM1-POPC-protein systems. Snapshots of the simulation systems of (A) WALP23 and (B) AQP1 are shown. WALP23 is shown in cyan, and subunits 1–4 of AQP1 are shown in red, light blue, orange, and light green, respectively. POPC lipids and GM1 lipids are shown in gray and red, respectively. Additionally, in panel B, the two GM1 lipids, which are stably bound to AQP1 during the last 1 μs simulation, are in blue and the 7 GM1 lipids, which have interactions with AQP1 or the 2 stably bound GM1 lipids, are shown in magenta. Atomistic details of the interactions between GM1 and (C) WALP23 and (D) subunit 1 of AQP1. (E) Side view and (F) top view of the binding site of GM1 on subunit 1 of AQP1. Panels C–F are based on the last snapshots of the simulations. GM1 and protein residues are shown in yellow and cyan, and the backbones of WALP23 and AQP1 are shown in cyan and white, respectively.

Interactions between AQP1 and GM1

Gangliosides interact with AQP1 directly and indirectly by forming clusters with other gangliosides that contact AQP1 (see Figure B). Our CG simulations suggest that PW does not affect the strength of direct GM1-AQP1 interactions, as indicated by the high GM1 densities with comparable amplitudes at similar positions around the AQP1 in the CG W and PW simulations (Figure A). However, in the CG PW simulations, the regions with high relative GM1 densities (regions in blue and white) are smaller compared with those in the CG W simulations, which means the indirect GM1-AQP1 interactions are decreased because the GM1 clusters become smaller and less stable.

The same picture emerges from the RDFs in Figure C. The GM1-AQP1 RDFs in the CG simulations show at least two peaks. The first peak at ∼3.4 nm is comparable for all three CG simulations and corresponds to direct GM1-AQP1 interactions. However, the second peak at ∼3.8 nm is decreased significantly in CG PW simulations and corresponds to the GM1 lipids associated with the immediate GM1 shell around AQP1. The RDFs of GM1 headgroups in the CG simulations show the same trend (Figure D). The protein-GM1 RDF in the atomistic simulation (Figure C) shows two peaks at approximately 3.4 and 4.2 nm, respectively, but the amplitudes of these two peaks are much smaller than those in the CG simulations, indicating weaker AQP1-GM1 interactions at the atomistic level. The positions of the peaks of the GM1 tail RDFs are shifted to the right compared with the headgroups in all of the CG and atomistic simulations (Figure D).

Free energy calculations suggest little effects of PW on GM3-AQP1 direct interactions (the binding energies are −19, −24, and −19 kJ/mol in the W, PW, and PW PME CG simulations, respectively (Figure B)), consistent with the equilibrium simulations. These calculations do not give estimation of the real GM3-AQP1 binding free energy, as GM3 was restrained to the specific orientation of AQP1 (Figure S1). However, they provide reliable estimation of the effects of polarizable Martini water.

Despite quantitative and qualitative differences, the CG and atomistic simulations show similar overall binding sites for GM1 on the AQP1 surface, as indicated by the positions of high relative GM1 densities in the 2D density map plots (Figure A). Examples of the GM1-AQP1 interactions in the atomistic simulation are shown in Figure B. There are two GM1 lipids stably bound to AQP1 for the last 1 μs of the simulation, each forming ∼3 hydrogen bonds with AQP1 on average (blue lipids in Figure B). There are other GM1 lipids (magenta lipids in Figure B) interacting with AQP1 indirectly by forming clusters with the stably bound lipids. These lipids may also have transient direct contacts with AQP1 during some of the simulation time. Most interactions between GM1 and AQP1 occur at a specific site, as shown by the density maps in Figure A. The binding site of GM1 lipids on subunit 1 of AQP1 in the atomistic simulation (Figure D–F) is located in the vicinity of the M4 and M7 helices and the M4-M5 loop. Seven residues (Thr205, Asn207, Ser115, Leu121, Asp123, Asn124, and Ser124; Figure D) are involved in hydrogen bonds between GM1 and AQP1, each accounting for ∼10% of the total hydrogen bond numbers. Most of these residues have polar side chains, which can act as hydrogen bond donors and/or acceptors.

Modifications to the Martini Ganglioside Force Field

As shown in section , ganglioside-ganglioside interactions are more prevalent in the CG W simulation than in the atomistic simulation. Therefore, we modified the Martini ganglioside force field to reduce their self-interactions. The aim of the reparameterization was to better reproduce the sizes and structures of ganglioside clusters observed in the atomistic model. Specifically, some of the headgroup bead types were reassigned to reduce the attraction between gangliosides (see Figure for the CG topology of gangliosides, and the mapping of atoms to CG beads). We also made minor modifications to the parameters used for bonded interactions in the headgroups and the linkage parts to improve the internal structure of the GM headgroup.

Figure 7

Modifications to the Martini ganglioside force field. (A) The topology of GM1 and GM3 in the optimized force field. The names, types of CG beads, and the names of the sugars are labeled. (B) Mapping of atoms to CG beads for the GM1 headgroup. (C) The RDF profiles between the GM lipids, headgroups and tails, as well as the RDFs between GM lipids and POPC in the GM-POPC binary mixture simulations using the optimized force field.

Reoptimization of Bonded Interactions

The atomistic trajectories were converted to the CG level. The bond, angle, and dihedral distributions of the resulting GM lipid pseudo-CG trajectories were analyzed and taken as references to optimize the bonded interactions of CG GM lipids (Figures S5 and S6). The original and new parameters for the bonds, angles, and dihedrals are compared in Table S3.

Equilibrium values of the bonds between the AM1-AM2 and AM1-T1A beads were changed to be the same as the corresponding sphingomyelin parameters in the Martini force field (Figure and Table S3) based on their distributions in the pseudo-CG trajectories. Constraints were applied to the bonds between the GM1-GM2 and GM1-GM3 beads of the Glc sugar in the optimized force field. More minor changes were introduced in some other bonds in the headgroups (Table S3).

Angles and dihedrals were modified to improve the orientation of the headgroups and the relative orientation of the sugars, as described in Table S3. The optimization improves the distributions of these angles and dihedrals but does not fit them perfectly to the atomistic simulations (Figures S5 and S6).

Reassignment of Ganglioside Headgroup Bead Types

For GM1, less “sticky” beads were assigned to Neu5Ac. Specifically, beads GM14 and GM15 were changed from type P1 to SP1, and bead GM17 from type P5 to P4. Bead GM6 of sugar Gal1 and GM7 of sugar GalNAc were changed from types Nda and P5 to SNda and P4, respectively. Sugars Glc and Gal2 were kept unchanged. The optimization of the GM1 lipid is based on the fact that hydrogen bond interactions mainly involve Gal2 and Neu5Ac. In the new GM1 force field, Gal2 and Neu5Ac are the most sticky sugars, whereas GalNAc is less sticky than Gal2 and Neu5Ac. Glc and Gal1 are not very sticky (Figure ). The same changes were made to Neu5Ac in GM3; however, beads GM2 and GM5 were changed to type P1 to mimic the tightly and parallel packed ganglioside headgroups in GM3 clusters (Figure ).

The ganglioside-POPC binary mixtures in section were simulated again using these new Martini ganglioside parameters (newFF). The average cluster sizes of GM1 and GM3 in the CG W simulations were now 3 and 4, respectively, whereas using PW reduced the cluster sizes of GM1 and GM3 to 2 and 3, respectively (Table ). The RDFs of ganglioside lipids, headgroups, and tails in Figure also became closer to the atomistic results. The GM1-GM1 RDF is shifted to the right of the GM3-GM3 RDF in the simulations with the newFF. The RDFs of each sugar of the headgroups and the RDFs of specific CG beads in the CG simulations using the original and newFF are shown in Figures S7 and S8. The optimized parameters result in a significant reduction of ganglioside interactions.

The bilayer properties (Table S1) and lipid density profiles (Figure S9) of the newFF CG simulations were also calculated. The bilayer thicknesses are comparable with the atomistic simulations; however, the area per lipid (APL) still differs with the values in the atomistic simulations (Table S1). The diffusion coefficients of the GM lipids are larger than the corresponding values in the simulations with original Martini ganglioside force field due to weaker ganglioside clustering. However, the CG lipid diffusion coefficients are much larger than the atomistic values because of smoother interaction energy at the CG level. The chain order parameters of GM lipids are smaller than those of the atomistic simulations, suggesting the necessity of further optimizing the Martini ceramide building blocks (Table S1; see the Supporting Information for a detailed discussion). The density profiles of lipids are improved in the newFF CG simulations, particular for the density profiles of Gal sugar, which are shifted toward the solution (Figures S2 and S9). Although the new parameters did not exactly reproduce the bilayer properties at the atomistic level, they improved the ganglioside headgroup conformations and the strengths of clustering.

Discussion

Insights into Ganglioside-Ganglioside and Protein-Ganglioside Interactions at the Atomistic Level

Atomistic simulations revealed only weak ganglioside-ganglioside and protein-ganglioside interactions. GM1 and GM3 aggregate by two different strategies. GM3 clusters are more ordered and slightly larger than the GM1 clusters, and the GM3 lipids packed very tightly and almost parallel to each other (Figure ). Sugars at the end of the headgroups (Neu5Ac and Gal2) contribute most to the hydrogen bond interactions, and the ceramide also makes contributions in the case of GM3 (Figure and Table ). This difference might be explained by the shapes and sizes of their headgroups. The two branches (Neu5Ac and GalNAc-Gal2, Figure ) of the GM1 headgroup extend in different directions and contact each other, thereby preventing other parts of GM1 lipids (Glc, Gal1, and the ceramide tails) from approaching each other. The GM3 headgroups adjust their conformations to planar structures, resulting in more tightly packed parallel clusters. Simulations using the CHARMM force field also revealed weak clustering of GM1 (personal communication with Dr. Jeffery Klauda and Dr. Wonpil Im). For protein-ganglioside interactions, although GM1 lipids are depleted from WALP23, the headgroup tilted and interacted with the N-terminus of WALP23 during some of the simulation time (Figure A and C), whereas extensive hydrogen bonding was observed between GM1 and AQP1 (Figure B and D).

Saccharide and protein-saccharide interactions are crucial for molecular/cell recognization;[81] however, investigating these interactions is challenging for both experiments[82] and computations[83] not only because of the diversity of the glycan structures[84] but also because these interactions in solution are weak.[81,83] Although atomistic simulations are the most detailed approach possible at the moment, sampling, in particular inside the bilayer environment, is a significant limitation. In addition, anomeric and exoanomeric effects[85] and CH-π interactions[86] are difficult to model and may affect the reliability of the atomistic simulations of saccharides.

Effects of Polarizable Martini Water (PW) on Ganglioside-Ganglioside and Protein-Ganglioside Interactions at the CG Level

In the CG simulations, both ganglioside (Figures and 3 and Table ) and protein-ganglioside interactions (Figure ) are stronger than in the atomistic simulations. The ganglioside clusters of GM1 and GM3 did not show apparent differences at the CG level.

Previous simulations indicate that the partitioning free energies of charged amino acid residues at the POPC/water interface[87] and the association energies of charged amino acid pairs[88] are improved significantly by PW. Phase transition of bilayers containing GM1 also appear to be modeled better by PW.[56] We found that PW reduced ganglioside and ganglioside-WALP23 interactions in equilibrium simulations and quantified the effects of PW and PW PME by PMF calculations (Figure C and Figure B). We included PW PME simulations, although this is not standard practice with Martini, to compare our results to those of Kociurzynski et al.[56] In simulations of GM1, GM3, and WALP23-GM1, the PW roughly reduces the binding free energies by half. The PW increased GM1-AQP1 binding free energies slightly. However, PME did not significantly affect ganglioside and protein-ganglioside interactions. The changes of GM1-GM1 and GM3-GM3 PMFs when using PW are different, although these lipids have the same charge.

In the protein-ganglioside simulations, PW does not affect the direct GM1-AQP1 interactions significantly (Figure ). This might be explained by the positions of the charged residues of AQP1. For each subunit of AQP1, there are only 2–3 charged residues located on the surface of the protein, whereas the others are buried in the protein and do not interact with the GM1 headgroups (see Figure S10). Thus, the interactions between GM1 and AQP1 are dominated by LJ interactions between the CG beads, and the changes in electrostatic interactions induced by the PW have a minor effect.

Ganglioside Molecular Dynamics

Because atomistic ganglioside simulations at 310 K resulted in the gangliosides entering a gel phase, we performed our simulations at 330 K to maintain the bilayer at the liquid disordered phase; this has the added benefit that lipid diffusion and simulation convergence are faster. A CG W simulation of the GM3-POPC mixture at 330 K revealed similar GM3 cluster sizes and numbers with the simulation at 310 K (Figure S11A compared to Figures B and 3A and B), indicating that the differences of ganglioside clustering at the CG and AA levels are not the result of different temperatures.

In vivo, gangliosides have high phase transition temperatures because of their saturated tails and are probably segregated from lipids with lower phase transition temperatures[24] to form highly ordered lipid microdomains such as lipid rafts with sphingomyelin and cholesterol. The temperature we used in our atomistic simulation prevents phase separation and reduces the effects of different phase transition temperatures of POPC and gangliosides on ganglioside segregation; therefore, interactions between the headgroups play the most important role in ganglioside aggregation in our simulations. The ganglioside clustering is probably decreased compared to the situation in vivo due to the high temperature we used. To evaluate the effects of both headgroup interactions and phase separation on ganglioside aggregation, additional simulations in different lipid environments and at lower temperatures using enhanced sampling techniques might be useful, although in atomistic simulations sampling remains a serious challenge.

One may expect a strain between leaflets as we only included GM lipids in one leaflet and the APLs of GM1, GM3, and POPC are not identical. However, in the CG simulations (using the original Martini ganglioside force field) of pure POPC, symmetrical GM1-POPC, and GM3-POPC bilayers, the averaged APLs are 64.2 ± 0.8, 64.7 ± 0.8, and 64.1 ± 0.8 Å2, respectively, suggesting only minor differences between the APLs of the GM lipids and POPC. Therefore, we ignore the effects of bilayer asymmetry on ganglioside clustering and bilayer properties.

In the CG simulations, the gangliosides tails we used have 3–4 CG tail beads (and one linker bead) each and correspond to GM (d18:1/18:0) lipids at the atomistic level according to the four to one mapping Martini philosophy. However, the GM lipids with longer tails (d20:1/d20:0) were used in our atomistic simulations. To address the inconsistency of the ganglioside tails at the CG and atomistic levels, we conducted control CG simulations using GM lipids with longer tails (AM1-AM2-T1A-C2A-C3A-C4A-C1B-C2B-C3B-C4B-C5B) at 330 K. We found that the lengths of the tails did not influence ganglioside clustering at the CG level in both simulations using the original and modified Martini ganglioside force fields (Figures S11B and C compared to Figures B and 3A and B).

Modifications to the Martini Ganglioside Force Field

We used the extent of ganglioside aggregation and the internal structures of the gangliosides as reference data to modify the Martini ganglioside force field. The detailed atomistic structures of ganglioside clusters were also considered in the reoptimization, but we do not aim to reproduce the different structures of GM1 and GM3 clusters exactly at the CG level. Ganglioside clustering is reduced significantly by the modified force field. This is achieved by assigning new types for the CG beads used for the ganglioside headgroups, suggesting that attraction of the headgroups is an important driving force for ganglioside clustering at the CG level. The parameters of the amide linkage at the CG level are changed to be the same as those used for Martini sphingomyelin for consistency.

We would like to note that different types of CG beads are used for the same building blocks of GM1 and GM3 lipids. Specifically, in GM3, a P1 type is used for bead GM2 in the Glc sugar and GM5 in the Gal1 sugar, but the corresponding types are SP1 in GM1. This is not consistent with the Martini building block philosophy but is needed in this case. Because the two branches of the GM1 headgroup point in different directions (Neu5Ac and GalNAc-Gal2; see Figure ), Glc and Gal1 may not be able to approach each other closely, and the interactions between them are shielded. However, in GM3, all three sugars packed tightly in the clusters. The interactions involving Glc and Gal1 are much stronger than in GM1. Therefore, we used different types of beads for Glc and Gal1 to mimic the cluster sizes and structures in the atomistic simulations. Slightly different dihedrals are used for GM1 and GM3 due to the different orientations of the sugars of their headgroups (Table S3).

It has been reported that membrane incorporation restricts the motion of ganglioside headgroups.[89] The Glc sugar was buried in the interfacial region of the phosphatidylcholine bilayers, whereas the other sugars were extended fully into solution.[40,90] The density profiles of the sugars in our atomistic and CG newFF simulations reproduced these experimental results (Figures S2 and S9). However, the GM clustering behavior in our simulations should ultimately be subjected to more rigorous testing, including experimental data that is currently lacking. The new GM parameters are available on the Martini Web site (cgmartini.nl) and in the CHARMM-GUI Martini maker.[91]

Conclusions

In atomistic simulations, the sugars at the end of the headgroups (Neu5Ac and Gal2) contribute most of the hydrogen bonds between GM headgroups. GM1 and GM3 pack differently in their respective clusters due to the different headgroup geometry with GM3 lipids packing more tightly, leading to GM3 clusters being more stable and a little larger than GM1 clusters. In atomistic simulations, GM1 binds stably to specific positions of AQP1 by hydrogen bonding between polar residues of the protein and the GM1 headgroups. The GM1 lipids are depleted from WALP23.

The ganglioside and protein-ganglioside interactions in the CG Martini simulations are more prevalent than in the atomistic simulations. Using the polarizable Martini water decreases ganglioside and protein-ganglioside interactions in most of the conditions. We proposed an optimized Martini force field for gangliosides that reproduces the strengths of ganglioside aggregations at the atomistic level.