Metadynamics convergence law in a multidimensional system

Expanded Functionality and Portability for the Colvars Library

Colvars is an open-source C++ library that provides a modular toolkit for collective-variable-based molecular simulations. It allows practitioners to easily create and implement descriptors that best fit a process of interest and to apply a wide range of biasing algorithms in collective variable space. This paper reviews several features and improvements to Colvars that were added since its original introduction. Special attention is given to contributions that significantly expanded the capabilities of this software or its distribution with major MD simulation packages. Collective variables can now be optimized either manually or by machine-learning methods, and the space of descriptors can be explored interactively using the graphical interface included in VMD. Beyond the spatial coordinates of individual molecules, Colvars can now apply biasing forces to mesoscale structures and alchemical degrees of freedom and perform simulations guided by experimental data within ensemble averages or probability distributions. It also features advanced computational schemes to boost the accuracy, robustness, and general applicability of simulation methods, including extended-system and multiple-walker adaptive biasing force, boundary conditions for metadynamics, replica exchange with biasing potentials, and adiabatic bias molecular dynamics. The library is made available directly within the main distributions of the academic software GROMACS, LAMMPS, NAMD, Tinker-HP, and VMD. The robustness of the software and the reliability of the results are ensured through the use of continuous integration with a test suite within the source repository.

Papain-Mediated Conjugation of Peptide Nucleic Acids to Delivery Peptides: A Density Functional Theory/Molecular Mechanics Metadynamics Study in Aqueous and Organic Solvent

Enzymatic peptide synthesis is a powerful alternative to solid-phase methods, as enzymes can have high regio- and stereoselectivity and high yield and require mild reaction conditions. This is beneficial in formulation research due to the rise of nucleic acid therapies. Peptide nucleic acids (PNAs) have a high affinity toward DNA and RNA, and their solubility and cellular delivery can be improved via conjugation to peptides. Here, we designed and assessed the viability of the papain enzyme to conjugate four PNA-peptide models in water and an organic solvent using QM/MM metadynamics. We found that the reactions in water yield better results, where three conjugates could potentially be synthesized by the enzyme, with the first transition state as the rate-limiting step, with an associated energy of 14.53 kcal mol-1, although with a slight endergonic profile. The results highlight the importance of considering the enzyme pockets and different substrate acceptivities and contribute to developing greener, direct, and precise synthetic routes for nucleic acid-based therapies. By exploring the enzyme's potential in conjunction with chemical synthesis, current protocols can be simplified for the synthesis of longer nucleic acids and peptide sequences (and, by extension, proteins) from smaller oligo or peptide blocks.

Improving the Efficiency of Variationally Enhanced Sampling with Wavelet-Based Bias Potentials

Collective variable-based enhanced sampling methods are routinely used on systems with metastable states, where high free energy barriers impede the proper sampling of the free energy landscapes when using conventional molecular dynamics simulations. One such method is variationally enhanced sampling (VES), which is based on a variational principle where a bias potential in the space of some chosen slow degrees of freedom, or collective variables, is constructed by minimizing a convex functional. In practice, the bias potential is taken as a linear expansion in some basis function set. So far, primarily basis functions delocalized in the collective variable space, like plane waves, Chebyshev, or Legendre polynomials, have been used. However, there has not been an extensive study of how the convergence behavior is affected by the choice of the basis functions. In particular, it remains an open question if localized basis functions might perform better. In this work, we implement, tune, and validate Daubechies wavelets as basis functions for VES. The wavelets construct orthogonal and localized bases that exhibit an attractive multiresolution property. We evaluate the performance of wavelet and other basis functions on various systems, going from model potentials to the calcium carbonate association process in water. We observe that wavelets exhibit excellent performance and much more robust convergence behavior than all other basis functions, as well as better performance than metadynamics. In particular, using wavelet bases yields far smaller fluctuations of the bias potential within individual runs and smaller differences between independent runs. Based on our overall results, we can recommend wavelets as basis functions for VES.

Differential ion dehydration energetics explains selectivity in the non-canonical lysosomal K+ channel TMEM175

Structures of the human lysosomal K⁺ channel transmembrane protein 175 (TMEM175) in open and closed states revealed a novel architecture lacking the canonical K⁺ selectivity filter motif present in previously known K⁺ channel structures. A hydrophobic constriction composed of four isoleucine residues was resolved in the pore and proposed to serve as the gate in the closed state, and to confer ion selectivity in the open state. Here, we achieve higher-resolution structures of the open and closed states and employ molecular dynamics simulations to analyze the conducting properties of the putative open state, demonstrating that it is permeable to K⁺ and, to a lesser degree, also Na⁺. Both cations must dehydrate significantly to penetrate the narrow hydrophobic constriction, but ion flow is assisted by a favorable electrostatic field generated by the protein that spans the length of the pore. The balance of these opposing energetic factors explains why permeation is feasible, and why TMEM175 is selective for K⁺ over Na⁺, despite the absence of the canonical selectivity filter. Accordingly, mutagenesis experiments reveal an exquisite sensitivity of the channel to perturbations that mitigate the constriction. Together, these data reveal a novel mechanism for selective permeation of ions by TMEM175 that is unlike that of other K⁺ channels.

Peptide Dynamics and Metadynamics: Leveraging Enhanced Sampling Molecular Dynamics to Robustly Model Long-Timescale Transitions

Molecular dynamics simulations can in theory reveal the thermodynamics and kinetics of peptide conformational transitions at atomic-level resolution. However, even with modern computing power, they are limited in the timescales they can sample, which is especially problematic for peptides that are fully or partially disordered. Here, we discuss how the enhanced sampling methods accelerated molecular dynamics (aMD) and metadynamics can be leveraged in a complementary fashion to quickly explore conformational space and then robustly quantify the underlying free energy landscape. We apply these methods to two peptides that have an intrinsically disordered nature, the histone H3 and H4 N-terminal tails, and use metadynamics to compute the free energy landscape along collective variables discerned from aMD simulations. Results show that these peptides are largely disordered, with a slight preference for α-helical structures.

Assessing the Intestinal Permeability of Small Molecule Drugs via Diffusion Motion on a Multidimensional Free Energy Surface

A protocol that accurately assesses the intestinal permeability of small molecule compounds plays an essential role in decreasing the cost and time in inventing a new drug. This manuscript presents a novel computational method to study the passive permeation of small molecule drugs based on the inhomogeneous solubility-diffusion model. The multidimensional free energy surface of the drug transiting through a lipid bilayer is computed with transition-tempered metadynamics that accurately captures the mechanisms of passive permeation. The permeability is computed by following the diffusion motion of the drug molecules along the minimal free energy path found on the multidimensional free energy surface. This computational method is assessed by studying the permeability of five small molecule drugs (ketoprofen, naproxen, metoprolol, propranolol, and salicylic acid). The results demonstrate a remarkable agreement between the computed permeabilities and those measured with the intestinal assay. The in silico method reported in this manuscript also reproduces the permeability measured from the intestinal assay (in vivo) better than the cell-based assays (e.g., PAMPA and Caco-2) do. In addition, the multidimensional free energy surface reveals the interplay between the structure of the small molecule and its permeability, shedding light on strategies of drug optimization.

Force-Correction Analysis Method for Derivation of Multidimensional Free-Energy Landscapes from Adaptively Biased Replica Simulations

A methodology is proposed for the calculation of multidimensional free-energy landscapes of molecular systems, based on analysis of multiple molecular dynamics trajectories wherein adaptive biases have been applied to enhance the sampling of different collective variables. In this approach, which we refer to as the Force-Correction Analysis Method (FCAM), local averages of the total and biasing forces are evaluated post hoc, and the latter are subtracted from the former to obtain unbiased estimates of the mean force across collective-variable space. Multidimensional free-energy surfaces and minimum free-energy pathways are then derived by integrating the mean-force landscape with a kinetic Monte Carlo algorithm. To evaluate the proposed method, a series of numerical tests and comparisons with existing approaches were carried out for small molecules, peptides, and proteins, based on all-atom trajectories generated with standard, concurrent, and replica-exchange metadynamics in collective-variable spaces ranging from one to six dimensional. The tests confirm the correctness of the FCAM formulation and demonstrate that calculated mean forces and free energies converge rapidly and accurately, outperforming other methods used to unbias this kind of simulation data.

Small molecules that inhibit TNF signalling by stabilising an asymmetric form of the trimer

Tumour necrosis factor (TNF) is a cytokine belonging to a family of trimeric proteins; it has been shown to be a key mediator in autoimmune diseases such as rheumatoid arthritis and Crohn’s disease. While TNF is the target of several successful biologic drugs, attempts to design small molecule therapies directed to this cytokine have not led to approved products. Here we report the discovery of potent small molecule inhibitors of TNF that stabilise an asymmetrical form of the soluble TNF trimer, compromising signalling and inhibiting the functions of TNF in vitro and in vivo. This discovery paves the way for a class of small molecule drugs capable of modulating TNF function by stabilising a naturally sampled, receptor-incompetent conformation of TNF. Furthermore, this approach may prove to be a more general mechanism for inhibiting protein–protein interactions.

While biologics have been successfully applied in TNF antagonist treatments, there are no clinically approved small molecules that target TNF. Here, the authors discover potent small molecule inhibitors of TNF, elucidate their molecular mechanism, and demonstrate TNF inhibition in vitro and in vivo.

Naphthalene crystal shape prediction from molecular dynamics simulations

The crystal shape of naphthalene grown from ethanol solution at constant supersaturation was predicted using state-of-the-art molecular dynamics simulations.

Interplay of cysteine exposure and global protein dynamics in small-molecule recognition by a regulator of G-protein signaling protein

Regulator of G protein signaling (RGS) proteins play a pivotal role in regulation of G protein-coupled receptor (GPCR) signaling and are therefore becoming an increasingly important therapeutic target. Recently discovered thiadiazolidinone (TDZD) compounds that target cysteine residues have shown different levels of specificities and potencies for the RGS4 protein, thereby suggesting intrinsic differences in dynamics of this protein upon binding of these compounds. In this work, we investigated using atomistic molecular dynamics (MD) simulations the effect of binding of several small-molecule inhibitors on perturbations and dynamical motions in RGS4. Specifically, we studied two conformational models of RGS4 in which a buried cysteine residue is solvent-exposed due to side-chain motions or due to flexibility in neighboring helices. We found that TDZD compounds with aromatic functional groups perturb the RGS4 structure more than compounds with aliphatic functional groups. Moreover, small-molecules with aromatic functional groups but lacking sulfur atoms only transiently reside within the protein and spontaneously dissociate to the solvent. We further measured inhibitory effects of TDZD compounds using a protein-protein interaction assay on a single-cysteine RGS4 protein showing trends in potencies of compounds consistent with our simulation studies. Thermodynamic analyses of RGS4 conformations in the apo-state and on binding to TDZD compounds revealed links between both conformational models of RGS4. The exposure of cysteine side-chains appears to facilitate initial binding of TDZD compounds followed by migration of the compound into a bundle of four helices, thereby causing allosteric perturbations in the RGS/Gα protein-protein interface.

Adsorption Mechanisms of Nucleobases on the Hydrated Au(111) Surface

The solution environment is of fundamental importance in the adsorption of molecules on surfaces, a process that is strongly affected by the capability of the adsorbate to disrupt the hydration layer above the surface. Here we disclose how the presence of interface water influences the adsorption mechanism of DNA nucleobases on a gold surface. By means of metadynamics simulations, we describe the distinctive features of a complex free-energy landscape for each base, which manifests activation barriers for the adsorption process. We characterize the different pathways that allow each nucleobase to overcome the barriers and be adsorbed on the surface, discussing how they influence the kinetics of adsorption of single-stranded DNA oligomers with homogeneous sequences. Our findings offer a rationale as to why experimental data on the adsorption of single-stranded homo-oligonucleotides do not straightforwardly follow the thermodynamics affinity rank.

Confidence Analysis of DEER Data and Its Structural Interpretation with Ensemble-Biased Metadynamics

Given its ability to measure multicomponent distance distributions between electron-spin probes, double electron-electron resonance (DEER) spectroscopy has become a leading technique to assess the structural dynamics of biomolecules. However, methodologies to evaluate the statistical error of these distributions are not standard, often hampering a rigorous interpretation of the experimental results. Distance distributions are often determined from the experimental DEER data through a mathematical method known as Tikhonov regularization, but this approach makes rigorous error estimates difficult. Here, we build upon an alternative, model-based approach in which the distance probability distribution is represented as a sum of Gaussian components, and use propagation of errors to calculate an associated confidence band. Our approach considers all sources of uncertainty, including the experimental noise, the uncertainty in the fitted background signal, and the limited time span of the data collection. The resulting confidence band reveals the most and least reliable features of the probability distribution, thereby informing the structural interpretation of DEER experiments. To facilitate this interpretation, we also generalize the molecular simulation method known as ensemble-biased metadynamics (EBMetaD). This method, originally designed to generate maximal-entropy structural ensembles consistent with one or more probability distributions, now also accounts for the uncertainty in those target distributions exactly as dictated by their confidence bands. After careful benchmarks, we demonstrate the proposed techniques using DEER results from spin-labeled T4 lysozyme.

Combining theoretical and experimental data to decipher CFTR 3D structures and functions

Cryo-electron microscopy (cryo-EM) has recently provided invaluable experimental data about the full-length cystic fibrosis transmembrane conductance regulator (CFTR) 3D structure. However, this experimental information deals with inactive states of the channel, either in an apo, quiescent conformation, in which nucleotide-binding domains (NBDs) are widely separated or in an ATP-bound, yet closed conformation. Here, we show that 3D structure models of the open and closed forms of the channel, now further supported by metadynamics simulations and by comparison with the cryo-EM data, could be used to gain some insights into critical features of the conformational transition toward active CFTR forms. These critical elements lie within membrane-spanning domains but also within NBD1 and the N-terminal extension, in which conformational plasticity is predicted to occur to help the interaction with filamin, one of the CFTR cellular partners.

Overfill Protection and Hyperdynamics in Adaptively Biased Simulations

Two problems associated with adaptively biased simulations are considered: overfilling and time scale estimation. First, a simple and computationally efficient procedure for limiting the bias fill-depth of any adaptive biasing potential is introduced. The resulting bias potential floods only up to a specified level and avoids bias accumulation in higher regions of free energy. Second, hyperdynamics can be invoked in combination with this depth-limited bias potential to estimate time scales. We argue that allowing the bias to equilibrate to a set depth is one crux in efficient hyperdynamics application. In both simple (alanine dipeptide) and complex (ligand residence) test cases, the useable boost factors are 6-8 times larger when the bias fill depth is directly controlled, as opposed to controlling the update frequency. Update frequency has emerged as a proxy for controlling the fill height of adaptive biases, with infrequent metadynamics being an example. The Kolmogorov-Smirnov test is shown to be insufficient for determining the trustworthiness of a hyperdynamics simulation, and a more robust strategy is described. Overfill limiting and hyperdynamics have been added to the fABMACS project.

Optimal updating magnitude in adaptive flat-distribution sampling

We present a study on the optimization of the updating magnitude for a class of free energy methods based on flat-distribution sampling, including the Wang-Landau (WL) algorithm and metadynamics. These methods rely on adaptive construction of a bias potential that offsets the potential of mean force by histogram-based updates. The convergence of the bias potential can be improved by decreasing the updating magnitude with an optimal schedule. We show that while the asymptotically optimal schedule for the single-bin updating scheme (commonly used in the WL algorithm) is given by the known inverse-time formula, that for the Gaussian updating scheme (commonly used in metadynamics) is often more complex. We further show that the single-bin updating scheme is optimal for very long simulations, and it can be generalized to a class of bandpass updating schemes that are similarly optimal. These bandpass updating schemes target only a few long-range distribution modes and their optimal schedule is also given by the inverse-time formula. Constructed from orthogonal polynomials, the bandpass updating schemes generalize the WL and Langfeld-Lucini-Rago algorithms as an automatic parameter tuning scheme for umbrella sampling.

van der Waals Contribution to the Relative Stability of Aqueous Zn(2+) Coordination States

Many properties of aqueous cations depend on their coordination state. However, the lack of long-range order and the dynamic character of aqueous solutions make it difficult to obtain information beyond average coordination parameters. A thorough understanding of the molecular-scale environment of aqueous cations usually requires a combination of experimental and theoretical approaches. In the case of Zn²⁺, significant discrepancies occur among theoretical investigations based on first-principles molecular dynamics (FPMD) or free-energy calculations, although experimental data consistently point to a dominant hexaaquo-zinc complex (Zn[H₂O]₆)²⁺ in pure water. In the present study, the aqueous speciation of zinc is theoretically investigated by combining FPMD simulations and free-energy calculations based on metadynamics and umbrella-sampling strategies. The simulations are carried out within the density functional theory (DFT) framework using for the exchange-correlation functional either a standard generalized gradient approximation (GGA) or a nonlocal functional (vdw-DF2) which includes van der Waals interactions. The theoretical environment of Zn is confronted to experiment by comparing calculated and measured X-ray absorption spectra. It is shown that the inclusion of van der Waals interactions is crucial for the correct modeling of zinc aqueous speciation, whereas GGA incorrectly favors tetraaquo- (Zn[H₂O]₄)²⁺ and pentaaquo-zinc (Zn[H₂O]₅)²⁺ complexes, results obtained with the vdW-DF2 functional show that the hexaaquo-zinc complex is more stable than the tetraaquo and pentaaquo-zinc complexes by 13 and by 4 kJ mol^-1, respectively. These results highlight the critical importance of even subtle interactions for the correct balance of different coordination states in aqueous solutions. However, for a given coordination state, GGA leads to a reasonable description of the geometry of the aqueous complex.

Free-Energy Landscape of the Amino-Terminal Fragment of Huntingtin in Aqueous Solution

The first exon of Huntingtin-a protein with multiple biological functions whose misfolding is related to Huntington's disease-modulates its localization, aggregation, and function within the cell. It is composed of a 17-amino-acid amphipathic segment (Htt17), an amyloidogenic segment of consecutive glutamines (QN), and a proline-rich segment. Htt17 is of fundamental importance: it serves as a membrane anchor to control the localization of huntingtin, it modulates huntingtin's function through posttranslational modifications, and it controls the self-assembly of the amyloidogenic QN segment into oligomers and fibrils. Experimentally, the conformational ensemble of the Htt17 monomer, as well as the impact of the polyglutamine and proline-rich segments, remains, however, mostly uncharacterized at the atomic level due to its intrinsic flexibility. Here, we unveil the free-energy landscape of Htt17, Htt17Q17, and Htt17Q17P11 using Hamiltonian replica exchange combined with well-tempered metadynamics. We characterize the free-energy landscape of these three fragments in terms of a few selected collective variables. Extensive simulations reveal that the free energy of Htt17 is dominated by a broad ensemble of configurations that agree with solution NMR chemical shifts. Addition of Q17 at its carboxy-terminus reduces the extent of the main basin to more extended configurations of Htt17 with lower helix propensity. Also, the aliphatic carbons of Q17 partially sequester the nonpolar amino acids of Htt17. For its part, addition of Q17P11 shifts the overall landscape to a more extended and helical Htt17 stabilized by interactions with Q17 and P11, which almost exclusively form a PPII-helix, as well as by intramolecular H-bonds and salt bridges. Our characterization of Huntingtin's amino-terminus provides insights into the structural origin of its ability to oligomerize and interact with phospholipid bilayers, processes closely linked to the biological functions of this protein.

Transition-Tempered Metadynamics Is a Promising Tool for Studying the Permeation of Drug-like Molecules through Membranes

Metadynamics is an important enhanced sampling technique in molecular dynamics simulation to efficiently explore potential energy surfaces. The recently developed transition-tempered metadynamics (TTMetaD) has been proven to converge asymptotically without sacrificing exploration of the collective variable space in the early stages of simulations, unlike other convergent metadynamics (MetaD) methods. We have applied TTMetaD to study the permeation of drug-like molecules through a lipid bilayer to further investigate the usefulness of this method as applied to problems of relevance to medicinal chemistry. First, ethanol permeation through a lipid bilayer was studied to compare TTMetaD with nontempered metadynamics and well-tempered metadynamics. The bias energies computed from various metadynamics simulations were compared to the potential of mean force calculated from umbrella sampling. Though all of the MetaD simulations agree with one another asymptotically, TTMetaD is able to predict the most accurate and reliable estimate of the potential of mean force for permeation in the early stages of the simulations and is robust to the choice of required additional parameters. We also show that using multiple randomly initialized replicas allows convergence analysis and also provides an efficient means to converge the simulations in shorter wall times and, more unexpectedly, in shorter CPU times; splitting the CPU time between multiple replicas appears to lead to less overall error. After validating the method, we studied the permeation of a more complicated drug-like molecule, trimethoprim. Three sets of TTMetaD simulations with different choices of collective variables were carried out, and all converged within feasible simulation time. The minimum free energy paths showed that TTMetaD was able to predict almost identical permeation mechanisms in each case despite significantly different definitions of collective variables.

All-Atom Structural Models of the Transmembrane Domains of Insulin and Type 1 Insulin-Like Growth Factor Receptors

The receptor tyrosine kinase superfamily comprises many cell-surface receptors including the insulin receptor (IR) and type 1 insulin-like growth factor receptor (IGF1R) that are constitutively homodimeric transmembrane glycoproteins. Therefore, these receptors require ligand-triggered domain rearrangements rather than receptor dimerization for activation. Specifically, binding of peptide ligands to receptor ectodomains transduces signals across the transmembrane domains for trans-autophosphorylation in cytoplasmic kinase domains. The molecular details of these processes are poorly understood in part due to the absence of structures of full-length receptors. Using MD simulations and enhanced conformational sampling algorithms, we present all-atom structural models of peptides containing 51 residues from the transmembrane and juxtamembrane regions of IR and IGF1R. In our models, the transmembrane regions of both receptors adopt helical conformations with kinks at Pro961 (IR) and Pro941 (IGF1R), but the C-terminal residues corresponding to the juxtamembrane region of each receptor adopt unfolded and flexible conformations in IR as opposed to a helix in IGF1R. We also observe that the N-terminal residues in IR form a kinked-helix sitting at the membrane-solvent interface, while homologous residues in IGF1R are unfolded and flexible. These conformational differences result in a larger tilt-angle of the membrane-embedded helix in IGF1R in comparison to IR to compensate for interactions with water molecules at the membrane-solvent interfaces. Our metastable/stable states for the transmembrane domain of IR, observed in a lipid bilayer, are consistent with a known NMR structure of this domain determined in detergent micelles, and similar states in IGF1R are consistent with a previously reported model of the dimerized transmembrane domains of IGF1R. Our all-atom structural models suggest potentially unique structural organization of kinase domains in each receptor.

Exploring Valleys without Climbing Every Peak: More Efficient and Forgiving Metabasin Metadynamics via Robust On-the-Fly Bias Domain Restriction

Metadynamics is an enhanced sampling method designed to flatten free energy surfaces uniformly. However, the highest-energy regions are often irrelevant to study and dangerous to explore because systems often change irreversibly in unforeseen ways in response to driving forces in these regions, spoiling the sampling. Introducing an on-the-fly domain restriction allows metadynamics to flatten only up to a specified energy level and no further, improving efficiency and safety while decreasing the pressure on practitioners to design collective variables that are robust to otherwise irrelevant high energy driving. This paper describes a new method that achieves this using sequential on-the-fly estimation of energy wells and redefinition of the metadynamics hill shape, termed metabasin metadynamics. The energy level may be defined a priori or relative to unknown barrier energies estimated on-the-fly. Altering only the hill ensures that the method is compatible with many other advances in metadynamics methodology. The hill shape has a natural interpretation in terms of multiscale dynamics, and the computational overhead in simulation is minimal when studying systems of any reasonable size, for instance proteins or other macromolecules. Three example applications show that the formula is accurate and robust to complex dynamics, making metadynamics significantly more forgiving with respect to CV quality and thus more feasible to apply to the most challenging biomolecular systems.

Secondary Structure of Rat and Human Amylin across Force Fields

The aggregation of human amylin has been strongly implicated in the progression of Type II diabetes. This 37-residue peptide forms a variety of secondary structures, including random coils, α-helices, and β-hairpins. The balance between these structures depends on the chemical environment, making amylin an ideal candidate to examine inherent biases in force fields. Rat amylin differs from human amylin by only 6 residues; however, it does not form fibrils. Therefore it provides a useful complement to human amylin in studies of the key events along the aggregation pathway. In this work, the free energy of rat and human amylin was determined as a function of α-helix and β-hairpin content for the Gromos96 53a6, OPLS-AA/L, CHARMM22/CMAP, CHARMM22*, Amberff99sb*-ILDN, and Amberff03w force fields using advanced sampling techniques, specifically bias exchange metadynamics. This work represents a first systematic attempt to evaluate the conformations and the corresponding free energy of a large, clinically relevant disordered peptide in solution across force fields. The NMR chemical shifts of rIAPP were calculated for each of the force fields using their respective free energy maps, allowing us to quantitatively assess their predictions. We show that the predicted distribution of secondary structures is sensitive to the choice of force-field: Gromos53a6 is biased towards β-hairpins, while CHARMM22/CMAP predicts structures that are overly α-helical. OPLS-AA/L favors disordered structures. Amberff99sb*-ILDN, AmberFF03w and CHARMM22* provide the balance between secondary structures that is most consistent with available experimental data. In contrast to previous reports, our findings suggest that the equilibrium conformations of human and rat amylin are remarkably similar, but that subtle differences arise in transient alpha-helical and beta-strand containing structures that the human peptide can more readily adopt. We hypothesize that these transient states enable dynamic pathways that facilitate the formation of aggregates and, eventually, amyloid fibrils.

Atomic-resolution dissection of the energetics and mechanism of isomerization of hydrated ATP-Mg(2+) through the SOMA string method

The atomic mechanisms of isomerization of ATP-Mg(2+) in solution are characterized using the recently developed String Method with Optimal Molecular Alignment (SOMA) and molecular-dynamics simulations. Bias-Exchange Metadynamics simulations are first performed to identify the primary conformers of the ATP-Mg(2+) complex and their connectivity. SOMA is then used to elucidate the minimum free-energy path (MFEP) for each transition, in a 48-dimensional space. Analysis of the per-atom contributions to the global free-energy profiles reveals that the mechanism of these transitions is controlled by the Mg(2+) ion and its coordinating oxygen atoms in the triphosphate moiety, as well as by the ion-hydration shell. Metadynamics simulations in path collective variables based on the MFEP demonstrate these isomerizations proceed across a narrow channel of configurational space, thus validating the premise underlying SOMA. This study provides a roadmap for the examination of conformational changes in biomolecules, based on complementary enhanced-sampling techniques with different strengths. © 2015 Wiley Periodicals, Inc.

Extended Hamiltonian approach to continuous tempering

We introduce an enhanced sampling simulation technique based on continuous tempering, i.e., on continuously varying the temperature of the system under investigation. Our approach is mathematically straightforward, being based on an extended Hamiltonian formulation in which an auxiliary degree of freedom, determining the effective temperature, is coupled to the physical system. The physical system and its temperature evolve continuously in time according to the equations of motion derived from the extended Hamiltonian. Due to the Hamiltonian structure, it is easy to show that a particular subset of the configurations of the extended system is distributed according to the canonical ensemble for the physical system at the correct physical temperature.

A molecular view of DNA-conjugated nanoparticle association energies

Nanoparticles functionalized with short sequences of DNA represent a promising platform for customizable self assembly. Though much recent research has focused on the phase behavior and assembly of these structures, little has been done to precisely characterize the pairwise interaction between particles. Here we present a detailed calculation of the association between DNA-nanoparticle conjugates using 3SPN.2, a coarse-grained model of DNA that accounts for molecular structure and base-pairing. We compare our results to those obtained experimentally using μm sized particles and analyze the free energy surfaces that characterize interparticle hybridization. Next, we study the importance of three-body effects and their impact on particle association and melting. Lastly, we explore the observation by Park et al. [Nature, 451, 553 (2008)] that DNA-nanoparticle crystallization can be inhibited by the deletion of a single nucleotide. Using our model, we suggest that the role of this nucleotide is to disrupt frustration.

Topological obstructions in the way of data-driven collective variables

Nonlinear dimensionality reduction (NLDR) techniques are increasingly used to visualize molecular trajectories and to create data-driven collective variables for enhanced sampling simulations. The success of these methods relies on their ability to identify the essential degrees of freedom characterizing conformational changes. Here, we show that NLDR methods face serious obstacles when the underlying collective variables present periodicities, e.g., arising from proper dihedral angles. As a result, NLDR methods collapse very distant configurations, thus leading to misinterpretations and inefficiencies in enhanced sampling. Here, we identify this largely overlooked problem and discuss possible approaches to overcome it. We also characterize the geometry and topology of conformational changes of alanine dipeptide, a benchmark system for testing new methods to identify collective variables.

Computational structural biology of opioid receptors

The publication of high-resolution structures for all of the opioid receptor subfamilies has unveiled exciting opportunities for mechanistic insight into the molecular mechanisms underlying the biology of nociception, reward, and higher cognitive functions, as well as promises for progress in several clinical areas such as pain management, physiological dependence, addiction, and mood disorders. To turn this promise into novel and improved therapeutic entities, however, this information needs to be supplemented with research strategies that explore the dynamic behavior of the proteins and their interactions with other receptors and ligands in their physiological environment.Here we describe state-of-the-art molecular dynamics computational protocols, based on all-atom and coarse-grained modeling techniques, designed to estimate crucial thermodynamic and kinetic parameters describing the binding of small-molecule ligands and the formation of supramolecular complexes.

Basis function sampling: a new paradigm for material property computation

Wang-Landau sampling, and the associated class of flat histogram simulation methods have been remarkably helpful for calculations of the free energy in a wide variety of physical systems. Practically, convergence of these calculations to a target free energy surface is hampered by reliance on parameters which are unknown a priori. Here, we derive and implement a method built upon orthogonal functions which is fast, parameter-free, and (importantly) geometrically robust. The method is shown to be highly effective in achieving convergence. An important feature of this method is its ability to attain arbitrary levels of description for the free energy. It is thus ideally suited to in silico measurement of elastic moduli and other material properties related to free energy perturbations. We demonstrate the utility of such applications by applying our method to calculate the Frank elastic constants of the Lebwohl-Lasher model of liquid crystals.

Structural determinants in prion protein folding and stability

Prions are responsible for a heterogeneous group of fatal neurodegenerative diseases, involving post-translational modifications of the cellular prion protein. Epidemiological studies on Creutzfeldt-Jakob disease, a prototype prion disorder, show a majority of cases being sporadic, while the remaining occurrences are either genetic or iatrogenic. The molecular mechanisms by which PrP(C) is converted into its pathological isoform have not yet been established. While point mutations and seeds trigger the protein to cross the energy barriers, thus causing genetic and infectious transmissible spongiform encephalopathies, respectively, the mechanism responsible for sporadic forms remains unclear. Since prion diseases are protein-misfolding disorders, we investigated prion protein folding and stability as functions of different milieus. Using spectroscopic techniques and atomistic simulations, we dissected the contribution of major structural determinants, also defining the energy landscape of prion protein. In particular, we elucidated (i) the essential role of the octapeptide region in prion protein folding and stability, (ii) the presence of a very enthalpically stable intermediate in prion-susceptible species, and (iii) the role of the disulfide bridge in prion protein folding.

Well-tempered metadynamics converges asymptotically

Metadynamics is a versatile and capable enhanced sampling method for the computational study of soft matter materials and biomolecular systems. However, over a decade of application and several attempts to give this adaptive umbrella sampling method a firm theoretical grounding prove that a rigorous convergence analysis is elusive. This Letter describes such an analysis, demonstrating that well-tempered metadynamics converges to the final state it was designed to reach and, therefore, that the simple formulas currently used to interpret the final converged state of tempered metadynamics are correct and exact. The results do not rely on any assumption that the collective variable dynamics are effectively Brownian or any idealizations of the hill deposition function; instead, they suggest new, more permissive criteria for the method to be well behaved. The results apply to tempered metadynamics with or without adaptive Gaussians or boundary corrections and whether the bias is stored approximately on a grid or exactly.

Statistical mechanics of the denatured state of a protein using replica-averaged metadynamics

The characterization of denatured states of proteins is challenging because the lack of permanent structure in these states makes it difficult to apply to them standard methods of structural biology. In this work we use all-atom replica-averaged metadynamics (RAM) simulations with NMR chemical shift restraints to determine an ensemble of structures representing an acid-denatured state of the 86-residue protein ACBP. This approach has enabled us to reach convergence in the free energy landscape calculations, obtaining an ensemble of structures in relatively accurate agreement with independent experimental data used for validation. By observing at atomistic resolution the transient formation of native and non-native structures in this acid-denatured state of ACBP, we rationalize the effects of single-point mutations on the folding rate, stability, and transition-state structures of this protein, thus characterizing the role of the unfolded state in determining the folding process.

Molecular Dynamics in the Multicanonical Ensemble: Equivalence of Wang-Landau Sampling, Statistical Temperature Molecular Dynamics, and Metadynamics

We show a direct formal relationship between the Wang-Landau iteration [PRL 86, 2050 (2001)], metadynamics [PNAS 99, 12562 (2002)], and statistical temperature molecular dynamics (STMD) [PRL 97, 050601 (2006)] that are the major work-horses for sampling from generalized ensembles. We demonstrate that STMD, itself derived from the Wang-Landau method, can be made indistinguishable from metadynamics. We also show that Gaussian kernels significantly improve the performance of STMD, highlighting the practical benefits of this improved formal understanding.

Human immunoglobulin E flexes between acutely bent and extended conformations

Crystallographic and solution studies have shown that IgE molecules are acutely bent in their Fc region. Crystal structures reveal the Cɛ2 domain pair folded back onto the Cɛ3-Cɛ4 domains, but is the molecule exclusively bent or can the Cɛ2 domains adopt extended conformations and even 'flip' from one side of the molecule to the other? We report the crystal structure of IgE-Fc captured in a fully extended, symmetrical conformation and show by molecular dynamics, calorimetry, stopped-flow kinetic, surface plasmon resonance (SPR) and Förster resonance energy transfer (FRET) analyses that the antibody can indeed adopt such extended conformations in solution. This diversity of conformational states available to IgE-Fc offers a new perspective on IgE function in allergen recognition, as part of the B-cell receptor and as a therapeutic target in allergic disease.

A boundary correction algorithm for metadynamics in multiple dimensions

Metadynamics is an efficient method for simulation of the free energy of many-particle systems. Over the last several years it has been applied to study a wide variety of systems, ranging from simple fluids to biological macromolecules. The method relies on uniform sampling along specified collective variables or order parameters. Such order parameters, however, are often bounded, and metadynamics algorithms as originally developed suffer from systematic errors at the corresponding boundaries. While several approaches have been proposed in the past to correct these errors for unidimensional systems, no method exists to fully correct these errors in multi-dimensional systems at points where multiple boundaries meet. Here we present a correction scheme that circumvents this limitation.

Computational study of peptide permeation through membrane: Searching for hidden slow variables

Atomically detailed molecular dynamics trajectories in conjunction with Milestoning are used to analyze the different contributions of coarse variables to the permeation process of a small peptide (N-acetyl-L-tryptophanamide, NATA) through a 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) membrane. The peptide reverses its overall orientation as it permeates through the biological bilayer. The large change in orientation is investigated explicitly but is shown to impact the free energy landscape and permeation time only moderately. Nevertheless, a significant difference in permeation properties of the two halves of the membrane suggests the presence of other hidden slow variables. We speculate, based on calculation of the potential of mean force, that a conformational transition of NATA makes significant contribution to these differences. Other candidates for hidden slow variables may include water permeation and collective motions of phospholipids.