Calculation of protein-ligand binding entropies using a rule-based molecular fingerprint.

The use of fast in silico prediction methods for protein-ligand binding free energies holds significant promise for the initial phases of drug development. Numerous traditional physics-based models (e.g., implicit solvent models), however, tend to either neglect or heavily approximate entropic contributions to binding due to their computational complexity. Consequently, such methods often yield imprecise assessments of binding strength. Machine learning models provide accurate predictions and can often outperform physics-based models. They, however, are often prone to overfitting, and the interpretation of their results can be difficult. Physics-guided machine learning models combine the consistency of physics-based models with the accuracy of modern data-driven algorithms. This work integrates physics-based model conformational entropies into a graph convolutional network. We introduce a new neural network architecture (a rule-based graph convolutional network) that generates molecular fingerprints according to predefined rules specifically optimized for binding free energy calculations. Our results on 100 small host-guest systems demonstrate significant improvements in convergence and preventing overfitting. We additionally demonstrate the transferability of our proposed hybrid model by training it on the aforementioned host-guest systems and then testing it on six unrelated protein-ligand systems. Our new model shows little difference in training set accuracy compared to a previous model but an order-of-magnitude improvement in test set accuracy. Finally, we show how the results of our hybrid model can be interpreted in a straightforward fashion.

Unraveling Complex Local Protein Environments with 4-Cyano-l-phenylalanine.

We present a multifaceted approach to effectively probe complex local protein environments utilizing the vibrational reporter unnatural amino acid (UAA) 4-cyano-l-phenylalanine (pCNPhe) in the model system superfolder green fluorescent protein (sfGFP). This approach combines temperature-dependent infrared (IR) spectroscopy, X-ray crystallography, and molecular dynamics (MD) simulations to provide a molecular interpretation of the local environment of the nitrile group in the protein. Specifically, a two-step enantioselective synthesis was developed that provided an 87% overall yield of pCNPhe in high purity without the need for chromatography. It was then genetically incorporated individually at three unique sites (74, 133, and 149) in sfGFP to probe these local protein environments. The incorporation of the UAA site-specifically in sfGFP utilized an engineered, orthogonal tRNA synthetase in E. coli using the Amber codon suppression protocol, and the resulting UAA-containing sfGFP constructs were then explored with this approach. This methodology was effectively utilized to further probe the local environments of two surface sites (sites 133 and 149) that we previously explored with room temperature IR spectroscopy and X-ray crystallography and a new interior site (site 74) featuring a complex local environment around the nitrile group of pCNPhe. Site 133 was found to be solvent-exposed, while site 149 was partially buried. Site 74 was found to consist of three distinct local environments around the nitrile group including nonspecific van der Waals interactions, hydrogen-bonding to a structural water, and hydrogen-bonding to a histidine side chain.

Development and Validation of AMBER-FB15-Compatible Force Field Parameters for Phosphorylated Amino Acids.

Phosphorylation of select amino acid residues is one of the most common biological mechanisms for regulating protein structures and functions. While computational modeling can be used to explore the detailed structural changes associated with phosphorylation, most molecular mechanics force fields developed for the simulation of phosphoproteins have been noted to be inconsistent with experimental data. In this work, we parameterize force fields for the phosphorylated forms of the amino acids serine, threonine, and tyrosine using the ForceBalance software package with the goal of improving agreement with experiments for these residues. Our optimized force field, denoted as FB18, is parameterized using high-quality ab initio potential energy scans and is designed to be fully compatible with the AMBER-FB15 protein force field. When utilized in MD simulations together with the TIP3P-FB water model, we find that FB18 consistently enhances the prediction of experimental quantities such as 3J NMR couplings and intramolecular hydrogen-bonding propensities in comparison to previously published models. As was reported with AMBER-FB15, we also see improved agreement with the reference QM calculations in regions at and away from local minima. We thus believe that the FB18 parameter set provides a promising route for the further investigation of the varied effects of protein phosphorylation.

Phosphoserine inhibits neighboring arginine methylation in the RKS motif of histone H3.

The effects of phosphorylation of histone H3 at serine 10 have been studied in the context of other posttranslational modifications such as lysine methylation. We set out to investigate the impact of phosphoserine-10 on arginine-8 methylation. We performed methylation reactions using peptides based on histone H3 that contain a phosphorylated serine and compared the extent of arginine methylation with unmodified peptides. Results obtained via fluorography indicate that peptides containing a phosphorylated serine-10 inhibit deposition of methyl groups to arginine-8 residues. To further explore the effects of phosphoserine on neighboring arginine residues, we physically characterized the non-covalent interactions between histone H3 phosphoserine-10 and arginine-8 using 31P NMR spectroscopy. A salt bridge was detected between the negatively charged phosphoserine-10 and the positively charged unmodified arginine-8 residue. This salt bridge was not detected when arginine-8 was symmetrically dimethylated. Finally, molecular simulations not only confirm the presence of a salt bridge but also identify a subset of electrostatic interactions present when arginine is replaced with alanine. Taken together, our work suggests that the negatively charged phosphoserine maximizes its interactions. By limiting its exposure and creating new contacts with neighboring residues, it will inhibit deposition of neighboring methyl groups, not through steric hindrance, but by forming intrapeptide interactions that may mask substrate recognition. Our work provides a mechanistic framework for understanding the role of phosphoserine on nearby amino acid residues and arginine methylation.

Non-bonded force field model with advanced restrained electrostatic potential charges (RESP2).

The restrained electrostatic potential (RESP) approach is a highly regarded and widely used method of assigning partial charges to molecules for simulations. RESP uses a quantum-mechanical method that yields fortuitous overpolarization and thereby accounts only approximately for self-polarization of molecules in the condensed phase. Here we present RESP2, a next generation of this approach, where the polarity of the charges is tuned by a parameter, δ, which scales the contributions from gas- and aqueous-phase calculations. When the complete non-bonded force field model, including Lennard-Jones parameters, is optimized to liquid properties, improved accuracy is achieved, even with this reduced set of five Lennard-Jones types. We argue that RESP2 with δ≈0.6 (60% aqueous, 40% gas-phase charges) is an accurate and robust method of generating partial charges, and that a small set of Lennard-Jones types is good starting point for a systematic re-optimization of this important non-bonded term.

Benchmarking Electronic Structure Methods for Accurate Fixed-Charge Electrostatic Models.

The accuracy of classical molecular mechanics (MM) force fields used for condensed phase molecular simulations depends strongly on the accuracy of modeling nonbonded interactions between atoms, such as electrostatic interactions. Some popular fixed-charge MM force fields use partial atomic charges derived from gas phase electronic structure calculations using the Hartree-Fock (HF) method with the relatively small 6-31G* basis set (HF/6-31G*). It is generally believed that HF/6-31G* generates fortuitously overpolarized electron distributions, as would be expected in the higher dielectric environment of the condensed phase. Using a benchmark set of 47 molecules, we show that HF/6-31G* overpolarizes molecules by just under 10% on average with respect to experimental gas phase dipole moments. The overpolarization of this method/basis set combination varies significantly though and, in some cases, even leads to molecular dipole moments that are lower than experimental gas phase measurements. We further demonstrate that using computationally inexpensive density functional theory (DFT) methods, together with appropriate augmented basis sets and a continuum solvent model, can yield molecular dipole moments that are both more strongly and more uniformly overpolarized. These data suggest that these methods-or ones similar to them-should be adopted for the derivation of accurate partial atomic charges for next-generation MM force fields.

Systematic Optimization of Water Models Using Liquid/Vapor Surface Tension Data.

In this work, we investigate whether experimental surface tension measurements, which are less sensitive to quantum and self-polarization corrections, are able to replace the usual reliance on the heat of vaporization as experimental reference data for fitting force field models of molecular liquids. To test this hypothesis, we develop the fitting protocol necessary to utilize surface tension measurements in the ForceBalance optimization procedure to determine revised parameters for both three-point and four-point water models TIP3P-ST and TIP4P-ST. We find that the incorporation of surface tension in the fit results in a rigid three-point model that reproduces the correct temperature of maximum density of water for the first time but also leads to overstructuring of the liquid and less accurate transport properties. The rigid four-point TIP4P-ST model is highly accurate for a broad range of thermodynamic and kinetic properties, with similar performance compared to recently developed four-point water models. The results show surface tension to be a useful fitting property in general, especially when self-polarization corrections or nuclear quantum corrections are not readily available for correcting the heat of vaporization as is the case for other molecular liquids.

New developments in force fields for biomolecular simulations.

Biomolecular force field development has been instrumental in improving the predictive power of molecular simulations over the past four decades. More recently, the era of large quantitative experimental datasets and ubiquitous high performance computing power has enabled rapid progress in the field. In this review we summarize recent developments in all-atom protein, nucleic acid, and small molecule force fields, paying specific attention to developments in parameterization methods and improvements in the representations of nonbonded interactions that are critical for solving the challenging biophysical problems of the present. We also sketch out new avenues for force field development and grand challenge applications for the near future.

Disordered structural ensembles of vasopressin and oxytocin and their mutants.

Vasopressin and oxytocin are intrinsically disordered cyclic nonapeptides belonging to a family of neurohypophysial hormones. Although unique in their functions, these peptides differ only by two residues and both feature a tocin ring formed by the disulfide bridge between first and sixth cysteine residues. This sequence and structural similarity are experimentally linked to oxytocin agonism at vasopressin receptors and vasopressin antagonism at oxytocin receptors. Yet single- or double-residue mutations in both peptides have been shown to have drastic impacts on their activities at either receptor, and possibly the ability to bind to their neurophysin carrier protein. In this study we perform molecular dynamics simulations of the unbound native and mutant sequences of the oxytocin and vasopressin hormones to characterize their structural ensembles. We classify the subpopulations of these structural ensembles on the basis of the distributions of radius of gyration and secondary structure and hydrogen-bonding features of the canonical tocin ring and disordered tail region. We then relate the structural changes observed in the unbound form of the different hormone sequences to experimental information about peptide receptor binding, and more indirectly, carrier protein binding affinity, receptor activity, and protease degradation. This study supports the hypothesis that the structural characteristics of the unbound form of an IDP can be used to predict structural or functional preferences of its functional bound form.

Optimizing Protein-Protein van der Waals Interactions for the AMBER ff9x/ff12 Force Field.

The quality of molecular dynamics (MD) simulations relies heavily on the accuracy of the underlying force field. In recent years, considerable effort has been put into developing more accurate dihedral angle potentials for MD force fields, but relatively little work has focused on the nonbonded parameters, many of which are two decades old. In this work, we assess the accuracy of protein-protein van der Waals interactions in the AMBER ff9x/ff12 force field. Across a test set of 44 neat organic liquids containing the moieties present in proteins, we find root-mean-square (RMS) errors of 1.26 kcal/mol in enthalpy of vaporization and 0.36 g/cm(3) in liquid densities. We then optimize the van der Waals radii and well depths for all of the relevant atom types using these observables, which lowers the RMS errors in enthalpy of vaporization and liquid density of our validation set to 0.59 kcal/mol (53% reduction) and 0.019 g/cm(3) (46% reduction), respectively. Limitations in our parameter optimization were evident for certain atom types, however, and we discuss the implications of these observations for future force field development.

Embedding Aß42 in heterogeneous membranes depends on cholesterol asymmetries.

Using a coarse-grained lipid and peptide model, we show that the free energy stabilization of amyloid-ß in heterogeneous lipid membranes is predicted to have a dependence on asymmetric distributions of cholesterol compositions across the membrane leaflets. We find that a highly asymmetric cholesterol distribution that is depleted on the exofacial leaflet but enhanced on the cytofacial leaflet of the model lipid membrane thermodynamically favors membrane retention of a fully embedded Aß peptide. However, in the case of cholesterol redistribution that increases concentration of cholesterol on the exofacial layer, typical of aging or Alzheimer's disease, the free energy favors peptide extrusion of the highly reactive N-terminus into the extracellular space that may be vulnerable to aggregation, oligomerization, or deleterious oxidative reactivity.

Optimizing solute-water van der Waals interactions to reproduce solvation free energies.

An accurate representation of solute-water interactions is necessary for molecular dynamics simulations of biomolecules that reside in aqueous environments. Modern force fields and advanced water models describe solute-solute and water-water interactions reasonably accurately but have known shortcomings in describing solute-water interactions, demonstrated by the large differences between calculated and experimental solvation free energies across a range of peptide and drug chemistries. In this work, we introduce a method for optimizing solute-water van der Waals interactions to reproduce experimental solvation free energy data and apply it to the optimization of a fixed charge force field (AMBER ff99SB/GAFF) and advanced water model (TIP4P-Ew). We show that, with these optimizations, the combination of AMBER ff99SB/GAFF and TIP4P-Ew is able to reproduce the solvation free energies of a variety of biologically relevant small molecules to within 1.0 k(B)T. We further validate these optimizations by examining the aggregation propensities of dipeptide-water solutions, the conformational preferences of short disordered peptides, and the native state stability and dynamics of a folded protein.

Homogeneous and heterogeneous tertiary structure ensembles of amyloid-ß peptides.

The interplay of modern molecular simulation and high-quality nuclear magnetic resonance (NMR) experiments has reached a fruitful stage for quantitative characterization of structural ensembles of disordered peptides. Amyloid-ß 1-42 (Aß42), the primary peptide associated with Alzheimer's disease, and fragments such as Aß21-30 are both classified as intrinsically disordered peptides (IDPs). We use a variety of NMR observables to validate de novo molecular dynamics simulations in explicit water to characterize the tertiary structure ensemble of Aß42 and Aß21-30 from the perspective of their classification as IDPs. Unlike the Aß21-30 fragment that conforms to expectations of an IDP that is primarily extended, we find that Aß42 samples conformations reflecting all possible secondary structure categories and spans the range of IDP classifications from collapsed structured states to highly extended conformations, making it an IDP with a far more heterogeneous tertiary ensemble.

Optimizing Protein-Solvent Force Fields to Reproduce Intrinsic Conformational Preferences of Model Peptides.

While most force field efforts in biomolecular simulation have focused on the parametrization of the protein, relatively little attention has been paid to the quality of the accompanying solvent model. These considerations are especially relevant for simulations of intrinsically disordered peptides and proteins, for which energy differences between conformations are small and interactions with water are enhanced. In this work, we investigate the accuracy of the AMBER ff99SB force field when combined with the standard TIP3P model or the more recent TIP4P-Ew water model, to generate conformational ensembles for disordered trialanine (Ala3), triglycine (Gly3), and trivaline (Val3) peptides. We find that the TIP4P-Ew water model yields significantly better agreement with experimentally measured scalar couplings-and therefore more accurate conformational ensembles-for both Ala3 and Gly3. For Val3, however, we find that the TIP3P and TIP4P-Ew ensembles are equivalent in their performance. To further improve the protein-water force field combination and obtain more accurate intrinsic conformational preferences, we derive a straightforward perturbation to the ϕ' backbone dihedral potential that shifts the ß-PPII equilibrium. We find that the revised ϕ' backbone dihedral potential yields improved conformational ensembles for a variety of small peptides and maintains the stability of the globular ubiquitin protein in TIP4P-Ew water.

Perturbing the folding energy landscape of the bacterial immunity protein Im7 by site-specific N-linked glycosylation.

N-linked glycosylation modulates protein folding and stability through a variety of mechanisms. As such there is considerable interest in the development of general rules to predict the structural consequences of site-specific glycosylation and to understand how these effects can be exploited in the design and development of modified proteins with advantageous properties. In this study, expressed protein ligation is used to create site-specifically glycosylated variants of the bacterial immunity protein Im7 modified with the chitobiose disaccharide (GlcNAc-GlcNAc). Glycans were introduced at seven solvent exposed sites within the Im7 sequence and the kinetic and thermodynamic consequences of N-linked glycosylation analyzed. The ΔΔG° values for glycan incorporation were found to range from +5.2 to -3.8 kJ·mol(-1). In several cases, glycosylation influences folding by modulating the local conformational preferences of the glycosylated sequence. These locally mediated effects are most prominent in the center of α-helices where glycosylation negatively effects folding and in compact turn motifs between segments of ordered secondary structure where glycosylation promotes folding and enhances the overall stability of the native protein. The studies also provide insight into why glycosylation is commonly identified at the transition between different types of secondary structure and when glycosylation may be used to elaborate protein structure to protect disordered sequences from proteolysis or immune system recognition.

The contribution of interchain salt bridges to triple-helical stability in collagen.

Studies on collagen and collagen-like peptides suggest that triple-helical stability can vary along the amino acid chain. In this regard, it has been shown that lysine residues in the Y position and acidic residues in the X' position of (GPO)(3)GXYGX'Y'(GPO)(3) peptides lead to triple-helical structures with melting temperatures similar to (GPO)(8) (where O is hydroxyproline), which is generally regarded as the most stable collagen-like sequence of this length. This enhanced stability has been attributed to the formation of salt bridges between adjacent collagen chains. In this study, we explore the relationship between interchain salt bridge formation and triple-helical stability using detailed molecular simulations. Although our results confirm that salt bridges promote triple-helical stability, we find that not all salt bridges are created equal. In particular, lysine-glutamate salt bridges are most stabilizing when formed between residues in the middle strand (B) and the trailing strand (C), whereas lysine-aspartate salt bridges are most stabilizing when formed between residues in the leading (A) and middle (B) strand-the latter observation being consistent with recent NMR data on a heterotrimeric model peptide. Overall, we believe these data clarify the role of salt bridges in modulating triple-helical stability and can be used to guide the design of collagen-like peptides that have specific interchain interactions.

Cleavage site specificity and conformational selection in type I collagen degradation.

Excessive degradation of type I collagen is associated with a variety of human diseases such as arthritis, tumor metastasis, and atherosclerosis. Methods that further our understanding of collagenolysis may therefore provide insights into the mechanism of several important disorders. Prior experiments suggest that cleavage of collagen in vitro requires intact full-length collagenase, a multidomain protein containing both a catalytic and a hemopexin-like domain. In this work we demonstrate that type I collagen can be degraded at room temperature, a temperature well below the melting temperature of type I collagen, by collagenase deletion mutants that only contain the catalytic domain of the enzyme. Furthermore, these mutant enzymes hydrolyze the same peptide bond that is recognized by the corresponding full-length enzymes. Hence enzyme specificity at room temperature is achieved without the hemopexin-like domain. We demonstrate that these findings can be explained in light of a conformational selection mechanism that dictates that collagenases preferentially recognize and cleave preformed partially unfolded states of collagen.

Physical basis of metal-binding specificity in Escherichia coli NikR.

In Escherichia coli and other bacteria, nickel uptake is regulated by the transcription factor NikR. Nickel binding at high-affinity sites in E. coli NikR (EcNikR) facilitates EcNikR binding to the nik operon, where it then suppresses transcription of genes encoding the nickel uptake transporter, NikABCDE. A structure of the EcNikR-DNA complex suggests that a second metal-binding site is also present when NikR binds to the nik operon. Moreover, this co-crystal structure raises the question of what metal occupies the second site under physiological conditions: K(+), which is present in the crystal structure, or Ni(2+), which has been proposed to bind to low- as well as high-affinity sites on EcNikR. To determine which ion is preferred at the second metal-binding site and the physical basis for any preference of one ion over another in both the second metal-binding site and the high-affinity sites, we conducted a series of detailed molecular simulations on the EcNikR structure. Simulations that place Ni(2+) at high-affinity sites lead to stable trajectories with realistic ion-ligand distances and geometries, while simulations that place K(+) at these sites lead to conformational changes in the protein that are likely unfavorable for ion binding. By contrast, simulations on the second metal site in the EcNikR-DNA complex lead to stable trajectories with realistic geometries regardless of whether K(+) or Ni(2+) occupies this site. Electrostatic binding free energy calculations, however, suggest that EcNikR binding to DNA is more favorable when the second metal-binding site contains K(+). An analysis of the energetic contributions to the electrostatic binding free energy suggests that, while the interaction between EcNikR and DNA is more favorable when the second site contains Ni(2+), the large desolvation penalty associated with moving Ni(2+) from solution to the relatively buried second site offsets this favorable interaction term. Additional free energy simulations that account for both electrostatic and non-electrostatic effects argue that EcNikR binding to DNA is most favorable when the second site contains a monovalent ion the size of K(+). Taken together, these data suggest that the EcNikR structure is most stable when Ni(2+) occupies high-affinity sites and that EcNikR binding to DNA is more favorable when the second site contains K(+).

Differential unfolding of alpha1 and alpha2 chains in type I collagen and collagenolysis.

Collagenolysis plays a central role in many disease processes and a detailed understanding of the mechanism of collagen degradation is of immense interest. While a considerable body of information about collagenolysis exists, the details of the underlying molecular mechanism are unclear. Therefore, to further our understanding of the precise mechanism of collagen degradation, we used molecular dynamics simulations to explore the structure of human type I collagen in the vicinity of the collagenase cleavage site. Since post-translational proline hydroxylation is an important step in the synthesis of collagen chains, we used the DNA sequence for the alpha1 and alpha2 chains of human type I collagen, and the known amino acid sequences for bovine and chicken type I collagen, to infer which prolines are hydroxylated in the vicinity of the collagenase cleavage site. Simulations of type I collagen in this region suggest that partial unfolding of the alpha2 chain is energetically preferred relative to unfolding of alpha1 chains. Localized unfolding of the alpha2 chain leads to the formation of a structure that has disrupted hydrogen bonds N-terminal to the collagenase cleavage site. Our data suggest that this disruption in hydrogen bonding pattern leads to increased chain flexibility, thereby enabling the alpha2 chain to sample different partially unfolded states. Surprisingly, our data also imply that alpha2 chain unfolding is mediated by the non-hydroxylation of a proline residue that is N-terminal to the cleavage site in alpha1 chains. These results suggest that hydroxylation on one chain (alpha1) can affect the structure of another chain (alpha2), and point to a critical role for the non-hydroxylation of proline residues near the collagenase cleavage site.