2. Why Do Molecular Dynamics Simulations Cannot Accurately Quantify the Aβ Structural Ensemble?
Experimental studies have been unable to determine the properties of Aβ peptide in solution due to the fast conformational changes and enhanced aggregation tendency. These studies have produced time- and space-average results that are difficult to map into a conformational state of folded and unfolded proteins. Computational simulations can make a time series at the atomic level that could help us explore the protein structure, dynamics, misfolding and aggregation mechanism, becoming a particularly suitable complement to experimental studies of conformational changes of Aβ. Several force fields (FFs) to study biomolecules have been developed in the last decades, such as AMBER, GROMOS, OPLS families, namely AMBER94, AMBER96, AMBER99, AMBER99SBildn, AMBER03, AMBER12SB, AMBER14SB, CHARMM22*, CHARMM36, CHARMM36m, OPLS, GROMOS43a1, GROMOS43a2, GROMOS43a3, GROMOS53a5, GROMOS53a6 and GROMOS54a7. Most of the existing FFs describe phenomena associated with well-structured proteins. However, Saravanan et al. 
concluded in a review study that the AMBER99SB-ILDN and CHARMM36m are highly optimized FFs and better choices for the characterization of IDPs such as Aβ peptide. This statement is supported because these FFs rendered the well agreement with experimental NMR chemical shift and β-sheet content, and the AMBER99SB-disp 
force field is also worth considering for the same purpose.
Five recent FFs Amber ff14SB, Amber ff14SB_idps, Amber ff99SB, CHARMM36, CHARMM36m have been used by Pawel et al. 
to explore the large conformational space of monomeric Aβ42 peptide during 10μs conventional molecular dynamics (MD) and 48 trajectories of replica exchange MD for 28.8μs. These FFs provided better results than their predecessor older versions. The potential energy can be described by Etotal = Ebonded + Enonbonded
where the bonded term (Ebonded
) consists of bond, angle, and dihedral-angle potentials, which explain the interactions of the atoms linked by covalent bonds, and the nonbonded term (Enonbonded
) is constituted by van der Waals(vdW) and electrostatic interactions. The electrostatic and vdW components are the primary contribution to nonbonded energy for monomeric Aβ1-42
. In the case of the CHARMM force field, the role of vdW interaction is reduced for Aβ1-42
peptide and enhanced for the Aβ1-42
-water-ions interaction, whereas, in the case of Amber ff99SB, nonbonded potential energy slightly level up by the higher domination of electrostatic interaction, resulting in additional stabilization of the Aβ1-42
peptide related an over-structured β sheet. The interaction with water molecules contributes to the dynamics, misfolded and self-assembly of the Aβ peptide. The stronger solute-solvent interaction leads Aβ1-42
to be less stable and more hydrophilic. In addition, MD simulation studies with CHARMM36m and FF14SB_IDPs show antiparallel β-sheets between residues 16–21 and 29–36 of monomeric Aβ1-42
, and short a β-strand in the C-terminal of the same monomer, which is in excellent agreement with NMR studies 
. AMBER_ff14SB and AMBER_ff99SB overestimated α-helical and β-contents, respectively. Pawel et al. 
strongly recommended using CHARMM36m force field for the study of the Aβ42-water-ion complex system over the AMBER FFs.
It is a big challenge to determine an accurate description of the structure of IDPs through MD simulations based only on FFs. In this perspective, Chong et al. 
reviewed advanced computational methods that employ protein configuration entropy and render a thermodynamic connection between structural disorder and protein properties. For example, the CHARMM and OPLS FFs exhibit lower average β-sheet content in dimers of Aβ42 than that obtained with GROMOS 53a6 force field 
. Subsequently, the average β-sheet content of the Aβ42 dimers was found to be greater in OPLSAA 
than in AMBERFF99SB 
. Interestingly, both AMBER99SB-ILDN and OPLS/L FFs have produced results of the average secondary structure of Aβ42 tetramer similar to each other 
The structural and thermodynamics properties of IDPs are susceptible to solute-solvent interaction compared to the folded protein. The choice of a reliable water model is necessary to characterize the Aβ42 peptide. Chong et al. 
performed an MD simulation to investigate the structural properties of Aβ1-42
peptide by employing AMBER ff99SB force field with different solvent water models. They demonstrated that TIP4P-Ew exposed more Aβ1-42
-water interaction than conventional TIP3P water model 
. They strongly encouraged using the TIP4P-Ew water model to investigate the Aβ1-42
peptide structural properties. In a review, Chong et al. 
reported that existing FFs are insufficient to expose Aβ protein to water. Recently, the same problem has been addressed by a couple of groups. The first group 
opted to scale the Lennard-Jones potential between atoms in proteins and oxygen atoms in water by factor 1.1 without disturbing water-water and water-protein interaction. The second group 
introduced a new water model, TIP4P-D, which included an additional parameter in the TIP4P water model to overcome the deficiencies in water dispersion interaction.
Recently developed FFs and their default water model are tabulated in Table 1
. Rahman et al. 
have evaluated the accuracy of recent developed FFs ff99IDPs, ff14IDPs, ff14IDPSFF, ff03w, CHARMM36m, and CHARMM22* by performing MD simulations for two short peptides (HEWL19 and RS), five IDPs (HIV-rev, Aβ40
, and pdE-γ) and two folded proteins (CspTm and ubiquitin) using trajectories of 1, 1.5, 5 or 10 μs for each system. They have compared J-coupling between MD simulation and NMR experiment for folded and disordered protein using different FFs. The J-coupling (J3-HNH2
and J3-HαC) parameter measures the secondary structure distribution based on φ backbone dihedral angle. Three IDPs FFs, ff99IDPs, ff14IDPSFF, ff14IDPs were in good agreement with the experimental J-coupling constant compared with tested FFs ff03w, CHARMM36m, and CHARMM22*. The balance between the local structural property (NMR chemical shift) and global structural property (Rg) is still a challenging issue for molecular simulations of IDPs. Rahman et al., 
noted two observations: (1) average Rg for Aβ421Z0Q
is 12.1 Å which is in close agreement with experimental Rg 12.4 Å, while ff03w showed Rg equal to 10.53 Å and CHARMM36m displayed Rg about 13 Å, which suggests highly divergence among those FFs; (2) the three IDPs FFs render a good balance between secondary structures contents for both Aβ40
. While Amberff03w, CHARMM36m and CHARMM22* overestimated the α-helical structure for IDPs, thus favouring folded protein structures. Therefore, the three specific IDPs FFs were developed by incorporating the changes made in the pre-existing FFs (Table 1
) to enable an accurate description of the folded and misfolded proteins 
Table 1. Latest developed force field for intrinsically disordered protein and water model.
||Updated from ff99SBildn by adding a set of backbone torsion parameters of eight disordered promoting amino acids.
||Wei Y et al.  2015
||Updated from ff14SB by embedding a set of backbone torsion parameters of eight disordered promoting amino acids.
||Song et al.  2017
||Updated from ff14SB by introducing a set of backbone torsion parameters for 20 amino acids
||Song et al.  2017
||Updated from ff03 by introducing a correction maps (CMAP)-optimized force field
of the TIP4P)
|Zhang et al.  2019
||Updated from ff99SB by improving the Accuracy of Protein Side Chain and Backbone Parameters
||Maier et al.  2015
||Updated from ff03 by adding slight backbone modification
||Best et al.  2010
||Update from a99SB-ILDN by an introducing small change in the protein and water vdW interaction terms
||Robustelli et al.  2018
||Updated from CHARMM36 by a refined backbone correction map potential
||Huang et al.  2017
||Updated from CHARMM36m by CMAP corrections made for all 20 naturally occurring amino acids
||Liu H et al.  2019
||Updated from CHARMM by introducing modifications in backbone torsion potential
||Stefano Piana et al.  2011
||Van der Waals interaction between protein and water are included in CHARMM36m
||Samantray et al.  2020
The inconsistency of empirical physical models used in MD techniques can impact FFs and water models that affect the simulation result’s accuracy. Researchers have strived hard to develop perfect FFs to improve IDPs description; they aim to describe the high flexibility of these proteins, thus enlarging the conformational ensemble and increase the possibility of locating them in different local minima. Mu et al. 
reported a couple of ideas to improve the accuracy of FFs for IDPs structural characterization, (1) Modification of force field parameters aided by global optimization, and (2) Maintaining a good balance between secondary structure via reparameterization (backbone dihedral parameters and vdW interaction between water-protein interaction) of existing FFs. One of the most common problems among the IDPs force field is over-stabilizing protein-protein interaction that impacts the aggregation mechanism of IDPs. Due to the IDPs force field’s inaccuracy, Mu et al. 
encouraged improving backbone dihedral parameters and Lennard-Jones potential parameter (protein-water interaction) in the existing IDPs FFs and obtained training data from experimental observation and quantum chemical calculation. Undoubtedly, both reparameterization and training strategies may assist in new FFs development.
Another option to the atomistic FFs is to employ state-of-the-art coarse-grained (CG) models that not only sample more efficiently the entire space of protein conformations for large systems but also allow simulations for longer time scales of hundreds of µs 
, generally forbidden by brute force all-atom MD simulation and very relevant for biological processes 
(e.g., folding, allosteric communications, conformational changes under mutations, self-assembly process, etc). Some popular CG FFs such as UNRES has been employed to study the fibril formations initiated by templates of Aβ40
capturing the dock-lock mechanism and similar the crowding effect of fragments was studied by PRIMO CG FF 
. MARTINI 2 was employed to unveil the aggregation and organization of short Aβ16-22 peptides in lipid membranes 
. The abovementioned CG FFs capture processes in time scales inaccessible atomic FFs, but yet they were restricted by system size or a lack of flexibility in secondary structure transitions. In this regard, the new release of MARTINI 3 force field aided by Gō-Like model 
can become a new tool for realistic exploration of full-length Aβ peptide aggregation in contact with complex lipid-cholesterol membranes. Marrink group has made large efforts to develop a library of different lipid species (i.e., about 63 types) consistent with human plasma cells 
. Following this idea, Poma et al. 
have employed a very simplified CG model 
to unveil the mechanical properties of Aβ40 and Aβ42 fibrils under different mechanical deformation process (e.g., tensile, indentation and shearing stresses) 
and more recently they investigated 
the change in mechanical stability between the oligomer and matured fibrils. The soluble oligomers are characterized by a length size about 3 to 5 nm, whereas Aβ fibrils typically reach hundreds of nm (see Figure 2
). It is well-known the inverse relationship between toxicity and the length of the Aβ assembly 
. Oligomers are considered more toxic than fibrils because of their high degree of flexibility toward low molecular weights, and the possibility of forming hydrophobic structures that may impair cell functions. Instead, fibrils are more thermodynamic stable and stiffer and less capable to undergo transitions to smaller and more toxic assemblies. Over the years, still some questions still remain open in terms of the mechanical characterization in amyloidogenesis and Aβ aggregate maturation: (a) what is the major role of the mechanical stability during Aβ oligomerization (e.g., tetramer, hexamer, etc)? (b) is the toxicity of oligomers which are closer to the membrane correlated by a minimum of mechanostability of Aβ peptide complexes, (c) how is the gain in mechanical stability from oligomers to fibrils involved in disease progression? and (d) Can we devise new strategies to reduce the mechanical stability of oligomer-to-fibril step through small molecular breakers (i.e., drugs recognition process)? To provide new answers to those open problems, new research combining versatile CG FFs and single molecule force spectroscopy is highly advised. For instance, in MD simulations, one can trace hydrophobic native interactions in equilibrium and simultaneously during a deformation process. Hence, the idea of hydrophobic structures can find support in molecular simulations as the main driving force for cell damage. Furthermore, molecular pathways could elucidate the critical conformation that maintains the mechanical stability of the Aβ assembly in the single-molecule force spectroscopy experiment. Figure 2
depicts the most relevant structure during aggregation and current methodologies used for its computational characterization 
Figure 2. Representation of Aβ structures during aggregation: (a) free Aβ peptide, (b) oligomers, (c) fibrils and (d) plaques. Typical lengths are given for peptide (<3 nm), oligomers (3–5 nm) and fibrils (20–100 nm). All-atom MD (AA-MD), coarse-grained MD (CG-MD) and continuum models have been employed to unveil the mechanical stability for each structure. Young modulus (Y) defined by the ratio of applied tensile stress to a given strain provides an idea of the elastic regime. Y values for each structure is taken from ref (50 and 51). Some images are modified with permission of BrightFocus Foundation.
Samantray et al. 
have examined recently developed IDPs FFs, namely AMBER99SB-disp, CHARMM36m, and CHARMM36 with enhanced protein-water interactions (CHARMM36mW) for the study of Aβ16-22
(wildtype) aggregation and its mutation F19L Aβ16−22
(mutation 1) and F19 V/F20 V Aβ16−22
(mutation2), as model systems for testing purpose. In AMBER99-disp, the peptide-water interactions are increased too much resulting in an inhibition of the Aβ16-22 aggregation. The same trend has been observed in the simulation with the AMBERFF03w force field 
. The difference between CHARMM36m and CHARMM36mW is a reparameterization of the protein-water interaction. In contrast, an experimental study 
reported the following aggregation order mutation1>wildtype>mutation2≈0. AMBER99-disp does not apply for Aβ aggregation process because the interactions between peptide and water are drastically increased, leading to inhibition of the aggregation pathway. The CHARMM36mW can provide aggregation rate in the order of wildtype>mutation1>mutation2. Thus, FFs cannot reproduce the aggregation of Aβ peptide observed in experiments, but they maintain a better balance between peptide-water and peptide-peptide interaction. These results imply that improving the force field significantly impacted the simulation aggregation pathway than modifying the protein sequence. Samantray et al. 
strongly encourage the use of CHARMM36mW for studying a full-length Aβ, even though it is not a perfect force field, since it yielded a promising result for aggregation benchmark. Nevertheless, reparameterization of this specific force field is still required.
Lockhart et al. 
performed REMD simulations to examine the impact of the three popular FFs, CHARMM22 (protein FF) with CHARMM36 (lipid FF), CHARMM36m (protein FF) with CHARMM36 (lipid FF), and Amber14SB (protein FF) with Lipid14 (lipid FF), for the binding mechanism between the Aβ10-40
and the Dimyristoylgylcerophosphocholine (DMPC) bilayers. These three FFs have shown similar results in subjects like (a) stable helix formed in C-terminal of the peptide, (b) C-terminal of the peptide inserted into the bilayer hydrophobic core, (c) the thickness of the bilayer induced by the peptide about 10 Å and d) the disordered effect induced by the peptide on the fatty acid tails in the DMPC lipids. Nevertheless, these three FFs yielded different conformation ensembles of the peptide and bilayer that do not disturb the binding of the peptide with the bilayers.
Coskuner et al. 
have reviewed several studies extensively and suggested that widely used FFs CHARMM, AMBER, GROMOS and OPLS, cannot produce accurate results for disordered entities. Even more, all existing computational techniques were designed to describe phenomena in ordered protein systems rather than in disordered protein structures. Important to mention that, Density Functional Theory (DFT) suffered from a number of errors that originated in the approximation of exchange-correlation functionals. These errors have been identified as the underestimation of barriers in describing chemical reactions, the band gaps, charge transfer excitation energies, and binding energies of charge transfer species in a biomolecule. DFT is the basis for constructing FFs such as CHARMM, AMBER, GROMOS and OPLS for intrinsically disordered protein and their complexes with ligands. These FFs are associated with an approximate exchange-correlation function that may lead to mistakes during prediction of the structural properties of IDPs. Notably, overcoming deficiencies of the exchange-correlation functional in DFT will help to improve the accuracy of IDP FFs.