Enhanced Conformational Sampling of N-Glycans in Solution with

Article pubs.acs.org/JCTC

Enhanced Conformational Sampling of N‑Glycans in Solution with Replica State Exchange Metadynamics Raimondas Galvelis,† Suyong Re,†,§ and Yuji Sugita*,†,#,¶,§ †

RIKEN Theoretical Molecular Science Laboratory, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan RIKEN iTHES, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan ¶ RIKEN Advanced Institute for Computational Science, Integrated Innovation Building 7F, 6-7-1 minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan § RIKEN Quantitative Biology Center, Integrated Innovation Building 7F, 6-7-1 minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan

Downloaded via KAOHSIUNG MEDICAL UNIV on September 3, 2018 at 07:35:12 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

#

S Supporting Information *

ABSTRACT: Molecular dynamics (MD) simulation of a Nglycan in solution is challenging because of high-energy barriers of the glycosidic linkages, functional group rotational barriers, and numerous intra- and intermolecular hydrogen bonds. In this study, we apply different enhanced conformational sampling approaches, namely, metadynamics (MTD), the replicaexchange MD (REMD), and the recently proposed replica state exchange MTD (RSE-MTD), to a N-glycan in solution and compare the conformational sampling efficiencies of the approaches. MTD helps to cross the high-energy barrier along the ω angle by utilizing a bias potential, but it cannot enhance sampling of the other degrees of freedom. REMD ensures moderate-energy barrier crossings by exchanging temperatures between replicas, while it hardly crosses the barriers along ω. In contrast, RSE-MTD succeeds to cross the high-energy barrier along ω as well as to enhance sampling of the other degrees of freedom. We tested two RSE-MTD schemes: in one scheme, 64 replicas were simulated with the bias potential along ω at different temperatures, while simulations of four replicas were performed with the bias potentials for different CVs at 300 K. In both schemes, one unbiased replica at 300 K was included to compute conformational properties of the glycan. The conformational sampling of the former is better than the other enhanced sampling methods, while the latter shows reasonable performance without spending large computational resources. The latter scheme is likely to be useful when a N-glycan-attached protein is simulated.

1. INTRODUCTION

These chemical features give slow interconversion of N-glycan conformations with different rotamer states. Molecular dynamics (MD) simulations are commonly used to examine conformational fluctuation or dynamics of biomolecules.8 However, slow biological processes with the time-scale of milliseconds or longer cannot be easily simulated with brute-force all-atom MD methods.9 One way to overcome the difficulty is to use enhanced conformational sampling approaches.10,11 The replica-exchange MD (REMD) was already applied to sample multiple conformers of a N-glycan in solution or gas phase.12−14 Re et al. simulated multiple conformers of a biantennary complex-type N-glycan (Bi9) in solution, such as two extended, half- and tight-backfolded conformations.13 Nishima et al. examined the conformations of other N-glycans (BiF10, BiB10, and BiBF11) by REMD to study the effect of bisecting GlcNAc and core fucosylation on

Glycans consist of rigid saccharide units connected together through flexible glycosidic linkages, enabling a variety of conformers in solution.1 In contrast to the classical “lock-andkey” mechanisms for ordinary ligand binding to proteins, glycans provide multiple “keys” and one of them is used for specific glycan-receptor recognition.2−5 The most abundant oligosaccharide found in cellular proteins is a N-glycan, which covalently bonds to a protein via the nitrogen atom of asparagine (Asn) side chain (Figure 1a and Figure S1). The 1− 6 glycosidic linkage is characterized by three dihedral angles (φ, ψ, and ω; Figure 1b), where the ω angle gives the three distinct rotamer states: gauche−gauche (gg), gauche−trans (gt), and trans−gauche (tg). The available NMR data suggest that typical N-glycans have mainly gg and gt conformations in solution.6 It has been estimated that these rotamer states are separated by significant energy barriers (∼7 kcal/mol).7 In addition, the multiple hydroxyl groups form various intra- and intermolecular hydrogen-bond networks stabilizing distinct conformations. © 2017 American Chemical Society

Received: January 25, 2017 Published: April 11, 2017 1934

DOI: 10.1021/acs.jctc.7b00079 J. Chem. Theory Comput. 2017, 13, 1934−1942

Article

Journal of Chemical Theory and Computation

and Figure S1) contains a branched penta-saccharide core s e q ue nc e , M anα 1 - 6 (M a n α1 - 3 ) M a n β 1-4 G lcN Ac β 14GlcNAcβ1-Asn-X-Ser/Thr. This form is shared by all Nglycans. The conformational space of N-glycans is, in general, determined by the flexibility of glycosidic linkages. In Bi9, the α1-6 linkage, which is characterized by three dihedral angles φ, ψ, and ω, largely determines the overall conformational space (Figure 1b). We constructed two initial configurations corresponding to the gg and gt conformations using GLYCAM Web Carbohydrate Builder23 and solvated each with 3437 and 3494 water molecules respectively using VMD.24 Bi9 consists of eight linkages connecting nine saccharides. To characterize its global motions, the eight pseudodihedral angles are defined: χ(12), 2O5−2C1−1C4−1C3; χ(23), 3O5−3C1−2C4−2C3; χ(34), 4O5−4C1−3C3−3C2; χ(45), 5O5−5C1−4C2−4C1; χ(56), 6O5−6C1−5C4−5C3; χ(34′), 4′O5−4′C1−3C6−3C5; χ(4′5′), 5′O5−5′C1−4′C2−4′C1; and χ(5′6′), 6′O5−6′C1− 5′C4−5′C3 (Figure S1). In addition to the ω angle, these pseudodihedral angles are used as collective variables and/or for the trajectory analysis. All the simulations were carried out using the GENESIS software package.25 The CHARMM36 additive force field for carbohydrates26−28 and the TIP3P model for water29 were used. The particle mesh Ewald (PME) method30 was used to evaluate long-range electrostatic interactions. Nonbonded interactions were truncated at 10.0 Å with an atom-based switching function that was effective at 8.0 Å. The pairlist for nonbonded interactions was truncated at 11.5 Å and updated every 10 steps. All bonds involving hydrogen atoms in Bi9 were constrained using SHAKE31 and water molecules were treated as rigid with SETTLE.32 Langevin thermostat33 with a damping coefficient of 5 ps−1 was used to simulate the NVT ensemble. Leapfrog integrator was used with an integration time-step of 2 fs. All rings in Bi9 were restrained to prevent inversions by adding the flat-bottom harmonic potentials to the six dihedral angles (O5−C1−C2−C3, C1−C2−C3−C4, C2−C3−C4−C5, C3− C4−C5−O5, C4−C5−O5−C1, and C5−O5−C1−C2):

Figure 1. (a) The topological structure of a N-glycan, Bi9, with the branch names. (b) The central part of Bi9 with the φ, ψ, and ω angles, and the three pseudodihedral angles, χ(23), χ(34), and χ(34′).

the biantennary complex-type N-glycans in solution.12 Compared to the conventional MD, these studies extended the conformational space of N-glycans, but the free-energy landscapes of N-glycans in solution were not perfectly converged within the simulation times. An alternative generalized-ensemble approach is the Hamiltonian replicaexchange molecular dynamics (H-REMD) (or the replicaexchange umbrella sampling (REUS)).15,16 Patel et al. used a set of two-dimensional grid-based bias potentials in H-REMD to improve the rotamer sampling of the glycosidic linkages.17,18 However, the separate umbrella sampling (US) simulations of model disaccharides are required for each glycosidic linkage to parametrize the bias potentials before the production simulations. In this study, we examine conformational sampling efficiency of several enhanced sampling methods, namely, metadynamics (MTD),19,20 temperature-REMD,21 and the recently proposed replica state exchange metadynamics (RSE-MTD)22 for Nglycan simulations in solution. Since each method enhances sampling of a conformational space in a different way, it is meaningful to compare their features in the case of N-glycans. In the next section, we introduce a target N-glycan system and then summarize the theories of MTD, REMD, and RSE-MTD, respectively. The simulation results are analyzed in terms of the time-series of collective variables (CVs), the conformational entropy, the free-energy landscapes, and structure clustering of internal coordinates, and so on. We propose two RSE-MTD schemes: one uses 64 replicas with the bias potentials for the same CV at different temperatures and another reduces the number of replicas to four with three bias potentials for different CVs at the room temperature. The former shows the best sampling efficiency, while the latter gives a reasonable performance and could be used in simulations of larger Nglycan-attached proteins.

⎧1 2 ⎪ k(ϕ − ϕlb) ϕ < ϕlb ⎪2 ϕlb ≤ ϕ≤ub Vrest(ϕ) = ⎨ 0 ⎪ 1 2 ϕ < ϕ ⎪ ub ⎪ k(ϕ − ϕub) ⎩2

(1)

where φlb = 30°, φub = 90°, and k = 20.0 (kcal/mol)/deg2. The two initial Bi9 conformations, gg and gt, with explicit water molecules were minimized and equilibrated with the two MD simulations: 1 ns NPT dynamics at 300 K and 1 atm, and another 1 ns NVT dynamics at 300 K before the production runs. We use the notations of s(gg) and s(gt) to specify different simulations starting from gg or gt conformations. Additionally, two independent conventional MD simulations (1 μs), s(gg) and s(gt), were performed for comparison. Four graphic processing units (GPUs) were used for each simulation.34 To further speed up the simulation, the multiple time step integrator (r-RESPA, reference system propagator algorithm)35 was used with 2.5 fs time step for the short-range interactions and 5 fs time step for the long-range interactions. 2.2. MTD, REMD, and RSE-MTD for Enhancing Conformational Samplings. 2.2.1. Metadynamics (MTD). In MTD, history-dependent bias potential, Vbias, is added to the system potential V0.

2. METHODS 2.1. A Target Molecule and Simulation Methods. A biantennary complex-type N-glycan, named as Bi9, (Figure 1a 1935


Article


Table 1. Replica Conditions (Temperature and Bias) for the MTD, REMD, RSE-MTD(PT), and RSE-MTD(BE) Methods MTD

a

REMD

RSE-MTD(PT)

RSE-MTD(BE)

replica condition

T [K]

CVa

T [K]

CVa

T [K]

CVa

T [K]

CVa

1 2 3 4 63 64

300.00 300.00

none ω

300.00 302.49 305.00 307.52 493.04 496.90

none none none none none none

300.00 300.00 302.49 305.00 489.20 493.04

none ω ω ω ω ω

300.00 300.00 300.00 300.00

none ω χ(34′) χ(23)

CV “none” refers to unbiased replica.

V ( r ⃗) = V0( r ⃗) + Vbias( s ⃗)

of replicas in this study was determined due to the resource limitations. 2.2.3. Replica State Exchange Metadynamics (RSE-MTD). RSE-MTD22 is one of the generalized-ensemble algorithms, which was recently developed to combine the merits of MTD and REMD algorithms. Two similar methods42,43 were already proposed by developers of MTD. In one method, multiple MTDs are performed at different temperatures, and the temperatures and bias potentials are exchanged with the Metropolis-Hasting algorithm such as temperature-REMD. In another method, MTD with the bias potentials for different CVs are carried out and the bias potentials are exchanged with the Metropolis-Hasting algorithm such as Hamiltonian-REMD (or REUS). The former is called PT-MTD,42 while the latter is named as BE-MTD.43 RSE-MTD generalizes these methods in the following two points. First, we allow the replicas with different bias potentials and at multiple temperatures. Second, we use the Suwa-Todo algorithm for exchanging all the replicas at once, by imposing the global balance condition. To emphasize the exchanges of replica states, we named it as the replica state exchange metadynamics (RSE-MTD). We tested the conformational sampling efficiency of this method with alanine polypeptides in gas and water.22 Note that similar approaches combining the different types of enhanced sampling methods have also been reported recently.44 The construction of bias potentials in RSE-MTD was almost the same as that described in the MTD section. We carried out simulations with two RSE-MTD schemes. One simulation (RSE-MTD(PT)) includes 64 replicas that cover a temperature range of 300−493 K. One replica is simulated without any bias, while another 63 replicas have the Vbias for the ω angle in Bi9. The parameters and construction in Vbias are the same as those in MTD. In the RSE-MTD(PT) scheme, the exchanges were performed with the partial Suwa-Todo algorithm, where 64 replicas are divided into 15 overlapping groups: A groups (1−8, 9−16, 17−24, 25−32, 33−40, 41−48, 49−56, and 57−64) and B groups (5−12, 13−20, 21−28, 29−36, 37−44, 45−52, and 53−60). The exchanges within the A and B groups were attempted in alternating order (i.e., A, B, A, B, ...). In another scheme of RSE-MTD (RSE-MTD(BE)), one unbiased and another three replicas with different Vbias for each collective variables (ω, χ(23), and χ(3′4′)) were simulated at 300 K. The exchanges of replicas were done using the Suwa-Todo algorithm. The replica conditions for all the methods used in this study are summarized in Table 1. 2.3. Conformational Analysis of the Simulation Trajectories. The simulation length for each replica in MTD, REMD, RSE-MTD(PT), and RSE-MTD(BE) was 50 ns. We analyzed 20000 snapshots taken from the last 40 ns of the unbiased replica in each simulation. The trajectories are

(2)

The bias potential is a function of a few collective variables (CV), s ⃗ that describe motions related with high-energy barriers on the potential energy surface. During the simulation, every time interval τ, Vbias is updated with a small Gaussian-shaped kernel: ⎛ V ( s ⃗) ⎞ ⎛ | s ⃗ − st⃗ | ⎞ ΔVbias( s ⃗) ⎟ = wexp⎜ − bias ⎟exp⎜ − ⎝ τ ΔT ⎠ ⎝ 2σ 2 ⎠

(3)

where st⃗ is an instantaneous CV value at time t; w and σ represent the height and width of the Gaussian kernels, respectively. The kernel height is rescaled with the welltempered MTD framework.20,36 We selected the ω angle in Bi9 as a CV, because it was identified as crucial for the conformational flexibility of Bi9.13 For constructing Vbias in MTD, we used the parameters in eq 2 and 3: τ = 1.0 ps (500 steps), w = 1.0 kcal/mol, σ = 0.1 rad, ΔT = 2700 K. The system temperature, T, was kept at 300 K. We performed two-replicas MTD simulations: one with the bias and another unbiased for trajectory analysis. These two replicas were exchanged with the Metropolis−Hastings algorithm.37,38 2.2.2. Replica-Exchange Molecular Dynamics (REMD). In REMD,21 copies of the original target systems, replicas, are simulated with canonical MD simulations at different temperatures. The temperatures are exchanged with the Metropolis− Hastings algorithm37,38 between the neighboring replicas or among all the replicas by the Suwa-Todo algorithm39 (this method is often called the replica-permutation method40). Theoretically, the former is designed to impose the detailed balance conditions for the replica exchanges and the latter satisfies only the global balance condition. MD simulations of replicas at high temperatures allow us to sample wider conformational spaces of the target molecules by avoiding from trapping at local energy minimum states, while those at low temperatures search stable conformations near the minimum energy states. By exchanging temperatures, we can realize a random walk in the temperature space, which, in turn, induces a random walk in the potential energy space to overcome energy barriers in the systems with rugged energy surfaces. In the current study, we used 64 replicas, which cover a temperature range of 300−497 K.41 Replica exchanges were attempted every 2.0 ps (1000 steps) based on the MetropolisHasting algorithm with alternating even/odd pairs of replicas. The simulations show 41−55% of the acceptance ratios of replica exchanges (Figure S2). We found random walks in both replica (Figure S3) and temperature (Figure S4) spaces. The highest temperature might not be sufficient to overcome highenergy barriers along with the ω angles in Bi9, but the number 1936


Article

Journal of Chemical Theory and Computation analyzed in terms of the time-series of collective variables (CVs), the conformational entropy,45 the free-energy landscapes, and structure clustering. Regarding CVs, we focus on the ω angle which is enhanced with a bias potential in MTD as well as RSE-MTDs, while the χ(23), and χ(3′4′) angles are selected as examples of unbiased CVs in MTD, REMD, and RSE-MTD(PT). The conformational entropy, S, is defined using eqs 4 and 5: Sk = −

∫0

2π

Pk(q) ln(Pk(q)) dq

(4) Figure 2. Root mean square deviations of the free energy along the ω angle with respect to the final values at 50 ns in the MTD s(gg), MTD s(gt), RSE-MTD(PT) s(gg), and RSE-MTD(PT) s(gt) simulations, respectively.

K

S=

∑ (Sk − S0) K =1

(5)

where Pk(q), Sk, and K are the probability distribution function of the kth degree of freedom, the conformational entropy originating from the kth degree of freedom, and the number of freedoms, respectively. We use all the flexible dihedral angles in Bi9 to computing the entropy. This quantity was previously used in protein folding simulations as an estimator of the simulation convergence.45 The free-energy landscapes of Bi9 are drawn using the spherical coordinates (η and θ angles), which were proposed by Nishima and co-workers.12 η represents the swing motion of the 1-6 branch, whereas θ shows the up−down motion with respect to the xy-plane defined with the core hexagon (see Figure 1 in ref 12). They can separate major conformers of Bi9. In our structure clustering, we use the ω angle and the eight pseudodihedral angles defined in section 2.1: χ(12), χ(23), χ(34), χ(45), χ(56), χ(34′), χ(4′5′), and χ(5′6′) for investigating the linkage flexibility. The density-based spatial clustering of applications with noise (DBSCAN) algorithm46 (as implemented in Scikit-learn47) is used, in which the minimum size of cluster and the spatial resolution (ε) are set to 1% of samples and 15°, respectively. We consider the periodicity of the dihedral and pseudodihedral angles in the structure clustering.

angle in the selected two replicas of the simulations (replicas 1 and 2 for MTD, REMD, and RSE-MTD(PT)). In all the cases, the transitions between the gg and gt conformations were observed frequently, implying good random walks in the space. These transitions rarely happen in the normal MD simulations at room temperatures. Because of the bias potentials, MTD and RSE-MTD(PT) show transitions more often. In contrast, the transitions are less frequent in REMD. For 10−20 ns, the simulation stayed in one of the states (either ω ≈ 60° or −60°) and suddenly jumped into the other stable state. This slow transition happened due to the random walks in temperature space, ensuring all the moderate-energy barrier crossings. Since the frequent structural transitions along the CV (ω) in MTD and RSE-MTD(PT) were expected in theory, we examine the other coordinates, which were not biased in all the three simulations. We selected the pseudodihedral angles, χ(34′) and χ(23), for this purpose and plot their time-series in Figure 4 (see also Figure S6). We used the same replicas for REMD and RSE-MTD(PT) as those in Figure 3. In MTD s(gg) and s(gt), those two angles did not show frequent transitions: χ(34′) stayed at −120° in most of the time and rarely jumped into 120°. In MTD s(gg), no transition was observed from χ(23) = 45°, while only one transition event occurred in MTD s(gt) for the two replicas. In contrast, REMD and RSE-MTD(PT) show frequent transitions in both χ(34′) and χ(23). The frequency of the transitions between REMD and RSE-MTD(PT) are similar to each other. In Figure 5, the conformational entropy45 shows that the MTD s(gg) and s(gt) converge slowly to different values, suggesting that the sampling is nonergodic. In contrast, REMD and RSE-MTD(PT) rapidly converge to the equilibrium values within the error of ∼1kBT just after the first 10 ns. The result suggests that REMD and RSE-MTD(PT) sample approximately the same number of microstates. Those differences result from the methodology of each enhanced sampling method. In MTD and RSE-MTD(PT), simulations were enhanced to lower the critical energy barriers of local coordinates, ω (Figure 3), with the bias potentials. However, as shown in Figure 4, other degrees of freedoms were not sufficiently sampled just by the bias potentials in MTD. REMD and RSE-MTD(PT) utilized the exchanges of temperatures among the replicas. The high temperature simulations in those methods make the fluctuations of global motions of Bi9 larger than those at low temperatures. Moderate-energy barriers along the pseudodihedral angles in Figure 4 could be overcome by the increase of fluctuations. The conformational entropy in Figure 5 shows that these two methods realize more ergodic

3. RESULTS AND DISCUSSION First, we examine the conformational sampling features of MTD, REMD, and RSE-MTD(PT). After investigating the results of RSE-MTD(PT) simulations, three essential CVs are found and used in RSE-MTD(BE). Finally, we discuss the sampling efficiency of all the four simulation methods and compare the results with experimental data. 3.1. Validation of MTD, REMD, and RSE-MTD. First, we verify each of the enhanced conformational sampling methods used in this study. For MTD and RSE-MTD, the convergence of the bias potential is necessary for estimating free-energy changes. In Figure 2, we plot the root-mean-square deviations (RMSDs) of the free-energy changes along the ω angle with respect to the final values at 50 ns. Large deviations were observed in the first 10 ns of MTD s(gg) and s(gt), while RSEMTD(PT) s(gg) and s(gt) show quick convergences toward the final value with small errors of less than 1 kcal/mol. After 10 ns, the RMSDs for RSE-MTD(PT) s(gg) and s(gt) are much smaller than those of MTD. Figure S5 shows the free-energy changes along the ω angle. The free energy is converged within ∼1 kcal/mol in all simulation methods. 3.2. Time Series of ω and Other Degrees of Freedom. Next, we study the conformational fluctuations of Bi9 in different simulations. Figure 3 shows the time series of the ω 1937


Article


Figure 3. Time series of the ω angle in MTD, REMD, RSE-MTD(PT), and RSE-MTD(BE) simulations. s(gg) and s(gt) are different simulations starting from gg and gt conformations of Bi9, respectively.

Figure 4. Time series of the two pseudodihedral angles, χ(34′) and χ(23) in MTD, REMD, RSE-MTD(PT), and RSE-MTD(BE) simulations. s(gg) and s(gt) are different simulations starting from gg and gt conformations of Bi9, respectively.

6 (in the plot of RSE-MTD(PT) s(gg)) and Figure 7, respectively. In MTD, noticeable differences between s(gg) and s(gt) were observed and TBF structures were rarely sampled. In REMD, the sampled spaces are more consistent between s(gg) and s(gt) and all the four conformations were observed as independent basins. However, the populations of local minima are different; that is, the s(gg) favored the gg conformations (θ ≈ 40°), but the s(gt) sampled the gt conformations (θ ≈ −20°) more often. In other words, REMD s(gg) oversampled Ex-a, because of the initial conformation dependency. The RSE-MTD(PT) results look the most converged. The populations of Ex-a and Ex-b seem to be the same regardless of the starting conformations. We further examine the equilibrium solution structures of Bi9 in terms of 9 internal coordinates: the ω angle and 8 pseudodihedral angles. If we neglect the internal conformational changes in each saccharide unit, this should be a complete description of Bi9 structure through all the linkage connectivity. The trajectories of unbiased replicas from the RSE -MTD(PT) s(gg) and s(gt) were merged for this clustering analysis based on the DBSCAN algorithm. DBSCAN detected 5 structures in Figure 7. The first four conformers were already assigned in our previous study,12 such as Ex-a, Ex-b, TBF, and BF, while the last one, Extended C (Ex-c), was found only in this study. In Figure 8, we project each conformer to the space spanning with the η and θ angles (same as the free energy map

Figure 5. Conformational entropy from the MTD, REMD, RSEMTD(PT), and RSE-MTD(BE) simulations.

samplings. In terms of both local and global conformational sampling features, RSE-MTD(PT) seems to be the best candidate for efficient simulations of N-glycans in solution. 3.3. Free-Energy Landscapes and Structural Clustering at Room Temperatures. For an ergodic simulation, the free energy estimate has to be reproducible and independent from the initial conformation. The free energy plot along the spherical coordinates (η and θ angles) has been shown that it can show distributions of the important conformers of Bi9, namely, extended A (Ex-a), extended B (Ex-b), tight backfold (TBF), and backfold (BF).12 Their positions in the free-energy plot and the representative conformations are shown in Figure 1938


Article


Figure 6. Free-energy plots in the spherical coordinates by Nishima et al. (see text for details). Each plot is obtained from the trajectory of MTD, REMD, RSE-MTD(PT), and RSE-MTD(BE). For each method, two independent simulations with different initial conformations were carried out as s(gg) and s(gt). The positions of major conformers found by the clustering analysis were indicated in the plot of RSE-MTD(PT) s(gg).

ω, χ(34′), and χ(23), are necessary to distinguish all the clusters (Figure 8 and Figure S7). The ω angle makes a distinction between the gg (Ex-a and BF) and gt (Ex-b, TBF, and Ex-c) conformers. χ(34′) distinguishes the extended (Ex-a, Ex-b, and Ex-c) with compact (TBF and BF) conformations, while χ(23) distinguishes Ex-a with Ex-c. For quantitative analysis, conformations from unbiased trajectory in each simulation are assigned to the clusters using the nearestneighbor classifier. The estimated cluster populations in Table 2 show the same trends with the free energy plots (Figure 6); that is, the RSE-MTD(PT) simulations are better converged and more consistent than MTD and REMD. 3.4. RSE-MTD(BE) Simulations. The key information from the clustering is that three inter-residue angles, ω, χ(34′), and χ(23), are able to account for the N-glycan flexibility and all the representative conformers of Bi9 in solution. Assuming that the linkage flexibility is the only source of high-energy barriers, we propose RSE-MTD(BE), where the three angles are used as CVs. In addition to one unbiased replica, the four replicas are

Figure 7. Five representative clusters are shown with the top 10 representative structures in each cluster. The core, 1−3 branch, and 1−6 branch are colored in blue, green, and red, respectively. Both top and front views are shown.

in Figure 6). We found that the distributions of Ex-b and Ex-c are overlapped in the space. Projection of the five representative clusters to the nine internal coordinates reveal that just three inter-residue angles,

Figure 8. Projections of the five clusters to (a) (ω, χ(23)), (b) (ω, χ(34′)), and (c) (η, θ) spaces. Red, green, blue, yellow, and magenta dots represent Ex-a, Ex-b, TBF, BF, and Ex-c conformers, respectively. 1939


Article


Table 2. Experimental and Simulated Cluster Populationsa and Scalar 3J-Coupling Constants. Each Row Corresponds to an Independent Simulation with the Specified Method (MTD, REMD, RSE-MTD(PT), or RSE-MTD(BE)) and Initial Condition (s(gg) or s(gt)) gg conformers [%] method/simulation experimental data MTD REMD RSE-MTD(PT) RSE-MTD(BE)

s(gg) s(gt) s(gg) s(gt) s(gg) s(gt) s(gg) s(gt)

3

gt conformers [%]

J-coupling [Hz]

Ex-a

BF

HBF

Ex-b

TBF

Ex-c

J56

J56′

45b 55.9 61.9 56.0 37.9 42.7 46.9 43.1 37.8

10b 4.4 3.9 8.5 4.6 6.9 6.7 10.5 7.8

10b 0.0e 0.0e 0.0e 0.0e 0.0e 0.0e 0.0e 0.0e

35b 32.9 30.0 23.0 36.9 30.3 33.9 32.1 44.2

0d 5.9 1.7 10.8 14.8 16.9 10.6 10.5 10.1

0 0.0 2.1 0.6 4.9 2.2 1.3 3.7 0.1

2.1c 2.9 2.8 2.9 3.0 3.0 2.9 2.9 2.9

5.8c 4.8 4.4 4.6 6.5 5.9 5.6 5.6 6.3

a

DBSCAN ref 46 and 47 classifies some samples as noise, so the cluster populations do not add up to 100%. bThe experimental cluster populations obtained from the PCS data fitting ref 48. cThe experimental NMR scalar 3J-coupling constants taken from ref 49. dEx-c was not taken into account in the PCS data fitting (ref 48). eIn our simulations, HBF was not found as a separate cluster, but rather as a part of BF.

coupling constants (J56 and J56′),49 which are related to the ω angle by the modified Karplus equations.12,50 The J56 values (Table 2) are not sensitive to the sampling protocols, that is, there is just ∼0.2 Hz variations and ∼1 Hz systematic overestimation. Meanwhile, there is a larger variation of the J56′ values in the different simulations. The RSE-PT-MTD simulations show the best agreement (∼0.2 Hz error) with the experiment. 3.6. Comparison with Long Time Conventional MD Simulation. Finally, we compare the enhanced sampling approaches with conventional MD simulations (1 μs). The transitions between the gg and gt conformations occur every few hundreds nanoseconds (three times in s(gg) and six times in s(gt) simulations, Figure S8) due to the high-energy barrier separating the two states. This underscores the importance of the enhanced sampling approach to obtain the converged conformational sampling of N-glycans in solution.

simulated at 300 K. Since the bias potentials could enhance sampling of the CVs, good random walks in the coordinates were observed in Figures 3 and 4. Furthermore, the conformational entropy obtained with RSE-MTD(BE) converges to the same values as REMD and RSE-MTD(PT), suggesting that these three methods sampled almost the same microstates. However, the convergences of the free-energy map (Figure 6) and populations of the representative clusters with RSE-MTD(BE) are slightly worse compared to those with the RSE-MTD(PT) simulation: the differences between s(gg) and s(gt) are larger than those of RSE-MTD(PT). Considering that the total number of replicas in RSE-MTD(BE) is significantly smaller than those in RSE-MTD(PT), this method can be a reasonable choice for simulating larger biological systems, such as singly or multiply glycosilated proteins in solution or in membranes. 3.5. Comparison between Simulations and Experiments. The comparison with experimental data cannot be done without touching force field accuracy and experimental data interpretation issues. Nevertheless, it serves as an indicator how realistic are these simulations. First, we compare the cluster populations with those estimated from the NMR pseudocontact shift (PCS) data.48 In general, there is a reasonably good agreement in terms of the Ex-a, Ex-b, and BF populations (Table 2), for example, in the case of RSEMTD(PT), the estimated errors are below 5%. However, we find a discrepancy in the HBF and TBF populations between the simulated and the experimentally derived results. RSEMTD(PT) gives the gg/gt ratio of 52:48 on average, while the experiment suggests a ratio of 65:35 (Figure 9). Second, the simulation results are compared with the NMR scalar 3J-

4. CONCLUSION AND FUTURE PERSPECTIVES In the current study, we examine the conformational sampling efficiency of enhanced sampling schemes, MTD, REMD, RSEMTD(PT), and RSE-MTD(BE). The simple MTD does not show good performance except for the sampling of the biased CV. Among the four methods, RSE-MTD(PT) shows the best sampling performance in terms of the frequency of structural transitions, the accuracy of free-energy landscapes, and the cluster populations. However, considering the computational cost, RSE-MTD(BE) could be a reasonable choice for larger Nglycan-attached proteins. The penta-saccharide part, which is crucial for overall conformation of Bi9, is shared with the most of N-glycans. Although we simulated only Bi9, RSE-MTD(BE) methods, applying the bias potential to the three inter-residue angles: ω, χ(34′), and χ(23), would be generally useful for all N-glycans. In other systems, including proteins, nucleic acids, or other biomolecules, if we could find essential CVs a proiri, or after short REMD simulations, it is possible to reduce the computational cost in the same way as we propose with RSEMTD(BE).

■

ASSOCIATED CONTENT

S Supporting Information *

Figure 9. Ratio of gg/gt conformations estimated from MTD, REMD, RSE-MTD(PT), and RSE-MTD(BE). The corresponding ratio derived from NMR PCS data (ref 48) is also shown for comparison.

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00079. 1940


Article


■

(11) Mitsutake, A.; Sugita, Y.; Okamoto, Y. Generalized-ensemble algorithms for molecular simulations of biopolymers. Biopolymers 2001, 60 (2), 96−123. (12) Nishima, W.; Miyashita, N.; Yamaguchi, Y.; Sugita, Y.; Re, S. Effect of bisecting GlcNAc and core fucosylation on conformational properties of biantennary complex-type N-glycans in solution. J. Phys. Chem. B 2012, 116 (29), 8504−8512. (13) Re, S.; Miyashita, N.; Yamaguchi, Y.; Sugita, Y. Structural diversity and changes in conformational equilibria of biantennary complex-type N-glycans in water revealed by replica-exchange molecular dynamics simulation. Biophys. J. 2011, 101 (10), L44−L46. (14) Yamaguchi, Y.; Nishima, W.; Re, S.; Sugita, Y. Confident identification of isomeric N-glycan structures by combined ion mobility mass spectrometry and hydrophilic interaction liquid chromatography. Rapid Commun. Mass Spectrom. 2012, 26 (24), 2877−2884. (15) Fukunishi, H.; Watanabe, O.; Takada, S. On the Hamiltonian replica exchange method for efficient sampling of biomolecular systems: Application to protein structure prediction. J. Chem. Phys. 2002, 116 (20), 9058−9067. (16) Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional replicaexchange method for free-energy calculations. J. Chem. Phys. 2000, 113 (15), 6042−6051. (17) Patel, D. S.; Pendrill, R.; Mallajosyula, S. S.; Widmalm, G.; MacKerell, A. D., Jr. Conformational properties of α- or β-(1→6)linked oligosaccharides: Hamiltonian replica exchange MD simulations and NMR experiments. J. Phys. Chem. B 2014, 118 (11), 2851−2871. (18) Yang, M.; MacKerell, A. D., Jr. Conformational sampling of oligosaccharides using Hamiltonian replica exchange with twodimensional dihedral biasing potentials and the weighted histogram analysis method (WHAM). J. Chem. Theory Comput. 2015, 11 (2), 788−799. (19) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12562−12566. (20) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100 (2), 020603. (21) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314 (1−2), 141− 151. (22) Galvelis, R.; Sugita, Y. Replica state exchange metadynamics for improving the convergence of free energy estimates. J. Comput. Chem. 2015, 36 (19), 1446−1455. (23) Woods, R. J., GLYCAM Web; Complex Carbohydrate Research Center, University of Georgia: Athens, GA, 2005−2011. (24) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graphics 1996, 14 (1), 33−38. (25) Jung, J.; Mori, T.; Kobayashi, C.; Matsunaga, Y.; Yoda, T.; Feig, M.; Sugita, Y. GENESIS: a hybrid-parallel and multi-scale molecular dynamics simulator with enhanced sampling algorithms for biomolecular and cellular simulations. Wiley Interdiscip Rev. Comput. Mol. Sci. 2015, 5 (4), 310−323. (26) Guvench, O.; Greene, S. N.; Kamath, G.; Brady, J. W.; Venable, R. M.; Pastor, R. W.; MacKerell, A. D., Jr. Additive empirical force field for hexopyranose monosaccharides. J. Comput. Chem. 2008, 29 (15), 2543−2564. (27) Guvench, O.; Hatcher, E. R.; Venable, R. M.; Pastor, R. W.; MacKerell, A. D., Jr. CHARMM Additive All-Atom Force Field for Glycosidic Linkages between Hexopyranoses. J. Chem. Theory Comput. 2009, 5 (9), 2353−2370. (28) Guvench, O.; Mallajosyula, S. S.; Raman, E. P.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T. J.; Jamison, F. W., II; MacKerell, A. D., Jr. CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling. J. Chem. Theory Comput. 2011, 7 (10), 3162−3180. (29) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926−935.

Structure of a N-glycan Bi9; time series data; other figures supporting the text (PDF)

AUTHOR INFORMATION

Corresponding Author

*Tel.: +81-48-462-1407. Fax: +81-48-467-4532. E-mail: [email protected]. ORCID

Yuji Sugita: 0000-0001-9738-9216 Funding

This research was supported in part by the RIKEN basic science interdisciplinary research projects (dynamic structural biology, molecular systems research, iTHES, and integrated lipidology research), Innovative Drug Discovery Infrastructure through Functional Control of Biomolecular Systems, Priority Issue 1 in Post-K Supercomputer Development (Project ID: hp150270), MEXT Grant-in-Aid for Scientific Research on Innovative Areas Grant No. 26119006 (to Y.S.), and MEXT/ JSPS KAKENHI Grant No. 26220807 (to Y.S.), and Grant No. 16K00415 (to S.R.), Center of innovation Program from Japan Science and Technology Agency, JST (to Y.S.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research used computational resources of the HPCI system provided by the University of Tokyo through the HPCI System Research Project (Project ID: hp140157), the RIKEN Integrated Cluster of Clusters (RICC), and MEXT SPIRE Supercomputational Life Science (SCLS).

■

REFERENCES

(1) Werz, D. B.; Ranzinger, R.; Herget, S.; Adibekian, A.; von der Lieth, C. W.; Seeberger, P. H. Exploring the structural diversity of mammalian carbohydrates (″glycospace″) by statistical databank analysis. ACS Chem. Biol. 2007, 2 (10), 685−691. (2) Gabius, H. J. The magic of the sugar code. Trends Biochem. Sci. 2015, 40 (7), 341. (3) Gabius, H. J.; Andre, S.; Jimenez-Barbero, J.; Romero, A.; Solis, D. From lectin structure to functional glycomics: principles of the sugar code. Trends Biochem. Sci. 2011, 36 (6), 298−313. (4) Hardy, B. J. The glycosidic linkage flexibility and time-scale similarity hypotheses. J. Mol. Struct.: THEOCHEM 1997, 395, 187− 200. (5) Weiss, A. A.; Iyer, S. S. Glycomics aims to interpret the third molecular language of cells. Microbe 2007, 2, 489−497. (6) Wormald, M. R.; Petrescu, A. J.; Pao, Y. L.; Glithero, A.; Elliott, T.; Dwek, R. A. Conformational studies of oligosaccharides and glycopeptides: complementarity of NMR, X-ray crystallography, and molecular modelling. Chem. Rev. 2002, 102 (2), 371−386. (7) Peric-Hassler, L.; Hansen, H. S.; Baron, R.; Hunenberger, P. H. Conformational properties of glucose-based disaccharides investigated using molecular dynamics simulations with local elevation umbrella sampling. Carbohydr. Res. 2010, 345 (12), 1781−1801. (8) Karplus, M.; McCammon, J. A. Molecular dynamics simulations of biomolecules (vol 9, pg 646, 2002). Nat. Struct. Biol. 2002, 9 (10), 788−788. (9) Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H.; Shaw, D. E. Biomolecular simulation: a computational microscope for molecular biology. Annu. Rev. Biophys. 2012, 41, 429−452. (10) Bernardi, R. C.; Melo, M. C.; Schulten, K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850 (5), 872−877. 1941


Article

Journal of Chemical Theory and Computation (30) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103 (19), 8577−8593. (31) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. NumericalIntegration of Cartesian Equations of Motion of a System with Constraints - Molecular-Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23 (3), 327−341. (32) Miyamoto, S.; Kollman, P. A. SETTLE - an Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13 (8), 952−962. (33) Adelman, S. A.; Doll, J. D. Generalized Langevin Equation Approach for Atom-Solid-Surface Scattering - General Formulation for Classical Scattering Off Harmonic Solids. J. Chem. Phys. 1976, 64 (6), 2375−2388. (34) Jung, J.; Naurse, A.; Kobayashi, C.; Sugita, Y. Graphics Processing Unit Acceleration and Parallelization of GENESIS for Large-Scale Molecular Dynamics Simulations. J. Chem. Theory Comput. 2016, 12 (10), 4947−4958. (35) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular-Dynamics. J. Chem. Phys. 1992, 97 (3), 1990− 2001. (36) Dama, J. F.; Parrinello, M.; Voth, G. A. Well-tempered metadynamics converges asymptotically. Phys. Rev. Lett. 2014, 112 (24), 240602. (37) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21 (6), 1087−1092. (38) Hastings, W. K. Monte-Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika 1970, 57 (1), 97−109. (39) Suwa, H.; Todo, S. Markov chain Monte Carlo method without detailed balance. Phys. Rev. Lett. 2010, 105 (12), 120603. (40) Itoh, S. G.; Okumura, H. Replica-Permutation Method with the Suwa-Todo Algorithm beyond the Replica-Exchange Method. J. Chem. Theory Comput. 2013, 9 (1), 570−581. (41) Patriksson, A.; van der Spoel, D. A temperature predictor for parallel tempering simulations. Phys. Chem. Chem. Phys. 2008, 10 (15), 2073−2077. (42) Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M. Free-energy landscape for beta hairpin folding from combined parallel tempering and metadynamics. J. Am. Chem. Soc. 2006, 128 (41), 13435−13441. (43) Piana, S.; Laio, A. A bias-exchange approach to protein folding. J. Phys. Chem. B 2007, 111 (17), 4553−4559. (44) Xie, L.; Shen, L.; Chen, Z. N.; Yang, M. Efficient free energy calculations by combining two complementary tempering sampling methods. J. Chem. Phys. 2017, 146 (2), 024103. (45) Li, D. W.; Bruschweiler, R. In silico relationship between configurational entropy and soft degrees of freedom in proteins and peptides. Phys. Rev. Lett. 2009, 102 (11), 118108. (46) Ester, M.; Kriegel, H.-p.; S, J.; Xu, X. In A density-based algorithm for discovering clusters in large spatial databases with noise, AAAI Press, 1996; pp 226−231. (47) Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; Duchesnay, E. Scikit-learn: Machine Learning in Python. J. Mach. Learn. Res. 2011, 12, 2825−2830. (48) Canales, A.; Mallagaray, A.; Perez-Castells, J.; Boos, I.; Unverzagt, C.; Andre, S.; Gabius, H. J.; Canada, F. J.; JimenezBarbero, J. Breaking Pseudo-Symmetry in Multiantennary Complex NGlycans Using Lanthanide-Binding Tags and NMR Pseudo-Contact Shifts. Angew. Chem., Int. Ed. 2013, 52 (51), 13789−13793. (49) Homans, S. W.; Dwek, R. A.; Boyd, J.; Mahmoudian, M.; Richards, W. G.; Rademacher, T. W. Conformational transitions in Nlinked oligosaccharides. Biochemistry 1986, 25 (20), 6342−6350. (50) Haasnoot, C. A. G.; de Leeuw, F. A. A. M.; Altona, C. The relationship between proton-proton NMR coupling constants and substituent electronegativitiesI. Tetrahedron 1980, 36 (19), 2783− 2792.

1942


Enhanced Conformational Sampling of N-Glycans in Solution with

Recommend Documents