- Email: [email protected]
Kinetic consistency in the protein folding process1
Journal of Molecular Structure ŽTheochem. 461᎐462 Ž1999. 503᎐521
Kinetic consistency in the protein folding process 夽 Seiji Saito1, Masaki SasaiU Graduate School of Human Informatics, Nagoya Uni¨ ersity, Nagoya 464-8601, Japan Received 11 June 1998; accepted 7 July 1998
Abstract The folding process of a small helical protein is simulated with the knowledge-based potential model. Starting from a stretched random conformation, the model peptide chain rapidly collapses to the compact globule state within a 0.1- s time-scale. Three different types of behaviors are found in the molecular dynamics ŽMD. trajectories after compaction: folding trajectories along which the chain reaches the native conformation within a s time-scale; misfolding trajectories along which the chain itinerates globally different conformations from the native conformation; and escaping trajectories along which the chain once takes conformations fairly close to the native conformation but escapes to the other misfolded conformations with s scale. Consistent growth in both the short-range structure and the long-range structure is an important feature found in the folding trajectories. This kinetic consistency is not fulfilled in the escaping trajectories. A simple spin-glass-like model is developed and it is shown that the design of kinetic connectivity among conformations is important to make the yield of the folding trajectories large. 䊚 1999 Elsevier Science B.V. All rights reserved. Keywords: Protein folding; Knowledge-based potential; Spin glass; Energy landscape; Consistency principle
1. Introduction Recent theoretical w1᎐4x and experimental w5,6x study has shown that energy landscape is a key
夽 Dedicated to Professor Keiji Morokuma in celebration of his 65th birthday. U Corresponding author. E-mail: [email protected] 1 Present address: Information Science Lab. Multimedia Systems Laboratories, FUJITSU LABORATORIES LTD. 9-3, Nakase 1-chome, Mihama-ku, Chiba-shi Chiba 261-8588, Japan.
concept to describe the folding process of proteins. Due to the frustration among interactions, the energy landscape of proteins should be rugged with many minima corresponding to many different conformations w1᎐4x. In order to avoid being trapped into these glassy minima, the energy landscape has to be designed in some coherent fashion. From the results of lattice models w3,7᎐9x and spin-glass models w1,2,10x it was suggested that energy of fast folding proteins tends to decrease as the chain approaches its native conformation over the wide range of the conformational space. Such a globally biased energy landscape
0166-1280r99r$ - see front matter 䊚 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 6 - 1 2 8 0 Ž 9 8 . 0 0 4 4 0 - 0
504
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
was termed the folding funnel w7x. The folding process on the funnel-like energy surface is not adequately described by a single definite pathway but is more suitably described by a statistical ensemble of multiple pathways converging to the native conformation w1᎐3x. Different pathways could pass through different conformations at the transition state w1,2,9x. A new kinetic theory of the multidimensional process, which is capable of treating a variety of different pathways in a unified framework is required to be developed. By using off-lattice models w11᎐14x and spinglass models w15,16x, it has been suggested that the funnel is not structureless but has inhomogeneous features such as anisotropy w14,15x or hierarchy w12,13,16x. We should ask, therefore, how the characteristics of each protein are reflected in the funnel structure. Such a question is indeed important to interpret the experimental data on fast folding. On the refolding experiments of cytochrome c, for example, Sosnick et al. reported that the fast folding process of ms Žs 10y3 s. order and the slow folding process of 10 0 s order are coexisting in the same refolding condition w17x. This coexistence of multiple time-scales could be explained by the hierarchical structure of the funnel as shown schematically in Fig. 1. Starting from the extended state, the chain collapses to the compact globular state. When the chain collapses to conformations within the same basin of potential energy as the native conformation belongs to, then the chain would rapidly fold within a ms time-scale. When the chain arrives at the conformations in a different basin, then seconds or longer is needed to guide the chain to the basin of the native conformation. Folding landscape, therefore, should be described by multiple scales of energy ᎏ ones which characterize the energy difference between basins, and others which characterize the ruggedness within basins. The similar scheme on the coexistence of fast and slow folding processes was discussed by Guo and Thirumalai w11x and by Thirumalai w12x. Guo and Thirumalai simulated the folding process of a simplified model protein and showed that 60% of independent simulation runs reached the native conformation within a s Žs 10y6 s. time-scale
Fig. 1. A schematic representation of the trajectories on the hierarchical energy landscape. Trajectories bifurcate into the fast folding process with the probability and into the slow folding process with the probability 1 y .
and much longer time was needed for 40% of runs: the probability of the fast folding trajectories in their model was f 0.6. Thus, important questions are how is related to the funnel structure and how is determined by the sequence design. In the present paper we simulate the folding process of a model protein, calcium-binding protein wProtein Data Bank ŽPDB. code: 4icbx, and analyze how is determined. The X-ray conformation of 4icb is shown in Fig. 2. It has two loop regions where Ca2q ions are bound. For this size of peptide, molecular dynamics ŽMD. simulation with the full atomistic details is still too demanding for the present-day computer power, so that a bold simplification is necessary. Instead of using the standard molecular mechanics force field, we here use a knowledge-based potential introduced by Sasai w13x. This potential is able to fold the small helical proteins such as calcium-binding protein w13x, engrailed homeodomain w18x, and the
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
designed peptide w19x into their native conformations. Here, the chain is regarded to be in the native conformation when the global structure of the chain resembles to the X-ray structure on the designed structure in the resolution of a few angstroms. Though the precise configuration at the loop position in Fig. 2 should be essential to determine the Ca2q binding activity, we do not focus on such structural details. We show that the MD trajectories of this simplified model of folding are classified into three types: ‘folding trajectory’; ‘misfolding trajectory’; and ‘escaping trajectory’. It is shown that the consistent growth of both the short- and long-range structural order is an important feature of ‘folding trajectories’. One way to understand the complex behavior of the model is to construct a simpler model which mimics the behavior of the original model. Here, we introduce a simple spin-glass-like model to examine the mechanism determining in the knowledge-based potential model. It is shown that is much decreased when the probability of appearance of escaping trajectories is large.
505
2. Knowledge-based potential model 2.1. Model peptide chain
We consider a model polymer chain consisting of N s 75 amino-acid residues. A schematic configuration of the chain is shown in Fig. 3, where Ni , C ␣ i , CXi , and Oi are nuclei in the main chain at the ith site. -Carbon ŽC  . at the ith site is written as C i . Other atomistic details of the residues are neglected. Three-dimensional positions of Ni , C ␣ i , CXi , Oi , and C i are denoted by ª N ª ␣ ª CX ª O ri , ri , ri , ri , and ª ri  , respectively. It is assumed that nuclei around the peptide bond are kept on a single plane, so that the following constraints are imposed w20x: ªN ␣ O ri s 0.49rªiy1 q 0.72ª ri␣ y 0.21rªiy1 , ªC X ␣ ri s 0.45rªi␣ q 0.25rªiq1 q 0.30rªi O .
Ž1.
In MD calculations ª ri␣ , ª ri O , and ª ri are treated as X ªN independent variables and ri and ª ri C are determined by Eq. Ž1.. Potential energy V of the polymer has the form: V s V Ž main chain . q V Ž residue᎐residue . ,
Ž2.
where V Žmain chain. represents the steric constraints in the main chain. V Žresidue᎐residue. is the residue᎐residue interaction energy which will be explained in Section 2.2. The functional form of V Žmain chain. is similar
Fig. 2. X-ray structure of the calcium-binding protein Ž4icb.. Conformation of the main chain is represented by C ␣ s connected by lines.
Fig. 3. A schematic configuration of the model peptide chain. Side chain is represented only by C  . Hydrogen nuclei in the main chain are not explicitly taken into account.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
506
to the one used by Friedrichs et al. w20x. V Žmain chain. is a sum of four terms: V Ž main chain . s Vs q V q V q V .
Ž3.
Vs represents the constraints on distances between nuclei: N
Vs s ⌳ s
X
Ý ½ Ž ªriC yªri N
y 2.47 .
2
is1
2
2
X
qŽ ª ri yª ri N y 2.47 . q Ž ª ri C yª ri y 2.53 . qŽ
2
ª␣ ri yª ri
y 1.53 . q Ž
ªN ri yª ri␣
y 1.47 .
2
X
2
qŽ ª ri␣ yª ri C y 1.53 . q Ž ª ri␣ yª ri O y 2.40 . ␣ qŽ ª ri␣ yª riq1 y 3.81 .
qŽ
ªO ␣ ri yª riq1
2
2
y 2.80 .
2
N
Ý
Ž
5,
Ž4.
Ž region i . Ž region v and vi. Ž region ii, iii, and iv. , Ž7.
V Ž i . s
½
sin Ž i . sin Ž i y 30⬚.
Ž region i, ii, iii, and v. Ž region iv and vi. . Ž8.
Ž i y 0 . 2 ,
X ª␣ ri yª riC
. = ŽªriN yªri␣ . ⭈ Žªri␣ yªri . ,
Ž5.
where 0 is the equilibrium value, 0 s y2.50 ˚3. A V and V are potentials for dihedral angles. Dihedral angles i and ⌿i are defined around C ␣ i as shown in Fig. 3. We write N
V s ⌳
¡ysin Ž y 30⬚. q 0.5 ~sin Ž y 30⬚. q 0.5 ¢0 i
is1
i s
V Ž i . s
i
˚ .. where distances are expressed in angstroms ŽA V is the potential which represents the chiral constraint around C ␣ : V s ⌳
Fig. 4. The i y i plane is decomposed into six regions Žregion i᎐region vi. to define V Ž i . and V Ž i ..
Ý g iV Ž i . ,
V and V roughly reproduce the constraints exhibited in the Ramachandran plot w21x. Parameters ⌳ s , ⌳ , ⌳ , and ⌳ were chosen so that V Žmain chain. produces large energy penalties for the unrealistic distortions of the main chain. These parameter values, however, have to be made smaller than certain threshold values to avoid the numerical instability during the MD integration. We used ⌳ s s 0.2, ⌳ s 0.0001, ⌳ s 0.1, and ⌳ s 0.1. 2.2. Knowledge-based potential
is1 N
V s ⌳
Ý g iV Ž i . ,
Ž6.
is1
where g i s 0 when the ith site is glycine and g i s 1, otherwise. The i ᎐ i plane is decomposed into six regions Žregion i᎐region vi. as shown in Fig. 4 and V Ž i . and V Ž i . are defined in each region as:
We use V Žresidue᎐residue. similar to the one introduced by Sasai w13x. V Žresidue᎐residue. is a sum of two terms: V Ž residue᎐residue . s V Ž knowledge based . qV Ž hydrophobic attraction . . Ž9.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Table 1 PDB code of proteins used as the library to construct the knowledge-based potentials
V Žknowledge based. has the form: N
V Ž knowledge based . s
Ý Ý ⌳ k Vkp q Ž ri j . ,
Ž 10.
is1 j)i
where p and q distinguish 20 types of residues at the ith and jth site and ri j is distance between C  s, ri j s <ª ri yª r j <. For the short-range interactions of 0 - j y i F 10, k in Eq. Ž10. is distance along the sequence k s j y i. The long-range interactions are assumed to be insensitive to the precise value of j y i: interactions of 11 F j y i F 30 are represented by the potential of k s 11, interactions of 31 F j y i F 60 are represented by the potential of k s 12, and interactions of 61 F j y i are represented by the potential of k s 13. Vkp q Ž r . was constructed from a library of Nlibrary s 75 mutually dissimilar proteins. List of the PDB code of proteins used in the library is shown in Table 1, which is a subset of the library shown in Table 2 of Hendlich et al. w22x. Notice that the target protein, 4icb, or its homologue is not included in the library. When C  s of the residues p and q at the ith and Ž i q k .th site of protein p q, are separated at the distance ri,iqk , then a Gaussp q, ian function whose center is at r s ri,iqk was summed into the potential: Vkp q Ž r . s
y1 Nk Ž 2 Ck .
1r2
Nlibrary
=
Ý Ý exp
s1
y1 V11p q Ž r . s N11
2
q , . y Ž r y rip,iqk r2Ck ,
i
for k F 10,
Ý 11Fl-30
1 Ž 2 Cl . 1r2
Nlibrary
=
Ý Ý exp
s1
y1 V12p q Ž r . s N12
2
q , . y Ž r y rip,iql r2Cl ,
i
for 11 F l F 30,
Ý 31Fl-60
1 Ž 2 Cl . 1r2
Nlibrary
=
Ý Ý exp
s1
507
2
q , . y Ž r y rip,iql r2Cl ,
i
for 31 - l F 60,
155c 1cse e 1hmq a 1paz 1pyp 1tim a 2cab 2hhb a 2pka a 2ssi 3adk 3fxn 3rp2 a 4pti 7api a
1abp 1cts 1lh4 1pcy 1rei a 1wrp r 2cdv 2hhb b 2pka b 2stv 3b5c 3gap b 4ape 4sbv a 7api b
V13p q Ž r . s
y1 N13
1acx 1est 1lzl 1pfk a 1rhd 2alp 2cna 2lhb 2sbt 2taa a 3cyt o 3gpd g 4cpv 4tln 7wga a
1bds 1gcr 1mbd 1phn 1rn3 2aza a 2cyp 2lzm 2sga 2tbv a 3ebx 3grs 4dfr a 5adh 8cat a
1bp2 1gp1 a 1ovo a 1pp2 r 1sn3 2c2c 2fb4 l 2pab a 2sod 351c 3fxc 3pgk 4mdh a 5cpa 9pap
1 Ž 2 Cl . 1r2
Ý 61Fl
Nlibrary
=
Ý Ý exp
s1
2
q , . y Ž r y rip,iql r2Cl ,
i
for 61 F l F 100,
Ž 11.
where Nk is the number of residue pairs found in the library. Smoothness of the potential is determined by Ck . When Ck is too small, the potential has many small minima due to the sparseness Žsmall Nk . of the library data. We chose large enough Ck to eliminate these artificial minima. When Ck is too large, on the other hand, the potential becomes featureless and does not bring useful information to distinguish different se˚2 was used. quences. Here Ck s 0.5k A V Žknowledge based. is similar to the energy function developed by Sippl and colleagues to study the inverse folding problem on the sequence-conformation consistency w22,23x. Indeed, V Žknowledge based. is able to find the sequenceconformation consistency among library proteins: when 75 library sequences are threaded on the conformation which is arbitrarily chosen from the library, then V Žknowledge based. can discriminate the native sequence of the chosen conformation from other irrelevant sequences Ždata not shown..
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
508
It should be noted that V Žknowledge based. is constructed using only the data of compactly folded structures. In order to simulate the movement of the extended chain in a more reliable way, the data of high energy conformations different from the native one should be incorporated into the knowledge-based potentials. The method to take account of such data, however, has not yet been established. Here, we use a naive way to bypass this issue ᎏ V Žhydrophobic attraction . which has the similar Gaussian functional form to Vkp q Ž r . is added to express the long-range interactions of j y i ) 10: V Ž hydrophobic attraction . s ⌳ Ý i
y
Ý j)iq10
½
1 Ž 2 r2 .
1r2
exp yri2jr2 r2
g Ž p,q . 2 exp y Ž ri j y r . r2C 1r2 Ž 2 C .
5
, Ž 12.
mechanism because the energy difference between the native conformation and the misfolded conformation with incorrect stacking of helices is not significantly large. When R is as small as 0.2, on the other hand, the chain collapses at first and then helices begin to develop later. It is, however, difficult to find the conformation with well developed helical structures with this mechanism. In the present model the chain folds most efficiently when R has a moderate value; 0.3- R - 3. With this range of R, helix formation and collapse of the chain start simultaneously. Dependence of the folding kinetics on the ratio between the long- and short-range interactions has been discussed by several authors by using lattice models w25,26x. We will discuss this important issue in detail in the further publication. In the present paper, we focus on the folding mechanism under the condition of the balanced long- and short-range interactions: we here use ⌳ k s 0.2 for k F 10, ⌳ k s 0.5 for k ) 10, and ⌳ s 0.9. 2.3. Molecular dynamics simulation
where p and q are types of residues at the ith and jth site. g Ž p,q . s 1 when the ith or the jth site is hydrophobic and g Ž p,q . s 0, otherwise. For large ri j , V Žhydrophobic attraction . is a shallow attractive potential, so that C was chosen to ˚2 . In order to express be a large value, C s 100 A the repulsive interaction for small ri j , r was set ˚ to be 4 A. ⌳ k should have similar value to ⌳ since interactions of Eqs. Ž11. and Ž12. have the similar physical origins. Parameters ⌳ k and ⌳ should not be too large, ⌳ k , ⌳ F 1, in order to make V Žresidue᎐residue. softer potential than V Žmain chain.. Folding kinetics is sensitively dependent on the ratio between strengths of the short- and long-range interactions, R s Ž ⌳ k for k F 10.rŽ ⌳ k for k ) 10.; starting from the extended chain conformation, the folding process is much affected by the value of R. When R is as large as 5, helical structures are formed at first and the chain collapses after that. The folding process is much like the one described by the diffusion-collision model w24x. In the present model, however, the native conformation is not efficiently searched with this
We assume that the chain obeys the following Langevin equation of motion: ª dp ⭸V sy ªy␥ª pq Ž t . , dt ⭸r ª dª r p s , dt m
Ž 13.
where ª r is a 9N dimensional vector, ª rs ª␣ ªO ª ri ,ri ,ri 4 and ª p is the momentum vector. For the simplicity, mass of all atoms are assumed to ª Ž t . is a vector of random force, be equal; m s 1. ª ª ª Ž t . s i Ž t .4 , i s 1 y 9N. Ž t . is assumed to be Gaussian white noise and its variance is related to temperature T by the fluctuation-dissipation relation; ² i Ž t 1 . j Ž t 2 .: s m␥ T␦ i j ␦ Ž t 1 y t 2 . .
Ž 14.
A model chain under the consideration is assumed to have the sequence of calcium-binding protein Ž4cib.. In Eq. Ž13. sequence specific inter-
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
actions are taken into account through V Žresidue᎐residue. in V. We denote units of spatial length, mass, and energy in Eq. Ž13. by r 0 , m 0 , and 0 , respectively. ˚ These units can be roughly assigned as r 0 f 1 A, m 0 f 100 grmol, and 0 f 10 kJrmol. A typical length of time, t 0 is t 0 s r 02rD f 10y10 s,
Ž 15.
where D is the diffusion constant of small molecules in aqueous solution, D f 10y1 0 m2rs. ␥ can be estimated from the Stokes relation, ␥ s 6 a, by assuming to be the viscosity of water at room temperature, f 10y3 kgrmrs and af 10y1 0 m, which yields ␥ f 10y1 2 kgrs f 1 ⭈ m 0rt 0 . We used a slightly underdamped condition, ␥ s 0.6 but the qualitative results did not depend on the precise value of ␥ . Temperature T is measured in the unit of 0rk B . T was slowly decreased as the MD calculation proceeds. Temperature at the nth MD step was T s 0.053rwlogŽ0.0005n q 10.x2 0rk B . MD simulation was performed by numerically integrating Eqs. Ž13. and Ž14. with the discretizing width ⌬ t s 0.2. Initial conformation of the MD simulation was obtained by putting V Žresidue᎐residue. s 0 and stretching the N and C terminals of the random conformation in the opposite direction. Thirty independent MD runs were simulated starting from different initial conformations ª and using different random numbers Ž t .. The discretized equation of motion was integrated up to 10 5 steps. The total length of the simulation is, therefore, 10 5⌬ t ⭈ t 0 f 10y6 s s s. 2.4. Results Difference between the simulated conformation and the X-ray conformation is measured by
509
10 5 steps and 90% of MD runs ended up with ˚ In conformations of conformations of ␦ ) 3.5 A. ˚ ␦ - 2.8 A, four helices are stacking in the same way as in the X-ray structure. Conformations of ˚ on the other hand, are the misfolded ␦ ) 3.5 A, structures with the wrong way of stacking of four helices or with bent or distorted helices. Thus, in the present model calculation, the probability that the MD trajectories reach the native conformation within a s time-scale is f 0.1. These trajectories are termed ‘folding trajectories’. We found that the other trajectories leading to the misfolded conformations are classified into two types: ‘misfolding trajectories’ Ž60% of MD runs.; and ‘escaping trajectories’ Ž30% of MD runs.. A change of ␦ along each type of trajectory is shown in Fig. 5. Along the folding trajectory of Fig. 5a, ␦ gradually decreases to the small value. Along the misfolding trajectories of Fig. 5c, on ˚ the other hand, ␦ does not decrease below 3.5 A throughout the simulation. In the misfolding trajectories, more diverse behaviors of ␦ are found than in the folding trajectories, which is due to diversity of misfolded conformations. It is interesting to see in Fig. 5b that along the escaping trajectories ␦ decreases more rapidly than along the folding trajectories until 4 = 10 4 steps. After ˚ ␦ increases this rapid decrease to ␦ f 3.0 A, ˚ Along the escaping trajecagain to be ␦ ) 3.7 A. tory chain once takes conformations fairly close to the native conformation but these conformations are destroyed in the later stage of the process and the trajectory ‘escapes’ from the native conformation. Overall shape of the conformation is measured by asphericity A w27,28x: As
Ž R1 y R 2 . 2 q Ž R 2 y R 3 . 2 q Ž R 3 y R1 . 2 , 2 Ž R12 q R 22 q R 23 . Ž 17.
2 ␦s Ž N N y 1.
N
N
Ý Ý
ri␣j y riX-ray j
Ž 16.
is1 j)1
where ri␣j s <ª ri␣ y ª r j␣ < and riX-ray is distance j between C ␣ s of the ith and jth site in the X-ray conformation. Ten percent of independent MD ˚ within runs reached conformations of ␦ - 2.8 A
where R12 , R 22 , and R 23 are eigenvalues of the matrix
xx xy xz
xy yy zy
xz yz zz
0
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
510
Fig. 5. Difference between the simulated conformation and the X-ray conformation, ␦ observed along the folding trajectories Ža., along the escaping trajectories Žb., and along the misfolding trajectories Žc.. Two examples Žreal and dotted lines. of trajectories are shown for each case.
and each element of the matrix is defined by xxs
1 N
N
Ý Ž x i␣ . is1
2
,
xy s
1 N
N
Ý x i␣ yi␣ , is1
and so on. Here, x i␣ , yi␣ , and z i␣ are x, y, and z components of the vector ª ri␣. A f 1 when the chain has the stretched rod-like conformation and
A f 0 when the chain takes spherical globular conformation. Typical behaviors of A along the trajectories are shown in Fig. 6. In all kinds of trajectories A rapidly decreases from A f 1 to A f 0.03 within the first 10 4 steps and then fluctuates around A f 0.02: chain collapses rapidly at the early stage of the process within a 0.1- s time-scale and then chain searches the low energy
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
conformation among the compact globular conformations. Changes of A along the folding ŽFig. 6a. and escaping ŽFig. 6b. trajectories show the similar behaviors. There is larger diversity, however, in behaviors of A along the misfolding trajectory ŽFig. 6c.. The folded and misfolded conformations which appear at the late stage of processes have similar asphericity to each other.
511
Results shown in Figs. 5 and 6 are consistent with the picture described in Fig. 1. Type of the trajectory is determined just after the chain collapses to the globule. Thus, a specific collapse is necessary to guide the trajectory to the native conformation in a short time scale. Such explanation based on the difference between specific and non-specific collapses, however, is not enough to
Fig. 6. Asphericity, A, observed along the folding trajectories Ža., along the escaping trajectories Žb., and along the misfolding trajectories Žc.. Two examples Žreal and dotted lines. of trajectories are shown for each case.
512
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 7. A schematic representation of the MD trajectories on the hierarchical energy landscape. There are three types of trajectories: folding trajectories; misfolding trajectories; and escaping trajectories.
clarify the reason why the trajectory ‘escapes’ with a fairly high probability from the native
conformation. Taking account of the escaping trajectories, scheme of Fig. 1 is altered as in Fig. 7. In the following part of this section, differences between the folding and escaping trajectories are more closely investigated. Number of residues involved in helices, n h , is defined by the number of residues whose dihedral angles are in the region y80 - - y25⬚ and y90⬚ - - 0⬚. Change of n h along the trajectories is shown in Fig. 8. Both in folding and escaping trajectories, n h increases rapidly up to n h f 20 at the early stage of process: almost 40% of helices are formed before chain completes to collapse. After compaction, n h increases steadily as the chain searches the low energy conformation; n h in the escaping trajectories, however, does not grow as large as in the folding trajectories: growth of helices is limited to the level of 80% of the native conformation in the escaping trajectories. We denote four helices in the native conformation by helix-I with I s 1᎐4 counting from the N terminal side to the C terminal side. When helix-I is composed of the residues i,i q 1,...,i q a and helix-J is composed of the residues j, j q 1,..., j q b, then distance between helix-I and helix-J is
Fig. 8. Helicity, n h , observed along the folding trajectories Ža. and along the escaping trajectories Žb.. A horizontal line shows the value in the X-ray conformation. Two examples Žreal and dotted lines. of trajectories are shown for each case.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
defined by hI J s
1 aq 1
iqa
Ý min Ž rk␣j ,rk␣jq1 ,...,rk␣jqb . ,
Ž 18.
ksi
where min Ž.... is the function to select the minimum value. When helix-I is placed in the direction at angle I J to the direction of helix-J, CI J s cos I J is defined by the inner product of two unit vectors; CI J s
ž ⭈
2 aŽ aq 1 .
ž
iqay1
2 Ž b bq 1 .
Ý ksi
ª␣ r k yª r l␣ ª␣ ª␣ lskq1 < r k y r l <
jqby1
Ý X
k sj
iqa
Ý
/
ª␣ r k X yª r l␣X ª␣ ª␣ X X l sk q1 < r k X y r l X <
jqb
Ý
/
.
Ž 20.
These definitions of h I J and CI J have clear meanings only in the conformations which have well developed helices. h I J and CI J monitored along the MD trajectory, however, should give information of the growth process of tertiary structure. Changes of h13 and h 23 are shown in Fig. 9 and changes of C14 and C34 are shown in Fig. 10. Results shown in Figs. 9 and 10 are consistent with the results shown in Fig. 5. Along the folding trajectories the tertiary structure is gradually formed during the structural search among the compact globular conformations. Along the escaping trajectories, on the other hand, the tertiary structure develops rather rapidly and the conformations which have features common to the native conformation are formed at the early stage of the process. These features, however, are destroyed as the MD steps proceed and the chain reaches conformations which have globally different structures from the native one. Differences between the folding and escaping trajectories suggest that consistency between the secondary structure formation and the tertiary structure formation is a necessary condition to guide the trajectories to the native conformation. When the strengths of the short- and long range interactions are balanced, approximately half of the secondary structures are formed when the chain collapses to the globule. Then, both the secondary and tertiary structures start to develop
513
cooperatively in the compact globular state. When the tertiary structural order is formed much earlier than the secondary structure, then the later growth of the secondary structure can bring the chain to the different tertiary conformation. To describe this escaping process, we could use an analogy of solving a puzzle like Rubik’s cube. When one side of the cube is arranged to have the correct order, then the later arrangement of the other sides often destroys the side which was arranged in a earlier step. In order to arrange multiple sides of the cube, orders of different sides have to be developed in a synergetic and consistent way. Present model calculations suggest that competition and cooperation of different structural formations require the consistent growth of different structural orders for the chain to fold fast with s scale. We should call this principle for fast folding ‘the consistency principle on the kinetic process’, which is a reminiscent of the consistency principle first discussed by Go based on the analyses of the 1st order nature of the folding᎐unfolding transition and on the analyses of the stability of the native conformation w29x. 3. Spin-glass-like model 3.1. Hierarchical dynamics In this section we examine the idea of the consistency principle on the kinetic process by constructing a simple spin-glass model. Two kinds of spin variables, Si for i s 1 y N and TI for I s 1 y M, are introduced to represent the coarse-grained secondary and tertiary structures of the protein. We assume S i and TI take values of "1. Si represents the secondary structure: when the helical structure is formed at the ith residue, then Si s 1. Si s y1, otherwise. TI represents the direction of helix-I. When TI s 1, helix-I is placed in the same direction as in the native conformation. When TI s y1, helix-I is in a different direction. Each helix is assumed to be composed of NrMs K residues: residues i s 1,...,K are the constituents of helix-1, residues i s Kq 1,...,2 K are the constituents of helix-2, and so on. We write i g I if i belongs to the
514
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 9. Distance between helix-1 and helix-3, h13 , observed along the folding trajectories Ža. and along the escaping trajectories Žb., and distance between helix-2 and helix-3, h 23 , observed along the folding trajectories Žc. and along the escaping trajectories Žd.. A horizontal line shows the value in the X-ray conformation. Two examples Žreal and dotted lines. of trajectories are shown for each case.
region of the sequence where helix-I is formed: residues i g 1 for i s 1,...,K, i g 2 for i s Kq 1,...,2 K, and so on Žsee Fig. 11.. In the native conformation Si s TI s 1 for all i and I. Energy H of the system is
Ny1
Hsy
Ý
Ž rJiiq1 q 1 . Si Siq1
is1
y
1 R st
My1
M
Ý Ý Ž r t JI J q 1. TI TJ . Is1 J)1
Ž 21.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
515
Fig. 10. Inner product of unit vector parallel to helix-1 and the one parallel to helix-4, C14 , observed along the folding trajectories Ža. and along the escaping trajectories Žb., and inner product of unit vector parallel to helix-3 and the one parallel to helix-4, C34 , observed along the folding trajectories Žc. and along the escaping trajectories Žd.. A horizontal line shows the value in the X-ray conformation. Two examples Žreal and dotted lines. of trajectories are shown for each case.
Here, the first term is energy of formation of helices, which is a sum of the short-range interactions. The second term is the helix᎐helix interactions which are the long-range interactions. Ji iq1 is Gaussian random number with mean 0 and
standard deviation 1r62. JI J is Gaussian random number with mean 0 and standard deviation 1r 'My 1 . r and r t are parameters which determine the extent of randomness. When r or r t is larger, interactions are more frustrated and
516
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 11. Si represents whether ith site is helical Ž Si s 1. or not Ž Si s y1. and TI represents whether helix-I is aligned in the same direction as in the native conformation ŽTI s 1. or not ŽTI s y1.. A case for N s 16 and Ms 4 is shown.
energy landscape is more rugged. R st is a parameter which determines the balance between the short-range interactions and the long-range interactions, which resembles to the parameter R in V Žresidue᎐residue. of Section 2. For most of the calculations we used r s 0.3, r t s 0.3, and R st s 1.67. N and M were set to be N s 30 and Ms 10 or N s 16 and Ms 4. Si and TI are changed with the Monte Carlo ŽMC. method. We assume that Si moves much faster than TI and that motion of TI is slaving to the motion of Si : to generate a trial configuration to update the MC step, sign of Si is flipped, at first, at the randomly chosen position i. When i g I, sign of TI is altered with the probability 0 - ␣ - 1 and is kept unchanged with the probability 1 y ␣ . Trial configuration prepared in this way is judged to be accepted or rejected by the Metropolis algorithm with temperature T. This MC motion is schematically shown in Fig. 12. It should be noted that S i and TI are uncoupled in the energy expression of Eq. Ž21., so that the equilibrium properties of the system are described just in terms of two independent spin systems. Si and TI , however, are kinetically strongly coupled to each other, so that when the MC trajectory is truncated with the relatively small number of steps, then the probability of the folding trajectory, , should be decisively affected by ␣ which controls the kinetic coupling between Si and TI . We here use this spin model to focus on such kinetic effects on the fast folding system. T dependence of is shown in Fig. 13. Starting from a random configuration, MC trajectories of 3000 steps were performed. Temperature T was kept constant through 3000 steps. From Fig.
Fig. 12. A schematic representation of folding trajectory in the spin-glass like system. When Si is flipped, TI is altered with the probability ␣ . Vetical columns are energy spectra for fixes configuration of TI . Arrows of real lines are folding trajectories and arrows of dotted lines are escaping trajectories. Both the secondary and tertiary structures have to develop in a cooperative way to reach the native conformation.
13, we can see that folding is difficult both in the high and low temperature regions. For T ) 0.5 random conformations are favored by the strong thermal fluctuation. For T - 0.2 system is trapped into one of many random minima in energy. Thus, folding is most efficient in the temperature range 0.2- T - 0.5. In the rest of this section, in order to compare the results to those of Section 2, T is gradually decreased as MC steps proceed: T at the nth step is T s T Ž0. y w T Ž0. y T Ž1.x nr15 000 for n - 15 000 and T s T Ž1. for n G 15 000. We used T Ž0. s 0.5 and T Ž1. s 0.2. 3.2. Results Distance between the simulated conformation and the native conformation is measured by d: ds
1 N
N
1
M
Ý Ž1 y Si . r4 q M Ý Ž1 y TI . r4,
is1
Ž 22.
is1
which resembles ␦ in Section 2. Difference
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 13. Dependence of T on the probability of the folding trajectories . Results of 2000 independent MC runs are averaged. N s 16, Ms 4 and ␣ s 0.2.
between the secondary structure of the simulated conformation and the one of the native conformation is measured by
dS s
1 N
N
Ý Ž1 y Si . r2.
Ž 23.
is1
N Ž1 y d S . corresponds to n h of Section 2. In Fig. 14 typical examples of H, d, and d S observed along the folding trajectory are shown. Features of the trajectory are similar to those of the folding trajectories shown in Figs. 6a and 9a. Both the short- and long-range structures develop steadily in a cooperative way. Trajectories which did not reach the native conformation within 3000 steps were further elongated to 30 000 steps. Most of these trajectories were escaping trajectories for this parameterization. Typical example of escaping trajectory is shown in Fig. 15. Features of the trajectory is similar to those of the escaping trajectories shown in Fig. 6b and Fig. 9b. When the long-range structural order is formed at too early stage, then the growth of short-range structure in the later stage destabilizes the early formed struc-
517
tural order. Thus, this simple spin-glass-like model is able to reproduce the distinction between the folding trajectories and the escaping trajectories found in the knowledge-based potential model. is a decreasing function of r and r t : s 1 when r s r t s 0 and f 0 when r s r t ) 1.0. For large r and r t , the system is frozen into the glassy random conformation and does not reach the native conformation. For small r and r t the ratio ␦ Er⌬ E should be large, where ␦ E is energy gap between energy of the native conformation and the average energy of other conformations, and ⌬ E is the average ruggedness in the energy landscape w1,2x. This design of r and r t for fast folding is the requirement that the funnel should be well developed around the native conformation. Such a design principle was called the principle of minimal frustration w1,2,10x. When R st s 5, secondary structures develop quickly but the long-range order does not grow within 30 000 steps. When R st s 1.25, on the other hand, the long-range structures develop quickly but it takes as long as 30 000 steps for secondary structures to develop. Folding is most efficient in the intermediate range, 1.5- R st - 2: balance between strengths of the short- and long-range interactions is necessary to make large. This design principle was called the consistency principle w29x. Alternation of R st changes the energy spectrum and thus R st affects ␦ Er⌬ E: the consistency principle has a great overlap on its content with the principle of minimal frustration. R st , however, affects not only the first moment, ␦ E, and the second moment, ⌬ E, of the energy spectrum but also the higher moments of the spectrum. Therefore, we should note that the consistency principle refers not only to the existence of the funnel but also to the rich phenomena relating to the inner structure of the funnel. Dependence of on ␣ is shown in Fig. 16. is largest in the range 0.2- ␣ - 0.4. For ␣ - 0.1, misfolding trajectories are dominant and for ␣ ) 0.5, escaping trajectories are dominant. In order to make large, ␣ has to be tuned to assure the consistent growth of both the short- and longrange structural orders. This design principle is
518
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 14. Energy, H, distance to the native conformation, d, and difference in the secondary structures from the native conformation, d S, observed along an example of the folding trajectory. N s 30, Ms 10 and ␣ s 0.2.
the consistency principle on the kinetic process discussed in the previous section. By changing ␣ , ␦ Er⌬ E is not altered but the kinetic connectivity among conformations is altered. Thus, the toy model in this section suggests that design on the kinetic consistency, that is, design of the inner structure of the funnel is important to design the fast folding sequences. Effects of the connectivity
among conformations on the folding speed were discussed by Betancourt and Onuchic w30x by using the lattice model. In the present paper we stress the idea that characteristics of the hierarchical structures which have not been taken account of in the lattice model should be an important factor to design the connectivity among conformations.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
519
Fig. 15. Energy, H, distance to the native conformation, d, and difference in the secondary structures from the native conformation, d S, observed along an example of the escaping trajectory. N s 30, Ms 10 and ␣ s 0.2.
4. Discussion
It was shown from the results of the knowledge-based potential model and a spinglass-like model that consistent and cooperative growth of the short-range and the long-range
structural orders is a necessary condition for the trajectory to reach the native conformation in a short time. Therefore, to design the fast folding peptides, sequence has to be selected to assure the kinetic consistency among different structural orders. It should be interesting to ask whether the
520
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 16. Dependence of kinetic parameter ␣ on the probability of the folding trajectories . Results of 2000 independent MC runs are averaged. N s 16, Ms 4.
kinetic consistency is always in harmony with the stabilization of the native conformation or these two requirements sometimes contradict each other. This point might be relevant to the difficulty encountered in de novo design of peptides w31x. For most of cases, de novo design of peptides is based on the stability condition of the native conformation. It should not be clear whether such design principle gives an optimal sequence from a view point of the kinetic consistency. The same question might be relevant to the reason why natural proteins are marginal in their stability w32x. Free energy difference between the unfolded and the folded states are, for most of small globular proteins, less than 10᎐20 kcalrmol. This value is quite small considering the fact that small proteins consist of more than several hundreds atoms. The native conformation of natural protein is, therefore, not thoroughly stabilized but is marginally stable. Kinetic requirements could be reasons for this marginal stability. Thus, further investigations on the relation between stability and kinetic consistency should be important. There is much room to improve the knowledge-based potential model used in this paper. Especially important is to incorporate the
data of high energy states into potentials. As shown in the present paper the way of compaction is a large factor to determine . Thus, potentials must have ability to reproduce the compaction process in a proper way. Statistical data of extended conformations and compact globular conformations other than the native conformation should be used to construct such potentials. It should be important to check whether consequences of the present paper, such as existence of three types of trajectories and necessity of kinetic consistency in the folding trajectories should hold in the improved models of folding. Important insights will be obtained from experimental data on the growth of the secondary structure and on the tertiary structure measured at the same time in s time scale. Such data can provide information on how the kinetic consistency is realized in natural proteins and how the degree of consistency can be altered by the protein engineering technique. Through the comparison between experiments and simulations, more precise understanding of the structure of the funnel will be obtained, which should provide deep insights on how the fast folding proteins have been designed in the evolutionary history. References w1x J.D. Bryngelson, J.N. Onuchic, N.D. Socci, P.G. Wolynes, Proteins 21 Ž1995. 167. w2x J.N. Onucic, Z. Luthey-Schulten, P.G. Wolynes, Annu. Rev. Phys. Chem. 48 Ž1997. 539. w3x K.A. Dill, H.S. Chan, Nature Struct. Biol. 4 Ž1997. 10. w4x E.I. Shakhnovich, Curr. Opin. Struct. Biol. 7 Ž1997. 29. w5x W. Eaton, P.A. Thompson, C.-K. Chan, S.J. Hagen, J. Hofrichter, Structure 4 Ž1996. 1133. w6x W. Eaton, et al., Curr. Opin. Struct. Biol. 7 Ž1997. 10. w7x P.E. Leopold, M. Montal, J.N. Onuchic, Proc. Natl. Acad. Sci. USA. 89 Ž1992. 8721. w8x C.J. Camacho, D. Thirumalai, Proc. Natl. Acad. Sci. U.S.A. 90 Ž1993. 6369. w9x A. Sali, E.I. Shakhnovich, M. Karplus, J. Mol. Biol. 235 Ž1994. 1614. w10x J.D. Bryngelson, P.G. Wolynes, Proc. Natl. Acad. Sci. U.S.A. 84 Ž1997. 7524. w11x Z. Guo, D. Thirumalai, Biopolymers 36 Ž1995. 83. w12x D. Thirumalai, J. Phys. I ŽParis. 5 Ž1995. 1457. w13x M. Sasai, Proc. Natl. Acad. Sci. USA. 92 Ž1995. 8438. w14x B.A. Shoemaker, J. Wang, P.G. Wolynes, Proc. Natl. Acad. Sci. USA. 94 Ž1997. 777.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521 w15x S. Saito, M. Sasai, T. Yomo, Proc. Natl. Acad. Sci. USA. 94 Ž1997. 11324. w16x S.S. Plotkin, J. Wang, P.G. Wolynes, J. Chem. Phys. 106 Ž1997. 2932. w17x T.R. Sosnick, L. Mayne, S.W. Englander, Proteins 24 Ž1996. 413. w18x T. Yomo, S. Saito, M. Sasai Žto be published.. w19x S. Kamtekar, J.M. Sciffer, H. Xiong, J.M. Babik, M.H. Hecht, Science 262 Ž1993. 1680. w20x M.S. Friedrichs, R.A. Goldstein, P.G. Wolynes, J. Mol. Biol. 222 Ž1991. 1013. w21x J.S. Richardson, Adv. Protein Chem. 34 Ž1981. 168. w22x M. Hendlich, P. Lackner, S. Weitckus, H. Floeckner, R. Froschauer, K. Gottsbacher, G. Casari, M.J. Sippl, J. Mol. Biol. 216 Ž1992. 167. w23x M.J. Sippl, J. Mol. Biol. 213 Ž1992. 859.
521
w24x M. Karplus, D.L. Weaver, Nature 260 Ž1976. 404. w25x N. Go, H. Taketomi, Proc. Natl. Acad. Sci. USA. 75 Ž1979. 559. w26x A. Rey, J. Skolnick, Proteins 16 Ž1993. 8. w27x J. Rudnick, G. Gaspari, J. Phys. A 19 Ž1986. L191. w28x O.F. Olaj, G. Zifferer, J. Chem. Phys. 100 Ž1994. 636. w29x N. Go, Ann. Rev. Biophys. Bioeng. 12 Ž1983. 183. w30x M.R. Betancourt, J.N. Onuchic, J. Chem. Phys. 103 Ž1995. 773. w31x K.A. Dill, S. Bromberg, K. Yue, K.M. Fiebig, D.P. Yee, P.D. Thomas, H.S. Chan, Protein Sci. 4 Ž1995. 561. w32x P.G. Wolynes, in: H. Bohr, S. Brunak ŽEds.., Proceedings on Symposium on Distance-Based Approaches to Protein Structure Determination II, CRC Press, Boca Raton, 1994, p. 3᎐17.
Kinetic consistency in the protein folding process 夽 Seiji Saito1, Masaki SasaiU Graduate School of Human Informatics, Nagoya Uni¨ ersity, Nagoya 464-8601, Japan Received 11 June 1998; accepted 7 July 1998
Abstract The folding process of a small helical protein is simulated with the knowledge-based potential model. Starting from a stretched random conformation, the model peptide chain rapidly collapses to the compact globule state within a 0.1- s time-scale. Three different types of behaviors are found in the molecular dynamics ŽMD. trajectories after compaction: folding trajectories along which the chain reaches the native conformation within a s time-scale; misfolding trajectories along which the chain itinerates globally different conformations from the native conformation; and escaping trajectories along which the chain once takes conformations fairly close to the native conformation but escapes to the other misfolded conformations with s scale. Consistent growth in both the short-range structure and the long-range structure is an important feature found in the folding trajectories. This kinetic consistency is not fulfilled in the escaping trajectories. A simple spin-glass-like model is developed and it is shown that the design of kinetic connectivity among conformations is important to make the yield of the folding trajectories large. 䊚 1999 Elsevier Science B.V. All rights reserved. Keywords: Protein folding; Knowledge-based potential; Spin glass; Energy landscape; Consistency principle
1. Introduction Recent theoretical w1᎐4x and experimental w5,6x study has shown that energy landscape is a key
夽 Dedicated to Professor Keiji Morokuma in celebration of his 65th birthday. U Corresponding author. E-mail: [email protected] 1 Present address: Information Science Lab. Multimedia Systems Laboratories, FUJITSU LABORATORIES LTD. 9-3, Nakase 1-chome, Mihama-ku, Chiba-shi Chiba 261-8588, Japan.
concept to describe the folding process of proteins. Due to the frustration among interactions, the energy landscape of proteins should be rugged with many minima corresponding to many different conformations w1᎐4x. In order to avoid being trapped into these glassy minima, the energy landscape has to be designed in some coherent fashion. From the results of lattice models w3,7᎐9x and spin-glass models w1,2,10x it was suggested that energy of fast folding proteins tends to decrease as the chain approaches its native conformation over the wide range of the conformational space. Such a globally biased energy landscape
0166-1280r99r$ - see front matter 䊚 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 6 - 1 2 8 0 Ž 9 8 . 0 0 4 4 0 - 0
504
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
was termed the folding funnel w7x. The folding process on the funnel-like energy surface is not adequately described by a single definite pathway but is more suitably described by a statistical ensemble of multiple pathways converging to the native conformation w1᎐3x. Different pathways could pass through different conformations at the transition state w1,2,9x. A new kinetic theory of the multidimensional process, which is capable of treating a variety of different pathways in a unified framework is required to be developed. By using off-lattice models w11᎐14x and spinglass models w15,16x, it has been suggested that the funnel is not structureless but has inhomogeneous features such as anisotropy w14,15x or hierarchy w12,13,16x. We should ask, therefore, how the characteristics of each protein are reflected in the funnel structure. Such a question is indeed important to interpret the experimental data on fast folding. On the refolding experiments of cytochrome c, for example, Sosnick et al. reported that the fast folding process of ms Žs 10y3 s. order and the slow folding process of 10 0 s order are coexisting in the same refolding condition w17x. This coexistence of multiple time-scales could be explained by the hierarchical structure of the funnel as shown schematically in Fig. 1. Starting from the extended state, the chain collapses to the compact globular state. When the chain collapses to conformations within the same basin of potential energy as the native conformation belongs to, then the chain would rapidly fold within a ms time-scale. When the chain arrives at the conformations in a different basin, then seconds or longer is needed to guide the chain to the basin of the native conformation. Folding landscape, therefore, should be described by multiple scales of energy ᎏ ones which characterize the energy difference between basins, and others which characterize the ruggedness within basins. The similar scheme on the coexistence of fast and slow folding processes was discussed by Guo and Thirumalai w11x and by Thirumalai w12x. Guo and Thirumalai simulated the folding process of a simplified model protein and showed that 60% of independent simulation runs reached the native conformation within a s Žs 10y6 s. time-scale
Fig. 1. A schematic representation of the trajectories on the hierarchical energy landscape. Trajectories bifurcate into the fast folding process with the probability and into the slow folding process with the probability 1 y .
and much longer time was needed for 40% of runs: the probability of the fast folding trajectories in their model was f 0.6. Thus, important questions are how is related to the funnel structure and how is determined by the sequence design. In the present paper we simulate the folding process of a model protein, calcium-binding protein wProtein Data Bank ŽPDB. code: 4icbx, and analyze how is determined. The X-ray conformation of 4icb is shown in Fig. 2. It has two loop regions where Ca2q ions are bound. For this size of peptide, molecular dynamics ŽMD. simulation with the full atomistic details is still too demanding for the present-day computer power, so that a bold simplification is necessary. Instead of using the standard molecular mechanics force field, we here use a knowledge-based potential introduced by Sasai w13x. This potential is able to fold the small helical proteins such as calcium-binding protein w13x, engrailed homeodomain w18x, and the
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
designed peptide w19x into their native conformations. Here, the chain is regarded to be in the native conformation when the global structure of the chain resembles to the X-ray structure on the designed structure in the resolution of a few angstroms. Though the precise configuration at the loop position in Fig. 2 should be essential to determine the Ca2q binding activity, we do not focus on such structural details. We show that the MD trajectories of this simplified model of folding are classified into three types: ‘folding trajectory’; ‘misfolding trajectory’; and ‘escaping trajectory’. It is shown that the consistent growth of both the short- and long-range structural order is an important feature of ‘folding trajectories’. One way to understand the complex behavior of the model is to construct a simpler model which mimics the behavior of the original model. Here, we introduce a simple spin-glass-like model to examine the mechanism determining in the knowledge-based potential model. It is shown that is much decreased when the probability of appearance of escaping trajectories is large.
505
2. Knowledge-based potential model 2.1. Model peptide chain
We consider a model polymer chain consisting of N s 75 amino-acid residues. A schematic configuration of the chain is shown in Fig. 3, where Ni , C ␣ i , CXi , and Oi are nuclei in the main chain at the ith site. -Carbon ŽC  . at the ith site is written as C i . Other atomistic details of the residues are neglected. Three-dimensional positions of Ni , C ␣ i , CXi , Oi , and C i are denoted by ª N ª ␣ ª CX ª O ri , ri , ri , ri , and ª ri  , respectively. It is assumed that nuclei around the peptide bond are kept on a single plane, so that the following constraints are imposed w20x: ªN ␣ O ri s 0.49rªiy1 q 0.72ª ri␣ y 0.21rªiy1 , ªC X ␣ ri s 0.45rªi␣ q 0.25rªiq1 q 0.30rªi O .
Ž1.
In MD calculations ª ri␣ , ª ri O , and ª ri are treated as X ªN independent variables and ri and ª ri C are determined by Eq. Ž1.. Potential energy V of the polymer has the form: V s V Ž main chain . q V Ž residue᎐residue . ,
Ž2.
where V Žmain chain. represents the steric constraints in the main chain. V Žresidue᎐residue. is the residue᎐residue interaction energy which will be explained in Section 2.2. The functional form of V Žmain chain. is similar
Fig. 2. X-ray structure of the calcium-binding protein Ž4icb.. Conformation of the main chain is represented by C ␣ s connected by lines.
Fig. 3. A schematic configuration of the model peptide chain. Side chain is represented only by C  . Hydrogen nuclei in the main chain are not explicitly taken into account.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
506
to the one used by Friedrichs et al. w20x. V Žmain chain. is a sum of four terms: V Ž main chain . s Vs q V q V q V .
Ž3.
Vs represents the constraints on distances between nuclei: N
Vs s ⌳ s
X
Ý ½ Ž ªriC yªri N
y 2.47 .
2
is1
2
2
X
qŽ ª ri yª ri N y 2.47 . q Ž ª ri C yª ri y 2.53 . qŽ
2
ª␣ ri yª ri
y 1.53 . q Ž
ªN ri yª ri␣
y 1.47 .
2
X
2
qŽ ª ri␣ yª ri C y 1.53 . q Ž ª ri␣ yª ri O y 2.40 . ␣ qŽ ª ri␣ yª riq1 y 3.81 .
qŽ
ªO ␣ ri yª riq1
2
2
y 2.80 .
2
N
Ý
Ž
5,
Ž4.
Ž region i . Ž region v and vi. Ž region ii, iii, and iv. , Ž7.
V Ž i . s
½
sin Ž i . sin Ž i y 30⬚.
Ž region i, ii, iii, and v. Ž region iv and vi. . Ž8.
Ž i y 0 . 2 ,
X ª␣ ri yª riC
. = ŽªriN yªri␣ . ⭈ Žªri␣ yªri . ,
Ž5.
where 0 is the equilibrium value, 0 s y2.50 ˚3. A V and V are potentials for dihedral angles. Dihedral angles i and ⌿i are defined around C ␣ i as shown in Fig. 3. We write N
V s ⌳
¡ysin Ž y 30⬚. q 0.5 ~sin Ž y 30⬚. q 0.5 ¢0 i
is1
i s
V Ž i . s
i
˚ .. where distances are expressed in angstroms ŽA V is the potential which represents the chiral constraint around C ␣ : V s ⌳
Fig. 4. The i y i plane is decomposed into six regions Žregion i᎐region vi. to define V Ž i . and V Ž i ..
Ý g iV Ž i . ,
V and V roughly reproduce the constraints exhibited in the Ramachandran plot w21x. Parameters ⌳ s , ⌳ , ⌳ , and ⌳ were chosen so that V Žmain chain. produces large energy penalties for the unrealistic distortions of the main chain. These parameter values, however, have to be made smaller than certain threshold values to avoid the numerical instability during the MD integration. We used ⌳ s s 0.2, ⌳ s 0.0001, ⌳ s 0.1, and ⌳ s 0.1. 2.2. Knowledge-based potential
is1 N
V s ⌳
Ý g iV Ž i . ,
Ž6.
is1
where g i s 0 when the ith site is glycine and g i s 1, otherwise. The i ᎐ i plane is decomposed into six regions Žregion i᎐region vi. as shown in Fig. 4 and V Ž i . and V Ž i . are defined in each region as:
We use V Žresidue᎐residue. similar to the one introduced by Sasai w13x. V Žresidue᎐residue. is a sum of two terms: V Ž residue᎐residue . s V Ž knowledge based . qV Ž hydrophobic attraction . . Ž9.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Table 1 PDB code of proteins used as the library to construct the knowledge-based potentials
V Žknowledge based. has the form: N
V Ž knowledge based . s
Ý Ý ⌳ k Vkp q Ž ri j . ,
Ž 10.
is1 j)i
where p and q distinguish 20 types of residues at the ith and jth site and ri j is distance between C  s, ri j s <ª ri yª r j <. For the short-range interactions of 0 - j y i F 10, k in Eq. Ž10. is distance along the sequence k s j y i. The long-range interactions are assumed to be insensitive to the precise value of j y i: interactions of 11 F j y i F 30 are represented by the potential of k s 11, interactions of 31 F j y i F 60 are represented by the potential of k s 12, and interactions of 61 F j y i are represented by the potential of k s 13. Vkp q Ž r . was constructed from a library of Nlibrary s 75 mutually dissimilar proteins. List of the PDB code of proteins used in the library is shown in Table 1, which is a subset of the library shown in Table 2 of Hendlich et al. w22x. Notice that the target protein, 4icb, or its homologue is not included in the library. When C  s of the residues p and q at the ith and Ž i q k .th site of protein p q, are separated at the distance ri,iqk , then a Gaussp q, ian function whose center is at r s ri,iqk was summed into the potential: Vkp q Ž r . s
y1 Nk Ž 2 Ck .
1r2
Nlibrary
=
Ý Ý exp
s1
y1 V11p q Ž r . s N11
2
q , . y Ž r y rip,iqk r2Ck ,
i
for k F 10,
Ý 11Fl-30
1 Ž 2 Cl . 1r2
Nlibrary
=
Ý Ý exp
s1
y1 V12p q Ž r . s N12
2
q , . y Ž r y rip,iql r2Cl ,
i
for 11 F l F 30,
Ý 31Fl-60
1 Ž 2 Cl . 1r2
Nlibrary
=
Ý Ý exp
s1
507
2
q , . y Ž r y rip,iql r2Cl ,
i
for 31 - l F 60,
155c 1cse e 1hmq a 1paz 1pyp 1tim a 2cab 2hhb a 2pka a 2ssi 3adk 3fxn 3rp2 a 4pti 7api a
1abp 1cts 1lh4 1pcy 1rei a 1wrp r 2cdv 2hhb b 2pka b 2stv 3b5c 3gap b 4ape 4sbv a 7api b
V13p q Ž r . s
y1 N13
1acx 1est 1lzl 1pfk a 1rhd 2alp 2cna 2lhb 2sbt 2taa a 3cyt o 3gpd g 4cpv 4tln 7wga a
1bds 1gcr 1mbd 1phn 1rn3 2aza a 2cyp 2lzm 2sga 2tbv a 3ebx 3grs 4dfr a 5adh 8cat a
1bp2 1gp1 a 1ovo a 1pp2 r 1sn3 2c2c 2fb4 l 2pab a 2sod 351c 3fxc 3pgk 4mdh a 5cpa 9pap
1 Ž 2 Cl . 1r2
Ý 61Fl
Nlibrary
=
Ý Ý exp
s1
2
q , . y Ž r y rip,iql r2Cl ,
i
for 61 F l F 100,
Ž 11.
where Nk is the number of residue pairs found in the library. Smoothness of the potential is determined by Ck . When Ck is too small, the potential has many small minima due to the sparseness Žsmall Nk . of the library data. We chose large enough Ck to eliminate these artificial minima. When Ck is too large, on the other hand, the potential becomes featureless and does not bring useful information to distinguish different se˚2 was used. quences. Here Ck s 0.5k A V Žknowledge based. is similar to the energy function developed by Sippl and colleagues to study the inverse folding problem on the sequence-conformation consistency w22,23x. Indeed, V Žknowledge based. is able to find the sequenceconformation consistency among library proteins: when 75 library sequences are threaded on the conformation which is arbitrarily chosen from the library, then V Žknowledge based. can discriminate the native sequence of the chosen conformation from other irrelevant sequences Ždata not shown..
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
508
It should be noted that V Žknowledge based. is constructed using only the data of compactly folded structures. In order to simulate the movement of the extended chain in a more reliable way, the data of high energy conformations different from the native one should be incorporated into the knowledge-based potentials. The method to take account of such data, however, has not yet been established. Here, we use a naive way to bypass this issue ᎏ V Žhydrophobic attraction . which has the similar Gaussian functional form to Vkp q Ž r . is added to express the long-range interactions of j y i ) 10: V Ž hydrophobic attraction . s ⌳ Ý i
y
Ý j)iq10
½
1 Ž 2 r2 .
1r2
exp yri2jr2 r2
g Ž p,q . 2 exp y Ž ri j y r . r2C 1r2 Ž 2 C .
5
, Ž 12.
mechanism because the energy difference between the native conformation and the misfolded conformation with incorrect stacking of helices is not significantly large. When R is as small as 0.2, on the other hand, the chain collapses at first and then helices begin to develop later. It is, however, difficult to find the conformation with well developed helical structures with this mechanism. In the present model the chain folds most efficiently when R has a moderate value; 0.3- R - 3. With this range of R, helix formation and collapse of the chain start simultaneously. Dependence of the folding kinetics on the ratio between the long- and short-range interactions has been discussed by several authors by using lattice models w25,26x. We will discuss this important issue in detail in the further publication. In the present paper, we focus on the folding mechanism under the condition of the balanced long- and short-range interactions: we here use ⌳ k s 0.2 for k F 10, ⌳ k s 0.5 for k ) 10, and ⌳ s 0.9. 2.3. Molecular dynamics simulation
where p and q are types of residues at the ith and jth site. g Ž p,q . s 1 when the ith or the jth site is hydrophobic and g Ž p,q . s 0, otherwise. For large ri j , V Žhydrophobic attraction . is a shallow attractive potential, so that C was chosen to ˚2 . In order to express be a large value, C s 100 A the repulsive interaction for small ri j , r was set ˚ to be 4 A. ⌳ k should have similar value to ⌳ since interactions of Eqs. Ž11. and Ž12. have the similar physical origins. Parameters ⌳ k and ⌳ should not be too large, ⌳ k , ⌳ F 1, in order to make V Žresidue᎐residue. softer potential than V Žmain chain.. Folding kinetics is sensitively dependent on the ratio between strengths of the short- and long-range interactions, R s Ž ⌳ k for k F 10.rŽ ⌳ k for k ) 10.; starting from the extended chain conformation, the folding process is much affected by the value of R. When R is as large as 5, helical structures are formed at first and the chain collapses after that. The folding process is much like the one described by the diffusion-collision model w24x. In the present model, however, the native conformation is not efficiently searched with this
We assume that the chain obeys the following Langevin equation of motion: ª dp ⭸V sy ªy␥ª pq Ž t . , dt ⭸r ª dª r p s , dt m
Ž 13.
where ª r is a 9N dimensional vector, ª rs ª␣ ªO ª ri ,ri ,ri 4 and ª p is the momentum vector. For the simplicity, mass of all atoms are assumed to ª Ž t . is a vector of random force, be equal; m s 1. ª ª ª Ž t . s i Ž t .4 , i s 1 y 9N. Ž t . is assumed to be Gaussian white noise and its variance is related to temperature T by the fluctuation-dissipation relation; ² i Ž t 1 . j Ž t 2 .: s m␥ T␦ i j ␦ Ž t 1 y t 2 . .
Ž 14.
A model chain under the consideration is assumed to have the sequence of calcium-binding protein Ž4cib.. In Eq. Ž13. sequence specific inter-
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
actions are taken into account through V Žresidue᎐residue. in V. We denote units of spatial length, mass, and energy in Eq. Ž13. by r 0 , m 0 , and 0 , respectively. ˚ These units can be roughly assigned as r 0 f 1 A, m 0 f 100 grmol, and 0 f 10 kJrmol. A typical length of time, t 0 is t 0 s r 02rD f 10y10 s,
Ž 15.
where D is the diffusion constant of small molecules in aqueous solution, D f 10y1 0 m2rs. ␥ can be estimated from the Stokes relation, ␥ s 6 a, by assuming to be the viscosity of water at room temperature, f 10y3 kgrmrs and af 10y1 0 m, which yields ␥ f 10y1 2 kgrs f 1 ⭈ m 0rt 0 . We used a slightly underdamped condition, ␥ s 0.6 but the qualitative results did not depend on the precise value of ␥ . Temperature T is measured in the unit of 0rk B . T was slowly decreased as the MD calculation proceeds. Temperature at the nth MD step was T s 0.053rwlogŽ0.0005n q 10.x2 0rk B . MD simulation was performed by numerically integrating Eqs. Ž13. and Ž14. with the discretizing width ⌬ t s 0.2. Initial conformation of the MD simulation was obtained by putting V Žresidue᎐residue. s 0 and stretching the N and C terminals of the random conformation in the opposite direction. Thirty independent MD runs were simulated starting from different initial conformations ª and using different random numbers Ž t .. The discretized equation of motion was integrated up to 10 5 steps. The total length of the simulation is, therefore, 10 5⌬ t ⭈ t 0 f 10y6 s s s. 2.4. Results Difference between the simulated conformation and the X-ray conformation is measured by
509
10 5 steps and 90% of MD runs ended up with ˚ In conformations of conformations of ␦ ) 3.5 A. ˚ ␦ - 2.8 A, four helices are stacking in the same way as in the X-ray structure. Conformations of ˚ on the other hand, are the misfolded ␦ ) 3.5 A, structures with the wrong way of stacking of four helices or with bent or distorted helices. Thus, in the present model calculation, the probability that the MD trajectories reach the native conformation within a s time-scale is f 0.1. These trajectories are termed ‘folding trajectories’. We found that the other trajectories leading to the misfolded conformations are classified into two types: ‘misfolding trajectories’ Ž60% of MD runs.; and ‘escaping trajectories’ Ž30% of MD runs.. A change of ␦ along each type of trajectory is shown in Fig. 5. Along the folding trajectory of Fig. 5a, ␦ gradually decreases to the small value. Along the misfolding trajectories of Fig. 5c, on ˚ the other hand, ␦ does not decrease below 3.5 A throughout the simulation. In the misfolding trajectories, more diverse behaviors of ␦ are found than in the folding trajectories, which is due to diversity of misfolded conformations. It is interesting to see in Fig. 5b that along the escaping trajectories ␦ decreases more rapidly than along the folding trajectories until 4 = 10 4 steps. After ˚ ␦ increases this rapid decrease to ␦ f 3.0 A, ˚ Along the escaping trajecagain to be ␦ ) 3.7 A. tory chain once takes conformations fairly close to the native conformation but these conformations are destroyed in the later stage of the process and the trajectory ‘escapes’ from the native conformation. Overall shape of the conformation is measured by asphericity A w27,28x: As
Ž R1 y R 2 . 2 q Ž R 2 y R 3 . 2 q Ž R 3 y R1 . 2 , 2 Ž R12 q R 22 q R 23 . Ž 17.
2 ␦s Ž N N y 1.
N
N
Ý Ý
ri␣j y riX-ray j
Ž 16.
is1 j)1
where ri␣j s <ª ri␣ y ª r j␣ < and riX-ray is distance j between C ␣ s of the ith and jth site in the X-ray conformation. Ten percent of independent MD ˚ within runs reached conformations of ␦ - 2.8 A
where R12 , R 22 , and R 23 are eigenvalues of the matrix
xx xy xz
xy yy zy
xz yz zz
0
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
510
Fig. 5. Difference between the simulated conformation and the X-ray conformation, ␦ observed along the folding trajectories Ža., along the escaping trajectories Žb., and along the misfolding trajectories Žc.. Two examples Žreal and dotted lines. of trajectories are shown for each case.
and each element of the matrix is defined by xxs
1 N
N
Ý Ž x i␣ . is1
2
,
xy s
1 N
N
Ý x i␣ yi␣ , is1
and so on. Here, x i␣ , yi␣ , and z i␣ are x, y, and z components of the vector ª ri␣. A f 1 when the chain has the stretched rod-like conformation and
A f 0 when the chain takes spherical globular conformation. Typical behaviors of A along the trajectories are shown in Fig. 6. In all kinds of trajectories A rapidly decreases from A f 1 to A f 0.03 within the first 10 4 steps and then fluctuates around A f 0.02: chain collapses rapidly at the early stage of the process within a 0.1- s time-scale and then chain searches the low energy
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
conformation among the compact globular conformations. Changes of A along the folding ŽFig. 6a. and escaping ŽFig. 6b. trajectories show the similar behaviors. There is larger diversity, however, in behaviors of A along the misfolding trajectory ŽFig. 6c.. The folded and misfolded conformations which appear at the late stage of processes have similar asphericity to each other.
511
Results shown in Figs. 5 and 6 are consistent with the picture described in Fig. 1. Type of the trajectory is determined just after the chain collapses to the globule. Thus, a specific collapse is necessary to guide the trajectory to the native conformation in a short time scale. Such explanation based on the difference between specific and non-specific collapses, however, is not enough to
Fig. 6. Asphericity, A, observed along the folding trajectories Ža., along the escaping trajectories Žb., and along the misfolding trajectories Žc.. Two examples Žreal and dotted lines. of trajectories are shown for each case.
512
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 7. A schematic representation of the MD trajectories on the hierarchical energy landscape. There are three types of trajectories: folding trajectories; misfolding trajectories; and escaping trajectories.
clarify the reason why the trajectory ‘escapes’ with a fairly high probability from the native
conformation. Taking account of the escaping trajectories, scheme of Fig. 1 is altered as in Fig. 7. In the following part of this section, differences between the folding and escaping trajectories are more closely investigated. Number of residues involved in helices, n h , is defined by the number of residues whose dihedral angles are in the region y80 - - y25⬚ and y90⬚ - - 0⬚. Change of n h along the trajectories is shown in Fig. 8. Both in folding and escaping trajectories, n h increases rapidly up to n h f 20 at the early stage of process: almost 40% of helices are formed before chain completes to collapse. After compaction, n h increases steadily as the chain searches the low energy conformation; n h in the escaping trajectories, however, does not grow as large as in the folding trajectories: growth of helices is limited to the level of 80% of the native conformation in the escaping trajectories. We denote four helices in the native conformation by helix-I with I s 1᎐4 counting from the N terminal side to the C terminal side. When helix-I is composed of the residues i,i q 1,...,i q a and helix-J is composed of the residues j, j q 1,..., j q b, then distance between helix-I and helix-J is
Fig. 8. Helicity, n h , observed along the folding trajectories Ža. and along the escaping trajectories Žb.. A horizontal line shows the value in the X-ray conformation. Two examples Žreal and dotted lines. of trajectories are shown for each case.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
defined by hI J s
1 aq 1
iqa
Ý min Ž rk␣j ,rk␣jq1 ,...,rk␣jqb . ,
Ž 18.
ksi
where min Ž.... is the function to select the minimum value. When helix-I is placed in the direction at angle I J to the direction of helix-J, CI J s cos I J is defined by the inner product of two unit vectors; CI J s
ž ⭈
2 aŽ aq 1 .
ž
iqay1
2 Ž b bq 1 .
Ý ksi
ª␣ r k yª r l␣ ª␣ ª␣ lskq1 < r k y r l <
jqby1
Ý X
k sj
iqa
Ý
/
ª␣ r k X yª r l␣X ª␣ ª␣ X X l sk q1 < r k X y r l X <
jqb
Ý
/
.
Ž 20.
These definitions of h I J and CI J have clear meanings only in the conformations which have well developed helices. h I J and CI J monitored along the MD trajectory, however, should give information of the growth process of tertiary structure. Changes of h13 and h 23 are shown in Fig. 9 and changes of C14 and C34 are shown in Fig. 10. Results shown in Figs. 9 and 10 are consistent with the results shown in Fig. 5. Along the folding trajectories the tertiary structure is gradually formed during the structural search among the compact globular conformations. Along the escaping trajectories, on the other hand, the tertiary structure develops rather rapidly and the conformations which have features common to the native conformation are formed at the early stage of the process. These features, however, are destroyed as the MD steps proceed and the chain reaches conformations which have globally different structures from the native one. Differences between the folding and escaping trajectories suggest that consistency between the secondary structure formation and the tertiary structure formation is a necessary condition to guide the trajectories to the native conformation. When the strengths of the short- and long range interactions are balanced, approximately half of the secondary structures are formed when the chain collapses to the globule. Then, both the secondary and tertiary structures start to develop
513
cooperatively in the compact globular state. When the tertiary structural order is formed much earlier than the secondary structure, then the later growth of the secondary structure can bring the chain to the different tertiary conformation. To describe this escaping process, we could use an analogy of solving a puzzle like Rubik’s cube. When one side of the cube is arranged to have the correct order, then the later arrangement of the other sides often destroys the side which was arranged in a earlier step. In order to arrange multiple sides of the cube, orders of different sides have to be developed in a synergetic and consistent way. Present model calculations suggest that competition and cooperation of different structural formations require the consistent growth of different structural orders for the chain to fold fast with s scale. We should call this principle for fast folding ‘the consistency principle on the kinetic process’, which is a reminiscent of the consistency principle first discussed by Go based on the analyses of the 1st order nature of the folding᎐unfolding transition and on the analyses of the stability of the native conformation w29x. 3. Spin-glass-like model 3.1. Hierarchical dynamics In this section we examine the idea of the consistency principle on the kinetic process by constructing a simple spin-glass model. Two kinds of spin variables, Si for i s 1 y N and TI for I s 1 y M, are introduced to represent the coarse-grained secondary and tertiary structures of the protein. We assume S i and TI take values of "1. Si represents the secondary structure: when the helical structure is formed at the ith residue, then Si s 1. Si s y1, otherwise. TI represents the direction of helix-I. When TI s 1, helix-I is placed in the same direction as in the native conformation. When TI s y1, helix-I is in a different direction. Each helix is assumed to be composed of NrMs K residues: residues i s 1,...,K are the constituents of helix-1, residues i s Kq 1,...,2 K are the constituents of helix-2, and so on. We write i g I if i belongs to the
514
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 9. Distance between helix-1 and helix-3, h13 , observed along the folding trajectories Ža. and along the escaping trajectories Žb., and distance between helix-2 and helix-3, h 23 , observed along the folding trajectories Žc. and along the escaping trajectories Žd.. A horizontal line shows the value in the X-ray conformation. Two examples Žreal and dotted lines. of trajectories are shown for each case.
region of the sequence where helix-I is formed: residues i g 1 for i s 1,...,K, i g 2 for i s Kq 1,...,2 K, and so on Žsee Fig. 11.. In the native conformation Si s TI s 1 for all i and I. Energy H of the system is
Ny1
Hsy
Ý
Ž rJiiq1 q 1 . Si Siq1
is1
y
1 R st
My1
M
Ý Ý Ž r t JI J q 1. TI TJ . Is1 J)1
Ž 21.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
515
Fig. 10. Inner product of unit vector parallel to helix-1 and the one parallel to helix-4, C14 , observed along the folding trajectories Ža. and along the escaping trajectories Žb., and inner product of unit vector parallel to helix-3 and the one parallel to helix-4, C34 , observed along the folding trajectories Žc. and along the escaping trajectories Žd.. A horizontal line shows the value in the X-ray conformation. Two examples Žreal and dotted lines. of trajectories are shown for each case.
Here, the first term is energy of formation of helices, which is a sum of the short-range interactions. The second term is the helix᎐helix interactions which are the long-range interactions. Ji iq1 is Gaussian random number with mean 0 and
standard deviation 1r62. JI J is Gaussian random number with mean 0 and standard deviation 1r 'My 1 . r and r t are parameters which determine the extent of randomness. When r or r t is larger, interactions are more frustrated and
516
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 11. Si represents whether ith site is helical Ž Si s 1. or not Ž Si s y1. and TI represents whether helix-I is aligned in the same direction as in the native conformation ŽTI s 1. or not ŽTI s y1.. A case for N s 16 and Ms 4 is shown.
energy landscape is more rugged. R st is a parameter which determines the balance between the short-range interactions and the long-range interactions, which resembles to the parameter R in V Žresidue᎐residue. of Section 2. For most of the calculations we used r s 0.3, r t s 0.3, and R st s 1.67. N and M were set to be N s 30 and Ms 10 or N s 16 and Ms 4. Si and TI are changed with the Monte Carlo ŽMC. method. We assume that Si moves much faster than TI and that motion of TI is slaving to the motion of Si : to generate a trial configuration to update the MC step, sign of Si is flipped, at first, at the randomly chosen position i. When i g I, sign of TI is altered with the probability 0 - ␣ - 1 and is kept unchanged with the probability 1 y ␣ . Trial configuration prepared in this way is judged to be accepted or rejected by the Metropolis algorithm with temperature T. This MC motion is schematically shown in Fig. 12. It should be noted that S i and TI are uncoupled in the energy expression of Eq. Ž21., so that the equilibrium properties of the system are described just in terms of two independent spin systems. Si and TI , however, are kinetically strongly coupled to each other, so that when the MC trajectory is truncated with the relatively small number of steps, then the probability of the folding trajectory, , should be decisively affected by ␣ which controls the kinetic coupling between Si and TI . We here use this spin model to focus on such kinetic effects on the fast folding system. T dependence of is shown in Fig. 13. Starting from a random configuration, MC trajectories of 3000 steps were performed. Temperature T was kept constant through 3000 steps. From Fig.
Fig. 12. A schematic representation of folding trajectory in the spin-glass like system. When Si is flipped, TI is altered with the probability ␣ . Vetical columns are energy spectra for fixes configuration of TI . Arrows of real lines are folding trajectories and arrows of dotted lines are escaping trajectories. Both the secondary and tertiary structures have to develop in a cooperative way to reach the native conformation.
13, we can see that folding is difficult both in the high and low temperature regions. For T ) 0.5 random conformations are favored by the strong thermal fluctuation. For T - 0.2 system is trapped into one of many random minima in energy. Thus, folding is most efficient in the temperature range 0.2- T - 0.5. In the rest of this section, in order to compare the results to those of Section 2, T is gradually decreased as MC steps proceed: T at the nth step is T s T Ž0. y w T Ž0. y T Ž1.x nr15 000 for n - 15 000 and T s T Ž1. for n G 15 000. We used T Ž0. s 0.5 and T Ž1. s 0.2. 3.2. Results Distance between the simulated conformation and the native conformation is measured by d: ds
1 N
N
1
M
Ý Ž1 y Si . r4 q M Ý Ž1 y TI . r4,
is1
Ž 22.
is1
which resembles ␦ in Section 2. Difference
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 13. Dependence of T on the probability of the folding trajectories . Results of 2000 independent MC runs are averaged. N s 16, Ms 4 and ␣ s 0.2.
between the secondary structure of the simulated conformation and the one of the native conformation is measured by
dS s
1 N
N
Ý Ž1 y Si . r2.
Ž 23.
is1
N Ž1 y d S . corresponds to n h of Section 2. In Fig. 14 typical examples of H, d, and d S observed along the folding trajectory are shown. Features of the trajectory are similar to those of the folding trajectories shown in Figs. 6a and 9a. Both the short- and long-range structures develop steadily in a cooperative way. Trajectories which did not reach the native conformation within 3000 steps were further elongated to 30 000 steps. Most of these trajectories were escaping trajectories for this parameterization. Typical example of escaping trajectory is shown in Fig. 15. Features of the trajectory is similar to those of the escaping trajectories shown in Fig. 6b and Fig. 9b. When the long-range structural order is formed at too early stage, then the growth of short-range structure in the later stage destabilizes the early formed struc-
517
tural order. Thus, this simple spin-glass-like model is able to reproduce the distinction between the folding trajectories and the escaping trajectories found in the knowledge-based potential model. is a decreasing function of r and r t : s 1 when r s r t s 0 and f 0 when r s r t ) 1.0. For large r and r t , the system is frozen into the glassy random conformation and does not reach the native conformation. For small r and r t the ratio ␦ Er⌬ E should be large, where ␦ E is energy gap between energy of the native conformation and the average energy of other conformations, and ⌬ E is the average ruggedness in the energy landscape w1,2x. This design of r and r t for fast folding is the requirement that the funnel should be well developed around the native conformation. Such a design principle was called the principle of minimal frustration w1,2,10x. When R st s 5, secondary structures develop quickly but the long-range order does not grow within 30 000 steps. When R st s 1.25, on the other hand, the long-range structures develop quickly but it takes as long as 30 000 steps for secondary structures to develop. Folding is most efficient in the intermediate range, 1.5- R st - 2: balance between strengths of the short- and long-range interactions is necessary to make large. This design principle was called the consistency principle w29x. Alternation of R st changes the energy spectrum and thus R st affects ␦ Er⌬ E: the consistency principle has a great overlap on its content with the principle of minimal frustration. R st , however, affects not only the first moment, ␦ E, and the second moment, ⌬ E, of the energy spectrum but also the higher moments of the spectrum. Therefore, we should note that the consistency principle refers not only to the existence of the funnel but also to the rich phenomena relating to the inner structure of the funnel. Dependence of on ␣ is shown in Fig. 16. is largest in the range 0.2- ␣ - 0.4. For ␣ - 0.1, misfolding trajectories are dominant and for ␣ ) 0.5, escaping trajectories are dominant. In order to make large, ␣ has to be tuned to assure the consistent growth of both the short- and longrange structural orders. This design principle is
518
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 14. Energy, H, distance to the native conformation, d, and difference in the secondary structures from the native conformation, d S, observed along an example of the folding trajectory. N s 30, Ms 10 and ␣ s 0.2.
the consistency principle on the kinetic process discussed in the previous section. By changing ␣ , ␦ Er⌬ E is not altered but the kinetic connectivity among conformations is altered. Thus, the toy model in this section suggests that design on the kinetic consistency, that is, design of the inner structure of the funnel is important to design the fast folding sequences. Effects of the connectivity
among conformations on the folding speed were discussed by Betancourt and Onuchic w30x by using the lattice model. In the present paper we stress the idea that characteristics of the hierarchical structures which have not been taken account of in the lattice model should be an important factor to design the connectivity among conformations.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
519
Fig. 15. Energy, H, distance to the native conformation, d, and difference in the secondary structures from the native conformation, d S, observed along an example of the escaping trajectory. N s 30, Ms 10 and ␣ s 0.2.
4. Discussion
It was shown from the results of the knowledge-based potential model and a spinglass-like model that consistent and cooperative growth of the short-range and the long-range
structural orders is a necessary condition for the trajectory to reach the native conformation in a short time. Therefore, to design the fast folding peptides, sequence has to be selected to assure the kinetic consistency among different structural orders. It should be interesting to ask whether the
520
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521
Fig. 16. Dependence of kinetic parameter ␣ on the probability of the folding trajectories . Results of 2000 independent MC runs are averaged. N s 16, Ms 4.
kinetic consistency is always in harmony with the stabilization of the native conformation or these two requirements sometimes contradict each other. This point might be relevant to the difficulty encountered in de novo design of peptides w31x. For most of cases, de novo design of peptides is based on the stability condition of the native conformation. It should not be clear whether such design principle gives an optimal sequence from a view point of the kinetic consistency. The same question might be relevant to the reason why natural proteins are marginal in their stability w32x. Free energy difference between the unfolded and the folded states are, for most of small globular proteins, less than 10᎐20 kcalrmol. This value is quite small considering the fact that small proteins consist of more than several hundreds atoms. The native conformation of natural protein is, therefore, not thoroughly stabilized but is marginally stable. Kinetic requirements could be reasons for this marginal stability. Thus, further investigations on the relation between stability and kinetic consistency should be important. There is much room to improve the knowledge-based potential model used in this paper. Especially important is to incorporate the
data of high energy states into potentials. As shown in the present paper the way of compaction is a large factor to determine . Thus, potentials must have ability to reproduce the compaction process in a proper way. Statistical data of extended conformations and compact globular conformations other than the native conformation should be used to construct such potentials. It should be important to check whether consequences of the present paper, such as existence of three types of trajectories and necessity of kinetic consistency in the folding trajectories should hold in the improved models of folding. Important insights will be obtained from experimental data on the growth of the secondary structure and on the tertiary structure measured at the same time in s time scale. Such data can provide information on how the kinetic consistency is realized in natural proteins and how the degree of consistency can be altered by the protein engineering technique. Through the comparison between experiments and simulations, more precise understanding of the structure of the funnel will be obtained, which should provide deep insights on how the fast folding proteins have been designed in the evolutionary history. References w1x J.D. Bryngelson, J.N. Onuchic, N.D. Socci, P.G. Wolynes, Proteins 21 Ž1995. 167. w2x J.N. Onucic, Z. Luthey-Schulten, P.G. Wolynes, Annu. Rev. Phys. Chem. 48 Ž1997. 539. w3x K.A. Dill, H.S. Chan, Nature Struct. Biol. 4 Ž1997. 10. w4x E.I. Shakhnovich, Curr. Opin. Struct. Biol. 7 Ž1997. 29. w5x W. Eaton, P.A. Thompson, C.-K. Chan, S.J. Hagen, J. Hofrichter, Structure 4 Ž1996. 1133. w6x W. Eaton, et al., Curr. Opin. Struct. Biol. 7 Ž1997. 10. w7x P.E. Leopold, M. Montal, J.N. Onuchic, Proc. Natl. Acad. Sci. USA. 89 Ž1992. 8721. w8x C.J. Camacho, D. Thirumalai, Proc. Natl. Acad. Sci. U.S.A. 90 Ž1993. 6369. w9x A. Sali, E.I. Shakhnovich, M. Karplus, J. Mol. Biol. 235 Ž1994. 1614. w10x J.D. Bryngelson, P.G. Wolynes, Proc. Natl. Acad. Sci. U.S.A. 84 Ž1997. 7524. w11x Z. Guo, D. Thirumalai, Biopolymers 36 Ž1995. 83. w12x D. Thirumalai, J. Phys. I ŽParis. 5 Ž1995. 1457. w13x M. Sasai, Proc. Natl. Acad. Sci. USA. 92 Ž1995. 8438. w14x B.A. Shoemaker, J. Wang, P.G. Wolynes, Proc. Natl. Acad. Sci. USA. 94 Ž1997. 777.
S. Saito, M. Sasai r Journal of Molecular Structure (Theochem) 461᎐462 (1999) 503᎐521 w15x S. Saito, M. Sasai, T. Yomo, Proc. Natl. Acad. Sci. USA. 94 Ž1997. 11324. w16x S.S. Plotkin, J. Wang, P.G. Wolynes, J. Chem. Phys. 106 Ž1997. 2932. w17x T.R. Sosnick, L. Mayne, S.W. Englander, Proteins 24 Ž1996. 413. w18x T. Yomo, S. Saito, M. Sasai Žto be published.. w19x S. Kamtekar, J.M. Sciffer, H. Xiong, J.M. Babik, M.H. Hecht, Science 262 Ž1993. 1680. w20x M.S. Friedrichs, R.A. Goldstein, P.G. Wolynes, J. Mol. Biol. 222 Ž1991. 1013. w21x J.S. Richardson, Adv. Protein Chem. 34 Ž1981. 168. w22x M. Hendlich, P. Lackner, S. Weitckus, H. Floeckner, R. Froschauer, K. Gottsbacher, G. Casari, M.J. Sippl, J. Mol. Biol. 216 Ž1992. 167. w23x M.J. Sippl, J. Mol. Biol. 213 Ž1992. 859.
521
w24x M. Karplus, D.L. Weaver, Nature 260 Ž1976. 404. w25x N. Go, H. Taketomi, Proc. Natl. Acad. Sci. USA. 75 Ž1979. 559. w26x A. Rey, J. Skolnick, Proteins 16 Ž1993. 8. w27x J. Rudnick, G. Gaspari, J. Phys. A 19 Ž1986. L191. w28x O.F. Olaj, G. Zifferer, J. Chem. Phys. 100 Ž1994. 636. w29x N. Go, Ann. Rev. Biophys. Bioeng. 12 Ž1983. 183. w30x M.R. Betancourt, J.N. Onuchic, J. Chem. Phys. 103 Ž1995. 773. w31x K.A. Dill, S. Bromberg, K. Yue, K.M. Fiebig, D.P. Yee, P.D. Thomas, H.S. Chan, Protein Sci. 4 Ž1995. 561. w32x P.G. Wolynes, in: H. Bohr, S. Brunak ŽEds.., Proceedings on Symposium on Distance-Based Approaches to Protein Structure Determination II, CRC Press, Boca Raton, 1994, p. 3᎐17.