Predicted Binding Free Energies Research Articles

Flavonoids as a scaffold for development of novel anti-angiogenic agents: An experimental and computational enquiry

Relationship between structural diversity and biological activities of flavonoids has remained an important discourse in the mainstream of flavonoid research. In the current study anti-angiogenic, cytotoxic, antioxidant and cyclooxygenase (COX) inhibitory activities of diverse class of flavonoids including hydroxyl and methoxy substituted flavones, flavonones and flavonols have been evaluated in the light of developing flavonoids as a potential scaffold for designing novel anti-antiangiogenic agents. We demonstrate anti-angiogenic potential of flavonoids using in vivo chorioallantoic membrane model (CAM) and further elaborate the possible structural reasoning behind observed anti-angiogenic effect using in silico methods. Additionally, we report antioxidant potential and kinetics of free radical scavenging activity using DPPH and SOR scavenging assays. Current study indicates that selected flavonoids possess considerable COX inhibition potential. Furthermore, we describe cytotoxicity of flavonoids against selected cancer cell lines using MTT cell viability assay. Structural analysis of in silico docking poses and predicted binding free energy values are not only in accordance with the experimental anti-angiogenic CAM values from this study but also are in agreement with the previously reported literature on crystallographic data concerning EGFR and VEGFR inhibition.

Problems of robustness in Poisson-Boltzmann binding free energies.

Although models based on the Poisson–Boltzmann (PB) equation have been fairly successful at predicting some experimental quantities, such as solvation free energies (ΔG), these models have not been consistently successful at predicting binding free energies (ΔΔG). Here we found that ranking a set of protein–protein complexes by the electrostatic component (ΔΔGel) of ΔΔG was more difficult than ranking the same molecules by the electrostatic component (ΔGel) of ΔG. This finding was unexpected because ΔΔGel can be calculated by combining estimates of ΔGel for the complex and its components with estimates of the ΔΔGel in vacuum. One might therefore expect that if a theory gave reliable estimates of ΔGel, then its estimates of ΔΔGel would be reliable. However, ΔΔGel for these complexes were orders of magnitude smaller than ΔGel, so although estimates of ΔGel obtained with different force fields and surface definitions were highly correlated, similar estimates of ΔΔGel were often not correlated.

Structural and dynamic basis of acid amido synthetase GH3.1: an investigation of substrate selectivity and major active site access channels.

Auxin/IAA (Indole-3-acetic acid) plays critical roles in many aspects of plant growth and development. Gretchen Hagen 3.1 (GH3.1) enzyme from grapevine (Vitis vinifera) catalyzes the ATP-dependent conjugation of aspartate to IAA suggesting a significant way to modulate levels of cellular auxins/IAA. It is reported that VvGH3.1 prefers IAA as the substrate than 1-naphthaleneacetic acid (NAA) and benzothiazole-2-oxyacetic acid (BTOA), whereas the detailed interaction mechanism of these substrates remains unclear. In this study, based on the recently reported crystal structure of VvGH3.1 and AtGH3.12, the open form of VvGH3.1 was built. Then combined computational techniques including molecular docking, molecular dynamics (MD) simulations, Molecular Mechanics Generalized Born Surface Area and Poisson-Boltzmann Surface Area (MM-GB/PBSA) methods and active site access channel analysis were utilized to investigate the binding mechanism of three substrates, IAA, NAA and BTOA, and dynamic behaviors of GH3.1 induced by substrate binding. The predicted binding free energy is in agreement with the experimental work with an order of IAA < NAA < BTOA. Key residues interacting with three substrates and residues specifically binding to one of them were identified. The strong hydrogen bond interaction formed between IAA and Ser108, Ser339, Gln560 makes it the best substrate over NAA and BTOA. The dynamic behaviors, especially in the dominant access channels of IAA and NAA bound GH3.1, were found to be different from those in the BTOA-GH3.1 system with larger bottleneck radius and shorter length. This study provided novel insight to understand the substrate selectivity mechanism of GH3.1 and also for the rational design of novel auxin-based herbicides and growth regulators.

The feasibility of an efficient drug design method with high-performance computers.

In this study, we propose a supercomputer-assisted drug design approach involving all-atom molecular dynamics (MD)-based binding free energy prediction after the traditional design/selection step. Because this prediction is more accurate than the empirical binding affinity scoring of the traditional approach, the compounds selected by the MD-based prediction should be better drug candidates. In this study, we discuss the applicability of the new approach using two examples. Although the MD-based binding free energy prediction has a huge computational cost, it is feasible with the latest 10 petaflop-scale computer. The supercomputer-assisted drug design approach also involves two important feedback procedures: The first feedback is generated from the MD-based binding free energy prediction step to the drug design step. While the experimental feedback usually provides binding affinities of tens of compounds at one time, the supercomputer allows us to simultaneously obtain the binding free energies of hundreds of compounds. Because the number of calculated binding free energies is sufficiently large, the compounds can be classified into different categories whose properties will aid in the design of the next generation of drug candidates. The second feedback, which occurs from the experiments to the MD simulations, is important to validate the simulation parameters. To demonstrate this, we compare the binding free energies calculated with various force fields to the experimental ones. The results indicate that the prediction will not be very successful, if we use an inaccurate force field. By improving/validating such simulation parameters, the next prediction can be made more accurate.

Inhibition of acetylcholinesterase by two genistein derivatives: kinetic analysis, molecular docking and molecular dynamics simulation

In this study two genistein derivatives (G1 and G2) are reported as inhibitors of acetylcholinesterase (AChE) and butyrylcholinesterase (BuChE), and differences in the inhibition of AChE are described. Although they differ in structure by a single methyl group, the inhibitory effect of G1 (IC50=264nmol/L) on AChE was 80 times stronger than that of G2 (IC50=21,210nmol/L). Enzyme-kinetic analysis, molecular docking and molecular dynamics (MD) simulations were conducted to better understand the molecular basis for this difference. The results obtained by kinetic analysis demonstrated that G1 can interact with both the catalytic active site and peripheral anionic site of AChE. The predicted binding free energies of two complexes calculated by the molecular mechanics/generalized born surface area (MM/GBSA) method were consistent with the experimental data. The analysis of the individual energy terms suggested that a difference between the net electrostatic contributions (ΔEele+ΔGGB) was responsible for the binding affinities of these two inhibitors. Additionally, analysis of the molecular mechanics and MM/GBSA free energy decomposition revealed that the difference between G1 and G2 originated from interactions with Tyr124, Glu292, Val294 and Phe338 of AChE. In conclusion, the results reveal significant differences at the molecular level in the mechanism of inhibition of AChE by these structurally related compounds.

Computational study on the drug resistance mechanism of HCV NS5B RNA-dependent RNA polymerase mutants V494I, V494A, M426A, and M423T to Filibuvir

Filibuvir, a potent non-nucleoside inhibitor of the hepatitis C virus (HCV) NS5B RNA-dependent RNA polymerase (RdRp), has shown great promise in phase IIb clinical trial. However, drug resistant mutations towards Filibuvir have been identified. In the present study, the drug resistance mechanism of wild-type (WT) and mutant NS5B polymerases (including V494I, V494A, M426A, and M423T) toward Filibuvir was investigated by molecular modeling methods. The predicted binding free energy of these five complexes is highly consistent with the experimental EC50 values of Filibuvir to the wild-type and mutant NS5B RdRps, V494I<WT<V494A<M426A<M423T. Analysis of the individual energy terms indicates that the loss of binding affinity is mainly attributed to the decrease in the van der Waals interaction contribution. Through detailed analysis of the interaction between FBV and RdRpV494I, RdRpWT, RdRpV494A, RdRpM426A, and RdRpM423T, several conclusions are made. Firstly, the smaller size of residue 494 side chain results in the smaller binding affinity between FBV and RdRp. Secondly, the poor inhibition capacity of Filibuvir toward RdRpM426A is mainly due to the decrease in the van der Walls interaction between Filibuvir and residue Leu-497M426A caused by the spatial structure change of Ala-426M426A. Thirdly, the decrease in the binding affinity in mutation M423T is attributed to the smaller binding cave and the cyclopentyl group of Filibuvir exposing outside the cave. Our computational results will provide valuable information for developing more potent and selective inhibitors toward HCV NS5B polymerase.

Revealing origin of decrease in potency of darunavir and amprenavir against HIV-2 relative to HIV-1 protease by molecular dynamics simulations.

Clinical inhibitors Darunavir (DRV) and Amprenavir (APV) are less effective on HIV-2 protease (PR2) than on HIV-1 protease (PR1). To identify molecular basis associated with the lower inhibition, molecular dynamics (MD) simulations and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) calculations were performed to investigate the effectiveness of the PR1 inhibitors DRV and APV against PR1/PR2. The rank of predicted binding free energies agrees with the experimental determined one. Moreover, our results show that two inhibitors bind less strongly to PR2 than to PR1, again in agreement with the experimental findings. The decrease in binding free energies for PR2 relative to PR1 is found to arise from the reduction of the van der Waals interactions induced by the structural adjustment of the triple mutant V32I, I47V and V82I. This result is further supported by the difference between the van der Waals interactions of inhibitors with each residue in PR2 and in PR1. The results from the principle component analysis suggest that inhibitor binding tends to make the flaps of PR2 close and the one of PR1 open. We expect that this study can theoretically provide significant guidance and dynamics information for the design of potent dual inhibitors targeting PR1/PR2.

Seeking potential anticonvulsant agents that target GABAA receptors using experimental and theoretical procedures.

The aim of this study was to identify compounds that possess anticonvulsant activity by using a pentylenetetrazol (PTZ)-induced seizure model. Theoretical studies of a set of ligands, explored the binding affinities of the ligands for the GABA(A) receptor (GABA(A)R), including some benzodiazepines. The ligands satisfy the Lipinski rules and contain a pharmacophore core that has been previously reported to be a GABA(A)R activator. To select the ligands with the best physicochemical properties, all of the compounds were analyzed by quantum mechanics and the energies of the highest occupied molecular orbital and lowest unoccupied molecular orbital were determined. Docking calculations between the ligands and the GABA(A)R were used to identify the complexes with the highest Gibbs binding energies. The identified compound D1 (dibenzo(b,f)(1,4)diazocine-6,11(5H,12H)-dione) was synthesized, experimentally tested, and the GABA(A)R-D1 complex was submitted to 12-ns-long molecular dynamics (MD) simulations to corroborate the binding conformation obtained by docking techniques. MD simulations were also used to analyze the decomposition of the Gibbs binding energy of the residues involved in the stabilization of the complex. To validate our theoretical results, molecular docking and MD simulations were also performed for three reference compounds that are currently in commercial use: clonazepam (CLZ), zolpidem and eszopiclone. The theoretical results show that the GABA(A)R-D1, and GABA(A)R-CLZ complexes bind to the benzodiazepine binding site, share a similar map of binding residues, and have similar Gibbs binding energies and entropic components. Experimental studies using a PTZ-induced seizure model showed that D1 possesses similar activity to CLZ, which corroborates the predicted binding free energy identified by theoretical calculations.

Ligand-receptor affinities computed by an adapted linear interaction model for continuum electrostatics and by protein conformational averaging.

Accurate calculations of free energies involved in small-molecule binding to a receptor are challenging. Interactions between ligand, receptor, and solvent molecules have to be described precisely, and a large number of conformational microstates has to be sampled, particularly for ligand binding to a flexible protein. Linear interaction energy models are computationally efficient methods that have found considerable success in the prediction of binding free energies. Here, we parametrize a linear interaction model for implicit solvation with coefficients adapted by ligand and binding site relative polarities in order to predict ligand binding free energies. Results obtained for a diverse series of ligands suggest that the model has good predictive power and transferability. We also apply implicit ligand theory and propose approximations to average contributions of multiple ligand-receptor poses built from a protein conformational ensemble and find that exponential averages require proper energy discrimination between plausible binding poses and false-positives (i.e., decoys). The linear interaction model and the averaging procedures presented can be applied independently of each other and of the method used to obtain the receptor structural representation.

Computational study on the molecular mechanisms of drug resistance of Narlaprevir due to V36M, R155K, V36M+R155K, T54A, and A156T mutations of HCV NS3/4A protease.

Narlaprevir is a novel NS3/4A protease inhibitor of hepatitis C virus (HCV), and it has been tested in a phase II clinical trial recently. However, distinct drug-resistance of Narlaprevir has been discovered. In our study, the molecular mechanisms of drug-resistance of Narlaprevir due to the mutations V36M, R155K, V36M+R155K, T54A, and A156T of NS3/4A protease have been investigated by molecular dynamics (MD) simulations, free energy calculations, and free energy decomposition analysis. The predicted binding free energies of Narlaprevir towards the wild-type and five mutants show that the mutations V36M, R155K, and T54A lead to low-level drug resistance and the mutations V36M+R155K and A156T lead to high-level drug resistance, which is consistent with the experimental data. The analysis of the individual energy terms indicates that the van der Waals contribution is important for distinguishing the binding affinities of these six complexes. These findings again show that the combination of different molecular modeling techniques is an efficient way to interpret the molecular mechanism of drug-resistance. Our work mainly elaborates the molecular mechanism of drug-resistance of Narlaprevir and further provides valuable information for developing novel, safer, and more potent HCV antiviral drugs in the near future.

Comparative Analysis of Various Electrostatic Potentials on Docking Precision Against Cyclin-Dependent Kinase 2 Protein: A Multiple Docking Approach

The fundamental of molecular modeling is the interaction and binding to form a complex, because it explains the action of most drugs to a receptor active site. In the present study, different semiempirical (RM1, AM1, PM3, MNDO) and ab initio (HF, DFT) charge models were investigated for their performance in prediction of docking pose against CDK2 proteins with their respective inhibitor. Further, multiple docking approaches and Prime/MM-GBSA calculations were applied to predict the binding mode with respective charge model against CDK2 inhibitors. A reliable docking result was obtained using RRD, which showed significance improvement on ligand binding poses and docking score accuracy to the IFD. The combined use of RRD and Prime/MM-GBSA method could give a high correlation between the predicted binding free energy and experimental biological activity. The preliminary results point out that AM1 could be a precious charge model for design of new drugs with enhanced success rate. As a very similar result was also found for a different system of the protein-ligand binding, the suggested scoring function based on AM1 method seems to be applicable in drug design. The results from this study can provide insights into highest success rate for design of potent and selective CDK2 inhibitors.

Binding Free Energy Calculations for Lead Optimization: Assessment of Their Accuracy in an Industrial Drug Design Context.

Correctly ranking compounds according to their computed relative binding affinities will be of great value for decision making in the lead optimization phase of industrial drug discovery. However, the performance of existing computationally demanding binding free energy calculation methods in this context is largely unknown. We analyzed the performance of the molecular mechanics continuum solvent, the linear interaction energy (LIE), and the thermodynamic integration (TI) approach for three sets of compounds from industrial lead optimization projects. The data sets pose challenges typical for this early stage of drug discovery. None of the methods was sufficiently predictive when applied out of the box without considering these challenges. Detailed investigations of failures revealed critical points that are essential for good binding free energy predictions. When data set-specific features were considered accordingly, predictions valuable for lead optimization could be obtained for all approaches but LIE. Our findings lead to clear recommendations for when to use which of the above approaches. Our findings also stress the important role of expert knowledge in this process, not least for estimating the accuracy of prediction results by TI, using indicators such as the size and chemical structure of exchanged groups and the statistical error in the predictions. Such knowledge will be invaluable when it comes to the question which of the TI results can be trusted for decision making.

Insights from QM/MM Modeling the 3D Structure of the 2009 H1N1 Influenza A Virus Neuraminidase and Its Binding Interactions with Antiviral Drugs.

Control of drug release through the inhibition of neuraminidase-1 has been identified as a potential target for the treatment of H1N1 influenza; however, the drug binding mode of neuraminidase is not yet completely understood. In this work, we propose a model for a neuraminidase-1 complex based on four known X-ray structures of drug/neuraminidase-1 complexes. Specifically, H1N1 neuraminidase-1 complexed with 4 drugs (zanamivir, laninamivir, laninamivir octanoate and oseltamivir) was first investigated using a combined quantum mechanical and molecular mechanical (QM/MM) approach. Based on these structures, a model for the H1N1 neuraminidase-1 complex was proposed and simulated using the same computational protocol. Implications to drug/H1N1 neuraminidase-1 binding modes are discussed. From our simulations, the predicted binding free energies of the four drugs are in good agreement with the experimental results, with the correlation coefficient being 0.84.

Molecular Recognition in Complexes of TRF Proteins with Telomeric DNA

Telomeres are specialized nucleoprotein assemblies that protect the ends of linear chromosomes. In humans and many other species, telomeres consist of tandem TTAGGG repeats bound by a protein complex known as shelterin that remodels telomeric DNA into a protective loop structure and regulates telomere homeostasis. Shelterin recognizes telomeric repeats through its two major components known as Telomere Repeat-Binding Factors, TRF1 and TRF2. These two homologous proteins are therefore essential for the formation and normal function of telomeres. Indeed, TRF1 and TRF2 are implicated in a plethora of different cellular functions and their depletion leads to telomere dysfunction with chromosomal fusions, followed by apoptotic cell death. More specifically, it was found that TRF1 acts as a negative regulator of telomere length, and TRF2 is involved in stabilizing the loop structure. Consequently, these proteins are of great interest, not only because of their key role in telomere maintenance and stability, but also as potential drug targets. In the current study, we investigated the molecular basis of telomeric sequence recognition by TRF1 and TRF2 and their DNA binding mechanism. We used molecular dynamics (MD) to calculate the free energy profiles for binding of TRFs to telomeric DNA. We found that the predicted binding free energies were in good agreement with experimental data. Further, different molecular determinants of binding, such as binding enthalpies and entropies, the hydrogen bonding pattern and changes in surface area, were analyzed to decompose and examine the overall binding free energies at the structural level. With this approach, we were able to draw conclusions regarding the consecutive stages of sequence-specific association, and propose a novel aspartate-dependent mechanism of sequence recognition. Finally, our work demonstrates the applicability of computational MD-based methods to studying protein-DNA interactions.

Bringing Clarity to the Prediction of Protein-Ligand Binding Free Energies via "Blurring"

We present a method to evaluate the free energies of ligand binding utilizing a Monte Carlo estimation of the configuration integrals concomitant with uncertainty quantification. Ensembles for integration are built through systematically perturbing an initial ligand conformation in a rigid binding pocket, which is optimized separately prior to incorporation of the ligand. We call the procedure producing the ensembles “blurring”, and it is carried out using an in-house developed code. The Boltzmann factor contribution of each pose to the configuration integral is computed and from there the free energy is obtained. Potential function uncertainties are estimated using a fragment-based error propagation method. This method has been applied to a set of small aromatic ligands complexed with T4 Lysozyme L99A mutant. Microstate energies have been determined with the force fields ff99SB and ff94, and the semiempirical method PM6DH2 in conjunction with continuum solvation models including Generalized Born (GB), the Conductor-like Screening Model (COSMO), and SMD. Of the methods studied, PM6DH2-based scoring gave binding free energy estimates, which yielded a good correlation to the experimental binding affinities (R2 = 0.7). All methods overestimated the calculated binding affinities. We trace this to insufficient sampling, the single static protein structure, and inaccuracies in the solvent models we have used in this study.

English

Dengue has become a major global health threat, especially in tropical and subtropical regions. The development of antiviral agent targeting viral replication is really needed at this time. NS5 methyltransferase presents as a novel antiviral target. This enzyme plays an important role in the methylation of 5'-cap mRNA. Inhibition of the NS5 methyltransferase could inhibit dengue virus replication. In this research, two sites of NS5 methyltransferase (S-Adenosyl methionine/SAM binding site and RNA-cap site) were used as targets for inhibition. As much as 300 commercial cyclic peptides were screened to these target sites by means of molecular docking. Analysis of ligand-enzyme binding free energy and pharmacological prediction revealed two best ligands, namely [Tyr123] Prepro Endothelin (110-130), amide, human and Urotensin II, human. According to molecular dynamic simulation, both ligands maintain a stable complex conformation between enzyme and ligand at temperature 310 K and 312 K. Hence, Urotensin II, human is more reactive at 312 K than at 310 K. However, both ligands can be used as potential inhibitor candidates against NS5 methyltransferase of dengue virus with Urotensin II, human exposes more promising activity at 312 K.

Theoretical Studies on Binding and Specificity Mechanisms of Farnesyltransferase (FTase) and Geranylgeranyltransferase Type-I (GGTase-I) Inhibitors by Molecular Modeling

Farnesyltransferase (FTase) and geranylgeranyltransferase type-I (GGTase-I) are two members of protein prenyltransferases, which play critical roles in lipid post-translational modifications. Potent inhibitors of FTase and GGTase-I have been confirmed to show favorable influence on the therapies of various diseases, such as cancers, malaria and Toxoplasmosis. However, designing highly specific inhibitors toward FTase or GGTase-I without influencing their binding affinity remains a big challenge. In this work, molecular docking, molecular dynamics (MD) simulations and MM/GBSA free energy calculations were employed to study the bindings of two highly selective inhibitors (lonafarnib and GGTI-2133) towards FTase or GGTase-I. The specificities of the studied inhibitors derived from the predicted binding free energies are consistent with the experimental data. The analysis of the energetic components illustrates that both the non-polar and polar interactions play critical roles in determining the specificity between FTase and GGTase-I. Moreover, the protein-inhibitor interaction spectra for the studied inhibitors were determined through the decomposition of the binding free energies, and the important residues for binding and specificity were highlighted. Our study provides useful information for the rational design of selective FTase or GGTase-I inhibitors.

Theoretical studies on beta and delta isoform-specific binding mechanisms of phosphoinositide 3-kinase inhibitors

Phosphoinositide 3-kinase (PI3K) is known to be closely related to tumorigenesis and cell proliferation, and controls a variety of cellular processes, including proliferation, growth, apoptosis, migration, metabolism, etc. The PI3K family comprises eight catalytic isoforms, which are subdivided into three classes. Recently, the discovery of inhibitors that block a single isoform of PI3K has continued to attract special attention because they may have higher selectivity for certain tumors and less toxicity for healthy cells. The PI3Kβ and PI3Kδ share fewer studies than α/γ, and therefore, in this work, the combination of molecular dynamics simulations and free energy calculations was employed to explore the binding of three isoform-specific PI3K inhibitors (COM8, IC87114, and GDC-0941) to PI3Kβ or PI3Kδ. The isoform specificities of the studied inhibitors derived from the predicted binding free energies are in good agreement with the experimental data. In addition, the key residues critical for PI3Kβ or PI3Kδ selectivity were highlighted by decomposing the binding free energies into the contributions from individual residues. It was observed that although PI3Kβ and PI3Kδ share the conserved ATP-binding pockets, individual residues do behave differently, particularly the residues critical for PI3Kβ or PI3Kδ selectivity. It can be concluded that the inhibitor specificity between PI3Kβ and PI3Kδ is determined by the additive contributions from multiple residues, not just a single one. This study provides valuable information for understanding the isoform-specific binding mechanisms of PI3K inhibitors, and should be useful for the rational design of novel and selective PI3K inhibitors.

Molecular dynamics simulation-based evaluation of the binding free energies of computationally designed drug candidates: importance of the dynamical effects.

The computational structure-based drug design (SBDD) mainly aims at generating or discovering new chemical compounds with sufficiently large binding free energy. In any de novo drug design methods and virtual screening methods, drug candidates are selected by approximately evaluating the binding free energy (or the binding affinity). This approximate binding free energy, usually called "empirical score," is critical to the success of the SBDD. The purpose of this work is to yield physical insight into the approximate evaluation method in comparison with an exact molecular dynamics (MD) simulation-based method (named MP-CAFEE), which can predict binding free energies accurately. We calculate the binding free energies for 58 selected drug candidates with MP-CAFEE. Here, the compounds are generated by OPMF, a novel fragment-based de novo drug design method, and the ligand-protein interaction energy is used as an empirical score. The results show that the correlation between the binding free energy and the interaction energy is not strong enough to clearly distinguish compounds with nM-affinity from those with µM-affinity. This implies that it is necessary to take into account the natural protein motion with explicitly surrounded by water molecules to improve the efficiency of the drug candidate selection procedure.

Simultaneous prediction of binding free energy and specificity for PDZ domain–peptide interactions

Interactions between protein domains and linear peptides underlie many biological processes. Among these interactions, the recognition of C-terminal peptides by PDZ domains is one of the most ubiquitous. In this work, we present a mathematical model for PDZ domain-peptide interactions capable of predicting both affinity and specificity of binding based on X-ray crystal structures and comparative modeling with ROSETTA. We developed our mathematical model using a large phage display dataset describing binding specificity for a wild type PDZ domain and 91 single mutants, as well as binding affinity data for a wild type PDZ domain binding to 28 different peptides. Structural refinement was carried out through several ROSETTA protocols, the most accurate of which included flexible peptide docking and several iterations of side chain repacking and backbone minimization. Our findings emphasize the importance of backbone flexibility and the energetic contributions of side chain-side chain hydrogen bonds in accurately predicting interactions. We also determined that predicting PDZ domain-peptide interactions became increasingly challenging as the length of the peptide increased in the N-terminal direction. In the training dataset, predicted binding energies correlated with those derived through calorimetry and specificity switches introduced through single mutations at interface positions were recapitulated. In independent tests, our best performing protocol was capable of predicting dissociation constants well within one order of magnitude of the experimental values and specificity profiles at the level of accuracy of previous studies. To our knowledge, this approach represents the first integrated protocol for predicting both affinity and specificity for PDZ domain-peptide interactions.