Harmonic modes as variables to approximately account for receptor flexibility in ligand–receptor docking simulations: Application to DNA minor groove ligand complex

An approach to approximately account for receptor flexibility in ligand–receptor docking simulations is described and applied to a DNA/Hoechst 33258 analogue complex. Harmonic modes corresponding to eigenvectors with small eigenvalues of the Hessian matrix of the potential energy function were used as independent variables to describe receptor flexibility. For the DNA minor groove ligand case most of the conformational difference between an energy minimized free DNA and ligand-bound structure could be assigned to 5–40 harmonic receptor modes with small eigenvalues. During docking, deformations of the DNA receptor structure in the subset of harmonic modes were limited using a simple penalty function that avoided the summation over all intrareceptor atom pairs. Significant improvement of the sterical fit between ligand and receptor was found upon relaxation of the DNA in the subset of harmonic modes after docking of the ligand at the position found in the known crystal structure. In addition, the harmonic mode relaxation resulted in DNA structures that were more similar to the energy minimized ligand-bound form. Although harmonic mode relaxation also leads to improved sterical fit for other ligand placements, the placement as observed in the crystal structure could still be identified as the site with the most favorable sterical interactions. Because relaxation in the harmonic modes is orders of magnitude faster than conventional energy minimization using all atom coordinates as independent variables, the approach might be useful as a preselection tool to recognize ligand binding sites accessible only upon small conformational changes of the receptor. The harmonic mode relaxed structures can only be considered as approximate structures because deformation of the receptor in the harmonic modes can lead to small perturbations of the stereochemical geometry of the molecule. Energy minimization of preselected ligand–DNA docking candidates in all atom coordinates is required to reduce these deviations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 287–300, 1999     

Free energy landscape of A-DNA to B-DNA conversion in aqueous solution

The interconversion between the well-characterized A- and B-forms of DNA is a structural transition for which the intermediate states and the free energy difference between the two endpoints are not known precisely. In the present study, the difference between the Root Mean Square Distance (RMSD) from canonical A-form and B-form DNA is used as an order parameter to characterize this free energy difference using umbrella sampling molecular dynamics (MD) simulations with explicit solvent. The constraint imposed along this order parameter allows relatively unrestricted evolution of the intermediate structures away from both canonical A- and B-forms. The free energy difference between the A- and B-forms for the hexamer DNA sequence CTCGAG in aqueous solution is conservatively estimated to be at least 2.8 kcal/mol. A continuum of intermediate structures with no well-defined local minima links the two forms. The absence of any major barriers in the free energy surface is consistent with spontaneous conversion of the A-form DNA to B-form DNA in unconstrained simulations. The extensive sampling in the MD simulations (>0.1 mus) also allowed quantitative energetic characterization of local backbone conformational variables such as sugar pseudorotation angles and BI/BII state equilibria and their dependence on base identity. The absolute minimum in the calculated free energy profile corresponds closely to the crystal structure of the hexamer sequence, indicating that the present method has the potential to identify the most stable state for an arbitrary DNA sequence in water.     

Torsion angle preference and energetics of small-molecule ligands bound to proteins

Small organic molecules can assume conformations in the protein-bound state that are significantly different from those in solution. We have analyzed the conformations of 21 common torsion motifs of small molecules extracted from crystal structures of protein-ligand complexes and compared them with their torsion potentials calculated by an ab initio DFT method. We find a good correlation between the potential energy of the torsion motifs and their conformational distribution in the protein-bound state: The most probable conformations of the torsion motifs agree well with the calculated global energy minima, and the lowest torsion-energy state becomes increasingly dominant as the torsion barrier height increases. The torsion motifs can be divided into 3 groups based on torsion barrier heights: high (>4 kcal/mol), medium (2-4 kcal/mol), and low (<2 kcal/mol). The calculated torsion energy profiles are predictive for the most preferred bound conformation for the high and medium barrier groups, the latter group common in druglike molecules. In the high-barrier group of druglike ligands, >95% of conformational torsions occur in the energy region <4 kcal/mol. The conformations of the torsion motifs in the protein-bound state can be modeled by a Boltzmann distribution with a temperature factor much higher than room temperature. This high-temperature factor, derived by fitting the theoretical model to the experimentally observed conformation occurrence of torsions, can be interpreted as the perturbation that proteins inflict on the conformation of the bound ligand. Using this model, it is calculated that the average strain energy of a torsion motif in ligands bound to proteins is approximately 0.6 kcal/mol, a result which can be related to the lower binding efficiency of larger ligands with more rotatable bonds. The above results indicate that torsion potentials play an important role in dictating ligand conformations in both the free and the bound states.     

Molecular recognition of aminoglycoside antibiotics by bacterial defence proteins: NMR study of the structural and conformational features of streptomycin inactivation by Bacillus subtilis aminoglycoside-6-adenyl transferase

The molecular recognition of streptomycin by Bacillus subtilis aminoglycoside-6-adenyl transferase has been analysed by a combination of NMR techniques and molecular dynamic simulations. This protein inactivates streptomycin by transferring an adenyl group to position six of the streptidine moiety. Our combined approach provides valuable information about the bioactive conformation for both the antibiotic and ATP and shows that the molecular recognition process for streptomycin involves a conformational selection phenomenon. The binding epitope for both ligands has also been analysed by 1D-STD experiments. Finally, the specificity of the recognition process with respect to the aminoglycoside and to the nucleotide has been studied.     

Binding-affinity predictions of HSP90 in the D3R Grand Challenge 2015 with docking,MM/GBSA,QM/MM,and free-energy simulations

We have estimated the binding affinity of three sets of ligands of the heat-shock protein 90 in the D3R grand challenge blind test competition. We have employed four different methods, based on five different crystal structures: first, we docked the ligands to the proteins with induced-fit docking with the Glide software and calculated binding affinities with three energy functions. Second, the docked structures were minimised in a continuum solvent and binding affinities were calculated with the MM/GBSA method (molecular mechanics combined with generalised Born and solvent-accessible surface area solvation). Third, the docked structures were re-optimised by combined quantum mechanics and molecular mechanics (QM/MM) calculations. Then, interaction energies were calculated with quantum mechanical calculations employing 970–1160 atoms in a continuum solvent, combined with energy corrections for dispersion, zero-point energy and entropy, ligand distortion, ligand solvation, and an increase of the basis set to quadruple-zeta quality. Fourth, relative binding affinities were estimated by free-energy simulations, using the multi-state Bennett acceptance-ratio approach. Unfortunately, the results were varying and rather poor, with only one calculation giving a correlation to the experimental affinities larger than 0.7, and with no consistent difference in the quality of the predictions from the various methods. For one set of ligands, the results could be strongly improved (after experimental data were revealed) if it was recognised that one of the ligands displaced one or two water molecules. For the other two sets, the problem is probably that the ligands bind in different modes than in the crystal structures employed or that the conformation of the ligand-binding site or the whole protein changes.     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

Determination of side‐chain‐rotamer and side‐chain and backbone virtual‐bond‐stretching potentials of mean force from AM1 energy surfaces of terminally‐blocked amino‐acid residues,for coarse‐grained simulations of protein structure and folding. I. The method

In this and the accompanying article, we report the development of new physics‐based side‐chain‐rotamer and virtual‐bond‐deformation potentials which now replace the respective statistical potentials used so far in our physics‐based united‐reside UNRES force field for large‐scale simulations of protein structure and dynamics. In this article, we describe the methodology for determining the corresponding potentials of mean force (PMF's) from the energy surfaces of terminally‐blocked amino‐acid residues calculated with the AM1 quantum‐mechanical semiempirical method. The approach is based on minimization of the AM1 energy for fixed values of the angles λ for rotation of the peptide groups about the C^α ··· C^α virtual bonds, and for fixed values of the side‐chain dihedral angles χ, which formed a multidimensional grid. A harmonic‐approximation approach was developed to extrapolate from the energy at a given grid point to other points of the conformational space to compute the respective contributions to the PMF. To test the applicability of the harmonic approximation, the rotamer PMF's of alanine and valine obtained with this approach have been compared with those obtained by using a Metropolis Monte Carlo method. The PMF surfaces computed with the harmonic approximation are more rugged and have more pronounced minima than the MC‐calculated surfaces but the harmonic‐approximation‐and MC‐calculated PMF values are linearly correlated. The potentials derived with the harmonic approximation are, therefore, appropriate for UNRES for which the weights (scaling factors) of the energy terms are determined by force‐field optimization for foldability. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Modification and optimization of the united-residue (UNRES) potential energy function for canonical simulations. I. Temperature dependence of the effective energy function and tests of the optimization method with single training proteins

We report the modification and parametrization of the united-residue (UNRES) force field for energy-based protein structure prediction and protein folding simulations. We tested the approach on three training proteins separately: 1E0L (beta), 1GAB (alpha), and 1E0G (alpha + beta). Heretofore, the UNRES force field had been designed and parametrized to locate native-like structures of proteins as global minima of their effective potential energy surfaces, which largely neglected the conformational entropy because decoys composed of only lowest-energy conformations were used to optimize the force field. Recently, we developed a mesoscopic dynamics procedure for UNRES and applied it with success to simulate protein folding pathways. However, the force field turned out to be largely biased toward -helical structures in canonical simulations because the conformational entropy had been neglected in the parametrization. We applied the hierarchical optimization method, developed in our earlier work, to optimize the force field; in this method, the conformational space of a training protein is divided into levels, each corresponding to a certain degree of native-likeness. The levels are ordered according to increasing native-likeness; level 0 corresponds to structures with no native-like elements, and the highest level corresponds to the fully native-like structures. The aim of optimization is to achieve the order of the free energies of levels, decreasing as their native-likeness increases. The procedure is iterative, and decoys of the training protein(s) generated with the energy function parameters of the preceding iteration are used to optimize the force field in a current iteration. We applied the multiplexing replica-exchange molecular dynamics (MREMD) method, recently implemented in UNRES, to generate decoys; with this modification, conformational entropy is taken into account. Moreover, we optimized the free-energy gaps between levels at temperatures corresponding to a predominance of folded or unfolded structures, as well as to structures at the putative folding-transition temperature, changing the sign of the gaps at the transition temperature. This enabled us to obtain force fields characterized by a single peak in the heat capacity at the transition temperature. Furthermore, we introduced temperature dependence to the UNRES force field; this is consistent with the fact that it is a free-energy and not a potential energy function. beta     

Conformational disorder of conjugated polymers

Conformational disorder of conjugated polymers is an important issue to be understood and quantified. In this paper we present a new method to assess the chain conformation of conjugated polymers based on measurements of intrachain energy transfer. The chain conformation is modeled on the basis of monomer-monomer interactions, such as torsion, bending, and stretching of the connecting bond. The latter two potentials are assumed to be harmonic, while the torsional potential was calculated by density functional theory using B3-LYP functional with the SVP basis set. The energy transfer dynamics of excitons on these chains are quantitatively simulated using Forster-type line-dipole energy transfer. This allows us to compare the simulated ground state conformation of single polymer chains to ultrafast depolarization experiments of poly [3-(2,5-dioctylphenyl)thiophene] in solution. We identify torsional rotation as the main contributor to conformational disorder and find that this disorder is mainly controlled by the energy difference between syn and anti bonds.     

Simulations of lipid crystals: Characterization of potential energy functions and parameters for lecithin molecules

We evaluate an empirical potential energy function and associated parameters for classical molecular dynamics simulations of lecithins, a common class of lipid. The physical accuracy of the force field was tested through its application to molecular dynamics simulations of the known crystal structures of lipid molecules. Average atomic positions and molecular conformation are well maintained during the simulations despite considerable thermal motion. Calculated isotropic temperature factors correlate highly with those from experiment.     

Deciphering common failures in molecular docking of ligand-protein complexes

Common failures in predicting crystal structures of ligand-protein complexes are investigated for three ligand-protein systems by a combined thermodynamic and kinetic analysis of the binding energy landscapes. Misdocked predictions in ligand-protein docking are classified as `soft' and `hard' failures. While a soft failure arises when the search algorithm is unable to find the global energy minimum corresponding to the crystal structure, a hard failure results from a flaw of the energy function to qualify the crystal structure as the predicted lowest energy conformation in docking simulations. We find that neither the determination of a single structure with the lowest energy nor finding the most common binding mode is sufficient to predict crystal structures of the complexes, which belong to the category of hard failures. In a proposed hierarchical approach, structural similarity clustering of the conformations, generated from equilibrium simulations with the simplified energy function, is followed by energy refinement with the AMBER force field. This protocol, that involves a hierarchy of energy functions, resolves some common failures in ligand-protein docking and detects crystallographic binding modes that were not found during docking simulations.     

A new method for side-chain conformation prediction using a Hopfield network and reproduced rotamers

We present a new side-chain prediction method based on energy minimization using a Hopfield network, focusing on the buried residues of proteins. In this method, the network is composed of automata assigned to each rotamer to restrict side-chain conformational space. We reproduced a rotamer library that enabled us to more widely cover the space for side-chain conformations than those previously produced. The accuracy of the side-chain modeling was estimated by three standards: root mean square deviations (rmsds) between the modeled and the crystal structures, the percentages of correctly predicted side-chain torsion angles, and the percentages of correctly predicted hydrogen bonds. The average rmsd for buried side chains of 21 proteins was 1.10 Å. The value was almost always improved relative to the previous works. The percentage of side-chain X₁ angles for buried residues was 87.3%. By considering the hydrogen bond energy, the average percentage of correctly predicted hydrogen bonds rose from 33% without hydrogen bond energy to 52% with the bond energy. We applied this method to homology modeling, where the protein backbone used to predict side-chain conformations deviates from the correct conformation, and could predict side-chain conformations as correctly as those using the correct backbones. © 1996 by John Wiley & Sons, Inc.     

Conformational transition and energy landscape of ErbB4 activated by neuregulin1β: one microsecond molecular dynamics simulations

ErbB4, a receptor tyrosine kinase of the ErbB family, plays crucial roles in cell growth and differentiation, especially in the development of the heart and nervous system. Ligand binding to its extracellular region could modulate the activation process. To understand the mechanism of ErbB4 activation induced by ligand binding, we performed one microsecond molecular dynamics (MD) simulations on the ErbB4 extracellular region (ECR) with and without its endogenous ligand neuregulin1β (NRG1β). The conformational transition of the ECR-ErbB4/NRG1β complex from a tethered inactive conformation to an extended active-like form has been observed, while such large and function-related conformational change has not been seen in the simulation on the ECR-ErbB4, suggesting that ligand binding is indeed the active inducing force for the conformational transition and further dimerization. On the basis of MD simulations and principal component analysis, we constructed a rough energy landscape for the conformational transition of ECR-ErbB4/NRG1β complex, suggesting that the conformational change from the inactive state to active-like state involves a stable conformation. The energy barrier for the tether opening was estimated as ~2.7 kcal/mol, which is very close to the experimental value (1-2 kcal/mol) reported for ErbB1. On the basis of the simulation results, an atomic mechanism for the ligand-induced activation of ErbB4 was postulated. The present MD simulations provide a new insight into the conformational changes underlying the activation of ErbB4.     

Conformational analysis: crystallographic,NMR, and molecular mechanics studies of flexible sulfonic esters

Two novel X-ray structures of the sulfonic ester derivatives 2-(6-iodo-1,3-benzodioxol-5-yl)ethyl 4-nitrobenzenesulfonate, 3, and 2-(6-iodo-1,3-benzodioxol-5-yl)ethyl 4-methylbenzenesulfonate, 4, have been obtained in a study aimed at analyzing the structures and conformations of sulfonic ester derivatives that are routinely used in alkaloid syntheses. The crystal structure of 4 is highly unusual, containing four independent molecules that belong to two distinct conformational types: (1) a hairpin conformation (stabilized mainly by intramolecular pi-stacking) and (2) a stepped conformation (stabilized mainly by intermolecular pi-stacking). Compound 3, on the other hand, crystallizes exclusively as the hairpin conformer. New MM+ force field parameters for sulfonic esters have been developed using the X-ray data, empirical rules, and DFT calculations to estimate the bond dipole parameters. Grid searches of conformational space for 3 and 4 using MM methods show that there are several gas-phase conformations within 5 kcal/mol of the global minimum and that the lowest energy conformations (by approximately 4.6 kcal/mol) are of the hairpin type. Analysis of the MM conformational energies suggests that the dominant intramolecular interaction stabilizing the hairpin conformations of 3 and 4 is van der Waals attraction. Moreover, the lattice energies for packing the hairpin conformations of 3 and 4 are approximately 4 kcal/mol more favorable than for the stepped conformations. Various intermolecular interactions contribute to the complexity of the observed crystal structures of 3 and 4, including electrostatic attraction between O and I atoms in neighboring molecules. Langevin dynamics (LD) simulations at several temperatures (6.0 ns, friction coefficient = 2.5 ps(-1)) indicate that the conformational exchange rates are approximately 10(10)-10(11) s(-1) over the temperature range 213-400 K, accounting for the temperature-independent (1)H NMR spectra of 3 and 4.     

Conformational heterogeneity observed in simulations of a pyrene-substituted DNA

NMR studies previously carried out for a DNA system with a pyrene-substituted base did not observe NOEs involving the adenine located 5' to the pyrene, and thus the conformation of the adenine was poorly defined in the resulting family of refined structures. However, chemical shift data suggested that an AT base pair may be present. We have carried out fully unrestrained molecular dynamics simulations starting from several members of the family of structures, and these simulations support the existence of an AT base pair for this region. Simulations in both explicit and implicit solvent were carried out, with each converging to either anti or syn conformation for adenine and base pairing in all cases. During these simulations, large and dramatic conformational changes are observed that suggest pathways for complex conformational changes in the highly packed DNA interior. Our analysis reveals little difference in the energies of these syn and anti conformations, in contrast to control calculations carried out for standard DNA (in the absence of a neighboring pyrene). While no interconversion between the conformations was observed in standard simulations, reversible anti/syn exchange was directly simulated using the locally enhanced sampling approach. No exchange was seen in the non-pyrene control sequence. Together, these results suggest that an increased flexibility is introduced as a consequence of the pyrene substitution, offering an explanation that is consistent with the available experimental data. These results increase our optimism that simulations in atomic detail may provide accurate models for experimental observations in complex systems.     

Cooperative Effect of a Classical and a Weak Hydrogen Bond for the Metal‐Induced Construction of a Self‐Assembled β‐Turn Mimic

A novel metal‐induced template for the self‐assembly of two independent phosphane ligands by means of unprecedented multiple noncovalent interactions (classical hydrogen bond, weak hydrogen bond, metal coordination, π‐stacking interaction) was developed and investigated. Our results address the importance and capability of weak hydrogen bonds (WHBs) as important attractive interactions in self‐assembling processes based on molecular recognition. Together with a classical hydrogen bond, WHBs may serve as promoters for the specific self‐assembly of complementary monomeric phosphane ligands into supramolecular hybrid structures. The formation of an intermolecular C? H???N hydrogen bond and its persistence in the solid state and in solution was studied by X‐ray crystal analysis, mass spectrometry and NMR spectroscopy analysis. Further evidence was demonstrated by DFT calculations, which gave specific geometric parameters for the proposed conformations and allowed us to estimate the energy involved in the hydrogen bonds that are responsible for the molecular recognition process. The presented template can be regarded as a new type of self‐assembled β‐turn mimic or supramolecular pseudo amino acid for the nucleation of β‐sheet structures when attached to oligopeptides.     

Molecular dynamics simulation of the A-DNA to B-DNA transition in aqueous RbCl solution

Unrestrained molecular dynamics (MD) simulations have been carried out to characterize the stability of DNA conformations and the dynamics of A-DNA→B-DNA conformational transitions in aqueous RbCl solutions. The PARM99 force field in the AMBER8 package was used to investigate the effect of RbCl concentration on the dynamics of the A→B conformational transition in the DNA duplex d(CGCGAATTCGCG)2 . Canonical Aand B-form DNA were assumed for the initial conformation and the final conformation had a length per complete turn that matched the canonical B-DNA. The DNA structure was monitored for 3.0 ns and the distances between the C5′ atoms were obtained from the simulations. It was found that all of the double stranded DNA strands of A-DNA converged to the structure of B-form DNA within 1.0 ns during the unrestrained MD simulations. In addition, increasing the RbCl concentration in aqueous solution hindered the A→B conformational transition and the transition in aqueous RbCl solution was faster than that in aqueous NaCl solution for the same electrolyte strength. The effects of the types and concentrations of counterions on the dynamics of the A→B conformational transition can be understood in terms of the variation in water activity and the number of accumulated counterions in the major grooves of A-DNA. The rubidium ion distributions around both fixed A-DNA and B-DNA were obtained using the restrained MD simulations to help explain the effect of RbCl concentration on the dynamics of the A→B conformational transition.     

Molecular mechanics parameters for the FapydG DNA lesion

FapydG is a common oxidative DNA lesion involving opening of the imidazole ring. It shares the same precursor as 8-oxodG and can be excised by the same enzymes as 8-oxodG. However, the loss of the aromatic imidazole in FapydG results in a reduction of the double bond character between C5 and N7, with an accompanying increase in conformational flexibility. Experimental characterization of FapydG is hampered by high reactivity, and thus it is desirable to investigate structural details through computer simulation. We show that the existing Amber force field parameters for FapydG do not reproduce X-ray structural data. We employed quantum mechanics energy profile calculations to derive new molecular mechanics parameters for the rotation of the dihedral angles in the eximidazole moiety. Using these parameters, all-atom simulations in explicit water reproduce the nonplanar conformation of cFapydG in the crystal structure of the complex with L. lactis glycosylase Fpg. We note that the nonplanar structure is stabilized by an acidic residue that is not present in most Fpg sequences. Simulations of the E-->S mutant, as present in E. coli, resulted in a more planar conformation, suggesting that the highly nonplanar form observed in the crystal structure may not have direct biological relevance for FapydG.