Using the local elevation method to construct optimized umbrella sampling potentials: Calculation of the relative free energies and interconversion barriers of glucopyranose ring conformers in water

The counteracting influence of 2-hydroxypropyl substitution and the presence of a guest molecule on the shape and size of the β-cyclodextrin cavity

The aqueous solubility of β-cyclodextrin (β-CD), a cyclic carbohydrate comprising seven α-D-glucose molecules, is enhanced by 2-hydroxypropyl (2-HP) substitution of the hydroxyl groups at the CD rims. Our thorough analysis of the structural and solvation properties with different degrees of 2-hydroxypropyl substitution on β-CD using molecular dynamics simulations reveals that the solubility is enhanced at the cost of the structural distortion of the CD cyclic structure. Substitution at the secondary rim predominantly enhances the favourable interactions between CD and water by decreasing CD-CD hydrogen bonding and promoting CD-water hydrogen bonding. However, the effect of substitution at the primary rim on the CD-water interactions is minimal; the hydrogen bonds between water and the primary hydroxyl group in native CD merely get replaced by those between water and 2-HP, since the substitution makes the primary hydroxyl oxygen (O6 atom) inaccessible to water. In contrast, substitution at the primary rim maintains the structural integrity of CD, while substitution at the secondary rim results in structural distortion due to the disruption of the intramolecular hydrogen bond belt, even leading to cavity closure. Certain strategic substitutions of the primary hydroxyl groups can help in the reduction of structural distortion, depending upon the degree of substitution at the secondary hydroxyl rim. A detailed inspection of the simulation trajectory revealed that the tilting of glucose units with the primary hydroxyl oxygen (O6) pointing inward is the primary driver for cavity closure. Even though the dynamics of glucose tilting can influence the kinetics of host-guest complex formation, once the guest is well incorporated into the cavity, glucose tilting is inhibited and the cavity opens up as in native β-CD.

Flexible Gaussian Accelerated Molecular Dynamics to Enhance Biological Sampling

Molecular dynamics simulations often struggle to obtain sufficient sampling to study complex molecular events due to high energy barriers separating the minima of interest. Multiple enhanced sampling techniques have been developed and improved over the years to tackle this issue. Gaussian accelerated molecular dynamics (GaMD) is a recently developed enhanced sampling technique that works by adding a biasing potential, lifting the energy landscape up, and decreasing the height of its barriers. GaMD allows one to increase the sampling of events of interest without the need of a priori knowledge of the system or the relevant coordinates. All required acceleration parameters can be obtained from a previous search run. Upon its development, several improvements for the methodology have been proposed, among them selective GaMD in which the boosting potential is selectively applied to the region of interest. There are currently four selective GaMD methods that have shown promising results. However, all of these methods are constrained on the number, location, and scenarios in which this selective boosting potential can be applied to ligands, peptides, or protein-protein interactions. In this work, we showcase a GROMOS implementation of the GaMD methodology with a fully flexible selective GaMD approach that allows the user to define, in a straightforward way, multiple boosting potentials for as many regions as desired. We show and analyze the advantages of this flexible selective approach on two previously used test systems, the alanine dipeptide and the chignolin peptide, and extend these examples to study its applicability and potential to study conformational changes of glycans and glycosylated proteins.

Computational Approaches for Organic Semiconductors: From Chemical and Physical Understanding to Predicting New Materials

While a complete understanding of organic semiconductor (OSC) design principles remains elusive, computational methods─ranging from techniques based in classical and quantum mechanics to more recent data-enabled models─can complement experimental observations and provide deep physicochemical insights into OSC structure-processing-property relationships, offering new capabilities for in silico OSC discovery and design. In this Review, we trace the evolution of these computational methods and their application to OSCs, beginning with early quantum-chemical methods to investigate resonance in benzene and building to recent machine-learning (ML) techniques and their application to ever more sophisticated OSC scientific and engineering challenges. Along the way, we highlight the limitations of the methods and how sophisticated physical and mathematical frameworks have been created to overcome those limitations. We illustrate applications of these methods to a range of specific challenges in OSCs derived from π-conjugated polymers and molecules, including predicting charge-carrier transport, modeling chain conformations and bulk morphology, estimating thermomechanical properties, and describing phonons and thermal transport, to name a few. Through these examples, we demonstrate how advances in computational methods accelerate the deployment of OSCsin wide-ranging technologies, such as organic photovoltaics (OPVs), organic light-emitting diodes (OLEDs), organic thermoelectrics, organic batteries, and organic (bio)sensors. We conclude by providing an outlook for the future development of computational techniques to discover and assess the properties of high-performing OSCs with greater accuracy.

Potential Energy Surfaces Sampled in Cremer-Pople Coordinates and Represented by Common Force Field Functionals for Small Cyclic Molecules

The complex conformations of the cyclic moieties impact the physical and chemical properties of molecules. In this work, we chose 22 molecules of four-, five-, and six-membered rings and performed a thorough conformational sampling using Cremer-Pople coordinates. With consideration of symmetries, we obtained a total of 1504 conformational structures for four-membered, 5576 for five-membered, and 13509 for six-membered rings. All well-known and many less well-known conformers for each molecule were identified. We represented the potential energy surfaces (PESs) by fitting the data to common analytical force field (FF) functional forms. We found that the general features of PESs can be described by the essential FF functional forms; however, the accuracy of representation can be improved remarkably by including the torsion-bond and torsion-angle coupling terms. The best fit yields R-squared (R²) values close to 1.0 and mean absolute errors in energy less than 0.3 kcal/mol.

Multireference Generalization of the Weighted Thermodynamic Perturbation Method

We describe the generalized weighted thermodynamic perturbation (gwTP) method for estimating the free energy surface of an expensive "high-level" potential energy function from the umbrella sampling performed with multiple inexpensive "low-level" reference potentials. The gwTP method is a generalization of the weighted thermodynamic perturbation (wTP) method developed by Li and co-workers [J. Chem. Theory Comput. 2018, 14, 5583-5596] that uses a single "low-level" reference potential. The gwTP method offers new possibilities in model design whereby the sampling generated from several low-level potentials may be combined (e.g., specific reaction parameter models that might have variable accuracy at different stages of a multistep reaction). The gwTP method is especially well suited for use with machine learning potentials (MLPs) that are trained against computationally expensive ab initio quantum mechanical/molecular mechanical (QM/MM) energies and forces using active learning procedures that naturally produce multiple distinct neural network potentials. Simulations can be performed with greater sampling using the fast MLPs and then corrected to the ab initio level using gwTP. The capabilities of the gwTP method are demonstrated by creating reference potentials based on the MNDO/d and DFTB2/MIO semiempirical models supplemented with the "range-corrected deep potential" (DPRc). The DPRc parameters are trained to ab initio QM/MM data, and the potentials are used to calculate the free energy surface of stepwise mechanisms for nonenzymatic RNA 2'-O-transesterification model reactions. The extended sampling made possible by the reference potentials allows one to identify unequilibrated portions of the simulations that are not always evident from the short time scale commonly used with ab initio QM/MM potentials. We show that the reference potential approach can yield more accurate ab initio free energy predictions than the wTP method or what can be reasonably afforded from explicit ab initio QM/MM sampling.

Multifaceted Computational Modeling in Glycoscience

Glycoscience assembles all the scientific disciplines involved in studying various molecules and macromolecules containing carbohydrates and complex glycans. Such an ensemble involves one of the most extensive sets of molecules in quantity and occurrence since they occur in all microorganisms and higher organisms. Once the compositions and sequences of these molecules are established, the determination of their three-dimensional structural and dynamical features is a step toward understanding the molecular basis underlying their properties and functions. The range of the relevant computational methods capable of addressing such issues is anchored by the specificity of stereoelectronic effects from quantum chemistry to mesoscale modeling throughout molecular dynamics and mechanics and coarse-grained and docking calculations. The Review leads the reader through the detailed presentations of the applications of computational modeling. The illustrations cover carbohydrate-carbohydrate interactions, glycolipids, and N- and O-linked glycans, emphasizing their role in SARS-CoV-2. The presentation continues with the structure of polysaccharides in solution and solid-state and lipopolysaccharides in membranes. The full range of protein-carbohydrate interactions is presented, as exemplified by carbohydrate-active enzymes, transporters, lectins, antibodies, and glycosaminoglycan binding proteins. A final section features a list of 150 tools and databases to help address the many issues of structural glycobioinformatics.

Pyranose Ring Puckering Thermodynamics for Glycan Monosaccharides Associated with Vertebrate Proteins

The conformational properties of carbohydrates can contribute to protein structure directly through covalent conjugation in the cases of glycoproteins and proteoglycans and indirectly in the case of transmembrane proteins embedded in glycolipid-containing bilayers. However, there continue to be significant challenges associated with experimental structural biology of such carbohydrate-containing systems. All-atom explicit-solvent molecular dynamics simulations provide a direct atomic resolution view of biomolecular dynamics and thermodynamics, but the accuracy of the results depends on the quality of the force field parametrization used in the simulations. A key determinant of the conformational properties of carbohydrates is ring puckering. Here, we applied extended system adaptive biasing force (eABF) all-atom explicit-solvent molecular dynamics simulations to characterize the ring puckering thermodynamics of the ten common pyranose monosaccharides found in vertebrate biology (as represented by the CHARMM carbohydrate force field). The results, along with those for idose, demonstrate that the CHARMM force field reliably models ring puckering across this diverse set of molecules, including accurately capturing the subtle balance between 4C1 and 1C4 chair conformations in the cases of iduronate and of idose. This suggests the broad applicability of the force field for accurate modeling of carbohydrate-containing vertebrate biomolecules such as glycoproteins, proteoglycans, and glycolipids.

Efficient Alchemical Intermediate States in Free Energy Calculations Using λ-Enveloping Distribution Sampling

Alchemical free energy calculations generally require intermediate states along a coupling parameter λ to establish sufficient phase space overlap for obtaining converged results. Such intermediate states can also be engineered to lower the energy barriers and, consequently, reduce the required sampling time. The recently introduced λ-enveloping distribution sampling (λ-EDS) scheme combines the properties of the minimum variance pathway and the EDS methods to improve sampling and allow for larger steps along the alchemical pathway compared to conventional approaches. This scheme also eliminates the need for soft-core potentials and retains the behavior of conventional λ-intermediate states as a limiting case. In this study, an automated procedure is developed to select the parameters of λ-EDS for optimal performance. The underlying theory is illustrated based on simulations of simple test systems (bond length changes in harmonic oscillators, mutations of dihedral angles, and charge creation in water), as well as on the calculation of the absolute hydration free energies of 12 small organic molecules.

Solvent-scaling as an alternative to coarse-graining in adaptive-resolution simulations: The adaptive solvent-scaling (AdSoS) scheme

A new approach termed Adaptive Solvent-Scaling (AdSoS) is introduced for performing simulations of a solute embedded in a fine-grained (FG) solvent region itself surrounded by a coarse-grained (CG) solvent region, with a continuous FG ↔ CG switching of the solvent resolution across a buffer layer. Instead of relying on a distinct CG solvent model, the AdSoS scheme is based on CG models defined by a dimensional scaling of the FG solvent by a factor s, accompanied by an s-dependent modulation of the atomic masses and interaction parameters. The latter changes are designed to achieve an isomorphism between the dynamics of the FG and CG models, and to preserve the dispersive and dielectric solvation properties of the solvent with respect to a solute at FG resolution. This scaling approach offers a number of advantages compared to traditional coarse-graining: (i) the CG parameters are immediately related to those of the FG model (no need to parameterize a distinct CG model); (ii) nearly ideal mixing is expected for CG variants with similar s-values (ideal mixing holding in the limit of identical s-values); (iii) the solvent relaxation timescales should be preserved (no dynamical acceleration typical for coarse-graining); (iv) the graining level N_G (number of FG molecules represented by one CG molecule) can be chosen arbitrarily (in particular, N_G = s³ is not necessarily an integer); and (v) in an adaptive-resolution scheme, this level can be varied continuously as a function of the position (without requiring a bundling mechanism), and this variation occurs at a constant number of particles per molecule (no occurrence of fractional degrees of freedom in the buffer layer). By construction, the AdSoS scheme minimizes the thermodynamic mismatch between the different regions of the adaptive-resolution system, leading to a nearly homogeneous scaled solvent density s³ρ. Residual density artifacts in and at the surface of the boundary layer can easily be corrected by means of a grid-based biasing potential constructed in a preliminary pure-solvent simulation. This article introduces the AdSoS scheme and provides an initial application to pure atomic liquids (no solute) with Lennard-Jones plus Coulomb interactions in a slab geometry.

Comprehensive virtual screening of 4.8 k flavonoids reveals novel insights into allosteric inhibition of SARS-CoV-2 MPRO

SARS-CoV-2 main protease is a common target for inhibition assays due to its high conservation among coronaviruses. Since flavonoids show antiviral activity, several in silico works have proposed them as potential SARS-CoV-2 main protease inhibitors. Nonetheless, there is reason to doubt certain results given the lack of consideration for flavonoid promiscuity or main protease plasticity, usage of short library sizes, absence of control molecules and/or the limitation of the methodology to a single target site. Here, we report a virtual screening study where dorsilurin E, euchrenone a11, sanggenol O and CHEMBL2171598 are proposed to inhibit main protease through different pathways. Remarkably, novel structural mechanisms were observed after sanggenol O and CHEMBL2171598 bound to experimentally proven allosteric sites. The former drastically affected the active site, while the latter triggered a hinge movement which has been previously reported for an inactive SARS-CoV main protease mutant. The use of a curated database of 4.8 k flavonoids, combining two well-known docking software (AutoDock Vina and AutoDock4.2), molecular dynamics and MMPBSA, guaranteed an adequate analysis and robust interpretation. These criteria can be considered for future screening campaigns against SARS-CoV-2 main protease.

Computational methods for exploring protein conformations

Proteins are dynamic molecules that can transition between a potentially wide range of structures comprising their conformational ensemble. The nature of these conformations and their relative probabilities are described by a high-dimensional free energy landscape. While computer simulation techniques such as molecular dynamics simulations allow characterisation of the metastable conformational states and the transitions between them, and thus free energy landscapes, to be characterised, the barriers between states can be high, precluding efficient sampling without substantial computational resources. Over the past decades, a dizzying array of methods have emerged for enhancing conformational sampling, and for projecting the free energy landscape onto a reduced set of dimensions that allow conformational states to be distinguished, known as collective variables (CVs), along which sampling may be directed. Here, a brief description of what biomolecular simulation entails is followed by a more detailed exposition of the nature of CVs and methods for determining these, and, lastly, an overview of the myriad different approaches for enhancing conformational sampling, most of which rely upon CVs, including new advances in both CV determination and conformational sampling due to machine learning.

Extension of the Variational Free Energy Profile and Multistate Bennett Acceptance Ratio Methods for High-Dimensional Potential of Mean Force Profile Analysis

We redevelop the variational free energy profile (vFEP) method using a cardinal B-spline basis to extend the method for analyzing free energy surfaces (FESs) involving three or more reaction coordinates. We also implemented software for evaluating high-dimensional profiles based on the multistate Bennett acceptance ratio (MBAR) method which constructs an unbiased probability density from global reweighting of the observed samples. The MBAR method takes advantage of a fast algorithm for solving the unbinned weighted histogram (UWHAM)/MBAR equations which replaces the solution of simultaneous equations with a nonlinear optimization of a convex function. We make use of cardinal B-splines and multiquadric radial basis functions to obtain smooth, differentiable MBAR profiles in arbitrary high dimensions. The cardinal B-spline vFEP and MBAR methods are compared using three example systems that examine 1D, 2D, and 3D profiles. Both methods are found to be useful and produce nearly indistinguishable results. The vFEP method is found to be 150 times faster than MBAR when applied to periodic 2D profiles, but the MBAR method is 4.5 times faster than vFEP when evaluating unbounded 3D profiles. In agreement with previous comparisons, we find the vFEP method produces superior FESs when the overlap between umbrella window simulations decreases. Finally, the associative reaction mechanism of hammerhead ribozyme is characterized using 3D, 4D, and 6D profiles, and the higher-dimensional profiles are found to have smaller reaction barriers by as much as 1.5 kcal/mol. The methods presented here have been implemented into the FE-ToolKit software package along with new methods for network-wide free energy analysis in drug discovery.

Swarm-CG: Automatic Parametrization of Bonded Terms in MARTINI-Based Coarse-Grained Models of Simple to Complex Molecules via Fuzzy Self-Tuning Particle Swarm Optimization

We present Swarm-CG, a versatile software for the automatic iterative parametrization of bonded parameters in coarse-grained (CG) models, ideal in combination with popular CG force fields such as MARTINI. By coupling fuzzy self-tuning particle swarm optimization to Boltzmann inversion, Swarm-CG performs accurate bottom-up parametrization of bonded terms in CG models composed of up to 200 pseudo atoms within 4-24 h on standard desktop machines, using default settings. The software benefits from a user-friendly interface and two different usage modes (default and advanced). We particularly expect Swarm-CG to support and facilitate the development of new CG models for the study of complex molecular systems interesting for bio- and nanotechnology. Excellent performances are demonstrated using a benchmark of 9 molecules of diverse nature, structural complexity, and size. Swarm-CG is available with all its dependencies via the Python Package Index (PIP package: swarm-cg). Demonstration data are available at: www.github.com/GMPavanLab/SwarmCG.

High-throughput screening of metal-organic frameworks for kinetic separation of propane and propene

We apply molecular simulations to screen a database of reported metal-organic framework structures from the computation-ready, experimental (CoRE) MOF database to identify materials potentially capable of separating propane and propene by diffusion. We report a screening workflow that uses descriptor analysis, conventional molecular dynamics (MD), and Nudged Elastic Band (NEB) energy barrier calculations at both classical force field and Density Functional Theory (DFT) levels. For the first time, the effects of framework flexibility on guest transport properties were fully considered in a screening process and led to the identification of candidate MOFs. The hits identified by this proof-of-concept workflow include ZIF-8 and ZIF-67 previously shown to have large differences in propane and propene diffusivities as well as two other materials that have not been tested experimentally yet. This work emphasises the importance of taking into account framework flexibility when studying guest transport in porous materials, demonstrates the potential of the data-driven identification of high-performance materials and highlights the ways of improving the predictive power of the screening workflow.

Chirality Effects in Biomolecular Systems: Calculation of the Relative Free Energies by Molecular Dynamics Simulations

Chirality plays an essential role in chemical and biological sciences. At the molecular level, the effects associated with this phenomenon can be studied by using the well-established technique of molecular dynamics simulations. In this work, we present several approaches suited for the molecular dynamics-based free energy calculation in chiral systems. In particular, we have proposed and tested the following strategies relying on the application of general, enhanced sampling methods: (i) biased sampling in the two-dimensional space, along the coordinates defined by the values of the selected torsional angles; (ii) biased sampling in the one- or two-dimensional space, along the path-based coordinate(s); (iii) rational alteration of the system's Hamiltonian in order to enable the interconversion between stereoisomers and reweighting the biased distribution of configurations; (iv) using the free energy landscape generated within approaches (i) or (ii) as time-independent bias in order to further improve sampling efficiency and simultaneously account for multiple chiral centers. All approaches have been tested on a set of model compounds (fenoterol, fructofuranose, and bromochlorofluoromethane), demonstrating the good performance but also some differences in the range of their applicabilities.

Computationally predicting binding affinity in protein-ligand complexes: free energy-based simulations and machine learning-based scoring functions

Accurately predicting protein-ligand binding affinities can substantially facilitate the drug discovery process, but it remains as a difficult problem. To tackle the challenge, many computational methods have been proposed. Among these methods, free energy-based simulations and machine learning-based scoring functions can potentially provide accurate predictions. In this paper, we review these two classes of methods, following a number of thermodynamic cycles for the free energy-based simulations and a feature-representation taxonomy for the machine learning-based scoring functions. More recent deep learning-based predictions, where a hierarchy of feature representations are generally extracted, are also reviewed. Strengths and weaknesses of the two classes of methods, coupled with future directions for improvements, are comparatively discussed.

Investigating Conformational Dynamics of Lewis Y Oligosaccharides and Elucidating Blood Group Dependency of Cholera Using Molecular Dynamics

Cholera is caused by Vibrio cholerae and is an example of a blood-group-dependent disease. Recent studies suggest that the receptor-binding B subunit of the cholera toxin (CT) binds histo-blood group antigens at a secondary binding site. Herein, we studied the conformational dynamics of Lewis Y (Le^Y) oligosaccharides, H-tetrasaccharides and A-pentasaccharides, in aqueous solution by conducting accelerated molecular dynamics (aMD) simulations. The flexible nature of both oligosaccharides was displayed in aMD simulations. Furthermore, aMD simulations revealed that for both oligosaccharides in the free form, ⁴C₁ and ¹C₄ puckers were sampled for all but GalNAc monosaccharides, while either the ⁴C₁ (GlcNAc, Gal, GalNAc) or ¹C₄ (Fuc2, Fuc3) pucker was sampled in the CT-bound forms. In aMD, the complete transition from the ⁴C₁ to ¹C₄ pucker was sampled for GlcNAc and Gal in both oligosaccharides. Further, we have observed a transition from the open to closed conformer in the case of A-pentasaccharide, while H-tetrasaccharide remains in the open conformation throughout the simulation. Both oligosaccharides adopted an open conformation in the CT binding site. Moreover, we have investigated the molecular basis of recognition of Le^Y oligosaccharides by the B subunit of the cholera toxin of classical and El Tor biotypes using the molecular mechanics generalized Born surface area (MM/GBSA) scheme. The O blood group determinant, H-tetrasaccharide, exhibits a stronger affinity to both biotypes compared to the A blood group determinant, A-pentasaccharide, which agrees with the experimental data. The difference in binding free energy between O and A blood group determinants mainly arises due to the increased entropic cost and desolvation energy in the case of A-pentasaccharide compared to that of H-tetrasaccharide. Our study also reveals that the terminal Fuc3 contributes most to the binding free energy compared to other carbohydrate residues as it forms multiple hydrogen bonds with CT. Overall, our study might help in designing glycomimetic drugs targeting the cholera toxin.

Molecular Conformations of Di-, Tri-, and Tetra-α-(2→8)-Linked Sialic Acid from NMR Spectroscopy and MD Simulations

By using molecular dynamics simulations with an efficient enhanced sampling technique and in combination with nuclear magnetic resonance (NMR) spectroscopy quantitative structural information on α-2,8-linked sialic acids is presented. We used a bottom-up approach to obtain a set of larger ensembles for tetra- and deca-sialic acid from model dimer and trimer systems that are in agreement with the available J-coupling constants and nuclear Overhauser effects. The molecular dynamic (MD) simulations with enhanced sampling are used to validate the force field used in this study for its further use. This empowered us to couple NMR observables in the MD framework via J-coupling and distance restraining simulations to obtain conformations that are supported by experimental data. We used these conformations in thermodynamic integration and one-step perturbation simulations to calculate the free-energy of suggested helical conformations. This study brings most of the available NMR experiments together and supplies information to resolve the conflict on the structures of poly-α-2,8-linked sialic acid.

Multidimensional virtual-system coupled canonical molecular dynamics to compute free-energy landscapes of peptide multimer assembly

An enhanced-sampling method termed multidimensional virtual-system coupled canonical molecular dynamics (mD-VcMD) method is developed. In many cases, generalized-ensemble methods realizing enhanced sampling, for example, adaptive umbrella sampling, apply an effective potential, which is derived from temporarily assumed canonical distribution as a function of one or more arbitrarily defined reaction coordinates. However, it is not straightforward to estimate the appropriate canonical distribution, especially for cases applying multiple reaction coordinates. The current method, mD-VcMD, does not rely on the form of the canonical distribution. Therefore, it is practically useful to explore a high-dimensional reaction-coordinate space. In this article, formulation of mD-VcMD and its evaluation with the simple molecular models consisting of three or four alanine peptides are presented. We confirmed that mD-VcMD efficiently searched 2D and 3D reaction-coordinate spaces defined as interpeptide distances. Direct comparisons with results of long-term canonical MD simulations revealed that mD-VcMD produces correct canonical ensembles. © 2019 Wiley Periodicals, Inc.

Distinct Fcα receptor N-glycans modulate the binding affinity to immunoglobulin A (IgA) antibodies

Human immunoglobulin A (IgA) is the most prevalent antibody class at mucosal sites with an important role in mucosal defense. Little is known about the impact of N-glycan modifications of IgA1 and IgA2 on binding to the Fcα receptor (FcαRI), which is also heavily glycosylated at its extracellular domain. Here, we transiently expressed human epidermal growth factor receptor 2 (HER2)-binding monomeric IgA1, IgA2m(1), and IgA2m(2) variants in Nicotiana benthamiana ΔXT/FT plants lacking the enzymes responsible for generating nonhuman N-glycan structures. By coinfiltrating IgA with the respective glycan-modifying enzymes, we generated IgA carrying distinct homogenous N-glycans. We demonstrate that distinctly different N-glycan profiles did not influence antigen binding or the overall structure and integrity of the IgA antibodies but did affect their thermal stability. Using size-exclusion chromatography, differential scanning and isothermal titration calorimetry, surface plasmon resonance spectroscopy, and molecular modeling, we probed distinct IgA1 and IgA2 glycoforms for binding to four different FcαRI glycoforms and investigated the thermodynamics and kinetics of complex formation. Our results suggest that different N-glycans on the receptor significantly contribute to binding affinities for its cognate ligand. We also noted that full-length IgA and FcαRI form a mixture of 1:1 and 1:2 complexes tending toward a 1:1 stoichiometry due to different IgA tailpiece conformations that make it less likely that both binding sites are simultaneously occupied. In conclusion, N-glycans of human IgA do not affect its structure and integrity but its thermal stability, and FcαRI N-glycans significantly modulate binding affinity to IgA.

Enhanced sampling in molecular dynamics

Although molecular dynamics simulations have become a useful tool in essentially all fields of chemistry, condensed matter physics, materials science, and biology, there is still a large gap between the time scale which can be reached in molecular dynamics simulations and that observed in experiments. To address the problem, many enhanced sampling methods were introduced, which effectively extend the time scale being approached in simulations. In this perspective, we review a variety of enhanced sampling methods. We first discuss collective-variables-based methods including metadynamics and variationally enhanced sampling. Then, collective variable free methods such as parallel tempering and integrated tempering methods are presented. At last, we conclude with a brief introduction of some newly developed combinatory methods. We summarize in this perspective not only the theoretical background and numerical implementation of these methods but also the new challenges and prospects in the field of the enhanced sampling.

Structural Aspects of the O-glycosylation Linkage in Glycopeptides via MD Simulations and Comparison with NMR Experiments

A powerful conformational searching and enhanced sampling simulation method, and unbiased molecular dynamics simulations have been used along with NMR spectroscopic observables to provide a detailed structural view of O-glycosylation. For four model systems, the force-field parameters can accurately predict experimental NMR observables (J couplings and NOE's). This enables us to derive conclusions based on the generated ensembles, in which O-glycosylation affects the peptide backbone conformation by forcing it towards to an extended conformation. An exception is described for β-GalNAc-Thr where the α content is increased and stabilized via hydrogen bonding between the sugar and the peptide backbone, which was not observed in the rest of the studied systems. These observations might offer an explanation for the evolutionary preference of α-linked GalNAc glycosylation instead of a β link.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Orientational and Folding Thermodynamics via Electric Dipole Moment Restraining

The projection of molecular processes onto a small set of relevant descriptors, the so-called reaction coordinates or collective variables (CVs), is a technique nowadays routinely employed by the biomolecular simulation community. In this work, we implemented two CVs to manipulate the orientation (i.e., angle) (μ⃗_a) and magnitude (|μ⃗|) of the electric dipole moment. In doing so, we studied the thermodynamics of water orientation under the application of external voltages and the folding of two polypeptides at zero-field conditions. The projection of the free-energy [potential of mean force (PMF)] along water orientation defined an upper limit of around 0.3 V for irrelevant thermodynamic effects. On the other hand, sufficiently strong μ⃗_a restraints applied on 12-alanine (Ala₁₂) triggered structural effects because of the alignment of local dipoles; for lower restraints, a full-body rotation is achieved. The manipulation of |μ⃗| produced strong perturbations on the secondary structure of Ala₁₂, promoting an enhanced sampling to its configurational space. Rigorous free-energy calculations in the form of 2-D PMFs for deca-alanine showed the utility of |μ⃗| as a reaction coordinate to study folding in small α helices. As a whole, we propose that the manipulation of both components of the dipole moment, μ⃗_a and |μ⃗|, provides thermodynamics insights into the structural conformation and stability of biomolecules. These new CVs are implemented in the Colvars module, available for NAMD and LAMMPS.

Using correlated motions to determine sufficient sampling times for molecular dynamics

Here we present a time-dependent correlation method that provides insight into how long a system takes to grow into its equal-time (Pearson) correlation. We also show a usage of an extant time-lagged correlation method that indicates the time for parts of a system to become decorrelated, relative to equal-time correlation. Given a completed simulation (or set of simulations), these tools estimate (i) how long of a simulation of the same system would be sufficient to observe the same correlated motions, (ii) if patterns of observed correlated motions indicate events beyond the timescale of the simulation, and (iii) how long of a simulation is needed to observe these longer timescale events. We view this method as a decision-support tool that will aid researchers in determining necessary sampling times. In principle, this tool is extendable to any multidimensional time series data with a notion of correlated fluctuations; however, here we limit our discussion to data from molecular-dynamics simulations.

Identification of Rare Lewis Oligosaccharide Conformers in Aqueous Solution Using Enhanced Sampling Molecular Dynamics

Determining the conformations accessible to carbohydrate ligands in aqueous solution is important for understanding their biological action. In this work, we evaluate the conformational free-energy surfaces of Lewis oligosaccharides in explicit aqueous solvent using a multidimensional variant of the swarm-enhanced sampling molecular dynamics (msesMD) method; we compare with multi-microsecond unbiased MD simulations, umbrella sampling, and accelerated MD approaches. For the sialyl Lewis A tetrasaccharide, msesMD simulations in aqueous solution predict conformer landscapes in general agreement with the other biased methods and with triplicate unbiased 10 μs trajectories; these simulations find a predominance of closed conformer and a range of low-occupancy open forms. The msesMD simulations also suggest closed-to-open transitions in the tetrasaccharide are facilitated by changes in ring puckering of its GlcNAc residue away from the ⁴C₁ form, in line with previous work. For sialyl Lewis X tetrasaccharide, msesMD simulations predict a minor population of an open form in solution corresponding to a rare lectin-bound pose observed crystallographically. Overall, from comparison with biased MD calculations, we find that triplicate 10 μs unbiased MD simulations may not be enough to fully sample glycan conformations in aqueous solution. However, the computational efficiency and intuitive approach of the msesMD method suggest potential for its application in glycomics as a tool for analysis of oligosaccharide conformation.

Molecular dynamics coupled with a virtual system for effective conformational sampling

An enhanced conformational sampling method is proposed: virtual-system coupled canonical molecular dynamics (VcMD). Although VcMD enhances sampling along a reaction coordinate, this method is free from estimation of a canonical distribution function along the reaction coordinate. This method introduces a virtual system that does not necessarily obey a physical law. To enhance sampling the virtual system couples with a molecular system to be studied. Resultant snapshots produce a canonical ensemble. This method was applied to a system consisting of two short peptides in an explicit solvent. Conventional molecular dynamics simulation, which is ten times longer than VcMD, was performed along with adaptive umbrella sampling. Free-energy landscapes computed from the three simulations mutually converged well. The VcMD provided quicker association/dissociation motions of peptides than the conventional molecular dynamics did. The VcMD method is applicable to various complicated systems because of its methodological simplicity. © 2018 Wiley Periodicals, Inc.

Validation and Comparison of Force Fields for Native Cyclodextrins in Aqueous Solution

Molecular dynamics simulations of native α-, β-, and γ-cyclodextrin in aqueous solution have been conducted with the goal to investigate the performance of the CHARMM36 force field, the AMBER-compatible q4md-CD force field, and five variants of the GROMOS force field. The properties analyzed are structural parameters derived from X-ray diffraction and NMR experiments as well as hydrogen bonds and hydration patterns, including hydration free enthalpies. Recent revisions of the torsional-angle parameters for carbohydrate systems within the GROMOS family of force fields lead to a significant improvement of the agreement between simulated and experimental NMR data. Therefore, we recommend using the variant 53A6_GLYC instead of 53A6 and 56A6_{CARBO_R} or 2016H66 instead of 56A6_CARBO to simulate cyclodextrins in solution. The CHARMM36 and q4md-CD force fields show a similar performance as the three recommended GROMOS parameter sets. A significant difference is the more flexible nature of the cyclodextrins modeled with the CHARMM36 and q4md-CD force fields compared to the three recommended GROMOS parameter sets.

Ligand-induced action of the W2866.48 rotamer toggle switch in the β2-adrenergic receptor

Studies focused on GPCRs, particularly on the β₂-adrenergic receptor (β₂-AR), have demonstrated the relationship between ligand structure, receptor conformational changes and the corresponding pharmacological outcomes. Herein, we studied the molecular details of the rotameric flip of the W286^6.48 sidechain, i.e. a presumed action switch that has not been reported in native β₂-AR thus far. It is believed that although both the 'active' and 'inactive' conformers of β₂-AR exhibit similar conformations of this switch, it may still play a substantial role in the ligand-induced activation of the receptor. By using both experimental methods (time-resolved fluorescence spectroscopy) and molecular modeling techniques (enhanced-sampling molecular dynamics), we characterized the conformational rearrangements of W286^6.48 in relation to the type of ligand present in the binding cavity and to the conformation of the receptor ('active' vs. 'inactive' β₂-AR). We found that the conformational behaviour of W286^6.48 is correlated with the pharmacological character of the ligand present in the binding cavity but not with the instantaneous conformation of the receptor. Namely, agonists promote the W286^6.48 conformations that facilitate the increase of the solvation within the inner receptor channel. In contrast, antagonists and inverse agonists act toward the decrease of the solvation in the inner channel. This creates an opportunity for using computational methodologies in determining the pharmacological properties of various ligands. The combination of the time-resolved fluorescence spectroscopy technique with the enhanced-sampling molecular dynamics simulations is shown to be a powerful tool for studying the ligand-induced conformational rearrangements in GPCRs.

Modeling of Oligosaccharides within Glycoproteins from Free-Energy Landscapes

In spite of the abundance of glycoproteins in biological processes, relatively little three-dimensional structural data is available for glycan structures. Here, we study the structure and flexibility of the vast majority of mammalian oligosaccharides appearing in N- and O-glycosylated proteins using a bottom up approach. We report the conformational free-energy landscapes of all relevant glycosidic linkages as obtained from local elevation simulations and subsequent umbrella sampling. To the best of our knowledge, this represents the first complete conformational library for the construction of N- and O-glycan structures. Next, we systematically study the effect of neighboring residues, by extensively simulating all relevant trisaccharides and one tetrasaccharide. This allows for an unprecedented comparison of disaccharide linkages in large oligosaccharides. With a small number of exceptions, the conformational preferences in the larger structures are very similar as in the disaccharides. This, finally, allows us to suggest several efficient approaches to construct complete N- and O-glycans on glycoproteins, as exemplified on two relevant examples.

Exploration, Sampling, And Reconstruction of Free Energy Surfaces with Gaussian Process Regression

Practical free energy reconstruction algorithms involve three separate tasks: biasing, measuring some observable, and finally reconstructing the free energy surface from those measurements. In more than one dimension, adaptive schemes make it possible to explore only relatively low lying regions of the landscape by progressively building up the bias toward the negative of the free energy surface so that free energy barriers are eliminated. Most schemes use the final bias as their best estimate of the free energy surface. We show that large gains in computational efficiency, as measured by the reduction of time to solution, can be obtained by separating the bias used for dynamics from the final free energy reconstruction itself. We find that biasing with metadynamics, measuring a free energy gradient estimator, and reconstructing using Gaussian process regression can give an order of magnitude reduction in computational cost.

Free-energy calculations of residue mutations in a tripeptide using various methods to overcome inefficient sampling

Previous free-energy calculations have shown that the seemingly simple transformation of the tripeptide KXK to KGK in water holds some unobvious challenges concerning the convergence of the forward and backward thermodynamic integration processes (i.e., hysteresis). In the current study, the central residue X was either alanine, serine, glutamic acid, lysine, phenylalanine, or tyrosine. Interestingly, the transformation from alanine to glycine yielded the highest hysteresis in relation to the extent of the chemical change of the side chain. The reason for that could be attributed to poor sampling of φ2 /ψ2 dihedral angles along the transformation. Altering the nature of alanine's Cβ atom drastically improved the sampling and at the same time led to the identification of high energy barriers as cause for it. Consequently, simple strategies to overcome these barriers are to increase simulation time (computationally expensive) or to use enhanced sampling techniques such as Hamiltonian replica exchange molecular dynamics and one-step perturbation. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Multistate λ-local-elevation umbrella-sampling (MS-λ-LEUS): method and application to the complexation of cations by crown ethers

An extension of the λ-local-elevation umbrella-sampling (λ-LEUS) scheme [ Bieler et al. J. Chem. Theory Comput. 2014 , 10 , 3006 ] is proposed to handle the multistate (MS) situation, i.e. the calculation of the relative free energies of multiple physical states based on a single simulation. The key element of the MS-λ-LEUS approach is to use a single coupling variable Λ controlling successive pairwise mutations between the states of interest in a cyclic fashion. The Λ variable is propagated dynamically as an extended-system variable, using a coordinate transformation with plateaus and a memory-based biasing potential as in λ-LEUS. Compared to other available MS schemes (one-step perturbation, enveloping distribution sampling and conventional λ-dynamics) the proposed method presents a number of important advantages, namely: (i) the physical states are visited explicitly and over finite time periods; (ii) the extent of unphysical space required to ensure transitions is kept minimal and, in particular, one-dimensional; (iii) the setup protocol solely requires the topologies of the physical states; and (iv) the method only requires limited modifications in a simulation code capable of handling two-state mutations. As an initial application, the absolute binding free energies of five alkali cations to three crown ethers in three different solvents are calculated. The results are found to reproduce qualitatively the main experimental trends and, in particular, the experimental selectivity of 18C6 for K(+) in water and methanol, which is interpreted in terms of opposing trends along the cation series between the solvation free energy of the cation and the direct electrostatic interactions within the complex.

Theoretical Investigation of Solvent Effects on Glycosylation Reactions: Stereoselectivity Controlled by Preferential Conformations of the Intermediate Oxacarbenium-Counterion Complex

The mechanism of solvent effects on the stereoselectivity of glycosylation reactions is investigated using quantum-mechanical (QM) calculations and molecular dynamics (MD) simulations, considering a methyl-protected glucopyranoside triflate as a glycosyl donor equivalent and the solvents acetonitrile, ether, dioxane, or toluene, as well as gas-phase conditions (vacuum). The QM calculations on oxacarbenium-solvent complexes do not provide support to the usual solvent-coordination hypothesis, suggesting that an experimentally observed β-selectivity (α-selectivity) is caused by the preferential coordination of a solvent molecule to the reactive cation on the α-side (β-side) of the anomeric carbon. Instead, explicit-solvent MD simulations of the oxacarbenium-counterion (triflate ion) complex (along with corresponding QM calculations) are compatible with an alternative mechanism, termed here the conformer and counterion distribution hypothesis. This new hypothesis suggests that the stereoselectivity is dictated by two interrelated conformational properties of the reactive complex, namely, (1) the conformational preferences of the oxacarbenium pyranose ring, modulating the steric crowding and exposure of the anomeric carbon toward the α or β face, and (2) the preferential coordination of the counterion to the oxacarbenium cation on one side of the anomeric carbon, hindering a nucleophilic attack from this side. For example, in acetonitrile, the calculations suggest a dominant B2,5 ring conformation of the cation with preferential coordination of the counterion on the α side, both factors leading to the experimentally observed β selectivity. Conversely, in dioxane, they suggest a dominant (4)H3 ring conformation with preferential counterion coordination on the β side, both factors leading to the experimentally observed α selectivity.

Revision of the GROMOS 56A6(CARBO) force field: Improving the description of ring-conformational equilibria in hexopyranose-based carbohydrates chains

This article describes a revised version 56A6(CARBO_R) of the GROMOS 56A6(CARBO) force field for hexopyranose-based carbohydrates. The simulated properties of unfunctionalized hexopyranoses are unaltered with respect to 56A6CARBO . In the context of both O1 -alkylated hexopyranoses and oligosaccharides, the revision stabilizes the regular (4) C1 chair for α-anomers, with the opposite effect for β-anomers. As a result, spurious ring inversions observed in α(1→4)-linked chains when using the original 56A6(CARBO) force field are alleviated. The (4) C1 chair is now the most stable conformation for all d-hexopyranose residues, irrespective of the linkage type and anomery, and of the position of the residue along the chain. The methylation of a d-hexopyranose leads to a systematic shift in the ring-inversion free energy ((4) C1 to (1) C4 ) by 7-8 kJ mol(-1), positive for the α-anomers and negative for the β-anomers, which is qualitatively compatible with the expected enhancement of the anomeric effect upon methylation at O1. The ring-inversion free energies for residues within chains are typically smaller in magnitude compared to those of the monomers, and correlate rather poorly with the latter. This suggests that the crowding of ring substituents upon chain formation alters the ring flexibility in a nonsystematic fashion. In general, the description of carbohydrate chains afforded by 56A6(CARBO_R) suggests a significant extent of ring flexibility, i.e., small but often non-negligible equilibrium populations of inverted chairs, and challenges the "textbook" picture of conformationally locked carbohydrate rings.

Combine umbrella sampling with integrated tempering method for efficient and accurate calculation of free energy changes of complex energy surface

Umbrella sampling is an efficient method for the calculation of free energy changes of a system along well-defined reaction coordinates. However, when there exist multiple parallel channels along the reaction coordinate or hidden barriers in directions perpendicular to the reaction coordinate, it is difficult for conventional umbrella sampling to reach convergent sampling within limited simulation time. Here, we propose an approach to combine umbrella sampling with the integrated tempering sampling method. The umbrella sampling method is applied to chemically more relevant degrees of freedom that possess significant barriers. The integrated tempering sampling method is used to facilitate the sampling of other degrees of freedom which may possess statistically non-negligible barriers. The combined method is applied to two model systems, butane and ACE-NME molecules, and shows significantly improved sampling efficiencies as compared to standalone conventional umbrella sampling or integrated tempering sampling approaches. Further analyses suggest that the enhanced performance of the new method come from the complemented advantages of umbrella sampling with a well-defined reaction coordinate and integrated tempering sampling in orthogonal space. Therefore, the combined approach could be useful in the simulation of biomolecular processes, which often involves sampling of complex rugged energy landscapes.

Computational Calorimetry: High-Precision Calculation of Host-Guest Binding Thermodynamics

We present a strategy for carrying out high-precision calculations of binding free energy and binding enthalpy values from molecular dynamics simulations with explicit solvent. The approach is used to calculate the thermodynamic profiles for binding of nine small molecule guests to either the cucurbit[7]uril (CB7) or β-cyclodextrin (βCD) host. For these systems, calculations using commodity hardware can yield binding free energy and binding enthalpy values with a precision of ∼0.5 kcal/mol (95% CI) in a matter of days. Crucially, the self-consistency of the approach is established by calculating the binding enthalpy directly, via end point potential energy calculations, and indirectly, via the temperature dependence of the binding free energy, i.e., by the van't Hoff equation. Excellent agreement between the direct and van't Hoff methods is demonstrated for both host-guest systems and an ion-pair model system for which particularly well-converged results are attainable. Additionally, we find that hydrogen mass repartitioning allows marked acceleration of the calculations with no discernible cost in precision or accuracy. Finally, we provide guidance for accurately assessing numerical uncertainty of the results in settings where complex correlations in the time series can pose challenges to statistical analysis. The routine nature and high precision of these binding calculations opens the possibility of including measured binding thermodynamics as target data in force field optimization so that simulations may be used to reliably interpret experimental data and guide molecular design.

Toward Improved Force-Field Accuracy through Sensitivity Analysis of Host-Guest Binding Thermodynamics

Improving the capability of atomistic computer models to predict the thermodynamics of noncovalent binding is critical for successful structure-based drug design, and the accuracy of such calculations remains limited by nonoptimal force field parameters. Ideally, one would incorporate protein-ligand affinity data into force field parametrization, but this would be inefficient and costly. We now demonstrate that sensitivity analysis can be used to efficiently tune Lennard-Jones parameters of aqueous host-guest systems for increasingly accurate calculations of binding enthalpy. These results highlight the promise of a comprehensive use of calorimetric host-guest binding data, along with existing validation data sets, to improve force field parameters for the simulation of noncovalent binding, with the ultimate goal of making protein-ligand modeling more accurate and hence speeding drug discovery.