Evaluation of ab Initio Charge Determination Methods for Use in Continuum Solvation Calculations

Efficient calculation of molecular configurational entropies using an information theoretic approximation

Accurate computation of free energy changes upon molecular binding remains a challenging problem, and changes in configurational entropy are especially difficult due to the potentially large numbers of local minima, anharmonicity, and high-order coupling among degrees of freedom. Here we propose a new method to compute molecular entropies based on the maximum information spanning tree (MIST) approximation that we have previously developed. Estimates of high-order couplings using only low-order terms provide excellent convergence properties, and the theory is also guaranteed to bound the entropy. The theory is presented together with applications to the calculation of the entropies of a variety of small molecules and the binding entropy change for a series of HIV protease inhibitors. The MIST framework developed here is demonstrated to compare favorably with results computed using the related mutual information expansion (MIE) approach, and an analysis of similarities between the methods is presented.

Systematic placement of structural water molecules for improved scoring of protein-ligand interactions

Structural water molecules are found in many protein-ligand complexes. They are known to be vital in mediating hydrogen-bonding interactions and, in some cases, key for facilitating tight binding. It is thus very important to consider water molecules when attempting to model protein-ligand interactions for cognate ligand identification, virtual screening and drug design. While the rigid treatment of water molecules present in structures is feasible, the more relevant task of treating all possible positions and orientations of water molecules with each possible ligand pose is computationally daunting. Current methods in molecular docking provide partial treatment for such water molecules, with modest success. Here we describe a new method employing dead-end elimination to place water molecules within a binding site, bridging interactions between protein and ligand. Dead-end elimination permits a thorough, though still incomplete, treatment of water placement. The results show that this method is able to place water molecules correctly within known complexes and to create physically reasonable hydrogen bonds. The approach has also been incorporated within an inverse molecular design approach, to model a variety of compounds in the process of de novo ligand design. The inclusion of structural water molecules, combined with ranking based on the electrostatic contribution to binding affinity, improves a number of otherwise poor energetic predictions.

APBSmem: a graphical interface for electrostatic calculations at the membrane

Electrostatic forces are one of the primary determinants of molecular interactions. They help guide the folding of proteins, increase the binding of one protein to another and facilitate protein-DNA and protein-ligand binding. A popular method for computing the electrostatic properties of biological systems is to numerically solve the Poisson-Boltzmann (PB) equation, and there are several easy-to-use software packages available that solve the PB equation for soluble proteins. Here we present a freely available program, called APBSmem, for carrying out these calculations in the presence of a membrane. The Adaptive Poisson-Boltzmann Solver (APBS) is used as a back-end for solving the PB equation, and a Java-based graphical user interface (GUI) coordinates a set of routines that introduce the influence of the membrane, determine its placement relative to the protein, and set the membrane potential. The software Jmol is embedded in the GUI to visualize the protein inserted in the membrane before the calculation and the electrostatic potential after completing the computation. We expect that the ease with which the GUI allows one to carry out these calculations will make this software a useful resource for experimenters and computational researchers alike. Three examples of membrane protein electrostatic calculations are carried out to illustrate how to use APBSmem and to highlight the different quantities of interest that can be calculated.

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.

HIV-1 protease inhibitors from inverse design in the substrate envelope exhibit subnanomolar binding to drug-resistant variants

The acquisition of drug-resistant mutations by infectious pathogens remains a pressing health concern, and the development of strategies to combat this threat is a priority. Here we have applied a general strategy, inverse design using the substrate envelope, to develop inhibitors of HIV-1 protease. Structure-based computation was used to design inhibitors predicted to stay within a consensus substrate volume in the binding site. Two rounds of design, synthesis, experimental testing, and structural analysis were carried out, resulting in a total of 51 compounds. Improvements in design methodology led to a roughly 1000-fold affinity enhancement to a wild-type protease for the best binders, from a Ki of 30-50 nM in round one to below 100 pM in round two. Crystal structures of a subset of complexes revealed a binding mode similar to each design that respected the substrate envelope in nearly all cases. All four best binders from round one exhibited broad specificity against a clinically relevant panel of drug-resistant HIV-1 protease variants, losing no more than 6-13-fold affinity relative to wild type. Testing a subset of second-round compounds against the panel of resistant variants revealed three classes of inhibitors: robust binders (maximum affinity loss of 14-16-fold), moderate binders (35-80-fold), and susceptible binders (greater than 100-fold). Although for especially high-affinity inhibitors additional factors may also be important, overall, these results suggest that designing inhibitors using the substrate envelope may be a useful strategy in the development of therapeutics with low susceptibility to resistance.

Defining Cdk5 ligand chemical space with small molecule inhibitors of tau phosphorylation

Cyclin-dependent kinase 5 (Cdk5) is widely viewed as a possible target for a wide variety of neurological disorders. One pathological role attributed to Cdk5 is the abnormal phosphorylation of tau that may lead to the neuronal inclusions known as neurofibrillary tangles. A high through-put screen for inhibitors of Cdk5-mediated phosphorylation of tau resulted in three compounds with distinct mechanisms of action. One compound is competitive with ATP and has a high affinity for the Cdk5 ATP binding pocket. The second compound also competes with ATP, is noncompetitive with tau, and (uniquely among this class of inhibitors) displaces adjacent amino acid residues to make room for the nitrophenyl group. A third compound did not compete with ATP, but did compete with tau at low concentrations of tau. The SAR and charge optimization derived from cocrystals of the two ATP competitors along with cocrystals of three other ATP competitors map out the importance of filling and properly charging different regions of the ATP binding pocket. Taken together, this analysis shows how the structure of Cdk5 constrains the space of potential inhibitors and reveals a pocket unfilled in all of the structures. These leads could be a starting point for structure-based drug design of more potent and selective inhibitors.

A Computational Study of Oxiranyllithium

[reaction: see text] Computational methods were used to determine the structure, bonding, and aggregation states of oxiranyllithium in the gas phase and in THF solution, at 200 and 298 K. THF solvation was modeled by microsolvation with explicit THF ligands, forming a supermolecule that includes the oxiranyllithium aggregate and its first solvation shell. Because oxiranyllithium has a chiral center, two diastereomeric dimers were formed, the RR and the RS, along with their enantiomers. Similarly, three diastereomers of the tetramer were formed, the RRRR, RRRS, and RRSS and their enantiomers. Oxiranyllithium was found to exist predominantly as the tetramer in the gas phase, while the dimer was the dominant species in THF solution. The relative concentrations of the different stereoisomers were calculated from equilibrium constants.

Design of improved protein inhibitors of HIV-1 cell entry: Optimization of electrostatic interactions at the binding interface

Continuum electrostatic methods are a powerful tool for the analysis and design of biomolecular complexes, with methodologies that allow for the detailed analysis of the electrostatic contributions to binding affinities and procedures for computing the properties of electrostatically optimal ligands. We have applied these methods to the design of improved inhibitors of HIV-1 cell entry. HIV infection of a cell requires viral-cell membrane fusion, an event partially driven by a large-scale conformational change in the viral membrane glycoprotein gp41. This transformation involves the docking of a helix from the C-terminal region of three gp41 chains against a pre-formed trimeric-coiled coil; several protein constructs that inhibit membrane fusion act by binding to an isolated C-terminal helix and blocking the formation of the fusogenic structure. A detailed analysis of the electrostatic contributions to the binding of one such inhibitor (5-Helix) to a C-terminal helix was performed using the X-ray crystal structure of the core of the HIV-1 gp41 ectodomain as a structural model, and several residues on 5-Helix that make substantial contributions to binding, both favorable and unfavorable, were identified. An electrostatic affinity optimization methodology was applied to the side chains of 5-Helix, with the results showing that significant improvements in binding affinity are possible if the electrostatic contributions to the binding free energy are optimized. Several mutations accessible by experimental methods are suggested, with calculated improvements in binding affinity of as much as 500-fold and greater.

The Mechanism of Proton Transfer between Adjacent Sites Exposed to Water

The surface of a protein, or a membrane, is spotted with a multitude of proton-binding sites, some of which are only a few angstroms apart. When a proton is released from one site, it propagates through the water by a random walk under the bias of the local electrostatic potential determined by the distribution of the charges on the protein. Some of the released protons disperse into the bulk, but during the first few nanoseconds, the released protons can be trapped by encounter with nearby acceptor sites. This process resembles a scenario which corresponds with the time-dependent Debye-Smoluchowski equation. In the present study, we investigated the mechanism of proton transfer between sites that are only a few angstroms apart, using as a model the proton exchange between sites on a small molecule, fluorescein, having two, spectrally distinguishable, proton-binding sites. The first site is the oxyanion on the chromophore ring structure. The second site is the carboxylate moiety on the benzene ring of the molecule. Through our experiments, we were able to reconstruct the state of protonation at each site and the velocity of proton transfer between them. The fluorescein was protonated by a few nanosecond long proton pulse under specific conditions that ensured that the dye molecules would be protonated only by a single proton. The dynamics of the protonation of the chromophore were measured under varying initial conditions (temperature, ionic strength, and different solvents (H(2)O or D(2)O)), and the velocity of the proton transfer between the two sites was extracted from the overall global analysis of the signals. The dynamics of the proton transfer between the two proton-binding sites of the fluorescein indicated that the efficiency of the site-to-site proton transfer is very sensitive to the presence of the screening electrolyte and has a very high kinetic isotope effect (KIE = 55). These two parameters clearly distinguish the mechanism from proton diffusion in bulk water. The activation energy of the reaction (E(a) = 11 kcal mol(-1)) is also significantly higher than the activation energy for proton dissociation in bulk water (E(a) approximately 2.5 kcal mol(-1)). These observations are discussed with respect to the effect of the solute on the water molecules located within the solvation layer.

The Electrostatic Origin of Abraham's Solute Polarity Parameter

A computational method was developed which relates the empirical linear solvation energy relationship (LSER) solute polarity parameter, S (formerly denoted ), to two more fundamental quantities: a polarizability term and a computed solvent-accessible-surface electrostatic term. Electrostatics computations were conducted explicitly or with dielectric field polarizable continuum models (PCM, SCIPCM, IPCM), employing a density functional theory (B3LYP/6-311G(2df,2p)) or efficient Hartree-Fock (HF/MIDI!) method for 90 polar and nonpolar organic solutes. Electrostatic parameters calculated at electron isodensity solute surfaces were found to produce significantly better correlations with empirical S values than the same electrostatic parameters deduced from a fixed Bondi atomic radii based surface. The best-fit expression was found employing SCIPCM/IPCM at the 0.0004 e(-)/bohr(3) solvent-accessible-surface: S(fit)() = 0.46E - 0.091SigmaV(s)()(2), with squared correlation coefficient = 0.96 and standard deviation = 0.10, where E is a measured solute excess polarizability scale and SigmaV(s)()(2) is a quantum-calculated solute electrostatic descriptor in kcal A/mol. The resulting model is more accurate than previously developed estimation approaches and relies on only two fitted coefficients; it has the potential advantage of applicability to any solute composed of C, H, N, O, S, F, Cl, and Br. Finally, this investigation offers quantitative insight into the relative contributions of solute polarity and solute polarizability to the empirical LSER polarity parameter, S.

Escherichia coli glutaminyl-tRNA synthetase is electrostatically optimized for binding of its cognate substrates

Natural evolution has resulted in protein molecules displaying a wide range of binding properties that include extremes of affinity and specificity. A detailed understanding of the principles underlying protein structure-function relationships, particularly with respect to binding properties, would greatly enhance molecular engineering and ligand design studies. Here, we have analyzed the interactions of an aminoacyl-tRNA synthetase for which strong evolutionary pressure has enforced high specificity for substrate binding and catalysis. Electrostatic interactions have been identified as one efficient mechanism for enhancing binding specificity; as such, the effects of charged and polar groups were the focus of this study. The binding of glutaminyl-tRNA synthetase from Escherichia coli to several ligands, including the natural substrates, was analyzed. The electrostatic complementarity of the enzyme to its ligands was assessed using measures derived from affinity optimization theory. The results were independent of the details of the calculational parameters, including the value used for the protein dielectric constant. Glutamine and ATP, two of the natural ligands, were found to be extremely complementary to their binding sites, particularly in regions seen to make electrostatic interactions in the structure. These data suggest that the optimization of electrostatic interactions has played an important role in guiding the evolution of this enzyme. The results also show that the enzyme is able to effectively select for high affinity and specificity for the same chemical moieties both in the context of smaller substrates, and in that of a larger reactive intermediate. The regions of greatest non-complementarity between the enzyme and ligands are the portions of the ligand that make few polar contacts with the binding site, as well as the sites of chemical reaction, where overly strong electrostatic binding interactions with the substrate could hinder catalysis. The results also suggest that the negative charge on the phosphorus center of glutaminyl-adenylate plays an important role in the tight binding of this intermediate, and thus that adenylate analogs that preserve the negative charge in this region may bind substantially tighter than analogs where this group is replaced with a neutral group, such as the sulfamoyl family, which can make similar hydrogen bonds but is uncharged.

Substantial Energetic Improvement with Minimal Structural Perturbation in a High Affinity Mutant Antibody

Here, we compare an antibody with the highest known engineered affinity (K(d)=270 fM) to its high affinity wild-type (K(d)=700 pM) through thermodynamic, kinetic, structural, and theoretical analyses. The 4M5.3 anti-fluorescein single chain antibody fragment (scFv) contains 14 mutations from the wild-type 4-4-20 scFv and has a 1800-fold increase in fluorescein-binding affinity. The dissociation rate is approximately 16,000 times slower in the mutant; however, this substantial improvement is offset somewhat by the association rate, which is ninefold slower in the mutant. Enthalpic contributions to binding were found by calorimetry to predominate in the differential binding free energy. The crystal structure of the 4M5.3 mutant complexed with antigen was solved to 1.5A resolution and compared with a previously solved structure of an antigen-bound 4-4-20 Fab fragment. Strikingly, the structural comparison shows little difference between the two scFv molecules (backbone RMSD of 0.6A), despite the large difference in affinity. Shape complementarity exhibits a small improvement between the variable light chain and variable heavy chain domains within the antibody, but no significant improvement in shape complementarity of the antibody with the antigen is observed in the mutant over the wild-type. Theoretical modeling calculations show electrostatic contributions to binding account for -1.2 kcal/mol to -3.5 kcal/mol of the binding free energy change, of which -1.1 kcal/mol is directly associated with the mutated residue side-chains. The electrostatic analysis reveals several mechanistic explanations for a portion of the improvement. Collectively, these data provide an example where very high binding affinity is achieved through the cumulative effect of many small structural alterations.

A new interpretation of chlorine leaving group kinetic isotope effects; a theoretical approach

The chlorine leaving group kinetic isotope effects (KIEs) for the S(N)2 reactions between methyl chloride and a wide range of anionic, neutral, and radical anion nucleophiles were calculated in the gas phase and, in several cases, using a continuum solvent model. In contrast to the expected linear dependence of the chlorine KIEs on the C(alpha)-Cl bond order in the transition state, the KIEs fell in a very small range (1.0056-1.0091), even though the C(alpha)-Cl transition state bond orders varied widely from approximately 0.32 to 0.78, a range from reactant-like to very product-like. This renders chlorine KIEs, and possibly other leaving-group KIEs, less useful for studies of reaction mechanisms than commonly assumed. A partial explanation for this unexpected relationship between the C(alpha)-Cl transition state bond order and the magnitude of the chlorine KIE is presented.