Is the Conformational Ensemble of Alzheimer’s Aβ10-40 Peptide Force Field Dependent?

By applying REMD simulations we have performed comparative analysis of the conformational ensembles of amino-truncated Aβ10-40 peptide produced with five force fields, which combine four protein parameterizations (CHARMM36, CHARMM22*, CHARMM22/cmap, and OPLS-AA) and two water models (standard and modified TIP3P). Aβ10-40 conformations were analyzed by computing secondary structure, backbone fluctuations, tertiary interactions, and radius of gyration. We have also calculated Aβ10-40 3JHNHα-coupling and RDC constants and compared them with their experimental counterparts obtained for the full-length Aβ1-40 peptide. Our study led us to several conclusions. First, all force fields predict that Aβ adopts unfolded structure dominated by turn and random coil conformations. Second, specific TIP3P water model does not dramatically affect secondary or tertiary Aβ10-40 structure, albeit standard TIP3P model favors slightly more compact states. Third, although the secondary structures observed in CHARMM36 and CHARMM22/cmap simulations are qualitatively similar, their tertiary interactions show little consistency. Fourth, two force fields, OPLS-AA and CHARMM22* have unique features setting them apart from CHARMM36 or CHARMM22/cmap. OPLS-AA reveals moderate β-structure propensity coupled with extensive, but weak long-range tertiary interactions leading to Aβ collapsed conformations. CHARMM22* exhibits moderate helix propensity and generates multiple exceptionally stable long- and short-range interactions. Our investigation suggests that among all force fields CHARMM22* differs the most from CHARMM36. Fifth, the analysis of 3JHNHα-coupling and RDC constants based on CHARMM36 force field with standard TIP3P model led us to an unexpected finding that in silico Aβ10-40 and experimental Aβ1-40 constants are generally in better agreement than these quantities computed and measured for identical peptides, such as Aβ1-40 or Aβ1-42. This observation suggests that the differences in the conformational ensembles of Aβ10-40 and Aβ1-40 are small and the former can be used as proxy of the full-length peptide. Based on this argument, we concluded that CHARMM36 force field with standard TIP3P model produces the most accurate representation of Aβ10-40 conformational ensemble.


Introduction
Aβ peptides linked to the development of Alzheimer's disease (AD) are produced, through a normal cellular proteolysis, in a variety of alloforms, which differ with respect to sequence length and the extent of amino-or C-terminal truncation [1][2][3]. The most abundant form is a 40-residue version Aβ1-40, which constitutes about 90% of all Aβ species in cerebrospinal fluid [4]. Virtually all Aβ peptides are highly amyloidogenic [5,6] and play a central role in amyloid cascade hypothesis, which explains AD pathogenesis on the basis of multi-stage aggregation of Aβ species. In this process, Aβ monomers represent initial species involved in spontaneous aggregation. Moreover, according to experimental studies fibril elongation is also largely driven by deposition of Aβ monomers to the fibril edges [7]. Generally, Aβ peptides display a high level of cytotoxicity [2,[8][9][10], which is related to their ability to readily bind to cellular lipid bilayers and disrupt their structure [11]. Although the mechanism of binding to lipid bilayers is likely to be concentration dependent, Aβ peptides predominantly bind as monomers rather than oligomers at nanomolar concentrations [12,13].
Aβ peptides belong to the class of intrinsically disordered proteins implying that they lack well defined native structure in aqueous environment. Indeed, experimental investigations, including solution NMR studies, have shown that the conformational ensemble of Aβ monomer in water is populated by heterogeneous coil-like conformations [14][15][16][17]. More recently, several NMR measurements, including chemical shifts, nuclear Overhauser effects, and J-couplings, have been used to confirm that Aβ peptides adopt generally random coil structures at neutral pH [18]. At the same time, careful analysis of electron paramagnetic resonance studies revealed that Aβ monomers still contain short structured regions (His14-Val18, Gly29-Ala30, and Gly38-Val40) under normal physiological conditions [19]. Similarly, several NMR studies have pointed to a formation of a turn or bend structures in the sequence region (Phe20-Ser26) between the central hydrophobic cluster (Leu17-Ala21) and the C-terminal (Ala30-Val40) [15,17]. Due to generally disordered state of Aβ in water, it is not surprising that Aβ conformations are highly dependent on solvent properties. For example, in the membrane-like environments Aβ conformational ensemble undergoes considerable reorganization manifested in the formation of helical structure in the sequence regions Glu15-Val24 and Gly29-Met35 [20][21][22]. Moreover, Aβ helix propensity is pH dependent illustrated by the observation that the Glu15-Val24 helix, but not the C-terminal helix, becomes destabilized at normal pH. Conformational plasticity of Aβ peptides is also consistent with mutagenesis studies. For example, Iowa (D23N) or Osaka deletion (E22Δ) mutants aggregate significantly faster than the wild-type [23,24]. Similarly, according to in vitro studies many single-point mutations mainly affecting hydrophobic or charged residues can either accelerate or reduce Aβ aggregation propensity [25,26]. One may expect that Aβ monomeric conformations, being the initial species involved in aggregation, are impacted by these single-point mutations.
A survey of experimental findings presented above suggests that computational characterization of the conformational ensemble formed by Aβ monomers is important for understanding its aggregation and cytotoxicity. Several previous molecular dynamics studies have probed Aβ monomers in water. Garcia and coworkers have studied the conformations of Aβ1-40 and Aβ1-42 peptides using replica exchange molecular dynamics (REMD) simulations extended to microsecond timescales [16,17,27]. Consistent with the experiments, their analysis revealed generally disordered Aβ conformational ensemble augmented by several structured regions, especially in the C-terminal of Aβ1-42, where a β-hairpin has been detected. Qualitatively similar conclusions have been reached in the recent REMD study conducted by Head-Gordon and coworkers [28]. In our previous studies, we have utilized REMD and all-atom CHARMM22 force field with CMAP corrections to investigate the conformations of amino-truncated Aβ10-40 peptide in water [29]. We found that, similar to the full-length peptide, Aβ10-40 samples predominantly turn and random coil structures, whereas helical and especially β-sheet propensities are low. Furthermore, the peptide almost completely lacks tertiary structure with most stable intrapeptide interactions forming between the amino acids adjacent along the sequence. In light of disordered state of Aβ monomer it is important to test the dependence of its conformational ensemble on the force fields employed in the simulations. Recently, such investigation has been carried out for Aβ1-40 monomer, which was probed using OPLS-AA/L, AMBERff99sb-ILDN, and CHARMM22 Ã protein force fields and several water models [30]. All the three force fields predicted the formation of β-structure in the Leu17-Ala21 and Ala30-Leu34 regions, but with markedly different β propensities. In particular, OPLS-AA/L and AMBERff99sb-ILDN simulations revealed stable β-structures, whereas suppressed β fraction has been observed in the CHARMM22 Ã force field. It is conceivable that the three force fields overestimate the β propensity in Aβ monomer, which is expected to be low according to the NMR studies [18]. Indeed, REMD study of two natively unfolded peptides, NTL9 (1-22) and NTL9 (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17), using AMBERff99sb-ILDN, CHARMM22/CMAP, and CHARMM36 force fields [31] has shown that both CHARMM force fields predict much smaller β fraction than AMBERff99sb-ILDN. Additionally, the parameterization of water may affect peptide conformational ensembles. This point has been recently demonstrated using REMD and CHARMM36 for two Ala-rich peptides and GB1 peptide, which form considerably more solvated and extended structures with modified TIP3P water model compared to its standard version [32].
Because previous investigations have emphasized the importance of force field parameterization, a natural question arises about how a force field affects sampling of amino-truncated Aβ10-40 peptide, which was extensively studied by us in the context of peptide-lipid bilayer interactions [33][34][35][36]. To address this question, we used all-atom REMD simulations and performed a systematic comparison of Aβ10-40 conformational ensembles in four protein force fields (CHARMM36, CHARMM22 Ã , CHARMM22/CMAP, OPLS-AA) and two water models (standard and modified TIP3P). We show that although all force fields are consistent in predicting the Aβ propensity to form turn and random coil structures, they strongly disagree on the extent and distribution of tertiary interactions or helix and β propensities in Aβ monomer.
We have also compared the J-coupling and residual dipolar coupling (RDC) constants computed from Aβ10-40 simulations with Aβ1-40 experimental data. Surprisingly, we found that in silico Aβ10-40 conformational ensemble produced by CHARMM36 force field agrees better with the experimental measurements than in silico Aβ1-40 ensemble simulated earlier [17,27]. Based on these comparisons we suggest that Aβ10-40 can serve as a proxy of the full-length Aβ1-40 peptide and the CHARMM36 force field with standard TIP3P water model provides most accurate reproduction of Aβ conformational ensemble.

Molecular model and simulations
To explore the impact of force field parameterizations, we have performed molecular dynamics (MD) simulations of Aβ10-40 peptide, which is an amino-truncated fragment of the fulllength peptide Aβ1-40 (Fig 1a). In all, we have investigated four all-atom protein force fields and two explicit water models [37,38] resulting in five simulation systems, which utilized CHARMM36 [39] with modified TIP3P water model (denoted as C36), CHARMM36 with standard TIP3P water model (C36s), CHARMM22 Ã with modified TIP3P water model (C22 Ã ) [40], CHARMM22 with CMAP corrections and modified TIP3P water model (C22cmap) [41], and OPLS-AA with modified TIP3P water model (OPLS-AA) [42]. The C22cmap system was already studied by us previously [29,43,44] and is used here to expand the force fields comparison. All simulation systems contained a single Aβ10-40 peptide, 4959 water molecules, and one sodium ion to set the net system charge to zero. In all, the simulation systems contained 15,354 atoms. Aβ termini were capped with acetylated and aminated groups. The charged states of amino acids corresponded to neutral pH (in particular, histidines were deprotonated). For all force fields we used periodic boundary conditions with the cubic unit cell having the edge dimension of about 53.8 Å and resulting in the mass density of 0.9848 g/cm 3 at 330 K. Non-bonded interactions were computed using a smooth switching functions acting within the interval of 8 and 12 Å. Electrostatic interactions were computed using particle mesh Ewald summation with the grid size of % 1Å. All hydrogen associated covalent bonds except in water molecules were treated as rigid by applying ShakeH algorithm. Water molecules were treated as rigid using SETTLE algorithm. All simulations used the integration step of 1 fs. Full electrostatic evaluation frequency was set to 4 integration steps.

Replica exchange protocol
To produce exhaustive sampling of Aβ10-40 conformational ensembles we have utilized canonical (NVT) replica exchange molecular dynamics (REMD) [45], which is implemented in NAMD MD program [46]. For all systems we used R = 40 replicas distributed exponentially in the temperature range from 300 to 440 K. Canonical ensembles in the replicas were generated by applying underdamped Langevin dynamics with the damping coefficient γ = 5 ps −1 . Replica exchanges were attempted every 2 ps between all neighboring replicas along a temperature scale resulting in the average acceptance rates ranging between 27 and 29%. For each simulation system we have produced four REMD trajectories collecting in total 3.2 μs of sampling (or 80 ns per replica). Each REMD trajectory has been initiated with random initial conformations equilibrated at 330 K with preliminary isothermal-isobaric simulations to set correct mass densities. In addition, each replica was equilibrated at its own temperature for an additional 1 ns. With this approach no sampling data from REMD trajectories needed to be discarded as non-equilibrated. Although probing REMD convergence is generally a difficult task (see, e.g., [47]), our analysis in S1 Text suggests that the resulting simulation times appear sufficient for sampling convergence.

Computation of structural probes
Secondary structures in Aβ were assigned using the STRIDE program [48]. A helical state includes α-, 3 10 -, or π-helix conformations, whereas a β-strand state includes extended conformations or isolated bridges. Tertiary interactions were probed by side chain contacts. A contact occurs if the distance between the geometric centers of heavy atoms in two side chains is  less than 6.5 Å. This cutoff approximately corresponds to the onset of hydration of side chains as their separation increases. We classified the contacts between residues i and j as long-range if |j − i| ! 5 or short-range otherwise. To explore the fluctuations in peptide backbone, we computed the standard deviations δϕ(i) and δψ(i) for backbone dihedral angles ϕ and ψ of a residue i. These standard deviations are referred to as root-mean-square fluctuations (RMSF). To evaluate peptide dimensions we computed the radius of gyration, R g , using the positions of side chain centers of mass and C α atoms.
To quantitatively evaluate the consistency between the experimental and computational conformational ensembles we have utilized two quantities. The first is the 3 J HNHα coupling constants associated with three-bond coupling interaction between H N and H α protons. These constants are sensitive to peptide's secondary structure [49]. In this study we used three sets of experimentally measured J-coupling constants, J exp , for Aβ1-40 peptide [16][17][18]. In silico J-coupling constants, J comp , were determined from the backbone dihedral angles using Karplus equation [50] where ϕ is a backbone dihedral angle and A, B, and C are the coefficients determined by fitting with the experimental data. We  [53]. The N-terminal amino acid was excluded from our computation of J-couplings as its ϕ angle may be distorted by the capping group.
As a second quantity we chose residual dipolar coupling (RDC) constants experimentally measured for Aβ1-40 peptide [54]. RDC measurements probe orientation of amide NH bonds in protein backbones with respect to external magnetic field, and can identify long-range structural correlations across biomolecular structure. To compute RDC constants from in silico structures we used the Prediction of ALignmEnt from Structure (PALES) program [55]. We processed our simulation structures by using default PALES settings, including the application of steric interaction model, which determines alignment orientation based on steric properties of a molecule. As described in S1 Text we applied global alignment of Aβ10-40 structures for RDC computations. Once the RDC constants were produced, they were multiplied by the scaling factors determined from the least-squares fitting that minimizes the deviation between the experimental and in silico data.
To assess the similarity between experimental and computational quantities, we used the Pearson's correlation coefficient (PCC), root mean square deviations (RMSD), and quality fac- , where D exp,k and D comp,k are the experimental measurements and their computed counterparts available for amino acids k. The sum is taken over all amino acids for which experimental and computed values are simultaneously available. We applied Q to evaluate agreement for J-coupling or RDC constants. To preserve consistency with our previous studies all ensemble averages were computed using weighted histogram analysis method (WHAM) [56] of REMD data at 330K. J-coupling and RDC constants were computed at 300K, which is the closest simulation temperature to experimental conditions (% 280K) [16][17][18].

Results
Conformational ensemble of Aβ10-40 monomer in CHARMM36 force field with modified TIP3P water model  Table 1 that its conformational ensemble is dominated by random coil (hRCi = 0.49 ± 0.04) and turn (hTi = 0.44 ± 0.03), which together account for 93% of all amino acid states. Indeed, the turn structure is stable (i.e., its fraction hT(i)i > 0.5) for few S1 residues (His13, His14) and within the long sequence span Phe19-Gly29 covering the part of hydrophobic region S2, the entire hydrophilic S3, and the beginning of C-terminal S4. Random coil structure dominates the N-and C-terminals and the Gln15-Val18 region. The occurrence of     fairly uniform throughout Aβ sequence. The most rigid backbone conformation is observed for hydrophobic Phe20. The values of δϕ and δψ averaged over the entire Aβ10-40 sequence are 55.2 ± 1.6˚and 81.9 ± 3.8˚, respectively. Relatively uniform distribution of the backbone fluctuations in Aβ sequence is consistent with the lack of well-defined secondary structure, such as helices or β-strands.
We next investigate the Aβ10-40 tertiary structure by probing the formation of intrapeptide interactions. Fig 5a presents the peptide contact map hC(i, j)i, which visualizes the probabilities of forming contacts between amino acid side chains. This figure suggests that there are very few stable interactions in Aβ peptide, i.e., those occurring with the probability hC(i, j)i > 0.35. According to Table 2 there are no stable long-range (|i − j| ! 5) contacts, and there are only two stable short-range (|i − j| < 5) contacts, namely, Leu17-Phe19 and Asp23-Ser26. The Leu17-Phe19 contact is likely to explain a rigid backbone conformation at Phe20. The probability of forming a salt bridge Asp23-Lys28, which is important for amyloid fibril assembly, is low being equal to hC(23, 28)i = 0.10 ± 0.03. However, weak electrostatic interactions are formed between Glu22 and Lys28 (0.21 ± 0.13) and between Glu11 and Lys16 (0.20 ± 0.11). Overall, the average number of all side chain contacts forming in Aβ monomer is hCi = 21.1 ± 1.4, of which 10.7 ± 0.9 (or 51%) are long-range.
Finally, we consider the probability distribution P(R g ) of Aβ10-40 radius of gyration R g . Fig  6 shows that for C36 P(R g ) is broad with the maximum at R g % 15 Å. The equilibrium value of R g , hR g i, is 16.9 ± 0.5 Å, which is the largest among all force field considered (see below). In addition, we have computed the average end-to-end distance hR 1N i = 24.6 ± 1.2Å, which is also the largest among other force fields. Thus, Aβ10-40 peptide in C36 force field lacks stable secondary or tertiary structure and forms expanded conformations as illustrated in Fig 1b. Conformational ensemble of Aβ10-40 monomer in CHARMM36 force field with standard TIP3P water model To check the impact of water model, we repeated CHARMM36 REMD simulations of Aβ10-40 peptide using standard TIP3P water (denoted as C36s). Following the analysis for C36 we have computed Aβ secondary and tertiary structure.  Table 1 demonstrate that in close agreement to C36 simulations Aβ10-40 primarily samples random coil (hRCi = 0.47 ± 0.04) or turn (hTi = 0.45 ± 0.03) conformations, which together represent 92% of all residue states. In contrast, helix or β-state occur rarely. Fig 3 further reveals that the residue-specific turn hT(i)i and helix hH(i)i propensities are almost identical to those observed in C36. Similar to C36, stable turn structure is present at His13-His14 and in the region Phe19-Gly29. Consequently, the root-mean-square deviations (RMSD) between hT(i)i and helix hH(i)i distributions computed from C36s and C36 simulations are low (0.02 and 0.03). Random coil distributions between C36s and C36 also nearly match. Furthermore, according to Fig 4a and 4b a very close agreement is observed in the RMSF distributions, δϕ(i) and δψ(i), as confirmed by the low respective RMSD values (3.9 and 3.2˚). As expected the average values hδϕi and hδψi computed using all amino acids are 54.6 ± 1.9˚and 81.1 ± 5.0˚, which are nearly identical to those of C36.
To analyze Aβ tertiary structure we have computed the contact map hC(i, j)i, which is presented in Fig 5b. In general, C36s tertiary interactions are very similar to those seen for C36. It follows from Table 3 that Aβ peptide has no stable long-range contacts and the same two stable short-range contacts as in C36 are formed. There are four common contacts among top five long-range interactions in C36s and C36 simulations (Ala21-Ser26, Glu22-Lys28, Glu11-Lys16, Val24-Gly29), whereas three short-range contacts are shared between the two simulations (Leu17-Phe19, Asp23-Ser26, Gly25-Asn27). As in C36 the probability of forming Conformational Ensemble of Alzheimer's Aβ10-40 Peptide the salt bridge Asp23-Lys28 is low (0.14 ± 0.07). Overall, the RMSD computed between C36 and C36 contact maps is small (0.02) confirming their similarity. It is noteworthy, however, that in C36s simulations the peptide forms, on an average, more side chain contacts than in C36, as their number reaches hCi = 23.5 ± 1.4, of which 12.8 ± 0.8 (or 54%) are classified as long-range. Furthermore, comparison of the probability distributions P(R g ) for Aβ10-40 radius of gyration in Fig 6 reveals that the C36s distribution is systematically shifted to smaller R g values suggesting that the C36s peptide is more compact. Indeed, from Fig 6 we determine that the equilibrium hR g i is 15.9 ± 0.4Å, which is smaller than the C36 value. Thus, although the secondary and tertiary structure propensities in C36s and C36 simulations are in close agreement, the former generates more compact structures with larger number of intrapeptide interactions (Fig 1c).  Fig 4c display distinctive features characteristic of C22 Ã force field. Specifically, the plot reveals two sequence regions with suppressed backbone fluctuations, Val17-Asp23 and Leu34-Val36, which approximately coincide with the formation of stable turn and marginally stable helix structures. Furthermore, in C22 Ã simulations for most sequence positions δϕ(i) is much smaller than in C36. Overall, the average RMSF hδϕi and hδψi observed in C22 Ã are 30.2 ± 1.5˚and 66.1 ± 4.6˚, i.e., compared to C36 a particularly significant decrease by a factor of 1.8 is seen in hδϕi. It is then not surprising that the RMSD for δϕ(i) and δψ(i) distributions computed between C22 Ã and C36 simulations are 31.1˚and 33.0˚, respectively, which are about an order of magnitude larger than those comparing C36 force fields.
The C22 Ã equilibrium contact map hC(i, j)i displayed in Fig 5c shows a formation of multiple side chain contacts, of which eight long-range and ten short-range contacts are stable. Moreover, Table 4 lists top five long-range contacts, one of which, the salt bridge Lys16-Asp23, is formed with exceptionally high probability of 0.85 ± 0.02. Two hydrophobic long-range contacts, Phe19-Val24 and Val24-Ala30, also occur with high probabilities (0.52 ± 0.04 and 0.46 ± 0.08). Interestingly, in C22 Ã simulations Asp23-Lys28 salt bridge is effectively disrupted (0.06 ± 0.03). Also, multiple very stable short-range contacts are observed, such as the helixlike contact Gly33-Gly37 (0.86 ± 0.03). Overall, the equilibrium number of side chain contacts in Aβ10-40 is hCi = 29.0 ± 0.8, of which hC LR i = 15.4 ± 0.3 (or 53%) are long-range. Compared to C36, the values of hCi and hC LR i are larger by 37% and 44%, respectively. Finally, we note that none of the top five C22 Ã long-or short contacts are observed in C36 or C36s simulations. Consequently, the RMSD value measuring the difference between C36 and C22 Ã contact maps is much larger (0.12) than the RMSD comparing C36 and C36s force fields. As seen in Fig 6 the probability distribution P(R g ) for Aβ10-40 monomer computed from C22 Ã simulations peaks at smaller values of R g and is more narrow than for either of C36 force fields. From Fig 6 we find hR g i = 14.5 ± 0.2Å, which is smaller than the respective values for C36 simulations (16.9 or 15.9 Å). The same conclusion applies to the end-to-end distance (hR 1N i = 20.0 ± 0.9 Å). In summary, compared to C36 simulations, C22 Ã force field significantly enhances turn and, particularly, helical propensities, makes several sections of Aβ backbone rigid, and dramatically strengthens tertiary interactions as illustrated in Fig 1d resulting in peptide compaction. Therefore, there is little consistency between the conformational ensembles mapped using C22 Ã and C36 force fields.
Conformational ensemble of Aβ10-40 monomer in CHARMM22 force field with CMAP corrections In our previous work [29,43], we have used REMD to study the conformational ensemble of the Aβ10-40 monomer using CHARMM22 force field with CMAP corrections (denoted as C22cmap). Below we extend our previous C22cmap analysis to provide comparison with other force fields. We begin by focusing on Aβ secondary structure. Fig 2 and Table 1 demonstrate that in C22cmap Aβ mainly forms turn (hTi = 0.49 ± 0.01) and random coil (hRCi = 0.38 ± 0.01) structures. These observations are qualitatively consistent with other force fields. However, compared to both C36 the population of helical structure in C22cmap (hHi = 0.12 ± 0.01) is elevated, whereas the formation of β structure is still rare. Fig 3 demonstrates residue-specific secondary structure propensities for C22cmap. Stable turn structure is observed in His13, Phe19-Gly25, Asn27-Gly29, and Met35-Gly37, which with the exception of the last region match well the C36 turn distribution. Random coil, hRC(i)i, is observed in Tyr10-Glu11, Gln15-Val18, and Val39-Val40, i.e., it follows the respective C36 propensities. Fig 3a also reveals an appearance of marginally stable helix in the C-terminal region S4 (Ile32-Val36), which approximately coincides with the distribution of C22 Ã helix. The RMSD values comparing C22cmap and C36 distributions of turn and helix structure are 0.11 and 0.12, respectively, implicating moderate differences between the two systems. Fig 4d demonstrates that the fluctuations of backbone dihedral angles, δϕ(i) and δψ(i), with the exception for Gly residues, agree generally well with those computed for C36. However, in contrast to C36 distributions, the backbone fluctuations are suppressed in a wider interval of Val18-Val24 compared to a single position Phe20 in C36. The averages hδϕi and hδψi are 49.5 ± 1.0˚and 74.0 ± 1.2˚, which are somewhat smaller than the respective averages computed from C36 simulations. The RMSD values comparing δϕ(i) and δψ(i) between C22cmap and C36 are 13.7˚and 21.1˚, respectively. These RMSD values are much larger than those comparing C36 and C36s force fields, but significantly smaller (by 2.3 or 1.6 times, respectively) than the RMSD probing the difference between C36 and C22 Ã . In line with these observations, the C22cmap fluctuations in Fig 4d, particularly δϕ(i), are generally higher than those observed for C22 Ã . Fig 5d presents the equilibrium contact map hC(i, j)i computed using C22cmap simulations. In all, we have detected the formation of nine short-range and one long-range (Lys16-Asp23) stable contacts that stands in contrast to C36 results, which showed only two stable short-range contacts. The Asp23-Lys28 salt-bridge, which is disrupted in C22 Ã or C36, has also low probability of occurrence (0.15 ± 0.03). Moreover, Table 5 demonstrates none of the top long-or short-range contacts formed in C22cmap force field are observed in any of C36 simulations. The overall number of side chain contacts formed in C22cmap is hCi = 26.2 ± 0.2, of which hC LR i = 12.3 ± 0.2 (or 47%) are long-range. Compared to C36 simulations, hCi and hC LR i are larger by 24 and 15%, respectively. As a result, the RMSD comparing the C22 and C36 contact maps is moderately large (0.07), exceeding the RMSD between C36 and C36s simulations (0.02), but still being smaller than the RMSD comparing C36 and C22 Ã (0.12).
Last, we examined the probability distribution P(R g ) plotted in Fig 6. The distribution observed for C22cmap shows surprisingly good agreement with C36 (particularly, C36s) force fields. From P(R g ) we found hR g i = 15.8 ± 0.2 Å, which is almost equal to that observed for C36s (hR g i = 15.9 Å), somewhat smaller than that of C36 (16.9Å), but larger than the C22 Ã result (14.5Å). Taken together, C22cmap simulations exhibit a weak helix propensity in the Cterminal consistent with the formation of large number of stable short-range contacts that positions the C22cmap conformational ensemble in-between C36 and C22 Ã ensembles (see Discussion for further analysis). The representative structure of Aβ10-40 peptide in C22cmap force field is displayed in Fig 1e.

Conformational ensemble of Aβ10-40 monomer in OPLS-AA force field
To evaluate Aβ10-40 conformational ensemble in the force field unrelated to any CHARMM version, we have performed OPLS-AA REMD simulations. Similar to all previously considered force fields Fig 2 and Table 1 show that Aβ10-40 in OPLS-AA simulations adopt mainly turn (hTi = 0.46 ± 0.02) or random coil (hRCi = 0.40 ± 0.03) conformations, which together represent 86% of all amino acid states. However, in contrast to other force fields (e.g., C36) β-state fraction is elevated four-fold to hSi = 0.12 ± 0.02, whereas α-helix occurrence is negligible. Fig  3, which displays residue-specific secondary structure propensities, further underscores the characteristic OPLS-AA feature-an elevated sampling of β-structure in the S2 and S4 hydrophobic regions, where for some positions i hS(i)i reaches %0.4. The helix propensity across Aβ sequence is uniformly weak, whereas the turn and random coil fractions generally follow the other force field trends. Accordingly, the RMSD comparing the β distributions hS(i)i in OPL-S-AA and C36 simulations has a large value of 0.13, whereas the RMSD for helix and turn propensities are much smaller (0.04 and 0.08). Fig 4e displays  clearly exhibits more extensive tertiary interactions forming in OPLS-AA comparing to C36 force field. Indeed, the numbers of all and long-range contacts are increased by about 50% and 100%, respectively, to hCi = 30.8 ± 0.5 and hC LR i = 21.1 ± 0.6 resulting in the largest fraction of long-range interactions of 69% among all force fields tested by us. Interestingly, even though OPLS-AA promotes tertiary interactions, very few of them qualify as stable. Specifically, in sharp contrast to C22 Ã Table 6 lists only one stable long-range (Val18-Leu34) and one stable short-range (Tyr10-Val12) contacts. Among top five long-range contacts four are hydrophobic and three link the sequence regions S2 and S4 (Val18-Leu34, Val18-Val36, Leu17-Met35) with the probabilities of occurrence close to the stability threshold. Also, unique to OPLS-AA simulations, the salt-bridge Asp23-Lys28 appears among the top five long-range contacts. Finally, no long-or short-range top five contacts are shared between OPLS-AA and C36 simulations. It is then not surprising that extensive but weak tertiary interactions formed in OPLS-AA simulations lead to a fairly large RMSD value measuring the difference between the contact maps obtained in OPLS-AA and C36 force fields (0.07).
Probability distribution P(R g ) for Aβ10-40 radius of gyration presented in Fig 6 is far more narrow and reaches maximum at the smallest R g compared to any other force field. We determined from Fig 6 that the equilibrium value hR g i is 13.5 ± 0.2Å, which is the smallest among the force fields tested. Similar observation holds for the end-to-end distance (hR 1N i = 15.6 ± 1.5 Å). Thus, OPLS-AA force field promotes extensive but flickering tertiary interactions resulting in Aβ10-40 collapse and moderate enhancement of β-structure as shown in Fig 1f.

Comparison of experimental and computational J-coupling and RDC constants
We have compared 3 J HNHα -coupling and RDC constants computed from our simulations and measured experimentally (see Materials and Methods and Fig 7a). (It is important to note that experimental measurements refer to the full-length peptide Aβ1-40, whereas our simulations have examined the amino-truncated fragment Aβ10-40. This point is elaborated in Discussion.) We first investigated the agreement between 3 J HNHα -coupling constants. As stated in the Materials and Methods we used three sets of experimental 3 J HNHα -coupling constants, J exp ,  [18]. This combination of Karplus equation coefficients and experimental data provides the best agreement between J comp (i) and J exp (i). Only amino measured by Garcia and coworkers [16,17] and, more recently, by Bax and coworkers [18]. Furthermore, to compute 3 J HNHα -coupling constants in silico, J comp , we applied Eq (1) with three different sets of Karplus equation coefficients [51][52][53]. Because a priori it is unclear which combination of experimental data and Karplus equation coefficients is more accurate, we considered all nine possible combinations and for each computed J comp (i)-coupling constants using Aβ10-40 sequence positions i with available experimental measurements. Consistency between J comp (i) and J exp (i) was evaluated by calculating the root-mean-square deviation (RMSD), quality function Q, and Pearson's correlation coefficient (PCC) as described in the Materials and Methods. These quantities were then averaged over nine datasets and are presented in Table 7. Judged by RMSD, PCC, and Q values as well as their errors the best agreement between experimental and computational J-coupling constants is observed for C36 and C36s force fields. The ranking of other force fields in the descending order of agreement with the experiment is OPLS-AA, C22cmap, and C22 Ã . Notably, C22 Ã shows very significant increase in RMSD (by 41%) compared to C36 and has effectively no correlation with J exp (i) as measured by PCC (0.14).
To provide additional probe of Aβ10-40 conformational ensemble we have computed the RDC constants, RDC comp (i), as described in the Materials and Methods and S1 Text. To compare them with their experimental counterparts, RDC exp (i), [54] shown in Fig 7b, we used the same metrics as for 3 J HNHα -coupling constants. The results presented in Table 7 indicate that the best agreement between experimental and computational RDC constants is observed for C36s force field, which has the smallest RMSD or Q and the largest PCC of 0.65. Compared to C36s force field, C36 exhibits larger RMSD and Q (by % 20% for both) and significantly lower PCC (a 25% decrease). With respect to RMSD and Q the worst agreement is seen for C22cmap and C22 Ã , for which the RMSD or Q values are about two-fold larger than for C36s. Their PCC values are also lower than for C36s, but higher, especially of C22cmap, than for C36. Measured by RMSD and Q the agreement with the experiment places OPLS-AA between C36 and C22cmap/C22 Ã . However, OPLS-AA demonstrates the worst correlation between computational and experimental RDC constants as measured by PCC. Thus, taken together our data suggest that the force field providing the best consistency with the experimental data is C36s. If we disregard a negligible (within the error) difference in J-coupling RMSD for C36s and C36 (Table 7), this conclusion is supported by all comparison metrics. With the exception of RDC PCC, other metrics identify C36 as the next "best" force field. Again, with the exception of acids with experimentally available data are considered. (b) Distributions of RDC constants computed for Aβ10-40 peptide using PALES program [55] and REMD sampling generated with five force fields. Superimposed are experimental RDC constants (black lines with circles) measured by Wang and coworkers [54].
doi:10.1371/journal.pcbi.1005314.g007 PCC for RDC all the metrics determine that C22 Ã produces the worst agreement with the experimental data. The implications of our findings for the selection of force field and water model and for the differences between Aβ1-40 and Aβ10-40 conformational ensembles are presented in the Discussion.

Discussion
Comparison of Aβ10-40 conformational ensembles in different force fields Using REMD we have performed a comparative analysis of Aβ10-40 conformational ensembles generated by employing five force fields, which combine four protein parameterizations (C36, C22 Ã , C22cmap, and OPLS-AA) and two water models (standard and modified TIP3P). As a reference in our analysis we took the recent modification of CHARMM force field, CHARMM36, coupled with modified TIP3P water model. Its selection was motivated by recent tests showing that this force field provides the best agreement with experimental NMR data collected for six proteins [57]. Taken together, our results suggest several observations. First, all force fields produce fairly consistent distributions of secondary structure. According to Fig 2 and Table 1 all of them predict largely similar fractions of turn conformations (varying between 0.44 and 0.52) and, to a lesser extent, random coil (varying between 0.30 and 0.49). In all force fields, the turn structure dominates within approximately the same sequence regions, His13 and His14 in S1 region and a sequence interval Phe19-Gly29, whereas the random coil occurs at Aβ10-40 termini. Nevertheless, there are also significant variations among the force fields. For example, a unique feature of C22 Ã force field is a considerable helix bias in the S3 and S4 regions resulting in hH(i)i % 0.4 for few positions. Similarly, OPLS-AA differs from other simulations by significant β-structure propensity in the S2 and S4 regions, where hS(i)i peaks at % 0.4. Second, additional insight into peptide conformational ensemble is provided by the fluctuations in backbone dihedral angles, δϕ(i) and δψ(i). Similar to secondary structure, backbone fluctuations in C36 and C36s are in excellent agreement as evidenced by their consistent average values and small RMSDs. C22cmap differs moderately from C36 by having slightly more rigid backbone (i.e., smaller average hδϕi and hδψi values), whereas OPLS-AA, in contrast, demonstrates enhanced backbone fluctuations. However, the force field clearly standing apart from others is C22 Ã , which predicts suppressed fluctuations in two sequence regions (Val17-Asp23 and Leu34-Val36). This feature plus generally small fluctuations in ϕ angles result in the lowest averages hδϕi and hδψi as well as elevated, by an order of magnitude, RMSD value between C22 Ã and C36 simulations as opposed to that between C36 and C36s. Therefore, when secondary structure and backbone fluctuations are considered together, the force fields can be ranked in the descending order of similarity to C36 as C36s, which is nearly identical to C36, C22cmap and OPLS-AA, which moderately differ from C36, and C22 Ã , which reveals considerably differences due to helix formation and rigid backbone.
It is important to consider our results in the context of other studies of force field propensities. C22cmap tendency to bias peptide conformational ensembles toward helical states has been documented for Ala pentamer [58], to correct which C36 force field has been developed [39]. Our results show that C22cmap indeed produces a slightly elevated helix fraction in unstructured Aβ10-40 compared to C36 or C36s increasing it to 0.12 from 0.04 or 0.06, respectively. However, this bias is weak compared to that of C22 Ã , which demonstrates much stronger helix propensity in the C-terminal. Small differences in helix propensities between C22cmap and C36 have also been noted for the unstructured fragments of NTL9 peptides [31].
Our third observation is related to the distribution of tertiary interactions. In Aβ peptide C36 force field produces no stable long-range interactions and only two stable short-range contacts and, overall, it leads to the smallest numbers of tertiary contacts (21.1) and long-range contacts (10.7). Aβ conformations in this force field are also the least compact being characterized by the largest radius of gyration (16.9 Å). Thus, we conclude that equilibrium C36 conformations are dominated by expanded structures lacking stable interactions. C36s force field, which differs from C36 solely by the water model, generates very similar conformational ensemble, which also lacks stable interactions and shares four or three common top long-and short-range contacts with C36. Overall, the contact map RMSD for C36 and C36s is very low (0.02). Nevertheless, the C36s numbers of contacts, including long-range, are slightly larger (by 10-20%) than in C36. Additionally, compared to C36 C36s features slighly smaller Aβ radius of gyration. These results are consistent with recent comparison of standard and modified TIP3P water models, which showed that the latter enhances hydration and generates more open peptide conformations [32].
Although Aβ structures produced with C22cmap and C22 Ã share some similarity (three long-and short-range top contacts are common), C22 Ã is by far unique in generating the largest number of stable long-range contacts (eight as opposed to one in C22cmap). Some of these contacts, such as salt-bridge Lys16-Asp23, are effectively always formed. A distinctive characteristic of C22cmap is the large number (9) of stable short-range contacts combined with very few (1) stable long-range interactions. According to contact map RMSD, C22cmap differs moderately from C36 (0.07), whereas C22 Ã deviates from C36 by far larger degree (0.12). There are no common top long-or short-range contacts between C36 and C22cmap or C22 Ã . The tertiary interactions also set OPLS-AA force field apart from all other simulations. OPL-S-AA generates the largest number of all contacts and, more importantly, the largest number of long-range contacts, which is increased almost two-fold compared to other simulations (except for C22 Ã , with respect to which h C LR i increases one-third). As a result OPLS-AA conformations exhibit extremely large fraction of long-range interactions reaching almost 70% and, accordingly, Aβ adopts most compact structures (hR g i 13.5 Å). Similar observation concerning Aβ1-40 peptide collapse has been made previously for OPLS-AA/L force field [30]. Interestingly, OPLS-AA tertiary contacts, although numerous, are weak suggesting that Aβ samples compact but still disordered states. OPLS-AA and C36 do not share common top interactions, whereas the contact map RMSD between the two is moderate (0.07). Thus, using contact map RMSD and C36 as reference we rank the force fields in the descending order of similarity as C36s, C22cmap, OPLS-AA (due to compact state), and C22 Ã .
It might be argued that the computed Aβ10-40 conformational ensembles are specific to 330K. To investigate this possibility we used our REMD sampling to recompute the secondary structure propensities for five force fields at 300K. Fig D and Table A in S1 Text represent the analogues of Fig 2 and Table 1 obtained at 300K. It is seen that the secondary structure propensities at 300K and 330K are qualitatively similar, differring by no more than 7% for the propensities > 0.1. The same conclusion applies to the comparisons of the distributions of secondary structure along Aβ10-40 sequence, which demonstrate, for instance, that C22 Ã helical propensity < H(i) > only slighly increases at 300K. Finally, decrease in temperature triggers no qualitative changes in tertiary interactions (see S1 Text). Therefore, we surmize that Aβ10-40 conformational ensembles at 300K and 330K appear qualitatively similar.
Combining the analysis of secondary and tertiary structures we make the following conclusions: • All force fields predict that Aβ10-40 adopts unfolded structure dominated by turn and random coil conformations; • Water model does not dramatically affect secondary or tertiary Aβ peptide structure with standard TIP3P model favoring slightly more compact states; • Although there are no significant differences in secondary structures observed in C36 and C22cmap simulations, little similarity is seen in their tertiary interactions.
• Unique features of OPLS-AA force field are moderate β-structure propensity and extensive, but flickering long-range tertiary interactions leading to Aβ collapse. There are no common top tertiary interactions between C36 and OPLS-AA force fields.
• Unique features of C22 Ã are moderate helix propensity and multiple, exceptionally stable long-and short-range interactions. Based on RMSD computations applied to secondary (helix, turn, and backbone fluctuations distributions) and tertiary (contact maps) structure, we conclude that this force field differs the most from C36.
Which force field describes Aβ10-40 most accurately?
Using REMD simulations we have computed the 3 J HNHα -coupling and RDC constants for the amino-truncated peptide Aβ10-40 and compared them to the experimental measurements performed for the full-length peptide Aβ1-40. Although Aβ1-40 and Aβ10-40 peptides are not entirely identical having 78% of sequence homology, their comparison can still be instructive as demonstrated below. We first evaluate the agreement between in silico Aβ10-40 and experimental Aβ1-40 data in light of similar assessments made in the literature for Aβ1-40 or Aβ1-42 peptides, in which identical peptide species were used both in the experiments and simulations [17,27]. According to our analysis (Table 7) RMSD and PCC between the in silico and experimental J-coupling constants for the "best" C36s force field is 0.  [27].
A recent study has compared in silico and experimental J-coupling and RDC constants for three peptides, Aβ1-42, Aβ1-40, and Aβ1-42-M35ox using OPLS-AA/L force field and TIP3P water model [17]. Using the coefficients of Vuister and Bax [53], the RMSD values for J-coupling distributions were 1. 21 ). The values of PCC calculated by us are about the same or slightly lower than those reported for the three peptides (our PCCs of 0.44, 0.44 or 0.48 vs "their" PCC in the range of 0.4-0.6). As shown above these conclusions hold irrespective of computing the RMSD and PCC using all Karplus equation coefficient sets or only specific sets. Similarly, as measured by RMSD the agreement between in silico Aβ10-40 and experimental Aβ1-40 RDC data is much better than the previous comparisons, which involved identical in silico and experimental peptide species, such as Aβ1-42, Aβ1-40, or Aβ1-42-M35ox (our RMSD of 1.36 vs "their" RMSD of 1.49, 1.6, or 2.2 Hz). The same conclusion is supported by PCC comparing RDC distributions, which is 0.65 in our study against the approximate range of 0.35 to 0.45 reported previously.
Taken together the analysis above suggests two conclusions. First, if in silico Aβ10-40 and experimental Aβ1-40 J-coupling and RDC constants are generally in better agreement than these quantities computed and measured for strictly identical peptides, then the differences in the conformational ensembles of Aβ10-40 and Aβ1-40 are likely to be small or, at least, not exceeding the force field errors in reproducing the conformations of a specific peptide (Aβ1-42, Aβ1-40, or Aβ1-42-M35ox). Therefore, guided by the previous validations of protein force fields against Aβ NMR data [16,17,27,28] we argue that our analysis supports using Aβ10-40 peptide as a proxy of the full-length Aβ1-40. This conjecture has been made earlier by our [33,59] and other [60] groups. In this context, we note that the OPLS-AA/L simulations of the fulllength Aβ1-40 have predicted stable β-structure in Leu17-Ala21 and Ile31-Val36 regions [30]. In line with our view of Aβ10-40 as a proxy of Aβ1-40, the elevated β-structure propensity is observed in the same Aβ10-40 regions when sampled in our OPLS-AA simulations.
Second, a good agreement between in silico Aβ10-40 and experimental Aβ1-40 J-coupling and RDC constants argues that CHARMM36 force field with standard TIP3P water model is possibly the best force field for reproducing Aβ conformational ensemble. Previous studies evaluating the force fields for their ability to reproduce Aβ experimental data have identified OPLS-AA with TIP3P water model as most accurate [16]. However, to our knowledge CHARMM force fields have never been directly evaluated against the distributions of Aβ Jcoupling and RDC constants. In this study we addressed this issue. Recently, eight different force fields were evaluated using REMD simulations for their ability to reproduce small-angle X-ray scattering and NMR data for five natively unstructured peptides [61]. The authors have determined that CHARMM22 Ã generates the conformational ensembles most consistent with the experiments. They also noted erroneous CHARMM36 propensity to sample left-handed αhelix conformations. However, our study did not reach the same conclusions for Aβ peptides suggesting that the selection of the "best" force field still depends on the peptide and details of simulations. Incidentally, an updated version of CHARMM36 force field has been recently released, which corrects left-handed α-helix bias [62]. However, in the case of Aβ10-40 peptides this modification appears as not critically necessary given the lack of Aβ left-handed αhelix in the original CHARMM36 force field.

Conclusion
By applying REMD simulations we have performed comparative analysis of the conformational ensembles of amino-truncated Aβ10-40 peptide produced with five force fields, which combine four protein parameterizations (CHARMM36, CHARMM22 Ã , CHARMM22/cmap, and OPLS-AA) and two water models (standard and modified TIP3P). Aβ10-40 conformations were characterized by the analysis of secondary structure, backbone fluctuations, tertiary interactions, and radius of gyration. In addition, using computed conformational ensembles we have calculated Aβ10-40 3 J HNHα -coupling and RDC constants and compared them with their experimental counterparts obtained for the full-length Aβ1-40 peptide. Taken together, our study led us to several conclusions. First, all force fields predict that Aβ adopts unfolded structure dominated by turn and random coil conformations. Second, specific TIP3P water model does not dramatically affect secondary or tertiary Aβ10-40 peptide structure, albeit standard TIP3P model favors slightly more compact states. Third, although the secondary structures observed in CHARMM36 and CHARMM22/cmap simulations are qualitatively similar, their tertiary interactions show little consistency. Fourth, two force fields have unique features setting them apart from CHARMM36 or CHARMM22/cmap. Specifically, OPLS-AA reveals moderate β-structure propensity coupled with extensive, but weak long-range tertiary interactions leading to Aβ collapse. CHARMM22 Ã exhibits moderate helix propensity and generates multiple, exceptionally stable long-and short-range interactions. There are no common frequent tertiary interactions between CHARMM36 and OPLS-AA or CHARMM22 Ã force fields. Our investigation suggests that among all force fields CHARMM22 Ã differs the most from CHARMM36. Fifth, the analysis of 3 J HNHα -coupling and RDC constants based on CHARMM36 force field with standard TIP3P model led us to an unexpected finding that in silico Aβ10-40 and experimental Aβ1-40 constants are generally in better agreement than these quantities computed and measured for identical (100% homologous) peptides, such as Aβ1-40 or Aβ1-42. On the basis of this observation we argued that the differences in the conformational ensembles of Aβ10-40 and Aβ1-40 are likely to be small and the former can be used as proxy of the full-length peptide. We also concluded that CHARMM36 force field with standard TIP3P model produces the most accurate representation of Aβ10-40 conformational ensemble.