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

Abstract

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.

Fig. 1

Illustration of the local coordinate system of a virtual-bond side chain. Dot-dashed lines indicate the virtual bonds. The ${\mathrm{C}}_i^{\alpha }$ atom is at the origin of the reference system, the x axis of the reference system is the bisector of the virtual-bond angle θ, the y axis lies in the plane of the three C^α atoms, is perpendicular to the x axis and directed from ${\mathrm{C}}_i^{\alpha }$ to ${\mathrm{C}}_{i+1}^{\alpha }$. All three axes (x, y, and z) form a right-handed reference system. r_SC is the vector pointing from ${\mathrm{C}}_i^{\alpha }$ to the geometric center of the side chain. The angle α′ is the planar angle between the bisector of the θ angle and the C^α ⋯ SC vector and the angle β′ is the angle of clockwise rotation of the C^α ⋯ SC virtual-bond axis about the bisector of the θ angle from the plane of the three C^α atoms, taking the position of SC in the plane closer to ${\mathrm{C}}_{i+1}^{\alpha }$ (with positive y) as reference (β′ = 0°).

Fig. 2

Side-chain-rotamer potential surfaces plotted in the angles α′ and β′ (cf. Figure 1) of selected amino-acid residues: Ala (a), Pro (b), Val (c), Phe (d), Glu (e), and Arg (f) computed from the respective non-adiabatic AM1 energy surfaces, using the harmonic-approximation procedure described in the accompanying paper (left panels), fitted to the AM1-derived surfaces using eq. (3) (middle panels), and derived from the PDB as statistical potentials (right panels) for the virtual-bond-valence angle θ = 90°. The “South Pole” (α′ = 180°) is the point in the center of each panel, and the “North Pole” (α′ = 0°) is in the middle of the left and of the right vertical side of the rectangle. The parallels (lines of constant α′) are the distorted circles centered about the “South Pole” (except the parallel corresponding to α′ = 90°, which is a square centered at the “South Pole”) and semicircles centered about the “North Pole”, except those corresponding to α′ = 90° (the “Equator”) which constitute two half-squares. The meridians are the lines intersecting the parallels and running between the “North Pole” and the “South Pole”. The parallels and the meridians are each spaced 15°. The horizontal half-line going from the center of each panel to the right corresponds to β′ = 0° and that going from the center to the left corresponds to β′ = 180°; β′ increases when rotating clockwise about the “South Pole”. For all residues, the data to calculate the statistical potentials were taken from residues not involved in regular secondary structures. The free-energy color scale is shown in the small left-bottom panel.

Fig. 3

Side-chain-rotamers potential surfaces plotted in angles α′ and β′ (cf. Figure 1) of selected amino-acid residues: Ala (a), Pro (b), Val (c), Phe (d), Glu (e), and Arg (f) for virtual-bond-angle θ = 140°. See Figure 2 for detailed description. The data to calculate the statistical potentials were taken from residues not involved in regular secondary structures except proline for which all data were taken because of the small number of proline residues with θ ≈ 140°.

Fig. 4

Representative rotamers of the Phe side chain for (a) θ ≈ 90° and (b) θ ≈ 140°. The rotamers of the side chain referred to in the text are color coded: red – rotamer 1, orange – rotamer 2, green – rotamer 3. The drawings were done with MOLMOL. The colors of the rotamers have been chosen to match their free energies, as shown in Figures 2d and 3d, middle panels.

Fig. 5

The U_bond curves of the side chains of selected amino-acid residues: Val (a), Ile (b), and Arg (c). Left panels: the potentials of mean force calculated from non-adiabatic AM1 energy surfaces with the use of the harmonic-approximation procedure described in the accompanying paper (filled circles) and by fitting eq. (8) to these surfaces (solid lines). Right panels: the corresponding statistical potentials (solid lines).

Fig. 6

Stereoscopic views of the C^α traces of the average structures of the most native-like clusters of the proteins studied obtained in MREMD simulations with the UNRES force field incorporating the physics-based U_b potentials introduced in ref. and the physics-based U_rot and U_bond potentials introduced in this work (grey thin sticks) superposed on the C^α traces of the corresponding experimental structures (black thick sticks): (a) 1ENH, (b) 1E0L, (c) 1BDD, (d) 1GAB, (e) 1LQ7, (f) 1E0G, (g) 1PGA. For 1PGA only the fragments encompassing the middle helix and C-terminal β-hairpin (from residues 20 to 56) are superposed. The RMSD's between the computed and the experimental structures are listed in the ρ_ave column of Table 4. The drawings were done with MOLMOL.

Fig. 6

Stereoscopic views of the C^α traces of the average structures of the most native-like clusters of the proteins studied obtained in MREMD simulations with the UNRES force field incorporating the physics-based U_b potentials introduced in ref. and the physics-based U_rot and U_bond potentials introduced in this work (grey thin sticks) superposed on the C^α traces of the corresponding experimental structures (black thick sticks): (a) 1ENH, (b) 1E0L, (c) 1BDD, (d) 1GAB, (e) 1LQ7, (f) 1E0G, (g) 1PGA. For 1PGA only the fragments encompassing the middle helix and C-terminal β-hairpin (from residues 20 to 56) are superposed. The RMSD's between the computed and the experimental structures are listed in the ρ_ave column of Table 4. The drawings were done with MOLMOL.

Fig. 7

Plots of total energy vs. number of MD steps in microcanonical simulations of decaalanine with the time step δt = 4.89 fs in a run with (a) old U_rot, U_bond, and U_b potentials and (b) and with the new potentials.