How accurate are your simulations? Effects of confined aqueous volume and AMBER FF99SB and CHARMM22/CMAP force field parameters on structural ensembles of intrinsically disordered proteins: Amyloid-β₄₂ in water

ABSTRACT

Amyloid-β₄₂ (Aβ₄₂) is an intrinsically disordered peptide intimately related to the pathogenesis of several neurodegenerative diseases. Molecular dynamics (MD) simulations are extensively utilized in the characterization of the structures and conformational dynamics of intrinsically disordered proteins (IDPs) including Aβ₄₂, with AMBER and CHARMM parameters being commonly used in these studies. Recently, comparison of the effects of force field parameters on the Aβ₄₂ structures has started to gain significant attention. In this study, the structures of Aβ₄₂ are simulated using AMBER FF99SB and CHARMM22/CMAP parameters via replica exchange MD simulations utilizing a widely used clustering algorithm. These analyses show that the structural properties (extent and positioning of the elements of secondary and tertiary structure), radius of gyration values, number and position of salt bridges are extremely dependent on the chosen force field parameters notably with the usage of clustering algorithms. For example, predicted secondary structure elements, which are of the great importance for better understanding of the molecular mechanisms of neurodegenerative diseases, deviate enormously in models generated using currently available force field parameters for proteins. Based on the derived models, chemical shift values are calculated and compared to the experimentally determined data. This comparison revealed that although both force field parameters yield results in agreement with experiments, the obtained structural properties were rather different using a clustering algorithm. In other words, these results show that the predicted structures depend heavily on the force field parameters. Importantly, since none of the force field parameters currently utilized in MD studies were developed specifically taking into account the disordered nature of IDPs, these findings clearly indicate that new force field parameters have to be developed for IDPs considering their rapid flexibility and dynamics with high amplitude. Furthermore, molecular simulations of IDPs are typically conducted using one water volume. We show that the confined aqueous volume impacts the predicted structural properties of Aβ₄₂ in water. Although up to date, confined aqueous volume effects have been ignored in the MD simulations of IDPs in water, our data indicate that these effects have to be taken into account in predicting the structural and thermodynamic properties of disordered proteins in solution.

KEYWORDS:

• force field parameters
• confined aqueous volume effect
• intrinsically disordered proteins
• molecular dynamics simulations

Introduction

Based on the century-old structure-function paradigm, proteins can perform their function when they have stable two- and three-dimensional structures. On the other hand, intrinsically disordered proteins (IDPs) and intrinsically disordered protein regions (IDPRs) lack stable structures, but perform important biological functions, often playing important roles in crucial biological processes, such as signaling, recognition, and regulation.^1,2 IDPs/IDPRs are commonly associated with the pathology of various maladies, such as amyloidoses, cancer, cardiovascular diseases diabetes, and neurodegenerative diseases.^3–7 Interactions with specific small molecules of organic and inorganic nature, or with other proteins, nucleic acids, and membranes can cause a specific folding at least in the part of an IDP.^8–10 In other words, the biological medium surrounding IDP not only influence its structural state, but can also trigger a folding process or promote misfolding. In addition, the presence of macromolecular species in the biological medium is believed to generate macromolecular crowding effects, which have been shown to modulate the conformations, functional properties, and aggregation kinetics of proteins.^11–18 Jana and Sengupta evaluated the adsorption behavior of small Aβ oligomers on the surface of a single walled carbon nanotube of high curvature.¹¹⁹ While the intrinsic self-assembly propensity of Aβ is markedly hindered by adsorption, the oligomeric units show high degrees of surface immobilization. Immobilized complexes are capable of oligomeric growth but with a shifted monomer-oligomer equilibrium compared to the free states. In the presence of an ionic solution and external electric fields, magnitudes of the current blockades were found to be sensitive to the oligomeric number of the adsorbed complex. In addition, we should mention here that Menon and Sengupta studied the impact of glucose crowding on Aβ dimerization using MD simulations.¹²⁰ However, they have not studied the impact of confined aqueous volume effects on glucose crowding. These crowders confine the aqueous volume surrounding target protein.^19–25 It is important to note that the confinement and crowding, although related, are still discernibly different phenomena due to the differences in the dispersion of available volume.^26–28 For example, a change in the confined aqueous volume around a protein results in a different concentration of that protein in its discrete environment. Therefore, effects associated with confined aqueous volume can influence the hydration properties of a protein, and, thus, can impact measured or modeled conformational changes and modulate the protein (mis)folding mechanisms.

The major step at the beginning of the fibrillation process of IDPs involves stabilization of some partially folded structures, instead of being initiated from the random coil-like highly unfolded species.²⁹ Therefore, conformational changes of an IDP leading to its partial folding are considered as crucial prerequisites for the fibrillation process.²⁹ One such disordered fibrillogenic and pathogenic peptide is the amyloid-β (Aβ) peptide. Aβ is naturally produced by the multistate proteolytic cleavage of the amyloid precursor protein (APP) as a 4 kDa peptide that exists predominantly in the forms of Aβ₄₀ and Aβ₄₂.^30–35 Among the various sequences derived from APP, Aβ₄₂ (see Scheme 1) has been proposed to be the principal component in parenchymal plaques in Alzheimer's disease (AD) patients (see refs.^30–35 and references therein). Drugs in clinical trials that decrease the levels of Aβ₄₂ have been proposed to be linked to a reduced risk of the AD development (see, for example, refs.^36–41). Various experimental studies have shown that Aβ₄₂ possesses an increased toxicity in comparison with other Aβ segments.^42–46 In addition, experimental analyses revealed that Aβ₄₂ fibrillate in a nucleation-dependent manner and can form amorphous aggregates, oligomers, and fibrils.^47–49

Scheme 1. Amino acid residue sequence of Aβ₄₂.

Experimental studies face challenges in the studies of Aβ monomers and oligomers due to the lack of a stable fold in this protein and its propensity to aggregate, as well as because of the crowding and solvation effects. Nuclear magnetic resonance (NMR) spectroscopy measurements presented little regular structure for the monomeric Aβ in solution, supporting its mostly disordered nature.^50–52 Parallel to several NMR measurements, region-specific endoprotease sensitivity studies revealed the presence of noticeable protection at the C-terminal region of Aβ₄₂.^53,54 Additionally, by studying the Ala21-Ala30 fragment of the full-length peptide it was shown that hydrophobic interactions might play crucial roles in the predictions of potential Aβ conformations and in the nucleation mechanism of the aggregation of Aβ peptides in solution. Specifically, it was established that the Ala21-Ala30 fragment adopts a turn structure due to the hydrophobic interactions and formation of salt bridges.^32,34,35 The tendency for hydrophobic interactions in water is readily understood in terms of the dependence of hydrophobic solvation on solute size. Experimental and theoretical studies investigating a short region of the peptide might not directly reflect the impact of the conformational changes of the full-length peptide on the hydrophobic interactions.³² In addition, it is known that the solution environment affects the secondary structure predisposition of Aβ. For example, helical structure was reported for Aβ in solutions that include water and SDS or fluorinated alcohols.⁵⁵ Furthermore, experiments have shown α-helix to β-sheet transition with an increase in the amount of water in a heterogeneous solution that includes hexafluoroisopropanol (HFIP).⁵⁶ The presence of HFIP in water may also effect the measured structural properties.

MD simulations have been conducted extensively for the structural analysis of IDPs (including Aβ₄₂) in solution (see refs.^57–61 and references therein). These studies complement experiments and provide information on the structural and thermodynamic properties at the molecular level with dynamics. However, recently, researchers started questioning the impact of different force field parameters on the structural properties predicted by MD simulations. Garcia et al. studied Aβ utilizing an explicit water model.⁶² Specifically, they compared Replica-Exchange Molecular Dynamics (REMD) simulation results with the outputs of the NMR experiments and found that the OPLS force field parameters along with the TIP3P model for water performed better than AMBER94, AMBER96 and GROMOS potential functions.⁶² In a recent study, these authors found that the combination of AMBER FF99SB along with TIP4P-Ew water model provides comparable accuracy to OPLS simulations utilizing the TIP3P model for water.⁶³ Based on these analyses, OPLS and AMBER FF99SB parameters were considered by Garcia and co-workers and Head-Gordon and co-workers as the best force field parameters for modeling IDPs.^62–64 Furthermore, most recently, Carballo-Pacheco and Strodel analyzed the outputs of using different force field parameters for Aβ₄₂.⁶⁵ They compared various REMD simulations with NMR experiments and tested AMBER FF99SB, AMBER99SB*ILDN, AMBER99SBBLDN-NMR, OPLS, and CHARMM22 parameters using 200 ns per replica REMD production runs. In the case of C_α chemical shifts, they observed the lowest RMSD values for CHARMM22*, whereas other force fields, including AMBER FF99SB, had comparable results, except AMBER99SBILDN-NMR, for which the simulated secondary C_α chemical shift values were much higher than the experimental values.⁶⁵ For C_β chemical shifts, these authors obtained the best results using AMBER FF99SB, while the worst results were obtained utilizing OPLS parameters. The accuracy of these force field parameters seems to depend on the choice of the solvent, while AMBER parameters yielded comparable results to OPLS parameters (in explicit water) utilizing an implicit water model. Most recently, we studied the impacts of AMBER 14FFSB and CHARMM22/CMAP parameters on the simulated structures of Aβ₄₂ using the modified TIP5P model for water.^35,66,67

Classical MD simulations of Aβ₄₂ using explicit instead of implicit water models have been conducted extensively. Specific water volumes have been utilized in these simulations (see, for example, refs.^68–72). One cannot exclude the possibility that the confined aqueous volume may impact the predicted structural properties and energetics of Aβ₄₂ in water.^66,67 To the best of our knowledge, the effects of the confined aqueous volume on the dynamic behavior and structural properties of Aβ₄₂ in aqueous media have not been studied before. In other words, current literature includes molecular dynamics studies using specific aqueous volume values without studying the impact of varying confined aqueous volumes on the predicted structures of Aβ₄₂. For instance, Man et al. tested most recently the impact of varying force field parameters on the structures of Aβ₄₂ using only one water volume utilizing an explicit water model (TIP3P) without considering confined aqueous volume effects.⁷³ Moreover, Stroedel and co-workers simulate the structures of Aβ₄₂ using only one water volume utilizing an explicit water model. In one of their recent studies, they used only 3,621 – 3,636 water molecules to solvate the structures of Aβ₄₂ and to study the impact of varying force field parameters on the simulated structures of Aβ₄₂.⁶⁵ The comparison of force field parameters is affected by the chosen water volume, i.e., confined aqueous volume effects. Furthermore, Baumketner and Shea predict the structures of Aβ₄₂ using only one water volume in their simulations.¹¹⁹

In this study, we compare two most commonly used force field parameters, AMBER FF99SB and CHARMM22/CMAP, utilizing extensive special sampling molecular dynamics simulations using a widely used clustering algorithm. In addition, we analyzed the effects of the confined aqueous volume on the predicted structures of Aβ₄₂ utilizing CHARMM22/CMAP parameters and an explicit water model. Our studies demonstrate that although both parameters yield results in agreement with experiments, the obtained structural properties of models are different from each other. It is important to emphasize here that none of these force field parameters were specifically developed taking into account the disordered nature of IDPs. Furthermore, we find that the chosen water volume (confined aqueous volume) affects the predicted structural properties of Aβ₄₂ in water. Therefore, MD simulation using only one specific volume for the solvent might lead to the misinterpretation of the structural properties of Aβ₄₂ (and other IDPs) in aqueous media. Our findings show that the amount of water around Aβ₄₂ impacts its structural properties and therewith its predicted aggregation propensities based on monomeric Aβ₄₂ structures.

Methods

We performed all-atom MD simulations coupled with thermodynamic calculations using different Aβ₄₂ concentrations in aqueous solution (varying confined aqueous volume) utilizing an explicit model for water that enables the simulations of inter-molecular hydrogen bonds between the solute and solvent molecules. Separate sets of long-time MD simulations (a total of 160 ns; and 60 ns production time after convergence) using different volumes of water, corresponding to monomeric Aβ₄₂ concentrations of 5.04 mM and 2.05 mM with a difference of 480 nm³ in the initial volume of water were performed. The experimental structure of the Aβ₄₂ peptide (PDB ID: 1Z0Q) was employed as an initial structure,⁷⁴ and all-atom MD simulations were performed with NAMD using the CHARMM22/CMAP parameters and the TIP3P model for water, respectively.^75–77 Most recently, the CHARMM22/CMAP parameters were shown to yield accurate results for intrinsically disordered proteins using the TIP3P water model but such studies were conducted using only one water volume without considering the effects of confined aqueous volume.⁶⁵ Therefore, we chose these force field parameters to gain insights onto the confined aqueous volume effects. CMAP corrections to the CHARMM22 force field parameters improve dynamical and structural properties of proteins in molecular dynamics simulations. For instance, the recently developed CMAP correction to the CHARMM22 force field was evaluated by MacKerell and co-workers using a model protein in their molecular dynamics simulations.¹¹⁵ They reported that substantial deviations from experimental backbone root mean-square fluctuations and N-H NMR order parameters obtained in the trajectories using CHARMM22 parameters (without CMAP corrections) are eliminated by the CMAP correction. For more details see Ref. 115 and references therein.

For modeling the impact of varying Aβ₄₂ monomer concentration, or varying confined aqueous volume, on the predicted structural and thermodynamic properties of the aqueous peptide, the Aβ₄₂ monomer was solvated in cubic boxes with layers of water of l_s = 20 Å and 30 Å, where l_s is the size of the water layer around the peptide (corresponding to aqueous volumes containing 10,314 and 25,777 explicit water molecules, respectively) in two separate simulations. The solvation process of monomeric Aβ₄₂ was performed in a way so that all water molecules overlapping with or located in close proximity (< 2.4 Å) to any solute atoms were deleted. The integration time step was set to 2 fs. An isobaric-isothermal ensemble was applied using Langevin dynamics.⁷⁸ The temperature was set to 310 K to correlate to the physiological temperature. The long-range interactions were treated using the Particle Mesh Ewald method, and the cutoff value for non-bonded interactions was set to 12 Å.⁷⁹ Usually a cutoff value between 8 Å – 12 Å is utilized in molecular simulations. For minimizing the effect of the chosen cutoff value, we placed the peptide (peptide center) back to the center of the box for avoiding any possible interactions between the peptide and its periodic images when the peptide attempted to do so. Therefore, we do not see any additional salt bridges formed between the peptide and its periodic images. However, we should mention here that the choice of the cutoff value is critical in molecular simulations. Kostiuk and co-workers studied the effects of cut-off distance used in MD simulations on fluids utilizing both NVT and NPT ensembles.¹¹⁶ They reported that the cutoff distance in the NVT ensemble plays a insignificant role in determining the equilibrium structure of the fluid when the ensemble has a high density. However, thermodynamic properties calculated using the NVT ensemble strongly depend on the chosen cutoff value. In the NPT ensemble, cutoff plays a role in determining fluid equilibrium structure and thermodynamic as well as kinetic properties. On the other hand, Beck et al. studied the impact of cutoff variations on polypeptides and specifically studied three spherical atom-based cutoff approaches for use with all-atom explicit solvent MD: abrupt truncation, a CHARMM-style electrostatic shift truncation, and their own force-shifted truncation.¹¹⁷ They compared the time-averaged helical content of an end-capped 17-residue alanine-based helical peptide with experiments. They also examined the effect of varying cutoff treatment and distance on energy conversation. They reported that the force shifted spherical cutoff method conserves energy, correctly predicts the experimental helical content, and shows convergence in simulation statistics as the cutoff is increased. Despite, Schreiber and Steinhauser reported in 1992 that the cutoff size strongly influences MD results on solvated peptides.¹¹⁸ Specifically, the behavior of a 17-residue model peptide was analyzed by means of MD simulations. Coulomb interactions were truncated for three values of the cutoff radius: 6, 10 and 14 Å. It was found that the stability of α-helix is a function of the cutoff size. Additionally, they reported that the helix is conserved in a simulation using a cutoff of 10 Å, it is lost within a very short period of 100 ps when the cutoff is increased to 14 Å. This finding demonstrated that the commonly used cutoff size of 10 Å is inappropriate because it does not ensure the convergence of Coulomb interactions. In parallel, they conducted simulations using the Ewald sum method. In contrast to the 14 Å cutoff trajectory, the Ewald technique simulation conserves the helical character of the peptide conformation. This demonstrates that even 14 Å is too short a cutoff. Due to the fundamental uncertainty introduced by the use of a simple cutoff, this truncation scheme seems questionable for molecular dynamics simulations of solvated biomolecules. The impact of varying cutoff values along with confined aqueous volume effects on Aβ₄₂ still remains to be investigated.

To maintain the system's neutrality, three Na⁺ ions were inserted in each simulation. Before performing detailed analyses of the trajectories gained from our MD simulations, we checked the convergence and assessed the physical relevance of our theoretical studies by comparing their outputs with the results of the NMR spectroscopy measurements, which provides the highest possible structural resolution for the Aβ₄₂ peptide. The corresponding experimental data were obtained from Dr. Michael Zagorski (CWRU). The calculated potential energy values and chemical shift values (see below) indicated that a reliable convergence was reached in both simulations (see Fig. 1). The chemical shift values were calculated using the DSSP program and compared to experimental values obtained from Dr. Michael Zagorski for further testing the convergence and assessing the physical relevance of our theoretical investigations. Fig. 1B and C present the comparison of the calculated and experimental chemical shift values for the Cα atoms in both systems. The correlation factors demonstrate the physical relevance of our simulation results. The minimal deviation of our theoretical data from the experimental values shows further the reached convergence in our simulations.

Figure 1. Calculated cumulative potential energy (U) for Aβ₄₂ in aqueous solution with a water layer of 20 Å (dashed line) and 30 Å (dotted line) for 160 ns, respectively (A); correlation of calculated and experimental Cα chemical shift values obtained from Dr. Michael Zagorski for Aβ₄₂ in an aqueous solution with a water layer of 20 Å (B) and 30 Å (C), respectively.

In general, the lack of corresponding microscopic analogue; i.e., a function of configurational space variables to be averaged to obtain accurate results, defines the main difficulty in free energy calculations of highly dense systems. Various methods have been developed for calculating the free energy values, such as insertion and deletion methods, and perturbation methods.⁷⁹ To study the conformational preference, following our previous studies, the Gibbs free energy change values were calculated using the Molecular Mechanics Poisson-Boltzmann/Surface Area (MM/PBSA) method (see equation 2(2) $G_{nonpolar}=0.00542\times SASA+0.92$(2) ), which enables the prediction of thermodynamic properties of large biocomplexes in dense solutions.^80–83(1) $G=E_{total}+G_{sol}-TS$(1)

G_sol is calculated from the sum of the electrostatic and nonpolar contributions. For the G_{solvation-electrostatic} calculations, the internal and external dielectric constant values were set to 1 and 80, respectively, and dipolar boundary conditions were applied.^80–83 Separate calculations on the proteins using 3000 linear iterations showed that the calculations were converged. The nonpolar Gibbs free energy contribution was approximated employing the solvent accessible surface area (SASA) based on the following relationship (equation 5).(2) $G_{nonpolar}=0.00542\times SASA+0.92$(2)

The normal mode analysis method was used to calculate the entropy.⁸² All trajectories gained from our MD simulations using an explicit model for water were utilized to calculate the thermodynamic properties described above. We also should mention here that the water tension is factored into the calculation of G_nonpolar.

To study the impact of AMBERFF99SB and CHARM22/CMAP parameters on the predicted structural and thermodynamic properties of Aβ₄₂ in an aqueous environment, we conducted separate sets of replica exchange molecular dynamics (REMD) simulations.^84,85 The NMR structure, PDB ID: 1Z0Q, of the wild-type Aβ₄₂ was used as an initial structure.⁷⁴ We should note here that the structure of Aβ₄₂has been used widely in molecular simulations. See for example, Refs.119–125. This structure, however, was measured in the presence of detergents and Aβ₄₂ prefers a helical conformation in the presence of detergents. Specifically, it was measured in a solution with 2,2,2-trifluoroethanol (TFE), which is an α-helix stabilizer. This experimentally determined helical structure was chosen to aid in understanding the impact of the confined aqueous volume and chosen force field parameters on the debated secondary structure transformations in the current literature including the α-helix to β-sheet transition. REMD simulations were performed for 4.8 μs using the AMBER FF99SB and CHARMM22/CMAP parameters. In addition, the physical relevance of the obtained data was evaluated by calculating the C^α and H^α chemical shift values for Aβ₄₂ and comparing these values to the experimental data. The correlation between the calculated (δ_sim.) and experimental (δ_exp.) values for the chemical shifts is large after convergence with a Pearson correlation coefficient of 0.98 and 0.93 for the C^α and H^α chemical shifts, respectively (see Figures S1 and S2 in Supporting Information section). These results support previous simulation studies, which indicated that 60 ns of simulation time are required to reach a convergence for Aβ₄₂.

The C^α, C^β, and H^α chemical shift values were calculated for the Aβ₄₂ models generated using different force field parameters (CHARMM22/CMAP and AMBER FF14SB) in an explicit water (TIP5P; Supporting Information section). For gaining more detailed insights into the impact of chosen force field parameters on the structures of the wild-type Aβ₄₂, we conducted additional 7.2 μs of REMD simulations utilizing our modified TIP5P water model. Figures S3 and S4 show the calculated α-helix and β-sheet contents in the structures of Aβ₄₂ using AMBER FF99SB and CHARMM22/CMAP parameters, respectively. The obtained results were then compared to the experimental data. To this end, the root mean square (RMSD) values were calculated for simulated and experimental values (see Figure S5 in Supporting Information section). Furthermore, we compared calculated intra-molecular interactions of Aβ₄₂ using the AMBER FF99SB and the CHARMM22/CMAP parameters utilizing the TIP5P model for water and clearly shows that the tertiary structures are affected by various force field parameters (see Supporting Information section).

The REMD simulations were conducted utilizing the AMBER10 software package.⁸⁶ The Onufriev-Bashford-Case implicit model was used in the implicit solvent simulations to represent the aqueous solution environment around the disordered peptide. This allowed avoiding the effects of confined aqueous volume and eliminated specific heat errors on the simulated solute structures (see ref.³² and references therein). Long range interactions were treated using the particle mesh Ewald method with a cut off value of 25 Å.⁷⁹ Temperature was controlled using Langevin dynamics with a collision frequency of 2 ps⁻¹.⁷⁹ In each simulation, we used 24 different replicas with temperatures distributed exponentially between 280 K and 400 K.⁸⁷ Suitability of these conditions was validated in our previous studies.

The integration time step of 2 fs was used in simulations, and trajectories were saved every 500 steps. At each replica, initial structures were equilibrated for 200 ps, and then the peptides were simulated for 100 ns per replica in different simulation sets. A total simulation time of 4,8 µs was achieved by exchange between replicas at every 5 ps. Exchange probability was set to 0.74 for the wild-type Aβ₄₂. After convergence, we conducted structural and thermodynamic calculations through the replica closest to physiological temperature (310 K). The molecular mechanics/generalized Poisson-Boltzmann surface area (MM/PBSA) method was used to calculate thermodynamic properties (see above). In the normal mode analysis, the potential function is expanded in a Taylor series about point x_o. If the potential gradient vanishes at this point and third and higher order derivatives are ignored, then the dynamics of the system can be described in terms of the directions (Q_i) and frequencies (ω_i) of normal mode, which satisfy:(6) $M^{-1/2}F{M}^{-1/2}{Q}_i={\omega }_i^2{Q}_{i}and{Q}_i{Q}_j={\delta }_{ij},$(6) where the matrix M diagonal contains atomic masses, and the Hessian F contains the second derivatives of the potential energy calculated at x_o.For a protein with N atoms, computation of normal modes involves numerical diagonalization of the 3N × 3N matrix.

The presence of the intra-molecular interactions was assumed when the centers of mass of two residues were within 9.0 Å of each other. The presence of a salt bridge was assumed when two hydrogen bonded atoms possessed opposite formal charges. The presence of a hydrogen bond was assumed when the distance between a donor and the acceptor atoms was ≤ 2.5 Å and when the hydrogen bond angle was > 113°.³² For calculating the prominence of secondary structure elements per residue with dynamics, we used the DSSP algorithm.⁸⁸ A summary of our simulation protocols is listed in Table 1.

Table 1. Protocols and standards used in our different simulations.

Simulation Type	Temperature (K)	Pressure (MPa)	Water Model	Potential functions for the Solute	Number of water molecules	Simulation Time	Convergence Time	Corresponding Simulation Results
Classical MD simulations using the NPT ensemble	310	0.1	explicit TIP3P	CHARMM22/CMAP	10314	160 ns	100 ns	Table 2
								Figure 1
								Figure 2
								Figure 3
								Figure 4
Classical MD simulations using the NPT ensemble	310	0.1	Explicit TIP3P	CHARMM22/CMAP	25777	160 ns	100 ns	Table 2
								Figure 1
								Figure 2
								Figure 3
								Figure 4
REMD simulations with special sampling using the NVT ensemble utilizing 24 replicas	310 (closest replica to 310 K)	0.1 (in average)	Implicit Onufriev-Bashford-Case	CHARMM22/CMAP	—	4.8 μs	60 ns	Figure 5
								Figure 6
								Figure 7
								Figure 8
								Figure 9
REMD simulations with special sampling using the NVT ensemble utilizing 24 replicas	310 (closest replica to 310 K)	0.1 (in average)	Implicit Onufriev Basford	AMBER FF99SB	—	4.8 μs	60 ns	Figure 5
								Figure 6
								Figure 7
								Figure 8
								Figure 9
REMD simulations using the NVT ensemble utilizing 32 replicas	310 (closest replica to 310 K)	0.1 (in average)	explicit TIP5P	CHARMM22/CMAP	10314	7.2 μs	60 ns	Figure S3
								Figure S4
								Figure S5
								Figure S6
REMD simulations using the NVT ensemble utilizing 32 replicas	310 (closest replica to 310 K)	0.1 (in average)	explicit TIP5P	AMBER FF99SB	25777	7.2 μs	60 ns	Figure S3
								Figure S4
								Figure S5
								Figure S6

Results and discussion

Effects of confined volume on structure and dynamics of monomeric Aβ₄₂

Classical MD simulations are conducted widely for the studies of IDPs in solution without considering the impact of the confined aqueous volume (see above). Based on our thermodynamic calculations, we find that the Aβ₄₂ monomer is able to form more favorable conformations in an aqueous solution with a larger confined aqueous volume, or a lower concentration of monomeric Aβ₄₂ (Fig. 2). The enthalpic contributions to the conformational preference (Gibbs free energy values) of monomeric Aβ₄₂ in both solutions are presented in Fig. 2A. As noted by the enthalpy and Gibbs free energy values, the conformational changes are spontaneous in both solutions with a lower concentration or a larger confined aqueous volume (l_s = 30 Å; l_s is the size of water layer around the peptide) and a higher concentration or a smaller confined aqueous volume (l_s = 20 Å) of monomeric Aβ₄₂ in water. However, the conformational change spontaneity through the enthalpic contribution (H) is more significant in the solution that includes a lower concentration of the monomer or a larger aqueous volume (l_s = 30 Å), since the negative enthalpy values are larger (Fig. 2A). These findings indicate that the energy required for the formation of the Aβ₄₂ conformations and the energy necessary for making the space available for these conformations of monomeric Aβ₄₂ are more preferred within a rather large available confined aqueous volume.

Figure 2. Enthalpic (A) and entropic (B) contributions to the Gibbs free energy values of monomeric Aβ₄₂ conformations in aqueous solution with Vs = 330 nm³ (blue) and Vs = 810 nm³ (red) where Vs presents the volume of water in the solution.

The entropic contribution (-TΔS) is also negative and contributes to the conformational spontaneity in both solutions of monomeric Aβ₄₂ (Fig. 2B). In parallel to the enthalpic contribution, the entropic contribution (-TΔS) to the Gibbs free energy is larger negative for the monomeric Aβ₄₂ solution with a larger confined water volume, or lower Aβ₄₂ monomer concentration in solution. This finding indicates that the disordered Aβ₄₂ peptide adopts more preferred conformational arrangements in larger confined aqueous volume (l_s = 30 Å). Therefore, Gibbs free energies and both enthalpic and entropic contributions present a larger spontaneity for the solution with a lower rather than a higher concentration of monomeric Aβ₄₂ in water or with a rather larger than a smaller confined aqueous volume. Overall, the differences in these thermodynamic properties suggest that the physical properties of monomeric Aβ₄₂ change with a difference in the confined aqueous volume. Therefore, the forces that delicately balance specific conformational preferences in monomeric Aβ₄₂ might also change depending on the confined water volume.

For gaining a deeper insight into the nature of these forces, detailed structural property analyses based on the free energy landscapes of the disordered monomeric Aβ₄₂ peptide in both solutions were performed (see below). Even though Aβ₄₂ is a disordered peptide with mostly random coil structure, the formations of α-helical and β-strand structures have been directly linked to the conformational favorability of the Aβ peptide as well as to the rigidity of its C-terminal region using experimental and theoretical tools (see, for example, refs.^{22,73,74,89–91}). Furthermore, β-strand formation, which was also detected in monomeric Aβ, was shown to play crucial roles in the oligomerization and fibrillation mechanisms of Aβ (see, for example, refs.^{8–10,32–35}). Whether a change in the size of the confined aqueous volume influences the formation, dynamics, and abundances of secondary structure elements in monomeric Aβ₄₂ needs to be addressed. Therefore, we calculated the difference in secondary structure formation per residue along with the abundances of the secondary structure elements for monomeric Aβ₄₂ in both solutions (l_s = 20 Å to l_s = 30 Å) (see Fig. 3). Based on these results, a change in the monomer concentration from high to low (or a change in the confined aqueous volume from small to large) starkly influences the formation and abundance of secondary structures per residue in part of the mid-domain and C-terminal regions of Aβ₄₂ (Fig. 3).

Figure 3. Differences in secondary structure abundances in monomeric Aβ₄₂ in aqueous solution with an increase in confined aqueous volume, corresponding to an increase from Vs = 330 nm³ to Vs = 810 nm³ (Vs is the volume of water).

The number of residues adopting β-sheet structure decreases in the C-terminal region with a decrease in monomeric Aβ₄₂ concentration in water, or the usage of a larger confined aqueous volume (l_s = 30 Å). We observe strong intra-molecular backbone hydrogen bonds in the C-terminal region between the residues Val36, Gly37, Val39, and Val40 in the solution with a smaller confined aqueous volume (l_s = 20 Å). Furthermore, these C-terminal residues interact with first shell water molecules via inter-molecular hydrogen bonds in an aqueous solution with a smaller monomeric Aβ₄₂ concentration (larger confined aqueous volume; l_s = 30 Å) and tend to adopt coil and α-helical conformations rather than β-sheet structure. Overall, these findings indicate that the residues of the monomeric Aβ₄₂ peptide are able to adopt different secondary structures with a change in the concentration of the monomer in water and with varying size of the confined aqueous volume. However, traditionally, MD of IDPs use one water volume to predict the general structural properties of aqueous disordered proteins. Our findings clearly show that such simulations using explicit water molecules and one water volume have to be questioned due to presence of significant confined aqueous volume effects.

Overall, our findings agree with those by Carballo-Pacheco and Strodel⁶⁵ as well as those by Derreumaux and co-workers⁷³ when we use a smaller confined aqueous volume for water. However, results change abruptly using a larger confined aqueous volume, indicating that the simulations of Aβ₄₂ in an aqueous solution environment utilizing an explicit water model need to be revised critically using varying confined aqueous volumes. Macromolecular crowding effects, even though indirectly related to the confined aqueous space (see above), have also been shown to have a differential effect on the secondary structure formation of the Aβ_10-35 fragment via MD simulations.⁹² Li and Mehler revealed that this peptide is able to adopt variable helical and β-sheet structures depending on the crowding effects; in simulations without crowding agents, Aβ_10-35 formed only helical structures, whereas simulations with crowders showed the β-sheet conformation, and in some cases it was even the only secondary structure that was formed.⁹² If we assume that the existence of crowders decreases or increases the confined aqueous volume around the Aβ₄₂ peptide, the predicted secondary structures for the peptide might change depending on the degree of changes in the confined aqueous volume around the peptide.

Not only the secondary structure elements and their corresponding abundances, but also the intra-molecular peptide interactions in solution that involve hydrogen bonds and hydrophobic interactions as well as salt bridge formations may vary with changes in the confined aqueous volume. Even though previous molecular simulations (using an explicit model for water) have been useful in determining the tertiary structures of Aβ monomers in water,¹⁸ the utilized confined aqueous volume is rather small (up to 2 nm) in these simulations and the impact of varying aqueous monomeric Aβ concentration on these generalized characteristics needs to be addressed (see, for example, ref.¹⁸).

Figure 4 presents the calculated intra-molecular peptide interaction maps for the monomeric Aβ₄₂ peptide in an aqueous solution with a change of the confined aqueous volume. A larger number of intra-peptide interactions between different regions of Aβ₄₂ occur in an aqueous solution with a higher concentration of the monomer or a smaller confined aqueous volume (Fig. 4A.I.; l_s = 20 Å). Specifically, mid-domain and N-terminal, N-terminal and C-terminal, and mid-domain and C-terminal interactions are more pronounced in water with a higher (Fig. 4A.I) rather than lower concentration of the monomer (Fig. 4B.I), or with a smaller (l_s = 20 Å) rather than a larger (l_s = 30 Å) confined aqueous volume. Furthermore, we note that the intra-molecular interactions within the hydrophobic region (Ala21-Gly30) are more abundant in an aqueous solution with a higher concentration of the monomer, or a smaller confined aqueous volume (l_s = 20 Å) as shown in Fig. 4B.I. Results of the previous simulation studies revealed that the N-terminal and mid-domain and C-terminal and mid-domain interactions serve as important general characteristics of the aqueous monomeric Aβ.^{9,16,22,28,35,57,62,93-105} Our findings for the monomeric Aβ₄₂ solution in a smaller confined aqueous space (l_s = 20 Å) support these earlier findings. However, our results clearly show that the predicted intra-molecular interactions strongly depend on the confined aqueous volume. Various intra-molecular interactions that occur with the usage of a smaller confined aqueous volume (l_s = 20Å) do not exist in a solution with a larger confined aqueous volume (l_s = 30 Å).

Figure 4. The intra-molecular peptide interaction (I), intra-molecular hydrogen bond (II), and hydrophobic interactions (III) maps for the Aβ₄₂ monomer in aqueous solution with Vs = 330 nm³ (A) and Vs = 810 nm³ (B) where Vs presents the volume of water.

The calculated intra-molecular hydrogen bond interactions and salt bridge formations in both solutions of the Aβ₄₂ peptide are presented in Fig. 4 and Table 2. We find that the number of residues forming intra-molecular hydrogen bonds is larger in the mid-domain and C-terminal regions of monomeric Aβ₄₂ in experiments with a smaller concentration of the monomer in aqueous solution and with a larger confined aqueous volume (l_s = 30 Å) as shown in Fig. 4A.II and 4B.II. Specifically, N-terminal residues, such as Glu3-Arg5, do not form intra-molecular hydrogen bonds with a larger confined aqueous volume (l_s = 30 Å) around the peptide but form inter-molecular hydrogen bonds with water molecules in the first hydration shell surrounding these N-terminal residues of the peptide (Fig. 4B.II). On the other hand, the intra-molecular peptide interactions via hydrogen bonds, such as the one between the residues Gly29 and Val40, occur with relatively high abundances within a larger confined aqueous volume (l_s = 30 Å) as shown in Fig. 4B.II.

Table 2. Formed salt bridges along with their abundances in monomeric Aβ₄₂ in aqueous solution with different confined aqueous volume effects, corresponding to water volumes of 330 nm³ and 810 nm³, respectively.

A combination of intra-molecular hydrogen bonds with electrostatic interactions forms salt bridges. We analyzed the salt bridges separately from the intra-molecular hydrogen bonds (see Table 1). Based on these results, the number of salt bridges is smaller in an aqueous solution with a larger confined aqueous volume (l_s = 30 Å). Salt bridges with high abundances are formed between the residues Glu3 and Arg5, Lys16 and Ala42, Asp23 and Lys28 in an aqueous solution with a smaller confined aqueous volume (l_s = 20 Å; larger monomeric Aβ₄₂ concentration). As shown in Table 1, the formation of salt bridges starkly depends on the available confined aqueous volume around the peptide. Experimental and theoretical results have indicated that the formation of a salt bridge between the residues Asp23 and Lys28 shifts the protein structure to a nucleation form that is prerequisite for aggregation.^{9,16,22,28,35,57,62,93–105} Therefore, the conformational changes elicited by the change in confined aqueous volume may result in the formation of disordered monomer conformations with a different propensity for aggregation.

Hydrophobic interactions are major driving factors in protein conformational changes, folding, misfolding, and aggregation mechanisms in solution. It is well known that hydrophobic interactions depend on many factors, including the solute size and the solution medium in which the solute exists.^66,67 Therefore, a change in the hydration phenomena through a variation in the confined aqueous volume in the solution of monomeric Aβ₄₂ may influence the hydrophobic interactions. For gaining deeper insights into the hydrophobic interactions in both solutions, we analyzed the hydrophobic interactions in both solutions with dynamics. The number of residues that form hydrophobic interactions and their corresponding abundances between the mid-domain and the C-terminal, within the mid-domain, and to some extent between the N-terminal and mid-domain regions of monomeric Aβ₄₂ are larger within a smaller confined aqueous volume (l_s = 20 Å) as shown in Fig. 4A.III and 4B.III. For instance, strong hydrophobic interactions between the residues Gly29 and Leu34 as well as Phe20 and Val36 occur in a solution with a larger monomeric Aβ₄₂ concentration, or in a solution with a smaller confined aqueous volume (l_s = 20 Å). Despite, the number of residues that form abundant hydrophobic interactions within the C-terminal region of the full-length peptide increase when a larger confined aqueous volume (l_s = 30 Å) is used. For example, Ile31 and Ala42 form strong hydrophobic interactions only in the solution with a smaller monomeric Aβ₄₂ concentration, or larger confined aqueous space. Overall, these results show that hydrophobic interactions depend on the concentration of the monomer in an aqueous solution and/or on the confined aqueous volume.

The calculated free energy landscapes as well as secondary and tertiary structures of the disordered Aβ₄₂ peptide in aqueous solution show that the outputs of measurements and simulations of these properties depend on the available confined aqueous space. In other words, hydration phenomena and associated water volume effect in a solution of a disordered peptide are able to impact the predicted chemical, physical, and biological characteristics of the peptide. Molecular simulations using only one water volume do not reflect the general structural properties and thermodynamic properties of intrinsically disordered proteins in water utilizing explicit water models.

Effects of AMBER FF99SB and CHARMM22/CMAP force field parameters on the structural properties of Aβ₄₂

AMBER FF99SB and CHARMM22/CMAP force field parameters are often used in the simulations of IDPs. As shown above, confined aqueous volume might have profound effects on such structural properties obtained in MD simulations. Therefore, we studied the impact of the chosen force field parameters (AMBER FF99SB and CHARMM22/CMAP) on the structural properties of Aβ₄₂ using an implicit model for water and the widely utilized RMSD clustering algorithm. In this approach, the root-mean-square deviations (RMSDs) between all structures are calculated and the cluster is identified as an ensemble containing largest number of conformations deviating from each other within 1 Å. These structures are removed from the list and a new search is performed to identify the second largest cluster. This procedure is repeated until all clusters counting 100 or more structures are found. We have tested the results presented below with a larger threshold of 1.2 Å and found them to be qualitatively the same. These findings are in agreement with earlier studies using the same general and widely-used clustering method.^{89–91,106–114} Representative structures and average radius of gyration values are depicted in Figure S7. We focus here on the first five cluster ensembles from the clustering algorithm. The convergence was tested in terms of α-helix content, random coil content, and β-turn structure content with time utilizing an implicit water model. Fig. 5 shows the convergence of these secondary structure elements. In agreement with previous studies,^8–10,32,62 we find that the convergence is reached within 60 ns of our simulation time. ^8–10,32,62 Fig. 5 shows the calculated secondary structure properties. In agreement with Strodel and co-workers,⁶⁵ we find that the secondary structure properties are affected by the choice of the force field parameters utilized in simulation (Fig. 5). We find that the 3₁₀-helix content is slightly higher (up to 3% – 8%) in Clusters 1–3 using the AMBER FF99SB parameters. A different trend is obtained for Clusters 4 and 5 with CHARMM22/CMAP parameters yielding a higher 3₁₀-helix abundance by 3% – 6%. These deviations become larger for α-helix (Fig. 5). Namely, the α-helix contents obtained using CHARMM22/CMAP parameters are larger than those obtained utilizing AMBER FF99SB parameters by up to 30% for all clusters of Aβ₄₂ ensembles in an aqueous medium. On the other hand, we find that the AMBER FF99SB parameters produce the largest content of the random coil structures in comparison to CHARMM22/CMAP parameters. The difference in random coil abundances varies between 10% and 22% (Fig. 5) using the AMBER FF99SB and CHARMM22/CMAP parameters, respectively. Furthermore, the turn structure content of Aβ₄₂ in an aqueous environment is – for all clusters – higher (by 10%-22%) when the AMBER FF99SB parameters. These results agree partially with those by Strodel and co-workers,⁶⁵ who used an explicit water model with 5.5 nm box length in their simulations. Please note the confined aqueous volume effects were not evaluated in these studies (see above). Furthermore, Derreumaux and co-workers reported the secondary structures for Aβ₄₂ using an explicit water model with only 21,858 atoms in their simulations.⁷³ They tested different force field parameters and found that the helix content was overestimated with the usage of AMBER FF14SB and CHARMM22/CMAP parameters. Again, their simulations suffer from the lack of the evaluation of the confined aqueous volume effects (see above).

Figure 5. Convergence of Aβ₄₂ simulations via REMD simulations tested using the secondary structure abundances with time.

The deviations in calculated radius of gyration values for the five clusters using the AMBER FF99SB and CHARMM22/CMAP (CHARMM27) are presented in Fig. 6. We find more compact structures using the AMBER FF99SB parameters (by 0.8 Å) in Cluster 1. The compactness of structures deviates by 0.2 Å using CHARMM22/CMAP and AMBER FF99SB parameters in Cluster 2. Same compactness is reached for cluster 3 using both parameters. More compact structures occur in cluster 4 using the AMBER FF99SB parameters (by 2.5 Å). A vice versa trend is obtained for the structures in cluster 5. Fig. 7 depicts the secondary structure elements and their abundances per residue with dynamics using AMBER FF99SB and CHARMM22/CMAP parameters in separate simulations. For the structures in cluster 1, we find that α-helix is abundantly present at His6-Tyr10 and Ala30-Val36 regions. More abundant 3₁₀-helix occurs at regions His6-Ser8 and Asn27-Val36 using CHARMM22/CMAP parameters (Fig. 7). β-Sheet structure formation is detected at Phe19, Ala21, Gly25, Ile31, Met35, Gly38, Ile41 using the AMBER FF99SB parameters that also detect more abundant turn structure formation at Arg5-Asp7, Glu11, Val12, His14-Phe20, Ile32 and Gly33. On the other hand, more abundant turn structure using the CHARMM22/CMAP parameters occurs at Ser8-Tyr10, Asn27-Ile31, Leu34-Gly37 in the structures located in cluster 1 (Fig. 7).

Figure 6. Radius of gyration values using AMBER FF99SB and CHARMM22/CMAP parameters in an implicit model of water. (A) Cluster 1, (B) Cluster 2, (C) Cluster 3, (D) Cluster 4, (E) Cluster 5.

Figure 7. Secondary structure properties simulated using AMBER FF99SB (Black) and CHARMM22/CMAP (Red) parameters in an implicit water environment; (A) Cluster 1, (B) Cluster 2, (C) Cluster 3, (D) Cluster 4, (E) Cluster 5.

More abundant α-helix exists in the N-terminal, mid-domain, and C-terminal regions of Aβ₄₂ in the structures located in cluster 2 when the CHARMM22/CMAP parameters are utilized. Highly abundant 3₁₀-helical structures occur at Glu3-His6, His14-Lys16, Asp23-Gly25, Leu34, and Met35 using AMBER FF99SB parameters in the structures assembled in cluster 2. However, more abundant 3₁₀-helix in the structures located in cluster 2 exists at regions Asp7-Tyr10 and Val36-Val40 utilizing the CHARMM/CMAP parameters (Fig. 7). For cluster 2 structures, we detect β-sheet structure at Ala21 and Gly37 and more prominent turn structure at His6, Asp7, Glu11, Val12, His14-Val18, Val24, Gly25, Ile32-Val36 using the AMBER FF99SB parameters. On the other hand, CHARMM22/CMAP parameters yield prominent turn structure formation in cluster 2 structures at Glu3-Arg5, Ser8, His13, Phe19, Phe20, Glu22, Asp23, Lys28, Gly29, Gly37-Val40.

For cluster 3, we observe more abundant α-helix at Glu3-His6 using AMBER FF99SB parameters. However, more abundant α-helix formation occurs at Asp7, Ser8, and Asn27-Val36 using the CHARMM/CMAP parameters. Utilizing AMBER FF99SB parameters, we detect prominent 3₁₀-helix formation at Ala2-Arg5, Glu22-Val24, Lys28-Ala30, Leu34-Val36 and β-sheet formation at Phe20, Ala21, Gly25-Lys28, Ala30 (Fig. 7). On the other hand, we note β-sheet formation at Phe4, Phe19, Gly25, and Gly37 using the CHARMM22/CMAP parameters. More prominent turn structure formation happens at Glu3, Phe4, Tyr10-Val12, His14-Lys16, Lys28-Ile31, Leu34-Val36 and Ile41 utilizing the AMBER FF99SB parameters for cluster 3 structures (Fig. 7). However, we detect more abundant turn structure at Arg5, His6, Ser8, Gly9, Glu22-Val24, Gly38, Val39 using the CHARMM22/CMAP parameters for the same cluster (Fig. 7).

For cluster 4, although more abundant α-helix structure occurs in the N-terminal and mid-domain regions using CHARMM22/CMAP parameters, more prominent α-helical structure occurs at Ala30-Leu34 utilizing the AMBER FF99SB parameters. Prominent 3₁₀-helix formation happens at Ser8-Lys16 using the AMBER FF99SB parameters in cluster 4 (Fig. 7). However, abundant 3₁₀-helix formation occurs in the C-terminal region utilizing the CHARMM22/CMAP parameters. We detect β-sheet formation at Val39 and Ala42 using the CHARMM22/CMAP parameters in cluster 4. More abundant turn structure formation exists in the N-terminal, mid-domain and C-terminal regions using the AMBER FF99SB parameters except for Ala30-Gly33 (Fig. 7).

For cluster 5, more abundant α-helix is detected in the N- and C-terminal regions using CHARMM22/CMAP. Prominent 3₁₀-helix occurs at Phe4-His6, Tyr10, Asp23, and Ile32-Leu34 utilizing the AMBER FF99SB parameters. However, more abundant 3₁₀-helix is formed at Glu11-His13, Val24-Ser26 using the CHARMM/CMAP parameters. β-Sheet occurs at Met35, Val36, Val40 and Ile41 using the AMBER FF99SB parameters, while only Val36 and Val40 form β-sheet structure utilizing the CHARMM22/CMAP parameters (Fig. 7). More prominent turn structure is formed at Glu3, Phe4, Ser8, Gly9, Val18-Phe20, Ser26 and Gly33 using the CHARMM22/CMAP parameters. However, more abundant turn structures exist at Arg5-Asp7, Tyr10-His13, Gln15-Leu17, Ala21-Gly25, Asn27-Ile31 and at Leu34 utilizing the AMBER FF99SB parameters (Fig. 7). Therefore, in agreement with Strodel and Derreumaux and their co-workers,^65,73 we find that the secondary structure properties are affected by the chosen force field parameters, however, we also find that this impact is larger with clustering techniques.

Figure 8 illustrates the simulated tertiary structure properties using both force field parameters. As depicted in Fig. 8, the calculated tertiary structure properties heavily depend on the choice of the force fields. For cluster 1, mid-domain interactions with the N- and C-terminal regions are more pronounced using the CHARMM22/CMAP parameters. However, interactions within the N- and C-terminus are more abundant utilizing the AMBER FF99SB parameters (Fig. 8). For cluster 2, N-terminal interactions with the mid-domain region are more pronounced using the AMBER FF99SB parameters. Central hydrophobic core (CHC) region interacts with the C-terminal region abundantly utilizing the AMBER FF99SB parameters. Furthermore, mid-domain interactions with the C-terminal region are more prominent utilizing the AMBER FF99SB parameters. For cluster 3, interactions between N-terminal and C-terminal regions are more prominent using the CHARMM22/CMAP parameters. Moreover, CHC region interacts with the mid-domain and C-terminal regions with a larger probability utilizing the CHARMM22/CMAP parameters. For cluster 4, interactions within the CHC region are more pronounced using the AMBER FF99SB parameters. N-terminal and C-terminal interactions occur utilizing the AMBER FF99SB parameters (Fig. 8). Using AMBER FF99SB parameters, we obtain more abundant interactions within the mid-domain region and between the mid-domain and C-terminal regions. For cluster 5, more frequent interactions occur between the C-terminal and N-terminal or mid-domain regions using the CHARMM22/CMAP parameters. CHC region interactions with the C-terminal region are more abundant utilizing the CHARMM22/CMAP parameters. However, CHC region interactions with the N-terminal region are more prominent using the AMBER FF99SB parameters (Fig. 8). Overall, these results show that the simulated tertiary structure properties depend on the choice of force field parameters. These finding are in agreement with the results of Strodel and co-workers and Derreaumaux and co-workers.^65,73 However, our results further demonstrate that conventional clustering techniques produce results that deviate from one another gigantically with the usage of different force field parameters.

Figure 8. Tertiary structure properties simulated using AMBER FF99SB and CHARMM22/CMAP parameters in an implicit water environment. (A) Cluster 1 using AMBER FF99SB parameters, (B) Cluster 1 using CHARMM22/CMAP parameters, (C) Cluster 2 using AMBER FF99SB parameters, (D) Cluster 2 using CHARMM22/CMAP parameters, (E) Cluster 3 using AMBER FF99SB parameters, (F) Cluster 3 using CHARMM22/CMAP parameters, (G) Cluster 4 using AMBER FF99SB parameters, (H) Cluster 4 using CHARMM22/CMAP parameters, (I) Cluster 5 using AMBER FF99SB parameters, (J) Cluster 5 using CHARMM22/CMAP parameters.

Furthermore, enormous deviations are noticeable for the salt bridge between Asp23 and Lys28 (Fig. 9), which is a critical component for the turn structure formation in the central hydrophobic core (CHC) region of Aβ₄₂. Although the highest peaks exist for the cluster 4 using both force field parameters, the peak location shifts to from 6 Å to 4 Å using CHARMM22/CMAP parameters (Fig. 9). A widely distributed broad curve is obtained for cluster 5 utilizing AMBER FF99SB parameters with a peak location at 9 Å, which shifts to 16.5 Å utilizing CHARMM22/CMAP parameters. Therefore, our results show that the probability distribution of the salt bridge between Asp23 and Lys28 is affected by the choice of force field parameters and that is effect is larger with clustering algorithms.

Figure 9. The calculated probability distribution of the distance between the Cγ atom of the Asp23 residue and the Nζ atom of the Lys28 residues for the Aβ₄₂ using AMBER FF99SB and CHARMM22/CMAP parameters in an implicit water environment. Simulations were conducted using AMBER FF99SB (A) and CHARMM22/CMAP parameters (B).

In addition to the studies mentioned above, Das et al. studied the impacts of A2T and A2V single point mutations on the structural ensembles of Aβ₄₂ in water by REMD simulations. They used only 5,600 water molecules to solvate the structure and utilized OPLS-AA parameters for the solute. Our parameters differ from theirs and a direct comparison cannot be provided at this stage but their studies lack the effects of confined aqueous volume. They also studied the impact of single point mutation on dimeric Aβ₄₂. The impact of confined aqueous volume effects on their predicted structural properties remains to be investigated. Moreover, Garcia and co-workers studied the effects of varying force field parameters on Aβ₄₀ (a different fragment size). They found that the OPLS-AA parameters with TIP3P water model, AMBER99sb-ILDN and TIP4P-Ew potential functions, and CHARMM22* and TIP3SP yield different structural properties for Aβ₄₀. Our peptide length, parameters and water models differ from theirs. However, we find also that varying force field parameters affect the simulated structural properties of Aβ₄₂. Confined aqueous volume effects were not studied by Garcia and co-workers. Therefore, we cannot provide a comparison to their studies.

Overall, traditional MD simulations suffer from severe limitations when applied for the analysis of structural and dynamic behavior of IDPs in water. Structural properties produced by the available force field parameters are very different. Furthermore, the usage of explicit water models in the simulations of IDPs in aqueous media are greatly impacted by the effects of confined aqueous volume. All in all, we are in need of new force field parameters and of dedicated hydration phenomena studies, which require new developments using physical models. Such basic research and development research activities have to be supported by currently active scientific agencies.

Conclusion

We have used all-atom MD simulations along with thermodynamic calculations – utilizing an explicit model for water – to identify whether and how a change of the confined aqueous volume impacts the free energy landscape and structural properties of the disordered fibrillogenic Aβ₄₂ peptide in aqueous solution. The insights obtained from the simulations lend support to some previous experimental investigations, and prove to be rather valuable in interpreting previously performed spectroscopic measurements.

Our Gibbs free energy calculations along with the analysis of enthalpic and entropic contributions suggest that the favorability of formed Aβ₄₂ conformations with dynamics in the solution depend on the confined aqueous volume. Specifically, we find that more favorable conformations of Aβ₄₂ are formed under the conditions of a larger aqueous volume. Another key finding from the simulations and calculations is the differences in formed structures. Significant differences in secondary structures occur in part of the mid-domain and C-terminal regions of the Aβ₄₂ peptide in an aqueous solution with varying confined volume. Based on our findings, C-terminal residues tend to adopt α-helical conformation with an increase of the confined aqueous volume, whereas abundant β-sheet structures are observed in the solution with a smaller confined aqueous volume.

Furthermore, predicted intra-molecular hydrogen bonds and hydrophobic interactions for monomeric Aβ₄₂ vary with a change in the confined aqueous volume around the peptide. Larger number of residues engaged in hydrophobic interactions are found in the solution with a smaller confined volume, and an increased number of residues that form abundant intra-molecular hydrogen bonds are seen in the solution with a larger confined volume. The formation of salt bridges depends also starkly on the confined aqueous volume around the Aβ₄₂ peptide. We note larger number of residues forming salt bridges in a solution with a smaller confined aqueous volume rather than the solution with a larger confined volume. To the best of our knowledge, this study presents the first systematic analysis of the effects of changes in the confined aqueous volume on the free energy landscape and structural properties of the full-length Aβ₄₂ peptide with dynamics.

Furthermore, the choice of force field parameters affects the predicted structural properties of IDPs in an aqueous solution medium. We tested two parameter sets, AMBER FF99SB and CHARMM22/CMAP parameters, and even though these parameters yield chemical shift values in agreement with the NMR experiments, the resulting conformational ensembles were characterized by different structural properties; i.e., secondary and tertiary structure properties, salt bridges and radius of gyration values.

Acknowledgment

This research was supported by an allocation and computing resources provided by the National Science Foundation (grant no. TG-CHE110044). The simulations were performed on Kraken at the National Institute for Computational Sciences. O.C. is grateful for the supports provided by the University of Texas at San Antonio and National Institute of Standards and Technology. The authors thank Michael G. Zagorski for the experimental data and helpful discussions. The authors also acknowledge Taizo Kitahara and Liang Xu for helpful discussions related to this project.