Overcoming the limitations of the MARTINI force field in Molecular Dynamics simulations of polysaccharides

TL;DR

This study evaluates the MARTINI force field for simulating polysaccharides in molecular dynamics, identifies its overestimation of aggregation, and proposes a scaling adjustment to improve accuracy.

Contribution

It demonstrates the limitations of MARTINI in polysaccharide simulations and introduces a simple scaling method to correct aggregation issues.

Findings

MARTINI overestimates polysaccharide aggregation in simulations.

Comparison with experiments shows non-physical aggregation.

Scaling Lennard-Jones interactions improves simulation reliability.

Abstract

Polysaccharides (carbohydrates) are key regulators of a large number of cell biological processes. However, precise biochemical or genetic manipulation of these often complex structures is laborious and hampers experimental structure-function studies. Molecular Dynamics (MD) simulations provide a valuable alternative tool to generate and test hypotheses on saccharide function. Yet, currently used MD force fields often overestimate the aggregation propensity of polysaccharides, affecting the usability of those simulations. Here we tested MARTINI, a popular coarse-grained (CG) force field for biological macromolecules, for its ability to accurately represent molecular forces between saccharides. To this end, we calculated a thermodynamic solution property, the second virial coefficient of the osmotic pressure ($B_{22}$). Comparison with light scattering experiments revealed a non-physical aggregation of a prototypical polysaccharide in MARTINI, pointing at an imbalance of the non-bonded solute-solute, solute-water, and water-water interactions. This finding also applies to smaller oligosaccharides which were all found to aggregate in simulations even at moderate concentrations, well below their solubility limit. Finally, we explored the influence of the Lennard-Jones (LJ) interaction between saccharide molecules and propose a simple scaling of the LJ interaction strength that makes MARTINI more reliable for the simulation of saccharides.

Figures

Tables

Table 1. Table 1: Experimental values of B 22 subscript 𝐵 22 B_{22} for diverse saccharides

Molecule	Mw [g/mol]	Solvent	T [K]	$B_{22}$ [L/mol]	$A_2$ [mol L/g2]	Ref.
glucose	180.16	water	298.15	0.117	$3.61\text{×}{10}^{-6}$	25
cellobiose	342.30	water	298.15	0.267	$2.28\text{×}{10}^{-6}$	64
sucrose	342.30	water	298.15	0.305	$2.60\text{×}{10}^{-6}$	25
trehalose	342.30	water	295	0.51	$4.3\text{×}{10}^{-6}$	65
$\alpha$-CD	972.85	water	298.15	0.830	$6.57\text{×}{10}^{-7}$	66
$\beta$-CD	1135	water	298.15	6.296	$4.89\text{×}{10}^{-6}$	67
A2 glycan	2224	0.1M NaCl	300	46	$9.3\text{×}{10}^{-6}$	this work
dextran	9000	0.01M NaN3	293.15	60.7	$7.49\text{×}{10}^{-7}$	63
dextran	37400	0.01M NaN3	293.15	590	$4.22\text{×}{10}^{-7}$	63
dextran	59000	0.01M NaN3	293.15	1590	$4.56\text{×}{10}^{-7}$	63
dextranT2000	$1.58\text{×}{10}^6$	0.03M NaCl	302.05	$7.5\text{×}{10}^5$	$3.0\text{×}{10}^{-7}$	68

Table 2. Table 2: Comparison of simulated and experimental B 22 subscript 𝐵 22 B_{22} values for the five studied saccharides. RDF-derived values (symmetric errors) represent means, PMF-derived values (asymmetric errors) medians. Errors always denote 95% confidence interval ( ± 2 plus-or-minus 2 \pm{2} SD).

	$B_{22}$ (L/mol)
	simulated		experimental
Saccharide	$\gamma =1.0$	$\gamma =0.5$
glucose	$-1.2\pm 0.3$	$0.22\pm 0.23$	0.117
sucrose	$-{11}_{-6}^{+4}$	$0.57\pm 0.55$	0.305
$\alpha$-CD	$-{8.8}_{-5.1}^{+3.3}$$\times$${10}^{07}$	$-3.0\pm 4.6$	0.830
$\beta$-CD	$-{1.1}_{-0.2}^{+0.2}$$\times$${10}^{03}$	$-4.2\pm 6.9$	6.296
A2 glycan	$-{1.3}_{-2.4}^{+0.8}$$\times$${10}^{06}$	$-3.0\pm 8.7$	$46\pm 11$

Table 3. Table S1: MARTINI bead types for coarse-grained carbohydrates

saccharide	bead no.a	bead namea,b	bead typec
$\beta$-D-glucose	1	B3	GP4
	2	B2	GP4
	3	B1	GP1
sucrose	1	B1	GP1
	2	B2	GP2
	3	B3	GP4
	4	B4	GP1
	5	B5	GP1
	6	B6	GP4
$\alpha$-cyclodextrin	1	B1	GP1
	2	B2	GP2
	3	B3	GP4
$\beta$-cyclodextrin	1	B1	GP1
	2	B2	GP2
	3	B3	GP4
A2 glycan	1	Gn11	GNa
	2	Gn12	GP4
	3	Gn13	GSP1
	4	Gn14	GSP1
	5	Gn21	GP5
	6	Gn22	GSP1
	7	Gn23	GP1
	8	bM1	GP1
	9	bM2	GSP1
	10	bM3	GNda
	11	aM31	GNda
	12	aM32	GP4
	13	aM33	GSP1
	14	aM61	GNda
	15	aM62	GP4
	16	aM63	GSP1
	17	Gn31	GP5
	18	Gn32	GNda
	19	Gn33	GSP1
	20	Ga31	GSP1
	21	Ga32	GSP1
	22	Ga33	GP1
	23	SA31	GQa
	24	SA32	GSP1
	25	SA33	GSP1
	26	SA34	GP5
	27	SA35	GP4
	28	Gn61	GP5
	29	Gn62	GNda
	30	Gn63	GSP1
	31	Ga61	GSP1
	32	Ga62	GSP1
	33	Ga63	GP1
	34	SA61	GQa
	35	SA62	GSP1
	36	SA63	GSP1
	37	SA64	GP5
	38	SA65	GP4

Table 4. Table S2: Force Field Parameters for coarse-grained carbohydrates a

saccharide	bond	$r_{eq}$	$k_r$	angle	${\theta }_{eq}$	$k_{\theta }$	dihedral	${\varphi }_{eq}$	$k_{\varphi }$
		[nm]	[$\frac{kJ}{mol\cdot n{m}^2}$]		[deg]	[$\frac{kJ}{mol\cdot {\left[rad\right]}^2}$]		[deg]	[$\frac{kJ}{mol\cdot {\left[rad\right]}^2}$]
$\beta$-d-glucose	B3-B2	0.323	$\infty$
	B3-B1	0.384	$\infty$
	B2-B1	0.331	$\infty$
sucrose	B1-B2	0.222	$\infty$	B1-B2-B4	130	10	B1-B2-B4-B5	130	25
	B2-B3	0.247	$\infty$	B3-B2-B4	110	150	B1-B2-B4-B6	80	2
	B2-B4	0.429	$\infty$	B5-B4-B2	20	50	B3-B2-B4-B5	-70	20
	B4-B5	0.293	$\infty$	B6-B4-B2	85	150
	B4-B6	0.372	$\infty$
$\alpha$-cyclodextrin	G11-G12	0.215	5500	G11-G12-G22	80	1
	G12-G13	0.220	2500	G12-G22-G32	120	28
	G12-G22	0.460	450	G23-G22-G12	94	100
$\beta$-cyclodextrin	G11-G12	0.215	10000	G11-G12-G22	70	5
	G12-G13	0.220	$\infty$	G12-G22-G32	135	30
	G12-G22	0.470	500	G23-G22-G12	90	120
A2 glycan	Gn12-Gn13	0.270	8000	Gn11-Gn12-Gn13	120	250	Gn23-Gn21-Gn13-Gn12	180	20
	Gn13-Gn14	0.340	12000	Gn11-Gn12-Gn14	147	500	bM3-bM1-Gn22-Gn21	80	10
	Gn12-Gn14	0.380	20000	Gn12-Gn13-Gn21	154	400	aM33-aM31-bM2-bM1	-10	40
	Gn13-Gn21	0.433	5000	Gn14-Gn13-Gn21	72	250	aM63-aM61-bM3-bM1	180	40
	Gn21-Gn22	0.378	10000	Gn13-Gn21-Gn22	86	600	Gn33-Gn31-aM31-aM33	170	40
	Gn21-Gn23	0.524	22000	Gn13-Gn21-Gn23	50	300	Ga33-Ga31-Gn32-Gn31	-160	10
	Gn22-bM1	0.338	15000	Gn21-Gn22-bM1	168	500	SA33-SA31-Ga32-Ga31	-50	5
	bM1-bM3	0.334	18000	Gn23-Gn22-bM1	90	50	Gn63-Gn61-aM61-aM63	160	35
	bM2-aM31	0.372	9000	Gn22-bM1-bM2	131	700	Ga63-Ga61-Gn62-Gn61	110	8
	aM32-aM33	0.322	16000	Gn22-bM1-bM3	78	200	SA63-SA61-Ga62-Ga61	20	30
	aM31-aM33	0.349	8500	bM1-bM2-aM31	105	200
	bM3-aM61	0.354	8000	bM3-bM2-aM31	169	600
	aM62-aM63	0.323	16000	bM2-aM31-aM32	110	300
	aM61-aM63	0.354	15000	bM2-aM31-aM33	86	200
	aM31-Gn31	0.340	4500	bM2-aM31-Gn31	122	60
	Gn31-Gn32	0.388	20000	aM33-aM31-Gn31	90	80
	Gn32-Gn33	0.315	20000	aM31-Gn31-Gn32	97	550
	Gn31-Gn33	0.516	15000	aM31-Gn31-Gn33	61	400
	Ga32-Ga33	0.312	12000	Gn31-Gn32-Ga31	167	700
	Ga32-SA31	0.335	2000	Gn33-Gn32-Ga31	78	250
	SA31-SA33	0.382	15000	Gn32-Ga31-Ga32	112	500
	SA32-SA34	0.357	15000	Gn32-Ga31-Ga33	71	350
	SA33-SA35	0.296	10000	Ga31-Ga32-SA31	108	180
	aM61-Gn61	0.340	4800	Ga33-Ga32-SA31	137	50
	Gn61-Gn62	0.388	20000	Ga32-SA31-SA32	100	250
	Gn62-Gn63	0.315	20000	Ga32-SA31-SA33	80	220
	Gn61-Gn63	0.516	15000	SA31-SA32-SA34	122	800
	Ga62-Ga63	0.320	12000	SA31-SA33-SA35	109	350
	Ga62-SA61	0.335	2000	SA32-SA33-SA35	157	400
	SA61-SA63	0.382	15000	SA33-SA32-SA34	70	80
	SA62-SA64	0.357	15000	bM1-bM3-aM61	90	150
	SA63-SA65	0.294	10000	bM2-bM3-aM61	117	150
	Gn11-Gn12	0.280	4000	bM3-aM61-aM62	118	400
	Gn22-Gn23	0.312	$\infty$	bM3-aM61-aM63	83	50
	bM1-bM2	0.276	$\infty$	bM3-aM61-Gn61	119	190
	bM2-bM3	0.307	$\infty$	aM61-Gn61-Gn62	80	100
	aM31-aM32	0.277	$\infty$	aM61-Gn61-Gn63	61	600
	aM61-aM62	0.277	$\infty$	Gn61-Gn62-Ga61	140	50
	Gn32-Ga31	0.365	$\infty$	Gn63-Gn62-Ga61	78	250
	Ga31-Ga32	0.269	$\infty$	Gn62-Ga61-Ga62	112	600
	Ga31-Ga33	0.399	$\infty$	Gn62-Ga61-Ga63	71	300
	SA31-SA32	0.337	$\infty$	Ga61-Ga62-SA61	107	230
	SA32-SA33	0.315	$\infty$	Ga63-Ga62-SA61	135	200
	Gn62-Ga61	0.366	$\infty$	Ga62-SA61-SA62	99	170
	Ga61-Ga62	0.269	$\infty$	Ga62-SA61-SA63	80	170
	Ga61-Ga63	0.399	$\infty$	SA61-SA62-SA64	122	700
	SA61-SA62	0.336	$\infty$	SA61-SA63-SA65	109	200
	SA62-SA63	0.315	$\infty$	SA62-SA63-SA65	157	350
				SA63-SA62-SA64	70	80
excluded non-bonded	Gn13-Gn23
interactions	Gn14-Gn21
	Gn22-bM3
	Gn23-bM1
	bM1-aM61
	bM3-aM63
	aM31-Gn33
	aM33-Gn31
	aM61-Gn63
	Gn32-Ga33
	SA33-SA34
	Gn62-Ga63
	SA63-SA64

Equations

$$V_{LJ}(r_{ij})=4{ϵ}_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]$$

$$\Pi (T,c)=RT(c+B_{22}{c}^2+B_{23}{c}^3+\dots )$$

$$B_{22}=-2\pi N_A{\int }_0^{\infty }[\exp (-\frac{W_2(r)}{RT})-1]r^{2}dr$$

$$g(r)\approx \exp \left(-\frac{W_2(r)}{RT}\right)$$

$$B_{22}(r^{\prime })=-2\pi N_A{\int }_0^{r^{\prime }}[g(r)-1]r^{2}dr$$

$$g(r)=\frac{\mathrm{d}N(V)}{\rho \mathrm{d}V}=\frac{1}{4\pi r^2\rho }\frac{\mathrm{d}N(r)}{\mathrm{d}r}\text{with }V(r)=\frac{4}{3}\pi r^3$$

$$B_{22}(r^{\prime })=N_A\left(\frac{2\pi }{3}{r}^{\prime 3}-\frac{N(r^{\prime })}{2\rho }\right)$$

$$B_{22}(r^{\prime })=-2\pi N_A{\int }_0^{r^{\prime }}[\exp (-\frac{W_2(r)}{RT})-1]r^{2}dr$$

$$g^{(2)}(\tau )=\frac{⟨I(t)I(t+\tau )⟩}{⟨I(t){⟩}^2}$$

$$g^{(2)}(\tau )=[g^{(1)}(\tau ){]}^2+\xi (\tau )+1$$

$$\begin{array}{rl}{g}^{(1)}(\tau )& ={\int }_0^{\infty }B(\Gamma ){\mathrm{e}}^{-\Gamma \tau }\mathrm{d}\Gamma ,\\ \text{with }\Gamma & =D_t{q}^2,\\ \text{and }q& =4\pi \frac{n}{{\lambda }_0}\sin \frac{\theta }{2}\end{array}$$

$$R_H=\frac{k_{B}T}{6\pi {\eta }_0{D}_t}=\frac{k_{B}T{q}^2}{6\pi {\eta }_0\Gamma }\iff \Gamma =\frac{k_{B}T{q}^2}{6\pi {\eta }_0{R}_H}$$

$$C_j\propto \frac{B_j(R_H,\theta )}{M_j(R_H)P(R_H,\theta )}$$

$$\begin{array}{rl}⟨R_{\theta }⟩\equiv \frac{⟨I_s(\theta )⟩}{⟨I_0(\theta )⟩}\frac{r^2}{V_s(\theta )}& =\frac{4{\pi }^2}{N_A}\frac{n_0^2}{{\lambda }_0^4}{\left(\frac{\mathrm{d}n}{\mathrm{d}C}\right)}^2{M}_{W}CS(\theta )\\ & \equiv K{M}_{W}CS(\theta )\end{array}$$

$$\frac{1}{S}=1+2A_2{M}_{W}C+3A_3{M}_W{C}^2+\dots$$

$$\frac{\Pi }{RT}=\frac{1}{M_W}C+A_2{C}^2+A_3{C}^3+\dots$$

$$\frac{KC}{R_{\theta }}=\frac{1}{M_W}+2A_{2}C+3A_3{C}^2+\dots$$

$$B_{22}=A_2{M}_W^2,B_{23}=A_3{M}_W^3,\dots$$

$${ϵ}_{ij,scaled}=2\text{ kJ/mol}+\gamma ({ϵ}_{ij,original}-2\text{ kJ/mol})$$

Full text

I INTRODUCTION

Polysaccharides

Polysaccharides are sugar polymers found in various biological contexts, e.g. glycoproteins, proteoglycans, bacterial lipopolysaccharides, and cell walls of plants and fungi, functioning as biomolecular interaction modulatorsJohnson et al. (2005); Langer et al. (2012); Moremen et al. (2012), structural elements and energy storage. Their structural diversity ranges from simple linear homo-polymers, like cellulose, to cyclic or branched structures composed of diverse sugars connected via glycosidic linkagesVarki and Sharon (2009) (fig. 1). This complexity and the often encountered microheterogeneity, i.e. the simultaneous occurrence of length and structure variants of a given polysaccharide, hamper experimental approaches to study the roles of polysaccharides in biological processes. This gap can be filled by molecular dynamics (MD) simulations if the representation of polysaccharides in the respective model accurately captures their biophysical properties. Many carbohydrate-specific force fields have been developed (see Foley et al.Foley et al. (2012) for a review) and applied to study the behavior of saccharides in solutions as well as their interactions with other biomoleculesFeng et al. (2015); Xiong et al. (2015); Sauter and Grafmüller (2017). In this regard, the non-bonded interactions are of particular importance as they determine the magnitude of intermolecular forces. In MD force fields they are typically represented by Lennard-Jones (LJ) and Coulomb potentials for van der Waals and electrostatic interactions, respectively.

The MARTINI coarse-grained force field

The goal of MD simulations is to extract biophysical properties from a system that is incrementally evolved over time. This essentially means averaging over stochastic events, which requires system size and simulation length to be sufficiently large. Simulations at atomistic level (all-atom, AA) require femtosecond time increments to model bond stretching and so they are typically limited to system sizes on the order of $10^{5}$ atoms and time scales below 1 $\mu$s. This is insufficient for many processes involving polysaccharides in terms of both size and time. A way to overcome these limitations is the use of coarse-grained (CG) force fields which sacrifice atomic-resolution detail to reduce computational cost. The usual strategy involves replacing groups of atoms with larger pseudo-atoms (beads) that retain averaged properties of the underlying atomistic system. This reduces the number of particles in the system and smooths the free energy landscape, increasing the effective time span covered by the simulation.

MARTINIMarrink et al. (2007) is a CG force field that maps groups of four neighbor ”heavy” atoms (C, O, N, P, S) onto pre-defined beads. This rather modest level of coarse-graining retains much structural detail while offering substantial speed-up compared to all-atom simulations. Although originally devised for lipid systems, MARTINI has been extended to proteinsMonticelli et al. (2008), polysaccharidesLópez et al. (2009) and nucleic acidsUusitalo et al. (2015). When transferring an atomistic system to MARTINI, parameters for bonded interactions between the beads have to be found empirically, e.g. by comparison to an AA simulation, while non-bonded interactions are predetermined by bead type and charge. MARTINI beads are not endowed with partial charges of the component heavy atoms, so non-bonded interactions of uncharged molecules are determined by LJ interactions only. The LJ potential $V_{LJ}$ between two beads $i$ and $j$ at a distance $r_{ij}$ takes the following shape:

[TABLE]

The finite distance at which the potential reaches zero, $\sigma_{ij}$, and the depth of the potential well, $\epsilon_{ij}$, vary with the bead type and have been fit to reproduce partition coefficients of small reference molecules in polar/apolar solvent systemsMarrink et al. (2007).

The problem of non-bonded interactions

This common strategy of simply extrapolating from small molecule interaction energies to the macromolecules of interest bears the risk of misrepresenting the strength of interactions between macromolecules due to multiplication of small systematic errors in the parametrization. Indeed, it has been found in recent years that AA simulations tend to overestimate the aggregation propensity of proteins Petrov and Zagrovic (2014); Best et al. (2014); Henriques et al. (2015) and polysaccharidesSauter and Grafmüller (2015). The widely accepted notion is that this is due to an imbalance of solute-solute, solute-water, and water-water interactions which are determined by the choice of the water model and the solute-solute interaction potentials defined in the force field. To alleviate this imbalance, it has been proposed to include experimentally addressable solution properties in the parametrization process, such as the osmotic pressureLuo and Roux (2010); Yoo and Aksimentiev (2012, 2016), Kirkwood-Buff integralsPloetz and Smith (2011); Karunaweera et al. (2012) or the osmotic coefficientMiller et al. (2016, 2017). All these parameters relate to the second virial coefficient, $B_{22}$, of the osmotic pressure, $\Pi$, which describes the deviation from ideal behavior of a solution with solute molar concentration $c$McMillan Jr. and Mayer (1945):

[TABLE]

Here, $T$ denotes temperature, $R$ is the gas constant and $B_{ij}$ are coefficients of the virial expansion, with $i=2$ referring to the solute in a binary mixture and $j$ enumerating consecutive coefficients. $B_{22}>0$ indicates net attraction and $B_{22}<0$ repulsion between molecules, whereas the magnitude corresponds to the aggregation propensity. Experimentally, $B_{22}$ can be determined e.g. through direct measurements of $\Pi$Stigter (1960), self-interaction chromatographyTessier et al. (2002), or diffraction experimentsGeorge et al. (1997). These methods are complemented by an established theoretical groundworkMcMillan Jr. and Mayer (1945) that allows for calculation of $B_{22}$ from MD simulations, making it a powerful tool for the refinement of intermolecular interactions. Such refinement has been recently performed for several force fieldsBlanco et al. (2013); Nikoubashman et al. (2015); Kim and Hummer (2008); Grünberger et al. (2013); Calero-Rubio et al. (2016), including a study from the Elcock groupStark et al. (2013), who found an abnormally strong aggregation propensity of proteins in MARTINI that could be remedied by reduction of $\epsilon_{ij}$ of solute-solute LJ interactions. We hypothesize that a similar problem exists for polysaccharides in MARTINI, and furthermore, due to a cumulative effect of many interacting atoms, this spurious aggregation propensity is expected to grow with the size of interacting molecules. Should this be the case, it would pose a threat that unrealistic aggregation behavior could influence outcomes of studies based on the present MARTINI model for saccharides (e.g.Ma et al. (2015); López et al. (2015); Kapla et al. (2016)). In this work we compute $B_{22}$ from MARTINI simulations of five different saccharides (1 to 11 residues) and compare these to experimental values. To this end we utilize light scattering experiments to determine $B_{22}$ of a complex branched polysaccharide prototypical for protein glycosylation. We demonstrate that MARTINI considerably overestimates the aggregation behavior of even small saccharides and show that uniform scaling of the LJ parameter $\epsilon_{ij}$ is sufficient to facilitate more realistic outcomes of MARTINI simulations.

II METHODS

II.1 Computational methods

II.1.1 MD simulations

All simulations in this work were performed using GROMACSAbraham et al. (2015) 5.0.4 for enhanced sampling simulations and 5.1 otherwise.

All atom simulations

All atom simulations of the carbohydrates glucose, sucrose, $\alpha$-cyclodextrin ($\alpha$-CD), $\beta$-cyclodextrin ($\beta$-CD), and A2 glycan were conducted using the GLYCAM06Kirschner et al. (2008) force field and initial structures and carbohydrate-specific force field parameters were obtained from the GLYCAM06j-1 carbohydrate builder suiteKirschner et al. (2008) except for CDs, which were manually prepared based on the corresponding amylohexa- and -heptaose structures. Topologies were converted to GROMACS format using ACPYPESousa da Silva and Vranken (2012). For each of the saccharides, the initial structure was placed in a cubic box large enough to prevent self-interactions, subjected to energy minimization and solvated with SPC waterBerendsen et al. (1981). Subsequently, a number of water beads was replaced with sodium and chloride ions in order to neutralize any charges and to obtain a salt concentration of approximately physiological 100 mM NaCl. Finally, the system was equilibrated in the canonical ensemble, (i.e. with constant volume and temperature) for 250 ps. Production runs were performed in the isothermal-isobaric ensemble, with Nosé-HooverNosé (1984) thermostat and Parrinello-RahmanParrinello and Rahman (1981) barostat, with long range electrostatic interactions calculated by the Particle-Mesh Ewald method and short-range electrostatics and van der Waals interactions cut off at 1.0 nm, as recommended for the GLYCAM06 force fieldKirschner et al. (2008). All simulations were performed at a temperature of 300 K using an integration step of 2 fs. Trajectories of 150 ns were used for mapping of bonded interactions.

Parametrization of saccharides in the MARTINI force field

We generally followed the proposed parametrization scheme for monosaccharides and linear oligo- and polysaccharidesLópez et al. (2009). A monosaccharide unit is therein represented by three ,,polar” beads which can be of type $\mathrm{P_{1}}$, $\mathrm{P_{2}}$, $\mathrm{P_{3}}$, or $\mathrm{P_{4}}$. The differences in bead polarity, denoted by the subscripts, are reflected in different strengths of the LJ interaction parameter $\epsilon$. The full interaction matrix is found in the original MARTINI force field publicationMarrink et al. (2007). For glucose and sucrose we adopted the published parameter setLópez et al. (2009). Glucose therein consists of a triangle of three beads, two more polar $\mathrm{P_{4}}$ beads and one less polar $\mathrm{P_{1}}$ bead (see fig. 1B and Tables S1 and S2 for topology and corresponding bead types). The two sugar units in sucrose are represented by a different topology. A central bead connects to two outer beads that are not inter-connected. The monosaccharide units are linked through a bond between their central beads. The employed bead types are $\mathrm{P_{1}}$, $\mathrm{P_{2}}$, and $\mathrm{P_{4}}$ (fig. 1C). For the cyclic oligosaccharides $\alpha$- and $\beta$-CD (fig. 1D and E) the parametersLópez et al. (2009) for maltose ($\alpha$-glucosyl-(1$\rightarrow$4)-$\alpha$-glucose) and previous work on $\beta$-CD López et al. (2013a) served as template for a mapping procedure using AA simulations as a reference.

To this end, we iterated through short 50 ns CG simulations (described below), extracted bond and angle distributions and made small adjustments to bond and angle (and dihedral, see below) force field parameters until satisfactory agreement with reference distributions obtained from a posteriori coarse grained AA simulations had been reachedLópez et al. (2009). The same procedure was used to parametrize A2 glycan (fig. 1A), for which the mapping rules put forward for oligosaccharidesLópez et al. (2009) together with the parameters of the glycolipid GM1López et al. (2013b) provided the starting point. Herein, the additional bead types $\mathrm{SP_{1}}$, $\mathrm{P_{5}}$, $\mathrm{N_{a}}$, $\mathrm{N_{da}}$, and $\mathrm{Q_{a}}$ were used. Unfortunately, the introduction of dihedral potentials for branched A2 glycan caused instability of simulations unless a time step of 5 fs instead of 30 fs was used, as has been noted beforeLópez et al. (2013b). The dihedral potentials were therefore used for sucrose, but disabled in A2 glycan parameters. Nevertheless the glycosidic dihedrals, being most instructive for the conformation of A2 glycan, in most cases show qualitatively similar behavior in MARTINI and AA simulations (see the supplementary files for comparison of AA to CG distributions with dihedral potentials switched on and off). The complete parameter set for the bonded interactions is listed in table S2.

Coarse-grained simulations

For all coarse-grained simulations the MARTINI force field version 2.2 was used, with modifications as described in the Results section. Briefly, the martinize script (v2.4), obtained from the MARTINI websitede Jong et al. (2013) and supplemented with mapping schemes and bonded parameters of saccharides studied in this work, was used to convert atomistic structures into coarse-grained ones and to generate the necessary topology files. Similarly to AA simulations, structures were energy-minimized, solvated with MARTINI water containing 100 mM NaCl and 10 mM $\mathrm{CaCl_{2}}$ and equilibrated in the canonical ensemble for a total of 1 ns. Additional tests showed that the inclusion of calcium ions, not thoroughly tested in MARTINI, had no influence on the aggregating properties of saccharides (see fig. S1B). Production runs were performed in the isothermal-isobaric ensemble using velocity-rescaleBussi et al. (2007) and Parrinello-Rahman coupling schemes to keep temperature and pressure constant. For all MARTINI simulations a time step of 30 fs was used. To allow for GPU-accelerated simulations, the recently described Verlet neighbor search algorithm was usedde Jong et al. (2016). Briefly, both LJ and Coulomb potentials were cut off beyond 1.1 nm. The LJ potential was in addition shifted to zero at the cutoff distance. For Coulomb interactions, the reaction-field potential was used. Details of the parameters were kept according to the recommended mdp files (http://cgmartini.nl). Water is modeled in MARTINI as uncharged $\mathrm{P_{4}}$ beads with each bead representing four water moleculesMarrink et al. (2007). Special ”antifreeze” beads have been introduced to avoid the freezing of water at temperatures around 300 KMarrink et al. (2007). Unless otherwise stated 10% of the water particles (W) are replaced by antifreeze particles (WF) (referred to as antifreeze water) in the course of this study.

II.1.2 Computation of $B_{22}$

Two methods were employed to calculate $B_{22}$ from MD simulations, based on either the cumulative solute-solute radial distribution function (RDF), or a reconstruction of the potential of mean force (PMF) between two solute molecules over their separation distance.

Cumulative RDF method

McMillan and MayerMcMillan Jr. and Mayer (1945) derived an expression for $B_{22}$ under the assumption that the total solute potential energy can be approximated as the sum of pairwise solute-solute interactions, i.e. the potential of mean force (PMF) $W_{2}(r)$ between two solute particles at distance $r$:

[TABLE]

where $N_{A}$ denotes Avogadro’s constant. For weakly interacting solute particles, $B_{22}$ can be calculated efficiently by means of the radial distribution function (RDF), $g(r)$, which in thermodynamic equilibrium is approximately Boltzmann-distributed according to $W_{2}(r)$:

[TABLE]

Inserting (4) into (8) gives:

[TABLE]

Equation (5) can be simplified by making use of the definition $g(r)$ as the solute molecule number increment d$N$ in the spherical shell d$V$ at distance $r$, normalized by the average particle density $\rho=n/V$:

[TABLE]

Inserting (6) into (5) yields $B_{22}$ in terms of the cumulative number distribution function $N(r^{\prime})$:

[TABLE]

which obviates the need for numerical integration.

To be of practical use, the RDF method requires sufficient solute molecule numbers over the separation distance of interest, a prerequisite that is only fulfilled when the free energy landscape, i.e. the PMF, does not contain wells so deep that the Boltzmann-distributed solute molecules are effectively depleted from other regions. We found this limit to be on the order of -1 $RT$ (-2.5 kJ/mol), therefore if trial MD simulations for a given condition led to PMF minima not lower than -2.5 kJ/mol we employed the cumulative RDF method and otherwise the HEUS method, as detailed below.

Where appropriate, RDFs were obtained from MARTINI simulations of $n=420$ (glucose, sucrose) or $n=105$ (CDs, A2 glycan) randomly placed saccharide molecules in water-filled cubic boxes of volume [$V_{b}=(19\;\mathrm{nm})^{3}$]. These numbers, corresponding to concentrations of 100 and 25 mM respectively, ensure a sufficient degree of solute-solute interaction while being well below the solubility limit (except for $\beta$-CD which is soluble up to 16.5 mM in water).

The resulting trajectories of 1 to 10 $\mu$s were split into 200 ns segments and for every segment the cumulative saccharide-saccharide RDF (in terms of centers of mass) was computed using the GROMACS *gmx rdf *programAbraham et al. (2015). The resulting curves were averaged to give the final RDF with standard deviation as an error estimate.

Direct calculation of PMF

In principle an MD simulation of two saccharide molecules could yield the PMF along the reaction coordinate of choice, that is the center of mass separation of the two molecules. However, in the case of large energy barriers along the reaction coordinate, (precluding efficient calculation of the RDF) mere Boltzmann sampling does not suffice to adequately sample the free energy landscape. In umbrella samplingTorrie and Valleau (1977) numerous individual MD simulations are run in parallel, each imposed with an additional biasing potential designed to restrain the reaction coordinate to a certain window with a given force constant. Afterwards the individual potentials are de-biased and used to reconstruct the final PMF. An additional improvement that mitigates entrapment in local minima is the stochastic exchange of the Hamiltonian, that is effectively the biasing potential, between neighboring windows (Hamiltonian exchange umbrella sampling, HEUS)Bussi (2014).

In order to perform HEUS simulations, GROMACS v5.0.4 was patched with the PLUMED v2.2.1 plug-inTribello et al. (2014). The initial setup consisted of two solute molecules placed at maximum distance along the longest axis of an $8\times 8\times 12$ nm3 simulation box. After minimization and solvation, as described above, molecules were pulled towards each other while recording intermediate positions every 0.15 nm, until sterical clashes prevented further motion. The centers of harmonic biasing potentials were set at these positions with uniform force constant of 500 $\mathrm{kJ/(mol\cdot nm^{2})}$, spanning typically distances between 0.2 and 5 nm and yielding up to 32 independent umbrella windows. Each window was further subjected to a short equilibration and subsequent production runs of 300 to 500 ns were carried out. Every 1000 steps an attempt was made to exchange Hamiltonians of neighboring windows. The center of mass positions and the corresponding value of biasing potential was recorded every 10 steps.

During pulling and production simulations, molecules were free to move along the long box axis and rotate, but the center of mass position along shorter axes was restrained with a force constant of $1\text{\times}{10}^{5}$ $\mathrm{kJ/(mol\cdot nm^{2})}$ to prevent dumbbell-like rotation of paired molecules across the box. Since this additional restraint was orthogonal to the direction in which the PMF was calculated and the rotational degrees of freedom were not affected, we expect this did not have an effect on the PMFs. The same restraint was applied in all dimensions during equilibration steps.

In order to construct the PMF, HEUS simulations were de-biased using the Weighted Histogram Analysis Method (WHAM)Kumar et al. (1992) implemented in GROMACSHub et al. (2010). Each simulation was repeated at least seven times with randomized molecule orientations in the starting configuration. Average and standard deviation were estimated from a bootstrapping procedure (100 cycles) using replica histograms to reconstruct random PMFs. The final PMFs were offset to zero energy at an intermolecular distance of 4 nm, where a plateau was reached for all molecules studied.

For the purpose of $B_{22}$ calculation, the integral in (3) has to be finite:

[TABLE]

and an $r^{\prime}=r^{\prime}_{0}$ has to be chosen large enough so that solute-solute interactions effectively vanish and $B_{22}(r^{\prime}_{0})\rightarrow B_{22}(\infty)$. We found this consistently to be the case for $r^{\prime}>4$ nm, therefore we settled at a (conservative) $r^{\prime}_{0}$ = 5 nm in the course of this study. Numerical integration of (8) was performed using the trapezoid rule and a set of homemade scripts. Asymmetric 95% confidence intervals were calculated, exploiting the fact that $B_{22}$ values obtained from bootstrapped PMFs were distributed log-normallyLimpert et al. (2001). PMFs derived with both RDF and HEUS methods under the same conditions were found remarkably similar (fig. 6).

II.2 Experimental methods

II.2.1 Biophysical characterization of A2 in aqueous solution

A2 glycan (1 mg) was obtained as a powder from Ludger Ltd, Oxfordshire, UK (cat. no. CN-A2-SPBULK, batch no. B623-01). The manufacturer had determined its purity, by hydrophilic interaction liquid chromatography of 2-aminobenzamide labeled A2 glycan, to be 92.0% with 3.7% and 0.6% contaminants resulting from loss of one or both sialic acid moieties, respectively.

II.2.2 Solution dispersity and hydrodynamic radius

To assess the size distribution of dissolved A2 glycan molecules, solutions of 4 to 10 g/L in 100 mM NaCl were prepared, filtered (Whatman Anotop, 0.02 $\mu$m pore size), and dynamic light scattering (DLS) intensities were recorded in a DynaPro NanoStar instrument equipped with a calibrated 1.25 $\mu$L quartz cuvette (Wyatt Technology Corporation, Santa Barbara, CA). Measurements were carried out at a single scattering angle $\theta=90^{\circ}$, temperature $T=300$ K, and incident vacuum laser wavelength $\lambda_{0}=658$ nm. The instrument’s digital correlator (512 channels) was employed to compute $g^{(2)}(\tau)$, the normalized second order autocorrelation function (ACF) of the scattered light intensity:

[TABLE]

The second order ACF $g^{(2)}(\tau)$ can be related to the normalized first order electric field ACF $g^{(1)}(\tau)$ throughGoldin (2002):

[TABLE]

with the experimental noise $\xi(\tau)$. The field ACF $g^{(1)}(\tau)$ of a non-uniform (polydisperse) solution can theoretically be describedØgendal (2016) as an integral of exponential decays corresponding to a distribution of different-sized scatterers:

[TABLE]

with the translational diffusion coefficient $D_{t}$, the magnitude of the scattering vector $q$, the solution refraction index $n$, the vacuum wavelength of the incident light $\lambda_{0}$, and the scattering angle $\theta$.

The DYNALS regularization method Goldin (2002) implemented in the DYNAMICS software (Wyatt Technology Corporation) was used to find a discretized approximation for $B(\Gamma)$ in (11) whose coefficients represent the intensity fractions of the scattered light attributed to the different scatterers. $B(\Gamma)$, and its discrete approximation, can be expressed in terms of hydrodynamic radii $R_{H}$ through the Stokes-Einstein relation:

[TABLE]

where $\eta_{0}$ denotes the solvent viscosity ($1.019\text{\times}{10}^{-03}$ Pa$\cdot$s).

The mass fraction attributable to every scatterer $j$ was estimated from its intensity contribution $B_{j}(R_{H},\theta)$ (eqns. (11), (12)) according to:

[TABLE]

where $C_{j}$ is the mass concentration of scatterer $j$,$B_{j}(R_{H},\theta)$ is the intensity of the scattered light for molecule $j$ with radius $R_{H}$ at scattering angle $\theta$, $M_{j}(R_{H})$ is the molar mass, and $P(R_{H},\theta)$ is the scattering function (form factor). The Rayleigh-Gans approximation for random coils was chosen to estimate the angular dependence of $P(R_{H},\theta)$, and it was assumed that $M_{j}(R_{H})\propto R_{H}^{2}$ (Wyatt Technology Corporation, Technical Note 2004).

II.2.3 Measurement of the second virial coefficient of the osmotic pressure, $B_{22}$

Static (i.e. time-averaged) light scattering (SLS) intensities recorded from a series of solutions with varying concentration allow for the calculation of $B_{22}$. SLS intensities from the filtered A2 glycan solutions described in the DLS experiment were acquired using the SLS detector of the DynaPro NanoStar instrument. Temperature, scattering angle and laser light vacuum wavelength were kept the same ($\theta=90^{\circ}$, $T=300$ K, $\lambda_{0}=658$ nm).

The time-averaged scattered light intensity $\langle I_{s}\rangle$ at scattering angle $\theta$ relates to the solute mass concentration $C$ and solute molecular mass $M_{W}$ through van Holde (1985); Øgendal (2016):

[TABLE]

with $\langle I_{0}\rangle$ as time-averaged scattered light intensity of the solvent alone, the distance $r$ between detector and scattering volume $V_{s}$, the solvent refractive index $n_{0}$, the vacuum wavelength of the incident light $\lambda_{0}$, the solution refractive index $n$, and the structure factor $S(\theta)$. The size of the refractive index increment $\frac{\mathrm{d}n}{\mathrm{d}C}$ can be assumed constant over the concentration range under examination, and here we used the empirically found value of 0.145 ($\pm$ 0.005) mL/g for polysaccharidesØgendal (2016). The structure factor $S(\theta)$ describes the phase relation of the light scattered from the molecules in $V_{s}$. Ideal solute molecules (no intermolecular forces) scatter with completely uncorrelated phases, and $S(\theta)$ becomes unity. Non-ideal solute molecules influence each other (excluded volume, attractive or repulsive intermolecular forces) thus having a structure factor $S(\theta)$ deviating from one. For non-ideal molecules that are much smaller than the wavelength of the incident light, as is the case for A2 glycan in the setup described, it has been shown Tanford (1961); Øgendal (2016) that

[TABLE]

with $A_{2}$, $A_{3}$, …as second, third, and higher coefficients of the virial expansion of the osmotic pressure $\Pi$ expressed in terms of the solute mass concentration $C$:

[TABLE]

Inserting eq. 15 in eq. 14 gives after rearrangement

[TABLE]

The third and higher virial coefficients can be neglected in the limit of low concentrations, resulting in a linear relationship of $C\to\frac{KC}{R_{\theta}}$ in (17) which yields the molar mass $M_{W}$ from the y-intercept and $A_{2}$ from the slope. The conversion of the mass concentration based virial coefficients $A_{i}$ into the mole concentration based coefficients $B_{2i}$ is achieved by comparing eqs. 2 and 16:

[TABLE]

III RESULTS

III.1 Selection of saccharides, simulation conditions and water model

We based our analysis on a complex biantennary N-glycan, hereafter referred to as A2 glycan (fig. 1A) which represents a structure commonly found as asparagine-linked N-glycan in vertebrate glycoproteinsRoyle et al. (2008); Hua et al. (2013); Gao et al. (2015) and can be thought of as a prototype of protein glycosylation. It consists of a conserved pentasaccharide core that is substituted by two identical trisaccharide antennae, each terminating with sialic acid residues. Sialic acid carries a carboxyl group which is deprotonated and thus negatively charged at physiological pH. A2 glycan is commercially available as a product purified from natural resources. In addition, to test the hypothesis that the aggregation behavior of polysaccharides in MARTINI depends on their size, we included in our analysis four extensively studied saccharides for which experimental values of $B_{22}$ are known: The monosaccharide glucose111Unless otherwise stated monosaccharides are referred to as in their D-pyranose form., the disaccharide sucrose ($\alpha$-glucosyl-(1$\rightarrow$2)-$\beta$-fructofuranose), and the cyclic oligosaccharides $\alpha$- and $\beta$-CD which are six- and seven-membered cyclic oligomers of $\alpha$1$\rightarrow$4-linked glucose.

Since the aggregation propensity of saccharides in aqueous condition depends on the balance of solute-solute, solute-water, and water-water interactions, the choice of the water model is of particular importance. Currently, three water models exist in MARTINI: First, the standard water modelMarrink et al. (2007), in which groups of four H2O molecules are represented by a single, uncharged P4 bead interacting solely through LJ potentials ($\sigma=0.47$ nm, $\epsilon=5.0$ kJ/mol towards other water particles). Second, antifreeze water, a mixture of standard water and antifreeze particles (usually 10%)Marrink et al. (2007) to prevent freezing at 280-300 K. The antifreeze particles behave as water particles except for interaction with the standard water particles ($\sigma=0.57$ nm, $\epsilon=5.6$ kJ/mol). Third, a polarizable water modelYesylevskyy et al. (2010) in which the central, uncharged bead is connected to two partially charged ($\pm{0.46\;e}$) non-LJ interacting beads and the water-water LJ interaction has been reduced compared to the standard water model ($\sigma=0.47$ nm, $\epsilon=4.0$ kJ/mol).

In our simulations we chose a temperature of 300 K to ensure comparability of calculated $B_{22}$ values with experimental results which have typically been acquired at 298 K. At this temperature, we frequently observed freezing of MARTINI simulations of very diluted systems containing standard water as a solvent. This precluded the HEUS approach to calculate the PMF (see Methods) and thus made calculations of $B_{22}$ for strongly interacting particles impossible with the standard water model. Nevertheless the standard water model in MARTINI is far more often used than the polarizable water model, which is why we chose to use antifreeze water in all our MARTINI simulations. The same water model has been employed by Stark et al.Stark et al. (2013) which makes their findings on the aggregation propensity of proteins directly comparable to ours.

III.2 Spurious aggregation of A2 glycan

We hypothesized that a potential discrepancy between MARTINI and experimental aggregation would increase with saccharide size; consequently, we chose the largest glycan in our selection, A2 glycan, for an initial test for intermolecular interaction in MARTINI. To this end, 105 A2 glycan molecules were placed in a cubic (19 nm)3 box filled with MARTINI water and 100 mM NaCl and 10 mM $\mathrm{CaCl_{2}}$. The A2 glycan is expected to be readily soluble at this concentration (25 mM or 55.6 g/L), as e.g. dextrans (branched polymers of mainly $\alpha$1$\rightarrow$6-linked glucose) of medium to high molecular weight easily dissolve in water up to 400 g/LIoan et al. (2000). In the simulation however, we observed a striking aggregation behavior, resulting in all A2 glycan molecules clumping together within a few tens of nanoseconds (fig. 2A). Furthermore, evolution of the system for a total of 1 $\mu$s revealed not a single dissociation event, which would limit A2 glycan solubility to less than one molecule per volume of the simulation box, i.e. $<$ 0.25 mM. As we demonstrate in the experimental results section below, A2 glycan is readily soluble at concentrations up to 4.5 mM, and we therefore suspected that, in line with our hypothesis, the balance of non-bonded forces was severely biased towards promoting saccharide attraction.

Interestingly, we observed very similar aggregation of A2 glycan molecules with the other two water models, i.e. standard water and polarizable waterYesylevskyy et al. (2010) (fig. S1C,D), suggesting that we found a general problem in MARTINI rather than an isolated issue with a particular water model.

To quantify our findings, we calculated the potential of mean force (PMF, see Methods) of A2 glycan using HEUS (Hamiltonian Exchange Umbrella Sampling) and computed $B_{22}$ (fig. 2B).

Corroborating our qualitative findings we found a pronounced well of -36 kJ/mol at a distance of 0.7-1 nm and concomitantly $B_{22}$ converged to $-1.3\text{\times}{10}^{06}$ L/mol. Experimental $B_{22}$ values for other saccharides have all been found positive (cf. table 1), so we sought to measure $B_{22}$ of A2 glycan experimentally to scrutinize the prediction from MARTINI.

III.3 Experimental determination of $B_{22}$ of A2 glycan

$B_{22}$ values have been determined for the smaller saccharides (glucose, cellobiose, sucrose, $\alpha$-, and $\beta$-CD)Stigter (1960); Herrington et al. (1983); Dagade et al. (2004); Terdale et al. (2006) and some high molecular weight dextrans Ioan et al. (2000); Schaink and Smit (2007), but to our knowledge not for A2 glycan or any similar structure. To elucidate the thermodynamic properties of A2 glycan in aqueous solution with only low amounts of available A2 glycan (1 mg), we performed dynamic and static light scattering experiments which can be conducted with sample volumes as small as 1.5 $\mu$L.

A2 glycan was dissolved in 100 mM NaCl to mimic physiological conditions. As a prerequisite to $B_{22}$ measurements, we determined the dispersity of the A2 glycan solution by dynamic light scattering. The hydrodynamic radius of the main species of A2 glycan solutions (4-10 g/L) was found to be 1.2-1.3 nm, which is in very good agreement with an A2 glycan monomer (fig. 3A, B). Only tiny amounts of larger aggregates were observed (fig. 3A), suggesting that A2 glycan dissolved in water into an essentially uniform (monodisperse) solution of monomers. Notably, A2 glycan dissolved entirely up to a concentration of 10 g/L (4.5 mM), well above the upper solubility limit of 0.25 mM suggested by the MARTINI simulation.

Static light scattering analysis of A2 glycan solutions of different concentrations ($C=4..10$ g/L in 100 mM NaCl) yielded a linear relationship (fig. 3C) between $KC/R_{\theta}$ (see Methods for details) and $C$. From the y-intercept of the linear fit in fig. 3C a molecular mass $M_{W}$ of $4500\pm{700}$ g/mol (95% confidence interval (CI)) was obtained which is around twice the theoretical $M_{W}$ of a monomer (2224 g/mol). From a chemical perspective it seems unlikely that A2 glycan spontaneously dimerizes in aqueous solution, which is why we attribute this discrepancy to systematic errors in $C$. The slope of the fit directly relates to $B_{22}$ of A2 glycan and was found to be ${46}\pm{11}$ L/mol (95% CI). It is noteworthy that this clearly positive value indicates a net repulsive interaction between A2 glycan molecules in solution which is strongly opposing the behavior of A2 glycan in the MARTINI simulation. Table 1 and fig. 3D show that the measured $B_{22}$ value for A2 glycan fits well into the series of known $B_{22}$ values for other saccharides. The $B_{22}$ of A2 glycan is found to be more positive than would be expected for an uncharged saccharide of the same mass, which can be explained by an additional repulsive interaction due to the two negative charges of A2 glycan. Notably, the physiological salt concentration employed in our MD simulations had negligible effect on $B_{22}$ (fig. S1).

Taken together, we found our hypothesis confirmed that at least for A2 glycan the non-bonded force balance was biased towards promoting aggregation in MARTINI. We next investigated ways to adjust the MARTINI force field parameters to make it better reproduce experimental $B_{22}$ values.

III.4 Ways to improve the MARTINI force field

Microscopically, $B_{22}$ reflects the relative strengths of solute-solute, solute-solvent, and solvent-solvent non-bonded interactions, which are modeled by a Coulomb and an LJ interaction potential in MARTINI. The overestimation of the aggregation propensity occurs irrespective of electrical charge (all beads are uncharged except for two negatively charged beads in A2 glycan), necessitating adjustments to the LJ interaction. Specifically, we sought a modification that would retain central force field properties, i.e. the mapping rules and the partitioning behavior between water and apolar solvents, to circumvent a complete reparametrization of MARTINI. Therefore the water-water and water-solute interactions were kept unchanged, leaving $\sigma_{ij}$ of the saccharide beads and $\epsilon_{ij}$ of the saccharide-saccharide interactions, both of which are expected to influence $B_{22}$, as tunable parameters. However, tests of scaling $\sigma_{ij}$ for all interactions involving saccharide beads either had insufficient effect on $B_{22}$ or, due to the potential shift of the LJ potentials, impacted $\epsilon_{ij}$ as well (fig. S3). We therefore left $\sigma_{ij}$ unmodified and instead analyzed the effect of a reduction of the interaction strength, $\epsilon_{ij}$, of saccharide-saccharide interactions on $B_{22}$.

It has to be emphasized that this approach is only possible because solute-solute interactions can be expected to have a minor effect on partitioning coefficients, which depend mostly on solute-solvent and solvent-solvent interactions and were a key feature in the parametrization of the MARTINI force fieldMarrink et al. (2007). It does not imply that these interactions are the only cause for an imbalance in the non-bonded interactions. In fact, the water model in MARTINI is somewhat problematic and a known imbalance exists between solute-water and water-water interactions (see discussion below). Other approaches like improving the water model or changing water-solute interactions are possible routes to establish a better force balance. However, this would essentially mean a complete reparametrization of the MARTINI force field which is beyond the scope of this study.

III.5 Scaling of solute-solute interactions in MARTINI: A2 glycan and glucose

In order to systematically investigate the effect of reducing $\epsilon$, we defined a scaling parameter $\gamma$:

[TABLE]

Thus, $\epsilon_{ij,scaled}$ changes linearly in the interval [0,1] from 2 kJ/mol, the weakest LJ interaction level in MARTINI, to its original value $\epsilon_{ij,original}$. This approach follows the work of Stark et al. who defined an analogous scaling factor $\alpha$ to modify protein-protein interaction levelsStark et al. (2013). Figure 4 shows the effect of varying $\gamma$ between 0 and 1 on the PMF and $B_{22}$ for glucose and A2 glycan. In the case of glucose, unscaled MARTINI ($\gamma=1$) yields a PMF with a potential well of -2.5 kJ/mol with a concomitant $B_{22}$ value of -1.1 L/mol. The implied net attractive intermolecular interaction disagrees with the experimental prediction of a weak repulsion ($B_{22}=0.117\;\mathrm{L/mol}$Stigter (1960)). The depth of the potential well in the PMF flattens with decreasing $\gamma$ (fig. 4A) and $B_{22}$ reaches good agreement with experiment at $\gamma=0.5$ (fig. 4B). Similarly, the deep potential well of the PMF of A2 glycan in unscaled MARTINI (-35 kJ/mol, fig. 4C) flattens quickly with decreasing $\gamma$, shifting $B_{22}$ by six orders of magnitude from $-1.3\text{\times}{10}^{06}$ L/mol to -3 L/mol at $\gamma=0.5$. However, the experimental $B_{22}$ value of 46 L/mol is never reached; even with lowest LJ interaction potentials ($\gamma=0$, $\epsilon_{ij}=2\;\mathrm{kJ/mol}$) $B_{22}$ reaches not more than $\approx$ 7 L/mol (fig. 4D). The value of $\gamma=0.5$ thus constitutes a compromise between reproduction of physical $B_{22}$ values and minimization of the changes to the original MARTINI force field parameters.

Again, we wish to emphasize that this is strictly valid for the antifreeze water model only. Qualitatively, the other two water models (standard water without antifreeze particles and polarizable water) seem to overestimate the aggregation propensity of A2 glycan, too (fig. S1C, D), however the optimal value of the scaling factor $\gamma$ will most likely be different. Quantitative tests with glucose and A2 glycan suggest for example that a scaling factor of $\gamma=0.8$ might suffice for standard water without antifreeze particles (fig. S1A, B).

III.6 Extension to sucrose and cyclodextrins

To determine what effect a reduced LJ interaction strength would have for intermediate sized saccharides, we compared unmodified MARTINI simulations for sucrose, $\alpha$-, and $\beta$-CD (at 25 mM) with their scaled counterparts ($\gamma$=0.5). In unscaled MARTINI nearly all CD molecules and most of the sucrose molecules formed clusters within the first hundreds of nanoseconds of simulation time (fig. 5). The discrepancy with the aqueous solubility for these saccharides (sucrose: 1.97 M, $\alpha$-CD: 121 mMPharr et al. (1989), $\beta$-CD: 16.3 mMPharr et al. (1989); 25 *∘*C) indicates already a clear overestimation of the aggregation propensity. The PMF well depth in unscaled MARTINI varied from -6 to -55 kJ/mol (fig. 6A), whereas the energy of thermal fluctuations (RT) at 300 K equals to 2.5 kJ/mol. Clearly, the most probable state for an assembly of such molecules is strongly bound. Consequently, the corresponding $B_{22}$ values, even though endowed with significant uncertainties, all point to a very strong aggregation propensity (fig. 6B). Conversely, scaling down the LJ interaction strength with $\gamma$=0.5 resulted in fully disperse solutions for all saccharides over a time course of up to 10 $\mu$s (fig. 5). The derived and experimental $B_{22}$ values have been gathered in table 2.

Clearly, unmodified MARTINI overestimates the aggregation propensity in all studied cases. The discrepancies generally escalate with growing $M_{w}$, corroborating our hypothesis about their size dependence. Interestingly, whereas experimental values of $B_{22}$ of saccharides indicate a direct correlation between size and increased intermolecular repulsion, MARTINI simulations point to a strong inverse relationship, fostering the idea that a small overestimation in individual $\epsilon_{ij}$ accumulates in large molecules, leading to the observed net attractive forces. Scaling of LJ interaction strengths led to a very good agreement of $B_{22}$ with corresponding experimental values for glucose and sucrose. For $\alpha$- and $\beta$-CD, agreement is still acceptable although $B_{22}$ values remained slightly negative, pointing to a weak tendency to aggregate. For A2 glycan, despite missing the experimental $B_{22}$ value, an improvement over orders of magnitude has been achieved. Notably, the obtained $B_{22}$ values for the CDs and A2 glycan seem to be identical, pointing to a possible size-independent upper limit of the proposed correction (see below for discussion).

The limitations notwithstanding, even if a precise match between predicted and experimental $B_{22}$ could not be met, the aggregation behavior between $\gamma=1.0$ and $\gamma=0.5$ is critically different. Shallow PMF well depths account for an observed lack of aggregate formation, allowing for more realistic simulation of saccharides and/or glycosylated macromolecules. .

III.7 Validation of the proposed modifications

The density of a solution directly relates to the partial volume of the solute which, on a molecular scale, can be interpreted as the volume increase upon addition of a single solute molecule to a given solution. In a MD simulation, solution density depends on the interplay of solute-solute, solute-solvent and solvent-solvent non-bonded interactions. To test if modified solute-solute interactions compromised the non-bonded force balance, we calculated the solution density for glucose, sucrose, $\alpha$- and $\beta$-CD in water over a range of concentrations and compared it with available experimental data. To this end, the system volume averaged over the last five nanoseconds of 15 ns trajectories was used to calculate the solution density which was plotted against the solute mole fraction (fig. 7). Comparison with experimental dataPaduano et al. (1990); Cerdeiriña et al. (1997); Fucaloro et al. (2007); Santos et al. (2016) shows that MARTINI tends to underestimate the density of saccharide solutions by up to 10% at high solute concentrations (fig. 7A, B). This finding is in line with previous MARTINI simulations for aqueous glucose solutionsLópez et al. (2009). At low concentrations (solute mole fraction $<$ 0.4%) this effect is overcompensated by a $\approx$ 1% overestimation of the density of pure water in MARTINI (fig. 7C, D). The use of antifreeze water leads to an additional systematic underestimation of the density by about 10%, as was stated beforeMarrink et al. (2007). Importantly, the densities obtained from simulations with scaled $\epsilon_{ij}$ ($\gamma=0.5$) closely followed the densities of the standard MARTINI model (deviations $<$ 1.5%), indicating that the proposed force field modification does not disturb the overall non-bonded force balance.

We furthermore tested the effect of the proposed modification in organic solvents by studying systems of A2 glycan and glucose in the apolar solvent hexadecane. The expected low solubility of saccharides under such conditionsMazzobre et al. (2005) was captured not only by the original but also by the scaled MARTINI force field (fig. S4), suggesting the validity of the proposed scaling also in organic solvents.

IV DISCUSSION

MARTINI strongly overestimates the aggregation propensity of saccharides

Every force field is an approximation and therefore its application is limited to cases it has been designed for. For MARTINI this is certainly the simulation of biological macromolecules (proteinsMonticelli et al. (2008), carbohydratesLópez et al. (2009), and DNAUusitalo et al. (2015)). In contrast to AA force fields that are mainly being used (and have been parametrized) for simulations of individual macromolecules, the advantage of a CG force field such as MARTINI lies in the possibility of simulating ensembles of macromolecules on microsecond time scales. This is exemplified by the simulation of a lipid membrane patch containing several membrane protein moleculesIngólfsson et al. (2016) which holds promise for future simulations at subcellular scale. An important aspect in these mesoscopic systems will be how macromolecules interact and thus the accurate representation of non-bonded forces between macromolecules is of crucial importance.

We found that the MARTINI force field strongly overestimates the aggregation propensity of saccharides in aqueous solution. We analyzed this quantitatively for a single water model (antifreeze water) and qualitative tests suggested similar problems for the other two water models (standard water without antifreeze particles and polarizable water; fig. S1C, D). The deviations between experimental and predicted $B_{22}$ values rapidly increase with saccharide size, suggesting that an imbalance in the parametrization of the LJ interaction potentials of the MARTINI beads amplifies as the molecules grow bigger. Hence it is perhaps not surprising that the propensity to aggregate of either proteinsStark et al. (2013) or saccharides, as shown in this study, is strongly overestimated in MARTINI. As a consequence recent MARTINI simulations involving saccharides suggested a tendency to aggregate that most likely does not reflect physical reality: Kapla et al.Kapla et al. (2016) found that trehalose, a disaccharide similar to sucrose, forms clusters in MARTINI polarizable water. The small positive experimental $B_{22}$ of trehalose (table 1), however, indicates net weak repulsion of trehalose molecules so that clustering at the reported concentrations (0.22 kg/kg water; solubility limit: 1.36 kg/kg water at 310 KLammert et al. (1998)) seems very unlikely to represent true physical behavior. Ma et al.Ma et al. (2015) simulated a model of the bacterial outer membrane in MARTINI containing lipopolysaccharide (LPS) molecules, whose polysaccharide components extend into the extracellular space. Similar to our observations with unscaled MARTINI, LPS molecules clustered and collapsed onto the membrane surface. Again, this disagrees with experiments describing LPS as a polymer brush extending tens to a few hundred nanometers into the extracellular spaceIvanov et al. (2011), suggesting instead an overestimation of the aggregation propensity of saccharides and/or lipids.

Validity of the proposed correction, other approaches

The proposed scaling of the solute-solute interaction is an ad hoc remedy to correct the imbalance of non-bonded forces in the MARTINI force field in combination with a particular water model (antifreeze water). It suffices to abrogate unrealistic aggregation behavior of saccharides in aqueous solution, retains the original solvent partitioning behavior, and is compatible with earlier findings for proteinsStark et al. (2013). The underlying cause for the non-bonded force imbalance, however, cannot be limited to the strength of solute-solute interactions alone as is evident by the failure to quantitatively reproduce $B_{22}$ values for larger saccharides. Two extreme tests illustrate the general limit of scaling solute-solute interactions: Neither doubling $\sigma_{solute-solute}$ (which, given the cutoff value and the potential shift, reduces $\epsilon_{solute-solute}$ to almost zero), nor the conversion of all beads in A2 glycan to special S beads (scaling $\epsilon_{solute-solute}$ to 75% and $\sigma_{solute-solute}$ to 0.43 nm, but violating the 4:1 mapping rule in MARTINIMarrink et al. (2007)), lifted $B_{22}$ beyond 15 L/mol (figs. S3 and S5). These findings are indicative of an existing imbalance between solute-water and water-water interactions. It is known from calculations of hydration free energies that the ratio of solute-water/water-water interaction strengths is too low with MARTINI standard waterMarrink et al. (2007). The ratio is even lower with antifreeze water whose average water-water interaction strength is slightly higher. It is not enough though to simply change the water model to standard water or polarizable water as shown by the rapid formation of saccharide aggregates also in these systems, albeit the magnitude of the imbalance is smaller (fig. S1C, D). Moreover, the water model in MARTINI is known to underestimate the water-water interaction strength in the liquid phaseMarrink et al. (2007). Hence, both solute-water and solute-solute interactions need to be reparametrized, possibly in combination with an improved water model, to reach experimental $B_{22}$ values and thus more realistic molecular interaction behavior. Also, we recommend to include not only single beads as chemical building blocks in the parametrization, but biological reference macromolecules (proteins, DNA, saccharides) to avoid an amplification of small errors as has become evident in this study.

Non-bonded interactions in AA force fields and conclusion

CG simulations, including MARTINI, are often verified based on agreement of selected observables with an atomistic approach and the same could be done, in principle, for solute-solute interactions. Recently, however, overestimation of protein aggregation propensity has been reported in practically all modern AA force fields Best et al. (2014); Henriques et al. (2015); Petrov and Zagrovic (2014); Abriata and Dal Peraro (2015); Miller et al. (2016). It has been argued that this is at least partially due to the way AA force fields are fine-tuned, i.e. to maintain protein native structure, leading to too strong protein-protein non-bonded interactions, and suggest scaling of solute-solute/solute-solvent interactions or partial charges adjustment as a remedy. In many of these cases this nonphysical behavior could also stem from the use of certain water models (SPC, TIP3P), which have been shownBest et al. (2014); Miller et al. (2016) to promote solute aggregation. Parenthetically, atomistic simulations of A2 glycan molecules with the GLYCAM06j force field and SPC water (for details see Supplementary Information and consult fig. S2) resulted in aggregating behavior very similar to unscaled MARTINI simulations, suggesting similar imperfections for saccharides and certain water models.

Indeed, the Grafmüller group found nonphysical aggregation of monosaccharides (glucose, mannose, xylose) Sauter and Grafmüller (2015) for GLYCAM06h with TIP3P waterJorgensen et al. (1983) which could largely be corrected by use of the TIP5P water modelMahoney and Jorgensen (2000). The Elcock groupLay et al. (2016) confirmed the finding for GLYCAM06/TIP3P with glucose and sucrose and found a very similar deficit for CHARMM36-TIP3P. They achieved correction by a substantial reduction of $\epsilon$ in LJ interactions between C-O and C-C atoms in the GLYCAM06 force field. It will be interesting to see if the newly developed TIP4P-D waterPiana et al. (2015), demonstrated to improve the compaction of disordered proteins, will also aid saccharide simulations.

We conclude that scaling saccharide-saccharide interactions provides an ad hoc solution to remedy nonphysical aggregation behavior in MARTINI 2.x. Due to the existing non-bonded force imbalance a reparametrization of the entire force field, with particular regard to the water model, will be necessary to facilitate quantitative predictions of aggregation propensities. Otherwise, bold simulation attempts of systems at micron scaleFeig et al. (2015); Yu et al. (2016) run the risk of producing drastically misleading results. Furthermore, macroscopic observables of macromolecular systems, like $B_{22}$, solubility or even mechanical properties of molecules, need to be taken into account in the reparametrization of MARTINI to ensure its suitability for MD simulations of ensembles of macromolecules.

Acknowledgment

The authors thank Antje Potthast, Marek Cieplak, Tomasz Włodarski, and Damien Thompson for fruitful discussions and the IST Austria Scientific Computing Facility for support. P.S.S. was supported by research fellowship 2811/1-1 from the German Research Foundation (DFG) and M.S. was supported by EMBO Long Term Fellowship ALTF 187-2013 and a grant no GC65-32 from the Interdisciplinary Centre for Mathematical and Computational Modelling (ICM), University of Warsaw, Poland.

Supplementary Tables

acf. fig. 1 in the main text

bbead names prefixed with ”B” follow López et al.

cbead types are according to Marrink et al.; prefix ”G” means LJ interaction with other ”G” beads is scaled according to (19). See main text for details.

aIn the Table S2 bond spring constants equal $\infty$ denote that the bond was turned into a constraint. Non-bonded interactions were excluded between selected beads to achieve good agreement of some of the angle distributions with AA simulations.

Supplementary Figures