Optimization of the GBMV2 implicit solvent force field for accurate simulation of protein conformational equilibria

Abstract

Accurate treatment of solvent environment is critical for reliable simulations of protein conformational equilibria. Implicit treatment of solvation, such as using the generalized Born (GB) class of models arguably provides an optimal balance between computational efficiency and physical accuracy. Yet, GB models are frequently plagued by a tendency to generate overly compact structures. The physical origins of this drawback are relatively well understood, and the key to a balanced implicit solvent protein force field is careful optimization of physical parameters to achieve a sufficient level of cancellation of errors. The latter has been hampered by the difficulty of generating converged conformational ensembles of non-trivial model proteins using the popular replica exchange sampling technique. Here, we leverage improved sampling efficiency of a newly developed multi-scale enhanced sampling technique to re-optimize the generalized-Born with molecular volume (GBMV2) implicit solvent model with the CHARMM36 protein force field. Recursive optimization of key GBMV2 parameters (such as input radii) and protein torsion profiles (via the CMAP torsion cross terms) has led to a more balanced GBMV2 protein force field that recapitulates the structures and stabilities of both helical and β-hairpin model peptides. Importantly, this force field appears to be free of the over-compaction bias, and can generate structural ensembles of several intrinsically disordered proteins of various lengths that seem highly consistent with available experimental data. © 2017 Wiley Periodicals, Inc.

Keywords: continuum electrostatics; enhanced sampling; generalized Born; intrinsically disordered proteins; nonpolar solvation; protein folding.

Figure 1

PMFs of the backbone hydrogen-bonding interaction between a modified alanine dipeptide dimer (see insert) in TIP3P water and GBMV2 implicit solvent with two different sets of atomic input radii (see Table 1). Details of the dimer and calculation of the TIP3P PMF were described previously.^,

Figure 2

(A) Optimized CMAP for GBMV2 for all non-Proline residues, and (B) the change made to the original CHARMM36 CMAP.

Figure 3

Average residue helicity of peptide (AAQAA)₃ at 270 K derived from various segments of the control and folding simulations. The 80–140 ns fragments were used for data analyses. The average helicities were shown in parentheses. The experimental values were taken from Shalongo et al. 1994 (See Ref 49).

Figure 4

Probability distributions of the number of native backbone hydrogen bonds for (a) GB1p, (b) GB1m1, and (c) GB1m3. These distributions were calculated from structure ensembles extracted from the various segments of MSES simulations at T = 270 K. The number in the parentheses is the ratio of the folded population, identified as structures with ≥5 native backbone hydrogen bonds formed.

Figure 5

Potential energy versus Cα RMSD for structure ensembles sampled at 270K from various MSES simulaitons of (A) GB1p, (B) GB1m1, and (C) GB1m3. Only conformationsl sampled during the last 30 ns are included.

Figure 6

Residue helicity profiles of (A) KID, (B) ACTR, and (C) the RS peptide, calculated from structures sampled at 300 K during the second half of MSES simulations.

Figure 7

Distributions of end-to-end distance for (A) KID, (B) ACTR, and (C) the RS peptide, calculated from structures sampled at 300 K during the second half of MSES simulations.

Figure 8

(A) R_g distribution of the RS peptide and (B) theoretical SAXS curves calculated from the structure ensemble. The experimental result is show in grey trace. The structure ensembles include all snapshots sampled at 298 K during the second half of MSES simulations.