Artificial protein cavities as specific ligand-binding templates: characterization of an engineered heterocyclic cation-binding site that preserves the evolved specificity of the parent protein

Hydrophobicity-Based Force Field In Enzymes

The biocatalysis process takes place with the participation of enzymes, which, depending on the reaction carried out, require, apart from the appropriate arrangement of catalytic residues, an appropriate external force field. It is generated by the protein body. The relatively small size of the part directly involved in the process itself is supported by the presence of an often complex structure of the protein body, the purpose of which is to provide an appropriate local force field, eliminating the influence of water. Very often, the large size of the enzyme is an expression of the complex form of this field. In this paper, a comparative analysis of arbitrarily selected enzymes, representatives of different enzyme classes, was carried out, focusing on the measurement of the diversity of the force field provided by a given protein. This analysis was based on the fuzzy oil drop model (FOD) and its modified version (FOD-M), which takes into account the participation of nonaqueous external factors in shaping the structure and thus the force field within the protein. The degree and type of ordering of the hydrophobicity distribution in the protein molecule is the result of the influence of the environment but also the supplier of the local environment for a given process, including the catalysis process in particular. Determining the share of a nonaqueous environment is important due to the ubiquity of polar water, whose participation in processes with high specificity requires control. It can be assumed that some enzymes in their composition have a permanently built-in part, the role of which is reduced to that of a permanent chaperone. It provides a specific external force field needed for the process. The proposed model, generalized to other types of proteins, may also provide a form of recording the environment model for the simulation of the in silico protein folding process, taking into account the impact of its differentiation.

LifeSoaks: a tool for analyzing solvent channels in protein crystals and obstacles for soaking experiments

Due to the structural complexity of proteins, their corresponding crystal arrangements generally contain a significant amount of solvent-occupied space. These areas allow a certain degree of intracrystalline protein flexibility and mobility of solutes. Therefore, knowledge of the geometry of solvent-filled channels and cavities is essential whenever the dynamics inside a crystal are of interest. Especially in soaking experiments for structure-based drug design, ligands must be able to traverse the crystal solvent channels and reach the corresponding binding pockets. Unsuccessful screenings are sometimes attributed to the geometry of the crystal packing, but the underlying causes are often difficult to understand. This work presents LifeSoaks, a novel tool for analyzing and visualizing solvent channels in protein crystals. LifeSoaks uses a Voronoi diagram-based periodic channel representation which can be efficiently computed. The size and location of channel bottlenecks, which might hinder molecular diffusion, can be directly derived from this representation. This work presents the calculated bottleneck radii for all crystal structures in the PDB and the analysis of a new, hand-curated data set of structures obtained by soaking experiments. The results indicate that the consideration of bottleneck radii and the visual inspection of channels are beneficial for planning soaking experiments.

Identifying Protein-Ligand Interactions via a Novel Distance Self-Feedback Biomolecular Interaction Network

Efficient and accurate characterizations of protein-ligand interactions are key to understanding biology at the molecular level. They are particularly useful in pharmaceutical industry applications. They are usually computationally demanding for those widely applied dynamics-based methods in identifying important residues or calculating ligand binding free energy. In this work, we proposed a graph deep learning (DL) framework, namely, the distance self-feedback biomolecular interaction network (DSBIN), in which the relationship between the complex structure and binding affinity can be established by means of a carefully designed distance self-feedback module and interaction layer. Our model can directly provide a quantitative evaluation of inhibitor binding affinities (pK_d). More importantly, the DSBIN model efficiently identifies key interactions for inhibitor binding and thus intrinsically bears the interpretability. Its generalization performance was further verified using 1405 unseen structures. The predicted binding free energies' deviations were calculated to be less than 1.37 kcal/mol for more than 55% structures. Moreover, we also compared the DSBIN model with a commonly used theoretical method in calculating the substrate binding free energy, MM/GBSA. Our results show that the current DL model has generally better performance in predicting the binding free energy. For a specific complex system, mannopentaose/TmCBM27, the DSBIN predicted binding free energy is -8.21 kcal/mol, which is very close to experimentally measured -7.76 kcal/mol and MM/GBSA calculated -7.16 kcal/mol. Meanwhile, all important aromatic residues around the binding pocket can be identified by our DL model. Considering the accuracy and efficiency of the newly developed DL model, it may be very helpful in the field of drug design and molecular recognition.

Why Do Most Aromatics Fail to Support Hole Hopping in the Cytochrome c Peroxidase-Cytochrome c Complex?

Electron transport through aromatic species (especially tryptophan and tyrosine) plays a central role in water splitting, redox signaling, oxidative damage protection, and bioenergetics. The cytochrome c peroxidase (CcP)-cytochrome c (Cc) complex (CcP:Cc) is used widely to study interprotein electron transfer (ET) mechanisms. Tryptophan 191 (Trp191) of CcP supports hole hopping charge recombination in the CcP:Cc complex. Experimental studies find that when Trp191 is substituted by tyrosine, phenylalanine, or redox-active aniline derivatives bound in the W191G cavity, enzymatic activity and charge recombination rates both decrease. Theoretical analysis of these CcP:Cc complexes finds that the ET kinetics depend strongly on the chemistry of the modified Trp site. The computed electronic couplings in the W191F and W191G species are orders of magnitude smaller than in the native protein, due largely to the absence of a hopping intermediate and the large tunneling distance. Small molecules bound in the W191G cavity are weakly coupled electronically to the Cc heme, and the structural disorder of the guest molecule in the binding pocket may contribute further to the lack of enzymatic activity. The couplings in W191Y are not substantially weakened compared to the native species, but the redox potential difference for tyrosine vs tryptophan oxidation accounts for the slower rate in the Tyr mutant. Thus, theoretical analysis explains why only the native Trp supports rapid hole hopping in the CcP:Cc complex. Favorable free energies and electronic couplings are essential for establishing an efficient hole hopping relay in this protein-protein complex.

Exploring Protein Cavities through Rigidity Analysis

The geometry of cavities in the surfaces of proteins facilitates a variety of biochemical functions. To better understand the biochemical nature of protein cavities, the shape, size, chemical properties, and evolutionary nature of functional and nonfunctional surface cavities have been exhaustively surveyed in protein structures. The rigidity of surface cavities, however, is not immediately available as a characteristic of structure data, and is thus more difficult to examine. Using rigidity analysis for assessing and analyzing molecular rigidity, this paper performs the first survey of the relationships between cavity properties, such as size and residue content, and how they correspond to cavity rigidity. Our survey measured a variety of rigidity metrics on 120,323 cavities from 12,785 sequentially non-redundant protein chains. We used VASP-E, a volume-based algorithm for analyzing cavity geometry. Our results suggest that rigidity properties of protein cavities are dependent on cavity surface area.

An Activator-Blocker Pair Provides a Controllable On-Off Switch for a Ketosteroid Isomerase Active Site Mutant

Control of enzyme activity is fundamental to biology and represents a long-term goal in bioengineering and precision therapeutics. While several powerful molecular strategies have been developed, limitations remain in their generalizability and dynamic range. We demonstrate a control mechanism via separate small molecules that turn on the enzyme (activator) and turn off the activation (blocker). We show that a pocket created near the active site base of the enzyme ketosteriod isomerase (KSI) allows efficient and saturable base rescue when the enzyme's natural general base is removed. Binding a small molecule with similar properties but lacking general-base capability in this pocket shuts off rescue. The ability of small molecules to directly participate in and directly block catalysis may afford a broad controllable dynamic range. This approach may be amenable to numerous enzymes and to engineering and screening approaches to identify activators and blockers with strong, specific binding for engineering and therapeutic applications.

Predicting Binding Free Energies: Frontiers and Benchmarks

Binding free energy calculations based on molecular simulations provide predicted affinities for biomolecular complexes. These calculations begin with a detailed description of a system, including its chemical composition and the interactions among its components. Simulations of the system are then used to compute thermodynamic information, such as binding affinities. Because of their promise for guiding molecular design, these calculations have recently begun to see widespread applications in early-stage drug discovery. However, many hurdles remain in making them a robust and reliable tool. In this review, we highlight key challenges of these calculations, discuss some examples of these challenges, and call for the designation of standard community benchmark test systems that will help the research community generate and evaluate progress. In our view, progress will require careful assessment and evaluation of new methods, force fields, and modeling innovations on well-characterized benchmark systems, and we lay out our vision for how this can be achieved.

Constraints on the Radical Cation Center of Cytochrome c Peroxidase for Electron Transfer from Cytochrome c

The tryptophan 191 cation radical of cytochrome c peroxidase (CcP) compound I (Cpd I) mediates long-range electron transfer (ET) to cytochrome c (Cc). Here we test the effects of chemical substitution at position 191. CcP W191Y forms a stable tyrosyl radical upon reaction with peroxide and produces spectral properties similar to those of Cpd I but has low reactivity toward reduced Cc. CcP W191G and W191F variants also have low activity, as do redox ligands that bind within the W191G cavity. Crystal structures of complexes between Cc and CcP W191X (X = Y, F, or G), as well as W191G with four bound ligands reveal similar 1:1 association modes and heme pocket conformations. The ligands display structural disorder in the pocket and do not hydrogen bond to Asp235, as does Trp191. Well-ordered Tyr191 directs its hydroxyl group toward the porphyrin ring, with no basic residue in the range of interaction. CcP W191X (X = Y, F, or G) variants substituted with zinc-porphyrin (ZnP) undergo photoinduced ET with Cc(III). Their slow charge recombination kinetics that result from loss of the radical center allow resolution of difference spectra for the charge-separated state [ZnP(+), Cc(II)]. The change from a phenyl moiety at position 191 in W191F to a water-filled cavity in W191G produces effects on ET rates much weaker than the effects of the change from Trp to Phe. Low net reactivity of W191Y toward Cc(II) derives either from the inability of ZnP(+) or the Fe-CcP ferryl to oxidize Tyr or from the low potential of the resulting neutral Tyr radical.

Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects

The calculation of a protein-ligand binding free energy based on molecular dynamics (MD) simulations generally relies on a thermodynamic cycle in which the ligand is alchemically inserted into the system, both in the solvated protein and free in solution. The corresponding ligand-insertion free energies are typically calculated in nanoscale computational boxes simulated under periodic boundary conditions and considering electrostatic interactions defined by a periodic lattice-sum. This is distinct from the ideal bulk situation of a system of macroscopic size simulated under non-periodic boundary conditions with Coulombic electrostatic interactions. This discrepancy results in finite-size effects, which affect primarily the charging component of the insertion free energy, are dependent on the box size, and can be large when the ligand bears a net charge, especially if the protein is charged as well. This article investigates finite-size effects on calculated charging free energies using as a test case the binding of the ligand 2-amino-5-methylthiazole (net charge +1 e) to a mutant form of yeast cytochrome c peroxidase in water. Considering different charge isoforms of the protein (net charges -5, 0, +3, or +9 e), either in the absence or the presence of neutralizing counter-ions, and sizes of the cubic computational box (edges ranging from 7.42 to 11.02 nm), the potentially large magnitude of finite-size effects on the raw charging free energies (up to 17.1 kJ mol(-1)) is demonstrated. Two correction schemes are then proposed to eliminate these effects, a numerical and an analytical one. Both schemes are based on a continuum-electrostatics analysis and require performing Poisson-Boltzmann (PB) calculations on the protein-ligand system. While the numerical scheme requires PB calculations under both non-periodic and periodic boundary conditions, the latter at the box size considered in the MD simulations, the analytical scheme only requires three non-periodic PB calculations for a given system, its dependence on the box size being analytical. The latter scheme also provides insight into the physical origin of the finite-size effects. These two schemes also encompass a correction for discrete solvent effects that persists even in the limit of infinite box sizes. Application of either scheme essentially eliminates the size dependence of the corrected charging free energies (maximal deviation of 1.5 kJ mol(-1)). Because it is simple to apply, the analytical correction scheme offers a general solution to the problem of finite-size effects in free-energy calculations involving charged solutes, as encountered in calculations concerning, e.g., protein-ligand binding, biomolecular association, residue mutation, pKa and redox potential estimation, substrate transformation, solvation, and solvent-solvent partitioning.

Blind prediction of charged ligand binding affinities in a model binding site

Predicting absolute protein-ligand binding affinities remains a frontier challenge in ligand discovery and design. This becomes more difficult when ionic interactions are involved because of the large opposing solvation and electrostatic attraction energies. In a blind test, we examined whether alchemical free-energy calculations could predict binding affinities of 14 charged and 5 neutral compounds previously untested as ligands for a cavity binding site in cytochrome c peroxidase. In this simplified site, polar and cationic ligands compete with solvent to interact with a buried aspartate. Predictions were tested by calorimetry, spectroscopy, and crystallography. Of the 15 compounds predicted to bind, 13 were experimentally confirmed, while 4 compounds were false negative predictions. Predictions had a root-mean-square error of 1.95 kcal/mol to the experimental affinities, and predicted poses had an average RMSD of 1.7Å to the crystallographic poses. This test serves as a benchmark for these thermodynamically rigorous calculations at predicting binding affinities for charged compounds and gives insights into the existing sources of error, which are primarily electrostatic interactions inside proteins. Our experiments also provide a useful set of ionic binding affinities in a simplified system for testing new affinity prediction methods.

Roles for ordered and bulk solvent in ligand recognition and docking in two related cavities

A key challenge in structure-based discovery is accounting for modulation of protein-ligand interactions by ordered and bulk solvent. To investigate this, we compared ligand binding to a buried cavity in Cytochrome c Peroxidase (CcP), where affinity is dominated by a single ionic interaction, versus a cavity variant partly opened to solvent by loop deletion. This opening had unexpected effects on ligand orientation, affinity, and ordered water structure. Some ligands lost over ten-fold in affinity and reoriented in the cavity, while others retained their geometries, formed new interactions with water networks, and improved affinity. To test our ability to discover new ligands against this opened site prospectively, a 534,000 fragment library was docked against the open cavity using two models of ligand solvation. Using an older solvation model that prioritized many neutral molecules, three such uncharged docking hits were tested, none of which was observed to bind; these molecules were not highly ranked by the new, context-dependent solvation score. Using this new method, another 15 highly-ranked molecules were tested for binding. In contrast to the previous result, 14 of these bound detectably, with affinities ranging from 8 µM to 2 mM. In crystal structures, four of these new ligands superposed well with the docking predictions but two did not, reflecting unanticipated interactions with newly ordered waters molecules. Comparing recognition between this open cavity and its buried analog begins to isolate the roles of ordered solvent in a system that lends itself readily to prospective testing and that may be broadly useful to the community.

Perspective: Alchemical free energy calculations for drug discovery

Computational techniques see widespread use in pharmaceutical drug discovery, but typically prove unreliable in predicting trends in protein-ligand binding. Alchemical free energy calculations seek to change that by providing rigorous binding free energies from molecular simulations. Given adequate sampling and an accurate enough force field, these techniques yield accurate free energy estimates. Recent innovations in alchemical techniques have sparked a resurgence of interest in these calculations. Still, many obstacles stand in the way of their routine application in a drug discovery context, including the one we focus on here, sampling. Sampling of binding modes poses a particular challenge as binding modes are often separated by large energy barriers, leading to slow transitions. Binding modes are difficult to predict, and in some cases multiple binding modes may contribute to binding. In view of these hurdles, we present a framework for dealing carefully with uncertainty in binding mode or conformation in the context of free energy calculations. With careful sampling, free energy techniques show considerable promise for aiding drug discovery.

Modeling regionalized volumetric differences in protein-ligand binding cavities

Identifying elements of protein structures that create differences in protein-ligand binding specificity is an essential method for explaining the molecular mechanisms underlying preferential binding. In some cases, influential mechanisms can be visually identified by experts in structural biology, but subtler mechanisms, whose significance may only be apparent from the analysis of many structures, are harder to find. To assist this process, we present a geometric algorithm and two statistical models for identifying significant structural differences in protein-ligand binding cavities. We demonstrate these methods in an analysis of sequentially nonredundant structural representatives of the canonical serine proteases and the enolase superfamily. Here, we observed that statistically significant structural variations identified experimentally established determinants of specificity. We also observed that an analysis of individual regions inside cavities can reveal areas where small differences in shape can correspond to differences in specificity.

Pollutant-induced modulation in conformation and β-lactamase activity of human serum albumin

Structural changes in human serum albumin (HSA) induced by the pollutants 1-naphthol, 2-naphthol and 8-quinolinol were analyzed by circular dichroism, fluorescence spectroscopy and dynamic light scattering. The alteration in protein conformational stability was determined by helical content induction (from 55 to 75%) upon protein-pollutant interactions. Domain plasticity is responsible for the temperature-mediated unfolding of HSA. These findings were compared to HSA-hydrolase activity. We found that though HSA is a monomeric protein, it shows heterotropic allostericity for β-lactamase activity in the presence of pollutants, which act as K- and V-type non-essential activators. Pollutants cause conformational changes and catalytic modifications of the protein (increase in β-lactamase activity from 100 to 200%). HSA-pollutant interactions mediate other protein-ligand interactions, such as HSA-nitrocefin. Therefore, this protein can exist in different conformations with different catalytic properties depending on activator binding. This is the first report to demonstrate the catalytic allostericity of HSA through a mechanistic approach. We also show a correlation with non-microbial drug resistance as HSA is capable of self-hydrolysis of β-lactam drugs, which is further potentiated by pollutants due to conformational changes in HSA.

Impact of Proximal and Distal Pocket Site-Directed Mutations on the Ferric/Ferrous Heme Redox Potential of Yeast Cytochrome-c-Peroxidase

Cytochrome-c-peroxidase (CCP) contains a five-coordinate heme active site. The reduction potential for the ferric to ferrous couple in CCP is anomalously low and pH dependent (Eo = ~-180 mV vs. S.H.E. at pH 7). The contribution of the protein environment to the tuning of the redox potential of this couple is evaluated using site directed mutants of several amino acid residues in the environment of the heme. These include proximal pocket mutation to residues Asp-235, Trp-191, Phe-202 and His-175, distal pocket mutation to residues Trp-51, His-52, and Arg-48; and a heme edge mutation to Ala-147. Where unknown, the structural changes resulting from the amino acid substitution have been studied by X-ray crystallography. In most cases, ostensibly polar or charged residues are replaced by large hydrophobic groups or alternatively by Ala or Gly. These latter have been shown to generate large, solvent filled cavities. Reduction potentials are measured as a function of pH by spectroelectrochemistry. Starting with the X-ray derived structures of CCP and the mutants, or with predicted structures generated by Molecular Dynamics (MD), predictions of redox potential changes are modeled using the Protein Dipoles Langevin Dipoles (PDLD) method. These calculations serve to model an electrostatic assessment of the redox potential change with simplified assumptions about heme iron chemistry, with the balance of the experimentally observed shifts in redox potential being thence attributed to changes in the ligand set and heme coordination chemistry, and/or other changes in the structure not directly evident in the X-ray structures (e.g. ionization states, specific roles played by solvent species, or conformationally flexible portions of the protein). Agreement between theory and experiment is good for all mutant proteins with the exception of the mutation Arg 48 to Ala, and His 52 to Ala. In the former case, the influence of phosphate buffer is adduced to account for the discrepancy, and measurements made in a bis-tris propane/2-(N-morpholino)ethanesulfonic acid buffer system agree well with theory. For the latter case, an unknown structural element relevant to His-52, and/or solvent influence in the mutant akin to anion binding in the distal pocket (though lacking proof that it is) manifests in this mutant. The use of exogenous (sixth) ligands in dissecting the contributions to control of redox potential are also explored as a pathway for model building.

Stereo-selectivity of human serum albumin to enantiomeric and isoelectronic pollutants dissected by spectroscopy, calorimetry and bioinformatics

1-naphthol (1N), 2-naphthol (2N) and 8-quinolinol (8H) are general water pollutants. 1N and 2N are the configurational enantiomers and 8H is isoelectronic to 1N and 2N. These pollutants when ingested are transported in the blood by proteins like human serum albumin (HSA). Binding of these pollutants to HSA has been explored to elucidate the specific selectivity of molecular recognition by this multiligand binding protein. The association constants (K(b)) of these pollutants to HSA were moderate (10(4)-10(5) M(-1)). The proximity of the ligands to HSA is also revealed by their average binding distance, r, which is estimated to be in the range of 4.39-5.37 nm. The binding free energy (ΔG) in each case remains effectively the same for each site because of enthalpy-entropy compensation (EEC). The difference observed between ΔC(p) (exp) and ΔC(p) (calc) are suggested to be caused by binding-induced flexibility changes in the HSA. Efforts are also made to elaborate the differences observed in binding isotherms obtained through multiple approaches of calorimetry, spectroscopy and bioinformatics. We suggest that difference in dissociation constants of pollutants by calorimetry, spectroscopic and computational approaches could correspond to occurrence of different set of populations of pollutants having different molecular characteristics in ground state and excited state. Furthermore, our observation of enhanced binding of pollutants (2N and 8H) in the presence of hemin signifies that ligands like hemin may enhance the storage period of these pollutants in blood that may even facilitate the ill effects of these pollutants.

Key factors for successful generation of protein-fragment structures requirement on protein, crystals, and technology

In the past two decades, fragment-based approaches have evolved as a predominant strategy in lead discovery. The availability of structural information on the interaction geometries of binding fragments is key to successful structure-guided fragment-to-lead evolution. In this chapter, we illustrate methodological advances for protein-fragment crystal structure generation in order to offer general lessons on the importance of fragment properties and the most appropriate crystallographic setup to evaluate them. We analyze elaborate protocols, methods, and clues applied to challenging complex formation projects. The results should assist medicinal chemists to select the most promising targets and strategies for fragment-based crystallography as well as provide a tutorial to structural biologists who attempt to determine protein-fragment structures.

Predictive power of molecular dynamics receptor structures in virtual screening

Molecular dynamics (MD) simulation is a well-established method for understanding protein dynamics. Conformations from unrestrained MD simulations have yet to be assessed for blind virtual screening (VS) by docking. This study presents a critical analysis of the predictive power of MD snapshots to this regard, evaluating two well-characterized systems of varying flexibility in ligand-bound and unbound configurations. Results from such VS predictions are discussed with respect to experimentally determined structures. In all cases, MD simulations provide snapshots that improve VS predictive power over known crystal structures, possibly due to sampling more relevant receptor conformations. Additionally, MD can move conformations previously not amenable to docking into the predictive range.

Rapid context-dependent ligand desolvation in molecular docking

In structure-based screens for new ligands, a molecular docking algorithm must rapidly score many molecules in multiple configurations, accounting for both the ligand's interactions with receptor and its competing interactions with solvent. Here we explore a context-dependent ligand desolvation scoring term for molecular docking. We relate the Generalized-Born effective Born radii for every ligand atom to a fractional desolvation and then use this fraction to scale an atom-by-atom decomposition of the full transfer free energy. The fractional desolvation is precomputed on a scoring grid by numerically integrating over the volume of receptor proximal to a ligand atom, weighted by distance. To test this method's performance, we dock ligands versus property-matched decoys over 40 DUD targets. Context-dependent desolvation better enriches ligands compared to both the raw full transfer free energy penalty and compared to ignoring desolvation altogether, though the improvement is modest. More compellingly, the new method improves docking performance across receptor types. Thus, whereas entirely ignoring desolvation works best for charged sites and overpenalizing with full desolvation works well for neutral sites, the physically more correct context-dependent ligand desolvation is competitive across both types of targets. The method also reliably discriminates ligands from highly charged molecules, where ignoring desolvation performs poorly. Since this context-dependent ligand desolvation may be precalculated, it improves docking reliability with minimal cost to calculation time and may be readily incorporated into any physics-based docking program.

Toward protein engineering for phytoremediation: possibilities and challenges

The combination of rational protein engineering and directed evolution techniques allow for the redesign of enzymes with tailored properties for use in environmental remediation. This review summarizes current molecular methods for either altering or improving protein function and highlights examples of how these methods can address bioremediation problems. Although much of the protein engineering applied to environmental clean-up employs microbial systems, there is great potential for and significant challenges to translating these approaches to plant systems for phytoremediation purposes. Protein engineering technologies combined with genomic information and metabolic engineering strategies hold promise for the design of plants and microbes to remediate organic and inorganic pollutants.

Quantum mechanics/molecular mechanics strategies for docking pose refinement: distinguishing between binders and decoys in cytochrome C peroxidase

We investigate the effect of systematically applying molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) to docked poses in an attempt to improve the correspondence between theoretical prediction and experimental observation. The proposed scheme involves running a short time scale MD simulation on a docked ligand pose (and any known structurally important crystal structure waters in the active site), followed by QM/MM minimization. Both of these steps are relatively fast for moderately sized ligands; longer time scale MD involving the protein is not found to improve the results. The final binding energy is given in terms of the QM/MM total energy, a van der Waals correction, and a term to account for desolvation effects. This methodology is first tested with a trypsin inhibitor, for which we establish the importance of running MD before reoptimizing with QM/MM. The method is then applied to cytochrome c peroxidase using a set of binders and decoys. In this example, the proposed methodology affords much better discrimination between binders and decoys than the traditional docking approach used. For both systems presented, application of this protocol results in a significantly better energetic ranking and a smaller root mean squared deviation from known crystallographic ligand poses. This work highlights the importance of including polarization effects through QM/MM and of sampling with MD to refine a set of initial docked poses.

Water in cavity-ligand recognition

We use explicit solvent molecular dynamics simulations to estimate free energy, enthalpy, and entropy changes along the cavity-ligand association coordinate for a set of seven model systems with varying physicochemical properties. Owing to the simplicity of the considered systems we can directly investigate the role of water thermodynamics in molecular recognition. A broad range of thermodynamic signatures is found in which water (rather than cavity or ligand) enthalpic or entropic contributions appear to drive cavity-ligand binding or rejection. The unprecedented, nanoscale picture of hydration thermodynamics can help the interpretation and design of protein-ligand binding experiments. Our study opens appealing perspectives to tackle the challenge of solvent entropy estimation in complex systems and for improving molecular simulation models.

Algorithm for rigid-body Brownian dynamics

We present an algorithm for performing rigid-body Brownian dynamics that can take into account the hydrodynamic properties (translational and rotational friction tensors and the coupling between them) of each rigid body. In the zero temperature limit, the error term scales as Delta;{4} for time step Delta , while at nonzero temperatures the error scaling is Delta;{5/2} . We test the algorithm by applying it to a molecule of four-aminopyridine in water. We intend to use the algorithm to model the interaction between biological ion channels and other channel blocker molecules, but it may also have applicability to modeling other small particles such as colloids or nanoparticles.

Evaluation of an inverse molecular design algorithm in a model binding site

Computational molecular design is a useful tool in modern drug discovery. Virtual screening is an approach that docks and then scores individual members of compound libraries. In contrast to this forward approach, inverse approaches construct compounds from fragments, such that the computed affinity, or a combination of relevant properties, is optimized. We have recently developed a new inverse approach to drug design based on the dead-end elimination and A* algorithms employing a physical potential function. This approach has been applied to combinatorially constructed libraries of small-molecule ligands to design high-affinity HIV-1 protease inhibitors (Altman et al., J Am Chem Soc 2008;130:6099-6013). Here we have evaluated the new method using the well-studied W191G mutant of cytochrome c peroxidase. This mutant possesses a charged binding pocket and has been used to evaluate other design approaches. The results show that overall the new inverse approach does an excellent job of separating binders from nonbinders. For a few individual cases, scoring inaccuracies led to false positives. The majority of these involve erroneous solvation energy estimation for charged amines, anilinium ions, and phenols, which has been observed previously for a variety of scoring algorithms. Interestingly, although inverse approaches are generally expected to identify some but not all binders in a library, due to limited conformational searching, these results show excellent coverage of the known binders while still showing strong discrimination of the nonbinders.

Replacement of the axial histidine heme ligand with cysteine in nitrophorin 1: spectroscopic and crystallographic characterization

To evaluate the potential of using heme-containing lipocalin nitrophorin 1 (NP1) as a template for protein engineering, we have replaced the native axial heme-coordinating histidine residue with glycine, alanine, and cysteine. We report here the characterization of the cysteine mutant H60C_NP1 by spectroscopic and crystallographic methods. The UV/vis, resonance Raman, and magnetic circular dichroism spectra suggest weak thiolate coordination of the ferric heme in the H60C_NP1 mutant. Reduction to the ferrous state resulted in loss of cysteine coordination, while addition of exogenous imidazole ligands gave coordination changes that varied with the ligand. Depending on the substitution of the imidazole, we could distinguish three heme coordination states: five-coordinate monoimidazole, six-coordinate bisimidazole, and six-coordinate imidazole/thiolate. Ligand binding affinities were measured and found to be generally 2-3 orders of magnitude lower for the H60C mutant relative to NP1. Two crystal structures of the H60C_NP1 in complex with imidazole and histamine were solved to 1.7- and 1.96-A resolution, respectively. Both structures show that the H60C mutation is well tolerated by the protein scaffold and suggest that heme-thiolate coordination in H60C_NP1 requires some movement of the heme within its binding cavity. This adjustment may be responsible for the ease with which the engineered heme-thiolate coordination can be displaced by exogenous ligands.

Replacement of an electron transfer pathway in cytochrome c peroxidase with a surrogate peptide

A proposed electron transfer pathway in cytochrome c peroxidase was previously excised from the structure by design. The engineered channel mutant was shown to bind peptide surrogates without restoration of cyt c oxidation. Here, we report the 1.6 A crystal structure of (N-benzimidazole-propionic acid)-Gly-Ala-Ala bound within the engineered channel. The peptide retains many features of the native electron transfer pathway: placement of benzimidazole at the position of the Trp-191 radical, hydrogen bonding to Asp235, and positioning of the C-terminus near the point where wild type CcP makes closest contact to cyt c. The inability of this surrogate pathway to restore function supports proposals that electron transfer requires the Trp-191 radical.

An improved relaxed complex scheme for receptor flexibility in computer-aided drug design

The interactions among associating (macro)molecules are dynamic, which adds to the complexity of molecular recognition. While ligand flexibility is well accounted for in computational drug design, the effective inclusion of receptor flexibility remains an important challenge. The relaxed complex scheme (RCS) is a promising computational methodology that combines the advantages of docking algorithms with dynamic structural information provided by molecular dynamics (MD) simulations, therefore explicitly accounting for the flexibility of both the receptor and the docked ligands. Here, we briefly review the RCS and discuss new extensions and improvements of this methodology in the context of ligand binding to two example targets: kinetoplastid RNA editing ligase 1 and the W191G cavity mutant of cytochrome c peroxidase. The RCS improvements include its extension to virtual screening, more rigorous characterization of local and global binding effects, and methods to improve its computational efficiency by reducing the receptor ensemble to a representative set of configurations. The choice of receptor ensemble, its influence on the predictive power of RCS, and the current limitations for an accurate treatment of the solvent contributions are also briefly discussed. Finally, we outline potential methodological improvements that we anticipate will assist future development.

Rescoring docking hit lists for model cavity sites: predictions and experimental testing

Molecular docking computationally screens thousands to millions of organic molecules against protein structures, looking for those with complementary fits. Many approximations are made, often resulting in low "hit rates." A strategy to overcome these approximations is to rescore top-ranked docked molecules using a better but slower method. One such is afforded by molecular mechanics-generalized Born surface area (MM-GBSA) techniques. These more physically realistic methods have improved models for solvation and electrostatic interactions and conformational change compared to most docking programs. To investigate MM-GBSA rescoring, we re-ranked docking hit lists in three small buried sites: a hydrophobic cavity that binds apolar ligands, a slightly polar cavity that binds aryl and hydrogen-bonding ligands, and an anionic cavity that binds cationic ligands. These sites are simple; consequently, incorrect predictions can be attributed to particular errors in the method, and many likely ligands may actually be tested. In retrospective calculations, MM-GBSA techniques with binding-site minimization better distinguished the known ligands for each cavity from the known decoys compared to the docking calculation alone. This encouraged us to test rescoring prospectively on molecules that ranked poorly by docking but that ranked well when rescored by MM-GBSA. A total of 33 molecules highly ranked by MM-GBSA for the three cavities were tested experimentally. Of these, 23 were observed to bind--these are docking false negatives rescued by rescoring. The 10 remaining molecules are true negatives by docking and false positives by MM-GBSA. X-ray crystal structures were determined for 21 of these 23 molecules. In many cases, the geometry prediction by MM-GBSA improved the initial docking pose and more closely resembled the crystallographic result; yet in several cases, the rescored geometry failed to capture large conformational changes in the protein. Intriguingly, rescoring not only rescued docking false positives, but also introduced several new false positives into the top-ranking molecules. We consider the origins of the successes and failures in MM-GBSA rescoring in these model cavity sites and the prospects for rescoring in biologically relevant targets.

Dynamics, Hydration, and Motional Averaging of a Loop-Gated Artificial Protein Cavity: The W191G Mutant of Cytochrome c Peroxidase in Water as Revealed by Molecular Dynamics Simulations

Five molecular dynamics simulations of the W191G cavity mutant of cytochrome c peroxidase in explicit water reveal distinct dynamic and hydration behavior depending on the closed or open state of the flexible loop gating the cavity, the binding of (K+ or small molecule) cations, and the system temperature. The conformational spaces sampled by the loop region and by the cavity significantly reduce upon binding. The largest ordering factor on water dynamics is the presence of the K+ ion occupying the gated cavity. Considerable water exchange occurs for the open-gate cavity when no ligand or cation is bound. In all cases, good correspondence is found between the calculated (ensemble-averaged) location of water molecules and the water sites determined by X-ray crystallography experiments. However, our simulations suggest that these sites do not necessarily correspond to the presence of bound water molecules. In fact, individual water molecules may repeatedly exchange within the cavity volume yet occupy on average these water sites. Four major conclusions emerge. First, it seems misleading to interpret the conformation of protein loop regions in terms of single dominant structures. Second, our simulations support the general picture of Pro 190 cis-trans isomerization as a determinant of the loop-opening mechanism. Third, receptor flexibility is fundamental for ligand binding and molecular recognition, and our results suggest its importance for the docking of small compounds to the artificial cavity. Fourth, after validation against the available experimental data, molecular dynamics simulations can be used to characterize the dynamics and exchange of water molecules and ions, providing atomic level and time-dependent information otherwise inaccessible to experiments.

Outliers in SAR and QSAR: is unusual binding mode a possible source of outliers?

A lead optimization is usually carried out by structure-activity relationship (SAR) and/or quantitative structure-activity relationship (QSAR) studies. One of the assumptions in SAR and QSAR studies is that similar analogs bind to the same binding site in a similar binding mode. One often observes that there are outliers, especially in QSAR. However, most QSAR studies are carried out focusing their attention to the development of QSAR and leave the outliers without much attention. We searched a number of ligand-bound X-ray crystal structures from the protein structure database to find evidences that could indicate a possible source of outliers in SAR or QSAR. Our results show that unusual binding mode could be a source of outliers.

Probing molecular docking in a charged model binding site

A model binding site was used to investigate charge-charge interactions in molecular docking. This simple site, a small (180A(3)) engineered cavity in cyctochrome c peroxidase (CCP), is negatively charged and completely buried from solvent, allowing us to explore the balance between electrostatic energy and ligand desolvation energy in a system where many of the common approximations in docking do not apply. A database with about 5300 molecules was docked into this cavity. Retrospective testing with known ligands and decoys showed that overall the balance between electrostatic interaction and desolvation energy was captured. More interesting were prospective docking scre"ens that looked for novel ligands, especially those that might reveal problems with the docking and energy methods. Based on screens of the 5300 compound database, both high-scoring and low-scoring molecules were acquired and tested for binding. Out of 16 new, high-scoring compounds tested, 15 were observed to bind. All of these were small heterocyclic cations. Binding constants were measured for a few of these, they ranged between 20microM and 60microM. Crystal structures were determined for ten of these ligands in complex with the protein. The observed ligand geometry corresponded closely to that predicted by docking. Several low-scoring alkyl amino cations were also tested and found to bind. The low docking score of these molecules owed to the relatively high charge density of the charged amino group and the corresponding high desolvation penalty. When the complex structures of those ligands were determined, a bound water molecule was observed interacting with the amino group and a backbone carbonyl group of the cavity. This water molecule mitigates the desolvation penalty and improves the interaction energy relative to that of the "naked" site used in the docking screen. Finally, six low-scoring neutral molecules were also tested, with a view to looking for false negative predictions. Whereas most of these did not bind, two did (phenol and 3-fluorocatechol). Crystal structures for these two ligands in complex with the cavity site suggest reasons for their binding. That these neutral molecules do, in fact bind, contradicts previous results in this site and, along with the alkyl amines, provides instructive false negatives that help identify weaknesses in our scoring functions. Several improvements of these are considered.

Grand canonical Monte Carlo simulation of ligand-protein binding

A new application of the grand canonical thermodynamics ensemble to compute ligand-protein binding is described. The described method is sufficiently rapid that it is practical to compute ligand-protein binding free energies for a large number of poses over the entire protein surface, thus identifying multiple putative ligand binding sites. In addition, the method computes binding free energies for a large number of poses. The method is demonstrated by the simulation of two protein-ligand systems, thermolysin and T4 lysozyme, for which there is extensive thermodynamic and crystallographic data for the binding of small, rigid ligands. These low-molecular-weight ligands correspond to the molecular fragments used in computational fragment-based drug design. The simulations correctly identified the experimental binding poses and rank ordered the affinities of ligands in each of these systems.

Defining the Structural Basis for Assembly of a Transmembrane Cytochrome

To define the structural basis for cofactor binding to membrane proteins, we introduce a manageable model system, which allows us, for the first time, to study the influence of individual transmembrane helices and of single amino acid residues on the assembly of a transmembrane cytochrome. In vivo as well as in vitro analyses indicate central roles of single amino acid residues for either interaction of the transmembrane helices or for binding of the cofactor. The results clearly show that interaction of the PsbF transmembrane helix is independent from binding of the heme cofactor. On the other hand, binding of the cofactor highly depends on helix-helix interactions. By site-directed mutagenesis critical amino acid residues were identified, which are involved in the assembly of a functional transmembrane cytochrome. Especially, a highly conserved glycine residue is critical for interaction of the transmembrane helices and assembly of the cytochrome. Based on the two-stage-model of alpha-helical membrane protein folding, the presented results clearly indicate a third stage of membrane protein folding, in which a cofactor binds to a pre-assembled transmembrane protein.

Automated docking of ligands to an artificial active site: augmenting crystallographic analysis with computer modeling

The W191G cavity of cytochrome c peroxidase is useful as a model system for introducing small molecule oxidation in an artificially created cavity. A set of small, cyclic, organic cations was previously shown to bind in the buried, solvent-filled pocket created by the W191G mutation. We docked these ligands and a set of non-binders in the W191G cavity using AutoDock 3.0. For the ligands, we compared docking predictions with experimentally determined binding energies and X-ray crystal structure complexes. For the ligands, predicted binding energies differed from measured values by +/- 0.8 kcal/mol. For most ligands, the docking simulation clearly predicted a single binding mode that matched the crystallographic binding mode within 1.0 A RMSD. For 2 ligands, where the docking procedure yielded an ambiguous result, solutions matching the crystallographic result could be obtained by including an additional crystallographically observed water molecule in the protein model. For the remaining 2 ligands, docking indicated multiple binding modes, consistent with the original electron density, suggesting disordered binding of these ligands. Visual inspection of the atomic affinity grid maps used in docking calculations revealed two patches of high affinity for hydrogen bond donating groups. Multiple solutions are predicted as these two sites compete for polar hydrogens in the ligand during the docking simulation. Ligands could be distinguished, to some extent, from non-binders using a combination of two trends: predicted binding energy and level of clustering. In summary, AutoDock 3.0 appears to be useful in predicting key structural and energetic features of ligand binding in the W191G cavity.

Trapping of peptide-based surrogates in an artificially created channel of cytochrome c peroxidase

As recently described, the deliberate removal of the proposed electron transfer pathway from cytochrome c peroxidase resulted in the formation of an extended ligand-binding channel. The engineered channel formed a template for the removed peptide segment, suggesting that synthetic surrogates might be introduced to replace the native electron transfer pathway. This approach could be united with the recent development of sensitizer-linked substrates to initiate and study electron transfer, allowing access to unresolved issues about redox mechanism of the enzyme. Here, we present the design, synthesis, and screening of a peptide library containing natural and unnatural amino acids to identify the structural determinants for binding this channel mutant. Only one peptide, (benzimidazole-propionic acid)-Gly-Ala-Ala, appeared to interact, and gave evidence for both reversible and kinetically trapped binding, suggesting multiple conformations for the channel protein. Notably, this peptide was the most analogous to the removed electron transfer sequence, supporting the use of a cavity-template strategy for design of specific sensitizer-linked substrates as replacements for the native electron transfer pathway.

Excision of a proposed electron transfer pathway in cytochrome c peroxidase and its replacement by a ligand-binding channel

A previously proposed electron transfer (ET) pathway in the heme enzyme cytochrome c peroxidase has been excised from the structure, leaving an open ligand-binding channel in its place. Earlier studies on cavity mutants of this enzyme have revealed structural plasticity in this region of the molecule. Analysis of these structures has allowed the design of a variant in which the specific section of protein backbone representing a previously proposed ET pathway is accurately extracted from the protein. A crystal structure verified the creation of an open channel that overlays the removed segment, extending from the surface of the protein to the heme at the core of the protein. A number of heterocyclic cations were found to bind to the proximal-channel mutant with affinities that can be rationalized based on the structures. It is proposed that small ligands bind more weakly to the proximal-channel mutant than to the W191G cavity due to an increased off rate of the open channel, whereas larger ligands are able to bind to the channel mutant without inducing large conformational changes. The structure of benzimidazole bound to the proximal-channel mutant shows that the ligand accurately overlays the position of the tryptophan radical center that was removed from the wild-type enzyme and displaces four of the eight ordered solvent molecules seen in the empty cavity. Ligand binding also caused a small rearrangement of the redesigned protein loop, perhaps as a result of improved electrostatic interactions with the ligand. The engineered channel offers the potential for introducing synthetic replacements for the removed structure, such as sensitizer-linked substrates. These installed "molecular wires" could be used to rapidly initiate reactions, trap reactive intermediates, or answer unresolved questions about ET pathways.