Single Point Mutation Alters the Microstate Dynamics of Amyloid β

May 10, 2013 - A number of single point mutations at residues A21, E22, D23, and M35 ... E22K, and E22Q) and D23N generate more hub-like microstates t...
1 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Single Point Mutation Alters the Microstate Dynamics of Amyloid β‑Protein Aβ42 as Revealed by Dihedral Dynamics Analyses Liang Xu,*,† Shengsheng Shan,† and Xicheng Wang‡ †

School of Chemistry, Dalian University of Technology, Dalian 116023, China Department of Engineering Mechanics, State Key Laboratory of Structural Analyses for Industrial Equipment, Dalian University of Technology, Dalian 116023, China

‡

S Supporting Information *

ABSTRACT: The aggregation of amyloid β-protein (Aβ) has been associated with the pathogenesis of Alzheimer’s disease. A number of single point mutations at residues A21, E22, D23, and M35 have been identified to show increased or decreased aggregation tendency. Although the effects of point mutations on the structural properties of Aβ peptides have been intensively studied, how single point mutation affects the kinetics of Aβ remains unknown. In this work, dihedral dynamics analyses, which combine dihedral principle component analysis (dPCA), potential of mean force (PMF) calculations, and Markov state models (MSMs), were proposed to elucidate the different global free energy landscapes (FELs), the PMF of individual dihedral angle, and microstates/macrostates for a number of Aβ42 mutants (Flemish A21G, Arctic E22G, Italian E22K, Dutch E22Q, Iowa D23N, Japanese E22Δ, and M35 oxidation Met35OX). Our simulation results show that one point mutation is sufficient to change the rugged FEL of Aβ42 by altering the energy barriers around basins. This alteration was also observed in the potential of each dihedral angle to varying degrees, although most minima of PMF do not shift. MSMs further reveal that E22 mutants (E22Δ, E22G, E22K, and E22Q) and D23N generate more hub-like microstates than wild type Aβ42, thus creating diverse alternative pathways for conformational transitions and increasing subsequent aggregation. In contrast, transitions are more preferred within the same microstate of A21G and Met35OX. Mapping MSM to FEL suggests that transitions between different sets of microstates are kinetically feasible but thermodynamically difficult.

1. INTRODUCTION Amyloid β (Aβ) peptides consisting of 39−43 amino acids have long been associated with the pathogenesis of Alzheimer’s disease (AD).1−7 Because Aβ peptides are intrinsically disordered proteins (IDPs), the lack of specific secondary and tertiary structures under physiological conditions provides great challenges in probing the structure−function relationships.2,6 Three NMR-based Aβ monomers (Aβ42/Aβ40) are available (PDB ID: 1IYT,8 1Z0Q,9 and 2LFM10); however, these structures were resolved either in nonpolar or varying polar aqueous environment, showing different propensities of αhelical structures. Numerous experimental and computational simulation studies have indicated that Aβ monomer primarily adopts a collapsed coil structure in aqueous solution close to physiological temperature and pH.2,11−20 Recent studies have pointed out that small Aβ oligomers are responsible for the cytotoxicity of Aβ in brain cells.21−26 The aggregation of Aβ monomers into oligomeric species is accompanied by conformational changes. Several pathogenic point mutations of Aβ including A21G (Flemish),27 E22G (Arctic),28 E22K (Italian),29 E22Q (Dutch),30 D23N (Iowa),31 and E22Δ (Japanese)32 have © 2013 American Chemical Society

been identified. These variants have been found to show different effects on the aggregation of Aβ.33−38 The A21G variant shows a slower aggregation rate than wild type Aβ. In contrast, the E22X (X = G, K, and Q) mutations show enhanced aggregation propensity and enhanced neurotoxicity but less fibril formation. The E22Δ has also been reported to accelerate oligomerization and synaptotoxicity but no cytotoxicity was found. Similar to E22X variants, the D23N exhibits accelerated aggregation tendency and thus is more neurotoxic. In addition, the side chain of Met35 has been demonstrated to be oxidized to sulfoxide and the oxidized species (Met35OX) has also been found in post-mortem senile plaques,39−41 a hallmark of AD. Further studies have shown that the Met35OX species prevent formation of protofibril and thus reduce the aggregation of Aβ42(Met35OX) compared to wild type Aβ42.39−42 These observations clearly show that the oligomerization is very sensitive to single amino acid substitution, leading to different pathogenicity of Aβ related to AD. Received: April 3, 2013 Revised: April 30, 2013 Published: May 10, 2013 6206

dx.doi.org/10.1021/jp403288b | J. Phys. Chem. B 2013, 117, 6206−6216

The Journal of Physical Chemistry B

Article

mutation does not shift most of the dihedral energy basins, it alters each dihedral energy surface to varying extents. Moreover, one point mutation significantly changes the internal dynamics between microstates, as revealed by MSMs. The different effects of point mutants on the aggregation tendency can be related to the increase or decrease of hub-like microstates compared to wild type Aβ42.

Elucidating the clear scenario of the Aβ misfolding process needs the combination of experiments and theory studies, with molecular dynamics (MD) simulations in particular yielding detailed structural information at an atomic level. Extensive MD simulations, especially replica-exchange and discrete MD simulations at both all-atom and coarse-grained level, have been performed to probe the effects of those familial AD mutations on the structure and dynamics of Aβ monomers and small oligomers.34,38,43−57 The simulation time scale varies from nanoseconds (ns) to hundreds of microseconds (μs), providing valuable insights into understanding the aggregation mechanism of Aβ. Although most of those simulations were validated by comparison with experimental observables, such as chemical shifts, J-coupling constants, residual dipolar couplings, and NOE values, consistent results on the structure propensity of Aβ seem not to be expected.12 For example, the secondary structure preferences are quite different between the 100 ns REMD simulations with the generalized Born implicit solvation model (100 ns trajectory used for analysis, and the total simulation time is 3.52 μs)17 and long-time-scale (>200 μs) classical MD simulations using the explicit TIP3P solvent model.52 The β-strand population of Aβ42 obtained from the later simulations are much less abundant (0−2.5%) compared to the results of REMD simulations (0−14%). On the other hand, the dynamics of Aβ are usually characterized by different clustering methods or statistical learning techniques in order to use the most discriminative conformations to represent the whole heterogeneous conformational space.14,16,17 In this way, the investigation of conformational conversions becomes amenable in a reduced space. However, such dynamics of Aβ are basically static because there is little information on the kinetics among different ensembles. Consequently, little is known about how single point mutation alters the pathway of oligomerization and the corresponding kinetics of conformational transitions. The Markov state model (MSM) provides a convenient way to model kinetic networks between different conformational states.58−61 In the same state, the free energy basin is quite flat such that conformations can interconvert easily, whereas different states are separated by large energy basins. The transition probability matrix can be calculated to group the kinetically relevant microstates into more coarse grained macrostates. The MSM approach has been intensively applied to analyze the folding/unfolding pathways of proteins such as the FiP35 WW domain,62 λ repressor protein,63 NTL9,64 and Wimley−White pentapeptides.65 Although MSM has been used in previous studies to compute the ensemble-averaged secondary structures of Aβ42 and its familial mutations (E22K, E22G, E22Q, and D23N),51,52 it remains unknown if Aβ42 and mutants behave as the Markov process. In this study, we systematically analyze the dihedral dynamics properties of Aβ42 and its single point variations (A21G, E22Δ, E22G, E22K, E22Q, D23N, Met35OX). On the basis of trajectories generated from REMD simulations, the free energy surface of each system is first constructed using dihedral principle component analysis (dPCA)66−68 to give a general picture of the effects of mutations on the energy landscapes of Aβ42. A closer examination of dihedral dynamics is carried out by calculation of the potential of mean force (PMF) for each dihedral angle.69 In order to build a MSM, conformations of each system are clustered according to ϕ and ψ dihedral angle pairs. Our results show that one point mutation is sufficient to affect the free energy landscape of Aβ42. Although one point

2. METHODS 2.1. REMD Simulations. In our previous studies, we have simulated Aβ42 with the GB implicit solvent model.15 Convergence of REMD simulations was checked by calculation of cumulative helix contents and comparisons with experimental chemical shifts.15 The same simulation protocol was used here. Briefly, the initial Aβ42 structure was taken from the NMR structure (PDB ID: 1Z0Q).9 The Amber ff99SB force field70 and the Onufriev−Bashford−Case GB implicit solvent model71 were applied to represent Aβ42 (including Aβ42 mutants) and solvent effects, respectively. The Amber force field parameters for the oxidized Met35 were taken from previous studies.47 Different from our previous studies, no cutoff was used for long-range interactions. The temperature ranges from 280 to 400 K and eight replicas were used for each system, resulting in an acceptance ratio of about 20%. The total simulation time was 800 ns for each system (with each replica being 100 ns). Previous REMD simulations have shown that 60 ns was sufficient to reach convergence.15,16 Without using a cutoff, it was found that converged simulation could be obtained after 40 ns (Figure S1, Supporting Information), allowing us to maximize the use of simulation trajectory to build reliable MSMs. Convergence was further validated by calculation of the radius of gyration of Aβ42 and its mutants at different time intervals (Figure S2, Supporting Information). All REMD simulations were performed using the AMBER12 software package.72 Trajectories collected at a physiological temperature of 310 K were used for data analysis. 2.2. Dihedral Dynamics Analyses. Compared to principle component analysis (PCA), dihedral PCA (dPCA) has been shown to be able to separate internal and overall motions of a biomolecule using backbone dihedral angles in PCA.66,67 Thus, dPCA can characterize more complex features of the free energy landscape of a flexible protein. We employ dPCA to construct the free energy surfaces of wild type Aβ42 and seven point variations (A21G, E22Δ, E22G, E22K, E22Q, D23N, Met35OX), aiming to show the general effects of the point mutations on the structure and energy properties of Aβ42. A closer examination of dihedral dynamics was carried out by calculation of the potential of mean force (PMF) for individual dihedral angle. Instead of using ϕ and ψ dihedral angle pairs as we have done in dPCA, we attempt to characterize the dihedral dynamics in an efficient and coarse way. In line with MARTINI coarse-grained models,73−75 a dihedral angle is defined by four consecutive beads, with each bead representing the center of mass of one residue. The distributions of all dihedral angles were calculated for each system. Then, the corresponding PMF of each dihedral angle was calculated from the probability distributions according to the Boltzmann inversion method, i.e., PMF = −kBT ln(Pi), where kB is the Boltzmann constant, T is the temperature, and Pi is the probability distribution of the ith dihedral angle.69 2.3. MSM Analysis. In most previous studies, conformations were clustered according to atom RMSD relative to the native protein structure.63−65 Thousands of microstates and a 6207

dx.doi.org/10.1021/jp403288b | J. Phys. Chem. B 2013, 117, 6206−6216

The Journal of Physical Chemistry B

Article

respectively. The high propensity of α-helix structure is different from numerous experimental and simulation results, which indicate that Aβ42 predominantly adopts a coil structure (60−80%), β-strand (10−20%), and negligible α-helix (80% of all conformations of Aβ42 and mutations, respectively (Figure 1). Such results

Figure 1. Clustering results based on the RMSD of backbone atoms. The Daura method was used with a cutoff of 3.0 Å.

suggest that the number of important macrostates cannot be too much. The Ward clustering method implemented in the MSMBuilder2 package76 was applied to cluster conformations of each system into different numbers of microstates. Then, the transition probability matrix was calculated to determine the implied time scales corresponding to the time scale required for transition between different microstates. The Perron cluster cluster analysis (PCCA) algorithm was used to group kinetically related microstates into different macrostates. MSMExplorer77 was used to visualize the results of MSM for each system.

3. RESULTS AND DISCUSSION 3.1. Structure Properties of Aβ42 and Its Mutants. The average radius of gyration of Aβ42 is 1.01 ± 0.07 nm (Figure S2, Supporting Information), consistent with recent experimental results (0.9 ± 0.1 nm)78 and several simulation studies.79−83 However, a larger radius of gyration (1.2−1.3 nm) has also been reported in other studies under different conditions.12,16,83 The radius of gyration does not change significantly for Aβ42 mutants, suggesting that single amino acid substitution may not alter the general conformations of Aβ42. The distribution of secondary structures of Aβ42 and its mutants was calculated and shown in Figure S3 (Supporting Information). We note that the average α-helix contents are 26.3, 31.5, 36.7, 23.6, 30.7, 38.4, 38.6, and 42.3% for Aβ42, A21G, D23N, E22Δ, E22G, E22K, E22Q, and Met35OX, 6208

dx.doi.org/10.1021/jp403288b | J. Phys. Chem. B 2013, 117, 6206−6216

The Journal of Physical Chemistry B

Article

Figure 2. Conformational free energy surfaces for Aβ42 and variants. Each was constructed by projecting its conformations onto the first two dihedral principle components (dPC1 and dPC2). The free energy values (in kcal/mol) were obtained by ΔG = −kBT[ln Pi − ln Pmax], where Pi and Pmax are the probability distributions calculated from a histogram of individual REMD simulation trajectory. ln Pi − ln Pmax was used to ensure ΔG = 0 for the lowest free energy points. The highlighted four regions in the free energy surface of Met35OX roughly correspond to four macrostates/ microstates identified in Met35OX MSM. The C- and N-terminals of each structure are represented as orange and cyan van der Waals spheres, respectively. The mutated residue is shown as a van der Waals sphere as well.

4A). Also note that another basin (−120°) was observed for E22Δ with comparable PMF (∼11 kJ/mol). Interestingly, the C terminal, which is far away from the mutation points, displays complex dihedral dynamics, suggesting that these point mutations have distinct long-rang effects on the dynamics of Aβ42 (Figure 4D). Moreover, besides the global basin, there are numerous local maxima and minima distributed along the dihedral potential of each peptide (Figure 4 and Figure S4, Supporting Information), suggesting that one point mutation could influence the misfolding pathways of Aβ42 dramatically.

average values of PCC for the same dihedral angle of Aβ42 mutants vary from 0.14 to 0.94, implying that the effects of point mutation on the dynamics of Aβ42 are different with respect to the distance from the mutated residue. Inspection of individual dihedral angle PMF (Figure S4, Supporting Information) indeed shows that, in the regions consisting of Leu17 to Asn27 and Lys29 to Gly33, the global minima of PMF do not shift (Figure 4B and C). Another region, in particular, from Asp1 to Lys16, known as the disordered segment of Aβ42, also demonstrates that the PMF values almost fall into the same global minima (30°) for residues from Val12 to Lys16 (Figure 6209

dx.doi.org/10.1021/jp403288b | J. Phys. Chem. B 2013, 117, 6206−6216

The Journal of Physical Chemistry B

Article

G21−E22 (D19) and E22−D23−V24−G25 (D22), are significantly affected (Figure 3). Moreover, the potential for dihedrals from D19 to D22 becomes lower and the surfaces are much flatter than Aβ42 (Figure 4B and Figure S2, Supporting Information), indicating that this region becomes more flexible. The same trend can also be observed for potentials of dihedral angles from D24 (V24−G25−S26−N27) to D28 (K28−G29− A30−I31). This finding is consistent with previous MD studies which show that A21G may decrease the aggregation propensity by increasing the flexibility of the central hydrophobic region (residues 17−21) and affecting the formation of salt bridges involving residues Glu22, Asp23, and Lys28.43 As discussed above, the average PCC values of E22X (X = Δ, G, K, and Q), D23N, and Met35OX vary from 0.60 to 0.72, indicating fair similarity of overall dynamics between Aβ42 variants and wild type. We observe that these mutations have minimal effects on the location of the global potential basin for most dihedral angles; however, the potential surfaces become quite rugged and vary with distinct mutants to different degrees (Figure 4 and Figure S2, Supporting Information). For example, E22Δ greatly influences the dynamics of the Cterminus of Aβ42 by increasing the potential of D37 to D39 in particular, and consequently increasing the rigidity of E22Δ. On the other hand, E22X and D23N exhibit quite distinctive behavior in the region of residues Ala30 to Ala42. A discernable difference can also be observed in other regions including Leu17

Figure 3. Pearson correlation coefficient (PCC) calculated for the PMF of the individual dihedral angle of Aβ42 mutants relative to wild type Aβ42.

The inverse correlation between Aβ42 and A21G for the dihedral angle consisting of residues Asp7 to Tyr10 can be primarily attributed to the fact that the N-terminus of Aβ is unstructured. Such negative correlation is also observed in E22K and E22Q mutants (Figure 3). Thus, lower PCC values are expected in this region. Although a high value (0.88) of PCC was observed for dihedral angle G21−E22−D23−V24, the potentials of adjacent dihedral angles, especially F19−F20−

Figure 4. Representative PMF for dihedral angles of Aβ42 and its mutants. The dihedral angle was defined by four consecutive coarse-grained beads with each bead representing the center of mass of each residue. 6210

dx.doi.org/10.1021/jp403288b | J. Phys. Chem. B 2013, 117, 6206−6216

The Journal of Physical Chemistry B

Article

Figure 5. Markov state model constructed for Aβ42. The most populated microstates with their equilibrium populations >1% are shown. The sizes are proportional to the corresponding equilibrium populations. Microstates belonging to the same macrostate are represented using the same color, whereas microstates without coloring indicate different macrostates (containing only one microstate). Transition probabilities (>1%) between microstates are also listed. The C- and N-terminals of each structure are represented as orange and cyan van der Waals spheres, respectively. The mutated residue is shown as a van der Waals sphere as well.

Figure 6. Markov state model constructed for A21G. The representations are the same as those shown in Figure 5.

to Ala21. From a structure point of view, these regions (residues 17−21 and 30−42) have higher propensities to form β-strands and thus are crucial to the aggregation of Aβ. Previous MD simulations of Aβ(21−30) assumed that E22X (X = G, K, and Q) mutants may affect long-range interactions outside the Aβ(21−30) region.34 Different aggregation mechanisms between E22X and D23N have also been suggested on the basis of structural comparisons.34 With respect to the salt bridge involving Glu22, Asp23, and Lys28, the corresponding dihedrals (D19 to D23, and D28 in particular, Figure S4, Supporting Information) display diverse dynamics, indicating different effects of mutations at residues E22 and D23 on the formation of a salt bridge. Compared to other mutants,

Met35OX is found to have higher potential energies for dihedrals consisting of residues Glu22 to Lys28 (Figure S4, Supporting Information), suggesting increased rigidity in this region. Overall, our results here demonstrate the capability of dihedral PMF to characterize the relevant dynamics of Aβ42 and its mutants, and reflect the corresponding structural properties as well. 3.4. MSM Analyses. Although dPCA and PMF of dihedral angles have shown that point mutation of Aβ42 has a distinct consequence for the overall dynamics behavior and subsequent aggregation mechanism, little is known about the kinetics of Aβ42 and mutants. In order to identify the most populated states and transitions between them, we construct the Markov 6211

dx.doi.org/10.1021/jp403288b | J. Phys. Chem. B 2013, 117, 6206−6216

The Journal of Physical Chemistry B

Article

Figure 7. Markov state model constructed for E22Q. The representations are the same as those shown in Figure 5.

Figure 8. Markov state model constructed for Met35OX. The representations are the same as those shown in Figure 5.

6212

dx.doi.org/10.1021/jp403288b | J. Phys. Chem. B 2013, 117, 6206−6216

The Journal of Physical Chemistry B

Article

probabilities. As expected, the largest state (S5) corresponds to the central hub. To elucidate the enhanced aggregation propensity of E22X (X = Δ, G, K, and Q) and D23N mutants, their structural properties have been intensively characterized using computational simulations. For instance, REMD simulations of Aβ(21− 30) suggested that E22 mutants destabilized this bend structure.34 The decreased propensity of turn structure at Ala21-Ala30 was also observed for full-length E22Δ from REMD simulations.53 Similar, discrete MD simulations suggested that E22G increased the propensity of the β-strand, and destabilized the main folding region of Ala21-Ala30.44 The oligomer structure of E22G more resembled Aβ42 and is thus potentially more toxic.46 Monte Carlo simulations also revealed that E22G reduced the propensity of the turn conformation at residues Asp23-Ser26, and thus accelerated its aggregation.54 On the other hand, MD simulations of Aβ(10−35) in aqueous solution indicated that the increased aggregation rate of E22Q was due to the differing charge state from wild type45 and enhanced hydrophobic solvation.55 REMD simulations also suggested that E22Q displayed an increased rigidity at its N-terminus.49 Moreover, recent long time MD simulations found that E22G, E22K, E22Q, and D23N decreased the propensity of helix in residues Gly33-Val36, but E22K and E22Q increased helix formation in residues Phe20-Val24.51,52 In addition, MD simulations suggested the higher internal stability of E22X oligomers than the wild type species.38 Coarse-grained MD simulations of D23N(Aβ40) monomer48 and dimer50 indicated that D23N(Aβ40) had higher β-strand propensities at residues Ala30-Val40. In contrast, all-atom MD simulations implied that A21G destabilized the β-strand structure, increased the flexibility of residues Leu17-Ala21, and affected the salt bridges between Glu22/Asp23 and Lys28 of Aβ40 and Aβ42 dimers differently.43 The oxidized Met35 has been suggested to stabilize the hydrophobic turn structures of residues Val24-Asn27 56 but destabilize cross-β fibrils according to MD simulations.47 Experimental studies revealed that oxidation of Met 35 significantly impeded the aggregation rate of Aβ42 mainly by increasing the backbone mobility at the C-terminal,40−42 and prevented the formation of a paranucleus and consequently produced indistinguishable polymorphic oligomers.39 Mapping the microstate MSMs to different basins in the freeenergy surface determined in terms of dPCA, we can find that each basin corresponds to one populated microstate, and one macrostate corresponds to adjacent basins separated by shallow barriers. Conformations within these basins are believed to have a higher transition probability, and therefore are kinetically relevant. For example, dPCA identifies five large basins in the free-energy surface of Met35OX, with four of them separated by rather high energy barriers whereas connected to the central basin by relatively low barriers (Figure 2). This finding is consistent with the corresponding MSM (Figure 8), which shows that, in addition to self-transitions within each microstate, transitions between the largest central state and other states are also preferred. This mapping strategy can also be applied to other systems, allowing us to determine the transition possibility (thermodynamics) from dPCA and transition probability from MSM (kinetics). Taking our dihedral dynamics analyses together, an alternative mechanism emerges from our results to explain the effects of Aβ42 mutants on the different propensity of aggregation. We summarize that all mutants alter the overall dynamics of Aβ42 differently, making corresponding free

state model (MSM) in terms of dihedral angles for each system. The implied time scales of the microstate MSM were examined and showed that a lag time of 3 ns could lead to Markovian behavior (Figure S5, Supporting Information). The slowest implied time falls in the range 10−20 ns (Figure S5, Supporting Information), showing that our simulation time is adequate to sample the kinetics of these peptides between different sets of microstates. NMR experiments with Lipari−Szabo model-free analysis have revealed that the internal motion of most residues has a relaxation time of about 4 ns for Aβ40 and Aβ42,93 confirming that our simulations are adequate to describe the global conformational dynamics on a ps to ns time scale. Lowresolution MSMs were created to facilitate interpretation of our models; i.e., the highly populated microstates with equilibrium probability larger than 1% and transition probability larger than 1% were analyzed. Figures 5−8 show MSMs for Aβ42, A21G, E22Q, and Met35OX, respectively. MSMs for D23N, E22Δ, E22G, and E22K are provided in Figures S6−S9 (Supporting Information), respectively. These distinct MSMs clearly show that point mutation can significantly alter the kinetic behavior of Aβ42. Since Aβ42 is intrinsically disordered, the so-called native state along the transition network is not expected. We assume that, if there are more nodes (microstates or macrostates) and higher transition probability between nodes, there should be more diverse misfolding pathways for the aggregation of Aβ42. Figure 5 shows the top 10 microstate MSMs built for Aβ42, which belong to five macrostate MSMs. The transition probabilities between microstates are also listed. It is observed that conformational transitions predominantly take place between the top four microstates which are further clustered into one kinetically relevant macrostate. Self-transitions in both highly (e.g., S1, S2, and S3) and less populated states (e.g., S6 and S7) are observed. In addition, two moderately populated states (no label) with relatively low transition probabilities (