Structural and Pathway Complexity of β-Strand Reorganization within

J. Phys. Chem. B 2007, 111, 5425-5433

5425

Structural and Pathway Complexity of β-Strand Reorganization within Aggregates of Human Transthyretin(105-115) Peptide Da-Wei Li,† Li Han,‡ and Shuanghong Huo*,† Gustaf H. Carlson School of Chemistry and Biochemistry and Department of Mathematics and Computer Science, Clark UniVersity, 950 Main Street, Worcester, Massachusetts 01610 ReceiVed: January 14, 2007; In Final Form: March 6, 2007

Interstrand conformational rearrangements of human transthyretin peptide (TTR(105-115)) within dimeric aggregates were simulated by means of molecular dynamics (MD) with implicit solvation model for a total length of 48 µs. The conformations sampled in the MD simulations were clustered to identify free energy minima without any projections of free energy surface. A connected graph was constructed with nodes ()clusters) and edges corresponding to free energy minima and transitions between nodes, respectively. This connected graph which reflects the complexity of the free energy surface was used to extract the transition disconnectivity graph, which reflects the overall free energy barriers between pairs of free energy minima but does not contain information on transition paths. The routes of transitions between important free energy minima were obtained by further processing the original graph and the MD data. We have found that both parallel and antiparallel aggregates are populated. The parallel aggregates with different alignment patterns are separated by nonnegligible free energy barriers. Multiroutes exist in the interstrand conformational reorganization. Most visited routes do not dominant the kinetics, while less visited routes contribute a little each but they are numerous and their total contributions are actually dominant. There are various kinds of reptation motions, including those through a β-bulge, side-chain aided reptation, and flipping or rotation of a hairpin formed by one strand.

Introduction Amyloid deposits formed by full-length amyloidogenic proteins or by fragments of large precursor proteins are associated with several types of amyloidoses. Emerging evidence has shown that soluble oligomers are primary cytotoxic species.1-4 Therefore, toward finding a cure for amyloid diseases, it is crucial to understand the mechanism of early formation of oligomers. Human transthyretin (TTR) is one of the 20 or so human amyloidogenic proteins. Although it belongs to the category of full-length amyloidogenic protein, TTR fibrils also contain fragments of protein.5 The solid-state NMR spectroscopy has revealed that TTR(105-115) peptide adopts an extended conformation in the fibril state.6,7 However, the analysis of the 15N, 13CO, 13C , and 13C chemical shifts only provides the R β information on backbone torsion angles6,7 without interstrand information, such as parallel vs antiparallel and in-registry vs off-registry. In addition, detailed experimental characterization of the mechanism of formation of oligomers is very difficult because of the solubility and the wide range distribution of conformations. Computational approaches that complements experimental methods have been used to investigate the early aggregation process, including reduced protein models8-17 and full atomic models.18-25 It has been proposed that the small oligomericaggregatesactasthenucleationsiteoffibrillogenesis26-28 and the kinetically trapped dimers/oligomers with a slow relaxation time might be involved in the structural evolution in the prenucleation stage.18 The formation of the initial oligomers * To whom correspondence should be addressed. Phone: 508-793-7533. Fax: 508-793-8861. E-mail: [email protected]. † Gustaf H. Carlson School of Chemistry and Biochemistry. ‡ Department of Mathematics and Computer Science.

and the distribution of conformations certainly affect the following aggregation process. Recently, a series of isotope-edited FTIR spectroscopic studies have demonstrated that there is a relaxation process from the initial β-sheet aggregates to the final product, involving strand realignment by the repeated detachment and annealing of strands in solution as well as reptation motion within an aggregate.29-31 The later has also been observed in simulations,14 but the underlying mechanisms for the reptation motion are not clear. To shed lights on the detailed mechanism of reptation motion, this work focuses on the conformational transitions within aggregates of TTR(105-115), which adopt parallel and antiparallel β-sheet arrangements displaced by various number of residues. Some intuitive and widely used approaches for the analysis of free energy surface rely on the projections of the surface on some progress variables, such as the number of native contacts, the radius of gyration, and the orientations of different chains and/or principal components. These approaches, however, might miss possible intermediate states and the transitions involved,32-34 even when multiple projections on different progress variables are used. In this work, we take a different approach to study the free energy surface as captured in our molecular dynamics (MD) simulations by analyzing the simulation data with some graph algorithms, as utilized by Krivov and Karplus on β-hairpin folding.32 In brief, we constructed a connected graph that reflects the complexity of the free energy surface and the transition pathways connecting free energy minima. This graph in turn was used to extract transition disconnectivity graph (TRDG)32 that contains the free energy minima and the overall free energy barriers between minima but no information on transition pathways. The routes of

10.1021/jp0703051 CCC: $37.00 © 2007 American Chemical Society Published on Web 04/14/2007

5426 J. Phys. Chem. B, Vol. 111, No. 19, 2007 transitions between important free energy minima were obtained by further processing the original graph and the MD data. Materials and Methods Molecular Dynamics Simulations. Two chains of TTR(105-115) peptide were simulated with the N and C terminals uncapped in a periodically replicated cubic box (side ) 50 Å). The sequence of this fragment is YTIAALLSPYS. The molecular dynamics simulations were performed using CHARMM and param19 force field.35 To enhance sampling and speed up the simulations, we employed the EEF1 implicit solvation model to describe the aqueous environment.36 EEF1 has been employed to study protein folding and unfolding successfully.32,37 Since we used Langevin dynamics in the low friction limit (friction coefficient ) 0.1 ps-1) to take into partial account the friction and the random collision between solute and solvent, the events took place in a faster time scale in our simulations than in experimental studies. However, the thermodynamics of the system is not affected by the low viscosity of the solvent. Nonbonded interactions were truncated with a 9 Å cutoff, as required by the EEF1 model. Initially, the two chains were in extended state and 25 Å apart, which was larger than the nonbonded interaction cutoff. The temperature of the system was raised to 360 K for an efficient sampling. All bonds involving hydrogen atoms were constrained.38 A time step of 2 fs was used and the nonbonded neighbor list was updated heuristically. Conformations were saved every 20 ps for analysis. Eight independent trajectories were collected for 6 µs each, totaling 48 µs simulations with 2.4 million conformations generated. Cluster Analysis. To study the ensemble of conformations generated by the MD simulations and the underlying free energy surface, we used clustering analysis to group similar conformations, with the resulting clusters corresponding to free energy minima. For conformation similarity measure, we used the allatom root-mean-square deviations (RMSDs) and the main-chain (intrastrand and interstrand) hydrogen bond patterns. We think that information on hydrogen bonds is particularly useful for the clustering analysis of the conformations in the dimerization study because these hydrogen bonds are sensitive to the displacement of two chains and their formation contributes a large amount of effective energy gain in aggregation. In our study, a hydrogen bond was considered to form if the corresponding hydrogen atom and the acceptor were within 2.4 Å, without taking into account the angle restriction. The hydrogen bonds between neighboring residues in sequence (di,j < 3) were excluded. We applied the clustering algorithm introduced by Daura et al.39 to generate a neighbor list for each conformation, where we considered two conformations to be neighbors if their RMSD was no greater than 4.0 Å and their hydrogen bond patterns were off at most by 1. Then the conformation with the largest number of neighbors was identified and used along with all its neighbors to form the first cluster. All the conformations of this cluster were then (conceptually) eliminated from the ensemble of conformations, and the neighbor list of all other conformations was updated. This process was repeated until the ensemble was empty. By this way, a series of non-overlapping clusters of conformations was obtained. For each cluster, the conformation with the most neighbors was used as the representative of the cluster and would be referred to as the cluster center in the rest of the paper. This algorithm, however, is computationally expensive because the required memory space and CPU hours scale up with n2, where n is the number of

Li et al. conformations. As an example, it needs hundreds of gigabytes of memory to store the neighbor list for 2.4 million conformations. To overcome these difficulties, we used the first 600 000 conformations as a training set and obtained the 300 most visited clusters. Afterward, the rest of the 2.4 million conformations were compared with these 300 clusters and a given conformation would join one of the 300 clusters if it was a neighbor of the cluster center. All of the conformations that did not belong to these 300 clusters would be processed by the clustering algorithm to obtain new clusters. Connected Graph and the Minimal Cut Tree. From the clustering result and the MD data, we then constructed a connected graph with nodes and edges corresponding to clusters and transitions between the clusters. We also assigned a weight to each node and edge; the weight of a node is the number of times that the cluster was visited and the weight of the edge between two nodes was the number of direct transitions between these two clusters as observed in the MD simulation. This graph well describes the thermodynamics and kinetics of the system and can be used to extract much useful information. For example, the free energy of cluster i is computed as

Fi ) -kBT ln Zi

(1)

where Zi is the partition function of this cluster and is taken to be the number of times it has been visited in the MD simulation and, thus, is the node weight. Again, we consider that each node represents a free energy minimum. Clearly, based on eq 1, the node with the largest weight corresponds to the global minimum of the free energy surface. This connected graph can also be used to extract the overall transitions between minima, which reflect the information on free energy barrier of the transitions. For this task, we used the Maxflow problem formulation40 and minimal cut tree algorithm41 from graph theory, as done by Krivov and Karplus on β-hairpin folding.32,42 The minimal cut tree algorithm was employed to compute the maximum flow values ()minimum cut values) for all pairs of graph nodes, resulted in a minimal cut tree (or transition disconnectivity graph, TRDG32,42), which was in turn used to compute the overall free energy barrier between clusters. More specifically, the free energy barrier between two clusters i and j is

Fij ) -kBT ln Zij

(2)

where Zij is the partition function of the barrier and is related to the minimum cut value (nij) by

1 h 1 Zij ) nij 2 kBT ∆t

(3)

where h is the Plank constant, kB is the Boltzmann constant, T ) 360 K, ∆t ) 20 ps is the sampling interval. The factor of 2 in eq 3 is due to the definition of nij as the sum of the numbers of transitions for ifj and jfi. Note that the minimal cut tree as shown in Figure 1 reflects the overall free energy barriers between clusters, including the contributions of all possible routes, direct or indirect. And the free energy barrier calculated from all possible routes is comparable to experimental results, whereas the free energy barrier from the direct route alone is not. To describe the interstrand side-chain contact, we defined a contact index. Two side chains are considered to be in contact (contact index ) 1) if the distance between one or more pair(s) of their heavy atoms is less than 4.5 Å; otherwise, they are not in contact (contact index ) 0). As shown in Figure 2, all

β-Strand Reorganization in TTR(105-115) Aggregates

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5427 in the connected graph have low weights. Thereby, the paths directly observed from the MD trajectories are shown in Figure 4. Chemical Shift Calculation. We used the typical conformations obtained in the MD trajectories to back-calculate the chemical shifts for N (of NH), C (of CO), CR, and Cβ and compared them with the experimental results.7 The chemical shifts were calculated using the ShiftX program.43 We used the total deviation (∆δtot.), which represents the standard deviation from the averaged chemical shift for residues 106-114, to measure the difference between the simulation and experimental results. The definition of ∆δtot. was taken from ref 22.

Figure 1. Minimal cut tree, showing only the first 300 lowest free energy minima. Six typical conformations corresponding to the centers of some important clusters are depicted. From left to right: parallel out-of-registry displaced by three residues (minima 16 and 19); antiparallel and in-register β-sheet (minimum 9); parallel and in-register (minima 10 and 11); parallel and out-of-registry displaced by one residue (minima 1, 2, 3, and 4); bent antiparallel in-registry (minimum 77 has the lowest average effective energy). Different residues are in different colors. The N-terminal is in blue, while the C-terminal is in red. The red circle indicates the free energy barrier between cluster 9 and the central basin including clusters 1-4.

conformations belonging to one cluster were calculated, and the gray scale map represents the mean value of the contact index. The minimal cut tree does not contain the information on transition paths. To illustrate the paths connecting the parallel in-register aggregates and the parallel off-register arrangements displaced by one residue, we extracted transition pathways connecting these free energy minima from the original connected graph. To highlight the routes with significant contributions to the transitions, we further processed the original connected graph by merging nodes that were connected by edges with weights greater than 700 (large edge weight means a low free energy barrier and fast transition rate) to form a basin, in which each basin is named after its most visited node. We also ignored the edges with weights smaller than 65. Only the through routes were preserved and any route that led to a dead end (due to the removal of low weight edge) was removed. The resulting graph that shows important basins as well as their transitions is given in Figure 3. Note that the number of edges between two basins in Figure 3 corresponds to the number of edges between the clusters of the basins in the original connected graph, but only the total weights are included in the figure. As mentioned above, the threshold of edge weight is 65. For the transitions between free energy minima with high free energy barriers, such as those between the antiparallel and parallel aggregates, we did not observe any transitions with (statistically) significant contributions because all involved edges

Results and Discussion In our simulation, dissociation (defined as zero interaction energy between the two chains) was not observed in any 6 µs trajectory, indicating that the hydrogen bonding and the hydrophobic interaction were too strong to allow dissociation within the simulation time. However, the systems could still transit among different conformations. The dimer visited parallel and antiparallel aggregates with different alignments between the two chains many times in each trajectory. Therefore, the following analyses focus on the β-strand reorganization within aggregates. Cluster and Minimal Cut Tree. Figure 1 shows the minimal cut tree including only the 300 lowest free energy minima. Other minima are not shown for the clarity of the figure. The minima are numbered from the lowest in free energy to the highest. The 20 lowest free energy minima can be grouped into five basins: (1) clusters 1 and 2 as well as the clusters connected to them with large edge weights; (2) clusters 3 and 4 as well as the clusters connected to them with large edge weights; (3) cluster 9; (4) clusters 10 and 11; (5) clusters 16 and 19. Even though cluster 77 is high in free energy, it has the lowest effective energy; therefore, we also depict it in Figure 1. As shown in Figure 1, the five basins plus cluster 77 fall into four categories of alignments: (1) parallel out-of-registry displaced by one residue (minima 1, 2, 3, and 4); (2) parallel in-registry (minima 10 and 11); (3) parallel out-of-registry displaced by three residues (minima 16 and 19); and (4) antiparallel inregistry (minima 9 and 77). Minima 1 and 3, 2 and 4, 10 and 11, and 16 and 19 are identical to each other, respectively, if the two chains are considered to be indistinguishable. The slight difference in their computed free energy stems from the small deviations in the number of visited times in the simulation data. The overall shape of free energy does not change, but the relative free energies and the barriers change if the two chains are treated as indistinguishable. The reason why we treated the two chains distinguishable in our clustering procedure is that it can provide more information on the kinetics (see the details in the section of Kinetics). Minima 1 and 2 have the same pattern of alignment but different numbers of hydrogen bonds and side-chain packing. The number of conformations in clusters 1, 2, 9, 10, 16, and 77 is 147 145, 119 536, 39 512, 30 923, 15 473, and 2918 (out of 2.4 million), respectively. In fact, the clusters connected to these clusters with large edge weight have similar conformations. Therefore, the total population of the clusters in the central basin (in Figure 1) is ca. 1.2 million that includes clusters 1, 2, 3, and 4 themselves and the clusters with the free energy barrier lower than 2.1 kcal/mol separated from them. In the clusters with low free energies, R-helical dimers are not found. However, one turn of R-helix formed in one strand has been observed (see detailed discussions in the section of Kinetics). Recall that in our clustering analysis, we considered two conformations to be neighbors if their RMSD was no greater

5428 J. Phys. Chem. B, Vol. 111, No. 19, 2007

Li et al.

Figure 2. Interstrand side-chain contact. Panels a-f show the interstrand side-chain contact map of the cluster 1, 2, 9,10, 16, and 77, respectively. The residues 105-115 are labeled with their letter names. The interstrand alignment pattern is illustrated in the x-axis. The conformations are described as parallel or antiparallel and in-register or off-register (displaced by number of residues). The arrow emphasizes the aromatic stacking of Y114-Y114 and side-chain hydrogen bonding of Y114-S115.

than 4.0 Å and their hydrogen bond patterns were off at most by 1. We chose such a similarity measure because it can correctly distinguish not only hydrogen-bonding pattern but also secondary structure distortion. In principle, the transitions between conformations within a single cluster should be significantly more frequent than those between conformations in different clusters. In practice, this requirement is not 100% fulfilled because no similarity measure available to date is 100% correlated with free energy. When one chooses clustering methods and thresholds, the methods that can satisfy the above requirement to the largest extent should be used. We have also tried other clustering algorithms as well as different thresholds and found that both free energy minima and barriers depend on the clustering method and threshold (see Supporting Information). However, the overall shape of free energy is similar to Figure 1 (Supporting Information Figures S3 and S4) if the clustering method and the threshold are reasonably good. As to the issue of sampling sufficiency, we have found that the conformations of the 100 most visited states and the free energy barriers between them based on the ensemble collected from eight trajectories are the same as those generated from only four or even two trajectories. Therefore, we believe our simulation is long enough to yield reasonable results near most visited states. However, the total number of clusters kept increasing when more trajectories were added. This suggests that the sampling in the less visited phase space is not enough, though the total simulation time is 48 µs. Most Visited Clusters. In Figure 2, we illustrate the interstrand side-chain contacts and the alignment pattern of the interstrand main-chain hydrogen bonding of most visited clusters. To summarize the intermolecular interaction and free energy of most visited nodes, we list the number of interstrand main-chain hydrogen bonds, interstrand side-chain contact index, effective energy, free energy, and free energy barrier separating them from the global minimum in Table 1. Node 1 is the global minimum of free energy because it is the most visited state. In

general, there are two possible interstrand main-chain hydrogen bonding patterns for the parallel off-register β-sheet displaced by one residue; one is that all the hydrogen bonds are formed between odd number residues; the other is that involving the even number residues only. Because residue 113 is proline, only the even number residues are involved in the interstrand mainchain hydrogen bonding in node 1(3). Although both node 1 and node 2 are parallel off-register aggregates displaced by one residue, on average node 1 has three more interstrand mainchain hydrogen bonds and one less side-chain contact index than node 2 (Table 1). As shown in Figure 2a, there is no interchain contact between side chains of Y114 in node 1(3) (contact index < 0.01), while in node 2(4) (Figure 2b), this contact index is 0.78. Besides the hydrophobic packing between Y114 and Y114 (from the other chain), S115 also interacts with Y114 via side-chain hydrogen bonding in node 2. As a result, the intermolecular interaction energy of node 2 is 3.4 kcal/mol lower than that of node 1, mainly attributed to the solvation free energy difference (Table 1). The conformations in node 9 are organized into an antiparallel in-register β-sheet (Figure 1). Even though the conformations in this cluster have the most interstrand main-chain hydrogen bonds, and its interstrand sidechain contact index is almost as high as that of node 2, but its intermolecular interaction is less favorable than that of node 2, giving rise to higher effective energy than that of node 2 but slightly lower energy than that of node 1. The aggregates in node 10 adopt a parallel in-register conformation with five interstrand main-chain hydrogen bonds (Figure 1). Their sidechain contact pattern is very similar to that of node 2 (Figure 2b,d) but with smaller contact index. Node 16 which is a parallel off-register aggregate displaced by three residues has almost the same number of hydrogen bonds as node 10 but the least interstrand contact index in Table 1. Node 77 is also antiparallel in-register as node 9, while the aggregates bend in “U” shape (Figure 1) to maximize intrastrand and interstrand side-chain contacts. As a result, besides the interstrand contact along the

β-Strand Reorganization in TTR(105-115) Aggregates

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5429

Figure 3. (a) Graph showing the paths that connect the parallel in-register aggregates (basins 10 and 11) and the parallel off-register arrangements displaced by one residue (basins 1 and 3). Four typical intermediate states (conformations a-d) corresponding to nodes 134, 90, 316, and 504, and their side-chain contact maps are also shown. Edge weights are labeled. If there are multiple edges between two basins, e.g., basins 1 and 90, only the total edge weight is labeled. The graph had good symmetry: the left part and the right part are related to the exchange of the two strands. (b) Illustration of the reptation motion along pathway 1(3)f134(122)f10(11) through a β-bulge. The strands are labeled as 1 and 2. For simplicity, the residues near C-terminal are not shown. Nodes 1 and 3, nodes 134 and 122, and nodes 10 and 11 are identical if the two strands are considered to be interexchangeable.

diagonal in Figure 2f, node 77 also has Y114-Y114, Y114S115, and S115-Y114 contacts as seen in the parallel aggregates (node 2(4), 10(11) in Figure 2b,2d). Actually, node 77 has the lowest average effective energy due to its compact intrastrand and interstrand side-chain packing. However, its free energy is larger than the other clusters listed above, attributed

to its low configurational entropy. Our observation on that stacking interaction of Y114-Y114 stabilizes the aggregates is consistent with the study of the role of side-chain interactions in the early steps of aggregation of the yeast prion Sup35 segment, in which the parallel β-sheet arrangement is favored due to the stacking of the tyrosine side chains.19

5430 J. Phys. Chem. B, Vol. 111, No. 19, 2007

Li et al.

Figure 4. Examples of pathways that connect the antiparallel (cluster 9) and parallel aggregates (cluster 1). The paths are numbered from I to III, and the conformations along the path are numbered from 1 to 9.

TABLE 1: Summary of Conformation, Interaction, and Free Energy of Most Visited Clusters (Nodes) node

conformna

Nhbondb

ncontactb

∆Gsolc

Eintrad

Einterd

∆Ee

∆Ff

FBg

1 2 9 10 16 77

|, 1 |, 1 anti-|, 0 |, 0 |, 3 anti-|, 0

7.4 4.6 7.7 5.0 5.2 7.4

8.6 9.7 9.6 8.5 7.8 9.8

0.0 -3.6 -0.6 -1.0 -4.8 6.1

-146.0, -145.8 -144.4, -141.8 -144.8, -144.8 -144.4, -143.2 -146.4, -145.2 -141.5, -141.7

-58.4 -67.4 -60.8 -60.9 -52.4 -79.3

0 -3.4 -0.2 1.7 6.2 -12.3

0.0 0.2 0.9 1.1 1.6 2.8

0.0 4.7 7.7 6.9 7.8 7.4

a | denotes parallel β-sheet and anti-| represents antiparallel arrangements. The Arabic numbers denote the number of residues that are offregister and 0 indicates in-register. b The reported number of interstrand main-chain hydrogen bonds (Nhbond) and the interstrand side-chain contact index (ncontact) are the mean values averaged over the conformations in a given cluster. For a given conformation (with 11 residues/chain), ncontact 121 ) ((∑i)1 ni)/(11 × 11)) × 100, where ni ) 1 if there is contact or, otherwise, ni ) 0. c The difference in solvation free energy (calculated with EEF1) between node 1 and the given node. d The effective energy (E) is the potential energy of the protein plus the solvation free energy. Eintra and Einter denote the intrastrand and interstrand effective energy, respectively. The first number of Eintra represents the intrastrand interaction for strand 1, and the second number is for strand 2. The effective energies of clusters 1, 2, 9, 10, 16, and 77 are -350.2, -353.6, -350.4, -348.5, -344.0, and -362.5 kcal/mol, respectively, which are the average values over all of the conformations in the cluster. e The difference in effective energy between node 1 and the given node. f The difference in free energy between node 1 and the given node. The free energy calculations follow eqs 1-3 in the text. g Free energy barrier (FB) from node 1 to the given node. All of the effective energy and free energy are in the unit of kilocalories per mole.

Kinetics. Figure 3a includes a graph showing the paths between important basins such as the parallel in-register aggregates and the parallel off-register arrangements displaced by one residue. We also depict four typical intermediate states (conformation a-d) and the corresponding side-chain contact maps. Again, node 1 and 3, 10 and 11, 122 and 134, 138 and 188, 229 and 316 are identical if the two chains are considered to be indistinguishable. The apparent symmetry in the graph with exchange of the two chains in Figure 3a indicates that the sampling around these most visited states is enough, and the clustering result based on our conformation similarity measure and the algorithm by Daura et al.39 is reasonable. There are two important pathways connecting the parallel offregister aggregates displaced by one residue in basin 1(3) and the parallel in-register arrangements in basin 10(11), one through basin 134(122) and the other through basin 90 (Figure 3a). Overall, both pathways involve a rotation of one chain (Rot1) by 180° around the strand axis (in other words, all of the peptide bonds in the strand flip 180°) and a reptation move of the strand with respect to the other by shift in one residue (Rep1) as described in the pathways of Aβ16-22;14 however, the detailed mechanisms of reptation motions are different. In basin 90, the broad and weak side-chain contact map (Figure 3a (contact map b)) indicates the side-chain packing is loose except that at the termini, and the two chains are much more flexible than those in minima 1, 3, 10, and 11. The strong side-chain contacts of

Y114-Y114 (contact index ) 0.80) and Y114-S115 (contact index ) 0.70) prevent the dimeric dissociation. However, these packing cannot hinder the reptation move of one chain with respect to the other. Almost all of the hydrogen bonds in the conformations of basin 1 are broken (the average number of interstrand main-chain hydrogen bonds of basin 90 is 0.75 vs 7.4 in basin 1); meanwhile, as seen in Figure 3a contact map b, most of the interstrand side chains repack so that the C-terminal adopts the Y114-Y114, Y114-S115, and S115-Y114 sidechain contact pattern as seen in node 10 (Figure 2d), whereas, the N-terminal packing of Y105-Y105, Y105-T106, and T106-Y105 still resembles node 1 (Figure 2a). In general, the reptation motion along the pathway of 1 f 90 f 10 is conducted mainly by adjusting the side-chain contacts. Instead, the route pass basin 134(122) takes an alternative mechanism: a parallel classic (PC) β-bulge44 occurs in A108 on one strand (X position in Figure 3 of ref 44) and A108 (position 1 in Figure 3 of ref 44) and A109 (position 2 in Figure 3 of ref 44) on the other strand; as a result, the reptation of one strand with respect to the other is carried out by forming a β-bulge as schematically illustrated in Figure 3b. The center of basin 134 was used to drawn Figure 3b, in which the N-terminal residues (105-107) of strand 2 have already rotated 180° around the strand axis or their peptide bonds have flipped 180°, whereas the residues after the β-bulge of strand 2 have not rotated, yet; therefore, the N-terminal of this aggregate is in-register while the C-terminal is off-register.

β-Strand Reorganization in TTR(105-115) Aggregates TABLE 2: Max-flow between Nodes a

wc

f(1T3)

1

6525 100% 4049 62.1% 2791 42.8% 2389 36.6% 1784 27.3% 1668 25.6% 1394 21.4% 1230 18.9%

2 6 10 50 80 100 400

f(1T10)

f(1T9)

2784 100% 1526 54.8% 789 28.3% 603 21.7% 236 8.5% 170 6.1% 0 0% 0 0%

981 100% 259 26.4% 36 3.7% 11 1.1% 0 0% 0 0% 0 0% 0 0%

a The max-flows (f) between nodes 1 and 3, between nodes 1 and 10, and between nodes 1 and 9 are listed as a function of the threshold of the edge weight (wc). The edges with their weights smaller than the threshold were excluded in the calculation of max-flow. The percentage is f/f(wc)1) × 100%.

The fact that basin 90 is a common intermediate state during the conformational transitions among minima 1, 3, 10, and 11 also suggests that the distribution of conformations in this basin is very broad. Similar to the pathway of 1 f 90 f 10, the reptation motion during the conformational reorganization from basin 1 through basin 90 to basin 3 is via side-chain rearrangements. We have also noticed another 1 f 3 route through a special β-bulge that includes more than two residues on one strand opposite a single residue on the other strand,44 but its corresponding edge weight is less than 65; therefore, this route is not illustrated in Figure 3a. Another important path from node 1 to node 3 is 1 f (188) f 316 f 504 f 229 f (138) f 3 (Figure 3a). The typical conformations of basin 316 and 504 are shown in Figure 3a conformations c and d. On average, there are five interstrand main-chain hydrogen bonds in the conformations of basin 316, mainly A108(NH)-T106(CO), A108(NH)-A108(CO), and L110(NH)-A108(CO), which resembles the hydrogen bond pattern of basin 1. Additionally, a pair of antiparallel in-register hydrogen bonds (S115(NH)P113(CO)) and new side-chain contacts, Y114 of one chain interacting with S112 and Y114 of the other chain and hydrogen bonding between the side chains of S115 from each chain, are seen near the C-terminal (Figure 3a (contact map c)), implying that the reptation move of one strand with respect to the other starts from the C-terminal. The C-terminal antiparallel in-register hydrogen bonding of S115(NH)-P113(CO) remains in the conformations of basin 504, whereas other interstrand mainchain hydrogen bonds in basin 316 are not present. The C-terminal side-chain contacts of basin 504 resemble those of basin 316, but the middle region in the contact map (Figure 3a (contact map d)) is wider and weaker than that of basin 316 (Figure 3a (contact map c)), suggesting that the side-chain contact is loose in basin 504 relative to basin 316. To estimate the contributions of edges with different weights to the conformational transitions, we computed the max-flow values from the original connected graph multiple times, when each time we only considered edges with weights no smaller than a certain threshold. The results of max-flow between basins as a function of thresholds, as summarized in Table 2, indicate that the contributions from less visited paths are quite significant. For example, among the transitions between basins 1 and 3, 37.9% ((6525 - 4049)/6525 × 100%) of the max-flow are attributed to the edges with weight equal to 1. Thereby, the

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5431 several most visited pathways in fact do not dominate the kinetics. We have noticed that our sampling might not be enough for the edges whose weights are very small, but the conclusion that the most visited paths do not dominate is unlikely to be changed even with more sampling of the less visited paths. The diversity of the pathways is consistent with the so-called ”new view” of protein folding;37 there are many ways to reach the native state, though it is actually an aggregation process studied here. The analytical studies on the statistics of kinetic pathways45,46 are also relevant to this work, in which our observations appear to be the case where the energy landscape is smooth. For pairs of minima with low free energy barriers, the MD simulations have captured large numbers of transitions including the transitions with significant contributions as reflected in the edges with large weights in the connected graph as shown in Figure 3a. As a result, some important statistical properties of the transitions between these clusters can be reliably deduced from the graph. However, for the minima separated by large barriers, such as the barrier between the parallel aggregates (cluster 1) and antiparallel (cluster 9) (barrier ) 7.7 kcal/mol in Table 1), the MD simulations generated fewer transition pathways, all involved edges with weights smaller than 50 (Table 2). It is not statistically reliable to apply the protocol in Figure 3a in the study of these transitions. Thus, instead we show several possible paths directly observed in the MD trajectories. Along pathway I in Figure 4, the two strands initially adopt an antiparallel conformation (cluster 9) and then the route passes through a perpendicular state (Figure 4 (conformation 2)), which was also observed in the study of Aβ16-22.14 The other two pathways connecting the antiparallel to parallel conformations are not reported in the same work of Aβ16-22.14 Along pathway II (Figure 4), the N terminal of one chain forms a turn of π-helix (Figure 4 (conformation 5)), and then this chain rotated 180° to become parallel to the other chain. The hydrogen bonds propagate from the C- to N-terminal (Figure 4 (conformation 6)). In contrast, a β-hairpin is formed by one strand along pathway III, followed by a rotation of 180° along the long axis of the hairpin (Figure 4 (conformation 8)) before transiting to the parallel aggregates. The time interval for the transitions shown in Figure 4 varies from tens of picoseconds to nanoseconds or even longer. Although these routes shown in Figure 4 connect the transitions between parallel and antiparallel aggregates, the basins are also intermediates for the transitions between parallel aggregates with different displacements. Furthermore, these intermediates may go either direction along the routes; thus they may end up forming either parallel or antiparallel aggregates with various patterns of registry. Comparison with Other Simulations and Experiments. Apparently, all of the simulation results depend on the force field and the solvation model. A recent comparison of several implicit solvation models including EEF1 with explicit water simulation on paired helical filament 6 (PHF6) has demonstrated that EEF1 can accurately reproduce the set of local energy minima arising from quenched dynamics simulations with explicit solvent and predict that PHF6 forms extended β-structures in solution, consistent with the observation that PHF6 initiates neurofibrillary tangle formation in patients with Alzheimer’s disease.47 Previous comparison has been done between EEF1, solvent accessible surface area (SASA) implicit solvation model, and explicit water on several proteins.48 Herein, we compare our results using EEF1 with those using the SASA model on the aggregation of the same peptide.22 These two models are consistent with each other in some aspects: (1) both

5432 J. Phys. Chem. B, Vol. 111, No. 19, 2007

Li et al. and Aβ16-22 have also suggested that the kinetically favored (initially formed) aggregates are likely a distribution of strand alignments, and the conformational rearrangement within the aggregates to transit to a thermodynamically more stable state is a general mechanism of aggregation.29-31 Summary

Figure 5. Number of interstrand main-chain hydrogen bonds and sidechain contacts averaged over the 80 lowest free energy minima.

simulations suggest that TTR(105-115) peptide predominantly adopts an extended conformation in the amyloid aggregates, which is supported by the solid-state NMR experiments;6,7 (2) both parallel and antiparallel, as well as in- and off-register, aggregates are possible (for parallel aggregates, off-register conformation is more frequently observed than fully in-register conformation); (3) ∆δtot. in chemical shift of the six typical conformations (in clusters 1, 2, 9, 10, 16, and 77) from the experimental results7 are in the same range as that using the SASA model22 (∆δtot. ) 1.85, 1.85, 1.88, 1.96, 2.12, and 1.95 for clusters 1, 2, 9, 10, 16, and 77, respectively). Note that no experimental restraints were used in our simulations. Meanwhile, the EEF1 and SASA results disagree with each other in three main aspects. First of all, the extended and aligned aggregates are very rare (dimer population ) 1.9%) in the unbiased SASA simulation,22 whereas these conformations are dominate in our simulations. Therefore, the SASA model has observed the prenucleation stage, while our EEF1 simulations have revealed the conformational reorganization within aggregates. Second, for the antiparallel aggregates, in contrast to our results, the SASA model has found that the aggregates displaced by at least one residue is more populated than the in-register conformation.22 Partly, the second disagreement may be due to the difference in the number of chains simulated in our work and theirs. Third, the stacking interactions of Y105-Y105 and Y114-Y114 are rare in the parallel aggregates of the SASA model, whereas we have found these stacking in our simulation, though the Y114-Y114 contact is not present in the global minimum of free energy. To compare the probability of interstrand main-chain hydrogen bonds formed during our simulations and SASA model, we present the average number of interstrand main-chain hydrogen bonds and side-chain contacts in Figure 5. Overall, the hydrogen bond pattern is similar to Figure 5 of ref 22 (except residues 108 and 109) in which the even numbered residues have more chance to form interstrand hydrogen bonds than the odd numbered residues. Since P113 does not have main-chain hydrogen, the dominant hydrogen bond pattern in the aggregated state is that formed between even numbered residues. For the interstrand side-chain contacts, we have found that Y105 and Y114 have the highest probability (Figure 5). The observations that the existence of several free energy minima slightly destabilized with respect to the global minimum and the multiple routes connecting different aggregates without dissociation are very similar to the simulations on the aggregations of Aβ16-22 by a coarse-grained model,14 suggesting these phenomena may be common in short peptide aggregation. The isotope-edited FTIR studies of Syrian hamster PrPC(109-122)

In this work, we have investigated the interstrand rearrangements of TTR(105-115) peptide within dimeric aggregates by means of molecular dynamics with EEF1 implicit solvation model. First, we have found that both parallel and antiparallel aggregates are populated. Second, parallel aggregates with different alignment patterns are separated by nonnegligible free energy barriers. Third, multiroutes exist in the interstrand conformational reorganization. Most visited routes do not dominate the kinetics, while less visited routes contribute a little each but they are numerous and their total contributions are actually dominant. Fourth, there are various kinds of reptation motions, including those through a β-bulge, side-chain-aided reptation, and flipping or rotation of a hairpin formed by one strand. Acknowledgment. This work is supported by NIH (Grant 1R15 AG025023-01) and NSF Major Research Instrumentation (Grant DBI-0320875). The molecular graphics images were generated using the Chimera49 package from the Computer Graphics Laboratory, University of California, San Francisco, CA (supported by NIH Grant P41 RR-01081). We think Prof. Ping Xuan for helpful discussions on clustering methods. We are very grateful to Prof. David Wales for the programs for creating and manipulating disconnectivity graphs. We also acknowledge the National Center of Supercomputing Applications and the Scientific Computing and Visualization group at Boston University for providing part of the computational resources. Supporting Information Available: Text discussing how free energy minima and barriers depend on the clustering method and threshold and figures showing minimal cut trees obtained with various RMSD thresholds and an illustration of two clusters with a short RMSD distance but separated by a high free energy barrier. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Caughey, B.; Lansbury, P. T. Annu. ReV. Neurosci. 2003, 26, 267. (2) Hardy, J.; Selkoe, D. J. Science 2002, 297, 353. (3) Kayed, R.; Head, E.; Thompson, J. L.; McIntire, T. M.; Milton, S. C.; Cotman, C. W.; Glabe, C. G. Science 2003, 300, 486. (4) Reixach, N.; Deechongkit, S.; Jiang, X.; Kelly, J. W.; Buxbaum, J. N. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 2817. (5) Hamilton, J. A.; Benson, M. D. Cell. Mol. Life Sci. 2001, 58, 1491. (6) Jaroniec, C. P.; MacPhee, C. E.; Bajaj, V. S.; McMahon, M. T.; Dobson, C. M.; Griffin, R. G. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 711. (7) Jaroniec, C. P.; MacPhee, C. E.; Astrof, N. S.; Dobson, C. M.; Griffin, R. G. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 16748. (8) Ding, F.; Borreguero, J. M.; Buldyrey, S. V.; Stanley, H. E.; Dokholyan, N. V. Proteins 2003, 53, 220. (9) Peng, S.; Ding, F.; Urbanc, B.; Buldyrev, S. V.; Cruz, L.; Stanley, H. E.; Dokholyan, N. V. Phys. ReV. E: Stat., Nonlinear, Soft Matter Phys. 2004, 69, 041908. (10) Borreguero, J. M.; Urbanc, B.; Lazo, N. D.; Buldyrev, S. V.; Teplow, D. B.; Stanley, H. E. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6015. (11) Nguyen, H. D.; Hall, C. K. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 16180. (12) Nguyen, H. D.; Hall, C. K. Biophys. J. 2004, 87, 4122. (13) Nguyen, H. D.; Hall, C. K. J. Biol. Chem. 2005, 280, 9074.

β-Strand Reorganization in TTR(105-115) Aggregates (14) Santini, S.; Wei, G.; Mousseau, N.; Derreumaux, P. Structure (London) 2004, 12, 1245. (15) Santini, S.; Mousseau, N.; Derreumaux, P. J. Am. Chem. Soc. 2004, 126, 11509. (16) Favrin, G.; Irback, A.; Mohanty, S. Biophys. J. 2004, 87, 3657. (17) Fawzi, N. L.; Chubukov, V.; Clark, L. A.; Brown, S.; Head-Gordon, T. Protein Sci. 2005, 14, 993. (18) Hwang, W.; Zhang, S.; Kamm, R. D.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 12916. (19) Gsponer, J.; Haberthur, U.; Caflisch, A. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 5154. (20) Zanuy, D.; Porat, Y.; Gazit, E.; Nussinov, R. Structure (London) 2004, 12, 439. (21) Tsai, H. H.; Reches, M.; Tsai, C. J.; Gunasekaran, K.; Gazit, E.; Nussinov, R. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 8174. (22) Paci, E.; Gsponer, J.; Salvatella, X.; Vendruscolo, M. J. Mol. Biol. 2004, 340, 555. (23) Wu, C.; Lei, H.; Duan, Y. J. Am. Chem. Soc. 2005, 127, 13530. (24) Wu, C.; Lei, H.; Duan, Y. Biophys. J. 2004, 87, 3000. (25) Klimov, D. K.; Thirumalai, D. Structure (London) 2003, 11, 295. (26) Lomakin, A.; Teplow, D. B.; Kirschner, D. A.; Benedek, G. B. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 7942. (27) Perutz, M. F.; Windle, A. H. Nature 2001, 412, 143. (28) van Gestel, J.; de Leeuw, S. W. Biophys. J. 2006, 90, 3134. (29) Petty, S. A.; Decatur, S. M. J. Am. Chem. Soc. 2005, 127, 13488. (30) Petty, S. A.; Decatur, S. M. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 14272. (31) Silva, R. A.; Barber-Armstrong, W.; Decatur, S. M. J. Am. Chem. Soc. 2003, 125, 13674.

J. Phys. Chem. B, Vol. 111, No. 19, 2007 5433 (32) Krivov, S. V.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 14766. (33) Rao, F.; Caflisch, A. J. Mol. Biol. 2004, 342, 299. (34) Caflisch, A. Curr. Opin. Struct. Biol. 2006, 16, 71. (35) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. Journal Comput. Chem. 1983, 4, 187. (36) Lazaridis, T.; Karplus, M. Proteins 1999, 35, 133. (37) Lazaridis, T.; Karplus, M. Science 1997, 278, 1928. (38) Ryckaert, J. P.; Ciccotti, G.; Berebdseb, H. J. C. J. Comput. Phys. 1997, 23, 327. (39) Daura, X.; Gademann, K.; Jaun, B.; Seebach, D.; van Gunsteren, W. F.; Mark, A. E. Angew. Chem., Int. Ed. 1999, 38, 236. (40) Cormen, T. H.; Leiserson, C. E.; Rivest, R. L. Introduction to Algorithms; MIT Press, McGraw-Hill: Cambridge, MA, New York, 1990. (41) Gomory, R. E.; Hu, T. C. SIAM J. Appl. Math. 1961, 9, 551. (42) Krivov, S. V.; Karplus, M. J. Chem. Phys. 2002, 117, 10894. (43) Neal, S.; Nip, A. M.; Zhang, H.; Wishart, D. S. J. Biomol. NMR 2003, 26, 215. (44) Chan, A. W.; Hutchinson, E. G.; Harris, D.; Thornton, J. M. Protein Sci. 1993, 2, 1574. (45) Wang, J. J. Chem. Phys. 2003, 118, 952. (46) Wang, J.; Onuchic, J. N.; Wolynes, P. G. Phys. ReV. Lett. 1996, 76, 4861. (47) Huang, A.; Stultz, C. M. Biophys. J. 2007, 92, 34. (48) Ferrara, P.; Apostolakis, J.; Caflisch, A. Proteins 2002, 46, 24. (49) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. J. Comput. Chem. 2004, 25, 1605.