Protein folding dynamics via quantification of kinematic energy landscape

0.19 MB·English

by  Sëma Kachalo

Protein folding dynamics via quantification of kinematic energy landscape
S¨ema Kachalo, Hsiao-Mei Lu, and Jie Liang∗
Department of Bioengineering, MC-063, University of Illinois at Chicago, Chicago, IL, 60305, U.S.A.
(Dated: February 4, 2008, Accepted byPhys. Rev. Lett.)
We study folding dynamics of protein-like sequences on square lattice using physical move set
that exhausts all possible conformational changes. By analytically solving themaster equation, we
followthetime-dependentprobabilitiesofoccupancyofall802,075conformationsof16-mersover7-
ordersoftimespan. Wefindthat(i)foldingratesoftheseprotein-likesequencesofsamelengthcan
differ by 4-orders of magnitude, (ii) folding rates of sequences of the same conformation can differ
6 by a factor of 190, and (iii) parameters of the native structures, designability, and thermodynamic
0 properties are weak predictors of the folding rates, rather, basin analysis of the kinematic energy
0 landscape defined bythe movescan provideexcellent account of the observed folding rates.
2

n
a The dynamics of protein folding has been studied ex- it was necessary to cluster conformations of larger sys-
J tensively[1,2]. Aremarkableobservationisthatprotein temsintomacrostatestoreducethesizeofthetransition

2 folding rates are well correlated with their native struc- matrix [8], therefore making the use of physical moves
1 turalproperties[1]. Anativecentricviewpostulatesthat infeasible. Here we directly solve the eigenvalue problem

proteinfoldingratesarelargelydeterminedbythetopol- of the 802,075 × 802,075 transition matrix and develop
]
M ogy of its native structure [3]. Theoretical models using a method to monitor the time-dependent probability of
G¯o potential where only native contacts contribute en- occupancy of all conformations simultaneously from the
B
ergetically are very successful in reproducing observed firstkinetic moveuntil reachinghalf-equilibriumconcen-
.
o folding rates [2, 4]. tration over 7-orders of time scale.
i However, the extent to which native structure deter- Our results show that the properties of native struc-
b
- minesfoldingrateremainsunclear. Bythenative-centric ture, designability, and thermodynamic properties are
q
view, different sequences for the same protein struc- inadequate to explain protein folding dynamics in our
[
tural fold would all have very similar folding rates, as modelsystems. Wefoundthatprotein-likesequencescan

1 they share essentially the same native structure topol- foldintothesamenativestructurewithfoldingratesdif-
v ogy. However, this is not consistent with experimental ferasmuchas190timesandsequencesofthesamelength
8
results. Asthe foldingratesofsimple single-domainpro- and energy gap can differ by 4-orders of magnitude in
1
0 teins differ by 6 orders of magnitude [3], folding rates folding rate. Instead of thermodynamic properties, we
1 maybeveryheterogeneous. Arecentexperimentalstudy showthatpropertiesofthe move-connectedenergyland-
0 showed that a designed artificial protein with no homol- scape defined by the connection graph of physical moves
6 ogoussequenceinnature thatadopts the same structure can provide excellent account for observed folding rates.
0
as a natural protein can fold 4,000 times faster [5]. A Model. We use the following energy model for dif-
/
o distinct possibility is that the empirical correlation be- ferent types of nonbonded HP contacts: E = 1,
HH
i
b tween properties of native proteinstructures and folding EHP = 0, and EPP = 0. By evaluating the energy
- ratesmayarisefrominadequatesamplinginthesequence level of all 216 sequences of 16-mers on all enumerated
q
space due to accumulated biased natural selection and |Ω| = 802,075 conformations, we have identified 26 se-
:
v limited genetic drift, rather than from intrinsic physical quences that all fold into the same ground state confor-
i
X properties of proteins. mation (Fig. 1). This set of sequences forms the largest
r Inthisletter,weusetwo-dimensionalhydrophobicand proteinfamily,whereeachsequenceadoptsthesamecon-
a polar (HP) lattice model [6] to study the relationship of formation, and all are connected by (a series of) point
folding rates, native structure topology, thermodynamic mutations. Altogether,thereare1,539foldablesequences
properties, and effects of sequence variation. We model with unique ground state conformations. There are 456
the physical movement of protein chains. Real protein conformations that are the unique ground state for 1 or
cannot immediately jump from one conformation to an- more foldable sequences.
other arbitrary conformation. Two conformations of the Wedevelopasetofphysicallypossibleprimitivemoves
same energy may be well separated kinetically. We re- (Fig. 1c). They are generalizations of corner move,
gard protein movement as a sequence of successive con- crankshaft move, and pivot move. We exhaust all possi-
formationalchanges,eachrepresentedbyaphysicallyre- bleoccurrenceofsuchmovesforeveryconformation. We
alizableprimitivemove. Thephysicalmovesetwedevel- verified that this move set is ergodic, i.e., all conforma-
oped exhausts all possible conformational changes for a tions areconnected to eachotherby a seriesofprimitive
structure. We use master equation to study the folding moves. With this move set, the simple energy scheme
dynamics offoldable sequences of length 16. While mas- of the HP model leads to a complex energy landscapes,
ter equationprovidesanexactsolution[6, 7],inthe past with numerous local minima for a foldable sequence.
2
infeasible,λ andthe correspondingeigenvectorsn and
1 1
v can be computed by an Arnoldi method.
1
A
B (a) (b)
C f
D h g
GEF a c b d e 1 e−02
H
KLJI 1 e−03
MONP 2 4 ding rate 1 e−04
Q ol
RST 1 (c) F 1 e−05
U
V
WXY 3 5 1 e−06
Z 5 6 7 8
Contact energy
FIG. 1: Protein-like sequences and the set of primitive moves. FIG. 2: The correlation of logk and ground state contact en-
f
Thelargestproteinfamilycontains(a)26sequencesthatallfold ergy. Acirclerepresentsoneofthe26sequencesshowninFig.1,
into(b)thesamestructure. HerefilledcirclesareHresidues. (c) acrossrepresentsoneofthe79singletonsequences,andthedi-
Themovesetincludes: among(1,2,and3),singlepointmoves amondrepresentstheG¯omodel. Nativeconformationsforafew
rotate around a single point; between (1 and 4), generalized sequences are also shown.
corner moves reflect around a diagonal axis connecting any two
residues;between(1and5),generalizedcrankshaftmovesreflect Thermodynamics and folding rates. Severalther-
aroundahorizontalorverticalaxis. Pointsofrotationareongray modynamic properties have been proposed to be deter-
background. Foragivenconformation,weexhaustivelysearchall minants of protein folding rates. We found that protein
possible position for point moves, all possible pairs of positions
stabilityasmeasuredbythetotalcontactenergyarecor-
forpossiblegeneralizedcornermovesandgeneralizedcrankshaft
relatedwithlogk (R2 =0.71),i.e.,morestableproteins
moves. f
fold slower in general (Fig. 2). Because stable proteins
WeuseMetropolis-typeofdynamicstoassignthetran- have lower ground state energy, some local minima will
sition rate r from conformation i to a neighbor con- also have relatively deep energy traps. As a result, more
ij
formation j connected by a move: r = 1 if E(j) ≤ stableproteinswillhaveslowerfoldingratesbecausethey
ij
E(i); r = e−[E(j)−E(i)]/T if E(j) > E(i); and r = can be trapped in such local minima. However,the fold-
ij ij
−P r , if j = i. For non-neighbors, r = 0. We ing rates of sequences of the same ground state energy
i6=k ik ij
assumethe effects ofviscosityandfrictionarenegligible. can still differ as much as 104. The heterogeneity of
We follow [7, 8] and use a master equation to study folding rate was already noted in an earlier study using
protein folding dynamics. Let p (t) be the probability macrostate approximation [6]. Here we found that even
i
that the HP molecule takes the i-th conformationat the sequencesthatfoldinto the sameconformationshownin
timet,thendp (t)/dt=P [r p (t)−r p (t)]. Written Fig.1ademonstrateawiderangeofrates,from1.1×10−3
i i6=j ji j ij i
in vector form, we have: dp(t)/dt = Rp(t), where R is to5.8×10−6,whichismuchlargerthanthedifferencebe-
the rate matrix whose entries are defined by the above tween the averagefolding rates for sequences of different
expression. We choose temperature T = 0.2 in unit of native state energies. Protein stability therefore provide
∆E /k , which is below the folding temperature T some but not the main explanation of the heterogeneity
HH B f
when50%ofmoleculestakethe nativeconformation. T of folding rates.
f
varies from ∼0.2 to ∼0.5 for different sequences. Energygapbetweengroundstateandexcitedstatewas
A generalsolutionofthe master equationcanbe writ- thought to be the necessary and sufficient determinant
ten as p(t) = PiCinie−λit with Ci = vTi p(0), where of folding rate [9]. For all 1,456 protein-like sequences
λ is the i-the eigenvalue of the rate matrix R, n the of N =16, the energy gap between the lowest state and
i i
corresponding right eigenvector, v the left eigenvector, the next state is ∆E = 1. The diversity in folding rate
i
and p(0) the initial vector of distribution of conforma- k shown in Fig. 2 clearly indicates that energy gap is
f
tions. In this study, we use the high temperature con- not a determining factor for the folding rate. The corre-
dition and assign p(0) = 1/|Ω|. Two eigenvalues are of lation R2 between logk and energy gap normalized by
f
particular interest: λ = 0 corresponds to the equilib- standard deviation is 0.01.
0
riumBoltzmanndistribution,andthesmallestnone-zero Anotherthermodynamicpropertythoughttobeanim-
eigenvalueλ determines the slowestmode ofrelaxation. portantdeterminantoffoldingratesisthecollapsecoop-
1
Following [8], we take λ as the folding rate k of the erativityσ =1−T /T [10],whereT isasdefinedearlier,
1 f f θ f
protein. Althoughthefullcomputationofalleigenvalues andT thetemperaturewhenheatcapacityC(T)reaches
θ
and eigenvectors for a 802,075×802,075 matrix M is its maximum. Fig. 3 shows that for the 26 sequences
3
that fold to the same native structure in Fig. 1b, there limitednumber(∼30)ofconformations. Wecanidentify
is a weak correlation(R2 =0.38)between collapse coop- states that are local minima, i.e., all states connected to
erativity and logk . Large variance in observed folding whichbymoveshavehigherenergy. Asimplecharacteri-
f
ratesexistforsequencesofsimilarcollapsecooperativity. zationofthekinematicenergylandscapeisthenthenum-
ber count n of the local minima. Fig. 3b shows that
min
an excellent correlation of logk and n (R2 = 0.85)
f min
−04 −04 can be found for the 26 HP sequences that fold into the
e 5 e e 5 e same conformation.
at at Our conclusions are not sensitive to temperature T.
folding r 5 e−05 folding r 5 e−05 TWhfreonmT3i0s0rKaisteod3fr1o5mK)0,.2w0etofo0u.n2d1(tehqautivtahleenfotltdoinrgairsaintge
k ofthe26sequencesallincreases. Althoughk forslow
5 e−06 0.1 0.2 0.3 0.4 0.5 5 e−06 100 200 300 400 500 fooffld1e.4r fionrcrfeaasstesfomldoerrse),(bky aatfaTct=or0o.f212.i0svweerlslu-csfoarrefalactteodr
collapse cooperativity number of local minima f
withk atT =0.20. Thecorrelationcoefficientsoflogk
f f
FIG. 3: Examples of the correlation of folding rate k with with the number of local minima, collapse cooperativity
f
thermodynamic properties and kinetic landscape properties. (a) (Fig.3),andotherthermodynamicparametersareessen-
k and collapse cooperativity σ have weak correlation (R2 = tially unchanged.
f
0.38). (b) kf has excellent correlation with the numberof local Time evolution and basin analysis. Monitoring
minima (R2 = 0.85), a property of the kinetic landscape. The the exact time evolution of all conformations simulta-
diamond represents the G¯o model.
neously until reaching equilibrium during folding is a
The number of sequences that take a specific confor- challenging task. Mathematically, the model of master
mationastheuniquegroundstateisthoughttobecorre- equation is equivalent to a Markov process, where the
latedwithoverallproteinstability andfolding rates[11]. population vector of conformations at time t + k∆t is
We calculated in addition k for a group of 79 single- given by p(t+k∆t)=Mkp(t), where M =I +R·∆t,
f
ton sequences with no sequence homologs that fold to I being the identity matrix. However, the k-time step
the same native conformations. The distribution of k s Markov matrix Mk rapidly becomes a dense matrix,
f
for the singleton sequences and the 26 sequences shown andfollowingthe time evolutionoffolding by a straight-
in Fig. 2 demonstrate similarly large variation. For our forward matrix multiplication of O(|Ω|3logk) steps be-
model, designability is not an important determinant of comesimpossibleforalargematrixofsize|Ω|=802,075
the folding rates. and k = 106 −1010. The analytical solution of p(t) =
TheInverseParticipationRatioI iscommonlyusedto PiCinie−λit through diagonalization is also impracti-
characterizethelocalizationofeigenvectors. Itisdefined cal, as it is only possible to calculate a few eigenvectors
as I = P v2, where v is the k-th coefficient of the and eigenvalues for a large matrix.
k k k
normalized eigenvector. The correlation between I for We seek an accurate solution without the approxima-
the equlibrium eigenvector and the folding rate for the tion of macrostates. Taking advantage of the sparsity
26 sequences is rather poor (R2 =2×10−3). of the rate matrix R, we follow the approach of Sidje
Kinematic determinants of folding landscape. [12] and use the analytical solution of matrix exponen-
Protein folding kinetics are intrinsically determined by tial: p(t)=eRtp(0), where eRt is defined by the Taylor
physical movement of molecules. Weak correlations of expansioneRt =I+tR+t2R2+···+tkRk+···. This
2 k!
the folding rate with thermodynamic properties are not expansion itself is impractical, as it also involves large
surprising. Thermodynamicpropertiesofasequencecan matrixproductofincreasingdensity. Plus,the entriesin
becalculatedifitscompletesetofconformationsareenu- the matrix terms may have alternating signs and hence
merated. Such properties are not affected by the kinetic cause numerical instability. A better approach is to ex-
connections between conformations. A smooth energy pand eRtp(0) in the Krylov subspace K defined as:
m
landscape ensuring fast folding can be easily permuted
into a rugged landscape by assuming different transi- Km(Rt,p(0))≡Span{p(0),··· ,(Rt)m−1p(0)}. (1)
tion rules between conformations. Both will have the
Denoting || · || as the 2-norm of a vector or
same thermodynamic properties, but the resulting fold- 2
matrix, our approximation then becomes p(t) ≈
ingratesforthesamesequencewillbeverydifferent. The
energylandscapeoffoldingisdictatedbythe connection ||p(0)||2Vm+1etHm+1e1, where e1 is the first unit ba-
graph of states defined by the move set. Characterizing sis vector, Vm+1 is a (m + 1) × (m + 1) matrix
suchkinematicenergylandscapeisthereforeessentialfor formed by the orthonormal basis of the Krylov sub-
studying protein folding dynamics. space, and Hm+1 the upper Heisenberg matrix, both
Although the energy landscape contains 802,075 con- computed from an Arnoldi algorithm. The error can
formations, each is connected by the move set to only a be bounded by O(em−t||R||2(t||R||2/m)m). We now
4
only need to compute explicitly eHm+1t. Because m Tounderstandwhethertheformationofcertainnative
is much smaller than 802,075, this is a simpler prob- contacts facilitate folding, we examine the time evolu-
lem. A special form of the Pad´e rational of poly- tion of each of the 8 native contacts (a–h) in Fig. 1(b)
nomials instead of Taylor expansion is used for this for the 26 sequences. We found that fast folders have
[N12]:(tHetHm)+1=≈PNppp(tcH(tmH+1)/N)pkp(a−ntdHcm+1=), cwhere· lloagrgkefr),amanodunctonntaatcivteccoantttahcet tdra(nRs2ie=nt0t.i7m4e−of05.801−w1i0th0
pp m+1 k=0 k m+1 k k−1 (Fig. 4), indicating that these contacts are critical for
p+1−k . In our calculation, we select m=30.
(2p+1−k)k foldingbyrestrictingfavorablytheconforamtionalsearch
Fig. 4 shows
space. The formation of other native contacts seem not
an example of
to be directly related to folding rates.
an HP sequence
Toconclude,westudiedproteinfoldingdynamicsusing
(sequence C in 1
0 amodelbasedondetailedphysicalmovesandexactsolu-
Fig. 1a) and the e−
1 tionof the master equation. We found that folding rates
time evolution
varyenormouslyforsequencesofthesamelength,energy,
ocofnfoirtmsationnative opulation1e−03 emnaetrigoyn.gaIpn,caonndtreavsetntoofththeethsaermmeogdryonuanmdicstpaateracmonetfeorrs-
and several local P
which are weak predictors of folding rates, properties of
minima confor-
5 the kinematic landscape defined by the physical moves
mations. The −0
1e provideexcellentcorrelationwithfoldingrates. Withthe
time evolution
computationoftimeevolutionofindividualconformation
of the native 1e−02 1e+00 1e+02 1e+04
Time from the first move to half-time of equilibrium, we show
conformation
that many transiently accumulating intermediate states
shows an initial
FIG.4: Thetimeevolutionofthenative can be identified by basin analysis.
fast phase upto
state and several local minima states. We thank Drs. Ken Dill, Bosco Ho, Xiaofan Li, Banu
t ∼ 50 time
The probability of occupation of na-
Ozkan, Dev Thirumalai, Jin Wang, and Weitao Yang
units. In prin- tive state conformation (top) increases
for helpful discussions. This work is supported by NSF
ciple, the local monotonically through a time span of
minima con- 10−2−105,butlocalminimaconforma- DBI0133856,NIH GM68958,and Whitaker TF-04-0023.
formations can tions go through transiently accumulat-
follow different ing intermediate states. ∗ Corresponding author. Email: [email protected]
kinetic processes: Some could be transiently accumu-
lating, and others either monotopically accumulating
or monotopically decreasing. Based on the computed
trajectories of time evolution, we find that the dynamic
behavioroflocalminimaconformationscanbepredicted
from basin analysis of the move-connected energy land- [1] K.W. Plaxco, K.T. Simons, and D. Baker, J. Mol. Biol.
scape. We define the size of the basin associated with 227, 985 (1998).
[2] O.V.GalzitskayaandA.V.Finkelstein,Proc.Natl.Acad.
each local minimum state i computationally by artifi-
Sci. USA 96, 11299 (1999). V. Munoz and W.A. Eaton,
cially making every local minimum an absorption state,
Proc.Natl.Acad.Sci.USA96,11311(1999).E.Almand
i.e., a sink of infinite depth, such that once reached,
D. Baker, Proc. Natl. Acad. Sci. USA 96, 11305 (1999).
no molecule can escape. This is achieved by assigning [3] B. Gillespie andK.W.Plaxco, Annu.Rev.Biochem. 73,
rij = 0 and rii = 1 for each local minimum state i 837 (2004).
[13]. p′(t = ∞) therefore reflects the size of the basin [4] T.R. WeiklandK.A.Dill, J.Mol. Biol. 329, 585(2003).
i
of the i-th local minimum. We define the accumulation [5] M. Scalley-Kim and D. Baker, J. Mol. Biol. 338, 573
ratio as ̺ = p′i(∞) . If ̺ > 1, state i is most (2004).
e−Ei/T/Pje−Ej/T [6] H.S.ChanandK.A.Dill,J.Chem.Phys.100(12), 9238
likely a transient accumulating state, i.e., the other (1994).
conformations in its basin first rapidly flow to state i [7] M. Cieplak et al., Phys.Rev.Lett. 80, 3654 (1998).
before transiting to conformations outside the basin. If [8] S.B. Ozkan, K.A. Dill, and I. Bahar, Protein Sci. 11,
̺ < 1, depending on its initial probability of occupancy 1958 (2002).
[9] A. Sˇali, E.I. Shakhnovich,and M. Karplus, Nature 369,
and the final Boltzmann factor, state i may be either
248 (1994).
a monotonically decaying or accumulating state. We
[10] D.K. Klimov and D. Thirumalai, Phys. Rev. Lett. 76,
find that among the 493 local minima states for this
4070 (1995).
sequence, all except 3 are transiently accumulating, [11] R. Melin et al.,J. Chem. Phys. 110, 1252 (1999).
indicating they are responsible for forming transient [12] R.B.Sidje,ACMTrans.Math.Softw.24(1),130(1998).
state ensemble of various time scale. [13] I. Chang et al., Protein Sci. 13, 2446 (2004).