PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Chem Theory Comput. Author manuscript; available in PMC 2010 August 11.
Published in final edited form as:
J Chem Theory Comput. 2009 August 11; 5(8): 2115–2124.
PMCID: PMC2725330
NIHMSID: NIHMS133289

A Transferable Coarse Grain Non-bonded Interaction Model For Amino Acids

Abstract

The large quantity of protein sequences being generated from genomic data has greatly outpaced the throughput of experimental protein structure determining methods and consequently brought urgency to the need for accurate protein structure prediction tools. Reduced resolution, or coarse grained (CG) models, have become a mainstay in computational protein structure prediction perfoming among the best tools available. The quest for high quality generalized CG models presents an extremely challenging yet popular endeavor. To this point, a CG based interaction potential is presented here for the naturally occurring amino acids. In the present approach, three to four heavy atoms and associated hydrogens are condensed into a single CG site. The parameterization of the site-site interaction potential relies on experimental data thus providing a novel approach that is neither based on all-atom (AA) simulations nor experimental protein structural data. Specifically, intermolecular potentials, which are based on Lennard-Jones (LJ) style functional forms, are parameterized using thermodynamic data including surface tension and density. Using this approach, an amino acid potential dataset has been developed for use in modeling peptides and proteins. The potential is evaluated here by comparing the solvent accessible surface area (SASA) to AA representations and ranking of protein decoy data sets provided by Decoys ’R’ Us. The model is shown to perform very well compared to other existing prediction models for these properties.

Keywords: coarse grain, molecular dynamics, amino acids, proteins, peptides, surface tension, decoys, solvent accessible surface area

I. INTRODUCTION

With the rate at which genomic data has become available over the last several years, the need to accurately predict protein structures from the amino acid sequence has become paramount. Protein structure prediction has employed many approaches to attempt this feat with a few noted here.116 Many of these models are based on a reduced resolution, or coarse grained (CG) representation of the protein.13,5,1730 Such models are used because they provide a means to expand the capabilities of existing computational resources. However, this enhanced efficiency comes necessarily at the cost of reduced resolution in the description of the system since typically the details of several atoms are essentially averaged into a single CG interaction site. When this reduction of resolution is introduced, it is crucial to accurately capture the phenomenon of interest for the specific system under investigation or the resulting CG description is simply a toy model capable of providing only generic behavior.

The history of CG models for amino acids dates back over three decades at least, with the pioneering work of Levitt & Warshel, where they used reduced resolution models to investigate protein dynamics.1 Since then, much work has been put forth to develop reduced resolution models to overcome the computational barriers of studying proteins. Much of this work has been built on knowledge based models which take advantage of the existing structural data acquired from experimental techniques including x-ray crystallography and nmr spectroscopy. Tanaka and Scheraga were the first to develop knowledge based contact potentials from the frequency of contacts between residues in known structures.2 This approach was then extended by many groups to include for example distance and direction dependence.3,5,1720 Others such as Huang et al. have broken from this paradigm to use an approach that does not rely on known protein structural data other than the general nature of the packing of hydrophobic and hydrophilic residues in protein structures.9,31

Further, the use of CG protein models for molecular dynamics (MD) simulations has helped extend the capabilities of current computational resources allowing access to much larger temporal and spatial scales.3247 In particular, the MARTINI model has become a widely used CG model with an array of applications to protein and membrane simulations.3335 The later approach employees a LJ potential for non-bonded interactions with parameters in a tabulated form based on the interaction types (hydrophobic/hydrophobic, hydrophilic/hydrophobic, etc.). Many other groups have pushed CG models to a still higher level with approaches such as force matching,3840 and other schemes which have enjoyed success.4153

Herein we present the application of a recently proposed methodology to the development of CG non-bonded interaction potentials for amino acids.44,4851 The current model is essentially a distance dependent potential which distinguishes the 20 naturally occurring residues to a high level such that only a few residues are modeled with the same parameters. Molecular simulations are used to parameterize CG sites that are used here in a non-dynamical representation of proteins. Experimental thermodynamic data is used as the target for parameterization of the CG site interactions. This approach builds on previous work by Nielsen et al.52 who used a similar approach to obtain CG parameters for a series of alkanes and Shinoda et al.49,50 who extended the CG parameters to polyethylene glycol surfactants. Related approaches have been reported recently but with noteable differences.54,55 The work by Basdevant et al.54 uses a hybrid 1/r6 repulsive term and Gaussian attractive term that is paramterized using AA force fields. The work by Han et al.55 uses a thermodynamic based parameterization approach but with a higher resolution mapping.

Herein, the current model is used to predict the native structures from several protein decoy sets with a level of success on par with knowledge based potentials. Finally, although not explicitly demonstrated here, the non-bonded CG model presented could be coupled with an appropriate intramolecular forcefield and used for MD simulations of peptides and proteins in a spirit similar to the MARTINI force field. Indeed, an initial application of the present approach was previously utilized in an MD simulation investigation of peptide nanoring self assembly.32 In Section II, the methods are discussed including the force field details and parameterization. Section III presents results and discussion including evaluation of the model through solvent accessible surface area (SASA) calculations and the ranking of protein decoys. Finally, Section IV closes with the conclusions.

II. METHODS

A. Coarse Grain Model

Mainly for ease of implementation, the present CG model employs LJ style non-bonded potential functions.49 For the models developed herein, all CG beads interact via a LJ(9-6) potential given in Equation 1. The CG sites are charged and this is also include in the electrostatic contribution given by Equation 2.

υ(r)96=274ε(σ9r9σ6r6)
(1)

υ(r)elec=14πε0q1q2r
(2)

The choice of prefactor for the LJ function is selected such that υ(σ) = 0 and ε is the minimum energy. The choice of LJ functional form (for example the 9-6 in Equation 1) is essentially an adjustable parameter used in the fitting procedure. We have previously explored various options including 6-4, 8-4, 10-4, 9-6 and 12-4.49 For alkane interactions the choice of the 9-6 functional form was validated by comparison of CG liquid structure to AA MD results. The adoption of the 9-6 functional form reflects a desire to maintain consistency with existing CG models.49,50 A simple truncation, implemented at a distance of rcut = 15Å, is employed for the long range cut-off, with no smoothing or shifting. This length is sufficient to avoid gross artifacts resulting from the discontinuity; however the cut-off length, rcut does affect the thermodynamic properties and so is treated as a CG fitting parameter. The parameters in the LJ function are fixed by reproducing thermodynamic data. For bulk solutions, both ε and σ can be unambiguously fixed with a combination of density and surface tension. The cross interactions arising from non-identical CG sites can be generated using the combining rules given by Equation 3 and Equation 4.

σab=σaa+σbb2
(3)

εab=εaaεbb
(4)

Here, εaa and σaa represent the self interaction values for ε and σ values and εAB and σAB represent those for the cross interactions. Finally, the electrostatic interactions are calculated using a effective dielectric constant of 80 and without employing a cutoff for electrostatic interactions. The dielectric constant is given below in the form of scaled charges with the value of 0.1118.

Although the intramolecular forcefield is not discussed or used in detail here, it was necessary to employ bonds and angles in the development of the non-bonded interactions for PHE, TRP and TYR (throughout the paper, the 3-letter amino acid labels will be used) sidechains. For those cases, the intramolecular interactions are modeled via a harmonic potentials given by Equation 5 and Equation 6,

V(r)bond=Kb(rr0)2
(5)

V(r)angle=Ka(θθ0)2,
(6)

where r0 and kb represent the equilibrium bond length and force constant for bonds and θ0 and ka represent the equilibrium angle value and force constant for bends.

B. Molecular Dynamics Simulations

CG MD simulations were performed using the LAMMPS code developed at Sandia National Laboratory and extended by our group (now a part of the standard LAMMPS release) to implement our CG models.56 An integration timestep of up to at least 25fs can be used to evaluate the non-bonded interactions. Intramolecular degrees of freedom have to be dealt with on a case by case basis as the force constant will determine the timestep necessary. For these cases, the multi-timestep integrator rRESPA can be used to seperate the various degrees of freedom.57

Electrostatic interactions were calculated using the Particle-Particle Particle-Mesh method implemented in the LAMMPS MD code.5861 Non-bonded interactions are excluded for bonded CG sites seperated by two or less bonds (1–2, 1–3). Interactions between CG sites seperated by 3 or more bonds are included without scaling. For the sake of evaluating these exclusions, the connectivity (bonds) of the model was considered although the bond energy was not calculated in the potential evaluation. Only the non-bonded van der Waals and electrostatic energy were used. The bulk system MD simulations included roughly 500 CG sites each. For surface tension simulations, roughly 50ns simulation time was used. Systems were set up by equilibrating with an isothermal-isobaric (NPT) simulation. The equilibrated system box dimension was then extended in the z-direction creating a vacuum region large enough such that no interactions extended through the vacuum via the periodic boundary conditions to the opposite side. The extended system was then used in a MD run in a canonical (NVT) simulation. The surface tension, γ was then calculated via

γ=Lz2[Pzz(Pxx+Pyy)2]
(7)

where Lz is the box dimension in the z-direction and Paa is the isotropic pressure tensor for each cartesian direction where a=x,y,z.62 Density calculations were performed with NPT simulations averaged over roughly 2ns simulation time. Finally, all MD simulations were performed at a temperature of 303.15K, unless otherwise specified.

C. Solvent Accessible Surface Area

One property of interest in proteins is the SASA which is proportional to protein solvation free energy and represents an important property in protein research and design. In the work of Lee & Richards, it was shown that the hydrophobic residues tend to have a larger reduction of SASA when going from extended chains to the native conformation compared to the polar residues.63 Succinctly, the SASA is the area of the surface traced out by rolling a probe sphere, representing the solvent, over the surface of the protein. Any area accessible to the probe that does not require the probe to overlap with neighboring atoms is considered solvent accessible. The SASA was calculated here via the implementation in VMD.64 For the AA level calculations, a probe radius of 1.4Å and VMD default atomic radii for the protein atoms were used. For the CG SASA calculations, a probe radius of 2.5Å was used. This value was calculated from the minimum of the potential energy function, r where V (r)/r = 0, for the CG water-water interaction. CG site radii for the amino acids were calculated in a similar fashion. Note that the minima of the LJ functional forms used herein are not equal to the standard LJ (12-6) minimum of r = 21/6σ, but are given by r = 1.5σ and r = 3σ for the LJ (9-6) and LJ (12-4), respectively.

D. Parameterization

In the CG model presented here, the mapping includes roughly three to four heavy atoms and adjacent hydrogens per CG site. Each amino acid is divided into a backbone and sidechain section. The mappings for all of the naturally occurring amino acids are given in Figure 1, Figure 2 and Figure 3. All amino acids in this model (including GLY) share the same backbone representation with the exception of ALA. The standard backbone CG site includes the carbonyl carbon and oxygen, nitrogen, alpha carbon and related hydrogens. GLY is of course represented by the single standard backbone CG site. Due to its small sidechain, the ALA amino acid is mapped to a single site containing the backbone and sidechain. The mappings for the backbone sites are shown in Figure 1 with standard CG backbone site (labeled GBB in Table III used for all amino acids except ALA on the left and the ALA CG model (label ABB in Table III given on the right.

FIG. 1
Shown here are the CG mapping for the backbone sites. The standard backbone CG site (GBB) is shown on the left and the ALA CG site (ABB) is shown on the right.
FIG. 2
Shown are the mappings for the amino acid sidechains. The label here correspond to the labels used in the interaction potential Table III.
FIG. 3
The CG mappings for the charged residues are shown here. ASP, GLU, LYS2 and ARG2 are all modeled with the same interaction parameters. LYS1 and LYS2 are modeled with the VAL CG site paramters. Note that HIS is included in this Figure although it is not ...

Most sidechains are represented by a single CG site. The exceptions to this are the PHE, TYR, LYS and ARG sidechains, which are modeled with two CG sites, and TRP which is modeled with three. Finally, Figure 2 and Figure 3 shows the mapping for all of the neutral and charged sidechains, respectively.

The sidechain analogue molecules used in the model development of the neutral amino acids are given in Table I. To parameterize the LJ9-6 potentials for the self interactions of the CG sites (SER-SER, LEU-LEU, etc.) the surface tension and density were reproduced in bulk solutions of the representative sidechain analogue where possible. For example, this was not possible for the TRP sidechain analogue 3-methylindole due to the complexity (too many parameters to fix simultaneously with only two observables) and the fact that it is a solid at room temperature with a melting point of around 370 K. Instead, the CG sites for this sidechain were developed in a stepwise process (see below). On the other hand, having slightly higher melting point did not alter the approach for ASN and GLN sidechains. A slightly higher temperature was used in the simulations for the parameterization of this CG sites. Strictly speaking, the CG models are parameterized at a specific temperature and thus not transferable to different temperatures; however, there is a small range of temperature which will not alter the properties to a large degree. Results for the surface tension and density of each sidechain analogue is compared to experimental values in Table II.65 The parameters for the VAL sidechain are from previous work.49 The parameters for LEU and ILE were developed using the analogues 2,5-dimethylhexane and 3,4-dimethylhexane respectively since the true analogues are gases at room temperature. Equilibrium bond lengths of 3.87 Å and 3.22 Å with force constant of 5.0 kcal/mol where used for each, respectively.

TABLE I
Sidechain analogue molecules.
TABLE II
Thermodynamic data from CG simulations compared with experiment.

For PHE, TYR and TRP, it was necessary to implement a bond between the CG sites in order to parameterize the overall CG assembly. For PHE, and TRP the chosen minimum energy bond lengths, 2.5Å and the force constant of 150 kcal/mol, were used for all bonds. For TYR, the minimum energy bond length, 2.9 Å and the force constant of 150 kcal/mol were used. The choice of a two-site CG model for PHE and TYR was used to maintain consistency with the mapping of roughly 3 to 4 heavy atoms per CG site. This is admittedly a gross overestimation of these sidechains and the implications arising are discussed below. The PHE sidechain was developed from the CG sites for p-xylene (XYL) and benzene (BEN). Symmetric two site models were made for each of these two molecules and parameterized to reproduce the respective surface tension and density. These sites were then combined to produce the PHE sidechain model. Figure 4 shows the CG sites and corresponding AA representations (without hydrogens for clarity) used in the stepwise process. The blue and red spheres represent the BEN and XYL CG sites, respectively. In the PHE sidechain CG model (far right of Figure 4), the XYL CG site, represented by the red sphere, is labeled PHE1 and the BEN CG site, represented by the blue sphere, is labeled PHE2 in the interaction database (see Figure 4). The combination rules, Equation 3 and Equation 4, from the methods section were used to develop the XYL-BEN cross interaction. The result of using this approach can be seen in the reproduction of the bulk properties for the PHE sidechain model in Table II. The XYL CG site was then combined with a new site to create the cresol (TYR sidechain) CG model in a similar fashion. The labels in the TYR CG model are TYR1, modeled with the XYL CG site and TYR2 modeled with the new −OH containing cresol CG site (see Figure 4). Finally, the benzene CG site was used to develop a CG model for the TRP sidechain.

FIG. 4
The models used to develop the phenyl based rings (PHE, TYR and TRP) are shown here. The CG model for benzene (left and blue) and p-xylene (middle and red) were combined to produce the PHE CG model (right). These sites were also used in the TYR and TRP ...

The analogue molecule 3-methylindole was used to represent the TRP sidechain. The CG mapping for the TRP sidechain is shown in Figure 2. The new CG sites TRP1 and TRP3 represent propene and allylamine, respectively. Each of these molecules is represented by a single distinct CG site that was parameterized in a similar fashion to the amino acids with sidechains represented by single CG sites. In a stepwise process the TRP sidechain was built up first by combining the BEN and ALL CG sites to create the CG model for aniline. The bond was modeled with the equilibrium bond length and force constant stated above. This process yielded the BEN-ALL cross interaction from the bulk aniline simulations. Then, the BEN and ALL beads were combined with the PPE CG site to create the TRP sidechain CG model. The combination rules were employed to fix the interactions BEN-PPE and ALL-PPE.

For the charged residues (ASP, GLU, LYS and ARG) a primitive model was developed including charges. Figure 3 shows the CG mapping for the charged residues. It should be noted that although HIS is included here, it is not modeled with a charge in the current model. The CG sites labeled ASP, GLU, LYS2 and ARG2 are identical in the current model. These sites were parameterized to reproduce the surface tension of guanidium salt solutions at various concentrations. A charge of 0.1118 (with ε0 = 1 in Equation 2 scaled down from unity to compensate for the solvent screening that would take place normally but is absent due to the CG water model being represented by a single charge-less non-polar spherical site) is included on the sites ASP, GLU, LYS2 and ARG2. All other sites are modeled without charge. The addition of a scaled charge is part of a recent development designed for modeling charged surfactants. Details will be part of a forthcoming manuscript. As stated above, the current model will not distinguish ASP from GLU or LYS from ARG. For LYS and ARG, the charged CG sites, labeled LYS2 and ARG2 (modeled with a charge), are combined with a CG site equivalent to the VAL sidechain CG site (and modeled neutral) and labeled LYS1 and ARG1, respectively in Figure 3, to create the two bead CG site for these sidechains.

III. RESULTS AND DISCUSSIONS

A. Surface Area

Evaluation of forcefields for complex molecules such as proteins is a difficult task. Indeed it is an ongoing process where many systems must be explored in order to provide a good representation of possible interactions that can arise. Nonetheless, here we set out to provide at least a preliminary level of evaluation. For protein-protein interactions, steric effects arising from sidechain packing can play an important role in determining the spatial arrangements. Therefore the relative sizes of amino acid sidechains is a key property to preserve in any model of amino acids. Further, it is understood that CG modeling often distorts the representation of molecules especially when mapping to spherical sites. To evaluate the effects of this mapping, the SASA was calculated for each amino acid sidechain and backbone unit from a set of PDB databank structures with less than 30% homology. For these calculations, each sidechain or backbone unit was calculated as a solitary group of atoms (excluding all other atoms in the structure) to give a SASA for that residue sidechain or backbone type (GLY or ALA) taking into account only the various conformations of each subunit. From these calculations, the radii were calculated for each subunit and are shown in Table III. The results are compared to the radii calculated for the CG model sidechain units as described in the methods section. As shown in Table III, the values compare well suggesting conservation of SASA with the mapping from AA to CG. It is important to note that the radii for the CG model were not based on all atom or PDB calculations, but arise from the parameterization process as described above. To further evaluate the SASA predictive ability of the CG model, the SASA for full proteins were calculated. Table IV shows the results of the SASA calculations for several proteins comparing the CG results to those from the AA level. The results show that despite the fact that the CG model uses a bulky solvent (combining 3 water molecules) the SASA is consistent in the CG and AA models. Finally, it should be mentioned that calculation of SASA is not straightforward and many things, such as the probe radius, can affect the results. Measurements made here were done as consistently as possible using standard techniques and radii for the probe and atomic radii. Therefore, it is assumed that this is at least a reasonable comparison of the SASA.

TABLE IV
Comparison of SASA for CG and AA models. The SASA was calculated for the CG and AA representations with a probe radius of 2.5Å and 1.4Å, respectively.

B. Predicting Native Structures

The ranking of structures in decoy sets such as those provided by Decoys R Us is a useful test for evaluating amino acid potentials and has become a somewhat standard evaluation tool.66 The decoy sets are composed of native protein structures with decoys for each native structure generated through various techniques. The primary use of the protein decoy sets is to test a models ability to distinguish the native structure from the non-native decoy structures. Here we evaluated the current CG model using five decoy sets (4state_reduced, fisa, fisa_casp3, lmds, lattice_ssfit) provided by Decoys R Us (http//dd.stanford.edu).14,6669 For comparison, the MARTINI force field was also used to rank the decoy sets.35 For each structure, the CG model was mapped onto the AA model by placing the CG sites at the center of mass of the representative heavy atoms. The non-bonded potential presented here (and the MARTINI non-bonded potential) was then used to calculate the potential energy of each structure. These energies were then ranked from lowest to highest. The observable is the ranking of the potential energy of the native structure compared to the decoys. A value of 1 indicates that the native structure was predicted to have the lowest energy (is the most stable) of the set. The higher the value the worse the ability of the model to identify the native structure from the decoys. The ranking of the native structures for each set are shown in Table V. The results for the model presented here are shown in the columns labeled CG and those for the MARTINI model are labeled MARTINI.35

TABLE V
Results of ranking native structures for decoy sets from Decoys ’R’ Us. The lower the value, the better the performance of the model. A value of “1” indicates that the model predicts the native structure to be the most ...

During the evaluation of the decoy sets, one issue had to be addressed in an ad-hoc fashion. Disulfide bonds arising between CYS residues had to be treated as a special case. The parameterization the CYS sidechain site used bulk methylmercaptan solutions leading to an effective size which is much larger than the center of mass distance that arises when the CYS residues are involved in a covalent bond. To address this issue, sigma for the CYS-CYS interaction was shortened to a value of 2.4Å to represent the CYS disulfide bond. However, it is important to note that the cross terms (CYS-LEU, CYS-SER, etc.) generated with the combination rules used the original bulk solution sigma value of 4.16Å. To test the effect of the CYS-CYS interactions, the structures were reevaluated with all of the CYS-CYS interactions removed. The results did not show significant differences which indicates that the use of the modified sigma value is reasonable for the structures evaluated here. However, to make a more general potential, this issue will have to be dealt with and the CYS-CYS will need to be categorized as either covalent or non-covalent. Typically CYS-CYS interactions involved in a disulfide bond would be excluded from the non-bonded interactions; however, the current setup does not take into account these bonds and thus they had to be dealt with with a reduced sigma value. For the MARTINI model, the calculation was repeated with all of the CYS-CYS interactions turned off and showed no significant difference in the results.

The results in Table V show that the potential is extremely effective at picking the native structure from the decoy sets. The present model ranked the native structure as the most stable in 68% of the decoy sets and ranked the native structure as one of the top 3 most stable structures in 77% of the decoy sets. For the sake of comparison, other methods are shown including MARITINI (mentioned above and not designed for this type of evaluation), DFIRE1, Rosetta, ModPipe-Pair, ModPipe-Comb and DOPE.14,67,7075 The values for these methods were taken directly from the references of Zhou & Zhou73,74 and Shen & Sali.13 Another evaluation of the quality of ranking is the Z-score which indicates the models discriminatory ability. We define the Z-score as:

Z=EENσ
(1)

where EN is the energy of the native structure, < E > is the the average energy of all structures in a given set (for example 4state_reduced : 1ctf) and σ is the standard deviation of the distribution of energy values for a given set. The Z-scores for each set are given in Table V. To prevent the Z-scores from being skewed by structures with unreasonably high energies, any decoy structure with a positive potential energy value was not included in the Z-score calculation. The percentage of structures that had negative potential energy values as determined by the current model are listed in Table V under the heading “%NEG”. It should be noted that the fisa_casp3 set had a high percentage of decoy structures with a positive potential energy according the current model and so the Z-scores from those sets must be judged accordingly. Further, the lattice_ssfit had a modest amount of structures evaluated with positive potential energies. Finally, all of the native structures had a negative potential energy. For the sets in which the native structure was ranked as the most stable structure, the average Z-score was 3.13. For comparison, Z-scores for the DFIRE model have been included in Table V. The average Z-score from the DFIRE model for all of the systems tested here was 5.21.

The 4state_reduced structures 3icb, 4pti, 4rxn and the 1dtk structure from the lmds set are ranked very poorly by the current model. However, closer examination of these structures indicates that a substantial portion of the bad contacts arise from interactions involving phenyl based residues. The problem results from the mapping of the PHE and TYR sidechains to two site CG models. The planar nature of these sidechains is critical for the spatial arrangements they occupy and must be conserved as much as possible in the CG model. The CG mapped structure for 1dtk native conformation is shown in Figure 5. In this figure, the TYR1 sidechain site is shown as a red sphere and the ALA backbone site is shown as a blue sphere. It is easy to see here how the two site spherical representations of the phenyl based sidechains leads to cases where bad contacts are generated by the CG model. This arrangement is not present in all decoys and indeed the decoys with the lowest energy values do not have this conformation in which the TYR residue is in close proximity to other sites. The choice of a two site mapping for these residues was made to stay consistent with the mapping of roughly 3 to 4 heavy atoms. However, this issue will need to be addressed in future refinement of the model.

FIG. 5
Shown is the CG representation of the 1dtk native structure from the lmds decoy set. The close contact is highlighted for the TYR2 site (red sphere) and the ALA site (blue sphere).

As mentioned above, the electrostatic residues were crudely parameterized. The current parameter set is unable to distinguish the residues ASP from GLU or LYS from ARG. This is demonstrated by the fact that the electrostatics make only small improvements to the results (4state_reduced-1ctf went from the number 3 ranked structure to number 1; lmds-1igd went from the number 2 ranked structure to number 1; others saw small or no improvement with no structures showing detriment with inclusion of electrostatics). However, it is possible that the structures studied here have conformations that are not largely dependent on electrostatic effects. Further, the percentage of charged residues is typically low in a protein, so a small impact on results is not totally unexpected.

Finally, the backbone atoms of proteins pose a particularly difficult problem as the directionality of hydrogen bonding makes the backbone an extremely anisotropic collection of atoms. This problem can be seen in the decoy sets evaluated here where high energy contacts are formed upon mapping to a CG representation. In these structure, the all atom representation shows hydrogen bonding interactions between the backbone sites leading to a stabilization of that conformation. Further, it should be mentioned that in the model used here, the backbone was not specifically parameterized but took advantage of the CG sites developed for the ASN and GLN sidechains. This is a good first approximation to the interactions of the standard backbone and ALA backbone interaction sites; however, there is room for improvement of the interaction parameters of these sites. Also, although a nonpolar spherical CG site is a crude approximation to the backbone, the results shown here suggest that the approximation can perform reasonably well. Nonetheless, improvements to the performance of the backbone CG sites could possibly begin by addition of a dipole term such as that recently employed by Cascella et al.76

IV. CONCLUSIONS

A systematic approach has been used to develop a CG non-bonded interaction potential for amino acids. The model uses a reduced representation in which roughly 3 to 4 heavy atoms (non-hydrogen) and adjacent hydrogens are mapped to a single CG. The conservation of size in the model is reassuring for the protein-protein interaction predictive ability and for the future possibility of inclusion of a solvation contribution to the fold stability. The work here also demonstrates the consistency that can be maintained with parameterization techniques that are not based directly on protein structures. Further, it should be reiterated that no protein structural data or preconceived protein structural characteristics were used in the parameterization of this model. The encouraging results shown in the ranking of the decoy sets suggest that the current approach can be useful as a protein structure predictive tool. The decoy analysis results also highlights the importance of accurately treating the morphology of sidechains especially the phenyl based residues (PHE,TRP & TYR). As stated above, although the model is not presented here as suitable for MD simulations, there is nothing to prevent its implementation into a MD framework. Work on an intramolecular potential is currently in progress and will be presented in a forthcoming manuscript. In addition to this, further refinement of the non-bonded interactions will be a focal point of future work. The results presented herein suggest that the present approach is capable of yielding useful CG models.

Acknowledgments

The authors thank NSF, NIH and Procter & Gamble for generous support. Also, the authors would like to thank Axel Kohlmeyer for implementation of the CG model into the LAMMPS MD code. Finally, thanks to Vincenzo Carnevale for useful discussion. R.D. was supported by a NSF Biological Informatics Postdoctoral Fellowship under Grant No. DBI-0532800.

References

1. Levitt M, Warshel A. Nature. 1975;253:694–698. [PubMed]
2. Tanaka S, Scheraga HA. Macromolecules. 1976;9:945–950. [PubMed]
3. Miyazawa S, Jernigan RL. Macromolecules. 1985;18:534–552.
4. Treptow WL, Barbosa MAA, Garcia LG, de Araujo AFP. Proteins: Struct. Funct. Genet. 2002;49 167-18. [PubMed]
5. Sippl MJ. J. Mol. Biol. 1990;213:859–883. [PubMed]
6. Bahar I, Atilgan AR, Erman B. Folding & Design. 1997;2:173–181. [PubMed]
7. Ueda Y, Taketomi H, Go N. Biopolymers. 1978;17:1531–1548.
8. Dill KA, Ozkan SB, Shell MS, Weikl TR. Ann. Rev. Biophys. 2008;37:289–316. [PMC free article] [PubMed]
9. Huang ES, Subbiah S, Levitt M. J. Mol. Biol. 1995;252:709–720. [PubMed]
10. Schueler-Furman O, Wang C, Bradley P, Misura K, Baker D. Science. 2005;310:638–642. [PubMed]
11. Moult J. Curr. Opin. Struct. Biol. 2005;15:285–289. [PubMed]
12. Zhang Y. Curr. Opin. Struct. Biol. 2008;18:342–348. [PMC free article] [PubMed]
13. Shen MY, Sali A. Protein Sci. 2006;15:2507–2524. [PubMed]
14. Simons KT, Kooperberg C, Huang E, Baker D. J. Mol. Biol. 1997;268:209–225. [PubMed]
15. Reva BA, Finkelstein AV, Sanner M, Olson AJ, Skolnick J. Protein Eng. 1997;10:1123–1130. [PubMed]
16. Petrey D, Honig B. Molecular Cell. 2005;20:811–819. [PubMed]
17. Miyazawa S, Jernigan RL. J. Mol. Biol. 1996;256:623–644. [PubMed]
18. Miyazawa S, Jernigan RL. Proteins: Struct. Funct. Genet. 1999;36:357–369. [PubMed]
19. Sippl MJ. J. Comput.-Aided Mol. Des. 1993;7:473–501. [PubMed]
20. Sippl MJ, Weitckus S. Proteins: Struct. Func. Genet. 1992;13:258–271. [PubMed]
21. Kolinski A, Skolnick J. Polymer. 2004;45:511–524.
22. Clementi C. Curr. Opin. Struct. Biol. 2008;18:10–15. [PubMed]
23. Heath AP, Kavraki LE, Clementi C. Proteins: Struct. Funct. Bioinf. 2007;68:646–661. [PubMed]
24. Matysiak S, Clementi C. J. Mol. Biol. 2006;363:297–308. [PubMed]
25. Liwo A, Pincus MR, Wawak RJ, Rackovsky S, Scheraga HA. Protein Science. 1993;2:1697–1714. [PubMed]
26. Liwo A, Oldziej S, Pincus MR, Wawak RJ, Rackovsky S, Scheraga HA. Journal of Computational Chemistry. 1997;18:849–873.
27. Liwo A, Pincus MR, Wawak RJ, Rackovsky S, Oldziej S, Scheraga HA. Journal of Computational Chemistry. 1997;18:874–887.
28. Liwo A, Kazmierkiewicz R, Czaplewski C, Groth M, Oldziej S, Wawak RJ, Rackovsky S, Pincus MR, Scheraga HA. Journal of Computational Chemistry. 1998;19:259–276.
29. Yap EH, Fawzi NL, Head-Gordon T. Proteins-Structure Function and Bioinformatics. 2008;70:626–638. [PMC free article] [PubMed]
30. Khatun J, Khare SD, Dokholyan NV. Journal of Molecular Biology. 2004;336:1223–1238. [PubMed]
31. Huang ES, Subbiah S, Tsai J, Levitt M. J. Mol. Biol. 1996;257:716–725. [PubMed]
32. Khurana E, DeVane R, Kohlmeyer A, Klein ML. Nano Lett. 2008;8:3626. [PMC free article] [PubMed]
33. Marrink SJ, de Vries AH, Mark AE. J. Phys. Chem. B. 2004;108:750–760.
34. Marrink SJ, Risselada HJ, Yefimov S, Tieleman DP, de Vries AH. J. Phys. Chem. B. 2007;111:7812–7824. [PubMed]
35. Monticelli L, Kandasamy SK, Periole X, Larson RG, Tieleman DP, Marrink SJ. J. Chem. Theory Comput. 2008;4:819–834.
36. Bond PJ, Holyoake J, Ivetac A, Khalid S, Sansom MSP. J. Struct. Biol. 2007;157:593–605. [PubMed]
37. Treptow W, Marrink SJ, Tarek M. J. Phys. Chem. B. 2008;112:3277–3282. [PubMed]
38. Izvekov S, Voth GA. J. Phys. Chem. B. 2005;109:2469–2473. [PubMed]
39. Noid WG, Chu JW, Ayton GS, Voth GA. J. Phys. Chem. B. 2007;111:4116–4127. [PMC free article] [PubMed]
40. Zhou J, Thorpe IF, Izvekov S, Voth GA. Biophys. J. 2007;92:4289–4303. [PubMed]
41. Tozzini V. Curr. Opin. Struct. Biol. 2005;15:144–150. [PubMed]
42. Tozzini V, McCammon JA. Chem. Phys. Lett.s. 2005;413:123–128.
43. Masella M, Borgis D, Cuniasse P. J. Comput. Chem. 2008;29:1707–1724. [PubMed]
44. Klein ML, Shinoda W. Science. 2008;321:798–800. [PubMed]
45. Aksimentiev A, Schulten K. Proceedings of the National Academy of Sciences of the United States of America. 2004;101:4337–4338. [PubMed]
46. Arkhipov A, Yin Y, Schulten K. Biophys. J. 2008;95:2806–2821. [PubMed]
47. Shih AY, Arkhipov A, Freddolino PL, Schulten K. J. Phys. Chem. B. 2006;110:3674–3684. [PMC free article] [PubMed]
48. Shelley JC, Shelley MY, Reeder RC, Bandyopadhyay S, Klein ML. J. Phys. Chem. B. 2001;105:4464–4470. 26.
49. Shinoda W, DeVane R, Klein ML. Mol. Simul. 2006;33:27–36.
50. Shinoda W, DeVane R, Klein ML. Soft Matter. 2008;4:2454–2462.
51. Bhargava BL, DeVane R, Klein ML, Balasubramanian S. Soft Matter. 2007;3:1395–1400.
52. Nielsen SO, Lopez CF, Srinivas G, Klein ML. J. Chem. Phys. 2003;119:7043–7049.
53. Shih AY, Freddolino PL, Arkhipov A, Schulten K. J. Struct. Biol. 2007;157:579–592. [PubMed]
54. Basdevant N, Borgis D, Ha-Duong T. J. Phys. Chem. B. 2007;111:9390–9399. [PubMed]
55. Han W, Wan CK, Wu YD. J. Chem. Theory Comput. 2008;4:1891–1901.
56. Plimpton S. J. Comput. Phys. 1995;117:1–19.
57. Tuckerman ME, Berne BJ, Martyna GJ. J. Chem. Phys. 1991;94:6811–6815.
58. Hockney R, Eastwood J. Computer Simulation Using Particles. Bristol: IOP; 1988.
59. Darden T, York D, Pedersen L. J. Chem. Phys. 1993;98:10089–10092.
60. Deserno M, Holm C. J. Chem. Phys. 1998;109:7678–7693.
61. Deserno M, Holm C. J. Chem. Phys. 1998;109:7694–7701.
62. Allen M, Tildesley D. Computer Simulation of Liquids. Oxford Science; 1987.
63. Lee B, Richards FM. J. Mol. Biol. 1971;55:379. [PubMed]
64. Humphrey W, Dalke A, Schulten K. J. Mol. Graphics. 1996;14:33. [PubMed]
65. L YC. Chemical Properties Handbook. McGraw-Hill; 1999.
66. Samudrala R, Levitt M. Protein Sci. 2000;9:1399–1401. [PubMed]
67. Simons KT, Ruczinski I, Kooperberg C, Fox BA, Bystroff C, Baker D. Proteins: Struct. Func. Genet. 1999;34:82–95. [PubMed]
68. Keasar C, Levitt M. J. Mol. Biol. 2003;329:159–174. [PMC free article] [PubMed]
69. Xia Y, Huang ES, Levitt M, Samudrala R. J. Mol. Biol. 2000;300:171–185. [PubMed]
70. Melo F, Sanchez R, Sali A. Protein Science. 2002;11:430–448. [PubMed]
71. Zhang C, Liu S, Zhou HY, Zhou YQ. Protein Sci. 2004;13:400–411. [PubMed]
72. Zhang C, Liu S, Zhou HY, Zhou YQ. Biophys. J. 2004;86:3349–3358. [PubMed]
73. Zhou HY, Zhou YQ. Protein Sci. 2002;11:2714–2726. [PubMed]
74. Zhou HY, Zhou YQ. Protein Sci. 2003;12:2121–2121.
75. Melo F, Feytmans E. J. Mol. Biol. 1997;267:207–222. [PubMed]
76. Cascella M, Neri MA, Carloni P, Dal Peraro M. J. Chem. Theory Comput. 2008;4:1378–1385.