The study of tau protein aggregation and interactions with other molecules or solvents using molecular dynamics simulations (MDs) is of interest to many researchers to propose new mechanism-based therapeutics for neurodegenerative diseases such as Alzheimer’s disease, Pick’s disease, chronic traumatic encephalopathy, and other tauopathies.
Tau is a microtubule (MT)-associated protein that is essential for microtubule stabilization 
. In the adult human brain, there are six different isoforms, of which 4R2N (2N4R), the longest, has 441 residues. This full-length human tau has a charged N-terminal region, a proline-rich region (spanning amino acid residues 208-244) followed by microtubule-binding repeats (MTBRs) designated R1 to R4 (R1 aa 244-274; R2 aa 275-305; R3 aa 306-336; and R4 aa 337-368), which bind to axonal microtubules under physiological conditions, and a C-terminal region 
. The accumulation of abnormal tau aggregates in neurons is an important pathological feature in several neurodegenerative diseases grouped under the term tauopathies. These include Alzheimer’s disease (AD), chronic traumatic encephalopathy (CTE), and frontotemporal dementias (FTD) such as Pick’s disease (PiD), corticobasal degeneration (CBD), etc. 
. Tau has been found to fibrillate in vitro in the presence of negatively charged co-factors such as heparin, DNA, or RNA 
. At least two amino acid sequences in the microtubule-binding domain (R1 to R4) are critical in the aggregation of tau 
. The hexapeptide segments 275
(PHF6*) and 306
(PHF6) are located in the MT-binding domain (MTBD) in peptides R2 and R3, respectively 
. K18 contains the MTBD (R1–R4) and it is located in the C-terminal region 
Alzheimer’s disease (AD) is the most common cause of dementia, accounting for 60–80% of all dementia cases. An estimated 6.5 million Americans aged 65 years and older are now living with AD 
. This number is expected to increase to 13.85 million by 2060. It was Alois Alzheimer who first described the disease using the brain of his patient, Auguste Deter, who died in 1906, five years after she was diagnosed with what is known as presenile dementia 
. The characteristic pathologies of AD are the accumulation of beta-amyloid (Aβ) plaques outside neurons and neurofibrillary tangles (NFTs), composed of the protein tau, inside neurons in the brain 
. Several investigators have shifted their focus from Aβ to tau as an alternative target of novel therapeutics 
, as previous studies have found a weak or even nonexistent correlation between β-amyloid plaques and cognitive decline in the symptomatic phase of dementia 
The R2 and R3 peptides in the MTBD have different properties than the R1 and R4 peptides, which is probably due to the presence of the β-structure driving hexapeptides PHF6* and PHF6 
. The major functions of the microtubule-binding domain are listed in Figure 1 
. Paired helical filaments (PHFs-tau) and straight filaments (SFs-tau) form the neurofibrillary tangles (NFTs), one of the pathological hallmarks of AD 
. The PHF-tau structure from an AD-diseased brain (PDB ID:5O3L) has a total of 11,360 atoms, while SFs-tau (PDB ID:5O3T) has a total of 5,570 atoms. There are a total of ten chains in PHFs-tau, each chain having an equal but opposite pair 
, and consisting of eight β-sheets (β1- β8) running along the length of the protofilament and encompassing amino acids 306-378 
. SFs form neurofibrillary tangles with a width of ∼150 Å in AD. Hybrid filaments of SFs and PHFs have been observed, implying that they have a similar C-shaped subunit but are different in structure 
Figure 1. The major functions of MTBD.
Tau proteins may misfold and aggregate, leading to the formation of NFTs in the brain. Studies have shown that microtubule-binding repeats undergo a conformational change when they bind to microtubules, which may affect their tendency to aggregate. In Ding and their colleagues’ study , using an all-atom discrete molecular dynamics (MDs) simulation, they observed that both repeats can aggregate into metastable β-sheet-rich dimers according to their respective conformational ensembles and dimerization kinetics, with R2 and R4 being highly comparable in this regard. The β-sheets between chains in R2 aggregates were driven by residues in the PHF6* regions, whereas β-hairpins predominated in the formation of R4 dimers. The greatest propensity for amyloid aggregation was observed in the R3 repeat. In R3, residues in Paired helical filaments the PHF6 regions rapidly self-assembled into intermolecular β-sheets and subsequently stimulated the formation of larger β-sheets by other residues. Both PHF6 in R3 and PHF6* in R2 generated a considerable number of intermolecular interactions and contributed significantly to the early aggregation of tau.
The tau fibril that correlates with Pick’s disease is also called narrow Pick’s filament (NPF). Pick’s disease (PiD) is a rare neurological disorder in which the so-called Pick’s bodies serve as diagnostic markers, similar to neurofibrillary tangles in AD 
. They consist exclusively of 3R tau isoforms (isoforms comprising repeats R1, R3, and R4 in the MTBD) (Figure 2
. The NPF contains the R1 repeat in the fibril core and does not have the steric zipper PHF6: 306
or PHF6*: 275
. As in AD 
, a fuzzy coat consisting of the disordered N- and C-terminal regions of tau surrounds the filament cores and is removed by pronase treatment 
Figure 2. Microtubule-binding repeat region primary sequence, showing color-coded repeats of R1 (pink), R2 (yellow), R3 (gray), R4 (violet), and R’ (teal). The narrow Pick’s filament (NPF) lacks the R2 repeat. The residues of interest are outlined in red. A strand from the NPF cryo-EM structure (PDB-ID: 6GX5), colored according to the repeats of the primary sequence. Residues of interest are outlined in red and shown as spheres. An asterisk indicates the strong curvature of the R1 repeat. The largest tau isoform hTau40 is represented by the residue indices.
The four most commonly used software packages for molecular dynamics simulations of tau–peptide–lipid membrane complexes such as the tau–lipid bilayer, are GROMACS 
, AMBER 
, CHARMM 
, and NAMD 
. GROMACS, for example, is a versatile package for performing molecular dynamics, i.e., simulating Newton’s equations of motion for systems with hundreds to millions of particles 
. The wide availability of experimentally determined protein structures from the Protein Data Bank (PDB), which can be used directly in whole-atom simulations (AA) or as templates in coarse-grained (CG) MDs, and the use of HPC (High-Performance Computing) 
, have contributed much to new relevant biological insights, especially in the study of protein membrane systems. The improvement of GPUs (Graphics Processing Units) hardware and better accessibility of software packages have generated tremendous interest in using GPUs for scientific computations 
. The Fast Multipole Method (FMM) was recently integrated into GROMACS as an alternative to Particle Mesh Ewald (PME), as a first step toward exascale 
. This highly efficient GPU (GPU-FMM) has the ability to enable efficient and scalable biomolecular protein–membrane complex simulations on future exascale supercomputers 
. In large-scale calculations, the most time-consuming computational steps are often those involving nonbonding interactions, especially electrostatic interactions due to their long-range nature and the complexity of the Coulomb potential. The PME algorithm is the usual method for dealing with such interactions. However, the PME algorithm still requires significant computational resources, and in systems with a large number of charged particles, electrostatic interactions can become a bottleneck in the simulation. Recent advances in using FMM in combination with GPU hardware have shown that it is possible to solve such problems and enable faster and more efficient simulations.
2. Multiscale Molecular Dynamics Simulations of Tau Filaments (NPF and PHF/SF) in a Solution and in the Surface of a Lipid Bilayer
A series of simulations has focused on elucidating the structures of stable tau fibrils and oligomers. Narrow Pick’s filament (NPFPS
) phosphorylated at three experimentally verified sites in the MTBD (S262, S324, and S356) and is used to study the effects of phosphorylation on the local conformation 
. In addition, NPFGG
, a double-mutant fibril system with mutations E264G and D358G, was selected to understand the role of E264 and D358 on the local conformations and to study the influence of the salt bridges that they form throughout the fibril architecture. The fibril systems were studied together with a wild-type fibril, NPFWT
(as a control), using a 1 μs-long conventional molecular dynamic. Calculations were performed using CHARMM36m FF with the TIP3P water model. The fibril was constructed with water and 150 mM NaCl to neutralize the charges and replicate the in vivo environment. To keep the temperature at 310 K, the velocities are assigned according to the Maxwell distribution. Their results show that the phosphorylated (NPFPS
) and mutant (NPFGG
) systems were found to diverge during simulation, indicating major changes in the fibril architecture compared to NPFWT
. The largest divergences in backbone conformation from the wild-type fibril system were observed during phosphorylation and mutation (Figure 3
) (although a 1 μs-long simulation is not sufficient to capture all fibril conformations). For the production run, assuming a constrained structure, the last 700 ns of the trajectory are used for further analysis 
Figure 3. The fibril systems’ root mean square deviation (RMSD) from 1 μs simulations.
The stability of the C-shaped structural motifs of the tau protein was studied by Nussinov and et al. 
. They used an all-atom explicit solvent MD calculation with the CHARMM FF and the TIP3P water model. They also included ~150 mM NaCl to represent the physiological ion concentration. All systems were heated from 0 to 300 K and held at equilibrium for 4 ns, except for K18, which was heated to 310 K. Their simulations showed that only the third and fourth repeat domains (R3-R4) retain the C-shaped conformation after optimization, while the first and second repeat domains (R1-R2) adopted a linear shape (Figure 4
). The β2-β3 and β6-β7 angles of PHF34
remained relatively the same during the simulations, but the β6-β7 angles of PHF13
increased, indicating that the structure tends to elongate. PHF23
adopted a V-shaped conformation, whereas a U-shape was preferred for PHF13
. The stability resulted from the interaction between the C-terminal residues and the N-terminus of the adjacent chain. Of all the PHFs studied, PHF34
exhibited the strongest intra-chain interactions.
Figure 4. (a) Initial structure of one strand of protofilaments and two strands of filaments of the R3-R4 combination. The RMSD (b), H-bond, (c) and RMSF (d) trajectories for the four PHF protofilaments. All protofilament’s time evolution of β6-β7 (e) and β2-3 (f) angles. (g) The PHF34 Q351-I371; PHF12 Q288-I308; and PHF13 and PHF23 T319-V339 distances. (h) PHF34 and PHF12 snapshots (t = 100 ns).
Bhargava et al. studied the tau protein on lipid bilayers by CG MD and AA MD simulations 
. The tau-SF structure (PDB ID:5O3T) was modeled with CG using Martini representations 
to model the lipids and proteins. Water and 0.15 M KCl were added when calculating AA MD. The Nosé–Hoover thermostat was used to maintain the temperature at 310 K. CHARMM-GUI was used to generate the initial configurations for the fibrils and lipids. The simulations show that the tau proteins interact differently with the zwitterionic compared to the charged lipid membranes. The negatively charged POPG lipid membranes increase the binding tendency of tau fibrils. Fourteen systems consisting of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylethanolamine (POPE), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylglycerol (POPG), and cholesterol (CHOL) in seven different compositions were simulated. Throughout the study, the symmetric composition of lipids in the upper and lower leaflets of the bilayer was used. In the case of the pure POPC system, 255 POPC molecules were randomly placed. Using the clustering algorithm, it was shown that the tau fibrils have different modes of interaction with the lipids. They found that the binding of the tau fibrils to the lipid bilayers coincides with the loss of β-sheet zones over the tau filament, which breaks the lipid bilayers. In another CG model, multiple atoms in proteins and lipids are approximated as a single bead and four water molecules are treated as a single particle (known as one bead 4:1 mapping) 
. The beads can be distinguished by their polarity or hydrophilicity