Flexible protein‐ligand docking using the Fleksy protocol

Considering protein plasticity is important in accurately predicting the three‐dimensional geometry of protein‐ligand complexes. Here, we present the first public release of our flexible docking tool Fleksy, which is able to consider both ligand and protein flexibility in the docking process. We describe the workflow and different features of the software and present its performance on two cross‐docking benchmark datasets. © 2012 Wiley Periodicals, Inc.     

Exploring hierarchical refinement techniques for induced fit docking with protein and ligand flexibility

We present a series of molecular‐mechanics‐based protein refinement methods, including two novel ones, applied as part of an induced fit docking procedure. The methods used include minimization; protein and ligand sidechain prediction; a hierarchical ligand placement procedure similar to a‐priori protein loop predictions; and a minimized Monte Carlo approach using normal mode analysis as a move step. The results clearly indicate the importance of a proper opening of the active site backbone, which might not be accomplished when the ligand degrees of freedom are prioritized. The most accurate method consisted of the minimized Monte Carlo procedure designed to open the active site followed by a hierarchical optimization of the sidechain packing around a mobile flexible ligand. The methods have been used on a series of 88 protein‐ligand complexes including both cross‐docking and apo‐docking members resulting in complex conformations determined to within 2.0 Å heavy‐atom RMSD in 75% of cases where the protein backbone rearrangement upon binding is less than 1.0 Å α‐carbon RMSD. We also demonstrate that physics‐based all‐atom potentials can be more accurate than docking‐style potentials when complexes are sufficiently refined. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Flexible CDOCKER: Development and application of a pseudo‐explicit structure‐based docking method within CHARMM

Protein‐ligand docking is a commonly used method for lead identification and refinement. While traditional structure‐based docking methods represent the receptor as a rigid body, recent developments have been moving toward the inclusion of protein flexibility. Proteins exist in an interconverting ensemble of conformational states, but effectively and efficiently searching the conformational space available to both the receptor and ligand remains a well‐appreciated computational challenge. To this end, we have developed the Flexible CDOCKER method as an extension of the family of complete docking solutions available within CHARMM. This method integrates atomically detailed side chain flexibility with grid‐based docking methods, maintaining efficiency while allowing the protein and ligand configurations to explore their conformational space simultaneously. This is in contrast to existing approaches that use induced‐fit like sampling, such as Glide or Autodock, where the protein or the ligand space is sampled independently in an iterative fashion. Presented here are developments to the CHARMM docking methodology to incorporate receptor flexibility and improvements to the sampling protocol as demonstrated with re‐docking trials on a subset of the CCDC/Astex set. These developments within CDOCKER achieve docking accuracy competitive with or exceeding the performance of other widely utilized docking programs. © 2015 Wiley Periodicals, Inc.     

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     

Efficient inclusion of receptor flexibility in grid-based protein-ligand docking

Accounting for receptor flexibility is an essential component of successful protein-ligand docking but still marks a major computational challenge. For many target molecules of pharmaceutical relevance, global backbone conformational changes are relevant during the ligand binding process. However, popular methods that represent the protein receptor molecule as a potential grid typically assume a rigid receptor structure during ligand-receptor docking. A new approach has been developed that combines inclusion of global receptor flexibility with the efficient potential grid representation of the receptor molecule. This is achieved using interpolation between grid representations of the receptor protein deformed in selected collective degrees of freedom. The method was tested on the docking of three ligands to apo protein kinase A (PKA), an enzyme that undergoes global structural changes upon inhibitor binding. Structural variants of PKA were generated along the softest normal mode of an elastic network representation of apo PKA. Inclusion of receptor deformability during docking resulted in a significantly improved docking performance compared with rigid PKA docking, thus allowing for systematic virtual screening applications at small additional computational cost.     

Fast and accurate grid representations for atom‐based docking with partner flexibility

Macromolecular docking methods can broadly be divided into geometric and atom‐based methods. Geometric methods use fast algorithms that operate on simplified, grid‐like molecular representations, while atom‐based methods are more realistic and flexible, but far less efficient. Here, a hybrid approach of grid‐based and atom‐based docking is presented, combining precalculated grid potentials with neighbor lists for fast and accurate calculation of atom‐based intermolecular energies and forces. The grid representation is compatible with simultaneous multibody docking and can tolerate considerable protein flexibility. When implemented in our docking method ATTRACT, grid‐based docking was found to be ∼35x faster. With the OPLSX forcefield instead of the ATTRACT coarse‐grained forcefield, the average speed improvement was >100x. Grid‐based representations may allow atom‐based docking methods to explore large conformational spaces with many degrees of freedom, such as multiple macromolecules including flexibility. This increases the domain of biological problems to which docking methods can be applied. © 2017 Wiley Periodicals, Inc.     

GalaxyDock2: Protein–ligand docking using beta‐complex and global optimization

In this article, an enhanced version of GalaxyDock protein–ligand docking program is introduced. GalaxyDock performs conformational space annealing (CSA) global optimization to find the optimal binding pose of a ligand both in the rigid‐receptor mode and the flexible‐receptor mode. Binding pose prediction has been improved compared to the earlier version by the efficient generation of high‐quality initial conformations for CSA using a predocking method based on a beta‐complex derived from the Voronoi diagram of receptor atoms. Binding affinity prediction has also been enhanced by using the optimal combination of energy components, while taking into consideration the energy of the unbound ligand state. The new version has been tested in terms of binding mode prediction, binding affinity prediction, and virtual screening on several benchmark sets, showing improved performance over the previous version and AutoDock, on which the GalaxyDock energy function is based. GalaxyDock2 also performs better than or comparable to other state‐of‐the‐art docking programs. GalaxyDock2 is freely available at http://galaxy.seoklab.org/softwares/galaxydock.html . © 2013 Wiley Periodicals, Inc.     

Including ligand‐induced protein flexibility into protein tunnel prediction

In proteins with buried active sites, understanding how ligands migrate through the tunnels that connect the exterior of the protein to the active site can shed light on substrate specificity and enzyme function. A growing body of evidence highlights the importance of protein flexibility in the binding site on ligand binding; however, the influence of protein flexibility throughout the body of the protein during ligand entry and egress is much less characterized. We have developed a novel tunnel prediction and evaluation method named IterTunnel, which includes the influence of ligand‐induced protein flexibility, guarantees ligand egress, and provides detailed free energy information as the ligand proceeds along the egress route. IterTunnel combines geometric tunnel prediction with steered molecular dynamics in an iterative process to identify tunnels that open as a result of ligand migration and calculates the potential of mean force of ligand egress through a given tunnel. Applying this new method to cytochrome P450 2B6, we demonstrate the influence of protein flexibility on the shape and accessibility of tunnels. More importantly, we demonstrate that the ligand itself, while traversing through a tunnel, can reshape tunnels due to its interaction with the protein. This process results in the exposure of new tunnels and the closure of preexisting tunnels as the ligand migrates from the active site. © 2014 Wiley Periodicals, Inc.     

Protein‐specific force field derived from the fragment molecular orbital method can improve protein–ligand binding interactions

Accurate computational estimate of the protein–ligand binding affinity is of central importance in rational drug design. To improve accuracy of the molecular mechanics (MM) force field (FF) for protein–ligand simulations, we use a protein‐specific FF derived by the fragment molecular orbital (FMO) method and by the restrained electrostatic potential (RESP) method. Applying this FMO‐RESP method to two proteins, dodecin, and lysozyme, we found that protein‐specific partial charges tend to differ more significantly from the standard AMBER charges for isolated charged atoms. We did not see the dependence of partial charges on the secondary structure. Computing the binding affinities of dodecin with five ligands by MM PBSA protocol with the FMO‐RESP charge set as well as with the standard AMBER charges, we found that the former gives better correlation with experimental affinities than the latter. While, for lysozyme with five ligands, both charge sets gave similar and relatively accurate estimates of binding affinities. © 2013 Wiley Periodicals, Inc.     

Insights into docking and scoring neuronal α4β2 nicotinic receptor agonists using molecular dynamics simulations and QM/MM calculations

A combined quantum mechanical (QM)‐polarized docking and molecular dynamics approach to study the binding mode and to predict the binding affinity of ligands acting at the α4β2‐nAChR is presented. The results obtained in this study indicate that the quantum mechanical/molecular mechanics docking protocol well describes the charge‐driven interactions occurring in the binding of nicotinic agonists, and it is able to represent the polarization effects on the ligand exerted by the surrounding atoms of the receptor at the binding site. This makes it possible to properly score agonists of α4β2‐nAChR and to reproduce the experimental binding affinity data with good accuracy, within a mean error of 2.2 kcal/mol. Moreover, applying the QM‐polarized docking to an ensemble of nAChR conformations obtained from MD simulations enabled us to accurately capture nAChR‐ligand induced‐fit effects. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Pharmacophore Modeling and in Silico Screening Studies to Design Potential KDR Kinase Inhibitors

A novel ligand‐based pharmacophore model for KDR kinase was generated on the basis of chemical features of 30 KDR kinase inhibitors. This pharmacophore model consists of one hydrogen‐bond acceptor, one hydrogen‐bond donor and two hydrophobic groups. Several methods have been used to validate the model, suggesting that it can serve as a reliable tool for virtual screening to facilitate the discovery of novel KDR inhibitors. The model was then used as database search query from the National Cancer Institute (NCI) database for the rational design to identify new hit compound.     

Identification of Novel Dopamine D2 Receptor Ligands—A Combined In Silico/In Vitro Approach

Diseases of the central nervous system are an alarming global problem showing an increasing prevalence. Dopamine receptor D₂ (D₂R) has been shown to be involved in central nervous system diseases. While different D₂R-targeting drugs have been approved by the FDA, they all suffer from major drawbacks due to promiscuous receptor activity leading to adverse effects. Increasing the number of potential D₂R-targeting drug candidates bears the possibility of discovering molecules with less severe side-effect profiles. In dire need of novel D₂R ligands for drug development, combined in silico/in vitro approaches have been shown to be efficient strategies. In this study, in silico pharmacophore models were generated utilizing both ligand- and structure-based approaches. Subsequently, different databases were screened for novel D₂R ligands. Selected virtual hits were investigated in vitro, quantifying their binding affinity towards D₂R. This workflow successfully identified six novel D₂R ligands exerting micro- to nanomolar (most active compound K_I = 4.1 nM) activities. Thus, the four pharmacophore models showed prospective true-positive hit rates in between 4.5% and 12%. The developed workflow and identified ligands could aid in developing novel drug candidates for D₂R-associated pathologies.     

Modeling bacteriophage amplification as a predictive tool for optimized MALDI‐TOF MS‐based bacterial detection

Matrix‐assisted laser desorption/ionization time‐of‐flight mass spectrometry (MALDI‐TOF MS) is a valuable tool for rapid bacterial detection and identification but is limited by the need for relatively high cell count samples, which have been grown under strictly controlled conditions. These requirements can be eliminated by the natural infection of a viable bacterial species of interest with a host‐specific phage. This produces a rapid increase in phage protein concentrations in comparison to bacterial concentrations, which can in turn be exploited as a method for signal amplification during MALDI‐TOF MS. One drawback to this approach is the requirement for repetitive, time‐consuming sample preparation and analysis applied over the course of a phage infection to monitor phage concentrations as a function of time to determine the MALDI‐TOF MS detection limit. To reduce the requirement for repeated preparation and analysis, a modified phage therapy model was investigated as a means for predicting the time during a given phage infection when a detectable signal would occur. The modified model used a series of three differential equations composed of predetermined experimental parameters including phage burst size and burst time to predict progeny phage concentrations as a function of time. Using Yersinia pestis with plague diagnostic phage ?A1122 and Escherichia coli with phage MS2 as two separate, well‐characterized model phage–host pairs, we conducted in silico modeling of the infection process and compared it with experimental infections monitored in real time by MALDI‐TOF MS. Significant agreement between mathematically calculated phage growth curves and those experimentally obtained by MALDI‐TOF MS was observed, thus verifying this method's utility for significant time and labor reduction. Copyright © 2012 John Wiley & Sons, Ltd.     

Twist‐dependent stacking energy of base‐pair steps in B‐DNA geometry: A density functional theory approach

Stacking energy of all the 10 unique DNA base‐pair steps (bp step) are calculated using density functional theory within the ultrasoft pseudopotential plane wave method and local density approximation for the exchange‐correlation functional. We have studied the dependence of stacking energy on twist angle, an aspect found difficult to explain using classical theory. We have found that the twist angle for different bp steps at stacking energy minimum matches extremely well with the values of average twist obtained from B‐DNA crystal structure data. This indicates that the use of a proper quantum chemical method to calculate the π‐π electronic interactions may explain stacking energy without incorporating hydrophobic interaction through solvent or effect of backbone through pseudobond. From the twist angle‐dependent stacking energy profile, we have also generated the probability distributions of twist for all the bp steps and calculated the variance of the distribution. Our calculated variances show similar trend to that of the experimental data for which sufficient numbers of data are available. The TA, AT, and CG doublets show large variances among the 10 possible bp steps, indicating their maximum flexibility. This might be the case of unusual deformation observed at the TATA‐box while binding to TBP protein. © 2008 Wiley Periodicals, Inc. Int J Quantum Chem, 2008     

Multidimensional umbrella sampling and replica‐exchange molecular dynamics simulations for structure prediction of transmembrane helix dimers

Structural information of a transmembrane (TM) helix dimer is useful in understanding molecular mechanisms of important biological phenomena such as signal transduction across the cell membrane. Here, we describe an umbrella sampling (US) scheme for predicting the structure of a TM helix dimer in implicit membrane using the interhelical crossing angle and the TM–TM relative rotation angles as the reaction coordinates. This scheme conducts an efficient conformational search on TM–TM contact interfaces, and its robustness is tested by predicting the structures of glycophorin A (GpA) and receptor tyrosine kinase EphA1 (EphA1) TM dimers. The nuclear magnetic resonance (NMR) structures of both proteins correspond to the global free‐energy minimum states in their free‐energy landscapes. In addition, using the landscape of GpA as a reference, we also examine the protocols of temperature replica‐exchange molecular dynamics (REMD) simulations for structure prediction of TM helix dimers in implicit membrane. A wide temperature range in REMD simulations, for example, 250–1000 K, is required to efficiently obtain a free‐energy landscape consistent with the US simulations. The interhelical crossing angle and the TM–TM relative rotation angles can be used as reaction coordinates in multidimensional US and be good measures for conformational sampling of REMD simulations. © 2013 Wiley Periodicals, Inc.     

Functional assessment of “in vivo” and “in silico” mutations in the guanine binding site of RNase T1: A DFT study

Experimentally, the functional assessment of amino acid side chains in proteins is carried out by comparing parameters such as binding constants for the wild‐type protein and a mutant protein in which the considered side chain is deleted. In the present study, we apply a density functional theory (DFT) methodology to obtain changes in binding energy upon mutations in the enzyme ribonuclease T₁. Mutant structures were either taken directly from crystallographic data (“in vivo”) allowing for conformational changes upon mutation, or derived from the wild‐type (“in silico”). Excluding entropic contributions, the computed interaction energy changes upon mutation in vivo correlate qualitatively well with experimental binding free energy changes. In contrast, the in silico approach does not perform as well, especially for residues that contribute largely to binding. Subsequently, we assessed the applicability of the in vivo approach by analyzing the functional cooperativity between pairs of side chains. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

The solid‐state structure of the β‐blocker metoprolol: a combined experimental and in silico investigation

Metoprolol {systematic name: (RS)‐1‐isopropylamino‐3‐[4‐(2‐methoxyethyl)phenoxy]propan‐2‐ol}, C₁₅H₂₅NO₃, is a cardioselective β₁‐adrenergic blocking agent that shares part of its molecular skeleton with a large number of other β‐blockers. Results from its solid‐state characterization by single‐crystal and variable‐temperature powder X‐ray diffraction and differential scanning calorimetry are presented. Its molecular and crystal arrangements have been further investigated by molecular modelling, by a Cambridge Structural Database (CSD) survey and by Hirshfeld surface analysis. In the crystal, the side arm bearing the isopropyl group, which is common to other β‐blockers, adopts an all‐trans conformation, which is the most stable arrangement from modelling data. The crystal packing of metoprolol is dominated by an O—H…N/N…H—O pair of hydrogen bonds (as also confirmed by a Hirshfeld surface analysis), which gives rise to chains containing alternating R and S metoprolol molecules extending along the b axis, supplemented by a weaker O…H—N/N—H…O pair of interactions. In addition, within the same stack of molecules, a C—H…O contact, partially oriented along the b and c axes, links homochiral molecules. Amongst the solid‐state structures of molecules structurally related to metoprolol deposited in the CSD, the β‐blocker drug betaxolol shows the closest analogy in terms of three‐dimensional arrangement and interactions. Notwithstanding their close similarity, the crystal lattices of the two drugs respond differently on increasing temperature: metoprolol expands anisotropically, while for betaxolol, an isotropic thermal expansion is observed.