Binding affinities for sulfonamide inhibitors with human thrombin using Monte Carlo simulations with a linear response method

The free energy landscape of small molecule unbinding

The spontaneous dissociation of six small ligands from the active site of FKBP (the FK506 binding protein) is investigated by explicit water molecular dynamics simulations and network analysis. The ligands have between four (dimethylsulphoxide) and eleven (5-diethylamino-2-pentanone) non-hydrogen atoms, and an affinity for FKBP ranging from 20 to 0.2 mM. The conformations of the FKBP/ligand complex saved along multiple trajectories (50 runs at 310 K for each ligand) are grouped according to a set of intermolecular distances into nodes of a network, and the direct transitions between them are the links. The network analysis reveals that the bound state consists of several subbasins, i.e., binding modes characterized by distinct intermolecular hydrogen bonds and hydrophobic contacts. The dissociation kinetics show a simple (i.e., single-exponential) time dependence because the unbinding barrier is much higher than the barriers between subbasins in the bound state. The unbinding transition state is made up of heterogeneous positions and orientations of the ligand in the FKBP active site, which correspond to multiple pathways of dissociation. For the six small ligands of FKBP, the weaker the binding affinity the closer to the bound state (along the intermolecular distance) are the transition state structures, which is a new manifestation of Hammond behavior. Experimental approaches to the study of fragment binding to proteins have limitations in temporal and spatial resolution. Our network analysis of the unbinding simulations of small inhibitors from an enzyme paints a clear picture of the free energy landscape (both thermodynamics and kinetics) of ligand unbinding.

Postprocessing of protein-ligand docking poses using linear response MM-PB/SA: application to Wee1 kinase inhibitors

Prediction of the binding strength of untested ligands is a central issue in structure-based drug design. In order to rapidly screen large compound databases, simple scoring schemes are often used in target-based virtual screening. The resulting scores often correlate poorly with biological affinities. More rigorous scoring methods, such as MM-PB/SA, correlate better with biological data by considering solvation effects and protein flexibility in the calculation of the binding free energy of a ligand. Here we describe the performance of a modified MM-PB/SA method on 222 Wee1 kinase inhibitors (48 pyridopyrimidine and 174 pyrrolocarbazole derivatives). Docking of these inhibitors into the available Wee1 kinase crystal structure yielded a consistent binding mode, and the derived MM-PB/SA models showed a significant correlation between calculated and experimental data (r(2) values between 0.64 and 0.67). Further study of these models on external test sets of Wee1 kinase inhibitors and structurally related decoys showed that a model based on a single kinase-inhibitor conformation can discriminate the active inhibitors from decoys. We also tested whether the linear interaction energy method with continuum electrostatics (LIECE) yields comparable results to MM-PB/SA and whether the LIECE and MM-PB/SA models can be applied for virtual screening of compound libraries.

Modeling Protein-Ligand Binding by Mining Minima

We present the first application of the mining minima algorithm to protein-small molecule binding. This end-point approach use an empirical force field and implicit solvent models, treats the protein binding-site as fully flexible and estimates free energies as sums over local energy wells. The calculations are found to yield encouraging agreement with experiment for three sets of HIV-1protease inhibitors and a set of phosphodiesterase 10a inhibitors. The contributions of various aspects of the model to its accuracy are examined, and the Poisson-Boltzmann correction is found to be the most critical. Interestingly, the computed changes in configurational entropy upon binding fall roughly along the same entropy-energy correlation previously observed for smaller host-guest systems. Strengths and weaknesses of the method are discussed, as are the prospects for enhancing accuracy and speed.

Computational design of protein-ligand binding: modifying the specificity of asparaginyl-tRNA synthetase

A method for computational design of protein-ligand interactions is implemented and tested on the asparaginyl- and aspartyl-tRNA synthetase enzymes (AsnRS, AspRS). The substrate specificity of these enzymes is crucial for the accurate translation of the genetic code. The method relies on a molecular mechanics energy function and a simple, continuum electrostatic, implicit solvent model. As test calculations, we first compute AspRS-substrate binding free energy changes due to nine point mutations, for which experimental data are available; we also perform large-scale redesign of the entire active site of each enzyme (40 amino acids) and compare to experimental sequences. We then apply the method to engineer an increased binding of aspartyl-adenylate (AspAMP) into AsnRS. Mutants are obtained using several directed evolution protocols, where four or five amino acid positions in the active site are randomized. Promising mutants are subjected to molecular dynamics simulations; Poisson-Boltzmann calculations provide an estimate of the corresponding, AspAMP, binding free energy changes, relative to the native AsnRS. Several of the mutants are predicted to have an inverted binding specificity, preferring to bind AspAMP rather than the natural substrate, AsnAMP. The computed binding affinities are significantly weaker than the native, AsnRS:AsnAMP affinity, and in most cases, the active site structure is significantly changed, compared to the native complex. This almost certainly precludes catalytic activity. One of the designed sequences has a higher affinity and more native-like structure and may represent a valid candidate for Asp activity.

Correlation analyses on binding affinity of substituted benzenesulfonamides with carbonic anhydrase using ab initio MO calculations on their complex structures

Quantitative structure-activity relationship analyses on the free energy change during complex formation between substituted benzenesulfonamides (BSAs) and bovine carbonic anhydrase II (bCA II) were performed using generilized Born/surface area (GB/SA) and ab initio fragment molecular orbital (FMO) calculations for the whole complex structures. The result shows that the overall free energy change is governed by the contribution from solvation and dissociation free energy changes accompanying by complex formation. The FMO-IFIE (interfragment interaction energy) analysis quantitatively revealed that the intrinsic interaction energy of bCA II with BSAs is mostly from interactions with amino acid residues in the active site of bCA II. The "Zn block" (Zn(2+) and three histidine residues coordinated to Zn(2+)) in the active site shows the lowest interaction energy and the greatest variance of interaction energy with BSAs through their coordination interaction. The proposed procedure was demonstrated to provide a quantitative basis for understanding a ligand-protein interaction at electronic and atomic levels.

Alchemical free energy simulations for biological complexes: powerful but temperamental

Free energy simulations compare multiple ligand:receptor complexes by "alchemically" transforming one into another, yielding binding free energy differences. Since their introduction in the 1980s, many technical and theoretical obstacles were surmounted, and the method ("MDFE," since molecular dynamics are often used) has matured into a powerful tool. We describe its current status, its effectiveness, and the challenges it faces. MDFE has provided chemical accuracy for many systems but remains expensive, with significant human overhead costs. The bottlenecks have shifted, partly due to increased computer power. To study diverse sets of ligands, force field availability and accuracy can be a major difficulty. Another difficulty is the frequent need to consider multiple states, related to sidechain protonation or buried waters, for example. Sophisticated, automated methods to sample these states are maturing, such as constant pH simulations. Meanwhile, combinations of MDFE and simpler approaches, like continuum dielectric models, can be very effective. As illustrations, we show how, with careful force field parameterization, MDFE accurately predicts binding specificities between complex tetracycline ligands and their targets. We describe substrate binding to the aspartyl-tRNA synthetase enzyme, where many distinct electrostatic states play a role, and a histidine and a Mg(2+) ion act as coupled switches that help enforce a strict preference for the aspartate substrate, relative to several analogs. Overall, MDFE has achieved a predictive status, where novel ligands can be studied and molecular recognition elucidated in depth. It should play an increasing role in the analysis of complex cellular processes and biomolecular engineering.

Molecular dynamics simulations show that conformational selection governs the binding preferences of imatinib for several tyrosine kinases

Tyrosine kinases transmit cellular signals through a complex mechanism, involving their phosphorylation and switching between inactive and active conformations. The cancer drug imatinib binds tightly to several homologous kinases, including Abl, but weakly to others, including Src. Imatinib specifically targets the inactive, so-called "DFG-out" conformation of Abl, which differs from the preferred, "DFG-in" conformation of Src in the orientation of a conserved Asp-Phe-Gly (DFG) activation loop. However, recent x-ray structures showed that Src can also adopt the DFG-out conformation and uses it to bind imatinib. The Src/Abl-binding free energy difference can thus be decomposed into two contributions. Contribution i measures the different protein-imatinib interactions when either kinase is in its DFG-out conformation. Contribution ii depends on the ability of imatinib to select or induce this conformation, i.e. on the relative stabilities of the DFG-out and DFG-in conformations of each kinase. Neither contribution has been measured experimentally. We use molecular dynamics simulations to show that contribution i is very small, 0.2 +/- 0.6 kcal/mol; imatinib interactions are very similar in the two kinases, including long range electrostatic interactions with the imatinib positive charge. Contribution ii, deduced using the experimental binding free energy difference, is much larger, 4.4 +/- 0.9 kcal/mol. Thus, conformational selection, easy in Abl, difficult in Src, underpins imatinib specificity. Contribution ii has a simple interpretation; it closely approximates the stability difference between the DFG-out and DFG-in conformations of apo-Src. Additional calculations show that conformational selection also governs the relative binding of imatinib to the kinases c-Kit and Lck. These results should help clarify the current framework for engineering kinase signaling.

Ligand and structure-based models for the prediction of ligand-receptor affinities and virtual screenings: Development and application to the beta(2)-adrenergic receptor

In this study, we evaluated the applicability of ligand-based and structure-based models to quantitative affinity predictions and virtual screenings for ligands of the beta(2)-adrenergic receptor, a G protein-coupled receptor (GPCR). We also devised and evaluated a number of consensus models obtained through partial least square regressions, to combine the strengths of the individual components. In all cases, the bioactive conformation of each ligand was derived from molecular docking at the crystal structure of the receptor. We identified the most effective models applicable to the different scenarios, in the presence or in the absence of a training set. For ranking the affinity of closely related analogs when a training set is available, a ligand-based consensus model (LI-CM) seems to be the best choice, while the structure-based MM-GBSA score seems the best alternative in the absence of a training set. For virtual screening purposes, the structure-based MM-GBSA score was found to be the method of choice. Consensus models consistently had performances superior or close to those of the best individual components, and were endowed with a significantly increased robustness. Given multiple models with no a priori knowledge of their predictive capabilities, constructing a consensus model ensures results very close to those that the best model alone would have yielded.

Design, asymmetric synthesis, and evaluation of pseudosymmetric sulfoximine inhibitors against HIV-1 protease

The HIV-1 protease is a validated drug target for the design of antiretroviral drugs to combat AIDS. We previously established the sulfoximine functionality as a valid transition state mimetic (TSM) in the HIV-1 protease inhibitors (PI) design and have identified a lead pseudosymmetric compound with nanomolar enzymatic inhibitory activity. Here, we report the asymmetric synthesis of this compound and its application in the synthesis of sulfoximine-based peptidomimetic HIV-1 protease inhibitors. Molecular modeling revealed the potential mode of binding of the sulfoximine inhibitor as a TSM. The predicted absolute binding free energies suggested similar inhibitory effect as observed in our enzymatic inhibitory studies.

Ligand and structure-based methodologies for the prediction of the activity of G protein-coupled receptor ligands

Accurate in silico models for the quantitative prediction of the activity of G protein-coupled receptor (GPCR) ligands would greatly facilitate the process of drug discovery and development. Several methodologies have been developed based on the properties of the ligands, the direct study of the receptor-ligand interactions, or a combination of both approaches. Ligand-based three-dimensional quantitative structure-activity relationships (3D-QSAR) techniques, not requiring knowledge of the receptor structure, have been historically the first to be applied to the prediction of the activity of GPCR ligands. They are generally endowed with robustness and good ranking ability; however they are highly dependent on training sets. Structure-based techniques generally do not provide the level of accuracy necessary to yield meaningful rankings when applied to GPCR homology models. However, they are essentially independent from training sets and have a sufficient level of accuracy to allow an effective discrimination between binders and nonbinders, thus qualifying as viable lead discovery tools. The combination of ligand and structure-based methodologies in the form of receptor-based 3D-QSAR and ligand and structure-based consensus models results in robust and accurate quantitative predictions. The contribution of the structure-based component to these combined approaches is expected to become more substantial and effective in the future, as more sophisticated scoring functions are developed and more detailed structural information on GPCRs is gathered.

The HSP90 binding mode of a radicicol-like E-oxime determined by docking, binding free energy estimations, and NMR 15N chemical shifts

We determine the binding mode of a macrocyclic radicicol-like oxime to yeast HSP90 by combining computer simulations and experimental measurements. We sample the macrocyclic scaffold of the unbound ligand by parallel tempering simulations and dock the most populated conformations to yeast HSP90. Docking poses are then evaluated by the use of binding free energy estimations with the linear interaction energy method. Comparison of QM/MM-calculated NMR chemical shifts with experimental shift data for a selective subset of backbone (15)N provides an additional evaluation criteria. As a final test we check the binding modes against available structure-activity-relationships. We find that the most likely binding mode of the oxime to yeast HSP90 is very similar to the known structure of the radicicol-HSP90 complex.

ATP binding to the C terminus of the Arabidopsis thaliana nitrate/proton antiporter, AtCLCa, regulates nitrate transport into plant vacuoles

Nitrate, one of the major nitrogen sources for plants, is stored in the vacuole. Nitrate accumulation within the vacuole is primarily mediated by the NO(3)(-)/H(+) exchanger AtCLCa, which belongs to the chloride channel (CLC) family. Crystallography analysis of hCLC5 suggested that the C-terminal domain, composed by two cystathionine beta-synthetase motifs in all eukaryotic members of the CLC family is able to interact with ATP. However, interaction of nucleotides with a functional CLC protein has not been unambiguously demonstrated. Here we show that ATP reversibly inhibits AtCLCa by interacting with the C-terminal domain. Applying the patch clamp technique to isolated Arabidopsis thaliana vacuoles, we demonstrate that ATP reduces AtCLCa activity with a maximum inhibition of 60%. ATP inhibition of nitrate influx into the vacuole at cytosolic physiological nitrate concentrations suggests that ATP modulation is physiologically relevant. ADP and AMP do not decrease the AtCLCa transport activity; nonetheless, AMP (but not ADP) competes with ATP, preventing inhibition. A molecular model of the C terminus of AtCLCa was built by homology to hCLC5 C terminus. The model predicted the effects of mutations of the ATP binding site on the interaction energy between ATP and AtCLCa that were further confirmed by functional expression of site-directed mutated AtCLCa.

QM/MM linear response method distinguishes ligand affinities for closely related metalloproteins

Design of selective ligands for closely related targets is becoming one of the most important tasks in the drug development. New tools, more precise than fast scoring functions and less demanding than sophisticated Free Energy Perturbation methods, are necessary to help accomplish this goal. The methods of intermediate complexity, characterizing individual contributions to the binding energy, have been an area of intense research in the past few years. Our recently developed quantum mechanical/molecular mechanical (QM/MM) modification of the Linear Response (LR) method describes the binding free energies as the sum of empirically weighted contributions of the QM/MM interaction energies and solvent-accessible surface areas for the time-averaged structures of hydrated complexes, obtained by molecular dynamics (MD) simulations. The method was applied to published data on 27 inhibitors of matrix metalloproteinase-3 (MMP-3). The two descriptors explained 90% of variance in the inhibition constants with RMSE of 0.245 log units. The QM/MM treatment is indispensable for characterization of the systems lacking suitable force-field expressions. In this case, it provided characteristics of H-bonds of the inhibitors to Glu202, charges of binding site atoms, and accurate coordination geometries of the ligands to catalytic zinc. The geometries were constrained during the MD simulations, which characterized conformational flexibility of the complexes and helped in the elucidation of the binding differences for related compounds. A comparison of the presented QM/MM LR results with those previously published for inhibition of MMP-9 by the same set of ligands showed that the QM/MM LR approach was able to distinguish subtle differences in binding affinities for MMP-3 and MMP-9, which did not exceed one order of magnitude. This precision level makes the approach a useful tool for design of selective ligands to similar targets, because the results can be safely extrapolated to maximize selectivity.

Binding free energy calculations of adenosine deaminase inhibitor and the effect of methyl substitution in inhibitors

The binding affinity of an inhibitor is often improved ten times or more by introducing a simple substituent, such as a methyl group or a chlorine atom. We have investigated this phenomenon in the case of adenosine deaminase (ADA) inhibitors using molecular dynamics (MD) simulations and binding free energy calculations, by the linear interaction energy (LIE) method. For MD simulations, the coordination bond parameters and partial charges of atoms around the zinc ion in ADA have been determined by referring to ab initio MO calculations. The calculated binding free energies for seven inhibitors agreed well with the experimental ones, with a maximum error of 1.2 kcal/mol. The effect of methyl substitution in inhibitor molecules was examined on the basis of MD trajectories. It is suggested that the increase in binding affinity is caused by both van der Waals stabilizations by amino acid residues in contact with the introduced methyl group and through favored overall interactions with surrounding residues in the binding pocket.

Coherent multidimensional vibrational spectroscopy of biomolecules: concepts, simulations, and challenges

The response of complex molecules to sequences of femtosecond infrared pulses provides a unique window into their structure, dynamics, and fluctuating environments. Herein we survey the basic principles of modern two-dimensional infrared (2DIR) spectroscopy, which analogous to those of multidimensional NMR spectroscopy. The perturbative approach for computing the nonlinear optical response of coupled localized chromophores is introduced and applied to the amide backbone transitions of proteins, liquid water, membrane lipids, and amyloid fibrils. The signals are analyzed using classical molecular dynamics simulations combined with an effective fluctuating Hamiltonian for coupled localized anharmonic vibrations whose dependence on the local electrostatic environment is parameterized by an ab initio map. Several simulation methods, (cumulant expansion of Gaussian fluctuation, quasiparticle scattering, the stochastic Liouville equations, direct numerical propagation) are surveyed. Chirality-induced techniques which dramatically enhance the resolution are demonstrated. Signatures of conformational and hydrogen-bonding fluctuations, protein folding, and chemical-exchange processes are discussed.

Synthesis and quantitative structure-activity relationship of fatty acid amide hydrolase inhibitors: modulation at the N-portion of biphenyl-3-yl alkylcarbamates

Alkylcarbamic acid biphenyl-3-yl esters are a class of fatty acid amide hydrolase (FAAH) inhibitors that comprises cyclohexylcarbamic acid 3'-carbamoylbiphenyl-3-yl ester (URB597), a compound with analgesic, anxiolytic-like and antidepressant-like properties in rat and mouse models. Here, we extended the structure-activity relationships (SARs) for this class of compounds by replacing the cyclohexyl ring of the parent compound cyclohexylcarbamic acid biphenyl-3-yl ester (URB524) (FAAH IC50 = 63 nM) with a selected set of substituents of different size, shape, flexibility, and lipophilicity. Docking experiments and linear interaction energy (LIE) calculations indicated that the N-terminal group of O-arylcarbamates fits within the lipophilic region of the substrate-binding site, mimicking the arachidonoyl chain of anandamide. Significant potency improvements were observed for the beta-naphthylmethyl derivative 4q (IC50 = 5.3 nM) and its 3'-carbamoylbiphenyl-3-yl ester 4z (URB880, IC50 = 0.63 nM), indicating that shape complementarity and hydrogen bonds are crucial to obtain highly potent inhibitors.

Towards the development of universal, fast and highly accurate docking/scoring methods: a long way to go

Accelerating the drug discovery process requires predictive computational protocols capable of reducing or simplifying the synthetic and/or combinatorial challenge. Docking-based virtual screening methods have been developed and successfully applied to a number of pharmaceutical targets. In this review, we first present the current status of docking and scoring methods, with exhaustive lists of these. We next discuss reported comparative studies, outlining criteria for their interpretation. In the final section, we describe some of the remaining developments that would potentially lead to a universally applicable docking/scoring method.

QSAR: dead or alive?

This perspective concerns the methods employed within the current drug discovery community to develop predictive quantitative structure-activity relationships (QSAR). Specifically, a number of cautions are provided which may circumvent misuse and misunderstanding of the technique. Ignorance of such caveats has led to a discouraging tendency of the methods to result in poorly predictive models. Among these pitfalls are the fondness with which we associate correlation with causation, the mesmerizing influence of large numbers of molecular descriptors, the incessant misuse of the leave-one-out paradigm, and finally, the QSAR enigma wherein model predictivity is not a necessary component of a model's usefulness.

Calculation of binding affinities of HIV-1 RT and beta-secretase inhibitors using the linear interaction energy method with explicit and continuum solvation approaches

The linear interaction energy (LIE) approach has been applied to estimate the binding free energies of representative sets of HIV-1 RT and beta-Secretase inhibitors, using both molecular dynamics (MD) and tethered energy minimization sampling protocols with the OPLS-AA potential, using a range of solvation methodologies. Generalized Born (GB), 'shell' and periodic boundary condition (PBC) solvation were used, the latter with reaction field (RF) electrostatics. Poisson-Boltzmann (PB) and GB continuum electrostatics schemes were applied to the simulation trajectories for each solvation type to estimate the electrostatic ligand-water interaction energy in both the free and bound states. Reasonable agreement of the LIE predictions was obtained with respect to experimental binding free energy estimates for both systems: for instance, 'PB' fits on MD trajectories carried out with PBC solvation and RF electrostatics led to models with standard errors of 1.11 and 1.03 kcal mol(-1) and coefficients of determination, r (2) of 0.76 and 0.75 for the HIV-1 RT and beta-Secretase sets. However, it was also found that results from MD sampling using PBC solvation provided only slightly better fits than from simulations using shell or Born solvation or tethered energy minimization sampling.

Quantum and molecular dynamics study for binding of macrocyclic inhibitors to human alpha-thrombin

Molecular dynamics simulations followed by quantum mechanical calculation and Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) analysis have been carried out to study binding of proline- and pyrazinone-based macrocyclic inhibitors (L86 and T76) to human alpha-thrombin. Detailed binding interaction energies between these inhibitors and individual protein fragments are calculated using DFT method based on a new quantum mechanical approach for computing protein-ligand interaction energy. The analysis of detailed interaction energies provides insight on the protein-ligand binding mechanism. Study shows that T76 and L86 bind to thrombin in a very similar "inhibition mode" except that T76 has relatively weaker binding interaction with Glu(217). The analysis from quantum calculation of binding interaction is consistent with the MM-PBSA calculation of binding free energy, and the calculated free energies for L86/T76-thrombin binding agree well with the experimental data.

Calculation of Protein-Ligand Binding Affinities

Accurate methods of computing the affinity of a small molecule with a protein are needed to speed the discovery of new medications and biological probes. This paper reviews physics-based models of binding, beginning with a summary of the changes in potential energy, solvation energy, and configurational entropy that influence affinity, and a theoretical overview to frame the discussion of specific computational approaches. Important advances are reported in modeling protein-ligand energetics, such as the incorporation of electronic polarization and the use of quantum mechanical methods. Recent calculations suggest that changes in configurational entropy strongly oppose binding and must be included if accurate affinities are to be obtained. The linear interaction energy (LIE) and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) methods are analyzed, as are free energy pathway methods, which show promise and may be ready for more extensive testing. Ultimately, major improvements in modeling accuracy will likely require advances on multiple fronts, as well as continued validation against experiment.

Binding diversity of the two binding sites of ricin B lectin

Ricin B is a galactose-binding protein, which contains two binding sites. We have compared the binding properties of the two binding sites of ricin B chain toward different mono- and disaccharide ligands. The free energies of binding are calculated using the free energy perturbation simulation (thermodynamic integration method) and linear interaction energy approach using CHARMM force field. The second binding site of the protein was found to be weaker compared to the first. The details of the hydrogen-bonding scheme suggested the origin of the epimeric specificity of the protein. The reason for the weaker binding capacity of the second binding site has been addressed.

Linear interaction energy models for β-secretase (BACE) inhibitors: Role of van der Waals, electrostatic, and continuum-solvation terms

Computing the binding affinity of a protein-ligand complex is one of the most fundamental and difficult tasks in computer-aided drug design. Many approaches for computing binding affinities can be classified as linear interaction energy (LIE) models as they rely on some type of linear fit of computed interaction energies between ligand and protein. We have examined the computed interaction energies of a series of beta-secretase (BACE) inhibitors in terms of van der Waals, coulombic, and continuum-solvation contributions to ligand binding. We have also systematically examined the effect of different protonation states of the protein and ligands. We find that the binding affinities are relatively insensitive to the protonation state of the protein when neutral ligands are considered. Inclusion of charged ligands leads to large deviations in the coulomb, solvation, and even van der Waals terms. The latter is due to increased repulsive van der Waals interactions in the complex due to the strong coulomb attraction found between oppositely charged functional groups in the protein and ligand. In general, we find that the best models are obtained when the protein is judiciously charged (e.g. Asp32-, Arg235+) and the potentially charged ligands are treated as neutral.

Combining docking and molecular dynamic simulations in drug design

A rational approach is needed to maximize the chances of finding new drugs, and to exploit the opportunities of potential new drug targets emerging from genomic and proteomic initiatives, and from the large libraries of small compounds now readily available through combinatorial chemistry. Despite a shaky early history, computer-aided drug design techniques can now be effective in reducing costs and speeding up drug discovery. This happy outcome results from development of more accurate and reliable algorithms, use of more thoughtfully planned strategies to apply them, and greatly increased computer power to allow studies with the necessary reliability to be performed. Our review focuses on applications and protocols, with the main emphasis on critical analysis of recent studies where docking calculations and molecular dynamics (MD) simulations were combined to dock small molecules into protein receptors. We highlight successes to demonstrate what is possible now, but also point out drawbacks and future directions. The review is structured to lead the reader from the simpler to more compute-intensive methods. Thus, while inexpensive and fast docking algorithms can be used to scan large compound libraries and reduce their size, more accurate but expensive MD simulations can be applied when a few selected ligand candidates remain. MD simulations can be used: during the preparation of the protein receptor before docking, to optimize its structure and account for protein flexibility; for the refinement of docked complexes, to include solvent effects and account for induced fit; to calculate binding free energies, to provide an accurate ranking of the potential ligands; and in the latest developments, during the docking process itself to find the binding site and correctly dock the ligand a priori.

Molecular modeling of organic and biomolecular systems using BOSS and MCPRO

An overview is provided of the capabilities for the current versions of the BOSS and MCPRO programs for molecular modeling of organic and biomolecular systems. Recent applications are noted, particularly for QM/MM studies of organic and enzymatic reactions and for protein-ligand binding.

A combination of docking, QM/MM methods, and MD simulation for binding affinity estimation of metalloprotein ligands

To alleviate the problems in the receptor-based design of metalloprotein ligands due to inadequacies in the force-field description of coordination bonds, a four-tier approach was devised. Representative ligand-metalloprotein interaction energies are obtained by subsequent application of (1) docking with metal-binding-guided selection of modes, (2) optimization of the ligand-metalloprotein complex geometry by combined quantum mechanics and molecular mechanics (QM/MM) methods, (3) conformational sampling of the complex with constrained metal bonds by force-field-based molecular dynamics (MD), and (4) a single point QM/MM energy calculation for the time-averaged structures. The QM/MM interaction energies are, in a linear combination with the desolvation-characterizing changes in the solvent-accessible surface areas, correlated with experimental data. The approach was applied to structural correlation of published binding free energies of a diverse set of 28 hydroxamate inhibitors to zinc-dependent matrix metalloproteinase 9 (MMP-9). Inclusion of steps 3 and 4 significantly improved both correlation and prediction. The two descriptors explained 90% of variance in inhibition constants of all 28 inhibitors, ranging from 0.08 to 349 nM, with the average unassigned error of 0.318 log units. The structural and energetic information obtained from the time-averaged MD simulation results helped understand the differences in binding modes of related compounds.

An Experimental and Molecular-Modeling Study of the Binding of Linked Sulfated Tetracyclitols to FGF-1 and FGF-2

The experimental binding affinities of a series of linked sulfated tetracyclitols [Cyc2N-R-NCyc2, where Cyc = C6H6(OSO3Na)3 and R = (CH2)n (n = 2-10), p-xylyl or (C2H4)2-Ncyc] for the fibroblast growth factors FGF-1 and FGF-2 have been measured by using a surface plasmon resonance assay. The KD values range from 7.0 nM to 1.1 microM for the alkyl-linked ligands. The binding affinity is independent of the flexibility of the linker, as replacement of the alkyl linker with a rigid p-xylyl group did not affect the KD. Calculations suggest that binding modes for the p-xylyl-linked ligand are similar to those calculated for the flexible alkyl-linked tetracyclitols. The possible formation of cross-linked FGF:cyclitol complexes was examined by determining KD values at increasing protein concentrations. No changes in KD were observed; this suggesting that only 1:1 complexes are formed under these assay conditions. Monte Carlo multiple-minima calculations of low-energy conformers of the FGF-bound ligands showed that all of the sulfated tetracyclitol ligands can bind effectively in the heparan sulfate-binding sites of FGF-1 and FGF-2. Binding affinities of these complexes were estimated by the Linear Interaction Energy (LIE) method to within a root-mean-square deviation of 1 kcal mol(-1) of the observed values. The effect of incorporating cations to balance the overall charge of the complexes during the LIE calculations was also explored.

Processing multimode binding situations in simulation-based prediction of ligand-macromolecule affinities

The linear response (LR) approximation and similar approaches belong to practical methods for estimation of ligand-receptor binding affinities. The approaches correlate experimental binding affinities with the changes upon binding of the ligand electrostatic and van der Waals energies and of solvation characteristics. These attributes are expressed as ensemble averages that are obtained by conformational sampling of the protein-ligand complex and of the free ligand by molecular dynamics or Monte Carlo simulations. We observed that outliers in the LR correlations occasionally exhibit major conformational changes of the complex during sampling. We treated the situation as a multimode binding case, for which the observed association constant is the sum of the partial association constants of individual states/modes. The resulting nonlinear expression for the binding affinities contains all the LR variables for individual modes that are scaled by the same two to four adjustable parameters as in the one-mode LR equation. The multimode method was applied to inhibitors of a matrix metalloproteinase, where this treatment improved the explained variance in experimental activity from 75% for the unimode case to about 85%. The predictive ability scaled accordingly, as verified by extensive cross-validations.

Design, synthesis and evaluation of graftable thrombin inhibitors for the preparation of blood-compatible polymer materials

Piperazinyl-amide derivatives of N-alpha-(3-trifluoromethyl-benzenesulfonyl)-L-arginine (1) were synthesized as graftable thrombin inhibitors. The possible disturbance of biological activity due to a variable spacer-arm fixed on the N-4 piperazinyl position was evaluated in vitro, against human alpha-thrombin, and in blood coagulation assay. Molecular modelling (in silico analysis) and X-ray diffraction studies of thrombin-inhibitor complexes were also performed. The fixation of bioactive molecules on poly(butylene terephthalate) (PBT) and poly(ethylene terephthalate) (PET) membranes was performed by wet chemistry treatment and evaluated by XPS analysis. Surface grafting of inhibitor 1d improved the membrane hemocompatibility by reducing blood clot formation on the modified surface.

HIV-1 reverse transcriptase variants: molecular modeling of Y181C, V106A, L100I, and K103N mutations with nonnucleoside inhibitors using Monte Carlo simulations in combination with a linear response method

The energies and physical descriptors for the binding of 21 novel 1-(2,6-difluorobenzyl)-2-(2,6-difluorophenyl)-benzimidazole (BPBI) analogs to HIV-1 reverse transcriptase (RT) variants Y181C, L100I, V106A, and K103N have been determined using Monte Carlo (MC) simulations. The crystallographic structure of the lead compound, 4-methyl BPBI, was used as a starting point to model the inhibitors in both the mutant bound and the unbound states. The energy terms and physical descriptors obtained from the calculations were reasonably correlated with the respective experimental EC50 values for the inhibitors against the various mutant RTs. Using the linear response correlations from the calculations, 2 novel BPBI inhibitors have been designed and simulations have been carried out. The results show the computed deltaG(binding) values match the experimental data for the analogs. Given the ongoing problem with drug resistance, the ability to predict the activity of novel analogs against variants prior to synthesis is highly advantageous.

Protein-ligand binding free energy estimation using molecular mechanics and continuum electrostatics. Application to HIV-1 protease inhibitors

A method is proposed for the estimation of absolute binding free energy of interaction between proteins and ligands. Conformational sampling of the protein-ligand complex is performed by molecular dynamics (MD) in vacuo and the solvent effect is calculated a posteriori by solving the Poisson or the Poisson-Boltzmann equation for selected frames of the trajectory. The binding free energy is written as a linear combination of the buried surface upon complexation, SASbur, the electrostatic interaction energy between the ligand and the protein, Eelec, and the difference of the solvation free energies of the complex and the isolated ligand and protein, deltaGsolv. The method uses the buried surface upon complexation to account for the non-polar contribution to the binding free energy because it is less sensitive to the details of the structure than the van der Waals interaction energy. The parameters of the method are developed for a training set of 16 HIV-1 protease-inhibitor complexes of known 3D structure. A correlation coefficient of 0.91 was obtained with an unsigned mean error of 0.8 kcal/mol. When applied to a set of 25 HIV-1 protease-inhibitor complexes of unknown 3D structures, the method provides a satisfactory correlation between the calculated binding free energy and the experimental pIC5o without reparametrization.

Binding affinity prediction with different force fields: examination of the linear interaction energy method

A systematic study of the linear interaction energy (LIE) method and the possible dependence of its parameterization on the force field and system (receptor binding site) is reported. We have calculated the binding free energy for nine different ligands in complex with P450cam using three different force fields (Amber95, Gromos87, and OPLS-AA). The results from these LIE calculations using our earlier parameterization give relative free energies of binding that agree remarkably well with the experimental data. However, the absolute energies are too positive for all three force fields, and it is clear that an additional constant term (gamma) is required in this case. Out of five examined LIE models, the same one emerges as the best for all three force fields, and this, in fact, corresponds to our earlier one apart from the addition of the constant gamma, which is almost identical for the three force fields. Thus, the present free energy calculations clearly indicate that the coefficients of the LIE method are independent of the force field used. Their relation to solvation free energies is also demonstrated. The only free parameter of the best model is gamma, which is found to depend on the hydrophobicity of the binding site. We also attempt to quantify the binding site hydrophobicity of four different proteins which shows that the ordering of gamma's for these sites reflects the fraction of hydrophobic surface area.

Application of the linear interaction energy method (LIE) to estimate the binding free energy values of Escherichia coli wild-type and mutant arginine repressor C-terminal domain (ArgRc)–l-arginine and ArgRc–l-citrulline protein–ligand complexes

Protein-ligand binding free energy values of wild-type and mutant C-terminal domain of Escherichia coli arginine repressor (ArgRc) protein systems bound to L-arginine or L-citrulline molecules were calculated using the linear interaction energy (LIE) method by molecular dynamics (MD) simulation. The binding behaviour predicted by the dissociation constant (K(d)) calculations from the binding free energy values showed preferences for binding of L-arginine to the wild-type ArgRc but not to the mutant ArgRc(D128N). On the other hand, L-citrulline do not favour binding to wild-type ArgRc but prefer binding to mutant ArgRc(D128N). The dissociation constant for the wild-type ArgRc-L-arginine complex obtained in this study is in agreement with reported experimental results. Our results also support the experimental data for the binding of L-citrulline to the mutant ArgRc(D128N). These showed that LIE method for protein-ligand binding free energy calculation could be applied to the wild-type and the mutant E. coli ArgRc-L-arginine and ArgRc-L-citrulline protein-ligand complexes and possibly to other transcriptional repressor-co-repressor systems as well.