Search

Research has been conducted on the haloalkane dehalogenase substrates from Xanthobacter autotrophicus GJ10. The results of the comparative binding energy analysis carried out for these substrates are described.

A study was conducted to analyze the computation of absolute protein-ligand binding free energy using an alchemistic application of the free energy perturbation technique. Camphor was exchanged with water molecules in molecular dynamics simulations in a protein region characterized by an active site and an associated free energy change. Results indicated support for the validity of the modified technique and the computation procedure.

The stability of the salt-link bridges in cytochrome P450cam was investigated through detailed numerical computation of the finite-difference linearized Poisson-Boltzmann equation. Results have shown that mutations to the hydrophobic isosteres determine the stability and instability of salt links. Favorable interactions between salt link center and protein moiety also play an additive role in the stability of the salt-link bonds.

The effects of substrate-induced changes in cytochrome P-450cam on protein conformation were investigated based on four criteria: spin-state equilibrium, substrate dissociation, thermal unfolding, and subconformer equilibria. The results showed that binding of the enantiomers (1R)-camphor and (1S)-camphor to cytochrome P-450cam enhanced ligand binding to the active site of the hemothiolate proteins. Data indicated that a high-spin state was more favored in (1R)-camphor-bound P-450 than in (1S)-camphor complex.

Four novel camphor analogues of the cytochrome P450cam were generated by replacing the camphor carbonyl group by an ether oxygen attached to a hydrocarbon tail of different lengths. Findings showed that the binding energies determined from molecular dynamics simulations of the protein-ligand interactions, protein-water interactions and ligand-water interactions closely agree with empirical binding enthalpies.

A detailed model of electron transfer (ET) events that was applied to all of the generated protein configurations was used to compute ET rates and understand Brownian dynamics simulations of the diffusional association of plastocyanin (PC) and cytochrome f (CytF) from the cyanobacterium Phormidium laminosum. The simulations offer insight into the strong electrostatic interactions between P. laminosum PC and CytF which influenced the ET rates from both diffusion and thermal activation processes.

A comparison of the binding and reactivity of plant and mammalian peroxidases to indole derivatives by computational docking is presented. The differentially lowered reaction rates of co II of horseradish peroxidase (HRP) and myeloperoxidase (MPO) with respect to those of co I might stem from as yet undetermined conformational or electrostatic differences between the co I and co II states of MPO, which are absent in HRP.

The classical quantitative structure-activity relationship (QSAR) technique and comparative binding energy (COMBINE) analysis is simultaneously applied to analyze substrate specificity in the haloalkene dehalogenase LinB enzyme from Sphingomonas paucimobilis UT26. The mechanism of molecular adaptation toward a different specificity is discussed on the basis of quantitative comparison of models derived for two protein family members.

Feed-forwardneural networks trained by crystallographic data to predict water binding sitesin proteins on the basis of amino acid sequence and secondary structure were described. The networks were used to determine in a test protein the residues bound to a water ligand and atoms in close contact with a water molecule. The performance of the networks can be expected to improve as more accurate proteinstructures become available.

An evaluation of the search capabilities of three widely used programs, namely, GRID, DOCK and AutoDock, for heparin binding sites was conducted using known crystal structures of the protein-heparin complexes antithrombin, and acidic and basic fibroblast growth factor. The heparin binding site for all test cases was correctly determined with all three programs. When the protocols were used to predict the heparin binding site on Interleukin-8 (IL-8), results suggested that His18, Lys20, Arg60, Lys64 and Arg68 in IL-8 all bind to heparin.

A study was conducted to analyze the YTGP tetrapeptide in aqueous solution using multiple molecular dynamics trajectories. Random coil values were utilized to determine the chemical shifts for the protons in an unstructured peptide. Trajectories were generated by propagation for 300 ps at 300 K an 1 atm. Experimental results indicated that the structural shift spectra for different conformation states correlated with the experimental spectrum.

High macromolecular concentrations are a distinguishingfeatureof living organisms. Understanding how the high concentration of solutesaffects the dynamic properties of biological macromolecules is fundamentalfor the comprehension of biological processes in living systems. Inthis paper, we describe the implementation of mean field models oftranslational and rotational hydrodynamic interactions into an atomicallydetailed many-protein Brownian dynamics simulation method. Concentratedsolutions (30–40% volume fraction) of myoglobin, hemoglobinA, and sickle cell hemoglobin S were simulated, and static structurefactors, oligomer formation, and translational and rotational self-diffusioncoefficients were computed. Good agreement of computed propertieswith available experimental data was obtained. The results show theimportance of both solvent mediated interactions and weak protein–proteininteractions for accurately describing the dynamics and the associationproperties of concentrated protein solutions. Specifically, they showa qualitative difference in the translational and rotational dynamicsof the systems studied. Although the translational diffusion coefficientis controlled by macromolecular shape and hydrodynamic interactions,the rotational diffusion coefficient is affected by macromolecularshape, direct intermolecular interactions, and both translationaland rotational hydrodynamic interactions. [ABSTRACT FROM AUTHOR]

The factors that determine the extent to which diffusion and thermal activation processes govern electron transfer (ET) between proteins are debated. The process of ET between plastocyanin (PC) and cytochrome f(CytF) from the cyanobacterium Phormidium laminosum was initially thought to be diffusion-controlled but later was found to be under activation control (Schlarb-Ridley, B. G.; et al. Biochemistry 2005, 44, 6232). Here we describe Brownian dynamics simulations of the diffusional association of PC and CyIF, from which El rates were computed using a detailed model of ET events that was applied to all of the generated protein configurations. The proteins were modeled as rigid bodies represented in atomic detail. In addition to electrostatic forces, which were modeled as in our previous simulations of protein-protein association, the proteins interacted by a nonpolar desolvation (hydrophobic) force whose derivation is described here. The simulations yielded close to realistic residence times of transient protein-protein encounter complexes of up to tens of microseconds. The activation barrier for individual ET events derived from the simulations was positive. Whereas the electrostatic interactions between P. laminosum PC and CytF are weak, simulations for a second cyanobacterial PC-CytF pair, that from Nostoc sp. PCC 7119, revealed El rates influenced by stronger electrostatic interactions. In both cases, the simulations imply significant contributions to ET from both diffusion and thermal activation processes. [ABSTRACT FROM AUTHOR]

There is a growing appreciation of the importance of drug-target binding kinetics for lead optimization. For G protein-coupled receptors (GPCRs), which mediate signaling over a wide range of time scales, the drug dissociation rate is often a better predictor of in vivo efficacy than binding affinity, although it is more challenging to compute. Here, we assess the ability of the τ-Random Acceleration Molecular Dynamics (τRAMD) approach to reproduce relative residence times and reveal dissociation mechanisms and the effects of allosteric modulation for two important membrane-embedded drug targets: the β2-adrenergic receptor and the muscarinic acetylcholine receptor M2. The dissociation mechanisms observed in the relatively short RAMD simulations (in which molecular dynamics (MD) simulations are performed using an additional force with an adaptively assigned random orientation applied to the ligand) are in general agreement with much more computationally intensive conventional MD and metadynamics simulations. Remarkably, although decreasing the magnitude of the random force generally reduces the number of egress routes observed, the ranking of ligands by dissociation rate is hardly affected and agrees well with experiment. The simulations also reproduce changes in residence time due to allosteric modulation and reveal associated changes in ligand dissociation pathways.

The binding kinetic properties of potential drugs may significantly influence their subsequent clinical efficacy. Predictions of these properties based on computer simulations provide a useful alternative to their expensive and time-consuming experimental counterparts, even at an early drug discovery stage. Herein, we perform scaled molecular dynamics (ScaledMD) simulations on a set of 27 ligands of HSP90 belonging to more than seven chemical series to estimate their relative residence times. We introduce two new techniques for the analysis and the classification of the simulated unbinding trajectories. The first technique, which helps in estimating the limits of the free energy well around the bound state, and the second one, based on a new contact map fingerprint, allow the description and the comparison of the paths that lead to unbinding. Using these analyses, we find that ScaledMD's relative residence time generally enables the identification of the slowest unbinders. We propose an explanation for the underestimation of the residence times of a subset of compounds, and we investigate how the biasing in ScaledMD can affect the mechanistic insights that can be gained from the simulations.

There is growing consensus that the optimization of the kinetic parameters for drug-protein binding leads to improved drug efficacy. Therefore, computational methods have been developed to predict kinetic rates and to derive quantitative structure-kinetic relationships (QSKRs). Many of these methods are based on crystal structures of ligand-protein complexes. However, a drawback is that each ligand-protein complex is usually treated as having a single structure. Here, we present a modification of COMparative BINding Energy (COMBINE) analysis, which uses the structures of ligand-protein complexes to predict binding parameters. We introduce the option of using multiple structures to describe each ligand-protein complex in COMBINE analysis and apply this to study the effects of protein flexibility on the derivation of dissociation rate constants ( k off ) for inhibitors of p38 mitogen-activated protein (MAP) kinase, which has a flexible binding site. Multiple structures were obtained for each ligand-protein complex by performing docking to an ensemble of protein configurations obtained from molecular dynamics simulations. Coefficients to scale ligand-protein interaction energies determined from energy-minimized structures of ligand-protein complexes were obtained by partial least squares regression, and they allowed for the computation of k off values. The QSKR model obtained using single, energy-minimized crystal structures for each ligand-protein complex had higher predictive power than the QSKR model obtained with multiple structures from ensemble docking. However, incorporation of ligand-protein flexibility helped to highlight additional ligand-protein interactions that lead to longer residence times, such as interactions with residues Arg67 and Asp168, which are close to the ligand in many crystal structures. These results show that COMBINE analysis is a promising method to guide the design of compounds that bind to flexible proteins with improved binding kinetics.

Accurate protein druggability predictions are important for the selection of drug targets in the early stages of drug discovery. Because of the flexible nature of proteins, the druggability of a binding pocket may vary due to conformational changes. We have therefore developed two statistical models, a logistic regression model (TRAPP-LR) and a convolutional neural network model (TRAPP-CNN), for predicting druggability and how it varies with changes in the spatial and physicochemical properties of a binding pocket. These models are integrated into TRAnsient Pockets in Proteins (TRAPP), a tool for the analysis of binding pocket variations along a protein motion trajectory. The models, which were trained on publicly available and self-augmented datasets, show equivalent or superior performance to existing methods on test sets of protein crystal structures and have sufficient sensitivity to identify potentially druggable protein conformations in trajectories from molecular dynamics simulations. Visualization of the evidence for the decisions of the models in TRAPP facilitates identification of the factors affecting the druggability of protein binding pockets.

A growing consensus is emerging that optimizing the drug-target affinity alone under equilibrium conditions does not necessarily translate into higher potency in vivo and that instead binding kinetic parameters should be optimized to ensure better efficacy. Therefore, in silico methods are needed to predict the kinetic parameters and the mechanistic determinants of drug-protein binding. Here we demonstrate the application of COMparative BINding Energy (COMBINE) analysis to derive quantitative structure-kinetics relationships (QSKRs) for the dissociation rate constants ( k off ) of inhibitors of heat shock protein 90 (HSP90) and HIV-1 protease. We derived protein-specific scoring functions by correlating k off rate constants with a subset of weighted interaction energy components determined from the energy-minimized structures of drug-protein complexes. As the QSKRs derived for these sets of chemically diverse compounds have good predictive ability and provide insights into important drug-protein interactions for optimizing k off , COMBINE analysis offers a promising approach for binding kinetics-guided lead optimization., Competing Interests: The authors declare no competing financial interest.

Drug-target residence time (τ), one of the main determinants of drug efficacy, remains highly challenging to predict computationally and, therefore, is usually not considered in the early stages of drug design. Here, we present an efficient computational method, τ-random acceleration molecular dynamics (τRAMD), for the ranking of drug candidates by their residence time and obtaining insights into ligand-target dissociation mechanisms. We assessed τRAMD on a data set of 70 diverse drug-like ligands of the N-terminal domain of HSP90α, a pharmaceutically important target with a highly flexible binding site, obtaining computed relative residence times with an accuracy of about 2.3τ for 78% of the compounds and less than 2.0τ within congeneric series. Analysis of dissociation trajectories reveals features that affect ligand unbinding rates, including transient polar interactions and steric hindrance. These results suggest that τRAMD will be widely applicable as a computationally efficient aid to improving drug residence times during lead optimization.

Pteridine reductase-1 (PTR1) is a promising drug target for the treatment of trypanosomiasis. We investigated the potential of a previously identified class of thiadiazole inhibitors of Leishmania major PTR1 for activity against Trypanosoma brucei ( Tb ). We solved crystal structures of several Tb PTR1-inhibitor complexes to guide the structure-based design of new thiadiazole derivatives. Subsequent synthesis and enzyme- and cell-based assays confirm new, mid-micromolar inhibitors of Tb PTR1 with low toxicity. In particular, compound 4m , a biphenyl-thiadiazole-2,5-diamine with IC 50 = 16 μM, was able to potentiate the antitrypanosomal activity of the dihydrofolate reductase inhibitor methotrexate (MTX) with a 4.1-fold decrease of the EC 50 value. In addition, the antiparasitic activity of the combination of 4m and MTX was reversed by addition of folic acid. By adopting an efficient hit discovery platform, we demonstrate, using the 2-amino-1,3,4-thiadiazole scaffold, how a promising tool for the development of anti- T. brucei agents can be obtained.

Computer simulations of molecular systems often make use of regular rectangular grids with equidistant spacing to store information on their molecular interaction fields, e.g., electrostatic potential. These grids provide an easy way to store the data as they do not require any particular specification of the structure of the data. However, such grids may easily become large, and the storage and memory demands may become so high that calculations become infeasible. To overcome this problem, we show here how the data structure DT-Grid can be adapted and applied to efficiently represent macromolecular interaction grids by exploiting the nonuniformity of information on the grid; at the same time, this data structure ensures fast random data access. We demonstrate use of the DT-Grid data structure for potential of mean force and Brownian dynamics simulations of protein-surface binding and macromolecular association with the SDA software. We further demonstrate that the DT-Grid structure enables systems of large size, such as a viral capsid, and high resolution grids to be handled that are beyond current computational feasibility.

Flavonoids represent a potential source of new antitrypanosomatidic leads. Starting from a library of natural products, we combined target-based screening on pteridine reductase 1 with phenotypic screening on Trypanosoma brucei for hit identification. Flavonols were identified as hits, and a library of 16 derivatives was synthesized. Twelve compounds showed EC50 values against T. brucei below 10 μM. Four X-ray crystal structures and docking studies explained the observed structure-activity relationships. Compound 2 (3,6-dihydroxy-2-(3-hydroxyphenyl)-4H-chromen-4-one) was selected for pharmacokinetic studies. Encapsulation of compound 2 in PLGA nanoparticles or cyclodextrins resulted in lower in vitro toxicity when compared to the free compound. Combination studies with methotrexate revealed that compound 13 (3-hydroxy-6-methoxy-2-(4-methoxyphenyl)-4H-chromen-4-one) has the highest synergistic effect at concentration of 1.3 μM, 11.7-fold dose reduction index and no toxicity toward host cells. Our results provide the basis for further chemical modifications aimed at identifying novel antitrypanosomatidic agents showing higher potency toward PTR1 and increased metabolic stability.

Simulations of the long-time scale motions of a ligand binding pocket in a protein may open up new perspectives for the design of compounds with steric or chemical properties differing from those of known binders. However, slow motions of proteins are difficult to access using standard molecular dynamics (MD) simulations and are thus usually neglected in computational drug design. Here, we introduce two nonequilibrium MD approaches to identify conformational changes of a binding site and detect transient pockets associated with these motions. The methods proposed are based on the rotamerically induced perturbation (RIP) MD approach, which employs perturbation of side-chain torsional motion for initiating large-scale protein movement. The first approach, Langevin-RIP (L-RIP), entails a series of short Langevin MD simulations, each starting with perturbation of one of the side-chains lining the binding site of interest. L-RIP provides extensive sampling of conformational changes of the binding site. In less than 1 ns of MD simulation with L-RIP, we observed distortions of the α-helix in the ATP binding site of HSP90 and flipping of the DFG loop in Src kinase. In the second approach, RIPlig, a perturbation is applied to a pseudoligand placed in different parts of a binding pocket, which enables flexible regions of the binding site to be identified in a small number of 10 ps MD simulations. The methods were evaluated for four test proteins displaying different types and degrees of binding site flexibility. Both methods reveal all transient pocket regions in less than a total of 10 ns of simulations, even though many of these regions remained closed in 100 ns conventional MD. The proposed methods provide computationally efficient tools to explore binding site flexibility and can aid in the functional characterization of protein pockets, and the identification of transient pockets for ligand design.

Fluorescent labels are often attached to proteins to monitor binding and adsorption processes. Docking simulations for native hen egg white lysozyme (HEWL) and HEWL labeled with fluorescein isothiocyanate show that these adsorb differently on charged surfaces. Attachment of even a small label can significantly change the interaction properties of a protein. Thus, the results of experiments with fluorescently labeled proteins should be interpreted by modeling the structures and computing the interaction properties of both labeled and unlabeled species.

Protein-nanoparticle associations have important applications in nanoscience and nanotechnology such as targeted drug delivery and theranostics. However, the mechanisms by which proteins recognize nanoparticles and the determinants of specificity are still poorly understood at the microscopic level. Gold is a promising material in nanoparticles for nanobiotechnology applications because of the ease of its functionalization and its tunable optical properties. Ubiquitin is a small, cysteine-free protein (ubiquitous in eukaryotes) whose binding to gold nanoparticles has been characterized recently by nuclear magnetic resonance (NMR). To reveal the molecular basis of these protein-nanoparticle interactions, we performed simulations at multiple levels (ab initio quantum mechanics, classical molecular dynamics and Brownian dynamics) and compared the results with experimental data (circular dichroism and NMR). The results provide a model of the ensemble of structures constituting the ubiquitin-gold surface complex, and insights into the driving forces for the binding of ubiquitin to gold nanoparticles, the role of nanoparticle surfactants (citrate) in the association process, and the origin of the perturbations in the NMR chemical shifts.

Folate analogue inhibitors of Leishmania major pteridine reductase (PTR1) are potential antiparasitic drug candidates for combined therapy with dihydrofolate reductase (DHFR) inhibitors. To identify new molecules with specificity for PTR1, we carried out a virtual screening of the Available Chemicals Directory (ACD) database to select compounds that could interact with L. major PTR1 but not with human DHFR. Through two rounds of drug discovery, we successfully identified eighteen drug-like molecules with low micromolar affinities and high in vitro specificity profiles. Their efficacy against Leishmania species was studied in cultured cells of the promastigote stage, using the compounds both alone and in combination with 1 (pyrimethamine; 5-(4-chlorophenyl)-6-ethylpyrimidine-2,4-diamine). Six compounds showed efficacy only in combination. In toxicity tests against human fibroblasts, several compounds showed low toxicity. One compound, 5c (riluzole; 6-(trifluoromethoxy)-1,3-benzothiazol-2-ylamine), a known drug approved for CNS pathologies, was active in combination and is suitable for early preclinical evaluation of its potential for label extension as a PTR1 inhibitor and antiparasitic drug candidate.

In order to study protein-inorganic surface association processes, we have developed a physics-based energy model, the ProMetCS model, which describes protein-surface interactions at the atomistic level while treating the solvent as a continuum. Here, we present an approach to modeling the interaction of a protein with an atomically flat Au(111) surface in an aqueous solvent. Protein-gold interactions are modeled as the sum of van der Waals, weak chemisorption, and electrostatic interactions, as well as the change in free energy due to partial desolvation of the protein and the metal surface upon association. This desolvation energy includes the effects of water-protein, water-surface, and water-water interactions and has been parametrized using molecular dynamics (MD) simulations of water molecules and a test atom at a gold-water interface. The proposed procedure for computing the energy terms is mostly grid-based and is therefore efficient for application to long-time simulations of protein binding processes. The approach was tested for capped amino acid residues whose potentials of mean force for binding to a gold surface were computed and compared with those obtained previously in MD simulations with water treated explicitly. Calculations show good quantitative agreement with the results from MD simulations for all but one amino acid (Trp), as well as correspondence with available experimental data on the adhesion properties of amino acids.

In protein unfolding simulations, elevated temperature, significantly exceeding the melting temperature Tm, provides an important means to accelerate unfolding to a computationally accessible time range. This procedure is based on the assumption that protein thermal unfolding has Arrhenius behavior and therefore that increasing temperature does not alter the protein unfolding pathways. However, in nature, proteins can show non-Arrhenius behavior and, in practice, overly fast unfolding in high-temperature simulations can result in difficulties in identifying unfolding intermediates and distinguishing their relative stabilities. In this paper, we describe simulations of two WW domains, small protein domains that have a three-stranded β-sheet structure. Simulations were carried out at several temperatures ranging from 300 K to 500 K, starting from folded structures. The results demonstrate the temperature dependence of the unfolding pathways, showing that to obtain unfolding pathways corresponding to those observed in experiments, the elevation of the simulation temperature has to be controlled. Based on trajectory analysis, we proposed a qualitative criterion for judging when an elevated temperature is acceptable or not, namely, that the temperature must be such that the native folded state is sampled substantially before protein unfolding begins. While depending on force field parameters and protein fold complexity, this criterion can be quantified to obtain the upper bound of an "acceptable elevated temperature", which was observed to be dependent on the thermostabilities of the two WW domain proteins.

The secondary structure propensities observed in protein simulations depend heavily on the force field parameters used. The existing empirical force fields often have difficulty in balancing the relative stabilities of helical and extended conformations. The resultant secondary structure bias may not be apparent in short simulations at room temperature starting from the native folded states. However, it can manifest itself dramatically at high temperatures and lead to large deviations from experimentally observed secondary structure propensities. Motivated by thermal unfolding simulations of several WW domains, which have a three-stranded β-sheet structure, we chose the FBP28 WW domain as a well-characterized system to investigate several AMBER force fields as well as parametrization of the NPSA (Neutralized, Polarized ionizable side chains with a solvent-accessible Surface Area-dependent term) implicit solvent model. The ff94 force field and two variants with altered parameters for the backbone torsion term were found to convert the native β-sheet structure directly to a single helix at high temperatures, whereas the ff96 force field produced significant non-native β-sheet content at high temperatures. The ff03 force field was able to reproduce the β-sheet-coil transition and experimentally observed unfolding pathways with both an explicit water solvent and the NPSA implicit solvent model at relatively low temperatures. However, the protein domain became predominantly helical after unfolding. Modification of the solvation parameter in the NPSA implicit solvent model was not sufficient to remedy this problem. The results imply that the intrinsic secondary structure bias in a force field cannot easily be solved by modifying a single parameter such as backbone torsion potential or a solvation parameter of a solvent model. Nevertheless, the results show that the AMBER ff03 force field together with an explicit solvent model or the NPSA implicit solvent model is a useful tool for studying the unfolding of both α- and β-sheet structure protein domains, and an integrative consideration of all force field parameters is likely to be necessary for a complete solution.

The activity of a biological compound is dependent both on specific binding to a target receptor and its ADME (Absorption, Distribution, Metabolism, Excretion) properties. A challenge to predict biological activity is to consider both contributions simultaneously in deriving quantitative models. We present a novel approach to derive QSAR models combining similarity analysis of molecular interaction fields (MIFs) with prediction of logP and/or logD. This new classification method is applied to a set of about 100 compounds related to the auxin plant hormone. The classification based on similarity of their interaction fields is more successful for the indole than the phenoxy compounds. The classification of the phenoxy compounds is however improved by taking into account the influence of the logP and/or the logD values on biological activity. With the new combined method, the majority (8 out of 10) of the previously misclassified derivatives of phenoxy acetic acid are classified in accord with their bioassays. The recently determined crystal structure of the auxin-binding protein 1 (ABP1) enabled validation of our approach. The results of docking a few auxin related compounds with different biological activity to ABP1 correlate well with the classification based on similarity of MIFs only. Biological activity is, however, better predicted by a combined similarity of MIFs + logP/logD approach.

The periplasmic oligopeptide binding component (OppA) of the oligopeptide permease found in Gram-negative bacteria acts as a receptor for peptide transport across the cell membrane and is a potential target for antibacterial drug design. OppA exhibits broad specificity, binding to diverse peptides of 2-5 amino acid residues length. Crystallographic and calorimetric measurements have been carried out by Tame et al. of the binding of 28 peptides of sequence K-X-K to OppA, where X is a natural or nonnatural amino acid. Despite this extensive experimental characterization, a clear relationship between structural and thermodynamic parameters could not be readily identified, with a complicating factor being the observation of varying numbers of water molecules at the binding interface in the different complexes. Consequently, we have applied COMparative BINding Energy (COMBINE) analysis to derive quantitative structure-activity relationships (QSARs) for these 28 OppA-tripeptide complexes. This is the first application of COMBINE analysis to predict binding enthalpies and entropies, and predictive QSAR models were obtained for these quantities as well as for binding free energies. These QSAR models highlight several protein residues and bound water molecules in the binding site, as well as the electrostatic desolvation energies of the protein and the peptides, as responsible for most of the differences in binding thermodynamics between the peptides studied. The QSAR models aid rationalization of the determinants of binding affinity of the OppA:peptide complexes and provide guides for further ligand design. This study also points to the general applicability of COMBINE analysis to estimating thermodynamic parameters for protein-peptide complexes.