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.

Hydration free energy (HFE) calculations are often used to assess the performance of biomolecular force fields and the quality of assigned parameters. The AMOEBA polarizable force field moves beyond traditional pairwise additive models of electrostatics and may be expected to improve upon predictions of thermodynamic quantities such as HFEs over and above fixed-point-charge models. The recent SAMPL4 challenge evaluated the AMOEBA polarizable force field in this regard but showed substantially worse results than those using the fixed-point-charge GAFF model. Starting with a set of automatically generated AMOEBA parameters for the SAMPL4 data set, we evaluate the cumulative effects of a series of incremental improvements in parametrization protocol, including both solute and solvent model changes. Ultimately, the optimized AMOEBA parameters give a set of results that are not statistically significantly different from those of GAFF in terms of signed and unsigned error metrics. This allows us to propose a number of guidelines for new molecule parameter derivation with AMOEBA, which we expect to have benefits for a range of biomolecular simulation applications such as protein-ligand binding studies.

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.

Water molecules play integral roles in the formation of many protein-ligand complexes, and recent computational efforts have been focused on predicting the thermodynamic properties of individual waters and how they may be exploited in rational drug design. However, when water molecules form highly coupled hydrogen-bonding networks, there is, as yet, no method that can rigorously calculate the free energy to bind the entire network or assess the degree of cooperativity between waters. In this work, we report theoretical and methodological developments to the grand canonical Monte Carlo simulation technique. Central to our results is a rigorous equation that can be used to calculate efficiently the binding free energies of water networks of arbitrary size and complexity. Using a single set of simulations, our methods can locate waters, estimate their binding affinities, capture the cooperativity of the water network, and evaluate the hydration free energy of entire protein binding sites. Our techniques have been applied to multiple test systems and compare favorably to thermodynamic integration simulations and experimental data. The implications of these methods in drug design are discussed.

We present an efficient all-atom/coarse-grained hybrid model and apply it to membrane processes. This model is an extension of the all-atom/ELBA model applied previously to processes in water. Here, we improve the efficiency of the model by implementing a multiple-time step integrator that allows the atoms and the coarse-grained beads to be propagated at different timesteps. Furthermore, we fine-tune the interaction between the atoms and the coarse-grained beads by computing the potential of mean force of amino acid side chain analogs along the membrane normal and comparing to atomistic simulations. The model was independently validated on the calculation of small-molecule partition coefficients. Finally, we apply the model to membrane peptides. We studied the tilt angle of the Walp23 and Kalp23 helices in two different model membranes and the stability of the glycophorin A dimer. The model is efficient, accurate, and straightforward to use, as it does not require any extra interaction particles, layers of atomistic solvent molecules or tabulated potentials, thus offering a novel, simple approach to study membrane processes.

Hydrogenases use complex metal 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 Fe(2+) 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.

The use of the Zwanzig equation in the calculation of single-step perturbations to provide first-principles (ab initio) quantum mechanics (QM) correction terms to molecular mechanics (MM) free energy cycles is well established. A rigorous test of the ability to converge such calculations would be very useful in this context. In this work, we perform a direct assessment of the convergence of the MM to QM perturbation, by attempting the reverse QM to MM perturbation. This required the generation of extensive QM molecular dynamics trajectories, using density functional theory (DFT), within the representative biological system of a DNA adenosine-thymidine dimer. Over 100 ps of dynamics with the PBE functional and 6.25 ps with the LDA functional were generated. We demonstrate that calculations with total potential energies are very poorly convergent due to a lack of overlap of phase space distributions between ensembles. While not theoretically rigorous, the use of interaction energies provides far superior convergence, despite the presence of nonclassical charge transfer effects within the DFT trajectories. The source of poor phase space overlap for total energies is diagnosed, the approximate quantification of overlaps suggesting that even for the comparatively simple system considered here convergence of total energy calculations within a reasonable simulation time is unfeasible.

A comprehensive conformational analysis of both 2,3-difluorobutane diastereomers is presented based on density functional theory calculations in vacuum and in solution, as well as NMR experiments in solution. While for 1,2-difluoroethane the fluorine gauche effect is clearly the dominant effect determining its conformation, it was found that for 2,3-difluorobutane there is a complex interplay of several effects, which are of similar magnitude but often of opposite sign. As a result, unexpected deviations in dihedral angles, relative conformational energies and populations are observed which cannot be rationalised only by chemical intuition. Furthermore, it was found that it is important to consider the free energies of the various conformers, as these lead to qualitatively different results both in vacuum and in solvent, when compared to calculations based only on the electronic energies. In contrast to expectations, it was found that vicinal syn-difluoride introduction in the butane and by extension, longer hydrocarbon chains, is not expected to lead to an effective stabilisation of the linear conformation. Our findings have implications for the use of the vicinal difluoride motif for conformational control., (© 2015 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.)

Lipoyl cofactors are essential for living organisms and are produced by the insertion of two sulfur atoms into the relatively unreactive C-H bonds of an octanoyl substrate. This reaction requires lipoyl synthase, a member of the radical S-adenosylmethionine (SAM) enzyme superfamily. In the present study, we solved crystal structures of lipoyl synthase with two [4Fe-4S] clusters bound at opposite ends of the TIM barrel, the usual fold of the radical SAM superfamily. The cluster required for reductive SAM cleavage conserves the features of the radical SAM superfamily, but the auxiliary cluster is bound by a CX4CX5C motif unique to lipoyl synthase. The fourth ligand to the auxiliary cluster is an extremely unusual serine residue. Site-directed mutants show this conserved serine ligand is essential for the sulfur insertion steps. One crystallized lipoyl synthase (LipA) complex contains 5'-methylthioadenosine (MTA), a breakdown product of SAM, bound in the likely SAM-binding site. Modelling has identified an 18 Å (1 Å=0.1 nm) deep channel, well-proportioned to accommodate an octanoyl substrate. These results suggest that the auxiliary cluster is the likely sulfur donor, but access to a sulfide ion for the second sulfur insertion reaction requires the loss of an iron atom from the auxiliary cluster, which the serine ligand may enable.

Water molecules are commonplace in protein binding pockets, where they can typically form a complex between the protein and a ligand or become displaced upon ligand binding. As a result, it is often of great interest to establish both the binding free energy and location of such molecules. Several approaches to predicting the location and affinity of water molecules to proteins have been proposed and utilized in the literature, although it is often unclear which method should be used under what circumstances. We report here a comparison between three such methodologies, Just Add Water Molecules (JAWS), Grand Canonical Monte Carlo (GCMC), and double-decoupling, in the hope of understanding the advantages and limitations of each method when applied to enclosed binding sites. As a result, we have adapted the JAWS scoring procedure, allowing the binding free energies of strongly bound water molecules to be calculated to a high degree of accuracy, requiring significantly less computational effort than more rigorous approaches. The combination of JAWS and GCMC offers a route to a rapid scheme capable of both locating and scoring water molecules for rational drug design.

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.

Schemes of increasing sophistication for obtaining free energies of binding have been developed over the years, where configurational sampling is used to include the all-important entropic contributions to the free energies. However, the quality of the results will also depend on the accuracy with which the intermolecular interactions are computed at each molecular configuration. In this context, the energy change associated with the rearrangement of electrons (electronic polarization and charge transfer) upon binding is a very important effect. Classical molecular mechanics force fields do not take this effect into account explicitly, and polarizable force fields and semiempirical quantum or hybrid quantum-classical (QM/MM) calculations are increasingly employed (at higher computational cost) to compute intermolecular interactions in free-energy schemes. In this work, we investigate the use of large-scale quantum mechanical calculations from first-principles as a way of fully taking into account electronic effects in free-energy calculations. We employ a one-step free-energy perturbation (FEP) scheme from a molecular mechanical (MM) potential to a quantum mechanical (QM) potential as a correction to thermodynamic integration calculations within the MM potential. We use this approach to calculate relative free energies of hydration of small aromatic molecules. Our quantum calculations are performed on multiple configurations from classical molecular dynamics simulations. The quantum energy of each configuration is obtained from density functional theory calculations with a near-complete psinc basis set on over 600 atoms using the ONETEP program.

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.

A recently developed coarse-grain model is applied to simulate hydrated membranes containing the lamellar lipid DOPC and the nonlamellar lipid DOPE. In a first series of simulations, DOPC-water and DOPE-water systems are shown to form respectively bilayers and inverse hexagonal phases, in agreement with the well-known behaviour observed experimentally. A second set of calculations is then run to investigate several fundamental physical features of mixed DOPC-DOPE bilayers at different relative compositions. In particular, a quantitative characterisation is obtained of the internal distributions (profiles) of lateral pressure and electrical potential. These two properties, very difficult to measure experimentally, are thought to underpin many key membrane phenomena, including nonspecific lipid-mediated mechanisms of protein regulation. The molecular origin of the distributions, and their dependence on changes in the DOPC: DOPE ratio, are explained through an analysis of separate contributions from individual interaction types and molecular groups.

The identification of novel binding-site conformations can greatly assist the progress of structure-based ligand design projects. Diverse pocket shapes drive medicinal chemistry to explore a broader chemical space and thus present additional opportunities to overcome key drug discovery issues such as potency, selectivity, toxicity, and pharmacokinetics. We report a new automated approach to diverse pocket selection, PocketAnalyzer(PCA), which applies principal component analysis and clustering to the output of a grid-based pocket detection algorithm. Since the approach works directly with pocket shape descriptors, it is free from some of the problems hampering methods that are based on proxy shape descriptors, e.g. a set of atomic positional coordinates. The approach is technically straightforward and allows simultaneous analysis of mutants, isoforms, and protein structures derived from multiple sources with different residue numbering schemes. The PocketAnalyzer(PCA) approach is illustrated by the compilation of diverse sets of pocket shapes for aldose reductase and viral neuraminidase. In both cases this allows identification of novel computationally derived binding-site conformations that are yet to be observed crystallographically. Indeed, known inhibitors capable of exploiting these novel binding-site conformations are subsequently identified, thereby demonstrating the utility of PocketAnalyzer(PCA) for rationalizing and improving the understanding of the molecular basis of protein-ligand interaction and bioactivity. A Python program implementing the PocketAnalyzer(PCA) approach is available for download under an open-source license ( http://sourceforge.net/projects/papca/ or http://cpclab.uni-duesseldorf.de/downloads ).

Background: Cyclin-dependent kinases 2, 4 and 6 (Cdk2, Cdk4, Cdk6) are closely structurally homologous proteins which are classically understood to control the transition from the G1 to the S-phases of the cell cycle by combining with their appropriate cyclin D or cyclin E partners to form kinase-active holoenzymes. Deregulation of Cdk4 is widespread in human cancer, CDK4 gene knockout is highly protective against chemical and oncogene-mediated epithelial carcinogenesis, despite the continued presence of CDK2 and CDK6; and overexpresssion of Cdk4 promotes skin carcinogenesis. Surprisingly, however, Cdk4 kinase inhibitors have not yet fulfilled their expectation as 'blockbuster' anticancer agents. Resistance to inhibition of Cdk4 kinase in some cases could potentially be due to a non-kinase activity, as recently reported with epidermal growth factor receptor., Results: A search for a potential functional site of non-kinase activity present in Cdk4 but not Cdk2 or Cdk6 revealed a previously-unidentified loop on the outside of the C'-terminal non-kinase domain of Cdk4, containing a central amino-acid sequence, Pro-Arg-Gly-Pro-Arg-Pro (PRGPRP). An isolated hexapeptide with this sequence and its cyclic amphiphilic congeners are selectively lethal at high doses to a wide range of human cancer cell lines whilst sparing normal diploid keratinocytes and fibroblasts. Treated cancer cells do not exhibit the wide variability of dose response typically seen with other anticancer agents. Cancer cell killing by PRGPRP, in a cyclic amphiphilic cassette, requires cells to be in cycle but does not perturb cell cycle distribution and is accompanied by altered relative Cdk4/Cdk1 expression and selective decrease in ATP levels. Morphological features of apoptosis are absent and cancer cell death does not appear to involve autophagy., Conclusion: These findings suggest a potential new paradigm for the development of broad-spectrum cancer specific therapeutics with a companion diagnostic biomarker and a putative functional site for kinase-unrelated activities of Cdk4.

Triclocarban and triclosan, two potent antibacterial molecules present in many consumer products, have been subject to growing debate on a number of issues, particularly in relation to their possible role in causing microbial resistance. In this computational study, we present molecular-level insights into the interaction between these antimicrobial agents and hydrated phospholipid bilayers (taken as a simple model for the cell membrane). Simulations are conducted by a novel 'dual-resolution' molecular dynamics approach which combines accuracy with efficiency: the antimicrobials, modelled atomistically, are mixed with simplified (coarse-grain) models of lipids and water. A first set of calculations is run to study the antimicrobials' transfer free energies and orientations as a function of depth inside the membrane. Both molecules are predicted to preferentially accumulate in the lipid headgroup-glycerol region; this finding, which reproduces corresponding experimental data, is also discussed in terms of a general relation between solute partitioning and the intramembrane distribution of pressure. A second set of runs involves membranes incorporated with different molar concentrations of antimicrobial molecules (up to one antimicrobial per two lipids). We study the effects induced on fundamental membrane properties, such as the electron density, lateral pressure and electrical potential profiles. In particular, the analysis of the spontaneous curvature indicates that increasing antimicrobial concentrations promote a 'destabilizing' tendency towards non-bilayer phases, as observed experimentally. The antimicrobials' influence on the self-assembly process is also investigated. The significance of our results in the context of current theories of antimicrobial action is discussed., (© 2010 The Royal Society)

We present a molecular simulation protocol to compute free energies of binding, which combines a QM/MM correction term with rigorous classical free energy techniques, thereby accounting for electronic polarization effects. Relative free energies of binding are first computed using classical force fields, Monte Carlo sampling, and replica exchange thermodynamic integration. Snapshots of the configurations at the end points of the perturbation are then subjected to DFT-QM/MM single-point calculations using the B3LYP functional and a range of basis sets. The resulting quantum mechanical energies are then processed using the Zwanzig equation to give free energies incorporating electronic polarization. Our approach is conceptually simple and does not require tightly coupled QM and MM software. The method has been validated by calculating the relative free energies of hydration of methane and water and the relative free energy of binding of two inhibitors of cyclooxygenase-2. Closed thermodynamic cycles are obtained across different pathways, demonstrating the correctness of the technique, although significantly more sampling is required for the protein-ligand system. Our method offers a simple and effective way to incorporate quantum mechanical effects into computed free energies of binding.

A new coarse-grain model for molecular dynamics simulation of lipid membranes is presented. Following a simple and conventional approach, lipid molecules are modeled by spherical sites, each representing a group of several atoms. In contrast to common coarse-grain methods, two original (interdependent) features are here adopted. First, the main electrostatics are modeled explicitly by charges and dipoles, which interact realistically through a relative dielectric constant of unity (ε(r) = 1). Second, water molecules are represented individually through a new parametrization of the simple Stockmayer potential for polar fluids; each water molecule is therefore described by a single spherical site embedded with a point dipole. The force field is shown to accurately reproduce the main physical properties of single-species phospholipid bilayers comprising dioleoylphosphatidylcholine (DOPC) and dioleoylphosphatidylethanolamine (DOPE) in the liquid crystal phase, as well as distearoylphosphatidylcholine (DSPC) in the liquid crystal and gel phases. Insights are presented into fundamental properties and phenomena that can be difficult or impossible to study with alternative computational or experimental methods. For example, we investigate the internal pressure distribution, dipole potential, lipid diffusion, and spontaneous self-assembly. Simulations lasting up to 1.5 microseconds were conducted for systems of different sizes (128, 512 and 1058 lipids); this also allowed us to identify size-dependent artifacts that are expected to affect membrane simulations in general. Future extensions and applications are discussed, particularly in relation to the methodology's inherent multiscale capabilities.

Microtubules are supramolecular structures that make up the cytoskeleton and strongly affect the mechanical properties of the cell. Within the cytoskeleton filaments, the microtubule (MT) exhibits by far the highest bending stiffness. Bending stiffness depends on the mechanical properties and intermolecular interactions of the tubulin dimers (the MT building blocks). Computational molecular modeling has the potential for obtaining quantitative insights into this area. However, to our knowledge, standard molecular modeling techniques, such as molecular dynamics (MD) and normal mode analysis (NMA), are not yet able to simulate large molecular structures like the MTs; in fact, their possibilities are normally limited to much smaller protein complexes. In this work, we developed a multiscale approach by merging the modeling contribution from MD and NMA. In particular, MD simulations were used to refine the molecular conformation and arrangement of the tubulin dimers inside the MT lattice. Subsequently, NMA was used to investigate the vibrational properties of MTs modeled as an elastic network. The coarse-grain model here developed can describe systems of hundreds of interacting tubulin monomers (corresponding to up to 1,000,000 atoms). In particular, we were able to simulate coarse-grain models of entire MTs, with lengths up to 350 nm. A quantitative mechanical investigation was performed; from the bending and stretching modes, we estimated MT macroscopic properties such as bending stiffness, Young modulus, and persistence length, thus allowing a direct comparison with experimental data., (Copyright © 2010 Biophysical Society. Published by Elsevier Inc. All rights reserved.)

Structure-based drug design could benefit greatly from computational methodologies that accurately predict the binding affinity of small compounds to target biomolecules. However, the current scoring functions used to rank compounds in virtual screens by docking are not sufficiently accurate to guide reliably the design of tight binding ligands. Thus, there is strong interest in methodologies based on molecular simulations of protein-ligand complexes which are perceived to be more accurate and, with advances in computing power, amenable to routine use. This report provides an overview of the technical details necessary to understand, execute and analyze binding free energy calculations, using free energy perturbation or thermodynamic integration methods. Examples of possible applications in structure-based drug design are discussed. Current methodological limitations are highlighted as well as a number of ongoing developments to improve the scope, reliability, and practicalities of free energy calculations. These efforts are paving the way for a more common use of free energy calculations in molecular design., (Copyright © 2010 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.)

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.

Docking into multiple receptor conformations ("ensemble docking") has been proposed, and employed, in the hope that it may account for receptor flexibility in virtual screening and thus provide higher enrichments than docking into single rigid receptor structures. The statistical analyses presented in this paper provide quantitative evidence that in some cases docking into a crystallographically derived conformational ensemble does indeed yield better enrichment than docking into any of the individual members of the ensemble. However, these "successful" ensembles account for only a minority of those examined and it would not have been possible to prospectively predict their identity using only protein structural information. A more frequently observed outcome is that the ensemble enrichment is higher than the mean of the enrichments provided by its individual members. An additional and promising finding is that, if a set of known active compounds is available, an approach based on induced-fit docking appears to be a reliable way to construct ensembles which provide relatively high enrichments.

Our recently developed coarse-grain model for dimyristoylphosphatidylcholine (DMPC) has been improved and extended to dioleylphosphatidylcholine (DOPC), a more typical constituent of real biological membranes. Single-component DMPC and DOPC bilayers have been simulated using microsecond-long molecular dynamics. We investigated properties that are difficult or impossible to access experimentally, such as the pressure distribution, the spontaneous curvature and the diffusion pattern of individual lipid molecules. Moreover, we studied the dipole potential, a basic physical feature of paramount biological importance that cannot be currently modelled by other coarse-grain approaches. In fact, a complete representation of the system electrostatics and a realistic description of the water component make our method unique amongst the existing coarse-grain membrane models. The spontaneous permeation of water, a phenomenon out of reach of standard atomistic models, was also observed and quantified; this was possible thanks to the efficiency of our model, which is about two orders of magnitude less computationally expensive than atomic-level counterparts. Results are generally in good agreement with the literature data. Further model extensions and future applications are proposed.

The transmembrane permeation of eight small (molecular weight <100) organic molecules across a phospholipid bilayer is investigated by multiscale molecular dynamics simulation. The bilayer and hydrating water are represented by simplified, efficient coarse-grain models, whereas the permeating molecules are described by a standard atomic-level force-field. Permeability properties are obtained through a refined version of the z-constraint algorithm. By constraining each permeant at selected depths inside the bilayer, we have sampled free energy differences and diffusion coefficients across the membrane. These data have been combined, according to the inhomogeneous solubility-diffusion model, to yield the permeability coefficients. The results are generally consistent with previous atomic-level calculations and available experimental data. Computationally, our multiscale approach proves 2 orders of magnitude faster than traditional atomic-level methods.

One of the factors preventing the general application of free energy methods in rational drug design remains the lack of sufficient computational resources. Many nonequilibrium (NE) free energy methods, however, are easily made embarrassingly parallel in comparison to equilibrium methods and may be conveniently run on desktop computers using distributed computing software. In recent years, there has been a proliferation of NE methods, but the general applicability of these approaches has not been determined. In this study, a subset including only those NE methods which are easily parallelised were considered for examination, with a view to their application to the prediction of protein-ligand binding affinities. A number of test systems were examined, including harmonic oscillator (HO) systems and the calculation of relative free energies of hydration of water-methane. The latter system uses identical potentials to the protein ligand case and is therefore an appropriate model system on which methods may be tested. As well as investigating existing protocols, a replica exchange NE approach was developed, which was found to offer advantages over conventional methods. It was found that Rosenbluth-based approaches to optimizing the NE work values used in NE free energy estimates were not consistent in the improvements in accuracy achieved and that, given their computational cost, the simple approach of taking each work value in an unbiased way is to be preferred. Of the two free energy estimators examined, Bennett's acceptance ratio was the most consistent and is, therefore, to be preferred over the Jarzynski estimator. The recommended protocols may be run very efficiently within a distributed computing environment and are of similar accuracy and precision to equilibrium free energy methods.

HIV-1 protease performs a vital step in the propagation of the HIV virus and is therefore an important drug target in the treatment of AIDS. It consists of a homodimer, with access to the active site limited by two protein flaps. NMR studies have identified two time scales of motions that occur in these flaps, and it is thought that the slower of these is responsible for a conformational change that makes the protein ligand-accessible. This motion occurs on a time scale outside that achievable using traditional molecular dynamics simulations. Reversible Digitally Filtered Molecular Dynamics (RDFMD) is a method that amplifies low frequency motions associated with conformational change and has recently been applied to, among others, E. coli dihydrofolate reductase, inducing a conformational change between known crystal structures. In this paper, the conformational motions of HIV-1 protease produced during MD and RDFMD simulations are presented, including movement between the known semiopen and closed conformations, and the opening and closing of the protein flaps.

The HIV-1 IN enzyme is one of three crucial virally encoded enzymes (HIV-1 IN, HIV-1 PR, and HIV-1 RT) involved in the life-cycle of the HIV-1 virus, making it an attractive target in the development of drugs against the AIDS virus. The structure and mechanism of the HIV-1 IN enzyme is the least understood of the three enzymes due to the lack of three-dimensional structural information. X-ray cystallographic studies have not yet been able to resolve the full-length structure, and studies have been mainly focused on the catalytic domain. This central domain possesses an important catalytic loop observed to overhang the active site, and experimental studies have shown that its dynamics affects the catalytic activity of mutant HIV-1 IN enzymes. In this study, the enhanced sampling technique, Reversible Digitally Filtered Molecular Dynamics (RDFMD), has been applied to the catalytic domain of the WT and G140A/G149A HIV-1 IN enzymes and has highlighted significant differences between the behavior of the catalytic loop which may explain the decrease of activity observed in experimental studies for this mutant.

We present the theory of thermal equivalence in the framework of the Peyrard-Bishop model and some of its anharmonic variants. The thermal equivalence gives rise to a melting index τ which maps closely the experimental DNA melting temperatures for short DNA sequences. We show that the efficient calculation of the melting index can be used to analyse the parameters of the Peyrard-Bishop model and propose an improved set of Morse potential parameters. With this new set we are able to calculate some of the experimental melting temperatures to ± 1.2 °C. We review some of the concepts of sequencing probe design and show how to use the melting index to explore the possibilities of gene coverage by tuning the model parameters.

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.

The identification of lead molecules using computational modeling often relies on approximate, high-throughput approaches, of limited accuracy. We show here that, with a methodology we recently developed, it is possible to predict the relative binding free energies of structurally diverse ligands of the estrogen receptor-alpha using a rigorous statistical thermodynamics approach. Predictions obtained from the simulations with an explicit solvation model are in good qualitative agreement with experimental data, while simulations with implicit solvent models or rank ordering by empirical scoring functions yield predictions of lower quality. In addition, it is shown that free energy techniques can be used to select the most likely binding mode from a set of possible orientations generated by a docking program. It is suggested that the free energy techniques outlined in this study can be used to rank-order, by potency, structurally diverse compounds identified by virtual screening, de novo design or scaffold hopping programs.

Background: Sequencing by hybridisation is an effective method for obtaining large amounts of DNA sequence information at low cost. The efficiency of SBH depends on the design of the probe library to provide the maximum information for minimum cost. Long probes provide a higher probability of non-repeated sequences but lead to an increase in the number of probes required whereas short probes may not provide unique sequence information due to repeated sequences. We have investigated the effect of probe length, use of reference sequences, and thermal filtering on the design of probe libraries for several highly variable target DNA sequences., Results: We designed overlapping probe libraries for a range of highly variable drug target genes based on known sequence information and develop a formal terminology to describe probe library design. We find that for some targets these libraries can provide good coverage of a previously unseen target whereas for others the coverage is less than 30%. The optimal probe length varies from as short at 12 nt to as large as 19 nt and depends on the sequence, its variability, and the stringency of thermal filtering. It cannot be determined from inspection of an example gene sequence., Conclusions: Optimal probe length and the optimal number of reference sequences used to design a probe library are highly target specific for highly variable sequencing targets. The optimum design cannot be determined simply by inspection of input sequences or of alignments but only by detailed analysis of the each specific target. For highly variable sequences, shorter probes can in some cases provide better information than longer probes. Probe library design would benefit from a general purpose tool for analysing these issues. The formal terminology developed here and the analysis approaches it is used to describe will contribute to the development of such tools.

A simplified particle-based computer model for hydrated phospholipid bilayers has been developed and applied to quantitatively predict the major physical features of fluid-phase biomembranes. Compared with available coarse-grain methods, three novel aspects are introduced. First, the main electrostatic features of the system are incorporated explicitly via charges and dipoles. Second, water is accurately (yet efficiently) described, on an individual level, by the soft sticky dipole model. Third, hydrocarbon tails are modeled using the anisotropic Gay-Berne potential. Simulations are conducted by rigid-body molecular dynamics. Our technique proves 2 orders of magnitude less demanding of computational resources than traditional atomic-level methodology. Self-assembled bilayers quantitatively reproduce experimental observables such as electron density, compressibility moduli, dipole potential, lipid diffusion, and water permeability. The lateral pressure profile has been calculated, along with the elastic curvature constants of the Helfrich expression for the membrane bending energy; results are consistent with experimental estimates and atomic-level simulation data. Several of the results presented have been obtained for the first time using a coarse-grain method. Our model is also directly compatible with atomic-level force fields, allowing mixed systems to be simulated in a multiscale fashion.

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.

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.

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 approximately 0.1 micros 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.

Continuum electrostatics is combined with rigorous free-energy calculations in an effort to deliver a reliable and efficient method for in silico lead optimization. The methodology is tested by calculation of the relative binding free energies of a set of inhibitors of neuraminidase, cyclooxygenase2, and cyclin-dependent kinase 2. The calculated free energies are compared to the results obtained with explicit solvent simulations and empirical scoring functions. For cyclooxygenase2, deficiencies in the continuum electrostatics theory are identified and corrected with a modified simulation protocol. For neuraminidase, it is shown that a continuum representation of the solvent leads to markedly different protein-ligand interactions compared to the explicit solvent simulations, and a reconciliation of the two protocols is problematic. Cyclin-dependent kinase 2 proves more challenging, and none of the methods employed in this study yield high quality predictions. Despite the differences observed, for these systems, the use of an implicit solvent framework to predict the ranking of congeneric inhibitors to a protein is shown to be faster, as accurate or more accurate than the explicit solvent protocol, and superior to empirical scoring schemes.

Contemporary structural biology has an increased emphasis on high-throughput methods. Biomolecular simulations can add value to structural biology via the provision of dynamic information. However, at present there are no agreed measures for the quality of biomolecular simulation data. In this Letter, we suggest suitable measures for the quality assurance of molecular dynamics simulations of biomolecules. These measures are designed to be simple, fast, and general. Reporting of these measures in simulation papers should become an expected practice, analogous to the reporting of comparable quality measures in protein crystallography. We wish to solicit views and suggestions from the simulation community on methods to obtain reliability measures from molecular-dynamics trajectories. In a database which provides access to previously obtained simulations [Formula: see text] for example BioSimGrid ( http://www.biosimgrid.org/ ) [Formula: see text] the user needs to be confident that the simulation trajectory is suitable for further investigation. This can be provided by the simulation quality measures which a user would examine prior to more extensive analyses.

The Generalized Born Surface Area theory (GBSA) has become a popular method to model the solvation of biomolecules. While efficient in the context of molecular dynamics simulations, GBSA calculations do not integrate well with Monte Carlo simulations because of the nonlocal nature of the Generalized Born energy. We present a method by which Monte Carlo Generalized Born simulations can be made seven to eight times faster on a protein-ligand binding free energy calculation with little or no loss of accuracy. The method can be employed in any type of Monte Carlo or Hybrid Monte Carlo-molecular dynamics simulation and should prove useful in numerous applications.

Present chemical data storage methodologies place many restrictions on the use of the stored data. The absence of sufficient high-quality metadata prevents intelligent computer access to the data without human intervention. This creates barriers to the automation of data mining in activities such as quantitative structure-activity relationship modelling. The application of Semantic Web technologies to chemical data is shown to reduce these limitations. The use of unique identifiers and relationships (represented as uniform resource identifiers, URIs, and resource description framework, RDF) held in a triplestore provides for greater detail and flexibility in the sharing and storage of molecular structures and properties.

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.

To reach their biological target, drugs have to cross cell membranes, and understanding passive membrane permeation is therefore crucial for rational drug design. Molecular dynamics simulations offer a powerful way of studying permeation at the single molecule level. Starting from a computer model proven to be able to reproduce the physical properties of a biological membrane, the behaviour of small solutes and large drugs in a lipid bilayer has been studied. Analysis of dihedral angles shows that a few nano seconds are sufficient for the simulations to converge towards common values for those angles, even if the starting structures belong to different conformations. Results clearly show that, despite their difference in size, small solutes and large drugs tend to lie parallel to the bilayer normal and that, when moving from water solution into biomembranes, permeants lose degrees of freedom. This explains the experimental observation that partitioning and permeation are highly affected by entropic effects and are size-dependent. Tilted orientations, however, occur when they make possible the formation of hydrogen bonds. This helps to understand the reason why hydrogen bonding possibilities are an important parameter in cruder approaches which predict drug absorption after administration. Interestingly, hydration is found to occur even in the membrane core, which is usually considered an almost hydrophobic region. Simulations suggest the possibility for highly polar compounds like acetic acid to cross biological membranes while hydrated. These simulations prove useful for drug design in rationalising experimental observations and predicting solute behaviour in biomembranes.

Several methods for ultra high-throughput DNA sequencing are currently under investigation. Many of these methods yield very short blocks of sequence information (reads). Here we report on an analysis showing the level of genome sequencing possible as a function of read length. It is shown that re-sequencing and de novo sequencing of the majority of a bacterial genome is possible with read lengths of 20-30 nt, and that reads of 50 nt can provide reconstructed contigs (a contiguous fragment of sequence data) of 1000 nt and greater that cover 80% of human chromosome 1.

We discuss models fit to data collected by Duffy and Jorgensen to predict solvation free energies and partition equilibria of drugs, organic molecules, aromatic heterocycles, and other molecules. These data were originally examined using linear regression, but here more recently developed statistical models are applied. The data set is complicated due to the presence of discrepant observations and also curvature in the response. In some cases it is possible to discard a small number of the observations to get good fit to the data, but, in others, discarding an increasing proportion of the observations does not improve the fit. Our general preference is to use robust parameter estimation which downweights to reduce the influence of discrepant observations on the fitted models. Models are selected for four responses using linear or more complicated representations of the explanatory variables, such as cubic polynomials, B-splines, or smoothers via generalized additive models (GAMs). Variables are chosen using the traditional approach of formal tests to assess contribution to the fit of a model, and resampling methods including bootstrap are also considered to assess the prediction error for given models. Results of our analysis indicate that GAMs are an improvement on linear models for describing the data and making predictions. In general robust regression models and GAMs have the smallest conditional expected loss of prediction over the four responses. In addition, robust regression models offer the advantage of identifying molecules that perform poorly in the fit. In general, models were identified that yielded an improvement of approximately 50% in the conditional expected loss of prediction compared with the original parametrization of Duffy and Jorgensen. It was also found that the use of cross-validation to compare models was unreliable, and bootstrapping is preferred.

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.

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.