COFFDROP: A Coarse-Grained Nonbonded Force Field for Proteins Derived from All-Atom Explicit-Solvent Molecular Dynamics Simulations of Amino Acids

Abstract

We describe the derivation of a set of bonded and nonbonded coarse-grained (CG) potential functions for use in implicit-solvent Brownian dynamics (BD) simulations of proteins derived from all-atom explicit-solvent molecular dynamics (MD) simulations of amino acids. Bonded potential functions were derived from 1 μs MD simulations of each of the 20 canonical amino acids, with histidine modeled in both its protonated and neutral forms; nonbonded potential functions were derived from 1 μs MD simulations of every possible pairing of the amino acids (231 different systems). The angle and dihedral probability distributions and radial distribution functions sampled during MD were used to optimize a set of CG potential functions through use of the iterative Boltzmann inversion (IBI) method. The optimized set of potential functions-which we term COFFDROP (COarse-grained Force Field for Dynamic Representation Of Proteins)-quantitatively reproduced all of the "target" MD distributions. In a first test of the force field, it was used to predict the clustering behavior of concentrated amino acid solutions; the predictions were directly compared with the results of corresponding all-atom explicit-solvent MD simulations and found to be in excellent agreement. In a second test, BD simulations of the small protein villin headpiece were carried out at concentrations that have recently been studied in all-atom explicit-solvent MD simulations by Petrov and Zagrovic (PLoS Comput. Biol.2014, 5, e1003638). The anomalously strong intermolecular interactions seen in the MD study were reproduced in the COFFDROP simulations; a simple scaling of COFFDROP's nonbonded parameters, however, produced results in better accordance with experiment. Overall, our results suggest that potential functions derived from simulations of pairwise amino acid interactions might be of quite broad applicability, with COFFDROP likely to be especially useful for modeling unfolded or intrinsically disordered proteins.

Figure 1

Translational diffusion coefficients and association kinetics of amino acids calculated from all-atom MD simulations. (A) Translational diffusion coefficients of 12 amino acids calculated from MD and compared to experimental data from Longsworth. Each symbol is labeled using the one letter amino acid code. (B) Plot showing the correlation between the calculated effective association rate constant (k_on) and the sum of the individual amino acid translational diffusion coefficients. Each symbol represents a different amino acid pair.

Figure 2

Radial distribution functions of amino acid pairs and amino acid dendrogram calculated from all-atom MD. (A) Plot showing 231 Cα–Cα g(r) functions. (B) Dendrogram created by performing agglomerative clustering on the calculated correlation coefficients of the 231 Cα–Cα g(r) functions shown in A. (C) Plot showing 231 g(r) functions calculated using only the closest distance between any pair of heavy atoms.

Figure 3

Derivation of COFFDROP bonded potential functions using the IBI method. (A) Plot showing the error in the angle probability distributions obtained from BD simulations as a function of IBI iteration number for the amino acids arginine, alanine and tryptophan. (B) Same as A but showing results for dihedral probability distributions. (C) Comparison of the angle probability distributions obtained from MD (lines) with those obtained from BD (circles) for tryptophan. Each color represents a different angle. (D) Same as C but showing results for dihedral probability distributions. (E) Comparison of an example angle potential function (Ace–Cα–Nme for tryptophan) obtained from using IBI (blue) with that obtained from noniterative Boltzmann inversion of the MD probability distribution (red). (F) Same as E but showing an example dihedral potential function (Cγ–Cβ–Cα–Nme for tryptophan).

Figure 4

Derivation of COFFDROP nonbonded potential functions using the IBI method. (A) Plot showing the error in the nonbonded g(r) functions obtained from BD simulations as a function of IBI iteration number for the ile–leu (green circles), glu–arg (yellow upward triangles), and tyr–trp (red downward triangles) systems. (B) Comparison of binding affinities calculated from the Cα–Cα g(r) functions from MD (x-axis) and BD (y-axis). The green, yellow, and red symbols represent the ile–leu, glu–arg, and tyr–trp systems, respectively; the blue symbols represent the other 228 systems.

Figure 5

Example nonbonded potential functions. (A) Plot comparing the Ace–Ace nonbonded potential function of the val–val system obtained from IBI (blue) with that obtained from noniterative Boltzmann inversion (red) of the MD g(r) function. (B) Same as A but for the Cα–Cα interaction; (C) same as A but for the Cβ–Cβ interaction. (D) Plot comparing the Ace–Ace g(r) of the val-val system obtained from MD (black circles) with that obtained from BD using COFFDROP (blue line). (E) Same as D but for for the Cα–Cα interaction; (F) same as D but for the Cβ–Cβ interaction.

Figure 6

Clustering of alanine, leucine, asparagine, and tryptophan solutions in MD and BD. The plots show the fraction of solute molecules that are members of clusters of various sizes. Blue circles represent results from MD, green upward triangles represent results from BD using COFFDROP, and red downward triangles represent results from BD using steric nonbonded potentials.

Figure 7

Clustering of alanine, leucine, asparagine, and tryptophan solutions at concentrations of 200 and 300 mg/mL in MD and BD. Same as Figure 6 but showing results for much higher concentrations.

Figure 8

Clustering of villin headpiece solutions at a 9.2 mM concentration in BD. (A) Plot shows the fraction of villin headpiece molecules that are members of clusters of various sizes. Blue circles represent results using a 1.0 scaling factor with COFFDROP’s nonbonded potential functions, green upward triangles represent results using a 0.9 scaling factor, yellow downward triangles represent results using a 0.8 scaling factor, and red squares represent results using a 0.8 scaling factor and starting from a structure in which the villin molecules were already aggregated into a trimer and pentamer. (B) Image showing aggregated villin molecules obtained at the end of a 200 ns BD simulation using a 1.0 scaling factor. Each color represents a different villin molecule. (C) Image showing villin molecules at the end of a 200 ns BD simulation using a 0.8 scaling factor.

Figure 9

Spatial disposition of the Cδ pseudoatom of tryptophan in MD and BD. (A–C) Red contours show preferred locations of the Cδ pseudoatom of a tryptophan molecule interacting with a second tryptophan molecule (shown in black) sampled from MD; each of the panels A–C shows the same image viewed from a different orientation. (D–F) Same as panels A–C, respectively, but showing results from BD.