Identifying promising druggable binding sites and their flexibility to target the receptor-binding domain of SARS-CoV-2 spike protein

The spike protein of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is crucial for viral infection. The interaction of its receptor-binding domain (RBD) with the human angiotensin-converting enzyme 2 (ACE2) protein is required for the virus to enter the host cell. We identified RBD binding sites to block its function with inhibitors by combining the protein structural flexibility with machine learning analysis. Molecular dynamics simulations were performed on unbound or ACE2-bound RBD conformations. Pockets estimation, tracking and druggability prediction were performed on a large sample of simulated RBD conformations. Recurrent druggable binding sites and their key residues were identified by clustering pockets based on their residue similarity. This protocol successfully identified three druggable sites and their key residues, aiming to target with inhibitors for preventing ACE2 interaction. One site features key residues for direct ACE2 interaction, highlighted using energetic computations, but can be affected by several mutations of the variants of concern. Two highly druggable sites, located between the spike protein monomers interface are promising. One weakly impacted by only one Omicron mutation, could contribute to stabilizing the spike protein in its closed state. The other, currently not affected by mutations, could avoid the activation of the spike protein trimer.

Molecular Modeling Study of Uncharged Oximes Compared to HI-6 and 2-PAM Inside Human AChE Sarin and VX Conjugates

The deleterious effects of nerve agents over the enzymes acetylcholinesterase (AChE) and butyrylcholinesterase (BChE) turned these compounds into the most dangerous chemical weapons known. Among the antidotes in use today against these agents, oximes in combination with other drugs are the only treatment with any action. HI-6 and 2-PAM are cationic oximes proved to be effective for the reactivation of AChE inhibited by the nerve agents VX and sarin (GB). However, when it comes to reactivation of AChE inside the central or peripheral nervous systems, charged molecules present low diffusion due to low penetration through the blood-brain barrier. Uncharged oximes appear as an interesting alternative to solve this problem, but the development and enhancement of more efficient uncharged oximes capable of reactivating human AChE is still necessary. Given the limitations for in vivo and in vitro experimental studies with nerve agents, modeling is an important tool that can contribute to a better understanding of factors that may affect the efficiency of uncharged oximes. In order to investigate the interaction and behavior of cationic and uncharged oximes, we performed here molecular docking, molecular dynamics simulations, and binding energies calculations of the known cationic oximes HI-6 and 2-PAM plus four uncharged oximes found in the literature, complexed with human AChE (HssACHE) conjugated with the nerve agents VX and GB. The uncharged oximes showed different behaviors, especially RS194B, which presented stability inside AChE-VX, but presented free binding energy lower than cationic oximes, suggesting that structural alterations could favor its interactions with these complexes. In contrast, HI-6 and 2-PAM showed higher affinities with more negative binding energy values and larger contribution of the amino acid Asp74, demonstrating the importance of the quaternary nitrogen to the affinity and interaction of oximes with AChE-GB and AChE-VX conjugates.

Investigating Phosphorylation-Induced Conformational Changes in WNK1 Kinase by Molecular Dynamics Simulations

The With-No-Lysine (WNK) kinase is considered to be a master regulator for various cation-chloride cotransporters involved in maintaining cell-volume and ion homeostasis. Here, we have investigated the phosphorylation-induced structural dynamics of the WNK1 kinase bound to an inhibitor via atomistic molecular dynamics simulations. Results from our simulations show that the phosphorylation at Ser³⁸² could stabilize the otherwise flexible activation loop (A-loop). The intrahelix salt-bridge formed between Arg²⁶⁴ and Glu²⁶⁸ in the unphosphorylated system is disengaged after the phosphorylation, and Glu²⁶⁸ reorients itself and forms a stable salt-bridge with Arg³⁴⁸. The dynamic cross-correlation analysis shows that phosphorylation diminishes anticorrelated motions and increases correlated motions between different domains. Structural network analysis reveals that the phosphorylation causes structural rearrangements and shortens the communication path between the αC-helix and catalytic loop, making the binding pocket more suitable for accommodating the ligand. Overall, we have characterized the structural changes in the WNK kinase because of phosphorylation in the A-loop, which might help in designing rational drugs.

Interactions of human butyrylcholinesterase with phenylvalerate and acetylthiocholine as substrates and inhibitors: kinetic and molecular modeling approaches

Phenyl valerate (PV) is a substrate for measuring the PVase activity of neuropathy target esterase (NTE), a key molecular event of organophosphorus-induced delayed neuropathy. A protein with PVase activity in chicken (model for delayed neurotoxicity) was identified as butyrylcholinesterase (BChE). Purified human butyrylcholinesterase (hBChE) showed PVase activity with a similar sensitivity to inhibitors as its cholinesterase (ChE) activity. Further kinetic and theoretical molecular simulation studies were performed. The kinetics did not fit classic competition models among substrates. Partially mixed inhibition was the best-fitting model to acetylthiocholine (AtCh) interacting with PVase activity. ChE activity showed substrate activation, and non-competitive inhibition was the best-fitting model to PV interacting with the non-activated enzyme and partial non-competitive inhibition was the best fitted model for PV interacting with the activated enzyme by excess of AtCh. The kinetic results suggest that other sites could be involved in those activities. From the theoretical docking analysis, we deduced other more favorable sites for binding PV related with Asn289 residue, situated far from the catalytic site ("PV-site"). Both substrates acethylcholine (ACh) and PV presented similar docking values in both the PV-site and catalytic site pockets, which explained some of the observed substrate interactions. Molecular dynamic simulations based on the theoretical structure of crystallized hBChE were performed. Molecular modeling studies suggested that PV has a higher potential for non-competitive inhibition, being also able to inhibit the hydrolysis of ACh through interactions with the PV-site. Further theoretical studies also suggested that PV could yet be able to promote competitive inhibition. We concluded that the kinetic and theoretical studies did not fit the simple classic competition among substrates, but were compatible with the interaction with two different binding sites.

Design of inhibitors of thymidylate kinase from Variola virus as new selective drugs against smallpox: part II

Acknowledging the importance of studies toward the development of measures against terrorism and bioterrorism, this study aims to contribute to the design of new prototypes of potential drugs against smallpox. Based on a former study, nine synthetic feasible prototypes of selective inhibitors for thymidylate kinase from Variola virus (VarTMPK) were designed and submitted to molecular docking, molecular dynamics simulations and binding energy calculations. The compounds are simplifications of two more complex scaffolds, with a guanine connected to an amide or alcohol through a spacer containing ether and/or amide groups, formerly suggested as promising for the design of selective inhibitors of VarTMPK. Our study showed that, despite the structural simplifications, the compounds presented effective energy values in interactions with VarTMPK and HssTMPK and that the guanine could be replaced by a simpler imidazole ring linked to a –NH₂ group, without compromising the affinity for VarTMPK. It was also observed that a positive charge in the imidazole ring is important for the selectivity toward VarTMPK and that an amide group in the spacer does not contribute to selectivity. Finally, prototype 3 was pointed as the most promising to be synthesized and experimentally evaluated.

Communicated by Ramaswamy H. Sarma

RNA Structural Dynamics As Captured by Molecular Simulations: A Comprehensive Overview

With both catalytic and genetic functions, ribonucleic acid (RNA) is perhaps the most pluripotent chemical species in molecular biology, and its functions are intimately linked to its structure and dynamics. Computer simulations, and in particular atomistic molecular dynamics (MD), allow structural dynamics of biomolecular systems to be investigated with unprecedented temporal and spatial resolution. We here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA-ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field.

Coupled molecular dynamics and continuum electrostatic method to compute the ionization pKa's of proteins as a function of pH. Test on a large set of proteins

A computational method, to predict the pKa values of the ionizable residues Asp, Glu, His, Tyr, and Lys of proteins, is presented here. Calculation of the electrostatic free-energy of the proteins is based on an efficient version of a continuum dielectric electrostatic model. The conformational flexibility of the protein is taken into account by carrying out molecular dynamics simulations of 10 ns in implicit water. The accuracy of the proposed method of calculation of pKa values is estimated from a test set of experimental pKa data for 297 ionizable residues from 34 proteins. The pKa-prediction test shows that, on average, 57, 86, and 95% of all predictions have an error lower than 0.5, 1.0, and 1.5 pKa units, respectively. This work contributes to our general understanding of the importance of protein flexibility for an accurate computation of pKa, providing critical insight about the significance of the multiple neutral states of acid and histidine residues for pKa-prediction, and may spur significant progress in our effort to develop a fast and accurate electrostatic-based method for pKa-predictions of proteins as a function of pH.

Behavior of uncharged oximes compared to HI6 and 2-PAM in the human AChE-tabun conjugate: a molecular modeling approach

Tabun is one of the most dangerous nerve agents because it has deleterious effects like inhibition of the essential enzymes acetylcholinesterase (AChE) and butyrylcholinesterase. Some oximes such HI6 as 2-PAM are nucleophiles that are capable to reactivate inhibited human AChE under some conditions. Zwitterionic and cationic species have the best chance of productive action on inhibited AChE. However uncharged oximes can give important interaction information. In order to investigate the interaction and behavior of cationic and uncharged oximes, we performed molecular docking simulations and molecular dynamics and calculated binding energies of complexes of these compounds with human AChE. The uncharged oximes of larger structure were more susceptible to the influence of the substituents on the phosphorus atom and presented low binding energies. In contrast, HI 6 and 2-PAM showed high binding energy values with great contribution of the amino acid Asp74, demonstrating the importance of the quaternary nitrogen to the affinity and interaction of the oximes/AChE tabun-inhibited complexes.

Assessment of Detection and Refinement Strategies for de novo Protein Structures Using Force Field and Statistical Potentials

De novo predictions of protein structures at high resolution are plagued by the problem of detecting the native conformation from false energy minima. In this work, we provide an assessment of various detection and refinement protocols on a small subset of the second-generation all-atom Rosetta decoy set (Tsai et al. Proteins 2003, 53, 76-87) using two potentials:  the all-atom CHARMM PARAM22 force field combined with generalized Born/surface-area (GB-SA) implicit solvation and the DFIRE-AA statistical potential. Detection schemes included DFIRE-AA conformational scoring and energy minimization followed by scoring with both GB-SA and DFIRE-AA potentials. Refinement methods included short-time (1-ps) molecular dynamics simulations, temperature-based replica exchange molecular dynamics, and a new computational unfold/refold procedure. Refinement methods include temperature-based replica exchange molecular dynamics and a new computational unfold/refold procedure. Our results indicate that simple detection with only minimization is the best protocol for finding the most nativelike structures in the decoy set. The refinement techniques that we tested are generally unsuccessful in improving detection; however, they provide marginal improvements to some of the decoy structures. Future directions in the development of refinement techniques are discussed in the context of the limitations of the protocols evaluated in this study.

Selective and brain-permeable polo-like kinase-2 (Plk-2) inhibitors that reduce α-synuclein phosphorylation in rat brain

Polo-like kinase-2 (Plk-2) has been implicated as the dominant kinase involved in the phosphorylation of α-synuclein in Lewy bodies, which are one of the hallmarks of Parkinson's disease neuropathology. Potent, selective, brain-penetrant inhibitors of Plk-2 were obtained from a structure-guided drug discovery approach driven by the first reported Plk-2-inhibitor complexes. The best of these compounds showed excellent isoform and kinome-wide selectivity, with physicochemical properties sufficient to interrogate the role of Plk-2 inhibition in vivo. One such compound significantly decreased phosphorylation of α-synuclein in rat brain upon oral administration and represents a useful probe for future studies of this therapeutic avenue toward the potential treatment of Parkinson's disease.

Importance of polar solvation and configurational entropy for design of antiretroviral drugs targeting HIV-1 protease

Both KNI-10033 and KNI-10075 are high affinity preclinical HIV-1 protease (PR) inhibitors with affinities in the picomolar range. In this work, the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method has been used to investigate the potency of these two HIV-1 PR inhibitors against the wild-type and mutated proteases assuming that potency correlates with the affinity of the drugs for the target protein. The decomposition of the binding free energy reveals the origin of binding affinities or mutation-induced affinity changes. Our calculations indicate that the mutation I50V causes drug resistance against both inhibitors. On the other hand, we predict that the mutant I84V causes drug resistance against KNI-10075 while KNI-10033 is more potent against the I84V mutant compared to wild-type protease. Drug resistance arises mainly from unfavorable shifts in van der Waals interactions and configurational entropy. The latter indicates that neglecting changes in configurational entropy in the computation of relative binding affinities as often done is not appropriate in general. For the bound complex PR(I50V)-KNI-10075, an increased polar solvation free energy also contributes to the drug resistance. The importance of polar solvation free energies is revealed when interactions governing the binding of KNI-10033 or KNI-10075 to the wild-type protease are compared to the inhibitors darunavir or GRL-06579A. Although the contributions from intermolecular electrostatic and van der Waals interactions as well as the nonpolar component of the solvation free energy are more favorable for PR-KNI-10033 or PR-KNI-10075 compared to PR-DRV or PR-GRL-06579A, both KNI-10033 and KNI-10075 show a similar affinity as darunavir and a lower binding affinity relative to GRL-06579A. This is because of the polar solvation free energy which is less unfavorable for darunavir or GRL-06579A relative to KNI-10033 or KNI-10075. The importance of the polar solvation as revealed here highlights that structural inspection alone is not sufficient for identifying the key contributions to binding affinities and affinity changes for the design of drugs but that solvation effects must be taken into account. A detailed understanding of the molecular forces governing binding and drug resistance might assist in the design of new inhibitors against HIV-1 PR variants that are resistant against current drugs.

Mutation-induced loop opening and energetics for binding of tamiflu to influenza N8 neuraminidase

Tamiflu, also known as oseltamivir (OTV), binds to influenza A neuraminidase (H5N1) with very high affinity (0.32 nM). However, this inhibitor binds to other neuraminidases as well. In the present work, a systematic computational study is performed to investigate the mechanism underlying the binding of oseltamivir to N8 neuraminidase (NA) in "open" and "closed" conformations of the 150-loop through molecular dynamics simulations and the popular and well established molecular mechanics Poisson-Boltzmann (MM-PBSA) free energy calculation method. Whereas the closed conformation is stable for wild type N8, it transforms into the open conformation for the mutants Y252H, H274Y, and R292K, indicating that bound to oseltamivir these mutants are preferentially in the open conformation. Our calculations show that the binding of wild type oseltamivir to the closed conformation of N8 neuraminidase is energetically favored compared to the binding to the open conformation. We observe water mediated binding of oseltamivir to the N8 neuraminidase in both conformations which is not seen in the case of binding of the same drug to the H5N1 neuraminidase. The decomposition of the binding free energy reveals the mechanisms underlying the binding and changes in affinity due to mutations. Considering the mutant N8 variants in the open conformation adopted during the simulations, we observe a significant loss in the size of the total binding free energy for the N8(Y252H)-OTV, N8(H274Y)-OTV, and N8(R292K)-OTV complexes compared to N8(WT)-OTV, mainly due to the decrease in the size of the intermolecular electrostatic energy. For R292K, an unfavorable shift in the van der Waals interactions also contributes to the drug resistance. The mutations cause a significant expansion in the active site cavity, increasing its solvent accessible surface compared to the crystal structures of both the open and closed conformations. Our study underscores the need to consider dynamics in rationalizing the structure-function relationships of various antiviral inhibitor-NA complexes.

Potential of mean force of water-proton bath and molecular dynamic simulation of proteins at constant pH

An advanced implicit solvent model of water-proton bath for protein simulations at constant pH is presented. The implicit water-proton bath model approximates the potential of mean force of a protein in water solvent in a presence of hydrogen ions. Accurate and fast computational implementation of the implicit water-proton bath model is developed using the continuum electrostatic Poisson equation model for calculation of ionization equilibrium and the corrected MSR6 generalized Born model for calculation of the electrostatic atom-atom interactions and forces. Molecular dynamics (MD) method for protein simulation in the potential of mean force of water-proton bath is developed and tested on three proteins. The model allows to run MD simulations of proteins at constant pH, to calculate pH-dependent properties and free energies of protein conformations. The obtained results indicate that the developed implicit model of water-proton bath provides an efficient way to study thermodynamics of biomolecular systems as a function of pH, pH-dependent ionization-conformation coupling, and proton transfer events.

Molecular modeling of the interaction between heparan sulfate and cellular growth factors: bringing pieces together

Heparan sulfate is a polysaccharide belonging to the glycaminoglycan family. It interacts with numerous proteins of the extracellular matrix, in particular cellular growth factors. The number of experimental protein-heparin sulfate complexes obtained by crystallography or nuclear magnetic resonance is limited. Alternatively, computational approaches can be employed. Generally, they restrain the conformation of the glycosidic rings and linkages in order to reduce the complexity of the problem. Modeling the interaction between protein and heparan sulfate is indeed challenging because of the large size of the fragment needed for a strong binding, the flexibility brought by the glycosidic rings and linkages and the high density of negative charges. We propose a two-step method based on molecular docking and molecular dynamics simulation. Molecular docking allows exploring the positioning of a rigid heparin sulfate fragment on the protein surface. Molecular dynamics refine selected docking models by explicitly representing solvent molecules and not restraining the polysaccharide backbone. The interaction of a hexamer of heparin sulfate was studied in interaction with fibroblast growth factor 2 and stromal cell-derived factor 1α. This approach shed light on the plasticity of the growth factors interacting with heparan sulfate. This approach can be extended to the study of other protein/glycosaminoglycan complexes.

Advances in implicit models of water solvent to compute conformational free energy and molecular dynamics of proteins at constant pH

Modern implicit solvent models for macromolecular simulations in water-proton bath are considered. The fundamental quantity that implicit models approximate is the solute potential of mean force, which is obtained by averaging over solvent degrees of freedom. The implicit solvent models suggest practical ways to calculate free energies of macromolecular conformations taking into account equilibrium interactions with water solvent and proton bath, while the explicit solvent approach is unable to do that due to the need to account for a large number of solvent degrees of freedom. The most advanced realizations of the implicit continuum models by different research groups are discussed, their accuracy are examined, and some applications of the implicit solvent models to macromolecular modeling, such as free energy calculations, protein folding, and constant pH molecular dynamics are highlighted.

Blind docking method combining search of low-resolution binding sites with ligand pose refinement by molecular dynamics-based global optimization

This study describes the development of a new blind hierarchical docking method, bhDock, its implementation, and accuracy assessment. The bhDock method uses two-step algorithm. First, a comprehensive set of low-resolution binding sites is determined by analyzing entire protein surface and ranked by a simple score function. Second, ligand position is determined via a molecular dynamics-based method of global optimization starting from a small set of high ranked low-resolution binding sites. The refinement of the ligand binding pose starts from uniformly distributed multiple initial ligand orientations and uses simulated annealing molecular dynamics coupled with guided force-field deformation of protein-ligand interactions to find the global minimum. Assessment of the bhDock method on the set of 37 protein-ligand complexes has shown the success rate of predictions of 78%, which is better than the rate reported for the most cited docking methods, such as AutoDock, DOCK, GOLD, and FlexX, on the same set of complexes.

Explicitly solvated ligand contribution to continuum solvation models for binding free energies: selectivity of theophylline binding to an RNA aptamer

To provide more accurate computational estimates of binding free energies in solution from molecular dynamics (MD) simulations, a separate solvation contribution for the binding ligand is determined from a linear response treatment. We use explicit water coordinates for this term and combine with MM-PBSA (molecular mechanics, Poisson-Boltzmann, and surface area contributions) in a new approach (MM-PB/LRA-SA). To assess this method, application to the binding between theophylline and its derivatives to an RNA aptamer was performed and compared with experimental binding affinities. Explicitly solvated MD trajectories were generated with the same parameter set used in the previous work by Gouda et al., who compared the relative binding of these molecules by both the MM-PBSA and thermodynamic integration methods. Substituting the linear response term for the ligand in the MM-PB/LRA-SA approach led to an improvement upon MM-PBSA when compared with experimental and thermodynamic integration results at approximately twice the computational cost. The balance between accuracy and computational expense achieved using this method suggests potential advantages in applying it in the virtual drug-screening process.

Identifying native-like protein structures with scoring functions based on all-atom ECEPP force fields, implicit solvent models and structure relaxation

Availability of energy functions which can discriminate native-like from non-native protein conformations is crucial for theoretical protein structure prediction and refinement of low-resolution protein models. This article reports the results of benchmark tests for scoring functions based on two all-atom ECEPP force fields, that is, ECEPP/3 and ECEPP05, and two implicit solvent models for a large set of protein decoys. The following three scoring functions are considered: (i) ECEPP05 plus a solvent-accessible surface area model with the parameters optimized with a set of protein decoys (ECEPP05/SA); (ii) ECEPP/3 plus the solvent-accessible surface area model of Ooi et al. (Proc Natl Acad Sci USA 1987;84:3086-3090) (ECEPP3/OONS); and (iii) ECEPP05 plus an implicit solvent model based on a solution of the Poisson equation with an optimized Fast Adaptive Multigrid Boundary Element (FAMBEpH) method (ECEPP05/FAMBEpH). Short Monte Carlo-with-Minimization (MCM) simulations, following local energy minimization, are used as a scoring method with ECEPP05/SA and ECEPP3/OONS potentials, whereas energy calculation is used with ECEPP05/FAMBEpH. The performance of each scoring function is evaluated by examining its ability to distinguish between native-like and non-native protein structures. The results of the tests show that the new ECEPP05/SA scoring function represents a significant improvement over the earlier ECEPP3/OONS version of the force field. Thus, it is able to rank native-like structures with C(alpha) root-mean-square-deviations below 3.5 A as lowest-energy conformations for 76% and within the top 10 for 87% of the proteins tested, compared with 69 and 80%, respectively, for ECEPP3/OONS. The use of the FAMBEpH solvation model, which provides a more accurate description of the protein-solvent interactions, improves the discriminative ability of the scoring function to 89%. All failed tests in which the native-like structures cannot be discriminated as those with low energy, are due to omission of protein-protein interactions. The results of this study represent a benchmark in force-field development, and may be useful for evaluation of the performance of different force fields.

Molecular basis for Bcl-2 homology 3 domain recognition in the Bcl-2 protein family: identification of conserved hot spot interactions

The proteins of the Bcl-2 family are important regulators of apoptosis, or programmed cell death. These proteins regulate this fundamental biological process via the formation of heterodimers involving both pro- and anti-apoptotic family members. Disruption of the balance between anti- and pro-apoptotic Bcl-2 proteins is the cause of numerous pathologies. Bcl-xl, an anti-apoptotic protein of this family, is known to form heterodimers with multiple pro-apoptotic proteins, such as Bad, Bim, Bak, and Bid. To elucidate the molecular basis of this recognition process, we used molecular dynamics simulations coupled with the Molecular Mechanics/Poisson-Boltzmann Surface Area approach to identify the amino acids that make significant energetic contributions to the binding free energy of four complexes formed between Bcl-xl and pro-apoptotic Bcl-2 homology 3 peptides. A fifth protein-peptide complex composed of another anti-apoptotic protein, Bcl-w, in complex with the peptide from Bim was also studied. The results identified amino acids of both the anti-apoptotic proteins as well as the Bcl-2 homology 3 (BH3) domains of the pro-apoptotic proteins that make strong, recurrent interactions in the protein complexes. The calculations show that the two anti-apoptotic proteins, Bcl-xl and Bcl-w, share a similar recognition mechanism. Our results provide insight into the molecular basis for the promiscuous nature of this molecular recognition process by members of the Bcl-2 protein family. These amino acids could be targeted in the design of new mimetics that serve as scaffolds for new antitumoral molecules.

Structure and folding of disulfide-rich miniproteins: insights from molecular dynamics simulations and MM-PBSA free energy calculations

The fold of small disulfide-rich proteins largely relies on two or more disulfide bridges that are main components of the hydrophobic core. Because of the small size of these proteins and their high cystine content, the cysteine connectivity has been difficult to ascertain in some cases, leading to uncertainties and debates in the literature. Here, we use molecular dynamics simulations and MM-PBSA free energy calculations to compare similar folds with different disulfide pairings in two disulfide-rich miniprotein families, namely the knottins and the short-chain scorpion toxins, for which the connectivity has been discussed. We first show that the MM-PBSA approach is able to discriminate the correct knotted topology of knottins from the laddered one. Interestingly, a comparison of the free energy components for kalata B1 and MCoTI-II suggests that cyclotides and squash inhibitors, although sharing the same scaffold, are stabilized through different interactions. Application to short-chain scorpion toxins suggests that the conventional cysteine pairing found in many homologous toxins is significantly more stable than the unconventional pairing reported for maurotoxin and for spinoxin. This would mean that native maurotoxin and spinoxin are not at the lowest free energy minimum and might result from kinetically rather than thermodynamically driven oxidative folding processes. For both knottins and toxins, the correct or conventional disulfide connectivities provide lower flexibilities and smaller deviations from the initial conformations. Overall, our work suggests that molecular dynamics simulations and the MM-PBSA approach to estimate free energies are useful tools to analyze and compare disulfide bridge connectivities in miniproteins.

Molecular modeling and characterization of the B. thuringiensis and B. thuringiensis LDC-9 cytolytic proteins

The Cyt toxins are able to lyse a wide range of cell types in vitro, unlike the Cry delta-endotoxins. It exerts its activity by the formation of pores within target cell membranes. The structural information available for Cyt2Aa (PDB id: 1CBY) consists of a single domain in which two outer layers of alpha-helix wrap around a mixed beta-sheet. Beta-barrel was suggested as a possible structure of the pores. Hence, this study seeks to investigate the structural properties of other Cytolytic proteins by predicting the three-dimensional (3D) model using Cyt2Aa as template. The predicted models are expected to be significantly more accurate as all the Cyt proteins showed significant similarity with the template (PDB id: 1CBY). The refined homology models revealed similar secondary structures (alpha-helices and beta-sheets) and tertiary features as Cyt2Aa. The variation in the loop regions of the tertiary structure accounts for the differential toxicity.

FAMBE-pH: a fast and accurate method to compute the total solvation free energies of proteins

A fast and accurate method to compute the total solvation free energies of proteins as a function of pH is presented. The method makes use of a combination of approaches, some of which have already appeared in the literature; (i) the Poisson equation is solved with an optimized fast adaptive multigrid boundary element (FAMBE) method; (ii) the electrostatic free energies of the ionizable sites are calculated for their neutral and charged states by using a detailed model of atomic charges; (iii) a set of optimal atomic radii is used to define a precise dielectric surface interface; (iv) a multilevel adaptive tessellation of this dielectric surface interface is achieved by using multisized boundary elements; and (v) 1:1 salt effects are included. The equilibrium proton binding/release is calculated with the Tanford-Schellman integral if the proteins contain more than approximately 20-25 ionizable groups; for a smaller number of ionizable groups, the ionization partition function is calculated directly. The FAMBE method is tested as a function of pH (FAMBE-pH) with three proteins, namely, bovine pancreatic trypsin inhibitor (BPTI), hen egg white lysozyme (HEWL), and bovine pancreatic ribonuclease A (RNaseA). The results are (a) the FAMBE-pH method reproduces the observed pK a's of the ionizable groups of these proteins within an average absolute value of 0.4 p K units and a maximum error of 1.2 p K units and (b) comparison of the calculated total pH-dependent solvation free energy for BPTI, between the exact calculation of the ionization partition function and the Tanford-Schellman integral method, shows agreement within 1.2 kcal/mol. These results indicate that calculation of total solvation free energies with the FAMBE-pH method can provide an accurate prediction of protein conformational stability at a given fixed pH and, if coupled with molecular mechanics or molecular dynamics methods, can also be used for more realistic studies of protein folding, unfolding, and dynamics, as a function of pH.

Use of decoys to optimize an all-atom force field including hydration

A novel method of parameter optimization is proposed. It makes use of large sets of decoys generated for six nonhomologous proteins with different architecture. Parameter optimization is achieved by creating a free energy gap between sets of nativelike and nonnative conformations. The method is applied to optimize the parameters of a physics-based scoring function consisting of the all-atom ECEPP05 force field coupled with an implicit solvent model (a solvent-accessible surface area model). The optimized force field is able to discriminate near-native from nonnative conformations of the six training proteins when used either for local energy minimization or for short Monte Carlo simulated annealing runs after local energy minimization. The resulting force field is validated with an independent set of six nonhomologous proteins, and appears to be transferable to proteins not included in the optimization; i.e., for five out of the six test proteins, decoys with 1.7- to 4.0-A all-heavy-atom root mean-square deviations emerge as those with the lowest energy. In addition, we examined the set of misfolded structures created by Park and Levitt using a four-state reduced model. The results from these additional calculations confirm the good discriminative ability of the optimized force field obtained with our decoy sets.

Can molecular dynamics simulations help in discriminating correct from erroneous protein 3D models?

BACKGROUND

Recent approaches for predicting the three-dimensional (3D) structure of proteins such as de novo or fold recognition methods mostly rely on simplified energy potential functions and a reduced representation of the polypeptide chain. These simplifications facilitate the exploration of the protein conformational space but do not permit to capture entirely the subtle relationship that exists between the amino acid sequence and its native structure. It has been proposed that physics-based energy functions together with techniques for sampling the conformational space, e.g., Monte Carlo or molecular dynamics (MD) simulations, are better suited to the task of modelling proteins at higher resolutions than those of models obtained with the former type of methods. In this study we monitor different protein structural properties along MD trajectories to discriminate correct from erroneous models. These models are based on the sequence-structure alignments provided by our fold recognition method, FROST. We define correct models as being built from alignments of sequences with structures similar to their native structures and erroneous models from alignments of sequences with structures unrelated to their native structures.

RESULTS

For three test sequences whose native structures belong to the all-alpha, all-beta and alphabeta classes we built a set of models intended to cover the whole spectrum: from a perfect model, i.e., the native structure, to a very poor model, i.e., a random alignment of the test sequence with a structure belonging to another structural class, including several intermediate models based on fold recognition alignments. We submitted these models to 11 ns of MD simulations at three different temperatures. We monitored along the corresponding trajectories the mean of the Root-Mean-Square deviations (RMSd) with respect to the initial conformation, the RMSd fluctuations, the number of conformation clusters, the evolution of secondary structures and the surface area of residues. None of these criteria alone is 100% efficient in discriminating correct from erroneous models. The mean RMSd, RMSd fluctuations, secondary structure and clustering of conformations show some false positives whereas the residue surface area criterion shows false negatives. However if we consider these criteria in combination it is straightforward to discriminate the two types of models.

CONCLUSION

The ability of discriminating correct from erroneous models allows us to improve the specificity and sensitivity of our fold recognition method for a number of ambiguous cases.

Can a physics-based, all-atom potential find a protein's native structure among misfolded structures? I. Large scale AMBER benchmarking

Recent work has shown that physics-based, all-atom energy functions (AMBER, CHARMM, OPLS-AA) and local minimization, when used in scoring, are able to discriminate among native and decoy structures. Yet, there have been only few instances reported of the successful use of physics based potentials in the actual refinement of protein models from a starting conformation to one that ends in structures, which are closer to the native state. An energy function that has a global minimum energy in the protein's native state and a good correlation between energy and native-likeness should be able to drive model structures closer to their native structure during a conformational search. Here, the possible reasons for the discrepancy between the scoring and refinement results for the case of AMBER potential are examined. When the conformational search via molecular dynamics is driven by the AMBER potential for a large set of 150 nonhomologous proteins and their associated decoys, often the native minimum does not appear to be the lowest free energy state. Ways of correcting the potential function in order to make it more suitable for protein model refinement are proposed.

Hydrophobic potential of mean force as a solvation function for protein structure prediction

We have developed a solvation function that combines a Generalized Born model for polarization of protein charge by the high dielectric solvent, with a hydrophobic potential of mean force (HPMF) as a model for hydrophobic interaction, to aid in the discrimination of native structures from other misfolded states in protein structure prediction. We find that our energy function outperforms other reported scoring functions in terms of correct native ranking for 91% of proteins and low Z scores for a variety of decoy sets, including the challenging Rosetta decoys. This work shows that the stabilizing effect of hydrophobic exposure to aqueous solvent that defines the HPMF hydration physics is an apparent improvement over solvent-accessible surface area models that penalize hydrophobic exposure. Decoys generated by thermal sampling around the native-state basin reveal a potentially important role for side-chain entropy in the future development of even more accurate free energy surfaces.

Implicit solvent models

Implicit solvent models for biomolecular simulations are reviewed and their underlying statistical mechanical basis is discussed. The fundamental quantity that implicit models seek to approximate is the solute potential of mean force, which determines the statistical weight of solute conformations, and which is obtained by averaging over the solvent degrees of freedom. It is possible to express the total free energy as the reversible work performed in two successive steps. First, the solute is inserted in the solvent with zero atomic partial charges; second, the atomic partial charges of the solute are switched from zero to their full values. Consequently, the total solvation free energy corresponds to a sum of non-polar and electrostatic contributions. These two contributions are often approximated by simple geometrical models (such as solvent exposed area models) and by macroscopic continuum electrostatics, respectively. One powerful route is to approximate the average solvent density distribution around the solute, i.e. the solute-solvent density correlation functions, as in statistical mechanical integral equations. Recent progress with semi-analytical approximations makes continuum electrostatics treatments very efficient. Still more efficient are fully empirical, knowledge-based models, whose relation to explicit solvent treatments is not fully resolved, however. Continuum models that treat both solute and solvent as dielectric continua are also discussed, and the relation between the solute fluctuations and its macroscopic dielectric constant(s) clarified.

Towards the high-resolution protein structure prediction. Fast refinement of reduced models with all-atom force field

Background

Although experimental methods for determining protein structure are providing high resolution structures, they cannot keep the pace at which amino acid sequences are resolved on the scale of entire genomes. For a considerable fraction of proteins whose structures will not be determined experimentally, computational methods can provide valuable information. The value of structural models in biological research depends critically on their quality. Development of high-accuracy computational methods that reliably generate near-experimental quality structural models is an important, unsolved problem in the protein structure modeling.

Results

Large sets of structural decoys have been generated using reduced conformational space protein modeling tool CABS. Subsequently, the reduced models were subject to all-atom reconstruction. Then, the resulting detailed models were energy-minimized using state-of-the-art all-atom force field, assuming fixed positions of the alpha carbons. It has been shown that a very short minimization leads to the proper ranking of the quality of the models (distance from the native structure), when the all-atom energy is used as the ranking criterion. Additionally, we performed test on medium and low accuracy decoys built via classical methods of comparative modeling. The test placed our model evaluation procedure among the state-of-the-art protein model assessment methods.

Conclusion

These test computations show that a large scale high resolution protein structure prediction is possible, not only for small but also for large protein domains, and that it should be based on a hierarchical approach to the modeling protocol. We employed Molecular Mechanics with fixed alpha carbons to rank-order the all-atom models built on the scaffolds of the reduced models. Our tests show that a physic-based approach, usually considered computationally too demanding for large-scale applications, can be effectively used in such studies.

Protein-protein recognition and interaction hot spots in an antigen-antibody complex: free energy decomposition identifies "efficient amino acids"

The molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) method was applied to the study of the protein-protein complex between a camelid single chain variable domain (cAb-Lys3) and hen egg white lysozyme (HEL), and between cAb-Lys3 and turkey egg white lysozyme (TEL). The electrostatic energy was estimated by solving the linear Poisson-Boltzmann equation. A free energy decomposition scheme was developed to determine binding energy hot spots of each complex. The calculations identified amino acids of the antibody that make important contributions to the interaction with lysozyme. They further showed the influence of small structural variations on the energetics of binding and they showed that the antibody amino acids that make up the hot spots are organized in such a way as to mimic the lysozyme substrate. Through further analysis of the results, we define the concept of "efficient amino acids," which can provide an assessment of the binding potential of a particular hot spot interaction. This information, in turn, can be useful in the rational design of small molecules that mimic the antibody. The implications of using free energy decomposition to identify regions of a protein-protein complex that could be targeted by small molecules inhibitors are discussed.

Structural refinement of protein segments containing secondary structure elements: Local sampling, knowledge-based potentials, and clustering

In this article, we present an iterative, modular optimization (IMO) protocol for the local structure refinement of protein segments containing secondary structure elements (SSEs). The protocol is based on three modules: a torsion-space local sampling algorithm, a knowledge-based potential, and a conformational clustering algorithm. Alternative methods are tested for each module in the protocol. For each segment, random initial conformations were constructed by perturbing the native dihedral angles of loops (and SSEs) of the segment to be refined while keeping the protein body fixed. Two refinement procedures based on molecular mechanics force fields - using either energy minimization or molecular dynamics - were also tested but were found to be less successful than the IMO protocol. We found that DFIRE is a particularly effective knowledge-based potential and that clustering algorithms that are biased by the DFIRE energies improve the overall results. Results were further improved by adding an energy minimization step to the conformations generated with the IMO procedure, suggesting that hybrid strategies that combine both knowledge-based and physical effective energy functions may prove to be particularly effective in future applications.

FALS mutations in Cu, Zn superoxide dismutase destabilize the dimer and increase dimer dissociation propensity: a large-scale thermodynamic analysis

Mutations in the dimeric enzyme Cu, Zn superoxide dismutase (SOD1) leading to its aggregation are implicated in the toxicity in familial amyotrophic lateral sclerosis (FALS). We and others have previously shown that aggregation occurs by a pathway involving dimer dissociation, metal-loss from monomers and multimeric assembly of apo-SOD1 monomers. We postulate that FALS mutations cause enhanced aggregation by affecting one or more steps in the pathway, and computationally test this postulate for 75 known mis-sense FALS mutants of SOD1. Based on an extensive thermodynamic analysis of the stability of apo-dimer and apo-monomer forms of these mutants, we classify the mutations into the following groups: 70 out of 75 mutations in SOD1 lead to (i) decreased dimer stability, and/or (ii) increased dimer dissociation, compared to wild type, and four mutations lead to (iii) decreased monomer stability compared to wild type. Our results suggest that enhanced aggregation of SOD1 in FALS occurs due to an increased population of mutant SOD1 apo-monomers compared to wild type. The dissociation of multimeric proteins induced by diverse mutations may be a common theme in several human diseases.

Conformational sampling with implicit solvent models: application to the PHF6 peptide in tau protein

Implicit solvent models approximate the effects of solvent through a potential of mean force and therefore make solvated simulations computationally efficient. Yet despite their computational efficiency, the inherent approximations made by implicit solvent models can sometimes lead to inaccurate results. To test the accuracy of a number of popular implicit solvent models, we determined whether implicit solvent simulations can reproduce the set of potential energy minima obtained from explicit solvent simulations. For these studies, we focus on a six-residue amino-acid sequence, referred to as the paired helical filament 6 (PHF6), which may play an important role in the formation of intracellular aggregates in patients with Alzheimer's disease. Several implicit solvent models form the basis of this work--two based on the generalized Born formalism, and one based on a Gaussian solvent-exclusion model. All three implicit solvent models generate minima that are in good agreement with minima obtained from simulations with explicit solvent. Moreover, free-energy profiles generated with each implicit solvent model agree with free-energy profiles obtained with explicit solvent. For the Gaussian solvent-exclusion model, we demonstrate that a straightforward ranking of the relative stability of each minimum suggests that the most stable structure is extended, a result in excellent agreement with the free-energy profiles. Overall, our data demonstrate that for some peptides like PHF6, implicit solvent can accurately reproduce the set of local energy minimum arising from quenched dynamics simulations with explicit solvent. More importantly, all solvent models predict that PHF6 forms extended beta-structures in solution, a finding consistent with the notion that PHF6 initiates neurofibrillary tangle formation in patients with Alzheimer's disease.

An Efficient, Path-Independent Method for Free-Energy Calculations

Classical free-energy methods depend on the definition of physical or nonphysical integration paths to calculate free-energy differences between states. This procedure can be problematic and computationally expensive when the states of interest do not overlap and are far apart in phase space. Here we introduce a novel method to calculate free-energy differences that is path-independent by transforming each end state into a reference state in which the vibrational entropy is the sole component of the total entropy, thus allowing direct computation of the relative free energy. We apply the method to calculate side-chain entropies of a beta-hairpin-forming peptide in a variety of backbone conformations, demonstrating its importance in determining structural propensities. We find that low-free-energy conformations achieve their stability through optimal trade off between enthalpic gains due to favorable interatomic interactions and entropic losses incurred by the same.

Reducing CDK4/6-p16(INK4a) interface. Computational alanine scanning of a peptide bound to CDK6 protein

The tumor suppressor gene p16INK4a is commonly found altered in numerous and different types of cancer. The encoded protein arrests cell cycle in G1 phase by binding to CDK4 and CDK6, inhibiting their kinase function. In 1995, a 20-residue peptide, extracted from p16INK4a protein sequence, was discovered that retains the cell cycle inhibition properties of the endogenous tumor suppressor. However, its structure has not been determined yet. In this article, the features of a theoretical structure of the peptide bound to CDK6 are reported. The complex was modeled from CDK6-p16INK4a X-ray crystal structure and through molecular dynamics. Final structure was assessed by comparing computed binding free energy changes, when single-alanine substitutions were brought about on the peptide, to experimental data. Better concordance was obtained when including a high level of solvation effects. Solute-solvent vdW energy and electrostatic energy between solute and first shells of water, computed through a force field and considering explicit waters, were also to be included to achieve reasonably good concordance between theoretical and experimental data.