Pairwise decomposition of an MMGBSA energy function for computational protein design

Experimental and computational study on anti-gastric cancer activity and mechanism of evodiamine derivatives

Introduction: Human topoisomerase 1 (TOP1) is an important target of various anticancer compounds. The design and discovery of inhibitors targeting TOP1 are of great significance for the development of anticancer drugs. Evodiamine and thieno [2,3-d] pyridine hybrids show potential antitumor activity. Herein, the anti-gastric cancer activities of these hybrids were investigated. Methods: The inhibitory effects of different concentrations of ten evodiamine derivatives on the gastric cancer cell line SGC-7901 were assessed using a methyl thiazolyl tetrazolium assay. Compounds EVO-1 and EVO-6 strongly inhibited gastric cancer cell proliferation, with inhibition rates of 81.17% ± 5.08% and 80.92% ± 2.75%, respectively. To discover the relationship between the structure and activity of these two derivatives, density functional theory was used to investigate their optimized geometries, natural population charges, frontier molecular orbitals, and molecular electrostatic potentials. To clarify their anti-gastric cancer mechanisms, molecular docking, molecular dynamics simulations, and binding free energy calculations were performed against TOP1. Results: The results demonstrated that these compounds could intercalate into the cleaved DNA-binding site to form a TOP1-DNA-ligand ternary complex, and the ligand remained secure at the cleaved DNA-binding site to form a stable ternary complex. As the binding free energy of compound EVO-1 with TOP1 (-38.33 kcal·mol-1) was lower than that of compound EVO-6 (-33.25 kcal·mol-1), compound EVO-1 could be a more potent anti-gastric cancer agent than compound EVO-6. Discussion: Thus, compound EVO-1 could be a promising anti-gastric cancer drug candidate. This study may facilitate the design and development of novel TOP1 inhibitors.

Selectivity Mechanism of Pyrrolopyridone Analogues Targeting Bromodomain 2 of Bromodomain-Containing Protein 4 from Molecular Dynamics Simulations

Bromodomain and extra-terminal domain (BET) proteins play an important role in epigenetic regulation and are linked to several diseases; therefore, they are interesting targets. BET has two bromodomains: bromodomain 1 (BD1) and BD2. Selective targeting of BD1 or BD2 may produce different activities and greater effects than pan-BD inhibitors. However, the selective mechanism of the specific core must be studied at the atomic level. This study determined the effectiveness of pyrrolopyridone analogues to selectively inhibit BD2 using a pan-BD inhibitor (ABBV-075) and a selective-BD2 inhibitor (ABBV-744). Molecular dynamics simulations and calculations of binding free energies were used to systematically study the selectivity of BD2 inhibition by the pyrrolopyridone analogues. Overall, the pyrrolopyridone analogue inhibitors targeting BD2 interacted mainly with the following amino acid pairs between bromodomain-containing protein 4 (BRD4)-BD1 and BRD4-BD2 complexes: I146/V439, N140/N433, D144/H437, P82/P375, V87/V380, D88/D381, and Y139/Y432. The pyrrolopyridone analogues targeting BRD4-BD2 were divided into five regions based on selectivity mechanism. These results suggest that the R3 and R5 regions of pyrrolopyridone analogues can be modified to improve the selectivity between BRD4-BD1 and BRD4-BD2. The selectivity of BD2 inhibition by pyrrolopyridone analogues can be used to design novel BD2 inhibitors based on a pyrrolopyridone core.

Redesigning methionyl-tRNA synthetase for β-methionine activity with adaptive landscape flattening and experiments

Amino acids (AAs) with a noncanonical backbone would be a valuable tool for protein engineering, enabling new structural motifs and building blocks. To incorporate them into an expanded genetic code, the first, key step is to obtain an appropriate aminoacyl-tRNA synthetase. Currently, directed evolution is not available to optimize AAs with noncanonical backbones, since an appropriate selective pressure has not been discovered. Computational protein design (CPD) is an alternative. We used a new CPD method to redesign MetRS and increase its activity towards β-Met, which has an extra backbone methylene. The new method considered a few active site positions for design and used a Monte Carlo exploration of the corresponding sequence space. During the exploration, a bias energy was adaptively learned, such that the free energy landscape of the apo enzyme was flattened. Enzyme variants could then be sampled, in the presence of the ligand and the bias energy, according to their β-Met binding affinities. Eighteen predicted variants were chosen for experimental testing; 10 exhibited detectable activity for β-Met adenylation. Top predicted hits were characterized experimentally in detail. Dissociation constants, catalytic rates, and Michaelis constants for both α-Met and β-Met were measured. The best mutant retained a preference for α-Met over β-Met; however, the preference was reduced, compared to the wildtype, by a factor of 29. For this mutant, high resolution crystal structures were obtained in complex with both α-Met and β-Met, indicating that the predicted, active conformation of β-Met in the active site was retained.

Identification of abemaciclib derivatives targeting cyclin-dependent kinase 4 and 6 using molecular dynamics, binding free energy calculation, synthesis, and pharmacological evaluation

CDK4/6 plays a crucial role in various cancers and is an effective anticancer drug target. However, the gap between clinical requirements and approved CDK4/6 drugs is unresolved. Thus, there is an urgent need to develop selective and oral CDK4/6 inhibitors, particularly for monotherapy. Here, we studied the interaction between abemaciclib and human CDK6 using molecular dynamics simulations, binding free energy calculations, and energy decomposition. V101 and H100 formed stable hydrogen bonds with the amine-pyrimidine group, and K43 interacted with the imidazole ring via an unstable hydrogen bond. Meanwhile, I19, V27, A41, and L152 interacted with abemaciclib through π-alkyl interactions. Based on the binding model, abemaciclib was divided into four regions. With one region modification, 43 compounds were designed and evaluated using molecular docking. From each region, three favorable groups were selected and combined with each other to obtain 81 compounds. Among them, C2231-A, which was obtained by removing the methylene group from C2231, showed better inhibition than C2231. Kinase profiling revealed that C2231-A showed inhibitory activity similar to that of abemaciclib; additionally, C2231-A inhibited the growth of MDA-MB-231 cells to a greater extent than did abemaciclib. Based on molecular dynamics simulation, C2231-A was identified as a promising candidate compound with considerable inhibitory effects on human breast cancer cell lines.

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.

Computational protein design repurposed to explore enzyme vitality and help predict antibiotic resistance

In response to antibiotics that inhibit a bacterial enzyme, resistance mutations inevitably arise. Predicting them ahead of time would aid target selection and drug design. The simplest resistance mechanism would be to reduce antibiotic binding without sacrificing too much substrate binding. The property that reflects this is the enzyme "vitality", defined here as the difference between the inhibitor and substrate binding free energies. To predict such mutations, we borrow methodology from computational protein design. We use a Monte Carlo exploration of mutation space and vitality changes, allowing us to rank thousands of mutations and identify ones that might provide resistance through the simple mechanism considered. As an illustration, we chose dihydrofolate reductase, an essential enzyme targeted by several antibiotics. We simulated its complexes with the inhibitor trimethoprim and the substrate dihydrofolate. 20 active site positions were mutated, or "redesigned" individually, then in pairs or quartets. We computed the resulting binding free energy and vitality changes. Out of seven known resistance mutations involving active site positions, five were correctly recovered. Ten positions exhibited mutations with significant predicted vitality gains. Direct couplings between designed positions were predicted to be small, which reduces the combinatorial complexity of the mutation space to be explored. It also suggests that over the course of evolution, resistance mutations involving several positions do not need the underlying point mutations to arise all at once: they can appear and become fixed one after the other.

Interactions between curcumin and human salt-induced kinase 3 elucidated from computational tools and experimental methods

Natural products are widely used for treating mitochondrial dysfunction-related diseases and cancers. Curcumin, a well-known natural product, can be potentially used to treat cancer. Human salt-induced kinase 3 (SIK3) is one of the target proteins for curcumin. However, the interactions between curcumin and human SIK3 have not yet been investigated in detail. In this study, we studied the binding models for the interactions between curcumin and human SIK3 using computational tools such as homology modeling, molecular docking, molecular dynamics simulations, and binding free energy calculations. The open activity loop conformation of SIK3 with the ketoenol form of curcumin was the optimal binding model. The I72, V80, A93, Y144, A145, and L195 residues played a key role for curcumin binding with human SIK3. The interactions between curcumin and human SIK3 were also investigated using the kinase assay. Moreover, curcumin exhibited an IC₅₀ (half-maximal inhibitory concentration) value of 131 nM, and it showed significant antiproliferative activities of 9.62 ± 0.33 µM and 72.37 ± 0.37 µM against the MCF-7 and MDA-MB-23 cell lines, respectively. This study provides detailed information on the binding of curcumin with human SIK3 and may facilitate the design of novel salt-inducible kinases inhibitors.

Evodiamine and Rutaecarpine as Potential Anticancer Compounds: A Combined Computational Study

Evodiamine (EVO) and rutaecarpine (RUT) are the main active compounds of the traditional Chinese medicinal herb Evodia rutaecarpa. Here, we fully optimized the molecular geometries of EVO and RUT at the B3LYP/6-311++G (d, p) level of density functional theory. The natural population analysis (NPA) charges, frontier molecular orbitals, molecular electrostatic potentials, and the chemical reactivity descriptors for EVO and RUT were also investigated. Furthermore, molecular docking, molecular dynamics simulations, and the analysis of the binding free energies of EVO and RUT were carried out against the anticancer target topoisomerase 1 (TOP1) to clarify their anticancer mechanisms. The docking results indicated that they could inhibit TOP1 by intercalating into the cleaved DNA-binding site to form a TOP1−DNA−ligand ternary complex, suggesting that they may be potential TOP1 inhibitors. Molecular dynamics (MD) simulations evaluated the binding stability of the TOP1−DNA−ligand ternary complex. The calculation of binding free energy showed that the binding ability of EVO with TOP1 was stronger than that of RUT. These results elucidated the structure−activity relationship and the antitumor mechanism of EVO and RUT at the molecular level. It is suggested that EVO and RUT may be potential compounds for the development of new anticancer drugs.

Molecular Dynamics Simulations of the Conformational Plasticity in the Active Pocket of Salt-Inducible Kinase 2 (SIK2) Multi-State Binding with Bosutinib

The kinase domain is highly conserved among protein kinases 'in terms of both sequence and structure. Conformational rearrangements of the kinase domain are affected by the phosphorylation of residues and the binding of kinase inhibitors. Interestingly, the conformational rearrangement of the active pocket plays an important role in kinase activity and can be used to design novel kinase inhibitors. We characterized the conformational plasticity of the active pocket when bosutinib was bound to salt-inducible kinase 2 (SIK2) using homology modeling and molecular dynamics simulations. Ten different initial complex models were constructed using the Morph server, ranging from open to closed conformations of SIK2 binding with bosutinib. Our simulation showed that bosutinib binds SIK2 with up or down conformations of the P-loop and with all the conformations of the activation loop. In addition, the αC-helix conformation was induced by the conformation of the activation loop, and the salt bridge formed only with its open conformation. The binding affinity of the models was also determined using the molecular mechanics generalized Born surface area method. Bosutinib was found to form a strong binding model with SIK2 and hydrophobic interactions were the dominant factor. This discovery may help guide the design of novel SIK2 inhibitors.

The binding mechanism of NHWD-870 to bromodomain-containing protein 4 based on molecular dynamics simulations and free energy calculation

Bromodomain and extra-terminal (BET) proteins (BRD2, BRD3, BRD4, and BRDT) are epigenetic readers with tandem bromodomains.

Perspective of structural flexibility on selective inhibition towards CYP1B1 over CYP1A1 by α-naphthoflavone analogs

Research on action selectivity between CYP1A1 and CYP1B1 is particularly valuable for cancer chemoprevention and chemotherapy. However, they share a very close similarity in their ligand-binding pockets that α-naphthoflavone (ANF) is the co-crystal ligand for both isoforms, which poses a major challenge in revealing their selectivity mechanism. Therefore, three selective CYP1B1 inhibitors derived from ANF were selected to illustrate the structural basis for the selectivity between the two isoforms via a comprehensive computational strategy. It was found that the sustainability of the π-π stacking interactions with the phenylalanine residues of the two isoforms, namely, Phe123, Phe224, and Phe258 for CYP1A1, and Phe134, Phe231, and Phe268 for CYP1B1, played a crucial role in determining the selectivity of ligands with a classic aromatic conjugation system like ANF and its derivatives for CYP1B1 versus CYP1A1. Of note, the structural flexibility of the corresponding protein domains mainly orchestrated the sustainability of the corresponding π-π stacking interactions, thereby determining the binding selectivity. Therefore, the structure modification of naphthoflavone lead compounds into preferable binding configurations to satisfy the π-π stacking interactions of the key phenylalanine residues within CYP1B1 would be an inspiring strategy devised to improve the inhibitory selectivity towards CYP1B1. Collectively, this study revealed valuable insight into understanding the selective mechanism between CYP1A1 and CYP1B1 from the perspective of structural flexibility, which sheds light on the future rational design of CYP1B1 selective inhibitors.

Exploring the stability of inhibitor binding to SIK2 using molecular dynamics simulation and binding free energy calculation

Salt inducible kinase 2 (SIK2) is a calcium/calmodulin-dependent protein kinase-like kinase that is implicated in a variety of biological phenomena, including cellular metabolism, growth, and apoptosis. SIK2 is the key target for various cancers, including ovarian, breast, prostate, and lung cancers. Although potent inhibitors of SIK2 are being developed, their binding stability and functional role are not presently known. In this work, we studied the detailed interactions between SIK2 and four of its inhibitors, HG-9-91-01, KIN112, MRT67307, and MRT199665, using molecular docking, molecular dynamics simulation, binding free energy calculation, and interaction fingerprint analysis. Intermolecular interactions revealed that HG-9-91-01 and KIN112 have stronger interactions with SIK2 than those of MRT199665 and MRT67307. The key residues involved in binding with SIK2 are conserved among all four inhibitors. Our results explain the detailed interaction of SIK2 with its inhibitors at the molecular level, thus paving the way for the development of targeted efficient anti-cancer drugs.

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.

Dasatinib-SIK2 Binding Elucidated by Homology Modeling, Molecular Docking, and Dynamics Simulations

Salt-inducible kinases (SIKs) are calcium/calmodulin-dependent protein kinase (CAMK)-like (CAMKL) family members implicated in insulin signal transduction, metabolic regulation, inflammatory response, and other processes. Here, we focused on SIK2, which is a target of the Food and Drug Administration (FDA)-approved pan inhibitor N-(2-chloro-6-methylphenyl)-2-(6-(4-(2-hydroxyethyl)piperazin-1-yl)-2-methylpyrimidin-4-ylamino)thiazole-5-carboxamide (dasatinib), and constructed four representative SIK2 structures by homology modeling. We investigated the interactions between dasatinib and SIK2 via molecular docking, molecular dynamics simulation, and binding free energy calculation and found that dasatinib showed strong binding affinity for SIK2. Binding free energy calculations suggested that the modification of various dasatinib regions may provide useful information for drug design and to guide the discovery of novel dasatinib-based SIK2 inhibitors.

"Big Data" Fast Chemoinformatics Model to Predict Generalized Born Radius and Solvent Accessibility as a Function of Geometry

The Generalized Born (GB) solvent model is offering the best accuracy/computing effort ratio yet requires drastic simplifications to estimate of the Effective Born Radii (EBR) in bypassing a too expensive volume integration step. EBRs are a measure of the degree of burial of an atom and not very sensitive to small changes of geometry: in molecular dynamics, the costly EBR update procedure is not mandatory at every step. This work however aims at implementing a GB model into the Sampler for Multiple Protein-Ligand Entities (S4MPLE) evolutionary algorithm with mandatory EBR updates at each step triggering arbitrarily large geometric changes. Therefore, a quantitative structure-property relationship has been developed in order to express the EBRs as a linear function of both the topological neighborhood and geometric occupancy of the space around atoms. A training set of 810 molecular systems, starting from fragment-like to drug-like compounds, proteins, host-guest systems, and ligand-protein complexes, has been compiled. For each species, S4MPLE generated several hundreds of random conformers. For each atom in each geometry of each species, its "standard" EBR was calculated by numeric integration and associated to topological and geometric descriptors of the atom neighborhood. This training set (EBR, atom descriptors) involving >5 M entries was subjected to a boot-strapping multilinear regression process with descriptor selection. In parallel, the strategy was repurposed to also learn atomic solvent-accessible areas (SA) based on the same descriptors. Resulting linear equations were challenged to predict EBR and SA values for a similarly compiled external set of >2000 new molecular systems. Solvation energies calculated with estimated EBR and SA match "standard" energies within the typical error of a force-field-based approach (a few kilocalories per mole). Given the extreme diversity of molecular systems covered by the model, this simple EBR/SA estimator covers a vast applicability domain.

Revealing the interaction modes of 5-HT2A receptor antagonists and the Structure-Based virtual screening from FDA and TCMNP database

5-hydroxytryptamine 2A (5-HT2A) receptor is emerging as an important target for numerous psychoactive drugs due to its imperative roles in psychological diseases. In fact, multiple 5-HT2A receptor antagonists were developed to treat numerous psychiatric disorders, however, their clinical outcome was far from ideal probably due to a blurry information of the exact interaction modes between the receptor and its antagonists. Impressively, with a recent release of its crystal structure, we carefully analyzed the receptor-ligand interactions with Protein Contacts Atlas, structure-based pharmacophore models, and molecular dynamics (MD) simulations to sum up the chemical features for antagonists interacting with 5-HT2A receptor. Moreover, the molecular docking-based virtual screening was applied to discover potential 5-HT2A receptor antagonists from FDA and TCMNP databases. Intriguingly, after a systematic assessment of the docking scores, binding modes and free energies, as well as their MD simulations performances, three compounds in TCMNP database were highlighted to be potential 5-HT2A receptor antagonists. Fascinatedly, these three hits also exhibited highly binding affinities with dopamine D2 receptor (D2R) due to the similarity of the ligand binding pockets of the receptors, indicating them to be promising dual target molecules that are of great benefit for anti-psychotic-drug research and development. In addition, ADME/Tox predictions were conducted for a primary evaluation of their developing potential. Together, this study not only revealed the exact interaction modes between 5-HT2A receptor and its antagonists, which shed a light on a better access for developing its novel antagonists, but also provided promising dual D2 and 5-HT2A receptor antagonists.Communicated by Ramaswamy H. Sarma.

Adaptive landscape flattening allows the design of both enzyme: Substrate binding and catalytic power

Designed enzymes are of fundamental and technological interest. Experimental directed evolution still has significant limitations, and computational approaches are a complementary route. A designed enzyme should satisfy multiple criteria: stability, substrate binding, transition state binding. Such multi-objective design is computationally challenging. Two recent studies used adaptive importance sampling Monte Carlo to redesign proteins for ligand binding. By first flattening the energy landscape of the apo protein, they obtained positive design for the bound state and negative design for the unbound. We have now extended the method to design an enzyme for specific transition state binding, i.e., for its catalytic power. We considered methionyl-tRNA synthetase (MetRS), which attaches methionine (Met) to its cognate tRNA, establishing codon identity. Previously, MetRS and other synthetases have been redesigned by experimental directed evolution to accept noncanonical amino acids as substrates, leading to genetic code expansion. Here, we have redesigned MetRS computationally to bind several ligands: the Met analog azidonorleucine, methionyl-adenylate (MetAMP), and the activated ligands that form the transition state for MetAMP production. Enzyme mutants known to have azidonorleucine activity were recovered by the design calculations, and 17 mutants predicted to bind MetAMP were characterized experimentally and all found to be active. Mutants predicted to have low activation free energies for MetAMP production were found to be active and the predicted reaction rates agreed well with the experimental values. We suggest the present method should become the paradigm for computational enzyme design.

Designed enzymes are of major interest. Experimental directed evolution still has significant limitations, and computational approaches are another route. Enzymes must be stable, bind substrates, and be powerful catalysts. It is challenging to design for all these properties. A method to design substrate binding was proposed recently. It used an adaptive Monte Carlo method to explore mutations of a few amino acids near the substrate. A bias energy was gradually “learned” such that, in the absence of the ligand, the simulation visited most of the possible protein mutations with comparable probabilities. Remarkably, a simulation of the protein:ligand complex, including the bias, will then preferentially sample tight-binding sequences. We generalized the method to design binding specificity. We tested it for the methionyl-tRNA synthetase enzyme, which has been engineered in order to expand the genetic code. We redesigned the enzyme to obtain variants with low activation free energies for the catalytic step. The variants proposed by the simulations were shown experimentally to be active, and the predicted activation free energies were in reasonable agreement with the experimental values. We expect the new method will become the paradigm for computational enzyme design.

Adaptive landscape flattening in amino acid sequence space for the computational design of protein:peptide binding

For the high throughput design of protein:peptide binding, one must explore a vast space of amino acid sequences in search of low binding free energies. This complex problem is usually addressed with either simple heuristic scoring or expensive sequence enumeration schemes. Far more efficient than enumeration is a recent Monte Carlo approach that adaptively flattens the energy landscape in sequence space of the unbound peptide and provides formally exact binding free energy differences. The method allows the binding free energy to be used directly as the design criterion. We propose several improvements that allow still more efficient sampling and can address larger design problems. They include the use of Replica Exchange Monte Carlo and landscape flattening for both the unbound and bound peptides. We used the method to design peptides that bind to the PDZ domain of the Tiam1 signaling protein and could serve as inhibitors of its activity. Four peptide positions were allowed to mutate freely. Almost 75 000 peptide variants were processed in two simulations of 10⁹ steps each that used 1 CPU hour on a desktop machine. 96% of the theoretical sequence space was sampled. The relative binding free energies agreed qualitatively with values from experiment. The sampled sequences agreed qualitatively with an experimental library of Tiam1-binding peptides. The main assumption limiting accuracy is the fixed backbone approximation, which could be alleviated in future work by using increased computational resources and multi-backbone designs.

Recent advances in automated protein design and its future challenges

INTRODUCTION

Protein function is determined by protein structure which is in turn determined by the corresponding protein sequence. If the rules that cause a protein to adopt a particular structure are understood, it should be possible to refine or even redefine the function of a protein by working backwards from the desired structure to the sequence. Automated protein design attempts to calculate the effects of mutations computationally with the goal of more radical or complex transformations than are accessible by experimental techniques. Areas covered: The authors give a brief overview of the recent methodological advances in computer-aided protein design, showing how methodological choices affect final design and how automated protein design can be used to address problems considered beyond traditional protein engineering, including the creation of novel protein scaffolds for drug development. Also, the authors address specifically the future challenges in the development of automated protein design. Expert opinion: Automated protein design holds potential as a protein engineering technique, particularly in cases where screening by combinatorial mutagenesis is problematic. Considering solubility and immunogenicity issues, automated protein design is initially more likely to make an impact as a research tool for exploring basic biology in drug discovery than in the design of protein biologics.

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.

Comparing pairwise-additive and many-body generalized Born models for acid/base calculations and protein design

Generalized Born (GB) solvent models are common in acid/base calculations and protein design. With GB, the interaction between a pair of solute atoms depends on the shape of the protein/solvent boundary and, therefore, the positions of all solute atoms, so that GB is a many-body potential. For compute-intensive applications, the model is often simplified further, by introducing a mean, native-like protein/solvent boundary, which removes the many-body property. We investigate a method for both acid/base calculations and protein design that uses Monte Carlo simulations in which side chains can explore rotamers, bind/release protons, or mutate. The fluctuating protein/solvent dielectric boundary is treated in a way that is numerically exact (within the GB framework), in contrast to a mean boundary. Its originality is that it captures the many-body character while retaining the residue-pairwise complexity given by a fixed boundary. The method is implemented in the Proteus protein design software. It yields a slight but systematic improvement for acid/base constants in nine proteins and a significant improvement for the computational design of three PDZ domains. It eliminates a source of model uncertainty, which will facilitate the analysis of other model limitations. © 2017 Wiley Periodicals, Inc.

Computational Design of the Tiam1 PDZ Domain and Its Ligand Binding

PDZ domains direct protein-protein interactions and serve as models for protein design. Here, we optimized a protein design energy function for the Tiam1 and Cask PDZ domains that combines a molecular mechanics energy, Generalized Born solvent, and an empirical unfolded state model. Designed sequences were recognized as PDZ domains by the Superfamily fold recognition tool and had similarity scores comparable to natural PDZ sequences. The optimized model was used to redesign the two PDZ domains, by gradually varying the chemical potential of hydrophobic amino acids; the tendency of each position to lose or gain a hydrophobic character represents a novel hydrophobicity index. We also redesigned four positions in the Tiam1 PDZ domain involved in peptide binding specificity. The calculated affinity differences between designed variants reproduced experimental data and suggest substitutions with altered specificities.

Computer-Aided Antibody Design: An Overview

The use of monoclonal antibody as the next generation protein therapeutics with remarkable success has surged the development of antibody engineering to design molecules for optimizing affinity, better efficacy, greater safety and therapeutic function. Therefore, computational methods have become increasingly important to generate hypotheses, interpret and guide experimental works. In this chapter, we discussed the overall antibody design by computational approches.

Simultaneous Optimization of Biomolecular Energy Functions on Features from Small Molecules and Macromolecules

Most biomolecular modeling energy functions for structure prediction, sequence design, and molecular docking have been parametrized using existing macromolecular structural data; this contrasts molecular mechanics force fields which are largely optimized using small-molecule data. In this study, we describe an integrated method that enables optimization of a biomolecular modeling energy function simultaneously against small-molecule thermodynamic data and high-resolution macromolecular structural data. We use this approach to develop a next-generation Rosetta energy function that utilizes a new anisotropic implicit solvation model, and an improved electrostatics and Lennard-Jones model, illustrating how energy functions can be considerably improved in their ability to describe large-scale energy landscapes by incorporating both small-molecule and macromolecule data. The energy function improves performance in a wide range of protein structure prediction challenges, including monomeric structure prediction, protein-protein and protein-ligand docking, protein sequence design, and prediction of the free energy changes by mutation, while reasonably recapitulating small-molecule thermodynamic properties.

Comparing three stochastic search algorithms for computational protein design: Monte Carlo, replica exchange Monte Carlo, and a multistart, steepest-descent heuristic

Computational protein design depends on an energy function and an algorithm to search the sequence/conformation space. We compare three stochastic search algorithms: a heuristic, Monte Carlo (MC), and a Replica Exchange Monte Carlo method (REMC). The heuristic performs a steepest-descent minimization starting from thousands of random starting points. The methods are applied to nine test proteins from three structural families, with a fixed backbone structure, a molecular mechanics energy function, and with 1, 5, 10, 20, 30, or all amino acids allowed to mutate. Results are compared to an exact, "Cost Function Network" method that identifies the global minimum energy conformation (GMEC) in favorable cases. The designed sequences accurately reproduce experimental sequences in the hydrophobic core. The heuristic and REMC agree closely and reproduce the GMEC when it is known, with a few exceptions. Plain MC performs well for most cases, occasionally departing from the GMEC by 3-4 kcal/mol. With REMC, the diversity of the sequences sampled agrees with exact enumeration where the latter is possible: up to 2 kcal/mol above the GMEC. Beyond, room temperature replicas sample sequences up to 10 kcal/mol above the GMEC, providing thermal averages and a solution to the inverse protein folding problem. © 2016 Wiley Periodicals, Inc.

Protein side chain conformation predictions with an MMGBSA energy function

The prediction of protein side chain conformations from backbone coordinates is an important task in structural biology, with applications in structure prediction and protein design. It is a difficult problem due to its combinatorial nature. We study the performance of an "MMGBSA" energy function, implemented in our protein design program Proteus, which combines molecular mechanics terms, a Generalized Born and Surface Area (GBSA) solvent model, with approximations that make the model pairwise additive. Proteus is not a competitor to specialized side chain prediction programs due to its cost, but it allows protein design applications, where side chain prediction is an important step and MMGBSA an effective energy model. We predict the side chain conformations for 18 proteins. The side chains are first predicted individually, with the rest of the protein in its crystallographic conformation. Next, all side chains are predicted together. The contributions of individual energy terms are evaluated and various parameterizations are compared. We find that the GB and SA terms, with an appropriate choice of the dielectric constant and surface energy coefficients, are beneficial for single side chain predictions. For the prediction of all side chains, however, errors due to the pairwise additive approximation overcome the improvement brought by these terms. We also show the crucial contribution of side chain minimization to alleviate the rigid rotamer approximation. Even without GB and SA terms, we obtain accuracies comparable to SCWRL4, a specialized side chain prediction program. In particular, we obtain a better RMSD than SCWRL4 for core residues (at a higher cost), despite our simpler rotamer library. Proteins 2016; 84:803-819. © 2016 Wiley Periodicals, Inc.

Proteus and the Design of Ligand Binding Sites

This chapter describes the organization and use of Proteus, a multitool computational suite for the optimization of protein and ligand conformations and sequences, and the calculation of pK α shifts and relative binding affinities. The software offers the use of several molecular mechanics force fields and solvent models, including two generalized Born variants, and a large range of scoring functions, which can combine protein stability, ligand affinity, and ligand specificity terms, for positive and negative design. We present in detail the steps for structure preparation, system setup, construction of the interaction energy matrix, protein sequence and structure optimizations, pK α calculations, and ligand titration calculations. We discuss illustrative examples, including the chemical/structural optimization of a complex between the MHC class II protein HLA-DQ8 and the vinculin epitope, and the chemical optimization of the compstatin analog Ac-Val4Trp/His9Ala, which regulates the function of protein C3 of the complement system.

Influence of Chirality of Crizotinib on Its MTH1 Protein Inhibitory Activity: Insight from Molecular Dynamics Simulations and Binding Free Energy Calculations

As a promising target for the treatment of lung cancer, the MutT Homolog 1 (MTH1) protein can be inhibited by crizotinib. A recent work shows that the inhibitory potency of (S)-crizotinib against MTH1 is about 20 times over that of (R)-crizotinib. But the detailed molecular mechanism remains unclear. In this study, molecular dynamics (MD) simulations and free energy calculations were used to elucidate the mechanism about the effect of chirality of crizotinib on the inhibitory activity against MTH1. The binding free energy of (S)-crizotinib predicted by the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) and Adaptive biasing force (ABF) methodologies is much lower than that of (R)-crizotinib, which is consistent with the experimental data. The analysis of the individual energy terms suggests that the van der Waals interactions are important for distinguishing the binding of (S)-crizotinib and (R)-crizotinib. The binding free energy decomposition analysis illustrated that residues Tyr7, Phe27, Phe72 and Trp117 were important for the selective binding of (S)-crizotinib to MTH1. The adaptive biasing force (ABF) method was further employed to elucidate the unbinding process of (S)-crizotinib and (R)-crizotinib from the binding pocket of MTH1. ABF simulation results suggest that the reaction coordinates of the (S)-crizotinib from the binding pocket is different from (R)-crizotinib. The results from our study can reveal the details about the effect of chirality on the inhibition activity of crizotinib to MTH1 and provide valuable information for the design of more potent inhibitors.

Protein:Ligand binding free energies: A stringent test for computational protein design

A computational protein design method is extended to allow Monte Carlo simulations where two ligands are titrated into a protein binding pocket, yielding binding free energy differences. These provide a stringent test of the physical model, including the energy surface and sidechain rotamer definition. As a test, we consider tyrosyl-tRNA synthetase (TyrRS), which has been extensively redesigned experimentally. We consider its specificity for its substrate l-tyrosine (l-Tyr), compared to the analogs d-Tyr, p-acetyl-, and p-azido-phenylalanine (ac-Phe, az-Phe). We simulate l- and d-Tyr binding to TyrRS and six mutants, and compare the structures and binding free energies to a more rigorous "MD/GBSA" procedure: molecular dynamics with explicit solvent for structures and a Generalized Born + Surface Area model for binding free energies. Next, we consider l-Tyr, ac- and az-Phe binding to six other TyrRS variants. The titration results are sensitive to the precise rotamer definition, which involves a short energy minimization for each sidechain pair to help relax bad contacts induced by the discrete rotamer set. However, when designed mutant structures are rescored with a standard GBSA energy model, results agree well with the more rigorous MD/GBSA. As a third test, we redesign three amino acid positions in the substrate coordination sphere, with either l-Tyr or d-Tyr as the ligand. For two, we obtain good agreement with experiment, recovering the wildtype residue when l-Tyr is the ligand and a d-Tyr specific mutant when d-Tyr is the ligand. For the third, we recover His with either ligand, instead of wildtype Gln.

Dead-End Elimination with a Polarizable Force Field Repacks PCNA Structures

A balance of van der Waals, electrostatic, and hydrophobic forces drive the folding and packing of protein side chains. Although such interactions between residues are often approximated as being pairwise additive, in reality, higher-order many-body contributions that depend on environment drive hydrophobic collapse and cooperative electrostatics. Beginning from dead-end elimination, we derive the first algorithm, to our knowledge, capable of deterministic global repacking of side chains compatible with many-body energy functions. The approach is applied to seven PCNA x-ray crystallographic data sets with resolutions 2.5-3.8 Å (mean 3.0 Å) using an open-source software. While PDB_REDO models average an Rfree value of 29.5% and MOLPROBITY score of 2.71 Å (77th percentile), dead-end elimination with the polarizable AMOEBA force field lowered Rfree by 2.8-26.7% and improved mean MOLPROBITY score to atomic resolution at 1.25 Å (100th percentile). For structural biology applications that depend on side-chain repacking, including x-ray refinement, homology modeling, and protein design, the accuracy limitations of pairwise additivity can now be eliminated via polarizable or quantum mechanical potentials.

Tunable Coarse Graining for Monte Carlo Simulations of Proteins via Smoothed Energy Tables: Direct and Exchange Simulations

Many commonly used coarse-grained models for proteins are based on simplified interaction sites and consequently may suffer from significant limitations, such as the inability to properly model protein secondary structure without the addition of restraints. Recent work on a benzene fluid (Lettieri S.; Zuckerman D. M.J. Comput. Chem.2012, 33, 268-275) suggested an alternative strategy of tabulating and smoothing fully atomistic orientation-dependent interactions among rigid molecules or fragments. Here we report our initial efforts to apply this approach to the polar and covalent interactions intrinsic to polypeptides. We divide proteins into nearly rigid fragments, construct distance and orientation-dependent tables of the atomistic interaction energies between those fragments, and apply potential energy smoothing techniques to those tables. The amount of smoothing can be adjusted to give coarse-grained models that range from the underlying atomistic force field all the way to a bead-like coarse-grained model. For a moderate amount of smoothing, the method is able to preserve about 70-90% of the α-helical structure while providing a factor of 3-10 improvement in sampling per unit computation time (depending on how sampling is measured). For a greater amount of smoothing, multiple folding-unfolding transitions of the peptide were observed, along with a factor of 10-100 improvement in sampling per unit computation time, although the time spent in the unfolded state was increased compared with less smoothed simulations. For a β hairpin, secondary structure is also preserved, albeit for a narrower range of the smoothing parameter and, consequently, for a more modest improvement in sampling. We have also applied the new method in a "resolution exchange" setting, in which each replica runs a Monte Carlo simulation with a different degree of smoothing. We obtain exchange rates that compare favorably to our previous efforts at resolution exchange (Lyman E.; Zuckerman D. M.J. Chem. Theory Comput.2006, 2, 656-666).