Parametrization of Backbone−Electrostatic and Multibody Contributions to the UNRES Force Field for Protein-Structure Prediction from Ab Initio Energy Surfaces of Model Systems

Optimization of parallel implementation of UNRES package for coarse-grained simulations to treat large proteins

We report major algorithmic improvements of the UNRES package for physics-based coarse-grained simulations of proteins. These include (i) introduction of interaction lists to optimize computations, (ii) transforming the inertia matrix to a pentadiagonal form to reduce computing and memory requirements, (iii) removing explicit angles and dihedral angles from energy expressions and recoding the most time-consuming energy/force terms to minimize the number of operations and to improve numerical stability, (iv) using OpenMP to parallelize those sections of the code for which distributed-memory parallelization involves unfavorable computing/communication time ratio, and (v) careful memory management to minimize simultaneous access of distant memory sections. The new code enables us to run molecular dynamics simulations of protein systems with size exceeding 100,000 amino-acid residues, reaching over 1 ns/day (1 μs/day in all-atom timescale) with 24 cores for proteins of this size. Parallel performance of the code and comparison of its performance with that of AMBER, GROMACS and MARTINI 3 is presented.

Modeling the Structure, Dynamics, and Transformations of Proteins with the UNRES Force Field

The physics-based united-residue (UNRES) model of proteins ( www.unres.pl ) has been designed to carry out large-scale simulations of protein folding. The force field has been derived and parameterized based on the principles of statistical-mechanics, which makes it independent of structural databases and applicable to treat nonstandard situations such as, proteins that contain D-amino-acid residues. Powered by Langevin dynamics and its replica-exchange extensions, UNRES has found a variety of applications, including ab initio and database-assisted protein-structure prediction, simulating protein-folding pathways, exploring protein free-energy landscapes, and solving biological problems. This chapter provides a summary of UNRES and a guide for potential users regarding the application of the UNRES package in a variety of research tasks.

Theory and Practice of Coarse-Grained Molecular Dynamics of Biologically Important Systems

Molecular dynamics with coarse-grained models is nowadays extensively used to simulate biomolecular systems at large time and size scales, compared to those accessible to all-atom molecular dynamics. In this review article, we describe the physical basis of coarse-grained molecular dynamics, the coarse-grained force fields, the equations of motion and the respective numerical integration algorithms, and selected practical applications of coarse-grained molecular dynamics. We demonstrate that the motion of coarse-grained sites is governed by the potential of mean force and the friction and stochastic forces, resulting from integrating out the secondary degrees of freedom. Consequently, Langevin dynamics is a natural means of describing the motion of a system at the coarse-grained level and the potential of mean force is the physical basis of the coarse-grained force fields. Moreover, the choice of coarse-grained variables and the fact that coarse-grained sites often do not have spherical symmetry implies a non-diagonal inertia tensor. We describe selected coarse-grained models used in molecular dynamics simulations, including the most popular MARTINI model developed by Marrink's group and the UNICORN model of biological macromolecules developed in our laboratory. We conclude by discussing examples of the application of coarse-grained molecular dynamics to study biologically important processes.

ESCASA: Analytical estimation of atomic coordinates from coarse-grained geometry for nuclear-magnetic-resonance-assisted protein structure modeling. I. Backbone and Hβ protons

A method for the estimation of coordinates of atoms in proteins from coarse-grained geometry by simple analytical formulas (ESCASA), for use in nuclear-magnetic-resonance (NMR) data-assisted coarse-grained simulations of proteins is proposed. In this paper, the formulas for the backbone H^α and amide (H^N ) protons, and the side-chain H^β protons, given the C^α -trace, have been derived and parameterized, by using the interproton distances calculated from a set of 140 high-resolution non-homologous protein structures. The mean standard deviation over all types of proton pairs in the set was 0.44 Å after fitting. Validation against a set of 41 proteins with NMR-determined structures, which were not considered in parameterization, resulted in average standard deviation from average proton-proton distances of the NMR-determined structures of 0.25 Å, compared to 0.21 Å obtained with the PULCHRA all-atom-chain reconstruction algorithm and to the 0.12 Å standard deviation of the average-structure proton-proton distance of NMR-determined ensembles. The formulas provide analytical forces and can, therefore, be used in coarse-grained molecular dynamics.

Scale-consistent approach to the derivation of coarse-grained force fields for simulating structure, dynamics, and thermodynamics of biopolymers

In this chapter the scale-consistent approach to the derivation of coarse-grained force fields developed in our laboratory is presented, in which the effective energy function originates from the potential of mean force of the system under consideration and embeds atomistically detailed interactions in the resulting energy terms through use of Kubo's cluster-cumulant expansion, appropriate selection of the major degrees of freedom to be averaged out in the derivation of analytical approximations to the energy terms, and appropriate expression of the interaction energies at the all-atom level in these degrees of freedom. Our approach enables the developers to find correct functional forms of the effective coarse-grained energy terms, without having to import them from all-atom force fields or deriving them on a heuristic basis. In particular, the energy terms derived in such a way exhibit correct dependence on coarse-grained geometry, in particular on site orientation. Moreover, analytical formulas for the multibody (correlation) terms, which appear to be crucial for coarse-grained modeling of many of the regular structures such as, e.g., protein α-helices and β-sheets, can be derived in a systematic way. Implementation of the developed theory to the UNIfied COarse-gRaiNed (UNICORN) model of biological macromolecules, which consists of the UNRES (for proteins), NARES-2P (for nucleic acids), and SUGRES-1P (for polysaccharides) components, and is being developed in our laboratory is described. Successful applications of UNICORN to the prediction of protein structure, simulating the folding and stability of proteins and nucleic acids, and solving biological problems are discussed.

Extension of the UNRES Coarse-Grained Force Field to Membrane Proteins in the Lipid Bilayer

The physics-based UNRES coarse-grained force field for the simulations of protein structure and dynamics has been extended to treat membrane proteins. The lipid bilayer has been modeled by introducing a continuous nonpolar phase with the water-interface region of appropriate thickness. The potentials for average electrostatic and correlation interactions of the peptide groups have been rescaled to account for the reduction of the dielectric permittivity compared to the water phase and new potentials for protein side-chain-side-chain interactions inside and across the lipid phase have been introduced. The model was implemented in the UNRES package for coarse-grained simulations of proteins, and the package with the new functionality was tested for total energy conservation and thermostat behavior in microcanonical and canonical molecular dynamics simulations runs, respectively. The method was validated by running unrestricted ab initio blind-prediction tests of 10 short α-helical membrane proteins, all runs started from the extended structures. The modified UNRES force field was able to predict correctly the overall folds of the membrane proteins studied.

A general method for the derivation of the functional forms of the effective energy terms in coarse-grained energy functions of polymers. III. Determination of scale-consistent backbone-local and correlation potentials in the UNRES force field and force-field calibration and validation

The general theory of the construction of scale-consistent energy terms in the coarse-grained force fields presented in Paper I of this series has been applied to the revision of the UNRES force field for physics-based simulations of proteins. The potentials of mean force corresponding to backbone-local and backbone-correlation energy terms were calculated from the ab initio energy surfaces of terminally blocked glycine, alanine, and proline, and the respective analytical expressions, derived by using the scale-consistent formalism, were fitted to them. The parameters of all these potentials depend on single-residue types, thus reducing their number and preventing over-fitting. The UNRES force field with the revised backbone-local and backbone-correlation terms was calibrated with a set of four small proteins with basic folds: tryptophan cage variant (TRP1; α), Full Sequence Design (FSD; α + β), villin headpiece (villin; α), and a truncated FBP-28 WW-domain variant (2MWD; β) (the NEWCT-4P force field) and, subsequently, with an enhanced set of 9 proteins composed of TRP1, FSD, villin, 1BDC (α), 2I18 (α), 1QHK (α + β), 2N9L (α + β), 1E0L (β), and 2LX7 (β) (the NEWCT-9P force field). The NEWCT-9P force field performed better than NEWCT-4P in a blind-prediction-like test with a set of 26 proteins not used in calibration and outperformed, in a test with 76 proteins, the most advanced OPT-WTFSA-2 version of UNRES with former backbone-local and backbone-correlation terms that contained more energy terms and more optimizable parameters. The NEWCT-9P force field reproduced the bimodal distribution of backbone-virtual-bond angles in the simulated structures, as observed in experimental protein structures.

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.

A general method for the derivation of the functional forms of the effective energy terms in coarse-grained energy functions of polymers. I. Backbone potentials of coarse-grained polypeptide chains

A general and systematic method for the derivation of the functional expressions for the effective energy terms in coarse-grained force fields of polymer chains is proposed. The method is based on the expansion of the potential of mean force of the system studied in the cluster-cumulant series and expanding the all-atom energy in the Taylor series in the squares of interatomic distances about the squares of the distances between coarse-grained centers, to obtain approximate analytical expressions for the cluster cumulants. The primary degrees of freedom to average about are the angles for collective rotation of the atoms contained in the coarse-grained interaction sites about the respective virtual-bond axes. The approach has been applied to the revision of the virtual-bond-angle, virtual-bond-torsional, and backbone-local-and-electrostatic correlation potentials for the UNited RESidue (UNRES) model of polypeptide chains, demonstrating the strong dependence of the torsional and correlation potentials on virtual-bond angles, not considered in the current UNRES. The theoretical considerations are illustrated with the potentials calculated from the ab initiopotential-energysurface of terminally blocked alanine by numerical integration and with the statistical potentials derived from known protein structures. The revised torsional potentials correctly indicate that virtual-bond angles close to 90° result in the preference for the turn and helical structures, while large virtual-bond angles result in the preference for polyproline II and extended backbone geometry. The revised correlation potentials correctly reproduce the preference for the formation of β-sheet structures for large values of virtual-bond angles and for the formation of α-helical structures for virtual-bond angles close to 90°.

The Role of Conformational Entropy in the Determination of Structural-Kinetic Relationships for Helix-Coil Transitions

Coarse-grained molecular simulation models can provide significant insight into the complex behavior of protein systems, but suffer from an inherently distorted description of dynamical properties. We recently demonstrated that, for a heptapeptide of alanine residues, the structural and kinetic properties of a simulation model are linked in a rather simple way, given a certain level of physics present in the model. In this work, we extend these findings to a longer peptide, for which the representation of configuration space in terms of a full enumeration of sequences of helical/coil states along the peptide backbone is impractical. We verify the structural-kinetic relationships by scanning the parameter space of a simple native-biased model and then employ a distinct transferable model to validate and generalize the conclusions. Our results further demonstrate the validity of the previous findings, while clarifying the role of conformational entropy in the determination of the structural-kinetic relationships. More specifically, while the global, long timescale kinetic properties of a particular class of models with varying energetic parameters but approximately fixed conformational entropy are determined by the overarching structural features of the ensemble, a shift in these kinetic observables occurs for models with a distinct representation of steric interactions. At the same time, the relationship between structure and more local, faster kinetic properties is not affected by varying the conformational entropy of the model.

Shielding effect in protein folding

One of the most important interactions responsible for protein folding and stability are hydrogen bonds between peptide groups. There is a constant competition between the water molecules and peptide groups in a hydrogen bond formation. Also side-chains take part in this process by reducing hydration of peptide group (shielding effect) that promotes the protein folding. In this paper, a new approach to take into account a shielding effect is presented. A modification of the energy function is derived and incorporated into the UNited RESidue (UNRES) force field. Canonical Molecular Dynamics and Replica Exchange Molecular Dynamics with UNRES force field is applied to study the influence of this effect on protein structure, folding kinetics and free energy landscapes. The results of test calculations suggest that even small contribution of this effect into energy function changes force field behavior as well as speeds up the folding process significantly.

Role of the sulfur to α-carbon thioether bridges in thurincin H

Thurincin H is a small protein produced by Bacillus thuringiensis SF361 with gram-positive antimicrobial properties. The toxins produced by B. thuringiensis are widely used in the agriculture as, e.g. natural preservatives in dairy products. The structure of thurincin H possesses four covalent sulfur to [Formula: see text]-carbon bonds that involve the cysteine side-chains; these bonds are probably responsible for the shape and stability of the protein and, thereby, for its antimicrobial properties. To examine the influence of the formation of the sulfur-carbon bonds on the folding pathways and stability of the protein, a series of canonical and multiplexed replica-exchange simulations with the coarse-grained UNRES force field was carried out without and with distance restraints imposed on selected S-C[Formula: see text] atom pairs. It was found that the order of the formation and breaking of the S-C[Formula: see text] thioether bonds significantly impacts on the foldability and stability of the thurincin H. It was also observed that thioether bridges play a major role in stabilizing the global fold of the protein, although it significantly diminishes the entropy of the system. The maximum foldability of thurincin H was observed in the presence of the optimal set of three out of four thioether bridges. Thus, the results suggest that the presence of ThnB enzyme and other agents that catalyze the formation of thioether bridges can be essential for correct folding of thurincin H and that the formation of the fourth bridge does not seem to facilitate folding; instead, it seems to rigidify the loop and prevent proteolysis.

Optimization of a Nucleic Acids united-RESidue 2-Point model (NARES-2P) with a maximum-likelihood approach

Coarse-grained models are useful tools to investigate the structural and thermodynamic properties of biomolecules. They are obtained by merging several atoms into one interaction site. Such simplified models try to capture as much as possible information of the original biomolecular system in all-atom representation but the resulting parameters of these coarse-grained force fields still need further optimization. In this paper, a force field optimization method, which is based on maximum-likelihood fitting of the simulated to the experimental conformational ensembles and least-squares fitting of the simulated to the experimental heat-capacity curves, is applied to optimize the Nucleic Acid united-RESidue 2-point (NARES-2P) model for coarse-grained simulations of nucleic acids recently developed in our laboratory. The optimized NARES-2P force field reproduces the structural and thermodynamic data of small DNA molecules much better than the original force field.

Performance of protein-structure predictions with the physics-based UNRES force field in CASP11

Participating as the Cornell-Gdansk group, we have used our physics-based coarse-grained UNited RESidue (UNRES) force field to predict protein structure in the 11th Community Wide Experiment on the Critical Assessment of Techniques for Protein Structure Prediction (CASP11). Our methodology involved extensive multiplexed replica exchange simulations of the target proteins with a recently improved UNRES force field to provide better reproductions of the local structures of polypeptide chains. All simulations were started from fully extended polypeptide chains, and no external information was included in the simulation process except for weak restraints on secondary structure to enable us to finish each prediction within the allowed 3-week time window. Because of simplified UNRES representation of polypeptide chains, use of enhanced sampling methods, code optimization and parallelization and sufficient computational resources, we were able to treat, for the first time, all 55 human prediction targets with sizes from 44 to 595 amino acid residues, the average size being 251 residues. Complete structures of six single-domain proteins were predicted accurately, with the highest accuracy being attained for the T0769, for which the CαRMSD was 3.8 Å for 97 residues of the experimental structure. Correct structures were also predicted for 13 domains of multi-domain proteins with accuracy comparable to that of the best template-based modeling methods. With further improvements of the UNRES force field that are now underway, our physics-based coarse-grained approach to protein-structure prediction will eventually reach global prediction capacity and, consequently, reliability in simulating protein structure and dynamics that are important in biochemical processes.

AVAILABILITY AND IMPLEMENTATION

Freely available on the web at http://www.unres.pl/ CONTACT: has5@cornell.edu.

Molecular dynamics of protein A and a WW domain with a united-residue model including hydrodynamic interaction

The folding of the N-terminal part of the B-domain of staphylococcal protein A (PDB ID: 1BDD, a 46-residue three-α-helix bundle) and the formin-binding protein 28 WW domain (PDB ID: 1E0L, a 37-residue three-stranded anti-parallel β protein) was studied by means of Langevin dynamics with the coarse-grained UNRES force field to assess the influence of hydrodynamic interactions on protein-folding pathways and kinetics. The unfolded, intermediate, and native-like structures were identified by cluster analysis, and multi-exponential functions were fitted to the time dependence of the fractions of native and intermediate structures, respectively, to determine bulk kinetics. It was found that introducing hydrodynamic interactions slows down both the formation of an intermediate state and the transition from the collapsed structures to the final native-like structures by creating multiple kinetic traps. Therefore, introducing hydrodynamic interactions considerably slows the folding, as opposed to the results obtained from earlier studies with the use of Gō-like models.

A Nonradial Coarse-Grained Potential for Proteins Produces Naturally Stable Secondary Structure Elements

We introduce a nonradial potential term for coarse-grained (CG) molecular simulations of proteins. This term mimics the backbone dipole-dipole interactions and accounts for the needed directionality to form stable folded secondary structure elements. We show that α-helical and β-sheet peptide chains are correctly described in dynamics without the need of introducing any a priori bias potentials or ad hoc parametrizations, which limit broader applicability of CG simulations for proteins. Moreover, our model is able to catch the formation of supersecondary structural motifs, like transitions from long single α-helices to helix-coil-helix or β-hairpin assemblies. This novel scheme requires the structural information of Cα beads only; it does not introduce any additional degrees of freedom to the system and has a general formulation, which allows it to be used in synergy with various CG protocols, leading to an improved description of the structural and dynamic properties of protein assemblies and networks.

Parametrization of Backbone Flexibility in a Coarse-Grained Force Field for Proteins (COFFDROP) Derived from All-Atom Explicit-Solvent Molecular Dynamics Simulations of All Possible Two-Residue Peptides

Recently, we reported the parametrization of a set of coarse-grained (CG) nonbonded potential functions, derived from all-atom explicit-solvent molecular dynamics (MD) simulations of amino acid pairs and designed for use in (implicit-solvent) Brownian dynamics (BD) simulations of proteins; this force field was named COFFDROP (COarse-grained Force Field for Dynamic Representations Of Proteins). Here, we describe the extension of COFFDROP to include bonded backbone terms derived from fitting to results of explicit-solvent MD simulations of all possible two-residue peptides containing the 20 standard amino acids, with histidine modeled in both its protonated and neutral forms. The iterative Boltzmann inversion (IBI) method was used to optimize new CG potential functions for backbone-related terms by attempting to reproduce angle, dihedral, and distance probability distributions generated by the MD simulations. In a simple test of the transferability of the extended force field, the angle, dihedral, and distance probability distributions obtained from BD simulations of 56 three-residue peptides were compared to results from corresponding explicit-solvent MD simulations. In a more challenging test of the COFFDROP force field, it was used to simulate eight intrinsically disordered proteins and was shown to quite accurately reproduce the experimental hydrodynamic radii (Rhydro), provided that the favorable nonbonded interactions of the force field were uniformly scaled downward in magnitude. Overall, the results indicate that the COFFDROP force field is likely to find use in modeling the conformational behavior of intrinsically disordered proteins and multidomain proteins connected by flexible linkers.

Structural studies of the C-terminal 19-peptide of serum amyloid A and its Pro → Ala variants interacting with human cystatin C

Serum amyloid A (SAA) is a multifunctional acute-phase protein whose concentration in serum increases markedly following a number of chronic inflammatory and neoplastic diseases. Prolonged high SAA level may give rise to reactive systemic amyloid A (AA) amyloidosis, where the N-terminal segment of SAA is deposited as amyloid fibrils. Besides, recently, well-documented association of SAA with high-density lipoprotein or glycosaminoglycans, in particular heparin/heparin sulfate (HS), and specific interaction between SAA and human cystatin C (hCC), the ubiquitous inhibitor of cysteine proteases, was proved. Using a combination of selective proteolytic excision and high-resolution mass spectrometry, a hCC binding site in the SAA sequence was determined as SAA(86-104). The role of this SAA C-terminal fragment as a ligand-binding locus is still not clear. It was postulated important in native SAA folding and in pathogenesis of AA amyloidosis. In the search of conformational details of this SAA fragment, we did its structure and affinity studies, including its selected double/triple Pro → Ala variants. Our results clearly show that the SAA(86-104) 19-peptide has rather unordered structure with bends in its C-terminal part, which is consistent with the previous results relating to the whole protein. The results of affinity chromatography, fluorescent ELISA-like test, CD and NMR studies point to an importance of proline residues on structure of SAA(86-104). Conformational details of SAA fragment, responsible for hCC binding, may help to understand the objective of hCC-SAA complex formation and its importance for pathogenesis of reactive amyloid A amyloidosis.

Folding kinetics of WW domains with the united residue force field for bridging microscopic motions and experimental measurements

To demonstrate the utility of the coarse-grained united-residue (UNRES) force field to compare experimental and computed kinetic data for folding proteins, we have performed long-time millisecond-timescale canonical Langevin molecular dynamics simulations of the triple β-strand from the Formin binding protein 28 WW domain and six nonnatural variants, using UNRES. The results have been compared with available experimental data in both a qualitative and a quantitative manner. Complexities of the folding pathways, which cannot be determined experimentally, were revealed. The folding mechanisms obtained from the simulated folding kinetics are in agreement with experimental results, with a few discrepancies for which we have accounted. The origins of single- and double-exponential kinetics and their correlations with two- and three-state folding scenarios are shown to be related to the relative barrier heights between the various states. The rate constants obtained from time profiles of the fractions of the native, intermediate, and unfolded structures, and the kinetic equations fitted to them, correlate with the experimental values; however, they are about three orders of magnitude larger than the experimental ones for most of the systems. These differences are in agreement with the timescale extension derived by scaling down the friction of water and averaging out the fast degrees of freedom when passing from all-atom to a coarse-grained representation. Our results indicate that the UNRES force field can provide accurate predictions of folding kinetics of these WW domains, often used as models for the study of the mechanisms of proein folding.

Application of information theory to a three-body coarse-grained representation of proteins in the PDB: insights into the structural and evolutionary roles of residues in protein structure

Knowledge-based methods for analyzing protein structures, such as statistical potentials, primarily consider the distances between pairs of bodies (atoms or groups of atoms). Considerations of several bodies simultaneously are generally used to characterize bonded structural elements or those in close contact with each other, but historically do not consider atoms that are not in direct contact with each other. In this report, we introduce an information-theoretic method for detecting and quantifying distance-dependent through-space multibody relationships between the sidechains of three residues. The technique introduced is capable of producing convergent and consistent results when applied to a sufficiently large database of randomly chosen, experimentally solved protein structures. The results of our study can be shown to reproduce established physico-chemical properties of residues as well as more recently discovered properties and interactions. These results offer insight into the numerous roles that residues play in protein structure, as well as relationships between residue function, protein structure, and evolution. The techniques and insights presented in this work should be useful in the future development of novel knowledge-based tools for the evaluation of protein structure.

A unified coarse-grained model of biological macromolecules based on mean-field multipole-multipole interactions

A unified coarse-grained model of three major classes of biological molecules—proteins, nucleic acids, and polysaccharides—has been developed. It is based on the observations that the repeated units of biopolymers (peptide groups, nucleic acid bases, sugar rings) are highly polar and their charge distributions can be represented crudely as point multipoles. The model is an extension of the united residue (UNRES) coarse-grained model of proteins developed previously in our laboratory. The respective force fields are defined as the potentials of mean force of biomacromolecules immersed in water, where all degrees of freedom not considered in the model have been averaged out. Reducing the representation to one center per polar interaction site leads to the representation of average site–site interactions as mean-field dipole–dipole interactions. Further expansion of the potentials of mean force of biopolymer chains into Kubo’s cluster-cumulant series leads to the appearance of mean-field dipole–dipole interactions, averaged in the context of local interactions within a biopolymer unit. These mean-field interactions account for the formation of regular structures encountered in biomacromolecules, e.g., α-helices and β-sheets in proteins, double helices in nucleic acids, and helicoidally packed structures in polysaccharides, which enables us to use a greatly reduced number of interacting sites without sacrificing the ability to reproduce the correct architecture. This reduction results in an extension of the simulation timescale by more than four orders of magnitude compared to the all-atom representation. Examples of the performance of the model are presented.

Figure

Revised Backbone-Virtual-Bond-Angle Potentials to Treat the l- and d-Amino Acid Residues in the Coarse-Grained United Residue (UNRES) Force Field

Continuing our effort to introduce d-amino-acid residues in the united residue (UNRES) force field developed in our laboratory, in this work the Cα ··· Cα ··· Cα backbone-virtual-bond-valence-angle (θ) potentials for systems containing d-amino-acid residues have been developed. The potentials were determined by integrating the combined energy surfaces of all possible triplets of terminally blocked glycine, alanine, and proline obtained with ab initio molecular quantum mechanics at the MP2/6-31G(d,p) level to calculate the corresponding potentials of mean force (PMFs). Subsequently, analytical expressions were fitted to the PMFs to give the virtual-bond-valence potentials to be used in UNRES. Alanine represented all types of amino-acid residues except glycine and proline. The blocking groups were either the N-acetyl and N',N'-dimethyl or N-acetyl and pyrrolidyl group, depending on whether the residue next in sequence was an alanine-type or a proline residue. A total of 126 potentials (63 symmetry-unrelated potentials for each set of terminally blocking groups) were determined. Together with the torsional, double-torsional, and side-chain-rotamer potentials for polypeptide chains containing d-amino-acid residues determined in our earlier work (Sieradzan et al. J. Chem. Theory Comput., 2012, 8, 4746), the new virtual-bond-angle (θ) potentials now constitute the complete set of physics-based potentials with which to run coarse-grained simulations of systems containing d-amino-acid residues. The ability of the extended UNRES force field to reproduce thermodynamics of polypeptide systems with d-amino-acid residues was tested by comparing the experimentally measured and the calculated free energies of helix formation of model KLALKLALxxLKLALKLA peptides, where x denotes any d- or l- amino-acid residue. The obtained results demonstrate that the UNRES force field with the new potentials reproduce the changes of free energies of helix formation upon d-substitution but overestimate the free energies of helix formation. To test the ability of UNRES with the new potentials to reproduce the structures of polypeptides with d-amino-acid residues, an ab initio replica-exchange folding simulation of thurincin H from Bacillus thuringiensis, which has d-amino-acid residues in the sequence, was carried out. UNRES was able to locate the native α-helical hairpin structure as the dominant structure even though no native sulfide-carbon bonds were present in the simulation.

Improvement of the treatment of loop structures in the UNRES force field by inclusion of coupling between backbone- and side-chain-local conformational states

The UNited RESidue (UNRES) coarse-grained model of polypeptide chains, developed in our laboratory, enables us to carry out millisecond-scale molecular-dynamics simulations of large proteins effectively. It performs well in ab initio predictions of protein structure, as demonstrated in the last Community Wide Experiment on the Critical Assessment of Techniques for Protein Structure Prediction (CASP10). However, the resolution of the simulated structure is too coarse, especially in loop regions, which results from insufficient specificity of the model of local interactions. To improve the representation of local interactions, in this work we introduced new side-chain-backbone correlation potentials, derived from a statistical analysis of loop regions of 4585 proteins. To obtain sufficient statistics, we reduced the set of amino-acid-residue types to five groups, derived in our earlier work on structurally optimized reduced alphabets, based on a statistical analysis of the properties of amino-acid structures. The new correlation potentials are expressed as one-dimensional Fourier series in the virtual-bond-dihedral angles involving side-chain centroids. The weight of these new terms was determined by a trial-and-error method, in which Multiplexed Replica Exchange Molecular Dynamics (MREMD) simulations were run on selected test proteins. The best average root-mean-square deviations (RMSDs) of the calculated structures from the experimental structures below the folding-transition temperatures were obtained with the weight of the new side-chain-backbone correlation potentials equal to 0.57. The resulting conformational ensembles were analyzed in detail by using the Weighted Histogram Analysis Method (WHAM) and Ward's minimum-variance clustering. This analysis showed that the RMSDs from the experimental structures dropped by 0.5 Å on average, compared to simulations without the new terms, and the deviation of individual residues in the loop region of the computed structures from their counterparts in the experimental structures (after optimum superposition of the calculated and experimental structure) decreased by up to 8 Å. Consequently, the new terms improve the representation of local structure.

Local vs global motions in protein folding

It is of interest to know whether local fluctuations in a polypeptide chain play any role in the mechanism by which the chain folds to the native structure of a protein. This question is addressed by analyzing folding and non-folding trajectories of a protein; as an example, the analysis is applied to the 37-residue triple β-strand WW domain from the Formin binding protein 28 (FBP28) (PDB ID: 1E0L). Molecular dynamics (MD) trajectories were generated with the coarse-grained united-residue force field, and one- and two-dimensional free-energy landscapes (FELs) along the backbone virtual-bond angle θ and backbone virtual-bond-dihedral angle γ of each residue, and principal components, respectively, were analyzed. The key residues involved in the folding of the FBP28 WW domain are elucidated by this analysis. The correlations between local and global motions are found. It is shown that most of the residues in the folding trajectories of the system studied here move in a concerted fashion, following the dynamics of the whole system. This demonstrates how the choice of a pathway has to involve concerted movements in order for this protein to fold. This finding also sheds light on the effectiveness of principal component analysis (PCA) for the description of the folding dynamics of the system studied. It is demonstrated that the FEL along the PCs, computed by considering only several critically-placed residues, can correctly describe the folding dynamics.

Interaction of serum amyloid A with human cystatin C--identification of binding sites

Serum amyloid A (SAA) is a multifunctional acute-phase protein whose natural role seems to be participation in many physiologic and pathological processes. Prolonged increased SAA level in a number of chronic inflammatory and neoplastic diseases gives rise to reactive systemic amyloid A amyloidosis, where the N-terminal 76-amino acid residue-long segment of SAA is deposited as amyloid fibrils. Recently, a specific interaction between SAA and the ubiquitous inhibitor of cysteine proteases--human cystatin C (hCC)--has been described. Here, we report further evidence corroborating this interaction, and the identification of the SAA and hCC binding sites in the SAA-hCC complex, using a combination of selective proteolytic excision and high-resolution mass spectrometry. The shortest binding site in the SAA sequence was determined as SAA(86-104), whereas the binding site in hCC sequence was identified as hCC(96-102). Binding specificities of both interacting sequences were ascertained by affinity experiments (ELISA) and by registration of mass spectrum of SAA-hCC complex.

MOLECULAR DYNAMICS SIMULATIONS OF HELIX BUNDLE PROTEINS USING UNRES FORCE FIELD AND ALL-ATOM FORCE FIELD

We have investigated the folding of two helix-bundle proteins, 36-residue Villin headpiece and 56-residue E-domain of Staphylococcal protein A, by combining molecular dynamics (MD) simulations with Coarse-Grained United-Residue (UNRES) Force Field and all-atom force field. Starting from extended structures, each of the proteins was folded to a stable structure within a short time frame using the UNRES model. However, the secondary structures of helices were not well formed. Further refinement using MD simulations with the all-atom force field was able to fold the protein structure into the native-like state with the smallest main-chain root-mean-square deviation of around 3 Å. Detailed analysis of the folding trajectories was presented and the performance of GPU-based MD simulations was also discussed.

Extension of UNRES force field to treat polypeptide chains with D-amino-acid residues

Coarse-grained force fields for protein simulations are usually designed and parameterized to treat proteins composed of natural L-amino-acid residues. However, D-amino-acid residues occur in bacterial, fungal (e.g., gramicidins), as well as human-designed proteins. For this reason, we have extended the UNRES coarse-grained force field developed in our laboratory to treat systems with D-amino-acid residues. We developed the respective virtual-bond-torsional and double-torsional potentials for rotation about the C α · · · C α virtual-bond axis and two consecutive C α · · · C α virtual-bond axes, respectively, as functions of virtual-bond-dihedral angles γ. In turn, these were calculated as potentials of mean force (PMFs) from the diabatic energy surfaces of terminally-blocked model compounds for glycine, alanine, and proline. The potential-energy surfaces were calculated by using the ab initio method of molecular quantum mechanics at the Møller-Plesset (MP2) level of theory and the 6-31G(d,p) basis set, with the rotation angles of the peptide groups about [Formula: see text] and [Formula: see text] used as variables, and the energy was minimized with respect to the remaining degrees of freedom. The PMFs were calculated by numerical integration for all pairs and triplets with all possible combinations of types (glycine, alanine, and proline) and chirality (D or L); however, symmetry relations reduce the number of non-equivalent torsional potentials to 13 and the number of double-torsional potentials to 63 for a given C-terminal blocking group. Subsequently, one- (for torsional) and two-dimensional (for double-torsional potentials) Fourier series were fitted to the PMFs to obtain analytical expressions. It was found that the torsional potentials of the x-Y and X-y types, where X and Y are Ala or Pro, respectively, and a lowercase letter denotes D-chirality, have global minima for small absolute values of γ, accounting for the double-helical structure of gramicidin A, which is a dimer of two chains, each possessing an alternating D-Tyr-L-Tyr sequence, and similar peptides. The side-chain and correlation potentials for D-amino-acid residues were obtained by applying the reflection about the [Formula: see text] plane to the respective potentials for the L-amino-acid residues.

Effects of mutation, truncation, and temperature on the folding kinetics of a WW domain

The purpose of this work is to show how mutation, truncation, and change of temperature can influence the folding kinetics of a protein. This is accomplished by principal component analysis of molecular-dynamics-generated folding trajectories of the triple β-strand WW domain from formin binding protein 28 (FBP28) (Protein Data Bank ID: 1E0L) and its full-size, and singly- and doubly-truncated mutants at temperatures below and very close to the melting point. The reasons for biphasic folding kinetics [i.e., coexistence of slow (three-state) and fast (two-state) phases], including the involvement of a solvent-exposed hydrophobic cluster and another delocalized hydrophobic core in the folding kinetics, are discussed. New folding pathways are identified in free-energy landscapes determined in terms of principal components for full-size mutants. Three-state folding is found to be a main mechanism for folding the FBP28 WW domain and most of the full-size and truncated mutants. The results from the theoretical analysis are compared to those from experiment. Agreements and discrepancies between the theoretical and experimental results are discussed. Because of its importance in understanding protein kinetics and function, the diffusive mechanism by which the FBP28 WW domain and its full-size and truncated mutants explore their conformational space is examined in terms of the mean-square displacement and principal component analysis eigenvalue spectrum analyses. Subdiffusive behavior is observed for all studied systems.

Determination of effective potentials for the stretching of C(α) ⋯ C(α) virtual bonds in polypeptide chains for coarse-grained simulations of proteins from ab initio energy surfaces of N-methylacetamide and N-acetylpyrrolidine

The potentials of mean force (PMF's) for the deformation of the C(α) ⋯ C(α) virtual bonds in polypeptide chains were determined from the diabatic energy surfaces of N-methylacetamide (modeling regular peptide groups) and N-acetylpyrrolidine (modeling the peptide groups preceding proline), calculated at the Møller-Plesset (MP2) ab initio level of theory with the 6-31G(d,p) basis set. The energy surfaces were expressed in the C(α) ⋯ C(α) virtual-bond length (d) and the H-N-C(α) ⋯ C' improper dihedral angle (α) that describes the pyramidicity of the amide nitrogen, or in the C(α)-C'(O)-N-C(α) dihedral angle (ω) and the angle α. For each grid point, the potential energy was minimized with respect to all remaining degrees of freedom. The PMF's obtained from the (d, α) energy surfaces produced realistic free-energy barriers to the trans-cis transition (10 kcal/mol and 13 kcal/mol for the regular and proline peptide groups, respectively, compared to 12.6 - 13.9 kcal/mol and 17.3 - 19.6 kcal/mol determined experimentally for glycylglycine and N-acylprolines, respectively), while those obtained from the (ω, α) energy maps produced either low-quality PMF curves when direct Boltzmann summation was implemented to compute the PMF's or too-flat curves with too-low free-energy barriers to the trans-cis transition if harmonic extrapolation was used to estimate the contributions to the partition function. An analytical bimodal logarithmic-Gaussian expression was fitted to the PMF's, and the potentials were implemented in the UNRES force field. Test Langevin-dynamics simulations were carried out for the Gly-Gly and Gly-Pro dipeptides, which showed a 10(6)-fold increase of the simulated rate of the trans-cis isomerization with respect to that measured experimentally; effectively the same result was obtained with the analytical Kramers theory of reaction rate applied to the UNRES representation of the peptide groups. Application of Kramers' theory to compute the rate constants from the all-atom ab initio energy surfaces of the model compounds studied resulted in isomerization rates close to the experimental values, which demonstrates that the increase of the isomerization rate in UNRES simulations results solely from averaging out the secondary degrees of freedom.

Hidden protein folding pathways in free-energy landscapes uncovered by network analysis

A network analysis is used to uncover hidden folding pathways in free-energy landscapes usually defined in terms of such arbitrary order parameters as root-mean-square deviation from the native structure, radius of gyration, etc. The analysis has been applied to molecular dynamics (MD) trajectories of the B-domain of staphylococcal protein A, generated with the coarse-grained united-residue (UNRES) force field in a broad range of temperatures (270K ≤ T ≤ 325K). Thousands of folding pathways have been identified at each temperature. Out of these many folding pathways, several most probable ones were selected for investigation of the conformational transitions during protein folding. Unlike other conformational space network (CSN) methods, a node in the CSN variant implemented in this work is defined according to the nativelikeness class of the structure, which defines the similarity of segments of the compared structures in terms of secondary-structure, contact-pattern, and local geometry, as well as the overall geometric similarity of the conformation under consideration to that of the reference (experimental) structure. Our previous findings, regarding the folding model and conformations found at the folding-transition temperature for protein A (Maisuradze et al., J. Am. Chem. Soc. 132, 9444, 2010), were confirmed by the conformational space network analysis. In the methodology and in the analysis of the results, the shortest path identified by using the shortest-path algorithm corresponds to the most probable folding pathway in the conformational space network.

Coarse-grained force field: general folding theory

We review the coarse-grained UNited RESidue (UNRES) force field for the simulations of protein structure and dynamics, which is being developed in our laboratory over the last several years. UNRES is a physics-based force field, the prototype of which is defined as a potential of mean force of polypeptide chains in water, where all the degrees of freedom except the coordinates of α-carbon atoms and side-chain centers have been integrated out. We describe the initial implementation of UNRES to protein-structure prediction formulated as a search for the global minimum of the potential-energy function and its subsequent molecular dynamics and extensions of molecular-dynamics implementation, which enabled us to study protein-folding pathways and thermodynamics, as well as to reformulate the protein-structure prediction problem as a search for the conformational ensemble with the lowest free energy at temperatures below the folding-transition temperature. Applications of UNRES to study biological problems are also described.

Simple physics-based analytical formulas for the potentials of mean force of the interaction of amino-acid side chains in water. VI. Oppositely charged side chains

The two-site coarse-grained model for the interactions of charged side chains, to be used with our coarse-grained UNRES force field for protein simulations proposed in the accompanying paper, has been extended to pairs of oppositely charged side chains. The potentials of mean force of four pairs of molecules modeling charged amino-acid side chains, i.e., propionate-n-pentylamine cation (for aspartic acid-lysine), butyrate-n-pentylamine cation (for glutamic acid-lysine), propionate-1-butylguanidine (for aspartic acid-arginine), and butyrate-1-butylguanidine (for glutamic acid-arginine) pairs were determined by umbrella-sampling molecular dynamics simulations in explicit water as functions of distance and orientation, and the analytical expression was fitted to the potentials of mean force. Compared to pairs of like-charged side chains discussed in the accompanying paper, an average quadrupole-quadrupole interaction term had to be introduced to reproduce the Coulombic interactions, and a multistate model of charge distribution had to be introduced to fit the potentials of mean force of all oppositely charged pairs well. The model reproduces all salt-bridge minima and, consequently, is likely to improve the performance of the UNRES force field.

Simple physics-based analytical formulas for the potentials of mean force of the interaction of amino-acid side chains in water. V. Like-charged side chains

A new model of side-chain-side-chain interactions for charged side-chains of amino acids, to be used in the UNRES force-field, has been developed, in which a side chain consists of a nonpolar and a charged site. The interaction energy between the nonpolar sites is composed of a Gay-Berne and a cavity term; the interaction energy between the charged sites consists of a Lennard-Jones term, a Coulombic term, a generalized-Born term, and a cavity term, while the interaction energy between the nonpolar and charged sites is composed of a Gay-Berne and a polarization term. We parametrized the energy function for the models of all six pairs of natural like-charged amino-acid side chains, namely propionate-propionate (for the aspartic acid-aspartic acid pair), butyrate-butyrate (for the glutamic acid-glutamic acid pair), propionate-butyrate (for the aspartic acid-glutamic acid pair), pentylamine cation-pentylamine cation (for the lysine-lysine pair), 1-butylguanidine cation-1-butylguanidine cation (for the arginine-arginine pair), and pentylamine cation-1-butylguanidine cation (for the lysine-arginine pair). By using umbrella-sampling molecular dynamics simulations in explicit TIP3P water, we determined the potentials of mean force of the above-mentioned pairs as functions of distance and orientation and fitted analytical expressions to them. The positions and depths of the contact minima and the positions and heights of the desolvation maxima, including their dependence on the orientation of the molecules were well represented by analytical expressions for all systems. The values of the parameters of all the energy components are physically reasonable, which justifies use of such potentials in coarse-grain protein-folding simulations.

PDZ binding to the BAR domain of PICK1 is elucidated by coarse-grained molecular dynamics

A key regulator of α-amino-3-hydroxy-5-methylisoxazole-4-propionic acid (AMPA) receptor traffic, PICK1 is known to interact with over 40 other proteins, including receptors, transporters and ionic channels, and to be active mostly as a homodimer. The current lack of a complete PICK1 structure determined at atomic resolution hinders the elucidation of its functional mechanisms. Here, we identify interactions between the component PDZ and BAR domains of PICK1 by calculating possible binding sites for the PDZ domain of PICK1 (PICK1-PDZ) to the homology-modeled, crescent-shaped dimer of the PICK1-BAR domain using multiplexed replica-exchange molecular dynamics (MREMD) and canonical molecular dynamics simulations with the coarse-grained UNRES force field. The MREMD results show that the preferred binding site for the single PDZ domain is the concave cavity of the BAR dimer. A second possible binding site is near the N-terminus of the BAR domain that is linked directly to the PDZ domain. Subsequent short canonical molecular dynamics simulations used to determine how the PICK1-PDZ domain moves to the preferred binding site on the BAR domain of PICK1 revealed that initial hydrophobic interactions drive the progress of the simulated binding. Thus, the concave face of the BAR dimer accommodates the PDZ domain first by weak hydrophobic interactions and then the PDZ domain slides to the center of the concave face, where more favorable hydrophobic interactions take over.

Determination of side-chain-rotamer and side-chain and backbone virtual-bond-stretching potentials of mean force from AM1 energy surfaces of terminally-blocked amino-acid residues, for coarse-grained simulations of protein structure and folding. I. The method

In this and the accompanying article, we report the development of new physics-based side-chain-rotamer and virtual-bond-deformation potentials which now replace the respective statistical potentials used so far in our physics-based united-reside UNRES force field for large-scale simulations of protein structure and dynamics. In this article, we describe the methodology for determining the corresponding potentials of mean force (PMF's) from the energy surfaces of terminally-blocked amino-acid residues calculated with the AM1 quantum-mechanical semiempirical method. The approach is based on minimization of the AM1 energy for fixed values of the angles lambda for rotation of the peptide groups about the C(alpha)...C(alpha) virtual bonds, and for fixed values of the side-chain dihedral angles chi, which formed a multidimensional grid. A harmonic-approximation approach was developed to extrapolate from the energy at a given grid point to other points of the conformational space to compute the respective contributions to the PMF. To test the applicability of the harmonic approximation, the rotamer PMF's of alanine and valine obtained with this approach have been compared with those obtained by using a Metropolis Monte Carlo method. The PMF surfaces computed with the harmonic approximation are more rugged and have more pronounced minima than the MC-calculated surfaces but the harmonic-approximation- and MC-calculated PMF values are linearly correlated. The potentials derived with the harmonic approximation are, therefore, appropriate for UNRES for which the weights (scaling factors) of the energy terms are determined by force-field optimization for foldability.

Determination of side-chain-rotamer and side-chain and backbone virtual-bond-stretching potentials of mean force from AM1 energy surfaces of terminally-blocked amino-acid residues, for coarse-grained simulations of protein structure and folding. II. Results, comparison with statistical potentials, and implementation in the UNRES force field

Using the harmonic-approximation approach of the accompanying article and AM1 energy surfaces of terminally blocked amino-acid residues, we determined physics-based side-chain rotamer potentials and the side-chain virtual-bond-deformation potentials of 19 natural amino-acid residues with side chains. The potentials were approximated by analytical formulas and implemented in the UNRES mesoscopic dynamics program. For comparison, the corresponding statistical potentials were determined from 19,682 high-resolution protein structures. The low free-energy region of both the AM1-derived and the statistical potentials is determined by the valence geometry and the L-chirality, and its size increases with side-chain flexibility and decreases with increasing virtual-bond-angle theta. The differences between the free energies of rotamers are greater for the AM1-derived potentials compared with the statistical potentials and, for alanine and other residues with small side chains, a region corresponding to the C(ax)(7) conformation has remarkably low free-energy for the AM1-derived potentials, as opposed to the statistical potentials. These differences probably result from the interactions between neighboring residues and indicate the need for introduction of cooperative terms accounting for the coupling between side-chain rotamer and backbone interactions. Both AM1-derived and statistical virtual-bond-deformation potentials are multimodal for flexible side chains and are topologically similar; however, the regions of minima of the statistical potentials are much narrower, which probably results from imposing restraints in structure determination. The force field with the new potentials was preliminarily optimized using the FBP WW domain (1E0L) and the engrailed homeodomain (1ENH) as training proteins and assessed to be reasonably transferable.