Accurate calculation of side chain packing and free energy with applications to protein molecular dynamics

Abstract

To address the large gap between time scales that can be easily reached by molecular simulations and those required to understand protein dynamics, we present a rapid self-consistent approximation of the side chain free energy at every integration step. In analogy with the adiabatic Born-Oppenheimer approximation for electronic structure, the protein backbone dynamics are simulated as preceding according to the dictates of the free energy of an instantaneously-equilibrated side chain potential. The side chain free energy is computed on the fly, allowing the protein backbone dynamics to traverse a greatly smoothed energetic landscape. This computation results in extremely rapid equilibration and sampling of the Boltzmann distribution. Our method, termed Upside, employs a reduced model involving the three backbone atoms, along with the carbonyl oxygen and amide proton, and a single (oriented) side chain bead having multiple locations reflecting the conformational diversity of the side chain's rotameric states. We also introduce a novel, maximum-likelihood method to parameterize the side chain interactions using protein structures. We demonstrate state-of-the-art accuracy for predicting χ1 rotamer states while consuming only milliseconds of CPU time. Our method enables rapidly equilibrating coarse-grained simulations that can nonetheless contain significant molecular detail. We also show that the resulting free energies of the side chains are sufficiently accurate for de novo folding of some proteins.

Fig 1. Six step inner loop of Upside calculation.

The side chain potential enters into the integration step simply as a complicated, many-body energy function that may be treated with standard techniques of molecular simulations.

Fig 2. Error in the position as a function of the number of side chain states, resulting from a decomposition of rotamer states into coarse-grained states.

The table summarizes the number of states chosen for each amino acid type. The relative uncertainty is the positional uncertainty for each number of states divided by the accuracy at three states. One, three, or six rotamer states are used, depending on the residue type. For residues without a rotatable χ₂, such as valine, only three states are needed. The time to compute the pairwise interactions and solve for the free energy scales roughly as the number of coarse rotamer states squared, so the use of fewer coarse states is preferred. Ile, Leu and Lys are the three residues with rotatable χ₂ where only 3 states are assigned.

Fig 3. Fragment of protein G with associated interaction graph (R_cutoff = 7Å).

A pair of residues is assigned a connection whenever their side chain beads are within R_cutoff for any side chain states.

Fig 4. Example of optimized coarse states for arginine overlaid on the PDB distribution of the rotamer angles χ₁ and χ₂.

Each of the six coarse states contains only a single fine state that has high probability, so that the variance of dihedral angles within each coarse state is small.

Fig 5. Coordinates and potentials used for side chain interactions.

Top. For each residue, a reference backbone structure is aligned to the N, C_α, and C atomic coordinates. This alignment creates a reference frame to establish the position and direction of the side chain bead. The two side chain beads x₁ and x₂ for a pair of residues establishes three coordinates, the distance r and angles θ₁ and θ₂. Bottom. Example of distance-dependent potential, unif(r₁₂), after training, between the side chain and backbone residues. These interactions cutoff at 5 Å while the side chain-side chain interactions cutoff at 7 Å. The thin lines describe the V_radial of oxygen (red) and hydrogen (blue), and the thick lines describe V_radial + V_angular for the same interactions.

Fig 6. Comparison of χ₁ prediction accuracy for Upside and SCWRL4, ordered by Upside accuracy.

The “PDB χ₁ frequency” line represents the accuracy of the NDRD rotamer library without any interactions; this library is used in both Upside and SCWRL4.

Fig 7. Comparison of the accuracy of predicting side chains as well as cpu running time.

For all programs, time spent reading the protein structure and writing the results is excluded from the running time to focus on the cost of solving for the side chain positions. For Upside (10 Å cutoff), all side chain interactions with backbone or other side chains are cutoff at 10 Å.

Fig 8. Accuracy of MD simulations for three proteins at variable backbone hydrogen bond strength.

Results with and without side chain-backbone interactions are presented.

Fig 9. Closest structure to native protein (lowest C_α-RMSD) at optimal hydrogen bonding strength, with and without backbone-side chain hydrogen bonds.

For alpha3D, BBA and homeodomain with backbone-side chain hydrogen bonds, the optimal hydrogen bond strength and lowest C_α-RMSD are -2.0, -2.0 and -1.9 RT, and 3.6, 1.7 and 3 Å, respectively. Without these hydrogen bonds, the corresponding values are -2.0, -1.8 and -1.6 RT; and 2.7, 1.8, and 2.5 Å. Blue is the native structure and red is the simulation.