Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral angles

Abstract

While the quality of the current CHARMM22/CMAP additive force field for proteins has been demonstrated in a large number of applications, limitations in the model with respect to the equilibrium between the sampling of helical and extended conformations in folding simulations have been noted. To overcome this, as well as make other improvements in the model, we present a combination of refinements that should result in enhanced accuracy in simulations of proteins. The common (non Gly, Pro) backbone CMAP potential has been refined against experimental solution NMR data for weakly structured peptides, resulting in a rebalancing of the energies of the α-helix and extended regions of the Ramachandran map, correcting the α-helical bias of CHARMM22/CMAP. The Gly and Pro CMAPs have been refitted to more accurate quantum-mechanical energy surfaces. Side-chain torsion parameters have been optimized by fitting to backbone-dependent quantum-mechanical energy surfaces, followed by additional empirical optimization targeting NMR scalar couplings for unfolded proteins. A comprehensive validation of the revised force field was then performed against data not used to guide parametrization: (i) comparison of simulations of eight proteins in their crystal environments with crystal structures; (ii) comparison with backbone scalar couplings for weakly structured peptides; (iii) comparison with NMR residual dipolar couplings and scalar couplings for both backbone and side-chains in folded proteins; (iv) equilibrium folding of mini-proteins. The results indicate that the revised CHARMM 36 parameters represent an improved model for the modeling and simulation studies of proteins, including studies of protein folding, assembly and functionally relevant conformational changes.

Figure 1

Helix formation in Ac-(AAQAA)₃-NH₂. (A) Fraction helix per residue estimated from experiment at 300 K (black lines) compared with average fraction helix calculated from REMD simulations with C22/CMAP^, (solid blue symbols), C36 (solid red symbols), Amber ff99SB (open blue symbols) and Amber ff99SB* (open red symbols) using the replica closest in temperature to 300 K. (B) Average carbonyl carbon chemical shifts calculated with SPARTA+ from the simulations in (A) compared with experimental shifts. Color scheme as in (A).

Figure 2

RMS Differences in relative energies between MM and QM potential energy surfaces. RMS Differences were calculated for relative energies less than 12.0 kcal/mol above the global minimum for the three 1D or 2D surfaces. Higher-energy cutoffs of 20.0 kcal/mol for Arg and Lys or 25.0 kcal/mol for Asp, Glu and Hsp were used for the charged amino acids.

Figure 3

Sampling of backbone φ/ψ torsion angles in the central residues of Ala₃, Ala₅, Ala₇, Val₃, and Gly₃. Colors indicate relative free energies according to color bar given on the right.

Figure 4

Sampling of φ/ψ torsion angles in the central residues of Ala₅ with different force fields: the new force field (C36), the previous C22/CMAP force field (C22), Amber ff99SB, OPLS-AA, Gromos 53a6 and Amber ff99SB*.

Figure 5

A) RMS difference from experimental structure and B) Inter-helical angle of the dimeric coiled-coil 1U0I as a function of simulation time. RMS differences are for the backbone (N, C^α, C, O) or non-hydrogen atoms of both helices following least squares alignment of the backbone atoms in both helices. The N- and C-terminal residues were excluded from the analysis. For the inter-helical angles vectors defining the helical axes were calculated using the non-terminal residue C^α atoms using the approach of Chothia et al.. Horizontal lines represent the inter-helical angles from the 20 models generated in the NMR study.

Figure 6

Ramanchandran plot for backbone sampling. Crystallographic φ/ψ values (triangles) are overlaid onto probability distributions obtained from the MD simulation of full unit cells. Probabilities calculated with snapshots between 5ns-40ns at 5ps intervals.

Figure 7

RMSD from experimental structure for long runs of folded proteins. (A) Ubiquitin, (B) Hen Lysozyme, (C) GB3 and (D) BPTI. Upper plots show results with C22/CMAP, lower plots with C36. Various dotted lines indicated RMSD values of 0.5, 1.0, 1.5 and 2.0 Å.

Figure 8

β-hairpin folding test. H^α chemical shifts from experiment (red symbols) and calculated from simulation using SPARTA+ (black symbols) for: (A) the GB1 hairpin (folded structure from simulation shown in (B)); and (C) the designed GB1m3 hairpin (structure shown in (D)). Error bars on calculated data are the typical RMSD between experiment and H^α shifts predicted by SPARTA+. Experimental data at 280 K from Fesinmeyer et al; simulation data taken from 278 K replica.