Comparison of multiple Amber force fields and development of improved protein backbone ...

14 Pages·2006·0.68 MB·English

PROTEINS:Structure,Function,andBioinformatics65:712–725(2006)
Comparison of Multiple Amber Force Fields and
Development of Improved Protein Backbone Parameters
ViktorHornak,1RobertAbel,2AsimOkur,3BentleyStrockbine,4AdrianRoitberg,2andCarlosSimmerling1,3,5*
1CenterforStructuralBiology,StonyBrookUniversity,StonyBrook,NewYork11794
2DepartmentofChemistryandQuantumTheoryProject,UniversityofFlorida,Gainesville,Florida32611
3DepartmentofChemistry,StonyBrookUniversity,StonyBrook,NewYork11794
4DepartmentofPharmacologicalSciences,StonyBrookUniversity,StonyBrook,NewYork11794
5ComputationalScienceCenter,BrookhavenNationalLaboratory,Upton,NewYork11973,
ABSTRACT The ff94 force field that is com- chanical (QM) wave function. Unfortunately, because of
monly associated with the Amber simulation pack- thelargesizesofbiologicalmacromolecules,suchcalcula-
age is one of the most widely used parameter sets tions are extremely time consuming and can only be
forbiomolecularsimulation.Afteradecadeofexten- applied to a limited number of conformations. Therefore,
sive use and testing, limitations in this force field, mostsimulationprojectsemployclassicalmolecularmechan-
such as over-stabilization of a-helices, were reported ics(MM) energy functions. Because of the need to evalu-
byusandotherresearchers.Thisledtoanumberof ate these functions a large number of times during the
attemptstoimprovetheseparameters,resultingina simulation, these functions are relatively simple and uti-
variety of ‘‘Amber’’ force fields and significant diffi- lize many adjustable empirical parameters. These are
cultyindeterminingwhichshouldbeusedforapar- most often obtained by fitting to data from experiments
ticular application. We show that several of these or from high level QM calculations. The energy function
continue to suffer from inadequate balance between togetherwiththesetofempiricalparametersisknownas
different secondary structure elements. In addition, aforcefield.
theapproachusedinmostofthesestudiesneglected Inspiteofpersistent yetslowemergenceofforcefields
to account for the existence in Amber of two sets of that explicitly account for charge polarization,1 the fixed
backbone u/w dihedral terms. This led to parameter charge additive force fields remain popular because of
sets that provide unreasonable conformational pref- computationalefficiency.Foramorein-depthoverviewof
erences for glycine. We report here an effort to im- currentforcefieldsandtrendsinforcefielddevelopment,
provetheu/w dihedraltermsintheff99energyfunc- the reader is referred to recent review articles.2–4 In this
tion. Dihedral term parameters are based on fitting report,wefocusontheexistingAmberfixedchargeaddi-
the energies of multiple conformations of glycine tiveforcefields.
and alanine tetrapeptides from high level ab initio The ‘‘Cornell et al.’’ force field5 (denoted as ff94 in
quantum mechanical calculations. The new parame- Amber) has been the most widely used with the Amber
ters for backbone dihedrals replace those in the suite of programs6,7 since its publication over a decade
existing ff99 force field. This parameter set, which ago.Itintroducedthesetofparameters forall-atomsim-
we denote ff99SB, achieves a better balance of sec- ulationssuitableforproteinsimulationsinthecondensed
ondary structure elements as judged by improved phase, largely inspiredby theoptimized potentialforliq-
distribution of backbone dihedrals for glycine and uid simulations (OPLS). Some characteristic features of
alanine with respect to PDB survey data. It also ff94includefixedpartialchargesonatomcenters,explicit
accomplishes improved agreement with published use of all hydrogen atoms, no specific functional form for
experimental data for conformational preferences of hydrogenbonding,anddihedralparametersfittorelative
shortalaninepeptidesandbetteraccordwithexper- QMenergiesofalternaterotamersofsmallmolecules.In
imentalNMRrelaxationdataoftestproteinsystems. particular, the protein u/w dihedrals have specific rota-
Proteins2006;65:712–725. VVC 2006Wiley-Liss,Inc.
Keywords: trialanine; dihedral parameters; mole- TheSupplementaryMaterialreferredtointhisarticlecanbefound
cular dynamics; molecular mechanics; athttp://www.interscience.wiley.com/jpages/0887-3585/suppmat/
decoy analysis; NMR order parame- Grant sponsor: National Institutes of Health; Grant number:
ters; a-helix GM6167803;Grantsponsor:DepartmentofEnergy;Grantnumber:
DE-AC02-98CH10886.
*Correspondenceto:CarlosSimmerling,DepartmentofChemistry,
INTRODUCTION StonyBrookUniversity,StonyBrook,NY11794.
E-mail:[email protected]
Molecular modeling studies of biologically important
Received 23 February 2006; Revised 4 May 2006; Accepted 15
molecules usually require the evaluation of potential en-
June2006
ergies for alternate conformations. The most accurate
Publishedonline15September2006inWileyInterScience(www.
approach would encompass the use of the quantum me- interscience.wiley.com).DOI:10.1002/prot.21123
VVC 2006WILEY-LISS,INC.
IMPROVEDAMBERPROTEINBACKBONEPARAMETERS 713
tional parameters that affect relative energies of alter- haveimprovedthatspecificweaknessbutamoresystem-
nate backbone conformations. These were fit to optimize atic revision was clearly necessary to quantify the prob-
agreement with QM relative energies for several confor- lem and address the issue with generality and transfer-
mations of glycine and alanine. The partial atomic abilityinmind.
charges were derived by fitting the (gas phase) electro- Much of the confusion and suboptimal performance of
static potential calculated at the Hartree-Fock 6-31G* recent Amber protein force field variants may have been
level. This approach intentionally ‘‘overpolarizes’’ bond caused by a somewhat nonintuitive parameterization of
dipolesas present in the gas phase, such that the result- protein backbone dihedral angles. In Amber, each dihe-
ing charge distribution approximates those occurring in dral profile is defined by a set of four atoms. The set of
aqueouscondensedphase.9 atoms used to define u and w for glycine is as expected,
Becauseoflimitedcomputationalresourcesatthetime, following u and w along the main chain (u ¼
dihedral parameters were fit to a small number of low- C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)C,w¼N(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)N).Importantly,forother
energyconformationsofglycineandalaninedipeptides.A amino acids that have a side chain, an additional set of
possible limitation of using dipeptides is that their gas dihedrals also influences rotation about the u/w bonds
phase energy surfaces do not have a local minimum in connectingtheCaatomtotheamideCandNatoms.This
thea-helicalregion,whichoccurswithhighfrequencyin extra set of terms corresponds to dihedral angles
branched out to the Cb carbon, which in this work we
protein structures. This was addressed in subsequent
modifications of ff94, such as ff9610 and, more recently, identifyasu0¼C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)Cbandw0¼Cb(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)N.
ff99.11Inff96,identicalparameterswereusedforuandw Thisdefinitionofdihedralswasoriginallyimplemented
in ff94, with u and w that fit independently to glycine
terms, which were empirically adjusted to reproduce the
data and u0 and w0 used to adjust the behavior for
energydifferencebetweenextendedandconstraineda-hel-
alanine. Later modifications of ff94, such as ff96, ff99,
icalenergiesforalaninetetrapeptide.ff99reflectedanother
Garcia’s13, Pande’s19 and ours12,18 only changed the first
attempt at refitting backbone dihedral parameters by
setofu/wterms,usingthemtoadjustbackbonepreferen-
including11representativestructuresofalaninetetrapep-
ces for alanine. Thus the new u and w parameters were
tidealongwiththealaninedipeptides.Duringthisprocess,
fitinthepresenceoftheff94alanine-basedu0 andw0 and
however, a new problem was introduced in the way ff96/
will only give their intended behavior, when they accom-
ff99u/w dihedralparameterswereoptimized(discussedin
panytheseu0 andw0 dihedrals.Importantly,u0 andw0 are
detaillater),resultinginincorrectconformationalpreferen-
not present in glycine. For example, ff99 modified the
cesforglycine.Inadditiontochangesinbackbonedihedral
backbone u and w terms in order to reproduce relative
parameters,ff99alsogeneralizedandextendedatomtypes
energies for alanine. However, the alanine residue also
tocompoundsbeyondaminoandnucleicacids. hadtheu0 andw0 dihedralparameterspresentduringthe
Eventhoughbothff99and,particularly,ff94parameter
u/w modification. These optimizations may perform well
setshavebeensuccessfullyusedformanyyears,improve- fornonglycineresiduesthatpossesstheu0/w0terms.How-
ment in conformational sampling due to increased com-
ever, when these modified backbone dihedral parameters
puterpowerandalgorithmicadvancesrevealedthatboth
areappliedtoglycine,theresulthaslittlephysicaljustifi-
force fields over-stabilize a-helical peptide conforma- cation, since they were fit in the presence of u0 and w0. A
tions.12,13 On the other hand, the changes introduced in
similar inconsistency is present in Garcia’s modification;
the less frequently used ff96 were observed to overesti- only the u and w terms were zeroed; thus, while glycine
mate b-strand propensity.14–17 Because the backbone di-
has no backbone dihedral potential all other amino acids
hedralparametersaresharedbyallaminoacids,regard- retaintheff94u0andw0terms,whichwerefitinthepres-
less of the type of side-chain, their cumulative effect is ence of the glycine-based u and w terms, and thus have
most likely responsible for the bias towards the specific nomeaningwhenappliedwithoutthem.
secondary structure. This recently motivated us12,18 and We show below that these modifications result in
others13,19 to revisit backbone dihedral parameterization unreasonablesamplingofdihedralspaceforseveralpost-
for ff94/ff99. Garcia and Sanbonmatsu13 simply zeroed ff94Amberforcefields.Theseproblemsareinadditionto
thetorsionpotentialforuandw(althoughasweexplain the overstabilization of helical conformations present in
later, only some of the u/w terms were removed) and ff94 and ff99. Amber and other popular force fields (e.g.
noted improved agreement with experimental helix-coil Charmm22 and Gromos96 variants and OPLS-AA) were
parameters.Morerecently,SorinandPande19,20modified previouslycompared21,22intermsofconformationalsam-
ff99byreplacingitsudihedralparameterswiththeones pling of blocked glycine and alanine dipeptide and were
fromff94.This improved the agreementwith experimen- noted todisagree tovariousdegrees(particularly forgly-
talkineticand thermodynamicmeasurementsforfolding cine21) with simulations employing combined QM/MM
of two 21-residue a-helical peptides. Our own modifica- forcefieldaswellaswithstatisticalanalysisofhighreso-
tionofff99,publishedpreviously,18specificallyaddressed lution protein crystal structures. In another study, all
stronghelicalbiasand,inthisrespect,succeededtoachieve force fields were also shown to perform very differently
betterbalanceofmajorsecondarystructuretypes. with respect to a number of properties of trialanine com-
Allthesepreviousmodificationswereaimedatcorrect- paredwithNMRandinfraredobservables.23
inganapparentproblemencounteredwithaspecificsys- We also show that our reparameterization of backbone
tem.Theheuristicapproachesthatwereemployedmight dihedral parameters improves conformational preferen-
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
714 V.HORNAKETAL.
Fig.1. Blockedalaninetetrapeptide(threealanineresiduesbutfour peptide bonds)usedintheoptimi-
zation procedure. Note the definition of dihedral angles: u ¼ C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)C, w ¼ N(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)N, u0 ¼
C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)Cb,w0 ¼Cb(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)N.Therearethreeu/wpairsandthreeadditionalu0/w0 pairsinthistetra-
peptide. Glycine tetrapeptide would only have the u/w dihedrals due to the absence of Cb carbon. As an
example,w andu0 areshowninbold.
1 3
ces for typical secondary structures. It should be noted Eventhoughitmayatfirstseemmorestraightforwardto
thatthisworkdoesnotattempttocreateyetanothervar- have a single set of backbone u/w parameters for all
iantofAmberforcefieldbutstrivestoimprovetheexist- aminoacids,thedifferentnatureofglycine(noCbcarbon)
ingff94andff99forcefields,aswellasspecificallyaddress motivatedthefollowingapproach:thereisindeedasingle
the problem with glycine sampling. Thus rather than con- set of backbone u/w parameters for glycine (defined as
tributing to the diverging set of Amber force fields, we u ¼ C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)C, w ¼ N(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)N), but any other
extendtheevolutionstartedwithff94,followedbyff99and aminoacidhasanadditionalsetofu/wparameterswhich
subsequentlyrefinedinthisreport.Thisisalsoreflectedby are added to the glycine ones (here defined as u0 ¼
naming this modified force field ff99SB. We attempted to C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)Cb, w0 ¼ Cb(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)N). In other words, for
maximallyimprovethemethodologythatweusedtoderive any nonglycine amino acid, the dihedral energy is the
ournewparameterswhilestillfollowingthesamegeneral sum of dihedral energies calculated for u/w and u0/w0.
philosophyaswasusedwithff94/ff99. Notehowever,thatu0/w0 iscalculatedforadihedralangle
A different approach was taken by Duan et al.,24 who shiftedby(cid:2)1208becausethatistheoffsetofthetwotor-
recentlyintroducedamoreextensivemodificationofff94/
sions(asfollowsfromtheirdefinitionsusingdifferentsets
ff99(calledff03),inwhichafundamentallydifferentcon-
of four atoms). Figure 1 demonstrates the definition of
cept to derivation of partial atomic charges was used.
thesetwosetsofdihedralanglesonAla tetrapeptide.
Instead of relying on the HF/6-31G* approach to provide The glycine u/w parameters were opt3imized first, usu-
aqueous-phasecharges,alow-dielectriccontinuummodel
allybasedonthebestreasonablefitofQMandMMener-
corresponding to an organic solvent environment was
giesforasetofglycineconformers.Thesearetotalpoten-
included directly in the QM calculation of the dihedral
tial energies for the molecule, where backbone dihedral
parameters and electrostatic potential (from which the
parameters are the adjustable variables in the optimiza-
charges are obtained). Because of these differences, ff03
tion.Forexample,ifweoptimizedtheabsolutedifference
should be considered a distinct force field model rather
of QM and MM energies for glycine tetrapeptide (Gly ),
thanextensionofpreviousAmberforcefields. 3
we would minimize the absolute error (AE) between the
Amberforcefieldsarenottheonlyonesundergoingsuch
MMandQMenergies[Eq.(1)]:
evolutionarychangesbroughtaboutbydiscoveryoflimita-
tions through improved conformational sampling capabil- XN (cid:1) (cid:1)
ities. Similar efforts have recently been reported by Mac- AE¼1=N (cid:1)E ðiÞ(cid:1)E ðiÞ(cid:1) ð1Þ
QM MM
Kerell et al.,25 improving the backbone parameters in the i¼1
Charmm force field, specifically removing the bias toward
p-helical peptide conformations.26,27 In this case, the where EQM(i) and EMM(i) correspond to the QM and MM
energies respectively for i-th glycine tetrapeptide con-
Charmm22 alanine dipeptide dihedral energy surface was
former,andNisthenumberofallglycineconformers.In
corrected using a grid-based difference map (CMAP) to
our case, MM energy for a given conformer is given by
achieve an almost perfect match with LMP2/cc-pVQZ(-g)
theAmberenergyfunction(seeff945):
energy surface. Likewise, dihedral terms in the original
OprPoLveS-tAhAeafgorreceemfeienltdw28itwhehrieghallseoverleQpMaradmateat.e2r9ized to im- EMM¼EbondþEangleþEnon-bondþEdihedral ð2Þ
wherethedihedralenergytermis
METHODS X3
Optimization ofBackbone Dihedral Parameters E ðuÞ¼ V ð1þcosðnu(cid:1)g ÞÞ ð3Þ
dihedral n n
Tobetterunderstandtheprocedureweusedforobtain- n¼1
ing a new set of backbone dihedral parameters, we first V is dihedral force constant (amplitude), n is dihedral
n
outline the general procedure originally used in ff94.5 periodicity, and g is a phase of the dihedral angle y
n
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
IMPROVEDAMBERPROTEINBACKBONEPARAMETERS 715
(which would be either u or w for backbone dihedral dralparameters,andnochangewasmadetotheff94par-
terms). The Fourier series in E is approximated tialatomiccharges.
dihedral
using a small number of terms. Our choice of three The optimization of dihedral parameters was done in
terms in the expansion is consistent with the AMBER twosteps.First,theu/wtorsionparametersforGly were
3
philosophy of using dihedral terms that can be physi- optimized as follows. We used an exhaustive grid search
cally rationalized30 and yields a reasonable number of in12-dimensionaldihedralparameterspace:threeampli-
terms persingle dihedral (six)thatneed to beoptimized. tudes (V ,V ,V ) and three phases(g ,g ,g ) foreachu
1 2 3 1 2 3
DihedralenergyparametersinEq.(3)correspondingto and w. Amplitudes were systematically varied between
backboneu/wdihedrals(i.e.,V andg )aretheadjustable 0 and 2 kcal/mol with a step of 0.2, and phases were set
n n
parameters. We seek values of these parameters that to either 0 or p radians. Thus, the total number of
minimizethefunctioninEq.(1).Oncetheglycineparam- grid points in the initial search was 116 3 26 (over 100
etersareoptimized,theyareheldfixedandusedinasec- million).
ond round of fitting where u0/w0 parameters are opti- Theinitialgridspacingwas0.2kcal/molforamplitudes
mized. This is carried out by fitting the parameters to (V ),whichwouldbethemaximumprecisionofoptimized
n
best reproduce QM and MM energies for a set of alanine parameters. To improve this, we employed ‘grid focusing’
tetrapeptide (Ala3) conformations. This second round of by repeating the grid search with finer grid spacing cen-
u0/w0 optimization was initially misunderstood by us and tered on the parameters resulting from the initial course
byotherswhomodifiedtheoriginalff94.Theuseofasin- grid. Our refined grid had a spacing of 0.01 kcal/mol for
gleroundfittinginwhichtheglycineu/wdihedralparam- amplitudes,andthus,thisistheultimateprecisionofour
eters were modified to improve agreement for alanine optimizedgridparameters.
data, where the ‘real’ u0/w0 alanine parameters were The function we used for optimization is somewhat
untouched, led to incorrect parameterization (primarily more complicated than that shown in Eq. (1). Since the
forglycine)andoverallconfusionabouttheeffectofdihe- zerooftheMMenergyfunctionisarbitrary,theoptimiza-
dralparametersonthebackboneconformation.
tion should be performed using energy differences be-
In the following, we describe the specifics of this gen-
tweenalternateconformations.Thiswasachievedinff99
eralprocedure thatwere employed forthe present study.
bysettingoneoftheenergies(thelowestone)tozero,and
All optimizations were carried out with blocked (acetyl
all other energies were assigned values relative to this
and N-methylamine groups at the N- and C-termini,
‘‘zero energy’’ reference. However, we observed that the
respectively) glycine and alanine tetrapeptides (Gly ,
3 optimized parameter set changes depending on which
Ala ) as shown in Figure 1. A set of Gly and Ala tetra-
3 3 3 conformer’s energy is used as a reference. Therefore, the
peptideconformerswaschosentorepresentlocalminima
function we optimized was calculated as an average of
on the ff94 energy surface (ff99 was not used here
QM and MM energy differences with each conformer’s
because of apparent flaws in the ff99 energy function
energy set as a reference in turn. This gives an average
for glycine). Because the six-dimensional dihedral space
absoluteerror(AAE)definedasfollows:
of tetrapeptides is too large to be represented exhaus-
tively (as would be possible for simple two-dimensional XN XN (cid:1) (cid:1)
surface of dipeptides), we chose to optimize the relative AAE¼1=NðN(cid:1)1Þ (cid:1)Eo;i ðjÞ(cid:1)Eo;i ðjÞ(cid:1) ð4Þ
QM MM
depths of local minima on the energy surface, because i¼1 j>i
these are most relevant to the thermodynamic stability
ofalternatesecondarystructuretypes.Eachpointinthis where Eo,i (j) is the QM energy of conformer j with con-
QM
dihedralspace corresponds toaspecific tetrapeptide con- former i as a reference, and Eo,i (j) is the MM energy of
MM
formation. Local minima were identified through a sto- conformerjwithconformeriasareference,Nisnumber
chastic search in u/w dihedral space followed by local ofconformers,whichis28forGly .Anotherpossibilityto
3
minimization.
define a function for optimization is to concentrate on
All conformations corresponding to the proposed min-
whatthemaximumabsoluteerror(MAE)mightbewhen
ima were optimized at the Hartree-Fock level with 6-
we consider all individual QM and MM differences and,
31G*basisset.After discardingtheconformers thatcon-
again,takeintoaccountthatanyconformermayserveas
vergedtothesamegeometry,weendedupwith28glycine
a reference ‘‘zero energy.’’ This would give a function of
and51alaninetetrapeptideconformers(withamaximum
thefollowingform:
range ff94 energy cutoff of 10 kcal/mol). The energies of
(cid:1) (cid:1)
these conformations were calculated using higher level MAE¼max(cid:1)Eo;i ðjÞ(cid:1)Eo;i ðjÞ(cid:1) ð5Þ
(gasphase)QMLMP2/cc-pVTZ(-f)inJaguar,31employing i;j>i QM MM
the methodology described by Beachy et al.,32 who
showedthatinclusionofelectroncorrelation(inthiscase Here, we would search for such set of dihedral parame-
using LMP2) is important and has a significant effect on ters that give the smallest MAE. As it turns out, there
relative energies of peptide conformations. Tables of all are many sets of dihedral parameters with very similar
Gly /Ala conformersandtheirQMenergiesareprovided AAE.However,inspectingMAErevealsthatthisparame-
3 3
in the supplementary material (Tables S1 and S2). Note ter is much more sensitive to different parameter sets
that we only used these energies to refit backbone dihe- and therefore the selection of the final set was based on
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
716 V.HORNAKETAL.
Fig.2. Ramachandranplotof28Gly (left)and51Ala (right)conformersasobtainedfromQMgeome-
3 3
tryoptimization of structures that were acquired using stochastic search with MM energies.All threedihe-
dralpairsareshownforeachconformer(blacksquares).Thegraypointsinthebackgroundareu/wvalues
ofglycineandalanineresiduescollectedfromasubsetofthePDB.Typicalsecondarystructureregionsare
outlinedbyboxesintheAla plot(seeMethodsforadditionaldetails).Ala conformernumber16(outlier)in
3 3
thelowerleftofAla plotisdesignatedbyanarrow.
3
this error estimate. In general, small values of MAE in the main text: ff94,5 ff99,11 Garcia’s modification of
alwayscorrelatewithsmallvaluesofAAE,suchthatthe ff9413withC(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)CandN(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)Ntermszeroed
best parameter set picked based on MAE will also have (denoted ff94gs in the text), Pande’s modified ff9919 with
oneofthelowestAAEvalues. C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)C term replaced by the one from ff94
For Gly optimization, an additional constraint was (denotedff99/),ff0324(aspresentinAmber8distribution)
3
introduced,requiringthefirstphaseforu(i.e.g foru)to and ff99SB developed as described earlier. The time step
1
be 0. This was necessary to obtain physically reasonable was 2 fs, and all bonds involving hydrogen were con-
dihedral functions with energy barrier located at u (cid:3) 08 strained by SHAKE with a tolerance of 10(cid:1)4-A˚. Glycine
on the glycine Ramachandran map. Even though glycine and alanine tetrapeptide systems were solvated by (cid:2)520
parameter sets without this constraint achieved slightly TIP3P33 water molecules in a periodic box. Simulations
better fit of QM and MM data, the resulting dihedral were carried out in the NPTensemble at 300 K. A cutoff
energy function would have a physically unreasonable of 8-A˚ was used for nonbonded interactions, and long-
minimum in the u (cid:3) 08 region of Ramachandran map range electrostatic interactions were treated with the
(where we also did not have any QM data). This arises particlemeshEwald(PME)method.34
from our use of only local minima to train the dihedral Replica exchange molecular dynamics (REMD) simula-
correctionterms. tions,35 as implemented in Amber8, were run for Ala
3
Once the Gly parameter set was obtained, the same using explicit water solvation and the ff99SB and ff03
3
procedure was repeated with 51 Ala conformers, and forcefields.REMDsimulationsweresetupundersimilar
3
parameter sets corresponding to u0/w0 dihedrals were ob- conditions as the corresponding standard MD runs. Ala
3
tained.Althoughtheresultingsetofparametersappeared was solvated in a truncated octahedral box using 595
satisfactory, a more detailed analysis of errors revealed TIP3Pwatermolecules.Twentysixreplicaswereusedto
that Ala conformer number 16 was repeatedly responsi- spanthe temperaturerangeof267–571K.Temperatures
3
ble for the largest fitting errors (as measured by AAE or wereoptimizedtogiveauniformexchangeacceptancera-
MAE). This conformer also falls into the unusual region tio of (cid:2)30%. Exchange between neighboring tempera-
oftheRamachandranmap(Fig.2right),andtherefore,it tures was attempted every picosecond and each REMD
was excluded from the optimization of Ala dihedral simulation was run for (cid:2)30,000 exchange attempts (30
3
parameters. ns). Because of the high temperatures of some replicas,
The resulting parameters in the form suitable for use theNVTensemblewasemployed.
with Amber are provided in supplementary material
(TableS3).
Calculation of NMROrder Parameters
Heneggwhitelysozyme(PDBcode6LYT36)andubiqui-
Simulationsin ExplicitWater
tin(PDBcode1UBQ37,38)simulationswererununderthe
All molecular dynamics simulations for Gly and Ala same conditions as the short peptides described earlier,
3 3
tetrapeptideswerecarriedoutwiththesandermodulein but with truncated octahedral periodic boundary condi-
Amber87 using several different force fields as discussed tions including 4350 and 3300 water molecules, respec-
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
IMPROVEDAMBERPROTEINBACKBONEPARAMETERS 717
tively. Equilibration was done in multiple steps, with the secondary structure basins were normalized by the
positional restraints on all heavy atoms in the crystal totalnumberofstructuresinjustthefourregions.
structuregraduallyreleasedfrom2to0kcal/molA˚.2The The convergence of Ala simulations was estimated by
3
productionruntrajectorieswere30nslong.Thesetrajec- comparing the first and second half of the trajectories in
tories were used to calculate the N(cid:1)(cid:1)H internuclear vec- termsofdifferencesofpopulationfractions.Thefractions
tor autocorrelation functions with ptraj module. Order were calculated independently for the two parts of the
parameters (S2) were obtained from a plateau region of trajectoryandtheresultingvaluesaregivenasaverages,
autocorrelation functions.39 Specifically, the autocorrela- withdeviationsreflectingthedifferencesofthetwoinde-
tion function was calculated up to time of half of trajec- pendentcalculations.
tory length and the mean of last 5 ns was taken as S2.
Standarddeviationofthemeanwasusedtoplottheerror
Decoy Analysis
bars. Since autocorrelation functions for a few backbone
amidestypicallydonotconvergewell(meaningthatthey Similar to the work that we previously published12 for
do not show clear plateau), the mean will be a less reli- ff94andff99,weused‘‘decoyanalysis’’totestthecapabil-
able estimate of S2 for these amides, which will be ity of the force field to identify the experimentally deter-
reflected by larger standard deviations (i.e., larger error mined native structure as that with the lowest potential
bars). energy. Three systems were used for decoy analysis:
trpzip242 (SWTWENGKWTWK-NH ), Baldwin type43 al-
2
anine based a-helical peptide44 (Ace-GGG(KAAAA) K-
3
NH ), and trpcage miniprotein45 (NLYIQWLKDGGP-
Dihedral PlotsandFree Energy Surfaces 2
SSGRPPPS). The decoys for the three systems were gen-
A Ramachandran plot showing the distribution of gly- eratedbyanumberofunrestrainedaswellasforcedMD
cine or alanine u/w angles in protein crystal structures simulations at 300 K as described previously.12 The sets
wasconstructedbasedonasubsetofproteinsfromaPDB cover structures generated in simulations started from
database, with high resolution (<1.6 A˚), small R-factor various initial structures, such as native, extended, and
(<0.25) and less than 20% sequence homology.40,41 Dihe- multiple unfolded or partially folded structures from
dralplotsbasedonthisstatisticalsurveyofthePDBpro- intermediary stages of folding trajectories. The decoy
vide experimental indication of the values of backbone ensembles thus contain local and global structural varia-
torsionsthatarecommonlyfoundinproteins. tion with a greater diversity of structures than could be
Gly3/Ala3simulationswereusedtoproducefreeenergy obtained from a single simulation. The sizes of the three
surfaces(orpotentialofmeanforce,PMF,maps)forback- decoy sets are as follows: 177,000 for trpzip, 270,000 for
bone dihedral angles, assuming that u/ware represented thehelicalpeptide,and118,000fortrpcage.Theenergies
in the ensembles according to a Boltzmann distribution. of these decoy sets were analyzed with different force
First,thevaluesofallu/wdihedralpairsovertheproduc- fields using a Generalized Born solvation model46 imple-
tion portions of simulations were collected. Then, two- mented in Amber. For each force field, the potential
dimensionalu/wnormalizedhistogramswereconstructed energyandrootmeansquaredifference(RMSD)fromthe
from these values and converted to free energies using experimentally determined structure were plotted for
Eq.(6): each of the decoy structures. The expectation for a
‘‘correct’’ force field is that the lowest energies are ob-
DG ¼(cid:1)RTlnðN=N Þ ð6Þ
i i 0 tained forstructures with low RMSDs, thatis, fornative
folds.
where N is the population of a particular histogram bin Because many different structures can have similar
i
for specified values of u/w, and N is the most populated RMSD values, and because different force fields may
0
bin. Thus, the global free energy minimum always has a favor slightly different conformations in each RMSD
valueof0kcal/mol. range, we performed a further simplification of the data.
Theu/whistogramswerealsousedtoevaluatepopula- Since we are primarily interested in structures with low
tions in different regions of secondary structures. The energies, we only plot one energy value, which is the av-
definitions of the four principal regions were as follows: erage of the 20 lowest energies in that particular RMSD
right-handed a-helix (a ), (u, w) (cid:3) ((cid:1)708 (cid:4) 308, (cid:1)458 (cid:4) bin.Thisresultsinacurvethatoutlinesthelowestener-
R
458);lefthandeda-helix(a ),þ60(cid:4)308/þ45(cid:4)458;poly- gies sampled for the RMSD range (‘‘lowest energy pro-
L
proline II (PP ), (cid:1)708 (cid:4) 308, þ 1508 (cid:4) 308; extended b- file’’), which can be more easily interpreted as being rep-
II
strandconformation (b),(cid:1)1508(cid:4)308,þ1508 (cid:4)308.Most resentativeoftheunderlyingenergysurfacefortheforce
ranges are (cid:4)308 except from the w range of a and a , field. This approach has several limitations, including
R L
which were slightly wider ((cid:4)458) to capture the shift in neglectofentropiccontributionstotheactualfreeenergy
theseregionsthatresultedfromusingtheff94gsparame- andalsousesanimplicitsolvationmodelwhichmayhave
ter set. The number of structures in individual regions itsowneffectsonthestabilityofdifferentconformations.
were summed and divided by the total number of struc- Nevertheless, this methodology is useful for fast qualita-
tures to obtain population fractions. Relative fractions tivescreeningofmanyvariationsofaforcefieldonlarger
are also reported; these differ in that the populations in peptides.
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
718 V.HORNAKETAL.
RESULTSAND DISCUSSION TABLEI. OptimizedDihedralParameters:
Amplitudes(V)andPhases(g)forBackbone
Theuniqueconformationsof28Gly and51Ala thatre-
3 3 DihedralAngles
sultedfromQMoptimizationareshowninFigure2.Ingen-
eral,theseconformersoverlaptheregionsoutlinedbyPDB V1 g1 V2 g2 V3 g3
survey data. There are several distinct clusters that fall / C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)C 0.00 0 0.27 0 0.42 0
withintheregionsoftypicalproteinsecondarystructures, w N(cid:1)(cid:1)Ca(cid:1)(cid:1)C(cid:1)(cid:1)N 0.45 p 1.58 p 0.55 p
such as right handed (a ) and left handed (a ) a-helix, /0 C(cid:1)(cid:1)N(cid:1)(cid:1)Ca(cid:1)(cid:1)Cb 2.00 0 2.00 0 0.40 0
R L
polyprolineII(PP ),andextendedb-strandconformations. w0 N(cid:1)(cid:1)C(cid:1)(cid:1)Ca(cid:1)(cid:1)Cb 0.20 0 0.20 0 0.40 0
II
FollowingtheproceduredescribedinMethods,weopti- Theu/wtermsarealwaysusedbyallresidues,u0/C0 areusedbyall
mized a set of dihedral parameters (Table I). The RMSD residues except glycine. Phases are in radians and amplitudes are
between the QM and MM energies for glycine tetrapep- inkcal/mol.
tideusingtheseparametersis1.17kcal/mol.Theaverage
absolute error and maximum error are 0.93 and 2.71
Onemightexpectthattheu/wregionssampledincrystal
kcal/mol, respectively. Alanine tetrapeptide calculations
structures should be also accessible to solvated Gly /Ala
used the glycine u/w backbone parameters as well as the 3 3
u0/w0 dihedralterms.Inthiscase,theRMSD,absoluteav- peptides. However,itwouldbewrong toassume thatthe
relative populations of u/w angles in static crystal struc-
erage error and maximum error are 1.31, 1.05, and 3.43
tures should match populations in small tetrapeptides
kcal/mol. This is comparable to similar fitting errors
solvated in water. Not only are the two environments
obtained in ff9911 and ff03.24 For comparison, the values
very different but also populations in crystal structures
ofbackbonedihedralparametersformanyothervariants
will certainly bebiased by occurrence of common second-
ofAmberforcefieldaregiveninthesupplementaryinfor-
ary structures, such as a-helices and b-sheets. Further-
mation (Table S4). It should be noted again that RMSD,
more, there is no obvious choice for temperature associ-
absolute average error (AAE) and maximum absolute
ated with the PDB distribution, which would make any
error (MAE) are calculated such that each conformer in
freeenergycalculationproblematic.
turnissetasa‘‘zeroenergy’’reference,andthusallpairs
We first discuss the plots for Gly . The plots obtained
of energy differences are used in the average. The sensi- 3
from the Gly simulations are symmetric with respect to
tivityoferrors(RMSD,AAE,MAE)tothechoiceofrefer- 3
theorigin,whilealanineplotsarenotduetothepresence
ence can be demonstrated by calculating these errors for
of the chiral center on the alanine a-carbon. A certain
all 51 Ala conformers with a reference conformation
3 amount of asymmetry in the glycine PDB data arises
fixedtoeachconformerinturn.Suchevaluationoferrors
from the influence of chiral centers in other residues,
produces best RMSD of 0.94 if conformer #31 is a refer-
whicharenotpresentintheGly simulations.
enceand worstRMSD of1.97for conformer #45as aref- 3
The analysis of glycine conformations in the PDB data
erence. Similarly, depending on the choice of conformer
reveals thatbothright-and left-handeda-helical regions
as a reference, AAE ranges between 0.77 and 1.82 and
are densely populated. PP and extended b-strand
MAE between 1.74 and 3.43 kcal/mol. Thus, averaging II
regions are also found with high frequency. One should
overallpair-wisedifferencesshouldprovideamorestrin-
expect similar conformational preferences resulting from
gent error estimation (e.g., MAE is always the largest of
MD simulations (Fig. 3). This is largely true for ff99SB
all pair wise differences when using all pairs) and yields
and ff03 and to a lesser extent for ff94 (which is missing
an independence of results with respect to which con-
PP andb,andhasstrongbiasfavoringa).However,ff99
formerischosenasareference. II
andff94gsbothshowunexpectedpatternsforglycinecon-
formations, sampling regions which are not represented
ForceField Validation
inPDBsurveydata.Aswenotedbefore,theselikelyarise
Explicit watersimulations ofGly /Ala
3 3 from lack of treatment of both sets of backbone parame-
We performed molecular dynamics simulations on the ters (u/w and u0/w0). The ff99 glycine parameters were fit
same systems we used for parameter optimization, that usingalaninedataandfavorconformationsthatshowno
is, Gly and Ala peptides, with different Amber force population inthePDBdata (u¼(cid:1)1808and w¼08).The
3 3
fields. All simulations were fully solvated with explicit ff94gs surface is also in significant disagreement with
water and extended to (cid:2)80 ns. Histogram analysis was PDB data, with a surface that has minima and low bar-
used to calculate relative free energies of different con- rier heights that appear highly dissimilar to the PDB
formers with respect to u/w backbone angles. These free dataandtoanyoftheforcefieldsthatdoincludeuandw
energy (PMF) surfaces are shown in Figure 3 for both terms. This may artificially increase the glycine confor-
Gly andAla . mationaltransitionrates.Theplotsforff96,ff99/aswell
3 3
The correspondence of the distributions obtained from astheplotsobtainedfromsimulationsthatusedourpre-
PDB and simulations has been used extensively in force vious modification of dihedral parameters12,18 are shown
fieldvalidationanddevelopment.21,24–26Thisisthebasis in Figure S6 (supplementary material) and, similarly to
forour first testand qualitatively evaluatesour dihedral ff99 and ff94gs, exhibit incorrect sampling for glycine
parameter set. PDB survey data are presented here as dihedralsduetothelackofexplicittreatmentofbothsets
simpleRamachandranmaps,notasfreeenergysurfaces. ofbackboneparameters.
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
IMPROVEDAMBERPROTEINBACKBONEPARAMETERS 719
Alanine free energy maps for the different force fields rapid thermalization of these basins during the simula-
(Fig. 3) show more similarity to each other, with major tions.However,thea regionisseparatedfromtheother
L
regions represented by right handed a-helix, PP and basins by a higher barrier near u ¼ 0, resulting in much
II
extended b-conformations. Once again, differences in the less frequent transitions from other regions (see Figure
relative free energies of the basins are readily apparent S5inthesupplementary informationforfurtherdetails).
in the different force fields. Although ff99SB and ff03 To ensure that the relative populations in that region
show a reasonable balance between a and PP , ff94 were determined reliably, we ran REMD for ff99SB and
R II
and ff99 show clear bias in favor of a . ff99 also shows ff03,whichwerethetwoforcefieldsthatprovidedconfor-
R
the same tendency to adopt too-low (close to (cid:1)1808) mational ensembles most consistent with the PDB data.
valuesofuthatwasseeninitscorrespondingGly data. All results are summarized in Table II. Comparison of
3
ff94gssamplesbotha andPP ,butlittlepopulationof standardMDandREMDdemonstratesthatrelativepop-
R II
extendedbispresent;a isalsoshiftedlowinw.Thea ulations in the four regions are essentially unchanged,
R L
conformation (with a positive backbone u value) is which indicates that the results are well-converged and
sampled most in ff99SB and least in ff03. This confor- thatthedifferencesdonotarisefrompoorsampling.
mationoccursinseveraltypesofb-turns,47particularly Simulations with ff99SB produce ensembles that are
those found in b-hairpins, and thus reproducing the largelyPP (38%)withsignificantfractionofb(20%)and
II
energy of a region may be important for modeling lower population of a (13%). Compared with other force
L R
b-hairpins. fields, this appears to be in best agreement with experi-
To better interpret such significant differences in rela- mental observations, in that PP is the most favorable
II
tivepopulationsofvarioussecondarystructureelements, conformationwithlowerpresenceofextendedbandeven
we turned to available experimental data that provide lowerfractionofa-helicalstructures.
more direct evidence concerning the dominant structure ff94 and ff99 show clear preference towards a (72%
R
ofshortalaninepeptidesinaqueoussolution.Allofthese and 36% respectively), with the relative fractions of a
R
observations have important implications for the varied whencompared withtheotherbasinsbeingover90%for
conformationalstatesofproteinsandthereforeshouldbe bothforcefields.Thelowerabsolutefractioninff99(36%)
capturedintheforcefieldifitshouldeverbeusedforthe arises from the unusually high population of structures
description of unfolded proteins or for protein folding around ((cid:1)1508,08) which donot correspond to any ofthe
studies. four canonical secondary structure basins and are also
A recent two-dimensional infrared spectroscopy study much less populated in other force fields. Overall, the
ofalaninedipeptideinaqueoussolution48estimatedback- data from ff94 and ff99 are both at variance with experi-
bone dihedral angles in the range ((cid:1)708 (cid:4) 258, þ1208 (cid:4) mentalobservations.
258) corresponding to PP -like conformation. Another ff03 behaves most similarly to ff99SB, but it may still
II
structural study of alanine dipeptide using 13C NMR49 be slightly over stabilizing a-helices as the most popu-
confirmed the presence of PP conformation, but also lated basin is the a conformation (34%). It is alsosome-
II R
suggests that a mixture of PP and a is more likely in what unsettling that ff03 samples very poorly the a
II R L
water.ThedominanceofPP intrialaninewasreiterated region that is important for b-turn conformations. We
II
by two-dimensional vibrational spectroscopy studies of obtained no sampling of a in standard 80-ns-long MD
L
Woutersen and Hamm50,51 and later modified by the simulations (0.0% fraction), indicating that this region is
same authors to include around 20% of a apart from at least 5 kcal/mol (this is based on our cutoff for free
R
PP .52TheexperimentalstudiesofSchweitzer-Stenner& energy plots, Fig. 3) less stable than the most populated
II
Eker on trialanines53,54 and tetraalanines55 in water a region. However, REMD simulations do sample a
R L
usingpolarizedRaman,FTIR,andVCDspectroscopycon- region with 0.06% population fraction. The sampling is
firmthedominanceofPP conformationintetraalanines, insufficienttomakeareliableestimate(seeuversustime
II
whilea 50:50mixture ofPP and extended b-strand-like plots in Figure S5 of supplementary material) of free
II
conformations was observed for trialanines. Other stud- energy difference between a and a , but based on the
R L
ies,56–58 mostly from Kallenbach and colleagues, have populations in those regions, a is at least 3.8 kcal/mol
L
confirmed that short alanine peptides form predomi- higher than a . Even though ff99SB only adopts 2–3%
R
nantly PP conformation, which is in temperature de- population in that region, it is still the largest compared
II
pendent equilibrium with extended b-strand conforma- withotherforcefields.
tions. Taken together, these experimental findings sug- Simulations with ff99/ result in essentially the same
gestthatshortalaninepeptidesarefoundpredominantly behavior as was obtained for ff94, with the relative ba-
in PP -like structures with varying degree of extended sinpopulationsoverwhelminglydominatedbya confor-
II R
b-strand, and possibly smaller fractions of a-helical con- mations (97 vs. 95% for ff94, see relative fractions of in-
formations. dividual secondary structures in Table II and free
We analyzed the Ala simulation data obtained from energy surfaces in Figure S6 of supplementary mate-
3
each force field in terms of fractional population of local rial). Thus, at least at the level of glycine/alanine tetra-
conformational basins corresponding to the four preva- peptide,thesmallchangeofasingleuterminff99/does
lent secondary structure elements. The frequent transi- not seem to result in significantly different or improved
tions between PP , a , and b regions result in relatively behavior.
II R
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
720 V.HORNAKETAL.
TABLEII. Fractions(Percentages)oftheFourSecondaryStructuresinAla MDSimulationsinExplicitWater
3
Region
Forcefield PP b a a
II R L
ff99SB 38.1(cid:4)1.2(51.9(cid:4)1.7) 20.1(cid:4)0.7(27.4(cid:4)0.5) 13.0(cid:4)1.6(17.8(cid:4)2.6) 2.1(cid:4)0.9(2.8(cid:4)1.2)
ff99SB(REMD) 36.9(cid:4)0.1(50.1(cid:4)0.3) 20.5(cid:4)0.8(27.9(cid:4)1.0) 12.7(cid:4)0.5(17.2(cid:4)0.7) 3.4(cid:4)0.1(4.6(cid:4)0.1)
ff03 28.0(cid:4)0.4(37.8(cid:4)0.4) 12.4(cid:4)0.5(16.7(cid:4)0.5) 33.6(cid:4)0.5(45.4(cid:4)0.9) 0.0(cid:4)0.0(0.0(cid:4)0.0)
ff03(REMD) 28.6(cid:4)0.6(38.2(cid:4)0.6) 12.6(cid:4)0.4(16.9(cid:4)0.5) 33.4(cid:4)0.7(44.7(cid:4)1.1) 0.1(cid:4)0.0(0.1(cid:4)0.0)
ff94 2.7(cid:4)0.5(3.5(cid:4)0.7) 0.8(cid:4)0.1(1.1(cid:4)0.2) 72.0(cid:4)0.5(94.6(cid:4)0.6) 0.5(cid:4)0.2(0.7(cid:4)0.3)
ff99 1.0(cid:4)0.0(2.5(cid:4)0.1) 2.4(cid:4)0.3(6.3(cid:4)0.8) 35.8(cid:4)0.4(90.8(cid:4)0.5) 0.1(cid:4)0.1(0.3(cid:4)0.3)
ff94gs 20.7(cid:4)3.1(30(cid:4)6.0) 2.3(cid:4)0.3(3.4(cid:4)0.7) 46.3(cid:4)6.8(66.0(cid:4)6.5) 0.4(cid:4)0.2(0.6(cid:4)0.2)
ff99/ 1.6(cid:4)0.7(2.2(cid:4)0.9) 0.4(cid:4)0.2(0.5(cid:4)0.2) 68.8(cid:4)0.5(96.6(cid:4)0.6) 0.4(cid:4)0.4(1.2(cid:4)1.2)
See Methods for definitions of secondary structure regions. Numbers in parentheses are relative fractions, only considering populations
withinthefourregions.
Decoyanalysis oflonger peptides much better than ff94/ff99 but, as we noted before, may
still be slightly too favorable for helices. ff94gs behaves
A critical property of the energy function is whether it
similarlytoff03.
can distinguish the native state such that it corresponds
On the other hand, all tested force fields successfully
to the lowest (free) energy. We previously employed
scorethea-conformationoftheBaldwinhelixwithlowest
‘‘decoy analysis’’12 as a rapid test of how well the energy
energy. As expected, the energy gap between the helical
function identifies the native structure. By plotting the
conformation and those with larger RMSD values is
potentialenergyagainstRMSDfromthenativestructure
greater for the force fields with strong helical bias (ff94,
one can readily see if the native state generally has low-
ff99,ff99/).Eventhoughff94gsperformedalmostaswell
estenergies.This hasthe advantagethatenergyprofiles
as ff99SB for trpzip2, it scores the Baldwin helix even
for multiple force fields can be compared without obtain-
more favorably than the ff94/ff99 force fields that are
ingconvergedensemblesforeachpeptide/forcefieldcom-
known to over stabilize a-helix conformations. Likewise,
bination.Asaresultofneglectingentropyandthe useof
all ofthe forcefieldscorrectly identifythe trpcagenative
an implicit water model, however, the results do not pro-
fold.
vide quantitative stability data for any of the peptides
Overall, both ff99SB and ff03 perform well in all three
simulated.
decoy cases. ff99/ follows very closely the behavior of
We show the ‘‘lowest profile’’ (see Methods) of energies
ff94.Sinceff94gsandff99/havenotbeenusedasexten-
plotted versus RMSD from the native structure (Fig. 4)
sively as the original Amber force fields upon which
with the expectation of seeing lowest potential energies
theyarebased(ff94/ff99),andbecausetheyhavenotper-
(including solvation free energy) for small RMSD values,
formed well on the tetrapeptide and decoy screens, we
reflecting that the native basin indeed can be distin-
willnotconsidertheminfurthervalidation.
guished based merely on the energy function. We per-
formed decoy screening using several peptides that have
Comparison withNMR parameters
been experimentally demonstrated to adopt stable sec-
for lysozyme and ubiquitin
ondary and/or tertiary structures. Trpzip2 represents
an exceptionally well-defined b-hairpin,42 stabilized by Hen egg white lysozyme has become a standard for
cross-strand pairs of indole rings. The Baldwin type ala- evaluating the quality of force fields by comparing inter-
nine-based peptide44 was designed to adopt a-helix. nal dynamics parameters calculated from MD simula-
Lastly, the trpcage mini-protein45 has a mixture of sec- tions to NMR relaxation experiments.59–63 The degree of
ondarystructureelements(a,310andpolyproline). backboneflexibilityisspecifiedbyexperimentallyderived
We previously used decoy analysis of the trpzip2 b- order parameters S2, which correspond to amide bond
hairpin to demonstrate the extreme helical bias of the N(cid:1)(cid:1)H librational motion. Lower S2 values reflect in-
ff94/ff99 force fields,12 showing that nonnative a-helical creasedbackboneflexibility.Wegenerated30-nsMDsim-
decoys had the lowest energies and were also the most ulationsoflysozymewithexplicitwater,usingfourdiffer-
populated in simulations. Therefore it comes as no sur- ent parameter sets:ff94, ff99, ff99SBand ff03.Order pa-
prisethatthepresentdecoyplotsagainclearlyrevealthe rameterscalculatedfromthesesimulationsarecompared
inadequacyoftheff94andff99forcefieldstoidentifythe with experimentally derived values in Figure 5(A). The
native state of the hairpin. Indeed the pronounced mini- RMSDandthecorrelationcoefficientbetweenexperimen-
mum around 4–5 A˚ [Fig. 4(A)] represents a large set of talandcalculatedS2aregiveninTableIII.Itcanbeseen
helicalstructures. Consistent withourobservationsfrom that the agreement of ff99SB derived values with the
the tetrapeptide simulations, ff99/ behaves much like experiment is very good, with the best correlation to ex-
ff94. The native trpzip2 structure scores as lowest in perimental data and only half of the RMSD compared
potential energy only with ff99SB. ff03 also performs withallotherforcefields.
PROTEINS:Structure,Function,andBioinformatics DOI10.1002/prot
Fig.3. Free energyu/w mapsfor Gly (toprow)andAla (bottom row)from 80-nssimulationswithexplicit TIP3Pwater.Theenergiesarecolor
3 3
codedfrom0upto5kcal/mol.PDBsurveydataisrepresentedbyasimpleRamachandranplot.Individualforcefieldsaredesignatedas:ff99SB(this
work),ff03,24ff94,5ff99,11andff94gs.13
Fig. 4. ‘‘Lowest energy’’ profiles for three prototypical decoy systems, each tested with six Amber force field variants (ff94, ff99, ff99SB, ff03,
ff94gs,ff99/):(A)trpzip2,whichisknowntoadoptab-hairpin,(B)aBaldwin-typesequenceasarepresentativeofana-helix,and(C)trpcageasa
representative of mixed secondary structures.RMSD values are calculated using theexperimentally determined structure asa reference.Ideally, a
forcefieldshouldshowlowestenergiesforthelowestRMSDvalues.
Fig.5. Order parameters(S2)derivedfromexperiment(blackline)andcalculatedfromMDsimulations inexplicitwaterwithdifferentforcefield
parametersets:ff99(green),ff99SB(red),ff94(blue)andff03(magenta).ErrorbarsreflecttheconvergenceofcalculatedS2values.(A)Lysozyme.
Secondary structures of lysozyme are labeled: helix A (HA: residues 4–15), loop 1 (L1: 16–23), helix B (HB: 24–36), strand 1 (S1: 41–45), turn 1
(T1:46–49),strand2(S2:50–53),strand3(S3:58–60),longloop2(L2:61–78),3 helix1(H1:80–84),loop3(L3:85–89),helixC(HC:89–99),
10
loop 4 (L4: 100–107), helix D (HD: 108–115), loop 5 (L5: 116–119) and 3 helix 2 (H2: 120–124). (B) Ubiquitin. Secondary structure elements in
10
ubiquitin:strand1(S1:residues1–7),turn1(7–10),strand2(S2:10–17),turn2(T2:18–21),helix1(H1:23–34),turn3(T3:37–40),strand3(S3:
41–44),turn4(T4:45–48),turn5(T5:51–54),helix2(H2:56–59),turn6(T6:62–65)andstrand4(S4:66–71).Differencesaredescribedinthe
text.