Theory of macromolecule-ligand interactions

A knowledge-guided strategy for improving the accuracy of scoring functions in binding affinity prediction

Background

Current scoring functions are not very successful in protein-ligand binding affinity prediction albeit their popularity in structure-based drug designs. Here, we propose a general knowledge-guided scoring (KGS) strategy to tackle this problem. Our KGS strategy computes the binding constant of a given protein-ligand complex based on the known binding constant of an appropriate reference complex. A good training set that includes a sufficient number of protein-ligand complexes with known binding data needs to be supplied for finding the reference complex. The reference complex is required to share a similar pattern of key protein-ligand interactions to that of the complex of interest. Thus, some uncertain factors in protein-ligand binding may cancel out, resulting in a more accurate prediction of absolute binding constants.

Results

In our study, an automatic algorithm was developed for summarizing key protein-ligand interactions as a pharmacophore model and identifying the reference complex with a maximal similarity to the query complex. Our KGS strategy was evaluated in combination with two scoring functions (X-Score and PLP) on three test sets, containing 112 HIV protease complexes, 44 carbonic anhydrase complexes, and 73 trypsin complexes, respectively. Our results obtained on crystal structures as well as computer-generated docking poses indicated that application of the KGS strategy produced more accurate predictions especially when X-Score or PLP alone did not perform well.

Conclusions

Compared to other targeted scoring functions, our KGS strategy does not require any re-parameterization or modification on current scoring methods, and its application is not tied to certain systems. The effectiveness of our KGS strategy is in theory proportional to the ever-increasing knowledge of experimental protein-ligand binding data. Our KGS strategy may serve as a more practical remedy for current scoring functions to improve their accuracy in binding affinity prediction.

Deciphering the ligand-binding sites in the Borrelia burgdorferi complement regulator-acquiring surface protein 2 required for interactions with the human immune regulators factor H and factor H-like protein 1

Borrelia burgdorferi, the etiologic agent of Lyme disease, employs sophisticated means to evade killing by its mammalian hosts. One important immune escape mechanism is the inhibition of complement activation mediated by interactions of the host-derived immune regulators factor H (CFH) and factor H-like protein 1 (CFHL1) with borrelial complement regulator-acquiring surface proteins (BbCRASPs). BbCRASP-2 is a distinctive CFH- and CFHL1-binding protein that is produced by serum-resistant B. burgdorferi strains. Here we show that binding of CFH by BbCRASP-2 is due to electrostatic as well as hydrophobic forces. In addition, 14 individual amino acid residues of BbCRASP-2 were identified as being involved in CFH and CFHL1 binding. Alanine substitutions of most of those residues significantly inhibited binding of CFH and/or CFHL1 by recombinant BbCRASP-2 proteins. To conclusively define the effects of BbCRASP-2 residue substitutions on serum sensitivity in the bacterial context, a serum-sensitive Borrelia garinii strain was transformed with plasmids that directed production of either wild-type or mutated BbCRASP-2 proteins. Critical amino acid residues within BbCRASP-2 were identified, with bacteria producing distinct mutant proteins being unable to bind either CFH or CFHL1, showing high levels of complement components C3, C6, and C5b-9 deposited on their surfaces and being highly sensitive to killing by normal serum. Collectively, we mapped a structurally sensitive CFH/CFHL1 binding site within borrelial BbCRASP-2 and identified single amino acid residues potentially involved in the interaction with both complement regulators.

Ligand-based assessment of factor Xa binding site flexibility via elaborate pharmacophore exploration and genetic algorithm-based QSAR modeling

The flexibility of activated factor X (fXa) binding site was assessed employing ligand-based pharmacophor modeling combined with genetic algorithm (GA)-based QSAR modeling. Four training subsets of wide structural diversity were selected from a total of 199 direct fXa inhibitors and were employed to generate different fXa pharmacophoric hypotheses using CATALYST software over two subsequent stages. In the first stage, high quality binding models (hypotheses) were identified. However, in the second stage, these models were refined by applying variable feature weight analysis to assess the relative significance of their features in the ligand-target affinity. The binding models were validated according to their coverage (capacity as a three-dimensional (3D) database search queries) and predictive potential as three-dimensional quantitative structure-activity relationship (3D-QSAR) models. Subsequently, GA and multiple linear regression (MLR) analysis were employed to construct different QSAR models from high quality pharmacophores and explore the statistical significance of combination models in explaining bioactivity variations across 199 fXa inhibitors. Three orthogonal pharmacophoric models emerged in the optimal QSAR equation suggesting they represent three binding modes accessible to ligands in the binding pocket within fXa.

Molecular dynamics simulations of the bacterial outer membrane protein FhuA: a comparative study of the ferrichrome-free and bound states

FhuA is one of the more complex members of the superfamily of bacterial outer membrane proteins. Its primary function is to provide a binding site on the outer membrane surface for siderophores, such as ferrichrome, and subsequently to facilitate their energy-dependent transport across the membrane, presumably powered by the TonB-ExbBD protein complex that resides in the cytoplasmic membrane. Crystal structures of FhuA with and without a bound ferrichrome molecule have provided some clues as to the initial stages of the siderophore transport mechanism. In the current study, we have employed 10-ns duration molecular dynamics simulations of FhuA and of the FhuA-ferrichrome complex, both embedded in a phospholipid bilayer, to probe the short timescale dynamics of this integral membrane protein, and to explore possible mechanistic implications of this dynamic behavior. Analysis of the dynamics of the protein suggests that the extracellular loops move as relatively rigid entities relative to the transmembrane beta-barrel. Comparison of the two simulations (with and without bound ferrichrome) revealed some ligand-induced changes in loop mobility. Specifically, loop L8 appears to be involved in a mechanism whereby the binding site is gated closed upon ligand binding. Analysis of the dynamics of water molecules within the core of the FhuA protein provided no evidence for a water-permeable protopore through which the ferrichrome might pass without a major perturbation of the FhuA protein. Overall, these simulations support the proposal that binding of ferrichrome initiates a signaling mechanism that ultimately leads to the TonB-mediated partial or total removal of the core domain from the beta-barrel, thus opening up a permeable pore. These simulations are among the longest that have been performed on a complex membrane protein. However, a simple analysis of sampling reveals that the description of protein motions is far from complete.

Solvation energy density occlusion approximation for evaluation of desolvation penalties in biomolecular interactions

In calculations involving many displacements of an interacting pair of biomolecules, such as brownian dynamics, the docking of a substrate/ligand to an enzyme/receptor, or the screening of a large number of ligands as prospective inhibitors for a particular receptor site, there is a need for rapid evaluation of the desolvation penalties of the interacting pair. Although continuum electrostatic treatments with distinct dielectric constants for solute and solvent provide an account of the electrostatics of solvation and desolvation, it is necessary to re-solve the Poisson equation, at considerable computational cost, for each displacement of the interacting pair. We present a new method that uses a formulation of continuum electrostatic solvation in terms of the solvation energy density and approximates desolvation in terms of the occlusion of this density. We call it the SEDO approximation. It avoids the need to re-solve the Poisson equation, as desolvation is now estimated by an integral over the occluded volume. Test calculations are presented for some simple model systems and for some real systems that have previously been studied using the Poisson equation approach: MHC class I protein-peptide complexes and a congeneric series of human immunodeficiency virus type 1 (HIV-1) protease--ligand complexes. For most of the systems considered, the trends and magnitudes of the desolvation component of interaction energies obtained using the SEDO approximation are in reasonable correlation with those obtained by re-solving the Poisson equation. In most cases, the error introduced by the SEDO approximation is much less than that of the often-used test-charge approximation for the charge-charge components of intermolecular interactions. Proteins 2001;43:12-27.

Energetics of drug-DNA interactions

Understanding the thermodynamics of drug binding to DNA is of both practical and fundamental interest. The practical interest lies in the contribution that thermodynamics can make to the rational design process for the development of new DNA targeted drugs. Thermodynamics offer key insights into the molecular forces that drive complex formation that cannot be obtained by structural or computational studies alone. The fundamental interest in these interactions lies in what they can reveal about the general problems of parsing and predicting ligand binding free energies. For these problems, drug-DNA interactions offer several distinct advantages, among them being that the structures of many drug-DNA complexes are known at high resolution and that such structures reveal that in many cases the drug acts as a rigid body, with little conformational change upon binding. Complete thermodynamic profiles (delta G, delta H, delta S, delta Cp) for numerous drug-DNA interactions have been obtained, with the help of high-sensitivity microcalorimetry. The purpose of this article is to offer a perspective on the interpretation of these thermodynamics parameters, and in particular how they might be correlated with known structural features. Obligatory conformational changes in the DNA to accommodate intercalators and the loss of translational and rotational freedom upon complex formation both present unfavorable free energy barriers for binding. Such barriers must be overcome by favorable free energy contributions from the hydrophobic transfer of ligand from solution into the binding site, polyelectrolyte contributions from coupled ion release, and molecular interactions (hydrogen and ionic bonds, van der Waals interactions) that form within the binding site. Theoretical and semiempirical tools that allow estimates of these contributions to be made will be discussed, and their use in dissecting experimental data illustrated. This process, even at the current level of approximation, can shed considerable light on the drug-DNA binding process.

Structure, modelling, and molecular dynamics studies of the inhibition of protein tyrosine phosphatase 1B by sulfotyrosine peptides

The protein tyrosine phosphatases comprise a class of enzymes that are crucial for the regulation of a number of cellular processes. Because of this, they are attracting increasing attention, not only as legitimate therapeutic targets, but also because of their relationship to many fundamental cellular processes. Certain sulfotyrosine peptides derived from casein are known to be good inhibitors of the protein tyrosine phosphatase, PTP1B. In this study, NMR transfer nuclear Overhauser effect studies have been used to ascertain the bound-state conformation adopted by the 12-amino acid residue casein-derived peptide, CAS200 (NANEEE(sY)SIGSA) and N-terminal truncated forms of this peptide, CAS203 and CAS205. Each of the peptides were found to bind in an extended beta-strand conformation. Extensive molecular modelling and molecular dynamics simulations of the PTP1B/peptide complexes, in a fully hydrated model, allowed a detailed description of the potential sources of the binding interactions to be developed. In agreement with the NMR studies, the modelling provided a picture of binding of CAS200 in which only the central (E203- I208) residues contributed significantly to the binding while the 3 N-terminal and 3 C-terminal residues were quite fluxional. Critical cationic surface residues, lying near to, but outside the active site pocket were the source of strong stabilizing forces that complemented the stabilizing interactions of the active site pocket. Electrostatic, hydrophobic, and hydrogen bonding interactions, in a residue specific manner, were all found to make significant contributions to the binding of these inhibitors.Key words: protein tyrosine phosphatase, PTP1B, casein peptide, inhibitor, NMR structure, molecular modelling, molecular dynamics.

Occluded molecular surface analysis of ligand-macromolecule contacts: application to HIV-1 protease-inhibitor complexes

Herein is described a method of quantifying and visualizing ligand-macromolecule contacts with the occluded surface algorithm by utilizing Connolly's van der Waals molecular surface dots together with the associated normals, to scoop out surrounding macromolecule atoms within a distance of 6.4 A from any ligand atom. On the basis of the intersections of surface normals with the van der Waals spheres of surrounding macromolecule atoms, the van der Waals molecular surface area for each atom is divided into occluded and nonoccluded surface areas. Also, we calculate a packing parameter for each occluded surface, measuring the closeness of the occluded surface against the macromolecule atom in contact. From the partial charges of ligand and macromolecule atoms and the occluded and nonoccluded surface areas due to the contact, we were able to identify favorable and unfavorable contacts. From the value of occluded surface constant, nonoccluded surface constant, and solvent-exposed constant for ligands, we qualitatively rank order the binding of ligands to the same target. From the individual parameters, group parameters for groups of atoms in a ligand or for each residue in a ligand-binding pocket of a macromolecule could be calculated. The group and the residue-based parameters could be used in structure-based ligand design and protein engineering experiments. In this paper, we present our analysis of ligand-macromolecule contacts, using five X-ray crystal structures of HIV-1 protease-ligand complexes.

How and why phosphotyrosine-containing peptides bind to the SH2 and PTB domains

BACKGROUND

Specific recognition of phosphotyrosine-containing protein segments by Src homology 2 (SH2) and phosphotyrosine-binding (PTB) domains plays an important role in intracellular signal transduction. Although many SH2/PTB-domain-containing receptor-peptide complex structures have been solved, little has been done to study the problem computationally. Prediction of the binding geometry and the binding constant of any peptide-protein pair is an extremely important problem.

RESULTS

A procedure to predict binding energies of phosphotyrosine-containing peptides with SH2/PTB domains was developed. The average deviation between experimentally measured binding energies and theoretical evaluations was 1.8 kcal/mol. Binding states of unphosphorylated peptides were also predicted reasonably well. Ab initio predictions of binding geometry of fully flexible peptides correctly identified conformations of two pentapeptides and a hexapeptide complexed with a v-Src SH2 domain receptor with root mean square deviations (rmsds) of 0.3 A, 1.2 A and 1.5 A, respectively.

CONCLUSIONS

The binding energies of phosphotyrosine-containing complexes can be effectively predicted using the procedure developed here. It was also possible to predict the bound conformations of flexible short peptides correctly from random starting conformations.

On the calculation of binding free energies using continuum methods: application to MHC class I protein-peptide interactions

This paper describes a methodology to calculate the binding free energy (delta G) of a protein-ligand complex using a continuum model of the solvent. A formal thermodynamic cycle is used to decompose the binding free energy into electrostatic and non-electrostatic contributions. In this cycle, the reactants are discharged in water, associated as purely nonpolar entities, and the final complex is then recharged. The total electrostatic free energies of the protein, the ligand, and the complex in water are calculated with the finite difference Poisson-Boltzmann (FDPB) method. The nonpolar (hydrophobic) binding free energy is calculated using a free energy-surface area relationship, with a single alkane/water surface tension coefficient (gamma aw). The loss in backbone and side-chain configurational entropy upon binding is estimated and added to the electrostatic and the nonpolar components of delta G. The methodology is applied to the binding of the murine MHC class I protein H-2Kb with three distinct peptides, and to the human MHC class I protein HLA-A2 in complex with five different peptides. Despite significant differences in the amino acid sequences of the different peptides, the experimental binding free energy differences (delta delta Gexp) are quite small (< 0.3 and < 2.7 kcal/mol for the H-2Kb and HLA-A2 complexes, respectively). For each protein, the calculations are successful in reproducing a fairly small range of values for delta delta Gcalc (< 4.4 and < 5.2 kcal/mol, respectively) although the relative peptide binding affinities of H-2Kb and HLA-A2 are not reproduced. For all protein-peptide complexes that were treated, it was found that electrostatic interactions oppose binding whereas nonpolar interactions drive complex formation. The two types of interactions appear to be correlated in that larger nonpolar contributions to binding are generally opposed by increased electrostatic contributions favoring dissociation. The factors that drive the binding of peptides to MHC proteins are discussed in light of our results.

Molecular docking to ensembles of protein structures

Until recently, applications of molecular docking assumed that the macromolecular receptor exists in a single, rigid conformation. However, structural studies involving different ligands bound to the same target biomolecule frequently reveal modest but significant conformational changes in the target. In this paper, two related methods for molecular docking are described that utilize information on conformational variability from ensembles of experimental receptor structures. One method combines the information into an "energy-weighted average" of the interaction energy between a ligand and each receptor structure. The other method performs the averaging on a structural level, producing a "geometry-weighted average" of the inter-molecular force field score used in DOCK 3.5. Both methods have been applied in docking small molecules to ensembles of crystal and solution structures, and we show that experimentally determined binding orientations and computed energies of known ligands can be reproduced accurately. The use of composite grids, when conformationally different protein structures are available, yields an improvement in computational speed for database searches in proportion to the number of structures.

On the bioactive conformation of a small peptide and its set of thermodynamically accessible conformations

In order to investigate the relationship between the bioactive conformation of a peptide and its set of thermodynamically accessible structures in solution, the conformational profile of the tetrapeptide Ac-Pro-Ala-Pro-Tyr-OH was characterized by computational methods. Search of the conformational space was performed within the molecular mechanics frame-work using the AMBER4.0 force field with an effective dielectric constant of 80. Unique structures of the peptide were compared with its bioactive conformation for the protein Streptomyces griseus Protease A, as taken from the crystal structure of the enzyme-peptide complex. The results show that the bound conformation is close to one of the unique conformations characterized in the conformational search of the isolated peptide. Moreover, the lowest energy minimum characterized in the conformational search exhibits large deviations when compared to the bound conformation of the crystal structure.

Low-resolution docking: prediction of complexes for underdetermined structures

One of the most fundamental questions concerning ligand-receptor interaction is whether such a process of intermolecular association is generally determined by local structural elements of the participating molecules, or whether there are also large-scale motifs in molecule structures that facilitate complex formation. From the point of view of practical docking computations, the elaborate character of local structural details in ligand-receptor interaction creates a large number of false-positive matches, which interfere with determination of the best fit. Another significant obstacle in protein docking is the problem of structural data inaccuracy (poor structure resolution, conformational changes upon complex formation, etc.). Our study [Vakser (1995) Protein Eng., 8, 371-377], based on ultralow (approximately 7 A resolution) representation of molecular structures, allowes to average all high-resolution structural details, and still predict most of the structural features of the ligand-receptor complex. The approach dramatically improves the signal-to-noise ratio in determination of the best fit, and moves the structure inaccuracy tolerance to the range of the macrostructure. In the present paper, we describe a further validation of the main principles of this approach and a detailed analysis of the low-resolution docking results. This includes clustering of ligand positions around the receptor molecule and cross-validation of ligands and receptors from different complexes. We also discuss the important implications of the approach to the multiple-minima problem and a possible role of different structural elements in the recognition mechanism.

Water molecules participate in proteinase-inhibitor interactions: crystal structures of Leu18, Ala18, and Gly18 variants of turkey ovomucoid inhibitor third domain complexed with Streptomyces griseus proteinase B

Crystal structures of the complexes of Streptomyces griseus proteinase B (SGPB) with three P1 variants of turkey ovomucoid inhibitor third domain (OMTKY3), Leu18, Ala18, and Gly18, have been determined and refined to high resolution. Comparisons among these structures and of each with native, uncomplexed SGPB reveal that each complex features a unique solvent structure in the S1 binding pocket. The number and relative positions of water molecules bound in the S1 binding pocket vary according to the size of the side chain of the P1 residue. Water molecules in the S1 binding pocket of SGPB are redistributed in response to the complex formation, probably to optimize hydrogen bonds between the enzyme and the inhibitor. There are extensive water-mediated hydrogen bonds in the interfaces of the complexes. In all complexes, Asn 36 of OMTKY3 participates in forming hydrogen bonds, via water molecules, with residues lining the S1 binding pocket of SGPB. For a homologous series of aliphatic straight side chains, Gly18, Ala18, Abu18, Ape18, and Ahp18 variants, the binding free energy is a linear function of the hydrophobic surface area buried in the interface of the corresponding complexes. The resulting constant of proportionality is 34.1 cal mol-1 A-2. These structures confirm that the binding of OMTKY3 to the preformed S1 pocket in SGPB involves no substantial structural disturbances that commonly occur in the site-directed mutagenesis studies of interior residues in other proteins, thus providing one of the most reliable assessments of the contribution of the hydrophobic effect to protein-complex stability.

Discriminating D1 and D2 agonists with a hydrophobic similarity index

Currently, methods for calculating molecular similarity indices have been developed for comparing steric, charge density, and molecular electrostatic potential (MEP) properties. Much of the existing technology may, however, be applied to the quantitative comparison of molecular hydrophobicities. In this article we present an empirical hydrophobic similarity index. We utilize atomic hydrophobic parameters derived from a quantum mechanical semiempirical wavefunction. Hydrophobicity at points on a grid is computed with a recently introduced "molecular lipophilicity potential." The overlap of pairs of molecules is calculated with the metric introduced by Carbó. This approach is applied to a case in which steric and electrostatic criteria have already been shown to be inadequate in rationalizing selectivity, namely, requirements for recognition at the dopamine D1 and D2 receptors. We demonstrate that, for a set of dopamine agonists, D1 ligands show higher similarity in this property that D2 analogs. This indicator of similarity is more successful at accounting for D1 selectivity than previous methods.

Structure-based drug design: progress, results and challenges

Protein structure-based drug design is rapidly gaining momentum. The new opportunities, developments and results in this field are almost unbelievable compared with the situation less than a decade ago.