United-residue Research Articles

Molecular Dynamics with the United-Residue Force Field: Ab Initio Folding Simulations of Multichain Proteins

The implementation of molecular dynamics with the united-residue (UNRES) force field is extended to treat multichain proteins. Constant temperature was maintained in the simulations with Berendsen or Langevin thermostats. The method was tested on three alpha-helical proteins (1G6U and GCN4-p1, each with two chains, and 1C94, with four chains). Simulations were carried out for both the isolated single chains and the multichain complexes. The proteins were folded by starting from the extended conformation with random initial velocities and with the chains parallel to each other. No symmetry constraints or structure information were included for the single chains or the multichain complexes. In the case of single-chain simulations, a high percentage of the trajectories (100% for 1G6U, 90% for GCN4-p1, and 80% for 1C94) converged to nativelike structures (assumed as the experimental structure of a monomer in the multichain complex), showing that, for the proteins studied in this work with the UNRES force field, the interactions between chains are not critical for stabilization of the individual chains. In the case of multichain simulations, the native structures of the 1G6U and GCN4-p1 complexes, but not that of 1C94, are predicted successfully. The association of the subunits does not follow a unique mechanism; the monomers were observed to fold both before and simultaneously with their association.

Implementation of a symplectic multiple-time-step molecular dynamics algorithm, based on the united-residue mesoscopic potential energy function

A symplectic multiple-time-step (MTS) algorithm has been developed for the united-residue (UNRES) force field. In this algorithm, the slow-varying forces (which contain most of the long-range interactions and are, therefore, expensive to compute) are integrated with a larger time step, termed the basic time step, and the fast-varying forces are integrated with a shorter time step, which is an integral fraction of the basic time step. Based on the split operator formalism, the equations of motion were derived. Separation of the fast- and slow-varying forces leads to stable molecular dynamics with longer time steps. The algorithms were tested with the Ala(10) polypeptide chain and two versions of the UNRES force field: the current one in which the energy components accounting for the energetics of side-chain rotamers (U(rot)) can lead to numerically unstable forces and a modified one in which the the present U(rot) was replaced by a numerically stable expression which, at present, is parametrized only for polyalanine chains. With the modified UNRES potential, stable trajectories were obtained even when extending the basic time step to 15 fs and, with the original UNRES potentials, the basic time step is 1 fs. An adaptive multiple-time-step (A-MTS) algorithm is proposed to handle instabilities in the forces; in this method, the number of substeps in the basic time step varies depending on the change of the magnitude of the acceleration. With this algorithm, the basic time step is 1 fs but the number of substeps and, consequently, the computational cost are reduced with respect to the MTS algorithm. The use of the UNRES mesoscopic energy function and the algorithms derived in this work enables one to increase the simulation time period by several orders of magnitude compared to conventional atomic-resolution molecular dynamics approaches and, consequently, such an approach appears applicable to simulating protein-folding pathways, protein functional dynamics in a real molecular environment, and dynamical molecular recognition processes.

Replica Exchange and Multicanonical Algorithms with the Coarse-Grained United-Residue (UNRES) Force Field

Three algorithms, namely a Replica Exchange method (REM), a Replica Exchange Multicanonical method (REMUCA), and Replica Exchange Multicanonical with Replica Exchange (REMUCAREM), were implemented with the coarse-grained united-residue force field (UNRES) in both Monte Carlo and Molecular Dynamics versions. The MD algorithms use the constant-temperature Berendsen thermostat, with the velocity Verlet algorithm and variable time step. The algorithms were applied to one peptide (20 residues of Alanine with free ends; ala(20)) and two small proteins, namely an α-helical protein of 46 residues (the B-domain of the staphylococal protein A; 1BDD), and an α+β-protein of 48 residues (the E. Coli Mltd Lysm Domain; 1E0G). Calculated thermodynamic averages, such as canonical average energy and heat capacity, are in good agreement among all simulations for poly-L-alanine, showing that the algorithms were implemented correctly, and that all three algorithms are equally effective for small systems. For protein A, all algorithms performed reasonably well, although some variability in the calculated results was observed whereas, for a more complicated α+β-protein (1E0G), only Replica Exchange was capable of producing reliable statistics for calculating thermodynamic quantities. Finally, from the Replica Exchange molecular dynamics results, we calculated free energy maps as functions of RMSD and radius of gyration for different temperatures. The free energy calculations show correct folding behavior for poly-L-alanine and protein A while, for 1E0G, the native structure had the lowest free energy only at very low temperatures. Hence, the entropy contribution for 1E0G is larger than that for protein A at the same temperature. A larger contribution from entropy means that there are more accessible conformations at a given temperature, making it more difficult to obtain an efficient coverage of conformational space to obtain reliable thermodynamic properties. At the same temperature, ala(20) has the smallest entropy contribution, followed by protein A, and then by 1E0G.

A Hierarchical Multiscale Approach to Protein Structure Prediction: Production of Low‐Resolution Packing Arrangements of Helices and Refinement of the Best Models with a United‐Residue Force Field

A hierarchical, two‐stage approach to ab initio protein structure prediction is presented and applied to four α‐helical proteins. In the first stage, a bank of low‐resolution models is generated using a highly simplified protein representation and energy function, coupled with a Conformation‐Family Monte Carlo (CFMC) search for the energy minimum. For helical proteins, this procedure (referred to as REPACK) produces a set of plausible packed arrangements of the helices, given their positions in the amino acid sequence. Secondary structure prediction methods such as JPRED can be used to provide the secondary structure assignment. In the second stage, these packing arrangements are used as starting points for a new search method (Local Search), based on the Monte Carlo‐with‐Minimization (MCM) algorithm and a united‐residue (UNRES) energy function. The focus of the Local Search is mainly on improving loop conformations and side‐chain positions, with minor modifications to the overall packing of the helices. ...

Kinetic Studies of Folding of the B-domain of Staphylococcal Protein A with Molecular Dynamics and a United-residue (UNRES) Model of Polypeptide Chains

Langevin dynamics is used with our physics-based united-residue (UNRES) force field to study the folding pathways of the B-domain of staphylococcal protein A (1BDD (alpha; 46 residues)). With 400 trajectories of protein A started from the extended state (to gather meaningful statistics), and simulated for more than 35 ns each, 380 of them folded to the native structure. The simulations were carried out at the optimal folding temperature of protein A with this force field. To the best of our knowledge, this is the first simulation study of protein-folding kinetics with a physics-based force field in which reliable statistics can be gathered. In all the simulations, the C-terminal alpha-helix forms first. The ensemble of the native basin has an average RMSD value of 4 A from the native structure. There is a stable intermediate along the folding pathway, in which the N-terminal alpha-helix is unfolded; this intermediate appears on the way to the native structure in less than one-fourth of the folding pathways, while the remaining ones proceed directly to the native state. Non-native structures persist until the end of the simulations, but the native-like structures dominate. To express the kinetics of protein A folding quantitatively, two observables were used: (i) the average alpha-helix content (averaged over all trajectories within a given time window); and (ii) the fraction of conformations (averaged over all trajectories within a given time window) with Calpha RMSD values from the native structure less than 5 A (fraction of completely folded structures). The alpha-helix content grows quickly with time, and its variation fits well to a single-exponential term, suggesting fast two-state kinetics. On the other hand, the fraction of folded structures changes more slowly with time and fits to a sum of two exponentials, in agreement with the appearance of the intermediate, found when analyzing the folding pathways. This observation demonstrates that different qualitative and quantitative conclusions about folding kinetics can be drawn depending on which observable is monitored.

Molecular dynamics with the united-residue model of polypeptide chains. II. Langevin and Berendsen-bath dynamics and tests on model alpha-helical systems.

The implementation of molecular dynamics (MD) with our physics-based protein united-residue (UNRES) force field, described in the accompanying paper, was extended to Langevin dynamics. The equations of motion are integrated by using a simplified stochastic velocity Verlet algorithm. To compare the results to those with all-atom simulations with implicit solvent in which no explicit stochastic and friction forces are present, we alternatively introduced the Berendsen thermostat. Test simulations on the Ala(10) polypeptide demonstrated that the average kinetic energy is stable with about a 5 fs time step. To determine the correspondence between the UNRES time step and the time step of all-atom molecular dynamics, all-atom simulations with the AMBER 99 force field and explicit solvent and also with implicit solvent taken into account within the framework of the generalized Born/surface area (GBSA) model were carried out on the unblocked Ala(10) polypeptide. We found that the UNRES time scale is 4 times longer than that of all-atom MD simulations because the degrees of freedom corresponding to the fastest motions in UNRES are averaged out. When the reduction of the computational cost for evaluation of the UNRES energy function is also taken into account, UNRES (with hydration included implicitly in the side chain-side chain interaction potential) offers about at least a 4000-fold speed up of computations relative to all-atom simulations with explicit solvent and at least a 65-fold speed up relative to all-atom simulations with implicit solvent. To carry out an initial full-blown test of the UNRES/MD approach, we ran Berendsen-bath and Langevin dynamics simulations of the 46-residue B-domain of staphylococcal protein A. We were able to determine the folding temperature at which all trajectories converged to nativelike structures with both approaches. For comparison, we carried out ab initio folding simulations of this protein at the AMBER 99/GBSA level. The average CPU time for folding protein A by UNRES molecular dynamics was 30 min with a single Alpha processor, compared to about 152 h for all-atom simulations with implicit solvent. It can be concluded that the UNRES/MD approach will enable us to carry out microsecond and, possibly, millisecond simulations of protein folding and, consequently, of the folding process of proteins in real time.

Molecular dynamics with the united-residue model of polypeptide chains. I. Lagrange equations of motion and tests of numerical stability in the microcanonical mode.

The Lagrange formalism was implemented to derive the equations of motion for the physics-based united-residue (UNRES) force field developed in our laboratory. The C(alpha)...C(alpha) and C(alpha)...SC (SC denoting a side-chain center) virtual-bond vectors were chosen as variables. The velocity Verlet algorithm was adopted to integrate the equations of motion. Tests on the unblocked Ala(10) polypeptide showed that the algorithm is stable in short periods of time up to the time step of 1.467 fs; however, even with the shorter time step of 0.489 fs, some drift of the total energy occurs because of momentary jumps of the acceleration. These jumps are caused by numerical instability of the forces arising from the U(rot) component of UNRES that describes the energetics of side-chain-rotameric states. Test runs on the Gly(10) sequence (in which U(rot) is not present) and on the Ala(10) sequence with U(rot) replaced by a simple numerically stable harmonic potential confirmed this observation; oscillations of the total energy were observed only up to the time step of 7.335 fs, and some drift in the total energy or instability of the trajectories started to appear in long-time (2 ns and longer) trajectories only for the time step of 9.78 fs. These results demonstrate that the present U(rot) components (which are statistical potentials derived from the Protein Data Bank) must be replaced with more numerically stable functions; this work is under way in our laboratory. For the purpose of our present work, a nonsymplectic variable-time-step algorithm was introduced to reduce the energy drift for regular polypeptide sequences. The algorithm scales down the time step at a given point of a trajectory if the maximum change of acceleration exceeds a selected cutoff value. With this algorithm, the total energy is reasonably conserved up to a time step of 2.445 fs, as tested on the unblocked Ala(10) polypeptide. We also tried a symplectic multiple-time-step reversible RESPA algorithm and achieved satisfactory energy conservation for time steps up to 7.335 fs. However, at present, it appears that the reversible RESPA algorithm is several times more expensive than the variable-time-step algorithm because of the necessity to perform additional matrix multiplications. We also observed that, because Ala(10) folds and unfolds within picoseconds in the microcanonical mode, this suggests that the effective (event-based) time unit in UNRES dynamics is much larger than that of all-atom dynamics because of averaging over the fast-moving degrees of freedom in deriving the UNRES potential.

Ab initio simulations of protein-folding pathways by molecular dynamics with the united-residue model of polypeptide chains

We report the application of Langevin dynamics to the physics-based united-residue (UNRES) force field developed in our laboratory. Ten trajectories were run on seven proteins [PDB ID codes 1BDD (alpha; 46 residues), 1GAB (alpha; 47 residues), 1LQ7 (alpha; 67 residues), 1CLB (alpha; 75 residues), 1E0L (beta; 28 residues), and 1E0G (alpha+beta; 48 residues), and 1IGD (alpha+beta; 61 residues)] with the UNRES force field parameterized by using our recently developed method for obtaining a hierarchical structure of the energy landscape. All alpha-helical proteins and 1E0G folded to the native-like structures, whereas 1IGD and 1E0L yielded mostly nonnative alpha-helical folds although the native-like structures are lowest in energy for these two proteins, which can be attributed to neglecting the entropy factor in the current parameterization of UNRES. Average folding times for successful folding simulations were of the order of nanoseconds, whereas even the ultrafast-folding proteins fold only in microseconds, which implies that the UNRES time scale is approximately three orders of magnitude larger than the experimental time scale because the fast motions of the secondary degrees of freedom are averaged out. Folding with Langevin dynamics required 2-10 h of CPU time on average with a single AMD Athlon MP 2800+ processor depending on the size of the protein. With the advantage of parallel processing, this process leads to the possibility to explore thousands of folding pathways and to predict not only the native structure but also the folding scenario of a protein together with its quantitative kinetic and thermodynamic characteristics.

A united residue force-field for calcium-protein interactions.

United-residue potentials are derived for interactions of the calcium cation with polypeptide chains in energy-based prediction of protein structure with a united-residue (UNRES) force-field. Specific potentials were derived for the interaction of the calcium cation with the Asp, Glu, Asn, and Gln side chains and the peptide group. The analytical expressions for the interaction energies for each of these amino acids were obtained by averaging the electrostatic interaction energy, expressed by a multipole series over the dihedral angles not considered in the united-residue model, that is, the side-chain dihedral angles chi and the dihedral angles lambda for the rotation of peptide groups about the C(alpha)...C(alpha) virtual-bond axes. For the side-chains that do not interact favorably with calcium, simple excluded-volume potentials were introduced. The parameters of the potentials were obtained from ab initio quantum mechanical calculations of model systems at the Restricted Hartree-Fock (RHF) level with the 6-31G(d,p) basis set. The energy surfaces of pairs consisting of Ca(2+)-acetate, Ca(2+)-propionate, Ca(2+)-acetamide, Ca(2+)-propionamide, and Ca(2+)-N-methylacetamide systems (modeling the Ca(2+)-Asp(-), Ca(2+)-Glu(-), Ca(2+)-Asn, Ca(2+)-Gln, and Ca(2+)-peptide group interactions) at different distances and orientations were calculated. For each pair, the restricted free energy (RFE) surfaces were calculated by numerical integration over the degrees of freedom lost when switching from the all-atom model to the united-residue model. Finally, the analytical expressions for each pair were fitted to the RFE surfaces. This force-field was able to distinguish the EF-hand motif from all potential binding sites in the crystal structures of bovine alpha-lactalbumin, whiting parvalbumin, calbindin D9K, and apo-calbindin D9K.

Optimization of the UNRES Force Field by Hierarchical Design of the Potential-Energy Landscape. 2. Off-Lattice Tests of the Method with Single Proteins

We describe the application of our recently proposed method of hierarchical optimization of the protein energy landscape to optimize our off-lattice united-residue (UNRES) force field using single training proteins. First, the IgG-binding domain from streptococcal protein G (PDB code 1IGD) was treated; earlier attempts to use this protein to optimize the force field by optimizing the energy gap and Z score between the nativelike and non-native structures failed. The structure of this protein consists of an N-terminal antiparallel β-hairpin, a middle α-helix, and a C-terminal antiparallel β-hairpin, these elements being referred to as β 1 , α 2 , and β 3 , respectively, with the two hairpins forming a parallel β-sheet packed against the α-helix. In our earlier study, one of these elements was assumed to form at level 1, two at level 2, and three at level 3, and higher levels corresponded to the proper packing of two or more elements. This approach resulted in a structure with the wrong packing of the β-sheet, and attempts at further optimization failed. We therefore tried a hierarchy scheme that corresponds to the sequence of folding events deduced from NMR experiments. In this scheme, level I corresponds to structures with either β 3 or α 2 , level 2 to structures with both β 3 and α 2 , level 3 to structures with β 3 , α 2 , and the N-terminal strand packed against α 2 (with β 1 still not fully formed), and level 4 to structures with β 1 , α 2 , and β 3 , with β 3 being packed to β 1 , which also implies the packing of β 1 and β 3 against α 2 . This optimization was successful and resulted in a reasonably transferable force field that led to well-foldable proteins. This corroborates the conclusion from our model on-lattice studies (Liwo, A.; Arlukowicz, P.; Oldziej, S.; Czaplewski, C.; Makowski, M.; Scheraga, H. A. J. Phys. Chem. B 2004, 108, 16918) that a proper design of the structural hierarchy is of crucial importance to the foldability with the resulting potential-energy function. Moreover, in the off-lattice approach, the design of the hierarchy also appears to be important to the success of the optimization procedure itself. The next series of calculations was carried out with the LysM domain from the E. coli 1E0G (a + β) protein, which is smaller than 1IGD. In this case, no experimental information about the folding pathway is available; nevertheless, we were able to deduce the appropriate hierarchy by a trial-and-error method. The resulting force field performed worse in tests on a + β- and β-proteins than that derived on the basis of 1IGD with a correct hierarchy, which suggests that the structure of the I IGD protein encodes more structure-determining interactions common to all proteins than the 1E0G protein does. For 1E0G, we also attempted to carry out a single energy gap and Z-score optimization; this effort resulted in an unsearchable force field. (The nativelike structures could not be found by a global search, although they were the lowest in energy). Technical details of the method, including the maintenance of proper secondary structure and a method to classify structure, are also described.

The protein folding problem: global optimization of the force fields.

The evolutionary development of a theoretical approach to the protein folding problem, in our laboratory, is traced. The theoretical foundations and the development of a suitable empirical all-atom potential energy function and a global optimization search are examined. Whereas the all-atom approach has thus far succeeded for relatively small molecules and for alpha-helical proteins containing up to 46 residues, it has been necessary to develop a hierarchical approach to treat larger proteins. In the hierarchical approach to single- and multiple-chain proteins, global optimization is carried out for a simplified united residue (UNRES) description of a polypeptide chain to locate the region in which the global minimum lies. Conversion of the UNRES structures in this region to all-atom structures is followed by a local search in this region. The performance of this approach in successive CASP blind tests for predicting protein structure by an ab initio physics-based method is described. Finally, a recent attempt to compute a folding pathway is discussed.

Ab initio structure prediction of two α‐helical oligomers with a multiple‐chain united‐residue force field and global search

A hierarchical methodology for ab initio structure prediction is extended to treat oligomeric proteins. Modifications are made to a united-residue (UNRES) force field and a Conformational Space Annealing (CSA) global search method. The computational cost of including additional chains and the increase in speed from symmetry optimizations are evaluated. The native structures of two oligomeric proteins from the CASP3 exercise, the retro-GCN4 leucine zipper and the synthetic domain-swapped dimer, were identified as the lowest-energy families resulting from the search of the proteins when rotational symmetry was imposed. Additional searches in different symmetries and oligomerization states were carried out, and the results indicate some problems in the thoroughness of the search and in the search of packing arrangements if symmetry constraints are not imposed.

Ab initio folding of multiple-chain proteins.

Our previous methodology for ab initio prediction of protein structure is extended here to treat multiple-chain proteins. This involved modification of our united-residue (UNRES) force field and our Conformational Space Annealing (CSA) Global Optimization procedure. Good results have been obtained for both a four- and a three-helix protein from the CASP3 exercise.

Development of Physics-Based Energy Functions that Predict Medium-Resolution Structures for Proteins of the α, β, and α/β Structural Classes

The development of three physics-based energy functions (force fields), designed to simulate the restricted free energy of proteins of the α, β, and α/β structural classes, is described. Each force field corresponds to a particular weighting of the united-residue (UNRES) interactions defined in earlier work.1-6 To find the optimal weights for the α, β, and α/β force fields, both the Z-score and energy gap of the native versus nonnative structures are minimized simultaneously for four benchmark proteins: 1pou (for the α force field), 1tpm (for the β force field), and 1bdd and betanova (for the α/β force field). The simultaneous minimization was carried out by using a novel Monte Carlo method, Vector Monte Carlo (VMC). For α-helical proteins, another weighting of the UNRES interactions (denoted as the α0 force field) was developed; this fourth force field is described in a companion publication (Lee, J. et al. J. Phys. Chem. B 2001, 105, 7291). The structural implications of the final weights of the four f...