Towards understanding the mechanisms of molecular recognition by computer simulations of ligand-protein interactions

Ligand binding free energy evaluation by Monte Carlo Recursion

The correct evaluation of ligand binding free energies by computational methods is still a very challenging active area of research. The most employed methods for these calculations can be roughly classified into four groups: (i) the fastest and less accurate methods, such as molecular docking, designed to sample a large number of molecules and rapidly rank them according to the potential binding energy; (ii) the second class of methods use a thermodynamic ensemble, typically generated by molecular dynamics, to analyze the endpoints of the thermodynamic cycle for binding and extract differences, in the so-called 'end-point' methods; (iii) the third class of methods is based on the Zwanzig relationship and computes the free energy difference after a chemical change of the system (alchemical methods); and (iv) methods based on biased simulations, such as metadynamics, for example. These methods require increased computational power and as expected, result in increased accuracy for the determination of the strength of binding. Here, we describe an intermediate approach, based on the Monte Carlo Recursion (MCR) method first developed by Harold Scheraga. In this method, the system is sampled at increasing effective temperatures, and the free energy of the system is assessed from a series of terms W(b,T), computed from Monte Carlo (MC) averages at each iteration. We show the application of the MCR for ligand binding with datasets of guest-hosts systems (N = 75) and we observed that a good correlation is obtained between experimental data and the binding energies computed with MCR. We also compared the experimental data with an end-point calculation from equilibrium Monte Carlo calculations that allowed us to conclude that the lower-energy (lower-temperature) terms in the calculation are the most relevant to the estimation of the binding energies, resulting in similar correlations between MCR and MC data and the experimental values. On the other hand, the MCR method provides a reasonable view of the binding energy funnel, with possible connections with the ligand binding kinetics, as well. The codes developed for this analysis are publicly available on GitHub as a part of the LiBELa/MCLiBELa project (https://github.com/alessandronascimento/LiBELa).

Comparative Analysis of Electrostatic Models for Ligand Docking

The precise modeling of molecular interactions remains an important goal among molecular modeling techniques. Some of the challenges in the field include the precise definition of a Hamiltonian for biomolecular systems, together with precise parameters derived from Molecular Mechanics Force Fields, for example. The problem is even more challenging when interaction energies from different species are computed, such as the interaction energy involving a ligand and a protein, given that small differences must be computed from large energies. Here we evaluated the effects of the electrostatic model for ligand binding energy evaluation in the context of ligand docking. For this purpose, a classical Coulomb potential with distance-dependent dielectrics was compared with a Poisson-Boltzmann (PB) model for electrostatic potential computation, based on DelPhi calculations. We found that, although the electrostatic energies were highly correlated for the Coulomb and PB models, the ligand pose and the enrichment of actual ligands against decoy compounds, were improved when binding energies were computed using PB as compared to the Coulomb model. We observed that the electrostatic energies computed with the Coulomb model were, on average, ten times larger than the energies computed with the PB model, suggesting a strong overestimation of the polar interactions in the Coulomb model. We also found that a slightly smoothed Lennard-Jones potential combined with the PB model resulted in a good compromise between ligand sampling and energetic scoring.

Towards a critical evaluation of an empirical and volume-based solvation function for ligand docking

Molecular docking is an important tool for the discovery of new biologically active molecules given that the receptor structure is known. An excellent environment for the development of new methods and improvement of the current methods is being provided by the rapid growth in the number of proteins with known structure. The evaluation of the solvation energies outstands among the challenges for the modeling of the receptor-ligand interactions, especially in the context of molecular docking where a fast, though accurate, evaluation is ought to be achieved. Here we evaluated a variation of the desolvation energy model proposed by Stouten (Stouten P.F.W. et al, Molecular Simulation, 1993, 10: 97–120), or SV model. The SV model showed a linear correlation with experimentally determined solvation energies, as available in the database FreeSolv. However, when used in retrospective docking simulations using the benchmarks DUD, charged-matched DUD and DUD-Enhanced, the SV model resulted in poorer enrichments when compared to a pure force field model with no correction for solvation effects. The data provided here is consistent with other empirical solvation models employed in the context of molecular docking and indicates that a good model to account for solvent effects is still a goal to achieve. On the other hand, despite the inability to improve the enrichment of retrospective simulations, the SV solvation model showed an interesting ability to reduce the number of molecules with net charge -2 and -3 e among the top-scored molecules in a prospective test.

Molecular Dynamics Assisted Mechanistic Study of Isoniazid-Resistance against Mycobacterium tuberculosis InhA

Examination of InhA mutants I16T, I21V, I47T, S94A, and I95P showed that direct and water mediated H-bond interactions between NADH and binding site residues reduced drastically. It allowed conformational flexibility to NADH, particularly at the pyrophosphate region, leading to weakening of its binding at dinucleotide binding site. The highly scattered distribution of pyrophosphate dihedral angles and chi1 side chain dihedral angles of corresponding active site residues therein confirmed weak bonding between InhA and NADH. The average direct and water mediated bridged H-bond interactions between NADH and mutants were observed weaker as compared to wild type. Further, estimated NADH binding free energy in mutants supported the observed weakening of InhA-NADH interactions. Similarly, per residue contribution to NADH binding was also found little less as compared to corresponding residues in wild type. This investigation clearly depicted and supported the effect of mutations on NADH binding and can be accounted for isoniazid resistance as suggested by previous biochemical and mutagenic studies. Further, structural analysis of InhA provided the crucial points to enhance the NADH binding affinity towards InhA mutants in the presence of direct InhA inhibitors to combat isoniazid drug resistance. This combination could be a potential alternative for treatment of drug resistant tuberculosis.

Can the energy gap in the protein-ligand binding energy landscape be used as a descriptor in virtual ligand screening?

The ranking of scores of individual chemicals within a large screening library is a crucial step in virtual screening (VS) for drug discovery. Previous studies showed that the quality of protein-ligand recognition can be improved using spectrum properties and the shape of the binding energy landscape. Here, we investigate whether the energy gap, defined as the difference between the lowest energy pose generated by a docking experiment and the average energy of all other generated poses and inferred to be a measure of the binding energy landscape sharpness, can improve the separation power between true binders and decoys with respect to the use of the best docking score. We performed retrospective single- and multiple-receptor conformation VS experiments in a diverse benchmark of 40 domains from 38 therapeutically relevant protein targets. Also, we tested the performance of the energy gap on 36 protein targets from the Directory of Useful Decoys (DUD). The results indicate that the energy gap outperforms the best docking score in its ability to discriminate between true binders and decoys, and true binders tend to have larger energy gaps than decoys. Furthermore, we used the energy gap as a descriptor to measure the height of the native binding phase and obtained a significant increase in the success rate of near native binding pose identification when the ligand binding conformations within the boundaries of the native binding phase were considered. The performance of the energy gap was also evaluated on an independent test case of VS-identified PKR-like ER-localized eIF2α kinase (PERK) inhibitors. We found that the energy gap was superior to the best docking score in its ability to more highly rank active compounds from inactive ones. These results suggest that the energy gap of the protein-ligand binding energy landscape is a valuable descriptor for use in VS.

Integrating ligand-based and protein-centric virtual screening of kinase inhibitors using ensembles of multiple protein kinase genes and conformations

The rapidly growing wealth of structural and functional information about kinase genes and kinase inhibitors that is fueled by a significant therapeutic role of this protein family provides a significant impetus for development of targeted computational screening approaches. In this work, we explore an ensemble-based, protein-centric approach that allows for simultaneous virtual ligand screening against multiple kinase genes and multiple kinase receptor conformations. We systematically analyze and compare the results of ligand-based and protein-centric screening approaches using both single-receptor and ensemble-based docking protocols. A panel of protein kinase targets that includes ABL, EGFR, P38, CDK2, TK, and VEGFR2 kinases is used in this comparative analysis. By applying various performance metrics we have shown that ligand-centric shape matching can provide an effective enrichment of active compounds outperforming single-receptor docking screening. However, ligand-based approaches can be highly sensitive to the choice of inhibitor queries. Employment of multiple inhibitor queries combined with parallel selection ranking criteria can improve the performance and efficiency of ligand-based virtual screening. We also demonstrated that replica-exchange Monte Carlo docking with kinome-based ensembles of multiple crystal structures can provide a superior early enrichment on the kinase targets. The central finding of this study is that incorporation of the template-based structural information about kinase inhibitors and protein kinase structures in diverse functional states can significantly enhance the overall performance and robustness of both ligand and protein-centric screening strategies. The results of this study may be useful in virtual screening of kinase inhibitors potentially offering a beneficial spectrum of therapeutic activities across multiple disease states.

Rational design of affinity peptide ligand by flexible docking simulation

Rational design of affinity peptide ligands of proteins by flexible docking simulation is performed using the SYBYL program package. This approach involves the use of experimental data to verify a scoring function that can be used to assess the affinity of a peptide for its target protein. The enzyme-linked immunosorbent assay (ELISA) data of several peptides displayed on phage surfaces for insulin and lysozyme, respectively, reported in literature are used for the purpose. It is found that the absolute values of the Dscore calculated from the docking correspond well to the ELISA data that relate to the affinity between the peptides and the target molecule. So, the Dscore function is used to assess the affinity of docked peptides in a pentapeptide library designed on the basis of protein (alpha-amylase) structure. As a result, a pentapeptide with a high Dscore value is selected and a hexapeptide (FHENWS) is built by linking serine to its C-terminal to lengthen the peptide. Molecular surface analysis with the MOLCAD program reveals that electrostatic interactions (including hydrogen bonds) and Van der Waals forces contribute to the affinity of the hexapeptide for alpha-amylase. Chromatographic experiments with the immobilized peptide have given further evidence for this observation. Adsorption isotherm described by the Langmuir equation indicates that the apparent binding constant of alpha-amylase to the immobilized hexapeptide was 2.5x10(5)L/mol. Finally, high affinity and specificity of the affinity adsorbent is exemplified by the purification of alpha-amylase from crude fermentation broth of Bacillus subtilis.

Activity of Dendrimer−Methotrexate Conjugates on Methotrexate-Sensitive and -Resistant Cell Lines

Dendritic nanostructures can play a key role in drug delivery, due to the high density and variety of surface functional groups that can facilitate and modulate the delivery process. We have investigated the effect of dendrimer end-functionality on the activity of polyamido amine (PAMAM) dendrimer-methotrexate (MTX) conjugates in MTX-sensitive and MTX-resistant human acute lymphoblastoid leukemia (CCRF-CEM) and Chinese hamster ovary (CHO) cell lines. Two amide-bonded PAMAM dendrimer-MTX conjugates were prepared using a dicyclohexylcarbodiimide (DCC) coupling reaction: one between a carboxylic acid-terminated G2.5 dendrimer and the amine groups of the MTX (conjugate A) and another between an amine-terminated G3 dendrimer and the carboxylic acid group of the MTX (conjugate B). Our studies suggest that conjugate A showed an increased drug activity compared to an equimolar amount of free MTX toward both sensitive and resistant cell lines, whereas conjugate B did not show significant activity on any of the cell lines. Despite substantially impaired MTX transport by MTX-resistant CEM/MTX and RII cells, conjugate A showed sensitivity increases of approximately 8- and 24-fold (based on IC50 values), respectively, compared to free MTX. Co-incubation of the cells with adenosine and thymidine along with either conjugate A or MTX resulted in almost complete protection, suggesting that the conjugate achieves its effect on dihyrofolate reductase (DHFR) enzyme through the same mechanism as that of MTX. The differences in cytotoxicity of these amide-bonded conjugates may be indicative of differences in the intracellular drug release from the cationic dendrimer (conjugate B) versus the anionic dendrimer (conjugate A), perhaps due to the differences in lysosomal residence times dictated by the surface functionality. These findings demonstrate the feasibility of using dendrimers as drug delivery vehicles for achieving higher therapeutic effects in chemotherapy, especially in drug-resistant cells.

In silico fragment-based discovery of DPP-IV S1 pocket binders

Dipeptidyl peptidase IV is a clinically validated target for type-2 diabetes and belongs to a family of peptidases with a quite unique post-proline cleavage specificity. Known inhibitors contain a limited number of molecular anchors occupying the small prototypical S1 pocket. A virtual screening approach for such S1-binding fragments was carried out using FlexX docking to evaluate its potential to confirm known and find novel compounds. Several low molecular weight inhibitors exhibiting activities in the micromolar range could be identified as starting points for structure-based design.

Computational analysis of ligand binding dynamics at the intermolecular hot spots with the aid of simulated tempering and binding free energy calculations

Equilibrium binding dynamics is studied for a panel of benzimidazole-containing compounds at the remodeled interface between human growth hormone (hGH) and the extracellular domain of its receptor (hGHbp), engineered by mutating to glycine hot spot residues T175 from the hormone and W104 from the receptor. Binding energetics is predicted in a good agreement with the experimental data for a panel of these small molecules that complement the engineered defect and restore the binding affinity of the wild-type hGH-hGHbp complex. The results of simulated tempering ligand dynamics at the protein-protein interface reveals a diversity of ligand binding modes that is consistent with the structural orientation of the benzimidazole ring which closely mimics the position of the mutated W104 hot spot residue in the wild-type hGH-hGHbp complex. This structural positioning of the benzimidazole core motif is shown to be a critical feature of the low-energy ligand conformations binding in the engineered cavity. The binding free energy analysis provides a plausible rationale behind the experimental dissociation constants and suggests a key role of ligand-protein van der Waals interactions in restoring binding affinity.

Computational detection of the binding-site hot spot at the remodeled human growth hormone-receptor interface

A hierarchical computational approach is used to identify the engineered binding-site cavity at the remodeled intermolecular interface between the mutants of human growth hormone (hGH) and the extracellular domain of its receptor (hGHbp). Multiple docking simulations are conducted with the remodeled hGH-hGHbp complex for a panel of potent benzimidazole-containing inhibitors that can restore the binding affinity of the wild-type complex, and for a set of known nonactive small molecules that contain different heterocyclic motifs. Structural clustering of ligand-bound conformations and binding free-energy calculations, using the AMBER force field and a continuum solvation model, can rapidly locate and screen numerous ligand-binding modes on the protein surface and detect the binding-site hot spot at the intermolecular interface. Structural orientation of the benzimidazole motif in the binding-site cavity closely mimics the position of the hot spot residue W104 in the crystal structure of the wild-type complex, which is recognized as an important structural requirement for restoring binding affinity. Despite numerous pockets on the protein surface of the mutant hGH-hGHbp complex, the binding-site cavity presents the energetically favorable hot spot for the benzimidazole-containing inhibitors, whereas for a set of nonactive molecules, the lowest energy ligand conformations do not necessarily bind in the engineered cavity. The results reveal a dominant role of the intermolecular van der Waals interactions in providing favorable ligand-protein energetics in the redesigned interface, in agreement with the experimental and computational alanine scanning of the hGH-hGHbp complex.

Importance of substrate and cofactor polarization in the active site of dihydrofolate reductase

By using a combined quantum-mechanical and molecular-mechanical potential in molecular dynamics simulations, we have investigated the effects of the enzyme electric field of dihydrofolate reductase on the electronic polarization of its 5-protonated dihydrofolate substrate at various stages of the catalyzed hydride transfer reaction. Energy decomposition of the total electrostatic interaction energy between the ligands and the enzyme shows that the polarization effect is 4% of the total electrostatic interaction energy, and, significantly, it accounts for 9kcal/mol of transition state stabilization relative to the reactant state. Therefore it is essential to take account of substrate polarization for quantitative interpretation of enzymatic function and for calculation of binding free energies of inhibitors to a protein. Atomic polarizations are calculated as the differences in the average atomic charges on the atoms in gas phase and in molecular simulations of the enzyme; this analysis shows that the glutamate tail and the pterin ring are the highly polarized regions of the substrate. Electron density difference plots of the reactant and product complexes at instantaneous configurations in the enzyme active center confirm the inferences made on the basis of partial atomic charges.

Monte Carlo simulations of the peptide recognition at the consensus binding site of the constant fragment of human immunoglobulin G: the energy landscape analysis of a hot spot at the intermolecular interface

Monte Carlo simulations of molecular recognition at the consensus binding site of the constant fragment (Fc) of human immunoglobulin G (Ig) protein have been performed to analyze structural and thermodynamic aspects of binding for the 13-residue cyclic peptide DCAWHLGELVWCT. The energy landscape analysis of a hot spot at the intermolecular interface using alanine scanning and equilibrium-simulated tempering dynamics with the simplified, knowledge-based energy function has enabled the role of the protein hot spot residues in providing the thermodynamic stability of the native structure to be determined. We have found that hydrophobic interactions between the peptide and the Met-252, Ile-253, His-433, and His-435 protein residues are critical to guarantee the thermodynamic stability of the crystallographic binding mode of the complex. Binding free energy calculations, using a molecular mechanics force field and a solvation energy model, combined with alanine scanning have been conducted to determine the energetic contribution of the protein hot spot residues in binding affinity. The conserved Asn-434, Ser-254, and Tyr-436 protein residues contribute significantly to the binding affinity of the peptide-protein complex, serving as an energetic hot spot at the intermolecular interface. The results suggest that evolutionary conserved hot spot protein residues at the intermolecular interface may be partitioned in fulfilling thermodynamic stability of the native binding mode and contributing to the binding affinity of the complex.

Complexity and simplicity of ligand-macromolecule interactions: the energy landscape perspective

The energy landscape approach has contributed to recent progress in understanding the complexity and simplicity of ligand-macromolecule interactions. Significant advances in computational structure prediction of ligand-protein complexes have been made using approaches that include the effects of protein flexibility and incorporate a hierarchy of energy functions. The results suggest that the complexity of structure prediction in molecular recognition may be determined by low-resolution properties of the underlying binding energy landscapes and by the nature of the energy funnels near the native structures of the complexes.

Hierarchy of simulation models in predicting structure and energetics of the Src SH2 domain binding to tyrosyl phosphopeptides

Structure and energetics of the Src Src Homology 2 (SH2) domain binding with the recognition phosphopeptide pYEEI and its mutants are studied by a hierarchical computational approach. The proposed structure prediction strategy includes equilibrium sampling of the peptide conformational space by simulated tempering dynamics with the simplified, knowledge-based energy function, followed by structural clustering of the resulting conformations and binding free energy evaluation of a single representative from each cluster, a cluster center. This protocol is robust in rapid screening of low-energy conformations and recovers the crystal structure of the pYEEI peptide. Thermodynamics of the peptide-SH2 domain binding is analyzed by computing the average energy contributions over conformations from the clusters, structurally similar to the predicted peptide bound structure. Using this approach, the binding thermodynamics for a panel of studied peptides is predicted in a better agreement with the experiment than previously suggested models. However, the overall correlation between computed and experimental binding affinity remains rather modest. The results of this study show that small differences in binding free energies between the Ala and Gly mutants of the pYEEI peptide are considerably more difficult to predict than the structure of the bound peptides, indicating that accurate computational prediction of binding affinities still remains a major methodological and technical challenge.

Hierarchy of simulation models in predicting molecular recognition mechanisms from the binding energy landscapes: structural analysis of the peptide complexes with SH2 domains

Computer simulations using the simplified energy function and simulated tempering dynamics have accurately determined the native structure of the pYVPML, SVLpYTAVQPNE, and SPGEpYVNIEF peptides in the complexes with SH2 domains. Structural and equilibrium aspects of the peptide binding with SH2 domains have been studied by generating temperature-dependent binding free energy landscapes. Once some native peptide-SH2 domain contacts are constrained, the underlying binding free energy profile has the funnel-like shape that leads to a rapid and consistent acquisition of the native structure. The dominant native topology of the peptide-SH2 domain complexes represents an extended peptide conformation with strong specific interactions in the phosphotyrosine pocket and hydrophobic interactions of the peptide residues C-terminal to the pTyr group. The topological features of the peptide-protein interface are primarily determined by the thermodynamically stable phosphotyrosyl group. A diversity of structurally different binding orientations has been observed for the amino-terminal residues to the phosphotyrosine. The dominant native topology for the peptide residues carboxy-terminal to the phosphotyrosine is tolerant to flexibility in this region of the peptide-SH2 domain interface observed in equilibrium simulations. The energy landscape analysis has revealed a broad, entropically favorable topology of the native binding mode for the bound peptides, which is robust to structural perturbations. This could provide an additional positive mechanism underlying tolerance of the SH2 domains to hydrophobic conservative substitutions in the peptide specificity region.

A genetic algorithm for structure-based de novo design

Genetic algorithms have properties which make them attractive in de novo drug design. Like other de novo design programs, genetic algorithms require a method to reduce the enormous search space of possible compounds. Most often this is done using information from known ligands. We have developed the ADAPT program, a genetic algorithm which uses molecular interactions evaluated with docking calculations as a fitness function to reduce the search space. ADAPT does not require information about known ligands. The program takes an initial set of compounds and iteratively builds new compounds based on the fitness scores of the previous set of compounds. We describe the particulars of the ADAPT algorithm and its application to three well-studied target systems. We also show that the strategies of enhanced local sampling and re-introducing diversity to the compound population during the design cycle provide better results than conventional genetic algorithm protocols.

Femtosecond studies of protein-ligand hydrophobic binding and dynamics: human serum albumin

In this contribution, we report studies of the nature of the dynamics and hydrophobic binding in protein-ligand complexes of human serum albumin with 2-(2'-hydroxyphenyl)-4-methyloxazole. With femtosecond time resolution, we examined the orientational motion of the ligand, its intrinsic nuclear motions, and the lifetime changes in the hydrophobic phase. For comparisons, with similar but chemical nanocavities, we also studied the same ligand in micelles and cyclodextrins. The hydrophobic interactions in the binding crevice are much stronger than those observed in cyclodextrins and micelles. The confined geometry restrains the nonradiative decay and significantly lengthens the excited-state lifetime. The observed dynamics over the femtosecond-to-nanosecond time scale indicate that the binding structure is rigid and the local motions of the ligand are nearly "frozen" in the protein. Another major finding is the elucidation of the directed dynamics by the protein. Proton transfer and intramolecular twisting of 2-(2'-hydroxyphenyl)-4-methyloxazole were observed to evolve along two routes: one involves the direct stretching motion in the molecular plane (approximately 200 fs) and is not sensitive to the environment; the second, less dominant, is related to the twisting motion (approximately 3 ps) of the two heterocyclic rings and drastically slows down in the protein hydrophobic pocket.