Computational optimization of the SARS-CoV-2 receptor-binding-motif affinity for human ACE2.

The coronavirus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is responsible for the coronavirus disease 2019 pandemic, and the closely related SARS-CoV coronavirus enter cells by binding at the human angiotensin converting enzyme 2 (hACE2). The stronger hACE2 affinity of SARS-CoV-2 has been connected with its higher infectivity. In this work, we study hACE2 complexes with the receptor-binding domains (RBDs) of the human SARS-CoV-2 and human SARS-CoV viruses, using all-atom molecular dynamics simulations and computational protein design with a physics-based energy function. The molecular dynamics simulations identify charge-modifying substitutions between the CoV-2 and CoV RBDs, which either increase or decrease the hACE2 affinity of the SARS-CoV-2 RBD. The combined effect of these mutations is small, and the relative affinity is mainly determined by substitutions at residues in contact with hACE2. Many of these findings are in line and interpret recent experiments. Our computational protein design calculations redesign positions 455, 493, 494, and 501 of the SARS-CoV-2 receptor binding motif, which contact hACE2 in the complex and are important for ACE2 recognition. Sampling is enhanced by an adaptive importance sampling Monte Carlo method. Sequences with increased affinity replace CoV-2 glutamine by a negative residue at position 493; serine by a nonpolar or aromatic residue or an asparagine at position 494; and asparagine by valine or threonine at position 501. Substitutions at positions 455 and 501 have a smaller effect on affinity. Substitutions suggested by our design are seen in viral sequences encountered in other species, including bat and pangolin. Our results might be used to identify potential virus strains with higher human infectivity and assist in the design of peptide-based or peptidomimetic compounds with the potential to inhibit SARS-CoV-2 binding at hACE2.

Recognition of LD motifs by the focal adhesion targeting domains of focal adhesion kinase and proline-rich tyrosine kinase 2-beta: Insights from molecular dynamics simulations.

The focal adhesion kinase (FAK) and the proline-rich tyrosine kinase 2-beta (PYK2) are implicated in cancer progression and metastasis and represent promising biomarkers and targets for cancer therapy. FAK and PYK2 are recruited to focal adhesions (FAs) via interactions between their FA targeting (FAT) domains and conserved segments (LD motifs) on the proteins Paxillin, Leupaxin, and Hic-5. A promising new approach for the inhibition of FAK and PYK2 targets interactions of the FAK domains with proteins that promote localization at FAs. Advances toward this goal include the development of surface plasmon resonance, heteronuclear single quantum coherence nuclear magnetic resonance (HSQC-NMR) and fluorescence polarization assays for the identification of fragments or compounds interfering with the FAK-Paxillin interaction. We have recently validated this strategy, showing that Paxillin mimicking polypeptides with 2 to 3 LD motifs displace FAK from FAs and block kinase-dependent and independent functions of FAK, including downstream integrin signaling and FA localization of the protein p130Cas. In the present work we study by all-atom molecular dynamics simulations the recognition of peptides with the Paxillin and Leupaxin LD motifs by the FAK-FAT and PYK2-FAT domains. Our simulations and free-energy analysis interpret experimental data on binding of Paxillin and Leupaxin LD motifs at FAK-FAT and PYK2-FAT binding sites, and assess the roles of consensus LD regions and flanking residues. Our results can assist in the design of effective inhibitory peptides of the FAK-FAT: Paxillin and PYK2-FAT:Leupaxin complexes and the construction of pharmacophore models for the discovery of potential small-molecule inhibitors of the FAK-FAT and PYK2-FAT focal adhesion based functions.

Physics-Based Computational Protein Design: An Update.

We describe methods for physics-based protein design and some recent applications from our work. We present the physical interpretation of a MC simulation in sequence space and show that sequences and conformations form a well-defined statistical ensemble, explored with Monte Carlo and Boltzmann sampling. The folded state energy combines molecular mechanics for solutes with continuum electrostatics for solvent. We usually assume one or a few fixed protein backbone structures and discrete side chain rotamers. Methods based on molecular dynamics, which introduce additional backbone and side chain flexibility, are under development. The redesign of a PDZ domain and an aminoacyl-tRNA synthetase enzyme were successful. We describe a versatile, adaptive, Wang-Landau MC method that can be used to design for substrate affinity, catalytic rate, catalytic efficiency, or the specificity of these properties. The methods are transferable to all biomolecules, can be systematically improved, and give physical insights.

Hybrid MC/MD for protein design.

Computational protein design relies on simulations of a protein structure, where selected amino acids can mutate randomly, and mutations are selected to enhance a target property, such as stability. Often, the protein backbone is held fixed and its degrees of freedom are modeled implicitly to reduce the complexity of the conformational space. We present a hybrid method where short molecular dynamics (MD) segments are used to explore conformations and alternate with Monte Carlo (MC) moves that apply mutations to side chains. The backbone is fully flexible during MD. As a test, we computed side chain acid/base constants or pKa's in five proteins. This problem can be considered a special case of protein design, with protonation/deprotonation playing the role of mutations. The solvent was modeled as a dielectric continuum. Due to cost, in each protein we allowed just one side chain position to change its protonation state and the other position to change its type or mutate. The pKa's were computed with a standard method that scans a range of pH values and with a new method that uses adaptive landscape flattening (ALF) to sample all protonation states in a single simulation. The hybrid method gave notably better accuracy than standard, fixed-backbone MC. ALF decreased the computational cost a factor of 13.

Simple models for nonpolar solvation: Parameterization and testing.

Implicit solvent models are important for many biomolecular simulations. The polarity of aqueous solvent is essential and qualitatively captured by continuum electrostatics methods like Generalized Born (GB). However, GB does not account for the solvent-induced interactions between exposed hydrophobic sidechains or solute-solvent dispersion interactions. These "nonpolar" effects are often modeled through surface area (SA) energy terms, which lack realism, create mathematical singularities, and have a many-body character. We have explored an alternate, Lazaridis-Karplus (LK) gaussian energy density for nonpolar effects and a dispersion (DI) energy term proposed earlier, associated with GB electrostatics. We parameterized several combinations of GB, SA, LK, and DI energy terms, to reproduce 62 small molecule solvation free energies, 387 protein stability changes due to point mutations, and the structures of 8 protein loops. With optimized parameters, the models all gave similar results, with GBLK and GBDILK giving no performance loss compared to GBSA, and mean errors of 1.7 kcal/mol for the stability changes and 2 Å deviations for the loop conformations. The optimized GBLK model gave poor results in MD of the Trpcage mini-protein, but parameters optimized specifically for MD performed well for Trpcage and three other small proteins. Overall, the LK and DI nonpolar terms are valid alternatives to SA treatments for a range of applications. © 2017 Wiley Periodicals, Inc.

Structure based inhibitor design targeting glycogen phosphorylase B. Virtual screening, synthesis, biochemical and biological assessment of novel N-acyl-ß-d-glucopyranosylamines.

Glycogen phosphorylase (GP) is a validated target for the development of new type 2 diabetes treatments. Exploiting the Zinc docking database, we report the in silico screening of 1888 N-acyl-ß-d-glucopyranosylamines putative GP inhibitors differing only in their R groups. CombiGlide and GOLD docking programs with different scoring functions were employed with the best performing methods combined in a 'consensus scoring' approach to ranking of ligand binding affinities for the active site. Six selected candidates from the screening were then synthesized and their inhibitory potency was assessed both in vitro and ex vivo. Their inhibition constants' values, in vitro, ranged from 5 to 377µM while two of them were effective at causing inactivation of GP in rat hepatocytes at low µM concentrations. The crystal structures of GP in complex with the inhibitors were defined and provided the structural basis for their inhibitory potency and data for further structure based design of more potent inhibitors.

Computational protein design: the Proteus software and selected applications.

We describe an automated procedure for protein design, implemented in a flexible software package, called Proteus. System setup and calculation of an energy matrix are done with the XPLOR modeling program and its sophisticated command language, supporting several force fields and solvent models. A second program provides algorithms to search sequence space. It allows a decomposition of the system into groups, which can be combined in different ways in the energy function, for both positive and negative design. The whole procedure can be controlled by editing 2-4 scripts. Two applications consider the tyrosyl-tRNA synthetase enzyme and its successful redesign to bind both O-methyl-tyrosine and D-tyrosine. For the latter, we present Monte Carlo simulations where the D-tyrosine concentration is gradually increased, displacing L-tyrosine from the binding pocket and yielding the binding free energy difference, in good agreement with experiment. Complete redesign of the Crk SH3 domain is presented. The top 10000 sequences are all assigned to the correct fold by the SUPERFAMILY library of Hidden Markov Models. Finally, we report the acid/base behavior of the SNase protein. Sidechain protonation is treated as a form of mutation; it is then straightforward to perform constant-pH Monte Carlo simulations, which yield good agreement with experiment. Overall, the software can be used for a wide range of application, producing not only native-like sequences but also thermodynamic properties with errors that appear comparable to other current software packages.

Novel compstatin family peptides inhibit complement activation by drusen-like deposits in human retinal pigmented epithelial cell cultures.

We have used a novel human retinal pigmented epithelial (RPE) cell-based model that mimics drusen biogenesis and the pathobiology of age-related macular degeneration to evaluate the efficacy of newly designed peptide inhibitors of the complement system. The peptides belong to the compstatin family and, compared to existing compstatin analogs, have been optimized to promote binding to their target, complement protein C3, and to enhance solubility by improving their polarity/hydrophobicity ratios. Based on analysis of molecular dynamics simulation data of peptide-C3 complexes, novel binding features were designed by introducing intermolecular salt bridge-forming arginines at the N-terminus and at position -1 of N-terminal dipeptide extensions. Our study demonstrates that the RPE cell assay has discriminatory capability for measuring the efficacy and potency of inhibitory peptides in a macular disease environment.

A computational protein design protocol for optimization of the SARS-CoV-2 receptor-binding-motif affinity for human ACE2.

The present protocol describes the computational design of the SARS-CoV-2 receptor binding motif (RBD) to identify mutations that can potentially improve binding affinity for the human ACE2 (hACE2) receptor. We focus on four positions located at the interface with the hACE2 receptor in the RBD:hACE2 complex. We conduct the design with a high-throughput computational protein design (CPD) program, Proteus, incorporating an adaptive Monte Carlo (MC) protocol that promotes the selection of sequences with good binding affinities. For complete details on the use and execution of this protocol, please refer to Polydorides and Archontis (2021).

Computational Design of Peptides with Improved Recognition of the Focal Adhesion Kinase FAT Domain.

We describe a two-stage computational protein design (CPD) methodology for the design of peptides binding to the FAT domain of the protein focal adhesion kinase. The first stage involves high-throughput CPD calculations with the Proteus software. The energies of the folded state are described by a physics-based energy function and of the unfolded peptides by a knowledge-based model that reproduces aminoacid compositions consistent with a helicity scale. The obtained sequences are filtered in terms of the affinity and the stability of the complex. In the second stage, design sequences are further evaluated by all-atom molecular dynamics simulations and binding free energy calculations with a molecular mechanics/implicit solvent free energy function.

Computational protein design with a generalized Born solvent model: application to Asparaginyl-tRNA synthetase.

Computational Protein Design (CPD) is a promising method for high throughput protein and ligand mutagenesis. Recently, we developed a CPD method that used a polar-hydrogen energy function for protein interactions and a Coulomb/Accessible Surface Area (CASA) model for solvent effects. We applied this method to engineer aspartyl-adenylate (AspAMP) specificity into Asparaginyl-tRNA synthetase (AsnRS), whose substrate is asparaginyl-adenylate (AsnAMP). Here, we implement a more accurate function, with an all-atom energy for protein interactions and a residue-pairwise generalized Born model for solvent effects. As a first test, we compute aminoacid affinities for several point mutants of Aspartyl-tRNA synthetase (AspRS) and Tyrosyl-tRNA synthetase and stability changes for three helical peptides and compare with experiment. As a second test, we readdress the problem of AsnRS aminoacid engineering. We compare three design criteria, which optimize the folding free-energy, the absolute AspAMP affinity, and the relative (AspAMP-AsnAMP) affinity. The sequences and conformations are improved with respect to our previous, polar-hydrogen/CASA study: For several designed complexes, the AspAMP carboxylate forms three interactions with a conserved arginine and a designed lysine, as in the active site of the AspRS:AspAMP complex. The conformations and interactions are well maintained in molecular dynamics simulations and the sequences have an inverted specificity, favoring AspAMP over AsnAMP. The method is not fully successful, since experimental measurements with the seven most promising sequences show that they do not catalyze at a detectable level the adenylation of Asp (or Asn) with ATP. This may be due to weak AspAMP binding and/or disruption of transition-state stabilization.

Design of a modified mouse protein with ligand binding properties of its human analog by molecular dynamics simulations: the case of C3 inhibition by compstatin.

The peptide compstatin and its derivatives inhibit the complement-component protein C3 in primate mammals and are potential therapeutic agents against the unregulated activation of complement in humans, but are inactive against C3 from lower mammals. Recent molecular dynamics (MD) simulations showed that the most potent compstatin analog comprised entirely of natural amino acids (W4A9) had a smaller affinity for rat C3, due to reproducible changes in the rat protein structure with respect to the human protein, which eliminated or weakened specific protein-ligand interactions seen in the human C3:W4A9 complex. Here, we study by MD simulations three W4A9 complexes with the mouse C3 protein, and two "transgenic" mouse derivatives, containing a small number (6-9) of human C3 substitutions. The mouse complex experiences the conformational changes and affinity reduction of the rat complex. In the "transgenic" complexes, the conformation remains closer to that of the human complex, the protein-ligand interactions are improved, and the affinity for compstatin becomes "human-like." The present work creates new avenues for a compstatin-sensitive animal model. A similar strategy, involving the comparison of a series of complexes by MD simulations, could be used to design "transgenic" sequences in other systems.

Kinetics, in silico docking, molecular dynamics, and MM-GBSA binding studies on prototype indirubins, KT5720, and staurosporine as phosphorylase kinase ATP-binding site inhibitors: the role of water molecules examined.

With an aim toward glycogenolysis control in Type 2 diabetes, we have investigated via kinetic experiments and computation the potential of indirubin (IC50 > 50 µM), indirubin-3'-oxime (IC50 = 144 nM), KT5720 (K(i) = 18.4 nM) and staurosporine (K(i) = 0.37 nM) as phosphorylase kinase (PhKγtrnc) ATP-binding site inhibitors, with the latter two revealed as potent inhibitors in the low nM range. Because of lack of structural information, we have exploited information from homologous kinase complexes to direct in silico calculations (docking, molecular dynamics, and MMGBSA) to predict the binding characteristics of the four ligands. All inhibitors are predicted to bind in the same active site area as the ATP adenine ring, with binding dominated by hinge region hydrogen bonds to Asp104:O and Met106:O (all four ligands) and also Met106:NH (for the indirubins). The PhKγtrnc-staurosporine complex has the greatest number of receptor-ligand hydrogen bonds, while for the indirubin-3'-oxime and KT5720 complexes there is an important network of interchanging water molecules bridging inhibitor-enzyme contacts. The MM-GBSA results revealed the source of staurosporine's low nM potency to be favorable electrostatic interactions, while KT5720 has strong van der Waals contributions. KT5720 interacts with the greatest number of protein residues either by direct or 1-water bridged hydrogen bond interactions, and the potential for more selective PhK inhibition based on a KT5720 analogue has been established. Including receptor flexibility in Schrödinger induced-fit docking calculations in most cases correctly predicted the binding modes as compared with the molecular dynamics structures; the algorithm was less effective when there were key structural waters bridging receptor-ligand contacts.

Species specificity of the complement inhibitor compstatin investigated by all-atom molecular dynamics simulations.

The development of compounds to regulate the activation of the complement system in non-primate species is of profound interest because it can provide models for human diseases. The peptide compstatin inhibits protein C3 in primate mammals and is a potential therapeutic agent against unregulated activation of complement in humans but is inactive against nonprimate species. Here, we elucidate this species specificity of compstatin by molecular dynamics simulations of complexes between the most potent natural compstatin analog and human or rat C3. The results are compared against an experimental conformation of the human complex, determined recently by X-ray diffraction at 2.4-A resolution. The human complex simulations provide information on the relative contributions to stability of specific C3 and compstatin residues. In the rat simulations, the protein undergoes reproducible conformational changes, which eliminate or weaken specific interactions and reduce the complex stability. The simulation insights can be used to design improved compstatin-based inhibitors for human C3 and active inhibitors against lower mammals.

Self-assembly of phenylalanine oligopeptides: insights from experiments and simulations.

Studies of peptide-based nanostructures provide general insights into biomolecular self-assembly and can lead material engineering toward technological applications. The diphenylalanine peptide (FF) self-assembles into discrete, hollow, well ordered nanotubes, and its derivatives form nanoassemblies of various morphologies. Here we demonstrate for the first time, to our knowledge, the formation of planar nanostructures with beta-sheet content by the triphenylalanine peptide (FFF). We characterize these structures using various microscopy and spectroscopy techniques. We also obtain insights into the interactions and structural properties of the FF and FFF nanostructures by 0.4-micros, implicit-solvent, replica-exchange, molecular-dynamics simulations of aqueous FF and FFF solutions. In the simulations the peptides form aggregates, which often contain open or ring-like peptide networks, as well as elementary and network-containing structures with beta-sheet characteristics. The networks are stabilized by polar and nonpolar interactions, and by the surrounding aggregate. In particular, the charged termini of neighbor peptides are involved in hydrogen-bonding interactions and their aromatic side chains form "T-shaped" contacts, as in three-dimensional FF crystals. These interactions may assist the FF and FFF self-assembly at the early stage, and may also stabilize the mature nanostructures. The FFF peptides have higher network propensities and increased aggregate stabilities with respect to FF, which can be interpreted energetically.

Glucose-based spiro-isoxazolines: a new family of potent glycogen phosphorylase inhibitors.

A series of glucopyranosylidene-spiro-isoxazolines was prepared through regio- and stereoselective [3+2]-cycloaddition between the methylene acetylated exo-glucal and aromatic nitrile oxides. The deprotected cycloadducts were evaluated as inhibitors of muscle glycogen phosphorylase b. The carbohydrate-based family of five inhibitors displays K(i) values ranging from 0.63 to 92.5 microM. The X-ray structures of the enzyme-ligand complexes show that the inhibitors bind preferentially at the catalytic site of the enzyme retaining the less active T-state conformation. Docking calculations with GLIDE in extra-precision (XP) mode yielded excellent agreement with experiment, as judged by comparison of the predicted binding modes of the five ligands with the crystallographic conformations and the good correlation between the docking scores and the experimental free binding energies. Use of docking constraints on the well-defined positions of the glucopyranose moiety in the catalytic site and redocking of GLIDE-XP poses using electrostatic potential fit-determined ligand partial charges in quantum polarized ligand docking (QPLD) produced the best results in this regard.

Toward Rational Design of Selective Molecularly Imprinted Polymers (MIPs) for Proteins: Computational and Experimental Studies of Acrylamide Based Polymers for Myoglobin.

Molecularly imprinted polymers (MIPs) have potential as alternatives to antibodies in the diagnosis and treatment of disease. However, atomistic level knowledge of the prepolymerization process is limited that would facilitate rational design of more efficient MIPs. Accordingly, we have investigated using computation and experiment the protein-monomer binding interactions that may influence the desired specificity. Myoglobin was used as the target protein and five different acrylamide-based monomers were considered. Protein binding sites were predicted using SiteMap and binding free energies of monomers at each site were calculated using MM-GBSA. Statistical thermodynamic analysis and study of atomistic interactions facilitated rationalization of monomer performance in MIP rebinding studies (% rebind; imprinting factors). CD spectroscopy was used to determine monomer effects on myoglobin secondary structure, with all monomers except the smallest monomer (acrylamide) causing significant changes. A complex interplay between different protein-monomer binding effects and MIP efficacy was observed. Validation of hypotheses for key binding features was achieved by rational selection of two different comonomer MIP combinations that produced experimental results in agreement with predictions. The comonomer studies revealed that uniform, noncompetitive binding of monomers around a target protein is favorable. This study represents a step toward future rational in silico design of MIPs for proteins.

Crystallographic and computational studies on 4-phenyl-N-(beta-D-glucopyranosyl)-1H-1,2,3-triazole-1-acetamide, an inhibitor of glycogen phosphorylase: comparison with alpha-D-glucose, N-acetyl-beta-D-glucopyranosylamine and N-benzoyl-N'-beta-D-glucopyranosyl urea binding.

4-Phenyl-N-(beta-D-glucopyranosyl)-1H-1,2,3-triazole-1-acetamide (glucosyltriazolylacetamide) has been studied in kinetic and crystallographic experiments with glycogen phosphorylase b (GPb), in an effort to utilize its potential as a lead for the design of potent antihyperglycaemic agents. Docking and molecular dynamics (MD) calculations have been used to monitor more closely the binding modes in operation and compare the results with experiment. Kinetic experiments in the direction of glycogen synthesis showed that glucosyltriazolylacetamide is a better inhibitor (K(i) = 0.18 mM) than the parent compound alpha-D-glucose (K(i) = 1.7 mM) or beta-D-glucose (K(i) = 7.4 mM) but less potent inhibitor than the lead compound N-acetyl-beta-D-glucopyranosylamine (K(i) = 32 microM). To elucidate the molecular basis underlying the inhibition of the newly identified compound, we determined the structure of GPb in complex with glucosyltriazolylacetamide at 100 K to 1.88 A resolution, and the structure of the compound in the free form. Glucosyltriazolylacetamide is accommodated in the catalytic site of the enzyme and the glucopyranose interacts in a manner similar to that observed in the GPb-alpha-D-glucose complex, while the substituent group in the beta-position of the C1 atom makes additional hydrogen bonding and van der Waals interactions to the protein. A bifurcated donor type hydrogen bonding involving O3H, N3, and N4 is seen as an important structural motif strengthening the binding of glucosyltriazolylacetamide with GP which necessitated change in the torsion about C8-N2 bond by about 62 degrees going from its free to the complex form with GPb. On binding to GP, glucosyltriazolylacetamide induces significant conformational changes in the vicinity of this site. Specifically, the 280s loop (residues 282-288) shifts 0.7 to 3.1 A (CA atoms) to accommodate glucosyltriazolylacetamide. These conformational changes do not lead to increased contacts between the inhibitor and the protein that would improve ligand binding compared with the lead compound. In the molecular modeling calculations, the GOLD docking runs with and without the crystallographic ordered cavity waters using the GoldScore scoring function, and without cavity waters using the ChemScore scoring function successfully reproduced the crystallographic binding conformation. However, the GLIDE docking calculations both with (GLIDE XP) and without (GLIDE SP and XP) the cavity water molecules were, impressively, further able to accurately reproduce the finer details of the GPb-glucosyltriazolylacetamide complex structure. The importance of cavity waters in flexible receptor MD calculations compared to "rigid" (docking) is analyzed and highlighted, while in the MD itself very little conformational flexibility of the glucosyltriazolylacetamide ligand was observed over the time scale of the simulations.

Recognition of ribonuclease A by 3'-5'-pyrophosphate-linked dinucleotide inhibitors: a molecular dynamics/continuum electrostatics analysis.

The proteins of the pancreatic ribonuclease A (RNase A) family catalyze the cleavage of the RNA polymer chain. The development of RNase inhibitors is of significant interest, as some of these compounds may have a therapeutic effect in pathological conditions associated with these proteins. The most potent low molecular weight inhibitor of RNase reported to date is the compound 5'-phospho-2'-deoxyuridine-3-pyrophosphate (P-->5)-adenosine-3-phosphate (pdUppA-3'-p). The 3',5'-pyrophosphate group of this compound increases its affinity and introduces structural features which seem to be unique in pyrophosphate-containing ligands bound to RNase A, such as the adoption of a syn conformation by the adenosine base at RNase subsite B(2) and the placement of the 5'-beta-phosphate of the adenylate (instead of the alpha-phosphate) at subsite P(1) where the phosphodiester bond cleavage occurs. In this work, we study by multi-ns molecular dynamics simulations the structural properties of RNase A complexes with the ligand pdUppA-3'-p and the related weaker inhibitor dUppA, which lacks the 3' and 5' terminal phosphate groups of pdUppA-3'-p. The simulations show that the adenylate 5'-beta-phosphate binding position and the adenosine syn orientation constitute robust structural features in both complexes, stabilized by persistent interactions with specific active-site residues of subsites P(1) and B(2). The simulation structures are used in conjunction with a continuum-electrostatics (Poisson-Boltzmann) model, to evaluate the relative binding affinity of the two complexes. The computed relative affinity of pdUppA-3'-p varies between -7.9 kcal/mol and -2.8 kcal/mol for a range of protein/ligand dielectric constants (epsilon(p)) 2-20, in good agreement with the experimental value (-3.6 kcal/mol); the agreement becomes exact with epsilon(p) = 8. The success of the continuum-electrostatics model suggests that the differences in affinity of the two ligands originate mainly from electrostatic interactions. A residue decomposition of the electrostatic free energies shows that the terminal phosphate groups of pdUppA-3'-p make increased interactions with residues Lys(7) and Lys(66) of the more remote sites P(2) and P(0), and His(119) of site P(1).

Computational sidechain placement and protein mutagenesis with implicit solvent models.

Structure prediction and computational protein design should benefit from accurate solvent models. We have applied implicit solvent models to two problems that are central to this area. First, we performed sidechain placement for 29 proteins, using a solvent model that combines a screened Coulomb term with an Accessible Surface Area term (CASA model). With optimized parameters, the prediction quality is comparable with earlier work that omitted electrostatics and solvation altogether. Second, we computed the stability changes associated with point mutations involving ionized sidechains. For over 1000 mutations, including many fully or partly buried positions, we compared CASA and two generalized Born models (GB) with a more accurate model, which solves the Poisson equation of continuum electrostatics numerically. CASA predicts the correct sign and order of magnitude of the stability change for 81% of the mutations, compared to 97% with the best GB. We also considered 140 mutations for which experimental data are available. Comparing to experiment requires additional assumptions about the unfolded protein structure, protein relaxation in response to the mutations, and contributions from the hydrophobic effect. With a simple, commonly-used unfolded state model, the mean unsigned error is 2.1 kcal/mol with both CASA and the best GB. Overall, the electrostatic model is not important for sidechain placement; CASA and GB are equivalent for surface mutations, while GB is far superior for fully or partly buried positions. Thus, for problems like protein design that involve all these aspects, the most recent GB models represent an important step forward. Along with the recent discovery of efficient, pairwise implementations of GB, this will open new possibilities for the computational engineering of proteins.