Folding simulation of Trp-cage utilizing a new AMBER compatible force ¯eld with coupled main ...

15 Pages·2014·5.53 MB·English

JournalofTheoreticalandComputationalChemistry
Vol.13,No.4(2014)1450026(15pages)
#.c WorldScienti¯cPublishingCompany
DOI:10.1142/S0219633614500266
FoldingsimulationofTrp-cageutilizinganewAMBERcompatible
force ¯eld with coupled main chain torsions
LirongMou*,XiangyuJia†,YaGao†,YongxiuLi†,JohnZ.H.Zhang†,‡andYeMei†,‡,§
*InstituteforAdvancedInterdisciplinaryResearch
m
o EastChinaNormalUniversity,Shanghai200062,P.R.China
c
ww.worldscientific.personal use only. Stat‡eNYKEUeay-sELtCaCIbnNh†oCsiUrtneiaatnCtuotNteerenryotorfeoomfSrfrhTaPfaLolhnrraeUegoscCnheriaireosviitmeaoi2crnnpa0sduil0Stt0yCapa,6netoi2cdSom,thnrCpPaoauonsl.tmcgCaRohtphp.aiuoeyiCtnma2ahat0inilisn0odtBa0rnyD6iao2leap,pSthPacyNri.setYiRmncUsc.eenCSthhoianfnaPghhayisics
m wFor §[email protected]
aded fro8/14/14. RAeccceeivpetded2120JManauracrhy22001144
o0
nln Published21May2014
wo
oY
DT
13. RSI A newly developed AMBER compatible force ¯eld with coupled backbone torsion potential
hem. 2014.K UNIVE trmeerpemclihcsaan(AiesxMmcBhwaEintRhg0ea32saDynn)deirsdgueirtteiicclitzfoemdldoiinnlegcaouflfaotrlhdeidnhygyndsairmmopiuchlsaotb(iiMocnDcoo)rfeas(imHmPuinlCai-t)pioarnnostd,etinhaeTm(cid:1)r-puh-letcila-ixgsteie.npTthhfoerlod¯uinngahgl
CR stageissuggested.Thenativestructurehasthelowestfreeenergyandthemeltingtemperature
mput. W YO psureredmicetendt.fTrohmisstthuedysp,etcoig¯ecthheerawtictahpoaucritpyreCvvioiussosntulydy1,2sKhohwisghtheratthAaMnBtEheRe0x3p2Deriismaennataclcumraeate-
oE
CN force¯eldthatcanbeusedforproteinfoldingsimulations.
or. by
he Keywords:Trp-cage;AMBERforce¯eld;coupledbackbonetorsion;foldingmechanism;melting
T
J. temperature.
1. Introduction
Understanding the mechanism of protein folding is critical for ¯ghting against
protein misfolding related diseases and for protein engineering. Much progress in
this¯eldhasbeenmadeinthepastfewdecades.1–3Asacomplementtoexperiment,
molecular dynamics (MD) is able to give an atomic description of the kinetic and
thermodynamic properties of proteins covering a wide range of time scales, which
are usually di±cult to observe in experiment. Some peptides and small proteins
have been studied4–7 by MD simulations bene¯ted from their small size and fast
folding rate. Among them, a 20-residue mini-protein Trp-cage8 (sequence: NLYIQ
§Correspondingauthor.
1450026-1
L.Mouetal.
WLKDG GPSSG RPPPS) is an ideal system that has been well investigated by
both theoretical9–22 and experimental means.23–29 Trp-cage consists of an (cid:1)-helix
from residue 2 to 8, a 310-helix from residue 11 to 14 and a polyproline II stretch.
The indole ring of Trp6 is buried in the center of a hydrophobic core (HPC), which
isformedbythesidechainsofTyr3andfourprolineresidues(Pro12,Pro17,Pro18,
and Pro19). The guanidine group of Arg16 and the carboxylic group of Asp9 can
form a salt bridge. Some suggest that salt bridge may also contribute to the
stability of protein structure.16,30–32 While some experimental and theoretical
studies showed that Trp-cage can maintain the folded state after the elimination of
this salt bridge.28,33 With continuous advancement of the force ¯elds for protein,
there is no di±culty in folding this short peptide. But to recapitulate some ex-
m
co perimental observables, such as melting temperature, is still a challenge. More
c.
worldscientifional use only. i(ocn1oft)nTr¯Wirgnpuhe-imcanatgegniestaitnsahd®aeerscmaitmtsoepr-tellhiemedtiwfito±iolndcgsuitnlsatgtteelpqamuniendosdsttceihoalenposperromoafcroeetrshesirsceoolffamoftloedplddilneixngt?gop?(r3io(t)2cs)HesWfoso?wlh3d4eti–nht3hg6eee(rm4nt)vheicWerhofanhonemldtisheimnnegtr.
ww.pers the secondary-like structures form before the protein collapsed to a compact
m wFor structure? Many studies have been devoted to ¯nd the answers. The NMR and
aded fro8/14/14. eCxDpeerximpeernitmsesnutgsg8e,2s7t,3e7dstuhgagtetshteedfoaldsinimgpmleecthwaon-isstmateofmtheechTarnpi-scma.geWphroilteeisnomwaesoetvheenr
o0
nln more complex and intermediates were involved before the occurrence of fold-
wo
DoTY ing.24,28,38After77MDsimulations,Duanetal.suggestedthatpackingoftheTrp6
13. RSI sidechaincouldbetherate-limitingstep.14,15Tianetal.34noticedthatthesize,the
4.E
m. 201UNIV ptooldaeritzeartmioinneotfhtehefocldonin¯gnepmatehnwtaayn.dZthhoeues®uegcgteosftesdolvtheantttchoem(cid:1)pe-hteedlixwwitahsefaocrhmoetdheinr
CheRK the ¯nal stage.39 On the contrary, some experimental results showed evidences of a
ut. YO helical structure in the denatured state of Trp-cage, suggesting an early formation
ompEW of the (cid:1)-helix.24
CN
or. by Thecontroversyisresultedfromthelimitationsofboththeexperimentsandthe
e
h computations. Theoretical prediction of protein folding pathway depends on the
T
J. ergodicityofphasespaceandthereliabilityofforce¯eldemployed.Duetotherapid
development of computational techniques, the urgent requirement for an accurate
force¯eldisbecomingmoreandmoresigni¯cant.Manystudieshavebeendevotedto
the re¯nement of force ¯elds. Best et al. optimized the additive AMBER40 and
CHARMM41 all-atom force ¯eld. Li et al.42 used to downhill simplex minimization
algorithm to optimize the dihedral angel potential with respect to NMR measure-
ments.Thisnewforce¯eld,AMBER99SBnmr,iscapableofgivingimprovedresults
compared to those of AMBER99SB. Lindor®-Larsen et al.43 optimized side-chain
dihedral parameters of certain residues by matching with the high-level quantum-
mechanicalcalculation.However,wedeemthatamajordefectincurrentforce¯elds
consistsinthesimplisticfunctionalformformainchaintorsions.Duethatthemain
chain torsions (’ and ) are not separable, employing coupled terms for the main
chain torsions is essential, especially for delineating the potential energy barriers,
1450026-2
FoldingsimulationofTrp-cage
which are crucial to large conformational change. A new set of backbone torsion
parametersutilizingcoupled2Dmainchaintorsiontermshasbeendevelopedinour
group.44 Preliminary examinations of its reliability for model peptides have shown
thatitiscapableofgeneratingmorebalancedsecondarystructuredistributionthan
original AMBER force ¯elds.
In this work, folding simulation of Trp-cage was carried out with this new
AMBER compatible force ¯eld incorporating re-optimized 2D main chain torsion
potential (AMBER032D). Other force ¯eld terms were directly extracted from
AMBER0345force¯eld.ThesimulatednativestatematchedverywellwiththeNMR
structureandthepredictedmeltingtemperaturewasclosetotheexperimentalvalue.
Meanwhile, a possible folding mechanism was proposed.
m
o
c
c.
m www.worldscientifiFor personal use only. 22Irdnai..h1tAeeM.ldyM2re.aDBtTlhsEhtoRfeoodrrpfmsooirteocednentb¯iayelldCe,n–teNhrge–yCmo(cid:1)af–inCeaccahhnadtinorNtsoi–orCnsi(cid:1)oi–snCsex–(pN’r,easnrsedesdp ea,cstwiavhe1iclDyh)Faaroreuerditeer¯renaetexedpdaanssesptiohane-
aded fro8/14/14. tpraurnacmaetteedriazfatteironthorfeteoorsriofonutrertmermmsa.inDlyuefotcousthoenlaimsmiteadllpnourmtiboenrooffsppaacraemareotuernsd, tthhee
o0
wnlon energybasis,whichisdeleteriousforthestudyoflargeconformationchangesuchas
oY
DT protein folding. Recently, a new functional form with 2D Fourier terms for main
13. RSI chain torsions has been proposed and parameterized by our group.44 The alanine
4.E
m. 201UNIV denipeergpytidmea(pAfDor)awsaesrciehsoosefn(’a;s t)hepamirosdwelassycsatelcmulfaotredpaartaMme0t6er2iXza/taioung.-cTch-pevptzo/te/nHtFia/l
CheRK 6-31G**levelusingGaussian09package.46ByequalizingtheQMandMMenergies,
ut. YO the main chain torsion energy can be expressed as:
mpW
CoNE E ¼E (cid:1)E þG (cid:1)G ; ð1Þ
or. by mct int oth PCM GB
e
h whereE andG aretheinternalenergyandpolarsolvationenergycalculatedat
T int PCM
J. QM level47 employing SMD, a continuous solvation model.48 It di®ers from our
previousstudyinwhichtheIEFPCMsolvationmodelwasemployed.E istheMM
oth
energy excluding the main chain torsions. G is the generalized Born solvation
GB
energy.49 E can be written as:
mct
XN’ XN
E ð’; Þ¼ Cðm;nÞeim’ein : ð2Þ
mct
i¼(cid:1)N’j¼(cid:1)N
The expansion coe±cient, the only parameters for the main chain torsion, can be
obtained by imposing transformation
1 XN’ XN
Cðm;nÞ¼ E ð’; Þe(cid:1)im’e(cid:1)in ð3Þ
4(cid:2)2 mct
i¼(cid:1)N’j¼(cid:1)N
1450026-3
L.Mouetal.
with Emctð’; Þ calculated from Eq. (1). Other torsion terms around N–C(cid:1) and
C(cid:1)–C bonds have their contributions also removed from the original AMBER force
¯eldandmergedinto2Dtorsionterms.Thisnewforce¯eldhasbeencodedintothe
sander and pmemd modules of AMBER 11.
2.2. Simulation protocol
The simulation began with a linear structure of trp-cage. It was ¯rst optimized by
1000 steps of steep decent method followed by 500 steps of conjugate gradient
method. The relaxed structure was gradually heated to 300K in 100ps. Then a
500ps MD simulation at constant temperature was conducted to further relax the
m whole system. The ¯nal structure was chosen as the initial structure for all the 12
o
c.c replicasinREMDsimulations.Temperaturesweresetinarangefrom261to542K.50
worldscientifional use only. GmbyeimnsueicrraftalhiczeeedasorBelvaoarttneiormmno.ed5®2eelI5cn1tt.weNgitrohanlaptniomlea®eressctotelipvvaewtsaiaoslntsecttoerntmoce1nwftasr.saTtaieopmnprpooefxr0iam.2tuaMrteewlwyaarsseuprterigelisuzelenadtteetddo
m www.For pers uristihnmg5B4ewreansdsuesnedthteormcoonstsatrta53inwaitlhl tahecocuopvlainlegnttimboencdosnsitnavnotlvoinfg1phsy.dSrHogAeKnEataolmgos-.
oaded fro08/14/14. Se5a0wc0ahnpssrewdpielrirececata.tSMtneDmappsstihmeodutselavwteireoyrne0s.s2aa5tvpe3ds0a0enKvderMwyeD0re.s2i5amlpsuosl.actBaioerrsniiseddweseorRuetEeMxtotDesntdsuiedmdyutoltah1tei6o0fnon,ldstiwfnoogr
nln
wo pathway.
oY
DT
3. SI
4.1ER 2.3. Analysis method
m. 201UNIV Parallel Tempering Weighted Histogram Analysis Method (PTWHAM)55 was used
heK for the calculation of free energy. The potential energy density of states can be
CR
ut. YO expressed as:
mpW P
heor. Coby NE (cid:1)m ¼ PNk¼T1g(cid:1)m1kNNkk¼(cid:2)T1Ug(cid:1)me1kxHp½mfkk(cid:1)(cid:3)kUm(cid:2); ð4Þ
T
J. where (cid:1) is the density of state and g is the statistical ine±ciency which usually
m mk
takesthevalueof1.(cid:3) istheBoltzmannconstantandf isthefreeenergyatthekth
k k
temperature. U is the energy and H is the number of conformations with the
m mk
potential energy valuesequal to U . N isall the conformations sampled atthe kth
m k
temperature.N isthetotalnumberofreplicas.Anyobservablesatthetemperature
T
of interest can be calculated from the density of state. The root mean square devi-
ation (RMSD) of the backbone atoms N, C(cid:1) and C and that of the HPC between
simulated and the NMR structure (the ¯rst model in PDB entry 1L2Y) were com-
puted.Otherreactioncoordinates(RCs),likeradiusofgyration(R ),nativecontact
g
(Q), heat capacity (C ) and (cid:1)-helix fraction were also analyzed. Native contact, a
v
widelyusedmetric,isde¯nedbytwonon-neighboringresiduesseparatedbylessthan
7.0(cid:3)A.56,57 There are 42 native contacts in the NMR structure. The formation of
(cid:1)-helixstructurewasdeterminedbymainchaindihedralangles((cid:1)100(cid:3) (cid:4)’(cid:4)(cid:1)30(cid:3)
1450026-4
FoldingsimulationofTrp-cage
and (cid:1)67(cid:3) (cid:4) (cid:4)(cid:1)7(cid:3)). To avoid numerical discontinuity problem, the helicity is
de¯ned as:
1
(cid:1) ¼ n h i on h i o: ð5Þ
h 1þ 1:5(cid:5)ð’þ60(cid:3)Þ 4 1þ 1:5(cid:5)ð þ45(cid:3)Þ 4
30 25
When ’ and equal (cid:1)60(cid:3) and (cid:1)45(cid:3), respectively, the helicity equals 1. When
dihedrals are away from ((cid:1)60(cid:3);(cid:1)45(cid:3)), the helicity decays rapidly. All the ensemble
statistical analysis was performed on the last 155ns trajectory of each simulation.
m 3. Results and Discussion
o
c
worldscientific.onal use only. Dcetiixoeucnnrhcisanynfgogfrotehrdreiap®RtheiraoEesMnebteDsttpweasmeicemepneusrlnaaaemttiuigoprhnleib,snogaarrianwengisddhteoetwmrhaenpnegcirnoeantFovufiergert.gseS.men1Tp.cheeTreahstppeueorreetedesanrwtweiaeearlrdeeeengqueuusreaagdrtya.enTdotivheseeterrdlieab±bpuys--
ww.pers between neighboring temperatures and the exchange ratios are around 0.3. All the
m wFor replicashavetravelledinthetemperaturespaceback-and-forthformanytimes,and
aded fro8/14/14. tahnedfSo3ld.eWdhsteantetshehafovledebdeesntreuvcetrurreeaicshreedacbhyedm,othsteopfetphteidreepstliacyass,inasthsheovwicniniintyFoigfst.hSa2t
o0
nln structure and at low temperature for several to tens of nanoseconds, indicating the
wo
DoTY potential energy superiority of the folded structure under this interaction potential.
3. SI The distribution of the backbone RMSD from the ¯rst model of the NMR
1R
4.E
m. 201UNIV sFtirgu.c1tu(ar)e.sA(Pcroont¯eginurDataiotan iBsacnoknsiednetrreyd1tLo2bYe)inshtohwesuntfhorledeedpestaaktse aifsitdsebpaicctkebdonine
heK RMSD is above 2.2(cid:3)A. This criterion is determined by the location of the distinct
CR
ut. YO troughinthedistributionofbackboneRMSDandisconsistentwiththatinthestudy
mpW by Day and coworkers.19 There are two narrow peaks in the folded domain, which
oE
CN
or. by cover 63.3% of the con¯gurations altogether. In the unfolded region, a broad dis-
he tributionpeaksat3.25(cid:3)Aandextendstoover8(cid:3)A.Thispeakpossiblycorrespondsto
T
J. an intermediate state with the peptide partially folded. Cluster analysis is an e®ec-
tive means to detect stable states. We used the kclust tool in MMTSB58 to classify
theconformationssampledat300K.Thecentroidstructuresofthetopthreeclusters
are shown in Fig. 1(b). The largest cluster holds 57% of the conformations. The
backbone RMSD of the centroid structure is only 1.1(cid:3)A, which is de¯nitely in the
folded region. The peptide adopts a U shape and both the (cid:1)-helix and the 310-helix
are well folded. Superimposition of the simulated conformation with the smallest
RMSD and the NMR structure is shown in Fig. 1(c). These two structures are well
alignedforboththetraceofthebackboneatomsandthepackofthesidechains.The
TRPresidueiswellburiedinaHPC.Thesecondandthethirdclusterscontain18%
and 15% of the conformations respectively. The centroid structures of these two
clustersareinanearUshapebutneitherofthemhas(cid:1)-helixcompletelyformed.We
also studied the positions of Trp-6 residue in the centroid structures of these three
1450026-5
L.Mouetal.
m
o
c
c.
worldscientifional use only.
ww.pers
m wFor
aded fro8/14/14. Ftriagj.ec1t.ory(aa)tP30r0obKa;b(icli)tiOesveorflatyheofRthMeSsDimduilsattreidbuctoinofnosrmatat3io0n0Kw;it(hbt)hMesamjoarllecsotnRfoMrmSaDtiaotn3c0l0uKste(rrsedo)fatnhde
o0
nln theNMRstructure(blue).
wo
oY
DT
13. RSI mainclusters.Inthe¯rstcluster,theHPChasformedandTrp-6isstable.Butinthe
4.E
m. 201UNIV odtrhaemractilcuasltleyr.s,TthhisessuidggeecstheadinthoefnTurcple-6atitoank-ecsonadwenrsoantgionormienecthatainoinsmanadnd°uincdtuicaatteeds
CheRK one possible folding pathway: in the early stage, the tertiary and second structures
ut. YO areonlypartiallyfolded;thenthefoldinggoesviathepackofTrp-6inHPCfollowed
ompEW by the full formation of (cid:1)-helix.
CN
or. by In order to con¯rm this conjecture, we plot the free energy landscapes (FEL)
e
h mapped to several RCs. The FEL shown in Fig. 2 also supports the collapse-¯rst
T
J. mechanism. The HPC collapses in the early stage of the folding with the RMSD of
theHPCdecreasesfrom16to5(cid:3)A.Thenthe(cid:1)-helixbeginstogrow.Aswecanseein
Fig. 2, there are two barriers separating this FEL into three regions. In region III,
TheRMSDofHPCisnear4.0(cid:3)Aand(cid:1)-helixfractiondistributesfrom0.25to0.65.In
regionII,theRMSDofHPCtranslatesfrom4.0(cid:3)Atonear2.5(cid:3)Aandthedistribution
of (cid:1)-helix varies from 0.4 to 0.65. As we can see, from region III to II, further
formationof(cid:1)-helixdoesnotoccuranditistheoptimizationofHPCthatdrivesthis
folding process. In region I, the RMSD of HPC is near 1.0(cid:3)A and (cid:1)-helix fraction
distributes from 0.5 to 0.8. So from region II to I, the growth of (cid:1)-helix and further
optimization of the HPC become synergetic. The (cid:1)-helix fraction of the native
structure is around 0.74, which is in the deepest free energy well.
WealsocarriedouttwodirectMDsimulationstostudythefoldingprocessofTrp-
cage. Each simulation ran for 500ns at 300K. For each trajectory, conformations
1450026-6
FoldingsimulationofTrp-cage
m
o
c
c.
worldscientifional use only.
ww.pers
m wFor
aded fro8/14/14. Fig.2. Fittedfreeenergylandscapeinkcal/molmappedtotheRMSDfromtheNMRstructureforthe
nlon 0 HPCandthefractionof(cid:1)-helixat300K.Theredarrowpointstothelocationofthenativestate.
wo
oY
DT
3. SI collected in the last 250ns are used for analysis, in which the native state has been
1R
4.E reachedandfolding/unfoldingprocesseshaveoccurred.InFig.3,theleftpanelshows
m. 201UNIV the variation of RMSD over simulation time for the ¯rst trajectory. After cluster
heK analysis,conformationsbelongingtothetop¯veclustersarecircledbydi®erentcolors.
CR
ut. YO
mpW
oE
CN
or. by
e
h
T
J.
Fig.3. TheRMSDvariesduringthelast250nsforthe¯rstdirectMDtrajectoryandeverycoloredcircle
encompassestheconformationsbelongtothesamecluster(leftpart).Thespecialstructuresnearestto
centersoftop¯veclustersandfoldingroutes(blackarrows)inthistrajectoryaredisplayed(rightpart).
1450026-7
L.Mouetal.
Eachclusterhasatypicalstructurenearesttoitscenter.These¯vetypicalstructures
arelistedintherightsideofFig.3.Aswecansee,forstructures3,4and5whoseRMSD
values are larger than 3(cid:3)A, their secondary structures are partially formed. But the
indoleringofTrp-6takesthewrongorientation.Thesestructuresmaycorrespondto
theconformationsinregionIIIofFig.2.Toarriveatthenativestate,thestructuresin
cluster 3, 4 and 5 must travel via cluster 2. RMSD of conformations in cluster 2
°uctuatelargelyfrom1.5(cid:3)Ato3(cid:3)A.Inthiscluster,Trp-6takesthecorrectorientation
whichfurtherpromotestheformationofHPCbutthe(cid:1)-helixisstillpartiallyformed.
ThiscorrespondstoregionIIinFig.2.HPCcanstabilizethesecondstructureandthe
second structure also contributes to the formation of HPC. This synergetic process
drives protein going from cluster 2 to cluster 1 (the native state). For the second
m
co trajectory,resultsaredescribedinFig.4.Inthistrajectory,structure1correspondsto
worldscientific.onal use only. pmcthoaumerstntpiaalgeltlotyievlfteyohrsrmtfooauertgdmeh..eFSdsrttorrbuumucctttuuunrTrfeeorlp22d-.e6hdIantssatsTakttrereuspct-t6otuhttrehaeekw5en,rnaoatntilhvlgeethosrteriagiesthneetc,taosottnrriidouenncst.ttuarSutrieoco,t3nut(rahRenesMsdehStaDhsveee>co(cid:1)a3n-lm:hd0eao(cid:3)Alrisxyt)
ww.pers structure elements cannot be stabilized by HPC and Trp-cage must leave this well
m wFor eventuallytosearchforamorestablestate.
aded fro8/14/14. witIhnstohmisesetxupdeyr,imaehnetlimxefraasgumreemnetnitssf.2o4u,5n9dIninsothmeeuontfhoeldrecdomstpauteta.tTiohnisalisstcuodniseiss,teXnut
o0
nln et al.17 used REMD with hybrid Hamiltonian to study the folding mechanism of
wo
DoTY Trp-cage and Chowdhury et al.14 carried out 77 100-ns direct MD simulations for
13. RSI Trp-cage,theybothfoundexistencesofhelixstructureintheunfoldedstateandthe
4.E
m. 201UNIV formation of HPC centered on Trp-6 was the rate-limiting step. Similar to their
heK
CR
ut. YO
mpW
oE
CN
or. by
e
h
T
J.
Fig.4. TheRMSDvariesduringthelast250nsfortheseconddirectMDtrajectoryandeverycolored
circleencompassestheconformationsbelongtothesamecluster(leftpart).Thespecialstructuresnearest
tocentersoftop¯veclustersandfoldingroutes(blackarrows)obtainedfromthistrajectoryaredisplayed
(rightpart).
1450026-8
FoldingsimulationofTrp-cage
m
o
c
c.
worldscientifional use only.
ww.pers
m wFor
aded fro8/14/14.
o0
nln
wo
oY
DT
3. SI
1R
4.E
m. 201UNIV
heK
CR
mput. W YO Fig.5. Fittedfreeenergylandscapeinkcal/molmappedtobackboneRMSDandRgat261K(upper
oE left),280K(upperright),300K(lowleft)and321K(lowright).
CN
or. by
e
h studies,wealsosuggestthatcorrectorientationofTrp-6andformationofHPCisan
T
J. important step in folding from the intermediate state to native state.
TheFELofTrp-cageatfourtemperaturesmappedtoR andRMSDareshownin
g
Fig. 5. The dominant free energy well at low temperature is located at the folded
region with the backbone RMSD and R around 1.0(cid:3)A and 7.0(cid:3)A, respectively. A
g
minorwellisseenintheunfoldedregionwiththebackboneRMSDaround3(cid:3)Aanda
more compact structure than the native structure. Population shift from the folded
region to the unfolded region can be detected with the increase of temperature,
indicating a thermal unfolding event. The conformations are mainly distributed in
the free energy well in the unfolded region at 321K. Plot of the population in the
folded region at all the temperatures is shown in Fig. 6(a) (see the black dots). The
predicted melting temperature, atwhich thefolded and unfoldedconformations are
equally populated, is 287K, which is much lower than the experimental measure-
ment(317K).26However,thepredictionmaybebiasedifonlyasingleRCischosen
1450026-9
L.Mouetal.
m
o (a) (b)
c
c.
worldscientifional use only. Ftteriimga.np6ge.lrea)t(ulaer)negVrtahanrogiafet;tiho(bne)snToafthitvehe¯e(cid:1)tpt-eohdpelusixplaertceiiog¯nicoonhfeafanotlddce(ubdrlvusete.astqeudaertee)rmthieneodccbuyrr(ebnlcaeckofdtohte)nRaMtivSeDco<n2ta:2ct(cid:3)Ai,n(trhede
ww.pers to delineate the thermodynamics.60 Therefore, the melting curves of the relative
m wFor lengthofthenative(cid:1)-helixregionandoftheoccurrenceofthenativecontactarealso
aded fro8/14/14. tshtuadnietdh.aTtohfetmheelgtilnogbatlemunpfeorladtinugre.Tofhtehnea(cid:1)t-ihveelicxonretgaicotnisise3v5e4nKle,sswlhaibcihlei,somfuwchhiclahrtgheer
o0
wnlon melting temperature is another 22K higher. Therefore, the predicted melting tem-
oY
DT perature shows strong dependence on the RC depicting the thermodynamics. An-
3. SI
1R other way to obtain the melting temperature is via the temperature-dependence of
4.E
m. 201UNIV sepneecrig¯ycahseatcapacityCv,whichcanbecalculatedthroughthe°uctuationofthetotal
heK
ut. CYOR C ¼ ðhE2i(cid:1)hEi2Þ:
mpW v RT2
oE
CN
or. by Asthetemperatureapproachesthemeltingtemperaturefrombothdirections,there
The willbeaclearspikeinCv theoretically.The¯ttedcurveofCv isshowninFig.6(b).
J. Optimized parameters through Monte Carlo indicate a melting temperature of
329K, which is only 12K higher than the experimental measurement. Zhou et al.39
used OPLS force ¯eld combined with explicit water model to study the melting
temperature of Trp-cage by REMD simulation. The reaction coordinate they chose
wasnativecontactandtheobtainedmeltingtemperaturewashigherthan400K.In
our study, melting temperature corresponding to native contact is about 339K.
Piteraetal.9usedAMBER94force¯eldtostudyTrp-cageandtheheatcapacitywas
used to study folding-unfolding transition. They found a clear pike of heat capacity
from 373K to 433K. Our calculation shows a better result by using heat capacity.
The FEL mapped to backbone RMSD and the distance between the titratable
groupsinArg16andAsp9isshowninFig.7.Itcanbeseenthatthesalt-bridgehas
anevendistributioninbothfoldedandunfoldedstateintheregionwithRMSDfrom
1(cid:3)A to 3(cid:3)A. Just as Zhou has pointed out,39 this salt bridge can stable not only the
1450026-10