Search

Advanced potential energy surfaces are defined as theoretical models that explicitly include many-body effects that transcend the standard fixed-charge, pairwise-additive paradigm typically used in molecular simulation. However, several factors relating to their software implementation have precluded their widespread use in condensed-phase simulations: the computational cost of the theoretical models, a paucity of approximate models and algorithmic improvements that can ameliorate their cost, underdeveloped interfaces and limited dissemination in computational code bases that are widely used in the computational chemistry community, and software implementations that have not kept pace with modern high-performance computing (HPC) architectures, such as multicore CPUs and modern graphics processing units (GPUs). In this Feature Article we review recent progress made in these areas, including well-defined polarization approximations and new multipole electrostatic formulations, novel methods for solving the mutual polarization equations and increasing the MD time step, combining linear-scaling electronic structure methods with new QM/MM methods that account for mutual polarization between the two regions, and the greatly improved software deployment of these models and methods onto GPU and CPU hardware platforms. We have now approached an era where multipole-based polarizable force fields can be routinely used to obtain computational results comparable to state-of-the-art density functional theory while reaching sampling statistics that are acceptable when compared to that obtained from simpler fixed partial charge force fields. [ABSTRACT FROM AUTHOR]

Cluster of differentiation 1c (CD1c)-dependent self-reactive T cells are abundant in human blood, but self-antigens presented by CD1c to the T-cell receptors of these cells are poorly understood. Here we present a crystal structure of CD1c determined at 2.4 Å revealing an extended ligand binding potential of the antigen groove and a substantially different conformation compared with known CD1c structures. Computational simulations exploring different occupancy states of the groove reenacted these different CD1c conformations and suggested cholesteryl esters (CE) and acylated steryl glycosides (ASG) as new ligand classes for CD1c. Confirming this, we show that binding of CE and ASG to CD1c enables the binding of human CD1c self-reactive T-cell receptors. Hence, human CD1c adopts different conformations dependent on ligand occupancy of its groove, with CE and ASG stabilizing CD1c conformations that provide a footprint for binding of CD1c self-reactive T-cell receptors. [ABSTRACT FROM AUTHOR]

Hydrogenases use complexmetal cofactors to catalyze the reversible formation of hydrogen. In [FeFe]-hydrogenases, the H-cluster cofactor includes a diiron subcluster containing azadithiolate, three CO, and two CN- ligands. During the assembly of the H cluster, the radical S-adenosyl methionine (SAM) enzyme HydG lyses the substrate tyrosine to yield the diatomic ligands. These diatomic products form an enzyme-bound Fe(CO)x(CN)y synthon that serves as a precursor for eventual H-cluster assembly. To further elucidate the mechanism of this complex reaction, we report the crystal structure and EPR analysis of HydG. At one end of the HydG (βα)8 triosephosphate isomerase (TIM) barrel, a canonical [4Fe-4S] cluster binds SAM in close proximity to the proposed tyrosine binding site. At the opposite end of the active-site cavity, the structure reveals the auxiliary Fe-S cluster in two states: one monomer contains a [4Fe-5S] cluster, and the other monomer contains a [5Fe-5S] cluster consisting of a [4Fe-4S] cubane bridged by a μ2-sulfide ion to a mononuclear Fe2+ center. This fifth iron is held in place by a single highly conserved protein-derived ligand: histidine 265. EPR analysis confirms the presence of the [5Fe-5S] cluster, which on incubation with cyanide, undergoes loss of the labile iron to yield a [4Fe-4S] cluster. We hypothesize that the labile iron of the [5Fe-5S] cluster is the site of Fe (CO)x(CN)y synthon formation and that the limited bonding between this iron and HydG may facilitate transfer of the intact synthon to its cognate acceptor for subsequent H-cluster assembly. [ABSTRACT FROM AUTHOR]

Water molecules play a crucial role in mediating the interaction between a ligand and a macromolecular receptor. An understanding of the nature and role of each water molecule in the active site of a protein could greatly increase the efficiency of rational drug design approaches: if the propensity of a water molecule for displacement can be determined, then synthetic effort may be most profitably applied to the design of specific ligands with the displacement of this water molecule in mind. In this paper, a thermodynamic analysis of water molecules in the binding sites of six proteins, each complexed with a number of inhibitors, is presented. Two classes of water molecules were identified: those conserved and not displaced by any of the ligands, and those that are displaced by some ligands. The absolute binding free energies of 54 water molecules were calculated using the double decoupling method, with replica exchange thermodynamic integration in Monte Carlo simulations. It was found that conserved water molecules are on average more tightly bound than displaced water molecules. In addition, Bayesian statistics is used to calculate the probability that a particular water molecule may be displaced by an appropriately designed ligand, given the calculated binding free energy of the water molecule. This approach therefore allows the numerical assessment of whether or not a given water molecule should be targeted for displacement as part of a rational drug design strategy. [ABSTRACT FROM AUTHOR]

Abstract: Comparative molecular dynamics (MD) simulations enable us to explore the conformational dynamics of the active sites of distantly related enzymes. We have used the BioSimGrid (http://www.biosimgrid.org) database to facilitate such a comparison. Simulations of four enzymes were analyzed. These included three hydrolases and a transferase, namely acetylcholinesterase, outer-membrane phospholipase A, outer-membrane protease T, and PagP (an outer-membrane enzyme which transfers a palmitate chain from a phospholipid to lipid A). A set of 17 simulations were analyzed corresponding to a total of ∼0.1μs simulation time. A simple metric for active-site integrity was used to demonstrate the existence of clusters of dynamic conformational behaviour of the active sites. Small (i.e. within a cluster) fluctuations appear to be related to the function of an enzymatically active site. Larger fluctuations (i.e. between clusters) correlate with transitions between catalytically active and inactive states. Overall, these results demonstrate the potential of a comparative MD approach to analysis of enzyme function. This approach could be extended to a wider range of enzymes using current high throughput MD simulation and database methods. [Copyright &y& Elsevier]

Abstract: A library of polyamine–peptide conjugates based around some previously identified inhibitors of trypanothione reductase was synthesised by parallel solid-phase chemistry and screened. Kinetic analysis of library members established that subtle structural changes altered their mechanism of action, switching between competitive and non-competitive inhibition. The mode of action of the non-competitive inhibitors was investigated in detail by a variety of techniques including enzyme kinetic analysis (looking at both NADPH and trypanothione disulfide substrates), gel filtration chromatography and analytical ultracentrifugation, leading to the identification of an allosteric mode of inhibition. [Copyright &y& Elsevier]

Investigates the synthesis and characterization of highly selective fluoride receptor and its binding properties. Use of Monte Carlo free energy perturbation technique and Poisson calculations; Effect of water on the calculated free energies of binding; Electrostatic interactions that the fluoride gains by sitting lower in the phenolic activity of the receptor.

XPO1 (Exportin-1/CRM1) is a nuclear export protein that is frequently overexpressed in cancer and functions as a driver of oncogenesis. Currently small molecules that target XPO1 are being used in the clinic as anticancer agents. We identify XPO1 as a target for natural killer (NK) cells. Using immunopeptidomics, we have identified a peptide derived from XPO1 that can be recognized by the activating NK cell receptor KIR2DS2 in the context of human leukocyte antigen-C. The peptide can be endogenously processed and presented to activate NK cells specifically through this receptor. Although high XPO1 expression in cancer is commonly associated with a poor prognosis, we show that the outcome of specific cancers, such as hepatocellular carcinoma, can be substantially improved if there is concomitant evidence of NK cell infiltration. We thus identify XPO1 as a bona fide tumor antigen recognized by NK cells that offers an opportunity for a personalized approach to NK cell therapy for solid tumors.

Fluorine substitution can have a profound impact on molecular conformation. Here, we present a detailed conformational analysis of how the 1,3-difluoropropylene motif (-CHF-CH 2 -CHF-) determines the conformational profiles of 1,3-difluoropropane, anti - and syn -2,4-difluoropentane, and anti - and syn -3,5-difluoroheptane. It is shown that the 1,3-difluoropropylene motif strongly influences alkane chain conformation, with a significant dependence on the polarity of the medium. The conformational effect of 1,3-fluorination is magnified upon chain extension, which contrasts with vicinal difluorination. Experimental evidence was obtained from NMR analysis, where polynomial complexity scaling simulation algorithms were necessary to enable J -coupling extraction from the strong second-order spectra, particularly for the large 16-spin systems of the difluorinated heptanes. These results improve our understanding of the conformational control toolkit for aliphatic chains, yield simple rules for conformation population analysis, and demonstrate quantum mechanical time-domain NMR simulations for liquid state systems with large numbers of strongly coupled spins.

Enhanced sampling algorithms are indispensable when working with highly disconnected multimodal distributions. An important application of these is the conformational exploration of particular internal degrees of freedom of molecular systems. However, despite the existence of many commonly used enhanced sampling algorithms to explore these internal motions, they often rely on system-dependent parameters, which negatively impact efficiency and reproducibility. Here, we present fully adaptive simulated tempering (FAST), a variation of the irreversible simulated tempering algorithm, which continuously optimizes the number, parameters, and weights of intermediate distributions to achieve maximally fast traversal over a space defined by the change in a predefined thermodynamic control variable such as temperature or an alchemical smoothing parameter. This work builds on a number of previously published methods, such as sequential Monte Carlo, and introduces a novel parameter optimization procedure that can, in principle, be used in any expanded ensemble algorithms. This method is validated by being applied on a number of different molecular systems with high torsional kinetic barriers. We also consider two different soft-core potentials during the interpolation procedure and compare their performance. We conclude that FAST is a highly efficient algorithm, which improves simulation reproducibility and can be successfully used in a variety of settings with the same initial hyperparameters., (© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).)

Hypoxia inducible factor (HIF) is a heterodimeric transcription factor composed of an oxygen-regulated α subunit and a constitutively expressed β subunit that serves as the master regulator of the cellular response to low oxygen concentrations. The HIF transcription factor senses and responds to hypoxia by significantly altering transcription and reprogramming cells to enable adaptation to a hypoxic microenvironment. Given the central role played by HIF in the survival and growth of tumors in hypoxia, inhibition of this transcription factor serves as a potential therapeutic approach for treating a variety of cancers. Here, we report the identification, optimization, and characterization of a series of cyclic peptides that disrupt the function of HIF-1 and HIF-2 transcription factors by inhibiting the interaction of both HIF-1α and HIF-2α with HIF-1β. These compounds are shown to bind to HIF-α and disrupt the protein-protein interaction between the α and β subunits of the transcription factor, resulting in disruption of hypoxia-response signaling by our lead molecule in several cancer cell lines.

We present a comparative study that evaluates the performance of a machine learning potential (ANI-2x), a conventional force field (GAFF), and an optimally tuned GAFF-like force field in the modeling of a set of 10 γ-fluorohydrins that exhibit a complex interplay between intra- and intermolecular interactions in determining conformer stability. To benchmark the performance of each molecular model, we evaluated their energetic, geometric, and sampling accuracies relative to quantum-mechanical data. This benchmark involved conformational analysis both in the gas phase and chloroform solution. We also assessed the performance of the aforementioned molecular models in estimating nuclear spin-spin coupling constants by comparing their predictions to experimental data available in chloroform. The results and discussion presented in this study demonstrate that ANI-2x tends to predict stronger-than-expected hydrogen bonding and overstabilize global minima and shows problems related to inadequate description of dispersion interactions. Furthermore, while ANI-2x is a viable model for modeling in the gas phase, conventional force fields still play an important role, especially for condensed-phase simulations. Overall, this study highlights the strengths and weaknesses of each model, providing guidelines for the use and future development of force fields and machine learning potentials.

Antigen processing is an immunological mechanism by which intracellular peptides are transported to the cell surface while bound to Major Histocompatibility Complex molecules, where they can be surveyed by circulating CD8+ or CD4+ T-cells, potentially triggering an immunological response. The antigen processing pathway is a complex multistage filter that refines a huge pool of potential peptide ligands derived from protein degradation into a smaller ensemble for surface presentation. Each stage presents unique challenges due to the number of ligands, the polymorphic nature of MHC and other protein constituents of the pathway and the nature of the interactions between them. Predicting the ensemble of displayed peptide antigens, as well as their immunogenicity, is critical for improving T cell vaccines against pathogens and cancer. Our predictive abilities have always been hindered by an incomplete empirical understanding of the antigen processing pathway. In this review, we highlight the role of computational and structural approaches in improving our understanding of antigen processing, including structural biology, computer simulation, and machine learning techniques, with a particular focus on the MHC-I pathway., (© 2023 The Author(s). Published by Portland Press Limited on behalf of the Biochemical Society.)

Water molecules play a key role in many biomolecular systems, particularly when bound at protein-ligand interfaces. However, molecular simulation studies on such systems are hampered by the relatively long time scales over which water exchange between a protein and solvent takes place. Grand canonical Monte Carlo (GCMC) is a simulation technique that avoids this issue by attempting the insertion and deletion of water molecules within a given structure. The approach is constrained by low acceptance probabilities for insertions in congested systems, however. To address this issue, here, we combine GCMC with nonequilibium candidate Monte Carlo (NCMC) to yield a method that we refer to as grand canonical nonequilibrium candidate Monte Carlo (GCNCMC), in which the water insertions and deletions are carried out in a gradual, nonequilibrium fashion. We validate this new approach by comparing GCNCMC and GCMC simulations of bulk water and three protein binding sites. We find that not only is the efficiency of the water sampling improved by GCNCMC but that it also results in increased sampling of ligand conformations in a protein binding site, revealing new water-mediated ligand-binding geometries that are not observed using alternative enhanced sampling techniques.

Antibody responses during infection and vaccination typically undergo affinity maturation to achieve high-affinity binding for efficient neutralization of pathogens 1,2 . Similarly, high affinity is routinely the goal for therapeutic antibody generation. However, in contrast to naturally occurring or direct-targeting therapeutic antibodies, immunomodulatory antibodies, which are designed to modulate receptor signalling, have not been widely examined for their affinity-function relationship. Here we examine three separate immunologically important receptors spanning two receptor superfamilies: CD40, 4-1BB and PD-1. We show that low rather than high affinity delivers greater activity through increased clustering. This approach delivered higher immune cell activation, in vivo T cell expansion and antitumour activity in the case of CD40. Moreover, an inert anti-4-1BB monoclonal antibody was transformed into an agonist. Low-affinity variants of the clinically important antagonistic anti-PD-1 monoclonal antibody nivolumab also mediated more potent signalling and affected T cell activation. These findings reveal a new paradigm for augmenting agonism across diverse receptor families and shed light on the mechanism of antibody-mediated receptor signalling. Such affinity engineering offers a rational, efficient and highly tuneable solution to deliver antibody-mediated receptor activity across a range of potencies suitable for translation to the treatment of human disease., (© 2023. The Author(s), under exclusive licence to Springer Nature Limited.)

Water molecules at protein-ligand interfaces are often of significant pharmaceutical interest, owing in part to the entropy which can be released upon the displacement of an ordered water by a therapeutic compound. Protein structures may not, however, completely resolve all critical bound water molecules, or there may be no experimental data available. As such, predicting the location of water molecules in the absence of a crystal structure is important in the context of rational drug design. Grand canonical Monte Carlo (GCMC) is a computational technique that is gaining popularity for the simulation of buried water sites. In this work, we assess the ability of GCMC to accurately predict water binding locations, using a dataset that we have curated, containing 108 unique structures of complexes between proteins and Food and Drug Administration (FDA)-approved small-molecule drugs. We show that GCMC correctly predicts 81.4% of nonbulk crystallographic water sites to within 1.4 Å. However, our analysis demonstrates that the reported performance of water prediction methods is highly sensitive to the way in which the performance is measured. We also find that crystallographic water sites with more protein/ligand hydrogen bonds and stronger electron density are more reliably predicted by GCMC. An analysis of water networks revealed that more than half of the structures contain at least one ligand-contacting water network. In these cases, displacement of a water site by a ligand modification might yield unexpected results if the larger network is destabilized. Cooperative effects between waters should therefore be explicitly considered in structure-based drug design.

Antibodies protect from infection, underpin successful vaccines and elicit therapeutic responses in otherwise untreatable cancers and autoimmune conditions. The human IgG2 isotype displays a unique capacity to undergo disulfide shuffling in the hinge region, leading to modulation of its ability to drive target receptor signaling (agonism) in a variety of important immune receptors, through hitherto unexplained molecular mechanisms. To address the underlying process and reveal how hinge disulfide orientation affects agonistic activity, we generated a series of cysteine to serine exchange variants in the hinge region of the clinically relevant monoclonal antibody ChiLob7/4, directed against the key immune receptor CD40. We report how agonistic activity varies with disulfide pattern and is afforded by the presence of a disulfide crossover between F(ab) arms in the agonistic forms, independently of epitope, as observed in the determined crystallographic structures. This structural "switch" affects directly on antibody conformation and flexibility. Small-angle x-ray scattering and ensemble modeling demonstrated that the least flexible variants adopt the fewest conformations and evoke the highest levels of receptor agonism. This covalent change may be amenable for broad implementation to modulate receptor signaling in an epitope-independent manner in future therapeutics.

The sampling problem is one of the most widely studied topics in computational chemistry. While various methods exist for sampling along a set of reaction coordinates, many require system-dependent hyperparameters to achieve maximum efficiency. In this work, we present an alchemical variation of adaptive sequential Monte Carlo (SMC), an irreversible importance resampling method that is part of a well-studied class of methods that have been used in various applications but have been underexplored in computational biophysics. Afterward, we apply alchemical SMC on a variety of test cases, including torsional rotations of solvated ligands (butene and a terphenyl derivative), translational and rotational movements of protein-bound ligands, and protein side chain rotation coupled to the ligand degrees of freedom (T4-lysozyme, protein tyrosine phosphatase 1B, and transforming growth factor β). We find that alchemical SMC is an efficient way to explore targeted degrees of freedom and can be applied to a variety of systems using the same hyperparameters to achieve a similar performance. Alchemical SMC is a promising tool for preparatory exploration of systems where long-timescale sampling of the entire system can be traded off against short-timescale sampling of a particular set of degrees of freedom over a population of conformers.

We have carried out a long-timescale simulation study on crystal structures of nine antibody-antigen pairs, in antigen-bound and antibody-only forms, using molecular dynamics with enhanced sampling and an explicit water model to explore interface conformation and hydration. By combining atomic level simulation and replica exchange to enable full protein flexibility, we find significant numbers of bridging water molecules at the antibody-antigen interface. Additionally, a higher proportion of interactions excluding bulk waters and a lower degree of antigen bound CDR conformational sampling are correlated with higher antibody affinity. The CDR sampling supports enthalpically driven antibody binding, as opposed to entropically driven, in that the difference between antigen bound and unbound conformations do not correlate with affinity. We thus propose that interactions with waters and CDR sampling are aspects of the interface that may moderate antibody-antigen binding, and that explicit hydration and CDR flexibility should be considered to improve antibody affinity prediction and computational design workflows., Competing Interests: MW’s studentship was partially funded by UCB. MW, SK, and RT are employees of UCB. XL and TB were former employees of UCB, and XL is a current employee of AngitiaBio. JE’s research is partially funded by UCB., (Copyright © 2022 Wong, Kelm, Liu, Taylor, Baker and Essex.)

Water molecules play important roles in all biochemical processes. Therefore, it is of key importance to obtain information of the structure, dynamics, and thermodynamics of water molecules around proteins. Numerous computational methods have been suggested with this aim. In this study, we compare the performance of conventional and grand-canonical Monte Carlo (GCMC) molecular dynamics (MD) simulations to sample the water structure, as well GCMC and grid-based inhomogeneous solvation theory (GIST) to describe the energetics of the water network. They are evaluated on two proteins: the buried ligand-binding site of a ferritin dimer and the solvent-exposed binding site of galectin-3. We show that GCMC/MD simulations significantly speed up the sampling and equilibration of water molecules in the buried binding site, thereby making the results more similar for simulations started from different states. Both GCMC/MD and conventional MD reproduce crystal-water molecules reasonably for the buried binding site. GIST analyses are normally based on restrained MD simulations. This improves the precision of the calculated energies, but the restraints also significantly affect both absolute and relative energies. Solvation free energies for individual water molecules calculated with and without restraints show a good correlation, but with large quantitative differences. Finally, we note that the solvation free energies calculated with GIST are ∼5 times larger than those estimated by GCMC owing to differences in the reference state., Competing Interests: The authors declare no competing financial interest., (© 2022 The Authors. Published by American Chemical Society.)

Water often plays a key role in protein structure, molecular recognition, and mediating protein-ligand interactions. Thus, free energy calculations must adequately sample water motions, which often proves challenging in typical MD simulation time scales. Thus, the accuracy of methods relying on MD simulations ends up limited by slow water sampling. Particularly, as a ligand is removed or modified, bulk water may not have time to fill or rearrange in the binding site. In this work, we focus on several molecular dynamics (MD) simulation-based methods attempting to help rehydrate buried water sites: BLUES, using nonequilibrium candidate Monte Carlo (NCMC); grand , using grand canonical Monte Carlo (GCMC); and normal MD. We assess the accuracy and efficiency of these methods in rehydrating target water sites. We selected a range of systems with varying numbers of waters in the binding site, as well as those where water occupancy is coupled to the identity or binding mode of the ligand. We analyzed the rehydration of buried water sites in binding pockets using both clustering of trajectories and direct analysis of electron density maps. Our results suggest both BLUES and grand enhance water sampling relative to normal MD and grand is more robust than BLUES, but also that water sampling remains a major challenge for all of the methods tested. The lessons we learned for these methods and systems are discussed.

Understanding the conformational ensembles of intrinsically disordered proteins and peptides (IDPs) in their various biological environments is essential for understanding their mechanisms and functional roles in the proteome, leading to a greater knowledge of, and potential treatments for, a broad range of diseases. To determine whether molecular simulation is able to generate accurate conformational ensembles of IDPs, we explore the structural landscape of the PLP peptide (an intrinsically disordered region of the proteolipid membrane protein) in aqueous and membrane-mimicking solvents, using replica exchange with solute scaling (REST2), and examine the ability of four force fields (ff14SB, ff14IDPSFF, CHARMM36 and CHARMM36m) to reproduce literature circular dichroism (CD) data. Results from variable temperature (VT) 1 H and Rotating frame Overhauser Effect SpectroscopY (ROESY) nuclear magnetic resonance (NMR) experiments are also presented and are consistent with the structural observations obtained from the simulations and CD. We also apply the optimum simulation protocol to TP2 and ONEG (a cell-penetrating peptide (CPP) and a negative control peptide, respectively) to gain insight into the structural differences that may account for the observed difference in their membrane-penetrating abilities. Of the tested force fields, we find that CHARMM36 and CHARMM36m are best suited to the study of IDPs, and accurately predict a disordered to helical conformational transition of the PLP peptide accompanying the change from aqueous to membrane-mimicking solvents. We also identify an α-helical structure of TP2 in the membrane-mimicking solvents and provide a discussion of the mechanistic implications of this observation with reference to the previous literature on the peptide. From these results, we recommend the use of CHARMM36m with the REST2 protocol for the study of environment-specific IDP conformations. We believe that the simulation protocol will allow the study of a broad range of IDPs that undergo conformational transitions in different biological environments., Competing Interests: JWE's research is part funded by AstraZeneca. LMR is currently employed by MedChemica Ltd., (This journal is © The Royal Society of Chemistry.)

Atomistic models provide a detailed representation of molecular systems, but are sometimes inadequate for simulations of large systems over long timescales. Coarse-grained models enable accelerated simulations by reducing the number of degrees of freedom, at the cost of reduced accuracy. New optimisation processes to parameterise these models could improve their quality and range of applicability. We present an automated approach for the optimisation of coarse-grained force fields, by reproducing free energy data derived from atomistic molecular simulations. To illustrate the approach, we implemented hydration free energy gradients as a new target for force field optimisation in ForceBalance and applied it successfully to optimise the un-charged side-chains and the protein backbone in the SIRAH protein coarse-grain force field. The optimised parameters closely reproduced hydration free energies of atomistic models and gave improved agreement with experiment.

Conformational analysis is of paramount importance in drug design: it is crucial to determine pharmacological properties, understand molecular recognition processes, and characterize the conformations of ligands when unbound. Molecular Mechanics (MM) simulation methods, such as Monte Carlo (MC) and molecular dynamics (MD), are usually employed to generate ensembles of structures due to their ability to extensively sample the conformational space of molecules. The accuracy of these MM-based schemes strongly depends on the functional form of the force field (FF) and its parametrization, components that often hinder their performance. High-level methods, such as ab initio MD, provide reliable structural information but are still too computationally expensive to allow for extensive sampling. Therefore, to overcome these limitations, we present a multilevel MC method that is capable of generating quantum configurational ensembles while keeping the computational cost at a minimum. We show that FF reparametrization is an efficient route to generate FFs that reproduce QM results more closely, which, in turn, can be used as low-cost models to achieve the gold standard QM accuracy. We demonstrate that the MC acceptance rate is strongly correlated with various phase space overlap measurements and that it constitutes a robust metric to evaluate the similarity between the MM and QM levels of theory. As a more advanced application, we present a self-parametrizing version of the algorithm, which combines sampling and FF parametrization in one scheme, and apply the methodology to generate the QM/MM distribution of a ligand in aqueous solution.

The influenza A M2 wild-type (WT) proton channel is the target of the anti-influenza drug rimantadine. Rimantadine has two enantiomers, though most investigations into drug binding and inhibition have used a racemic mixture. Solid-state NMR experiments using the full length-M2 WT have shown significant spectral differences that were interpreted to indicate tighter binding for ( R )- vs ( S )-rimantadine. However, it was unclear if this correlates with a functional difference in drug binding and inhibition. Using X-ray crystallography, we have determined that both ( R )- and ( S )-rimantadine bind to the M2 WT pore with slight differences in the hydration of each enantiomer. However, this does not result in a difference in potency or binding kinetics, as shown by similar values for k on , k off , and K d in electrophysiological assays and for EC 50 values in cellular assays. We concluded that the slight differences in hydration for the ( R )- and ( S )-rimantadine enantiomers are not relevant to drug binding or channel inhibition. To further explore the effect of the hydration of the M2 pore on binding affinity, the water structure was evaluated by grand canonical ensemble molecular dynamics simulations as a function of the chemical potential of the water. Initially, the two layers of ordered water molecules between the bound drug and the channel's gating His37 residues mask the drug's chirality. As the chemical potential becomes more unfavorable, the drug translocates down to the lower water layer, and the interaction becomes more sensitive to chirality. These studies suggest the feasibility of displacing the upper water layer and specifically recognizing the lower water layers in novel drugs.

Binding free energy calculations using alchemical free energy (AFE) methods are widely considered to be the most rigorous tool in the computational drug discovery arsenal. Despite this, the calculations suffer from accuracy, precision, and reproducibility issues. In this publication, we perform a high-throughput study of more than a thousand AFE calculations, utilizing over 220 μs of total sampling time, on three different protein systems to investigate the impact of the initial crystal structure on the resulting binding free energy values. We also consider the influence of equilibration time and discover that the initial crystal structure can have a significant effect on free energy values obtained at short timescales that can manifest itself as a free energy difference of more than 1 kcal/mol. At longer timescales, these differences are largely overtaken by important rare events, such as torsional ligand motions, typically resulting in a much higher uncertainty in the obtained values. This work emphasizes the importance of rare event sampling and long-timescale dynamics in free energy calculations even for routinely performed alchemical perturbations. We conclude that an optimal protocol should not only concentrate computational resources on achieving convergence in the alchemical coupling parameter (λ) space but also on longer simulations and multiple repeats.

Current computer architectures, coupled with state-of-the-art molecular dynamics simulation software, facilitate the in-depth study of large biomolecular systems at high levels of detail. However, biological phenomena take place at various time and length scales and as a result a multiscale approach must be adopted. One such approach is coarse-graining, where biochemical accuracy is sacrificed for computational efficiency. Here, we present a practical guide to setting up and carrying out coarse-grained molecular dynamics simulations.

Networks of water molecules can play a critical role at the protein-ligand interface and can directly influence drug-target interactions. Grand canonical methods aid in the sampling of these water molecules, where conventional molecular dynamics equilibration times are often long, by allowing waters to be inserted and deleted from the system, according to the chemical potential. Here, we present our open source Python module, grand (https://github.com/essex-lab/grand), which allows molecular dynamics simulations to be performed in conjunction with grand canonical Monte Carlo sampling, using the OpenMM simulation engine. We demonstrate the accuracy of this module by reproducing the density of bulk water observed from constant pressure simulations. Application of this code to the bovine pancreatic trypsin inhibitor protein reproduces three buried crystallographic water sites that are poorly sampled using conventional molecular dynamics.

Preorganization of large, directionally oriented, electric fields inside protein active sites has been proposed as a crucial contributor to catalytic mechanism in many enzymes, and it may be efficiently investigated at the atomistic level with molecular dynamics simulations. Here, we evaluate the ability of the AMOEBA polarizable force field, as well as the additive Amber ff14SB and Charmm C36m models, to describe the electric fields present inside the active site of the peptidyl-prolyl isomerase cyclophilin A. We compare the molecular mechanical electric fields to those calculated with a fully first-principles quantum mechanical (QM) representation of the protein, solvent, and ions, and find that AMOEBA consistently shows far greater correlation with the QM electric fields than either of the additive force fields tested. Catalytically relevant fields calculated with AMOEBA were typically smaller than those observed with additive potentials, but were generally consistent with an electrostatically driven mechanism for catalysis. Our results highlight the accuracy and the potential advantages of using polarizable force fields in systems where accurate electrostatics may be crucial for providing mechanistic insights.

ProtoCaller is a Python library distributed through Anaconda which automates relative protein-ligand binding free energy calculations in GROMACS. It links a number of popular specialized tools used to perform protein setup and parametrization, such as PDB2PQR, Modeller, and AmberTools. ProtoCaller supports commonly used AMBER force fields with additional cofactor parameters, and AM1-BCC is used to derive ligand charges. ProtoCaller also comes with an extensive PDB parser, an enhanced maximum common substructure algorithm providing improved ligand-ligand mapping, and a light GROMACS wrapper for running multiple molecular dynamics simulations. ProtoCaller is highly relevant to most researchers in the field of biomolecular simulation, allowing a customizable balance between automation and user intervention.

The intrinsic sampling method (ISM) is a powerful tool that allows the exploration of interfacial properties from molecular simulations by fitting a function that represents the local boundary between two phases. However, owing to the non-physical nature of an "intrinsic" surface, there remains an ambiguity surrounding the comparison of theoretical properties with the physical world. It is therefore important that the ISM remains internally consistent when reproducing simulated properties which match experiments, such as the surface tension or interfacial density distribution. We show that the current ISM procedure causes an over-fitting of the surface to molecules in the interface region, leading to a biased distribution of curvature at these molecular coordinates. We assert that this biased distribution is a cause of the disparity between predicted interfacial densities upon convolution to a laboratory frame, an artefact which has been known to exist since the development of the ISM. We present an improvement to the fitting procedure of the ISM in an attempt to alleviate the ambiguity surrounding the true nature of an intrinsic surface. Our "surface reconstruction" method is able to amend the shape of the interface so as to reproduce the global curvature distribution at all sampled molecular coordinates. We present the effects that this method has on the ISM predicted structure of a simulated Lennard-Jones fluid air-liquid interface. Additionally, we report an unexpected relationship between surface thermodynamic predictions of our reconstructed ISM surfaces and those of extended capillary wave theory, which is of current interest.

Computational methods to calculate ligand binding affinities are increasing in popularity, due to improvements in simulation algorithms, computational resources, and easy-to-use software. However, issues can arise in relative ligand binding free energy simulations if the ligands considered have different active site water networks, as simulations are typically performed with a predetermined number of water molecules (fixed N ensembles) in preassigned locations. If an alchemical perturbation is attempted where the change should result in a different active site water network, the water molecules may not be able to adapt appropriately within the time scales of the simulations-particularly if the active site is occluded. By combining the grand canonical ensemble (μVT) to sample active site water molecules, with conventional alchemical free energy methods, the water network is able to dynamically adapt to the changing ligand. We refer to this approach as grand canonical alchemical perturbation (GCAP). In this work we demonstrate GCAP for two systems; Scytalone Dehydratase (SD) and Adenosine A 2 A receptor. For both systems, GCAP is shown to perform well at reproducing experimental binding affinities. Calculating the relative binding affinities with a naı̈ve, conventional attempt to solvate the active site illustrates how poor results can be if proper consideration of water molecules in occluded pockets is neglected. GCAP results are shown to be consistent with time-consuming double decoupling simulations. In addition, by obtaining the free energy surface for ligand perturbations, as a function of both the free energy coupling parameter and water chemical potential, it is possible to directly deconvolute the binding energetics in terms of protein-ligand direct interactions and protein binding site hydration.

Covalent modification of outer membrane lipids of Gram-negative bacteria can impact the ability of the bacterium to develop resistance to antibiotics as well as modulating the immune response of the host. The enzyme LpxR from Salmonella typhimurium is known to deacylate lipopolysaccharide molecules of the outer membrane; however, the mechanism of action is unknown. Here, we employ molecular dynamics and Monte Carlo simulations to study the conformational dynamics and substrate binding of LpxR in representative outer membrane models as well as detergent micelles. We examine the roles of conserved residues and provide an understanding of how LpxR binds its substrate. Our simulations predict that the catalytic H122 must be Nε-protonated for a single water molecule to occupy the space between it and the scissile bond, with a free binding energy of -8.5 kcal mol -1 . Furthermore, simulations of the protein within a micelle enable us to predict the structure of the putative "closed" protein. Our results highlight the need for including dynamics, a representative environment, and the consideration of multiple tautomeric and rotameric states of key residues in mechanistic studies; static structures alone do not tell the full story., (Copyright © 2018 The Authors. Published by Elsevier Inc. All rights reserved.)

Complementary strategies of small-angle x-ray scattering (SAXS) and crystallographic analysis are often used to determine atomistic three-dimensional models of macromolecules and their variability in solution. This combination of techniques is particularly valuable when applied to macromolecular complexes to detect changes within the individual binding partners. Here, we determine the x-ray crystallographic structure of a F(ab) fragment in complex with CD32b, the only inhibitory Fc-γ receptor in humans, and compare the structure of the F(ab) from the crystal complex to SAXS data for the F(ab) alone in solution. We investigate changes in F(ab) structure by predicting theoretical scattering profiles for atomistic structures extracted from molecular dynamics (MD) simulations of the F(ab) and assessing the agreement of these structures to our experimental SAXS data. Through principal component analysis, we are able to extract principal motions observed during the MD trajectory and evaluate the influence of these motions on the agreement of structures to the F(ab) SAXS data. Changes in the F(ab) elbow angle were found to be important to reach agreement with the experimental data; however, further discrepancies were apparent between our F(ab) structure from the crystal complex and SAXS data. By analyzing multiple MD structures observed in similar regions of the principal component analysis, we were able to pinpoint these discrepancies to a specific loop region in the F(ab) heavy chain. This method, therefore, not only allows determination of global changes but also allows identification of localized motions important for determining the agreement between atomistic structures and SAXS data. In this particular case, the findings allowed us to discount the hypothesis that structural changes were induced upon complex formation, a significant find informing the drug development process. The methodology described here is generally applicable to deconvolute global and local changes of macromolecular structures and is well suited to other systems., (Copyright © 2018 Biophysical Society. Published by Elsevier Inc. All rights reserved.)

We present an unexpected finite size effect affecting interfacial molecular simulations that is proportional to the width-to-surface-area ratio of the bulk phase L l /A. This finite size effect has a significant impact on the variance of surface tension values calculated using the virial summation method. A theoretical derivation of the origin of the effect is proposed, giving a new insight into the importance of optimising system dimensions in interfacial simulations. We demonstrate the consequences of this finite size effect via a new way to estimate the surface energetic and entropic properties of simulated air-liquid interfaces. Our method is based on macroscopic thermodynamic theory and involves comparing the internal energies of systems with varying dimensions. We present the testing of these methods using simulations of the TIP4P/2005 water forcefield and a Lennard-Jones fluid model of argon. Finally, we provide suggestions of additional situations, in which this finite size effect is expected to be significant, as well as possible ways to avoid its impact.

Tuberculosis (TB), caused by Mycobacterium tuberculosis , remains a major human pandemic. Germline-encoded mycolyl lipid-reactive (GEM) T cells are donor-unrestricted and recognize CD1b-presented mycobacterial mycolates. However, the molecular requirements governing mycolate antigenicity for the GEM T cell receptor (TCR) remain poorly understood. Here, we demonstrate CD1b expression in TB granulomas and reveal a central role for meromycolate chains in influencing GEM-TCR activity. Meromycolate fine structure influences T cell responses in TB-exposed individuals, and meromycolate alterations modulate functional responses by GEM-TCRs. Computational simulations suggest that meromycolate chain dynamics regulate mycolate head group movement, thereby modulating GEM-TCR activity. Our findings have significant implications for the design of future vaccines that target GEM T cells., Competing Interests: The authors declare no conflict of interest.

We present blind predictions submitted to the SAMPL5 challenge on calculating distribution coefficients. The predictions were based on estimating the solvation free energies in water and cyclohexane of the 53 compounds in the challenge. These free energies were computed using alchemical free energy simulations based on a hybrid all-atom/coarse-grained model. The compounds were treated with the general Amber force field, whereas the solvent molecules were treated with the Elba coarse-grained model. Considering the simplicity of the solvent model and that we approximate the distribution coefficient with the partition coefficient of the neutral species, the predictions are of good accuracy. The correlation coefficient, R is 0.64, 82 % of the predictions have the correct sign and the mean absolute deviation is 1.8 log units. This is on a par with or better than the other simulation-based predictions in the challenge. We present an analysis of the deviations to experiments and compare the predictions to another submission that used all-atom solvent.

We present our predictions for the SAMPL4 hydration free energy challenge. Extensive all-atom Monte Carlo simulations were employed to sample the compounds in explicit solvent. While the focus of our study was to demonstrate well-converged and reproducible free energies, we attempted to address the deficiencies in the general Amber force field force field with a simple QM/MM correction. We show that by using multiple independent simulations, including different starting configurations, and enhanced sampling with parallel tempering, we can obtain well converged hydration free energies. Additional analysis using dihedral angle distributions, torsion-root mean square deviation plots and thermodynamic cycles support this assertion. We obtain a mean absolute deviation of 1.7 kcal mol(-1) and a Kendall's τ of 0.65 compared with experiment.

Recent efforts in the computational evaluation of the thermodynamic properties of water molecules have resulted in the development of promising new in silico methods to evaluate the role of water in ligand binding. These methods include WaterMap, SZMAP, GRID/CRY probe, and Grand Canonical Monte Carlo simulations. They allow the prediction of the position and relative free energy of the water molecule in the protein active site and the analysis of the perturbation of an explicit water network (WNP) as a consequence of ligand binding. We have for the first time extended these approaches toward the prediction of kinetics for small molecules and of relative free energy of binding with a focus on the perturbation of the water network and application to large diverse data sets. Our results support a qualitative correlation between the residence time of 12 related triazine adenosine A(2A) receptor antagonists and the number and position of high energy trapped solvent molecules. From a quantitative viewpoint, we successfully applied these computational techniques as an implicit solvent alternative, in linear combination with a molecular mechanics force field, to predict the relative ligand free energy of binding (WNP-MMSA). The applicability of this linear method, based on the thermodynamics additivity principle, did not extend to 375 diverse A(2A) receptor antagonists. However, a fast but effective method could be enabled by replacing the linear approach with a machine learning technique using probabilistic classification trees, which classified the binding affinity correctly for 90% of the ligands in the training set and 67% in the test set.

Many limitations of current computer-aided drug design arise from the difficulty of reliably predicting the binding affinity of a small molecule to a biological target. There is thus a strong interest in novel computational methodologies that claim predictions of greater accuracy than current scoring functions, and at a throughput compatible with the rapid pace of drug discovery in the pharmaceutical industry. Notably, computational methodologies firmly rooted in statistical thermodynamics have received particular attention in recent years. Yet free energy calculations can be daunting to learn for a novice user because of numerous technical issues and various approaches advocated by experts in the field. The purpose of this article is to provide an overview of the current capabilities of free energy calculations and to discuss the applicability of this technology to drug discovery.

Nonequilibrium (NE) free energy methods are embarrassingly parallel and may be very conveniently run on desktop computers using distributed computing software. In recent years there has been a proliferation of NE methods, but these approaches have barely, if at all, been used in the context of calculating protein-ligand binding free energies. In a recent study by these authors, different combinations of NE methods with various test systems were compared and protocols identified which yielded results as accurate as replica exchange thermodynamic integration (RETI). The NE approaches, however, lend themselves to extensive parallelization through the use of distributed computing. Here the best performing of those NE protocols, a replica exchange method using Bennett's acceptance ratio as the free energy estimator (RENE), is applied to two sets of congeneric inhibitors bound to neuraminidase and cyclooxygenase-2. These protein-ligand systems were originally studied with RETI, giving results to which NE and RENE simulations are compared. These NE calculations were carried out on a large, highly distributed group of low-performance desktop computers which are part of a Condor pool. RENE was found to produce results of a predictive quality at least as good as RETI in less than half the wall clock time. However, non-RE NE results were found to be far less predictive. In addition, the RENE method successfully identified a localized region of rapidly changing free energy gradients without the need for prior investigation. These results suggest that the RENE protocol is appropriate for use in the context of predicting protein-ligand binding free energies and that it can offer advantages over conventional, equilibrium approaches.

A pattern recognition algorithm for the alignment of drug-like molecules has been implemented. The method is based on the calculation of quantum mechanical derived local properties defined on a molecular surface. This approach has been shown to be very useful in attempting to derive generalized, non-atom based representations of molecular structure. The visualization of these surfaces is described together with details of the methodology developed for their use in molecular overlay and similarity calculations. In addition, this paper also introduces an additional local property, the local curvature (C (L)), which can be used together with the quantum mechanical properties to describe the local shape. The method is exemplified using some problems representing common tasks encountered in molecular similarity.

A methodology for the calculation of the free energy difference between a pair of molecules of arbitrary topology is proposed. The protocol relies on a dual-topology paradigm, a softening of the intermolecular interactions, and a constraint that prevents the perturbed molecules from drifting away from each other at the end states. The equivalence and the performance of the methodology against a single-topology approach are demonstrated on a pair of harmonic oscillators, the calculation of the relative solvation free energy of ethane and methanol, and the relative binding free energy of two congeneric inhibitors of cyclooxygenase 2. The stability of two alternative binding modes of an inhibitor of cyclin-dependent kinase 2 is then investigated. Finally, the relative binding free energy of two structurally different inhibitors of cyclin-dependent kinase 2 is calculated. The proposed methodology allows the study of a range of problems that are beyond the reach of traditional relative free energy calculation protocols and should prove useful in drug design studies.

The Schools Malaria Project (http://emalaria.soton.ac.uk/) brings together school students with university researchers in the hunt for a new antimalaria drug. The design challenge being offered to students is to use a distributed drug search and selection system to design potential antimalaria drugs. The system is accessed via a Web interface. This e-science project displays the results of the trials in an accessible manner, giving students an opportunity for discussion and debate both with peers and with the university contacts. The project has been implemented by using distributed computing techniques, spreading computer load over a network of machines that cross institutional boundaries, forming a grid. This provides access to greater computing power and allows a much more complex and detailed formulation of the drug design problem to be tackled for research, teaching, and learning.

Biomolecular computer simulations are now widely used not only in an academic setting to understand the fundamental role of molecular dynamics on biological function, but also in the industrial context to assist in drug design. In this paper, two applications of Grid computing to this area will be outlined. The first, involving the coupling of distributed computing resources to dedicated Beowulf clusters, is targeted at simulating protein conformational change using the Replica Exchange methodology. In the second, the rationale and design of a database of biomolecular simulation trajectories is described. Both applications illustrate the increasingly important role modern computational methods are playing in the life sciences.

Conventional paper-based scientific publications are limited in the amount of data and interaction they can provide. Simply placing an electronic copy of such a publication on the web makes for easier distribution, but more can be achieved by making use of the technology to navigate through the paper, to show the links between calculated parameters and the data, and provide access to the chain of calculated results all the way back to the raw data and ultimately the laboratory notebook. While the paper version of this publication is in a relatively conventional form, an interactive demonstrator version (available at http://epaper.combe.chem.soton.ac.uk and as an installable package) illustrates the concepts of publication@source, whereby all the figures and data presented in the paper are linked back to the original raw data together with a description of the processes by which the raw data were analysed. This level of interactivity is achieved using semantic technologies, which have the additional advantage of making the final document subsequently available and navigable by automated techniques. We present the combined information from experimental studies of surface tension and second harmonic generation (SHG) on the behaviour of benzo-15-crown-5 at the solution/air interface, together with a molecular dynamics computer simulation, to demonstrate how the simulation aids the interpretation of the SHG experiment. The adsorption isotherm was determined using SHG and fitted to a Langmuir isotherm giving Delta ads G0=26 kJ mol-1.

Reversible Digitally Filtered Molecular Dynamics (RDFMD) is a method of amplifying or suppressing motions in a molecular dynamics simulation, through the application of a digital filter to the simulation velocities. RDFMD and its derivatives have been previously used to promote conformational motions in liquid-phase butane, the Syrian hamster prion protein, alanine dipeptide, and the pentapeptide, YPGDV. The RDFMD method has associated with it a number of parameters that require specification to optimize the desired response. In this paper methods for the systematic analysis of these parameters are presented and applied to YPGDV with the specific emphasis of ensuring a gentle and progressive method that produces maximum conformation change from the energy put into the system. The portability of the new parameter set is then shown with an application to the M20 loop of E-coli dihydrofolate reductase. A conformational change is induced from a closed to an open structure similar to that seen in the DHFR-NADP(+) complex.