α-helices and β-sheets are the two main secondary structures stabilized by backbone hydrogen bonds in proteins. The propensity of each residue to form an α-helix or a β-sheet depends on the amino-acid sequence and on the fold (Chou and Fasman, 1974; Smith et al., 1994; Pace and Scholtz, 1998; Bhattacharjee and Biswas, 2010). The formation of long β-sheets differs significantly from α-helices as the former necessarily involves long-distance contacts in the sequence and depends on the mean curvature of the fold (Koh et al., 2006; Bhattacharjee and Biswas, 2010). The content of these two secondary structure elements is a criterion for classification of protein native structures and characterization of protein folding kinetics and pathways (Konagurthu et al., 2020). In intrinsically disordered proteins (IDPs) (Wright and Dyson, 1999; Uversky et al., 2000; Uversky, 2019; Deiana et al., 2019), α-helices and β-sheets are metastable secondary structures. As for folded proteins, a natural extension of protein classification of IDP can be based on the content of these two secondary structural elements in a statistical sense. As an IDP has no native state, the relative content of α-helices (α) and of β-sheets (β) possibly defines a fingerprint of each conformational state in a two-dimensional (α, β) propensity map (Ullman et al., 2011). Here, we aim to build such an (α, β) effective free-energy map for α-synuclein (α-syn), wild-type (WT), and A30P, A53T, and E46K mutants from first principles by using coarse-grained molecular dynamics (MD) simulations with the UNRES (UNited RESidue) force field (Liwo et al., 2001; Maisuradze et al., 2010).

α-syn is a 140 amino-acid protein abundant in the brain (Jakes et al., 1994; Mollenhauer et al., 2008). It plays a central role in the onset of the Parkinson disease (PD) and other neurological disorders, named synucleopathies (Soto, 2003; Chiti and Dobson, 2006; Stefanis, 2012; Chiti and Dobson, 2017; Tanudjojo et al., 2021). In these diseases, α-syn is found in high concentration, as filaments, in spherical inclusions (Lewy bodies) located in the neuron cytoplasm (Spillantini et al., 1997; Breydo et al., 2012; Lashuel, 2020; Trinkaus et al., 2021). Familial cases of PD are induced either by an overexpression of WT α-syn due to α-syn gene triplication or by pathogenic mutations in α-syn gene corresponding to single amino-acid substitution, namely, H50Q, G51D, A53E, A30P, A53T, and E46K (Polymeropoulos et al., 1997; Krüger et al., 1998; Zarranz et al., 2004; Fuchs et al., 2008; Appel-Cresswell et al., 2013; Pasanen et al., 2014; Petrucci et al., 2016). These inherited forms of PD have phenotypes that are similar to the sporadic PD, but patients with the A53T mutant exhibit an early onset of PD (Petrucci et al., 2016). The role of α-syn in PD development and the mechanisms of its aggregation and of the aggregate propagation are still debatable (Henrich et al., 2020; Sang et al., 2021) and involve α-syn-lipid interactions (Galvagnion et al., 2015; Suzuki et al., 2018; Antonschmidt et al., 2021), a liquid-liquid phase transition (Ray et al., 2020), and a secondary nucleation and fibril fragmentation (Cremades et al., 2012; Cascella et al., 2021; Kumari et al., 2021). Moreover, α-syn aggregation might disregulate the mitochondrial function, and increasing the level of free radicals and alterations in this pathway may be involved in the pathogenesis of PD (Hsu et al., 2000; Devi et al., 2008).

From a structural point of view, α-syn is a chameleon protein (Ullman et al., 2011; Mor et al., 2016): it is disordered (IDP) in solution under physiological conditions, in equilibrium with a minor α-helical tetrameric form in the cytoplasm (Bartels et al., 2011; Wang et al., 2011) and α-helical when bounded to a cell membrane (Fusco et al., 2018). When incubated under physiological conditions in vitro, α-syn aggregates into fibrils with polymorphic cross-β-sheet conformations, in which a core of β-strands is aligned perpendicular to the fibril axis forming extended regular β-sheets with different arrangements (Tuttle et al., 2016; Guerrero-Ferreira et al., 2019; Guerrero-Ferreira et al., 2020). In addition to cylindrical fibrils, ribbon aggregates have been also observed (Bousset et al., 2013). The different polymorphs characterized in vitro are believed to mimic the α-syn filament structures in synucleopathies.

Three main regions were identified in the primary sequence of α-syn regarding their role in α-syn conformational dynamics and aggregation. The N-terminal region (residues 1–60) contains a number of imperfect repeats, with the consensus motif KTKEGV, strongly similar to that found in the amphipathic helices and responsible for membrane binding (Perrin et al., 2000). High-resolution NMR structures revealed a broken helix featuring two curved α-helices of α-syn (residues 2—37 and 45—92) bound to micelles (Ulmer et al., 2005), whereas electron paramagnetic resonance characterizations feature an extended helix for α-syn bound to a lipid membrane (first 97 residues) (Cheng et al., 2013). In solution, exposure of this N-terminal region to solvents is correlated to aggregation propensity (Stephens et al., 2020). The central hydrophobic region of α-syn (residues 61—95), called NAC (non-amyloid component), contains a hydrophobic stretch of 12 residues 71VTGVTAVAQKTV82 necessary for the aggregation (Giasson et al., 2001). Fibrils of the NAC region are also found in plaques of Alzheimer’s disease (Uéda et al., 1993). The flexible acidic C-terminal region regulates fibril formation in vitro seeding experiments (Murray et al., 2003) and plays a role in the secondary nucleation process of amyloids via electrostatic interactions with the lysine-rich N-tail (Kumari et al., 2021). Transient, long-range interactions between the negatively charged C-terminus (residues 120—140) with the positively charged N-terminus and NAC (residues 30—100) were observed by paramagnetic relaxation enhancement (PRE) for WT α-syn (Dedmon et al., 2005). It is hypothesized that these transient interactions are responsible for the larger compactness of α-syn compared to a disordered chain (Dedmon et al., 2005). Structures with long-range contacts between the N-terminal and C-terminal represent about 14% of the conformational ensemble, and in a significant fraction of these structures, residues 68—78 (NAC) are exposed to a solvent, in contrast to the hypothesis that such long-range transient interactions prevent aggregation (Ullman et al., 2011).

Mutations affect the arrangement and growth of the fibrils in vitro (Guerrero-Ferreira et al., 2020). Compared to WT, A53T and E46K mutations aggregate faster, whereas A30P aggregates more slowly than WT (Tosatto et al., 2015; Stephens et al., 2020). In addition, WT, A53T, E46K, and A30P differ in the formation of the different early oligomeric moieties (Tosatto et al., 2015). The mutants A30P and A53T show a greater propensity to form non-fibrillar aggregates than WT (Li et al., 2002), and A53T promotes seeded aggregation in human neurons (Tanudjojo et al., 2021). The rate of lipid-induced aggregation and secondary nucleation have been found to differ by multiple orders of magnitude depending on which missense mutation is involved (Flagmeier et al., 2016). NMR studies have shown an increase flexibility of α-syn in nanosecond–microsecond time scales and a reduction of contacts between C- and N-terminals in mutants (Bertoncini et al., 2005). All single amino-acid substitutions have thus both a kinetic and structural effect on the formation of oligomeric structures. The purpose of the present theoretical study is the characterization of the differences between the conformational ensembles of monomeric WT and mutants, which may contribute to our understanding of the early steps of the aggregation process in solution.

The majority of monomeric WT α-syn conformations have no secondary structures (Ullman et al., 2011). The average helical and β-strand contents of WT α-syn deduced from restrained MD of α-syn fragments and NMR data are about 3 and 11%, respectively (Ullman et al., 2011). The maximum fraction of residues with helical and β-sheet secondary structures was found to be 20% (28 residues) and 28% (39 residues), respectively (Ullman et al., 2011). Little is known on how the missense mutations modify the propensity of α-syn to form secondary structures. As the α-helical/β-sheet equilibrium is central to the folding polymorphism of α-syn, we decided to characterize the propensities of these two secondary structures in WT, A30P, A53T, and E46K monomers using MD.

Because the huge conformational space of α-syn is out of range for all-atom MD simulations in the explicit solvent, we applied coarse-grained UNRES MD (Maisuradze et al., 2010) to sample the structures of WT and the most studied A30P, A53T, and E46K monomers. The effective time scale of UNRES is 3 orders of magnitude larger than the all-atom time scale (Khalili et al., 2005). To ensure the convergence of the conformational sampling of the monomeric states, we applied replica exchange MD (see the Material and Methods section). The total effective time scale of the present simulations is 29.7 ms (72 replicas of 412 μs each) for each protein studied. One of the gold standards to quantify the secondary structure elements of a protein from its structure is the Dictionary of Secondary Structure of Proteins (DSSP) (Kabsch and Sander, 1983; Touw et al., 2015) algorithm based on a simplified model of hydrogen bonds. Application of DSSP to coarse-grained structures simulated by UNRES requires to build a compatible all-atom structure from the C^α coordinates of the UNRES model using reconstruction programs (Feig et al., 2004; Rotkiewicz and Skolnick, 2008). To avoid the high computational cost of all-atom reconstructions from coarse-grained coordinates, we developed here an algorithm which assigns an α-helix or a β-sheet secondary structure to each residue based on the C^α -C^α distances and on the coarse-grained angles formed by C^α -C^α pseudobonds, which correspond to the local curvature and torsion of the protein main chain (Grassein et al., 2020). The accuracy of the present algorithm, named CUTABI (CUrvature and Torsion based of Alpha helix and Beta-sheet Identification), to quantify the β-sheet content of proteins is improved compared to an existing algorithm based on C^α coordinates [P-SEA (Labesse et al., 1997)] and is comparable to the accuracy of DSSP (see the Material and Methods section). For each structure of the conformational ensemble of WT and mutants monomers, the number of residues in α-helix or in β-sheet was computed with CUTABI. The probability density of this two-dimensional descriptor was computed by using the conformations at 300 and 310 K. Analysis of these maps and of the conformational ensembles of the WT protein and mutants revealed subtle effects of the single amino-acid substitutions which are possibly related to the differences observed in oligomerization between WT and mutants (see the Results and Discussion section).

The present analysis is based on a first-principle (no bias or constraints applied) very large conformational sampling of the WT α-syn and mutants. As in any MD simulations, the sampling is never complete and each force-field has some bias. Next, we compare the predictions of UNRES for WT and mutants to available experimental data and previous theoretical studies.

To illustrate the difficulties of producing a conformational ensemble of α-syn, we compared (see Figure 10) the helix and β-sheet probabilities computed in four different works using MD simulations for WT α-syn (Ullman et al., 2011; Yu et al., 2015; Ramis et al., 2019; Coskuner and Wise-Scira, 2013). Compared to the UNRES simulations and Ref. Ullman et al., 2011, the calculations of Ref. Yu et al., 2015 largely overestimate the helical properties of α-syn, whereas the simulations of Ref. Ramis et al., 2019 largely underestimate them. The simulations of Ref. Coskuner and Wise-Scira, 2013 also predict a larger helical propensity than the one found in the present work (Figure 8) and in Ref. Ullman et al., 2011 (Figure 10). It is worth noting that the helical region nearby residue 60 is found in the UNRES simulations (Figure 8) and in Refs. Ullman et al., 2011; Coskuner and Wise-Scira, 2013 (Figure 10).

Comparison with our predictions shown in Figure 8 indicates better agreement with Ref. Ullman et al., 2011 for the helical part, except that UNRES predicts a larger probability to form an α-helix at residue A78. In the present work, MD simulations are performed for the entire protein. In contrast, the method of Ref. Ullman et al., 2011 consists of a combination of MD simulations of peptides of α-syn with sequences which overlap each other and are extracted from the α-syn sequence. The conformation of the protein is rebuilt by a combination of these fragments with weights using NMR data. The formation of transient helical segments was observed experimentally in the WT α-syn tetramer in Ref. Wang et al., 2011 between residues 4—43 and 50—103. Larger probabilities of helical tendency were found experimentally in regions 20—23 and 48—90 (Wang et al., 2011). These findings are rather in good agreement with the transient existence of helices in the α-syn monomer, as shown in Figure 8 (top panel), including the peak at A78. The stability of the tetrameric form is supported by previous all-atom MD simulations showing that the stability of this helical oligomer is due to a hydrophobic core formed by non-polar residues pertaining to the second α-helix of each chain, along with salt bridges formed by positively charged lysine residues and negatively charged aspartate and glutamate residues (Cote et al., 2018).

A detailed quantitative comparison between the UNRES β-sheet propensities (Figure 8) and those predicted in Refs. Ullman et al., 2011; Coskuner and Wise-Scira, 2013; Yu et al., 2015; Ramis et al., 2019 is difficult. However, we observe that there is a significant tendency to form β-sheet conformations over the entire sequence in Refs. Ullman et al., 2011; Ramis et al., 2019. On the contrary, simulations of Refs. Coskuner and Wise-Scira, 2013; Yu et al., 2015 predicted much localized and lower β-sheet probabilities.

Although the probabilities represented in Figure 8 for intra-molecular β-sheets of an isolated monomer cannot be simply compared to the formation of intra-molecular contacts in fibrils, the peaks in UNRES predictions for the β-sheet propensities (Figure 8) agree with the pairs of residues forming intra-molecular contacts in protofilaments, i.e., 47–79, 48—78, 92—71, 93—70, 94—69, 95—68 (PDB ID: 2n0a), and 47—79, 48—78, 92—69, 93—68 (PDB ID 6h6b). In addition, the maxima in Figure 8 for β-sheet propensity are located at or close to valine residues [P _V3 = 0.52, P _V26 = 0.46, P _V40 = 0.57, P _V49 = 0.4, P _V55 = 0.69, P _V63 = 0.64, P _V70 = 0.32, P _V82 = 0.31, P _V95 = 0.35], which is the most frequently found amino acid in β-sheets (Chou and Fasman, 1974). The peak at A90 is in the hydrophobic stretch 88IAAA91. Present simulations show a significant propensity to form a β-sheet in the region 71—83 necessary for aggregation (Giasson et al., 2001). Except for G73, the probability is indeed about 40% in this region (Figure 8, bottom panel). However, as mentioned above, in the present simulations, the probability of helical tendency at A78 is also about 40%. The results shown in Figure 8 are thus compatible with the probability to form both the β-sheet and helix in the region 71—83. The helix can be stabilized in the tetrameric oligomer (Bartels et al., 2011; Wang et al., 2011; Cote et al., 2018).

Although the propensities along the sequence are difficult to compare to experimental data, more global metrics can be used. Indeed, an interesting experimental parameter is the average content in helix and β-sheet conformations of WT α-syn in solution. Circular dichroism (CD) data reported an average of 2 ± 3% and 11 ± 7% for helix and β-sheet contents, respectively (Rekas et al., 2010). It is worth noting that the algorithms to extract β-sheet conformations from CD spectra are not as accurate as those for helices (Micsonai et al., 2015) and, as shown for α-syn, the CD spectrum is also dependent on the buffer and concentration (Araki et al., 2016). Other experimental values extracted from CD for the SSE content of WT monomers were reported: < 2% for the helix and 30% for β-sheet contents in Ref. Weinreb et al., 1996 and 3 ± 1% for the helix and 23 ± 8% for β-sheet fractions in Ref. Davidson et al., 1998. In their construction of the WT conformational ensemble based on MD simulations of α-syn fragments using the CHARMM force field constrained by NMR data, Ullman et al. reported the values 2% for the helix and 11% for β-sheet fractions (Ullman et al., 2011) (orange curve in Figure 10). In the present work, the global proportions of residues in the helix and β-sheet computed from the complete converged conformational WT ensemble are 10 and 31%, respectively. There are about 3 times more residues in the β-sheet than in the helix, a ratio which is in agreement with the one estimated from CD in Ref. Rekas et al., 2010. Interestingly, the theoretical values reported by coarse-grained simulations of WT α-syn with a different force field (red curve in Figure 10) are 20 ± 4% and 26.8 ± 6.8% for helix and β-sheet fractions, respectively (Ramis et al., 2019). In the present simulations, one notes that half of the helix fraction arises from the short helices located in the C-terminal region. The helix fraction without the C-terminal region is only 5% here. The fractions of the helix and β-sheet are on the same order of magnitude for the mutants and WT: A30P (helix = 10%, β-sheet = 34%), A53T (helix = 13%, β-sheet = 26%), and E46K (helix = 10%, β-sheet = 29%). Again, A53T shows a larger deviation compared to WT in agreement with a less populated B state (Figures 4 and 5). Most likely, the present force field overestimates the formation of SSE but predicts the correct equilibrium between the two main SSEs.

Another global structural parameter is R _g . The average radius of gyration R _g measured by small-angle X-ray scattering for WT α-syn in solution depends on the protein concentration, the presence of dimers or trimers, buffer type, pH, acetylation, and the source of proteins (Araki et al., 2016). For recombinant α-syn in 10 mM ammonium acetate pH 7.4, R _g = 27.2 ± 0.44 Å extrapolated at infinite dilution, which is comparable to the average computed value (R _g = 24.7 Å, Figure 7). Addition of HCl and/or NaCl increases R _g significantly to 33—40 Å (Li et al., 2002; Allison et al., 2009; Araki et al., 2016), but these effects cannot be tested with the present model. It is interesting to note that R _g for dimers in the present simulations are about 10 Å larger than that of monomers, as will be discussed elsewhere. Therefore, any mixture of monomers and dimers increases the effective radius of gyration of the solution.

The effects of the missense mutations on the structural properties can be summarized as follows: for the SSE tendency (Figure 8), the most significant effects of the single mutations are for A30P and A53T. Residue P30 has a huge effect on neighboring residues for the formation of the β-sheet. Residue T53 has a great influence on the helical region centered at K58. For the B state, structures with low R _g (around 18 Å) are more compact for A30P and E46K due to the large average number of contacts between the residues 1—20 and 96—140. For the entire ensemble of conformations, the average number of contacts between these two regions is rather similar for WT (1.84), A30P (2.12), and E46K (1.61) but significantly less for A53T (0.75). The mutant A53T is thus expected to be more flexible and less compact as shown by its larger R _g . The analysis is slightly different if the average number of contacts is computed between the entire N-terminal (1—60) and C-terminal (96—140) regions for which the average number of contacts is 4.54 (WT), 4.22 (A30P), 2.71 (E46K), and 2.61 (A53T). One finds that both E46K and A53T are much more flexible than WT. NMR studies of mutants A30P and A53T showed a reduction of contacts between C- and N-terminals for both mutants compared to WT (Bertoncini et al., 2005). The present MD results fully agree with these data for A53T.

A major finding in the present work is that an isolated α-synuclein in solution occurs in two phases, which are clearly visible in Figures 4 and 5. Unfortunately, it is difficult to construct experimentally a two-dimensional map of SSE propensities such as in Figure 4. A possible road toward such experimental analysis could be the use of Raman single-molecule spectroscopy (Leray et al., 2016; Dai et al., 2021) as α-helices and β-sheets have been well described by Raman fingerprints (Dai et al., 2021). About 40% of conformations are in the B state for WT, A30P, and E46K and 25% for A53T, and such a significant fraction of the conformational ensemble might be detectable by spectroscopy.