Probing Protein Aggregation Using the Coarse-Grained UNRES Force Field.

Protein aggregation is the cause of many, often lethal, diseases, including the Alzheimer's, Parkinson's, and Huntington's diseases, and familial amyloidosis. Theoretical investigation of the mechanism of this process, including the structures of the oligomeric intermediates which are the most toxic, is difficult because of long time scale of aggregation. Coarse-grained models, which enable us to extend the simulation time scale by three or more orders of magnitude, are, therefore, of great advantage in such studies. In this chapter, we describe the application of the physics-based UNited RESidue (UNRES) force field developed in our laboratory to study protein aggregation, in both free simulations and simulations of aggregation propagation from an existing template (seed), and illustrate it with the examples of Aß-peptide aggregation and Aß-peptide-assisted aggregation of the peptides derived from the repeat domains of tau (TauRD).

Investigation of Phosphorylation-Induced Folding of an Intrinsically Disordered Protein by Coarse-Grained Molecular Dynamics.

Apart from being the most common mechanism of regulating protein function and transmitting signals throughout the cell, phosphorylation has an ability to induce disorder-to-order transition in an intrinsically disordered protein. In particular, it was shown that folding of the intrinsically disordered protein, eIF4E-binding protein isoform 2 (4E-BP2), can be induced by multisite phosphorylation. Here, the principles that govern the folding of phosphorylated 4E-BP2 (pT37pT46 4E-BP218-62) are investigated by analyzing canonical and replica exchange molecular dynamics trajectories, generated with the coarse-grained united-residue force field, in terms of local and global motions and the time dependence of formation of contacts between Cαs of selected pairs of residues. The key residues involved in the folding of the pT37pT46 4E-BP218-62 are elucidated by this analysis. The correlations between local and global motions are identified. Moreover, for a better understanding of the physics of the formation of the folded state, the experimental structure of the pT37pT46 4E-BP218-62 is analyzed in terms of a kink (heteroclinic standing wave solution) of a generalized discrete nonlinear Schrödinger equation. It is shown that without molecular dynamics simulations the kinks are able to identify not only the phosphorylated sites of protein, the key players in folding, but also the reasons for the weak stability of the pT37pT46 4E-BP218-62.

The structure of protein dynamic space.

We use a bioinformatic description of amino acid dynamic properties, based on residue-specific average B factors, to construct a dynamics-based, large-scale description of a space of protein sequences. We examine the relationship between that space and an independently constructed, structure-based space comprising the same sequences. It is demonstrated that structure and dynamics are only moderately correlated. It is further shown that helical proteins fall into two classes with very different structure-dynamics relationships. We suggest that dynamics in the two helical classes are dominated by distinctly different modes--pseudo-one-dimensional, localized helical modes in one case, and pseudo-three-dimensional (3D) global modes in the other. Sheet/barrel and mixed-α/ß proteins exhibit more conventional structure-dynamics relationships. It is found that the strongest correlation between structure and dynamic properties arises when the latter are represented by the sequence average of the dynamic index, which corresponds physically to the overall mobility of the protein. None of these results are accessible to bioinformatic methods hitherto available.

Curvature and Torsion of Protein Main Chain as Local Order Parameters of Protein Unfolding.

Thermal protein unfolding resembles a global (two-state) phase transition. At the local scale, protein unfolding is, however, heterogeneous and probe dependent. Here, we consider local order parameters defined by the local curvature and torsion of the protein main chain. Because chemical shifts (CS's) measured by NMR spectroscopy are extremely sensitive to the local atomic environment, CS has served as a local probe of thermal unfolding of proteins by varying the position of the atomic isotope along the amino acid sequence. The variation of the CS of each Cα atom along the sequence as a function of the temperature defines a local heat-induced denaturation curve. We demonstrate that these local heat-induced denaturation curves mirror the local protein nativeness defined by the free energy landscape of the local curvature and torsion of the protein main chain described by the Cα-Cα virtual bonds. Comparison between molecular dynamics simulations and CS data of the gpW protein demonstrates that some local native states defined by the local curvature and torsion of the main chain, mainly located in secondary structures, are coupled to each other whereas others, mainly located in flexible protein segments, are not. Consequently, CS's of some residues are faithful reporters of global protein unfolding, with heat-induced denaturation curves similar to the average global one, whereas other residues remain silent about the protein unfolded state. For the latter, the local deformation of the protein main chain, characterized by its local curvature and torsion, is not cooperatively coupled to global unfolding.

New Insights into Folding, Misfolding, and Nonfolding Dynamics of a WW Domain.

Intermediate states in protein folding are associated with formation of amyloid fibrils, which are responsible for a number of neurodegenerative diseases. Therefore, prevention of the aggregation of folding intermediates is one of the most important problems to overcome. Recently, we studied the origins and prevention of formation of intermediate states with the example of the Formin binding protein 28 (FBP28) WW domain. We demonstrated that the replacement of Leu26 by Asp26 or Trp26 (in ∼15% of the folding trajectories) can alter the folding scenario from three-state folding, a major folding scenario for the FBP28 WW domain (WT) and its mutants, toward two-state or downhill folding at temperatures below the melting point. Here, for a better understanding of the physics of the formation/elimination of intermediates, (i) the dynamics and energetics of formation of ß-strands in folding, misfolding, and nonfolding trajectories of these mutants (L26D and L26W) is investigated; (ii) the experimental structures of WT, L26D, and L26W are analyzed in terms of a kink (heteroclinic standing wave solution) of a generalized discrete nonlinear Schrödinger equation. We show that the formation of each ß-strand in folding trajectories is accompanied by the emergence of kinks in internal coordinate space as well as a decrease in local free energy. In particular, the decrease in downhill folding trajectory is ∼7 kcal/mol, while it varies between 31 and 48 kcal/mol for the three-state folding trajectory. The kink analyses of the experimental structures give new insights into formation of intermediates, which may become a useful tool for preventing aggregation.

PMFF: Development of a Physics-Based Molecular Force Field for Protein Simulation and Ligand Docking.

The physics-based molecular force field (PMFF) was developed by integrating a set of potential energy functions in which each term in an intermolecular potential energy function is derived based on experimental values, such as the dipole moments, lattice energy, proton transfer energy, and X-ray crystal structures. The term "physics-based" is used to emphasize the idea that the experimental observables that are considered to be the most relevant to each term are used for the parameterization rather than parameterizing all observables together against the target value. PMFF uses MM3 intramolecular potential energy terms to describe intramolecular interactions and includes an implicit solvation model specifically developed for the PMFF. We evaluated the PMFF in three ways. We concluded that the PMFF provides reliable information based on the structure in a biological system and interprets the biological phenomena accurately by providing more accurate evidence of the biological phenomena.

Assessing the One-Bond Cα-H Spin-Spin Coupling Constants in Proteins: Pros and Cons of Different Approaches.

In the present work, we explore three different approaches for the computation of the one-bond spin-spin coupling constants (SSCC) 1JCαH in proteins: density functional theory (DFT) calculations, a Karplus-like equation, and Gaussian process regression. The main motivation of this work is to select the best method for fast and accurate computation of the 1JCαH SSCC, for its use in everyday applications in protein structure validation, refinement, and/or determination. Our initial results showed a poor agreement between the DFT-computed and observed 1JCαH SSCC values. Further analysis leads us to the understanding that the model chosen for the DFT computations is inappropriate and that more complex models will require a higher, if not prohibitively, computational cost. Finally, we show that the Karplus-like equation and Gaussian Process regression provide faster and more accurate results than DFT-based calculations.

Outline of an experimental design aimed to detect a protein A mirror image in solution.

There is abundant theoretical evidence indicating that a mirror image of Protein A may occur during the protein folding process. However, as to whether such mirror image exists in solution is an unsolved issue. Here we provide outline of an experimental design aimed to detect the mirror image of Protein A in solution. The proposal is based on computational simulations indicating that the use of a mutant of protein A, namely Q10H, could be used to detect the mirror image conformation in solution. Our results indicate that the native conformation of the protein A should have a pKa, for the Q10H mutant, at ≈6.2, while the mirror-image conformation should have a pKa close to ≈7.3. Naturally, if all the population is in the native state for the Q10H mutant, the pKa should be ≈6.2, while, if all are in the mirror-image state, it would be ≈7.3, and, if it is a mixture, the pKa should be largerthan 6.2, presumably in proportion to the mirror population. In addition, evidence is provided indicating the tautomeric distribution of H10 must also change between the native and mirror conformations. Although this may not be completely relevant for the purpose of determining whether the protein A mirror image exists in solution, it could provide valuable information to validate the pKa findings. We hope this proposal will foster experimental work on this problem either by direct application of our proposed experimental design or serving as inspiration and motivation for other experiments.

A new protein nucleic-acid coarse-grained force field based on the UNRES and NARES-2P force fields.

Based on the coarse-grained UNRES and NARES-2P models of proteins and nucleic acids, respectively, developed in our laboratory, in this work we have developed a coarse-grained model of systems containing proteins and nucleic acids. The UNRES and NARES-2P effective energy functions have been applied to the protein and nucleic-acid components of a system, respectively, while protein-nucleic-acid interactions have been described by the respective coarse-grained potentials developed in our recent work (Yin et al., J. Chem Theory Comput. 2015, 11, 1792). The Debye-Hückel screening has been applied to the electrostatic-interaction energy between the phosphate groups and charged amino-acid side chains. The model has been integrated into the UNRES package for coarse-grained molecular dynamics simulations of proteins and the implementation has been tested for energy conservation in microcanonical molecular dynamics runs and for temperature conservation in canonical molecular dynamics runs. Two case studies were performed: (i) the dynamics of the Ku protein heterodimer bound to DNA, for which it was found that the Ku70/Ku80 protein complex plays an active role in DNA repairing and (ii) conformational changes of the multiple antibiotic resistance (MarA) protein occurring during DNA binding, for which the functionally important motions occurring during this process were identified. © 2018 Wiley Periodicals, Inc.

Dependence of the Formation of Tau and Aß Peptide Mixed Aggregates on the Secondary Structure of the N-Terminal Region of Aß.

One of the hallmarks of Alzheimer's disease is the formation of aggregates of the tau protein, a process that can be facilitated by the presence of fibrils formed by the amyloid ß peptide (Aß). However, the mechanism that triggers tau aggregation is still a matter of debate. The effect of Aß40 fibrils on the aggregation of the repeat domain of tau (TauRD) is investigated here by employing coarse-grained molecular dynamics simulations. The results indicate that the repeat domain of tau has a high affinity for Aß40 fibrils, with the 261GSTENLK267 fragment of tau driving TauRD toward the 16KLVFFA21 fragment in Aß40. Monomeric Aß40, in which the 16KLVFFA21 fragment is rarely found in an extended conformation (as in the fibril), has a low affinity for the TauRD, indicating that the ability of Aß40 fibrils to bind to the TauRD depends on the 16KLVFFA21 fragment of Aß adopting an extended conformation.

Statistical Model To Decipher Protein Folding/Unfolding at a Local Scale.

Protein folding/unfolding can be analyzed experimentally at a local scale by monitoring the physical properties of local probes as a function of the temperature, for example, the distance between fluorophores or the values of chemical shifts of backbone atoms. Here, the analytical Lifson-Roig model for the helix-coil transition is modified to analyze local thermal unfolding of the fast-folder W protein of bacteriophage lambda (gpW) simulated by all-atom molecular dynamics (MD) simulations in explicit solvent at 15 different temperatures. The protein structure is described by the coarse-grained dihedral angles (γ) and bond angles (θ) built between successive Cα-Cα virtual bonds. Each (γ,θ) pair serves as a local probe of protein unfolding. Local native/non-native states are defined for each pair of (γ,θ) angles by analyzing the free-energy landscapes Δ G(γ,θ) computed from MD trajectories. The three local elementary equilibrium constants of the model are extracted for each (γ,θ) pair along the sequence from MD simulations, and the model predictions are compared to the MD data. Using only the local equilibrium constants as an input, we show that the local denaturation curves computed from the model partition function fit their MD simulated counterparts in a satisfying manner without any adjustment. In the model and MD simulations, gpW unfolds gradually between 320 and 340 K, with an average native percentage decreasing from 0.8 (320 K) to 0.2 (340 K). In the prism of the model, there is no stable structure at the local scale in this 20 K unfolding temperature range. The enthalpy change upon local unfolding computed from the model and from MD trajectories suggests that the unfolded state between 320 and 340 K corresponds to a dynamical equilibrium between a large ensemble of constantly changing structures. The present results confirm the downhill unfolding of gpW, which does not obey a two-state global folding/unfolding model, and shed light on the interpretation of local denaturation curves.

From a Highly Disordered to a Metastable State: Uncovering Insights of α-Synuclein.

α-Synuclein (αS) is a major constituent of Lewy bodies, the insoluble aggregates that are the hallmark of one of the most prevalent neurodegenerative disorders, Parkinson's disease (PD). The vast majority of experiments in vitro and in vivo provide extensive evidence that a disordered monomeric form is the predominant state of αS in water solution, and it undergoes a large-scale disorder-to-helix transition upon binding to vesicles of different types. Recently, another form, tetrameric, of αS with a stable helical structure was identified experimentally. It has been shown that a dynamic intracellular population of metastable αS tetramers and monomers coexists normally; and the tetramer plays an essential role in maintaining αS homeostasis. Therefore, it is of interest to know whether the tetramer can serve as a means of preventing or delaying the start of PD. Before answering this very important question, it is, first, necessary to find out, on an atomistic level, a correlation between tetramers and monomers; what mediates tetramer formation and what makes a tetramer stable. We address these questions here by investigating both monomeric and tetrameric forms of αS. In particular, by examining correlations between the motions of the side chains and the main chain, steric parameters along the amino-acid sequence, and one- and two-dimensional free-energy landscapes along the coarse-grained dihedral angles γ and δ and principal components, respectively, in monomeric and tetrameric αS, we were able to shed light on a fundamental relationship between monomers and tetramers, and the key residues involved in mediating formation of a tetramer. Also, the reasons for the stability of tetrameric αS and inability of monomeric αS to fold are elucidated here.

Lysosomal enzyme tripeptidyl peptidase 1 destabilizes fibrillar Aß by multiple endoproteolytic cleavages within the ß-sheet domain.

Accumulation of amyloid-beta (Aß), which is associated with Alzheimer's disease, can be caused by excess production or insufficient clearance. Because of its ß-sheet structure, fibrillar Aß is resistant to proteolysis, which would contribute to slow degradation of Aß plaques in vivo. Fibrillar Aß can be internalized by microglia, which are the scavenger cells of the brain, but the fibrils are degraded only slowly in microglial lysosomes. Cathepsin B is a lysosomal protease that has been shown to proteolyze fibrillar Aß. Tripeptidyl peptidase 1 (TPP1), a lysosomal serine protease, possesses endopeptidase activity and has been shown to cleave peptides between hydrophobic residues. Herein, we demonstrate that TPP1 is able to proteolyze fibrillar Aß efficiently. Mass spectrometry analysis of peptides released from fibrillar Aß digested with TPP1 reveals several endoproteolytic cleavages including some within ß-sheet regions that are important for fibril formation. Using molecular dynamics simulations, we demonstrate that these cleavages destabilize fibrillar ß-sheet structure. The demonstration that TPP1 can degrade fibrillar forms of Aß provides insight into the turnover of fibrillar Aß and may lead to new therapeutic methods to increase degradation of Aß plaques.

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.

A comprehensive analysis of the computed tautomer fractions of the imidazole ring of histidines in Loligo vulgaris.

A recently introduced electrostatic-based method to determine the pKa values of ionizable residues and fractions of ionized and tautomeric forms of histidine (His) and acid residues in proteins, at a given fixed pH, is applied here to the analysis of a His-rich protein, namely Loligo vulgaris (pdb id 1E1A), a 314-residue all-ß protein. The average tautomeric fractions for the imidazole ring of each of the six histidines in the sequence were computed using an approach that includes, but is not limited to, molecular dynamic simulations coupled with calculations of the ionization states for all 94 ionizable residues of protein 1E1A in water at pH 6.5 and 300 K. The electrostatic-calculated tautomeric fractions of the imidazole ring of His were compared with predictions obtained from an existent NMR-based methodology. Our results indicate that: (i) the averaged electrostatic-based tautomeric predictions for the imidazole ring of all histidines of Loligo vulgaris are dominated by the Nε2-H rather than the Nδ1-H form, although such preferences from the NMR-based methodology are not so well defined; (ii) the computed average absolute difference between the electrostatic- and the NMR-based tautomeric predictions among all six histidines vary among 0% to 17%; (iii) for the His showing the largest fraction of the neutral form (81%), the absolute difference between the NMR- and electrostatic-based computed tautomeric predictions is only 3%; and (iv) the tautomeric predictions for the imidazole ring of His computed with the NMR-based methodology are stable within a certain, well-defined, range of variations of a tautomer-related parameter.

Dynamics of Disulfide-Bond Disruption and Formation in the Thermal Unfolding of Ribonuclease A.

Ribonuclease A (RNase A) is the most studied member of ribonucleases - a group of enzymes responsible for catalyzing RNA degradation. RNase A contains four disulfide bonds, which were found to be necessary for the native structure of the protein to form. In this work the kinetics and thermodynamics of RNase unfolding were studied by means of a series of coarse-grained canonical molecular-dynamics simulations, run at various temperatures, and replica-exchange molecular dynamics simulations with the UNRES force field, which is capable of dynamic formation and breaking the disulfide bonds during the course of the simulations. It was found that the Cys40-Cys95 bond was the first to break, while the Cys26-Cys84 bond was the last to break, independent of temperature, in agreement with available experimental data. Except for the Cys40-Cys95 bond, all disulfide bonds were relatively stable during the simulations. The formation/disruption of disulfide bonds was found to be temperature dependent for three out of four disulfide bonds in RNase A, except for the most stable disulfide bond between Cys65 and Cys72. A stable intermediate without the Cys40-Cys95 disulfide bond, with structure similar to that of the most common folding intermediate observed experimentally, was found in simulated unfolding. In agreement with experiment, non-native disulfide bonds were also observed. By analyzing residue-position fluctuations, it was found that native disulfide bonds are located in the highly flexible regions of the protein, which is probably why their presence is necessary for the stability of RNase A.

Limiting Values of the one-bond C-H Spin-Spin Coupling Constants of the Imidazole Ring of Histidine at High-pH.

Assessment of the relative amounts of the forms of the imidazole ring of Histidine (His), namely the protonated (H+) and the tautomeric Nε2-H and Nδ1-H forms, respectively, is a challenging task in NMR spectroscopy. Indeed, their determination by direct observation of the 15N and 13C chemical shifts or the one-bond C-H, 1JCH, Spin-Spin Coupling Constants (SSCC) requires knowledge of the "canonical" limiting values of these forms in which each one is present to the extent of 100%. In particular, at high-pH, an accurate determination of these "canonical" limiting values, at which the tautomeric forms of His coexist, is an elusive problem in NMR spectroscopy. Among different NMR-based approaches to treat this problem, we focus here on the computation, at the DFT level of theory, of the high-pH limiting value for the 1JCH SSCC of the imidazole ring of His. Solvent effects were considered by using the polarizable continuum model approach. The results of this computation suggest, first, that the value of 1JCε1H = 205 ± 1.0 Hz should be adopted as the canonical high-pH limiting value for this SSCC; second, the variation of 1JCε1H SSCC during tautomeric changes is minor, i.e., within ±1Hz; and, finally, the value of 1JCδ2H SSCC upon tautomeric changes is large (15 Hz) indicating that, at high-pH or for non-protonated His at any pH, the tautomeric fractions of the imidazole ring of His can be predicted accurately as a function of the observed value of 1JCδ2H SSCC.

Maximum Likelihood Calibration of the UNRES Force Field for Simulation of Protein Structure and Dynamics.

By using the maximum likelihood method for force-field calibration recently developed in our laboratory, which is aimed at achieving the agreement between the simulated conformational ensembles of selected training proteins and the corresponding ensembles determined experimentally at various temperatures, the physics-based coarse-grained UNRES force field for simulations of protein structure and dynamics was optimized with seven small training proteins exhibiting a variety of secondary and tertiary structures. Four runs of optimization, in which the number of optimized force-field parameters was gradually increased, were carried out, and the resulting force fields were subsequently tested with a set of 22 α-, 12 ß-, and 12 α + ß-proteins not used in optimization. The variant in which energy-term weights, local, and correlation potentials, side-chain radii, and anisotropies were optimized turned out to be the most transferable and outperformed all previous versions of UNRES on the test set.

Sequence-, structure-, and dynamics-based comparisons of structurally homologous CheY-like proteins.

We recently introduced a physically based approach to sequence comparison, the property factor method (PFM). In the present work, we apply the PFM approach to the study of a challenging set of sequences-the bacterial chemotaxis protein CheY, the N-terminal receiver domain of the nitrogen regulation protein NT-NtrC, and the sporulation response regulator Spo0F. These are all response regulators involved in signal transduction. Despite functional similarity and structural homology, they exhibit low sequence identity. PFM sequence comparison demonstrates a statistically significant qualitative difference between the sequence of CheY and those of the other two proteins that is not found using conventional alignment methods. This difference is shown to be consonant with structural characteristics, using distance matrix comparisons. We also demonstrate that residues participating strongly in native contacts during unfolding are distributed differently in CheY than in the other two proteins. The PFM result is also in accord with dynamic simulation results of several types. Molecular dynamics simulations of all three proteins were carried out at several temperatures, and it is shown that the dynamics of CheY are predicted to differ from those of NT-NtrC and Spo0F. The predicted dynamic properties of the three proteins are in good agreement with experimentally determined B factors and with fluctuations predicted by the Gaussian network model. We pinpoint the differences between the PFM and traditional sequence comparisons and discuss the informatic basis for the ability of the PFM approach to detect physical differences between these sequences that are not apparent from traditional alignment-based comparison.