Digital Commons @ Trinity Digital Commons @ Trinity

: The self-aggregation of tau, a microtubule-binding protein, has been linked to the onset of Alzheimer’s Disease. Recent studies indicate that the disordered tau aggregates, or oligomers, are more toxic than the ordered ﬁbrils found in the intracellular neuroﬁbrillary tangles of tau. At present, details of tau oligomer interactions with lipid rafts, a model of neuronal membranes, are not known. Using molecular dynamics simulations, the lipid-binding events, membrane-damage, and protein folding of tau oligomers on various lipid raft surfaces were investigated. Tau oligomers preferred to bind to the boundary domains (Lod) created by the coexisting liquid-ordered (Lo) and liquid-disordered (Ld) domains in the lipid rafts. Additionally, stronger binding of tau oligomers to the ganglioside (GM1) and phosphatidylserine (PS) domains, and subsequent protein-induced lipid chain order disruption and beta-sheet formation were detected. Our results suggest that GM1 and PS domains, located exclusively in the outer and inner leaﬂets, respectively, of the neuronal membranes, are speciﬁc membrane domain targets, whereas the Lod domains are non-speciﬁc targets, of tau oligomers binding to neurons. The molecular details of these speciﬁc and non-speciﬁc tau bindings to lipid rafts may provide new insights into understanding membrane-associated tauopathies leading to Alzheimer’s Disease.


Introduction
The misfolding and self-aggregation of tau, an intracellular microtubule-binding protein, has been linked to the progression of various neurodegenerative diseases, e.g., Alzheimer's Disease [1,2]. The presence of intracellular neurofibrillary tangles (NFTs) consisting of microscopic, highly ordered and beta-sheet rich tau fibrils [3] has been established as a histopathological marker in tauopathies [1]. However, recent studies have demonstrated that the nanoscopic and highly disordered tau oligomers are more neurotoxic than the mature tau fibrils, and that the interactions of these toxic tau oligomers with neuronal membranes may play a key role in the onset of tauopathies [4]. In addition, these tau oligomers, originally created in the cytosol, may exit the neurons and translocate to other areas of the brain, including the space between the synaptically associated neurons [5]. Therefore, a molecular level understanding of tau-membrane interactions related to tau oligomers binding to the cytosolic (inner) and extracellular (outer) membrane leaflets of neurons will provide useful insights in understanding the early onset of tauopathies [1,2,6].
As a class of intrinsically disordered protein (IDP), the molecular mechanisms of misfolding and aggregation events of tau oligomers in solution and their associations with membranes are challenging to study computationally and experimentally due to the highly dynamic and heterogeneous structures of IDP [7]. Recent experimental studies have established that the lateral (along the membrane surface) and transbilayer (across the lipid The construction of three laterally and transversely heterogeneous, phase-separated lipid raft systems, CO-raft, GM-raft, and PS-raft were created using CG MD simulations [14]. Our raft systems contain saturated phosphatidylcholine (PC), unsaturated PC, cholesterol, monosialotetrahexosylganglioside GM1, and PS lipids, and their design was based on a fully hydrated and equilibrated 3-component coarse-grained (CG) lipid raft model [15]. Details of constructing the two rafts, CO-raft and GM-raft, are described elsewhere [16,17]. The construction of the third raft, PS-raft, was very similar to that of GM-raft in our previous studies [16,17]. Briefly, the CO-raft contains 828 saturated dipalmitoyl-PC (DPPC), 540 un- Previous studies [15][16][17] have established that the CO-raft contains ordered DPPCrich and CHOL-rich (Lo) domains, disordered DLPC-rich (Ld) domains, DPPC-DLPC (Lod) domains. In the presence of the asymmetrically distributed GM1 in the GM-raft, GM1-clusters were formed on one leaflet of the bilayer, with Lo, Ld, and Lod domains on both leaflets [16,17]. Similarly, the new PS-raft in this study also exhibited asymmetric transbilayer distribution of PS-clusters located on one leaflet, and Lo, Ld, and Lod were present in both leaflets as well. After 15 µs of CG MD simulations, the symmetrically distributed Lo, Ld, and Lod domains as well as the asymmetrically distributed GM-cluster and PS-cluster domains were still present, and these lipid rafts were used as the initial membrane structures for the protein binding studies.

Tau-K18 Oligomers
Wild-type and mutated Tau-K18 monomers were constructed using MD simulation procedures. Tau-K18 is a truncated Human Tau peptide containing four microtubule binding repeats [6]. The creation of an all-atom (AA) 130-residue long Tau-K18 monomer was based on a cryoEM-derived, fibrillar tau pentamer structure [3] as shown in Figure 1A.
Here, a 73-residue long peptide Tau 308-372 was first extracted from the chain A of the pentamer structure. Thereafter, a 57-residue long random coil Tau 243-307 was attached to the N-terminus of Tau 308-372 using a homology modeling algorithm [18] to create the final 130 residue-long Tau 243-372 , or the wild-type Tau-K18 (WT-K18) monomer. A membranebinding deficient (MBD) mutant [13] of Tau-K18, or MBD-K18, was created by mutating three highly hydrophobic residues at V287, I308, and V318 to three negative residues at E287, E308, and E318. The mutations at V287E, I308E, and V318E are highlighted in both the primary and AA structures of WT-K18 and MBD-K18 ( Figure 1A-C).
A hydropathy index vs. residue number plot [19,20] was also used to characterize the hydrophobicity profile of the amino acids of the two K18 monomers. The hydropathy plot revealed significant drops in the hydropathy indices from~0 to −3, 2 to −2, and 0 to −3 at the three mutation sites, respectively, as shown in Figure 1E. As shown in Figure 1A, the amino acid charge counts of WT-K18 are +21e and −11e, with a net charge of +10e, while the counts of MBD-K18 are +21e and −14e, with a net charge of +7e. Therefore, WT-K18 also carries more positive charge than MBD-K18. Supplementary Materials Figure SA1 further demonstrates the R1, R2, R3 and R4 repeats of WT-K18 and MBD-K18. A close examination of the primary sequences reveals that the first mutation V287E is at the middle of R2 and the other two mutations, I308E and V281E, are near the beginning and the middle of R3 repeat.
CG monomers were created from the corresponding AA structures from above using the AA-to-CG resolution transformation program martinize.py [21] based on the MARTINI-2.20 CG force fields [22]. Since no secondary structures were present in the AA monomer, no external dynamic network or constraints were imposed on the atoms of the monomers. The CG monomers were solvated in 0.1 M NaCl, and underwent energy minimization and position-restraints to reduce the local high energy structure created during the solvation steps. Finally, the monomer in solution underwent a 5 µs-long CG MD simulation in the NPT ensemble using the Martini CG force fields [22] and ran on the GROMACS-5.1.5 MD simulation program [23]. Detailed CG-simulation procedures for simulating disordered monomers in solution can be found [16,17]. The final 5 µs monomer structures of WT-K18 and MBD-K18 were used to create dimer and tetramer. The chain A (blue ribbon) from the CryoEM all-atom (AA) tau pentamer structure (A) was used as a template (underlined in the primary sequence) to build the AA wild type K18 or WT-K18 (B) and mutated K18, MBD-K18 (C) with the three mutation sites identified. The coarse-grained (CG) structures, yellow beads for side chains and pink beads for backbones, of both monomers are overlaid to the AA structures (in licorice). The backbone structures of monomeric K18 peptide after 5 μs of CG simulations in 0.1 M NaCl, with water shown as small dot, Na + in dark blue sphere, Cl − in light blue sphere, and the residues before (blue) and after (red) mutations are shown (D). The hydropathy plot of K18 monomers, i.e., hydropathy index vs. residue or sequence number with a 5-point moving average fit, for WT-K18 (blue) and MBD-K18 (red), is given (E). The self-aggregation process of four WT-K18 monomers, with each chain in a different color (blue for chain A, red for chain B, gray for chain C, and orange for chain D) within 3 μs of CG simulation in 0.1 M NaCl is demonstrated (F). The chain A (blue ribbon) from the CryoEM all-atom (AA) tau pentamer structure (A) was used as a template (underlined in the primary sequence) to build the AA wild type K18 or WT-K18 (B) and mutated K18, MBD-K18 (C) with the three mutation sites identified. The coarse-grained (CG) structures, yellow beads for side chains and pink beads for backbones, of both monomers are overlaid to the AA structures (in licorice). The backbone structures of monomeric K18 peptide after 5 µs of CG simulations in 0.1 M NaCl, with water shown as small dot, Na + in dark blue sphere, Cl − in light blue sphere, and the residues before (blue) and after (red) mutations are shown (D). The hydropathy plot of K18 monomers, i.e., hydropathy index vs. residue or sequence number with a 5-point moving average fit, for WT-K18 (blue) and MBD-K18 (red), is given (E). The self-aggregation process of four WT-K18 monomers, with each chain in a different color (blue for chain A, red for chain B, gray for chain C, and orange for chain D) within 3 µs of CG simulation in 0.1 M NaCl is demonstrated (F). To generate the self-aggregated CG Tau-K18 dimer and tetramer, the above monomers were self-replicated along different x-, yand z-directions using the molecule replication tool genconf of GROMACS [23], followed by the same MD simulation procedure as above to model the self-aggregation process. The procedure was identical to the self-aggregation process published recently [17]. Briefly, a dimer was created by replicating the simulation box containing a single monomer, or monomer-self-replication, along the x-direction, to generate two separated monomers lined up along the x-direction in a larger simulation box (see Figure SC1). Similarly, a tetramer was generated by monomer-self-replication along three directions, x, y and diagonal x-y directions, as shown in Figure 1F.
To evaluate the protein conformations of oligomers, 3D protein residue-contact maps describing the color-coded residue-residue minimum distances, among all protein residues, along the xand y-axes, were calculated using the g_mdmat tool of GROMACS [23] and a statistical and molecular interaction analysis tool, CONAN [24]. Time-averaged protein residue-contact maps with standard deviation were generated. Details of residue-contact map generation were described in our previous studies [17,25].

CG Simulations of Oligomer-Lipid Raft Complexes
Incorporation of the CG oligomer (monomer, dimer, or tetramer) into the equilibrated CO-raft, GM-raft, and PS-raft was accomplished by placing each oligomer above the center of the upper lipid leaflet of each lipid raft using the genbox tool of GROMACS [23]. Some solvent atoms were removed to accommodate the newly added protein. In this study, three independent simulation replicates corresponding to three different initial positions of each oligomer placed above the surface of the lipid raft were created. The first replicate was placed at 5 nm above the center of the lipid leaflet, and the other two replicates were subsequently created with the protein position shifted ±2 nm along the x-direction relative to the position of the first replicate.
The use of three independent simulation replications for each oligomer-raft complex is important to improve phase sampling of protein-membrane binding events. Upon generating the initial structures of the three replicates, the same CG MD simulation procedures for oligomer simulations were performed for up to 15 µs or longer.
In this study, a total of 54 simulation replicates involving oligomers of 3 aggregation sizes, 3 lipid rafts, and 3 simulation replicates for each oligomer-lipid raft complex with an accumulated CG simulation time over 800 µs were generated. Since the CG simulation time based on the Martini force field is roughly four times faster than the real time based on the CG diffusion rate of water [15,22], this study may represent more than 3 milliseconds of biological processes that are relevant to the early aggregation of proteins on membrane surfaces.

AA Simulations of Oligomer-Lipid Raft Complexes
After the CG simulations of oligomer to lipid raft binding, each replicate of the CG tetramer raft system was converted to the AA structure using a CG-to-AA resolution transformation program, backward.py [26]. The newly transformed AA structure was equilibrated with similar energy minimization and position restraining procedures as in the simulations of the CG oligomer-raft complexes. However, instead of the Martini CG force fields, the atomistic AMBER99SB [27] for proteins and SLIPID [28][29][30][31] for lipids force fields were used in all AA simulations. and DLPC molecules that are not included in the Ld domain. Finally, Lo-CHOL, Ld-CHOL, and Lod-CHOL represent distinct groups of CHOL, for which at least one CHOL atom is within 0.5 nm of the PC lipid atoms in Lo, Ld, and Lod domains, respectively. The same g_select tool as above was also used to characterize annular lipid shells from each oligomer/raft system. Here, if an atom of any lipid (DPPC, DLPC, CHOL, or GM1) is within a certain threshold, e.g., 0.5 or 1.2 nm from an atom of an oligomer, that lipid is assigned to the 0.5 nm or 1.2 nm annular lipid shell, accordingly.

Characterizations of Phase-Separated Lipid Domains and Annular Lipids
The time-and replicate-averaged number of lipids in each lipid domain or annular lipid shell over the last 5 µs for CG and 50 ns for AA simulation were performed to assess the domain or annular shell compositions of our raft systems in the absence and presence of membrane-bound K18 oligomers. Additionally, the lipid orientational order parameter, a measurement of the tilt of three sequentially connected carbon atoms along the lipid acyl chains of DPPC, DLPC, POPS and GM1 with respect to the normal of the bilayer, was calculated as a function of the carbon number of lipid acyl chain. Detailed steps of characterizations of composition and lipid order parameters of the lipid domains and annular lipids are described elsewhere [16,17,25].

Membrane Binding Kinetics, Domains, and Energetics of Oligomers
Examination and visualization of the binding kinetics related to the time events of protein-membrane binding and the lipid types surrounding the membrane-bound oligomers for each replicate were performed using a molecular visualization program, VMD [32].
Statistical analysis of the lipid-binding kinetics and locations of the protein binding were obtained using the minimum-distance analysis tool, g_mindist, from GROMACS [23]. Briefly, a protein-lipid or protein-water minimum distance (mindist), defined as the minimum distance between any protein atom and the atom of its binding lipid or water neighbors, was determined as a function of simulation time. In addition, the number of contacts of the mindist within an interaction threshold (2 nm) vs. simulation time was also determined. Finally, the time-averaged mindist vs. protein residue number over the last 5 µs of the 15 µs long CG simulation or the last 50 ns of the 200 ns long AA simulation was calculated. These three parameters, mindist vs. time (upper panel), number of contacts vs. time (mid panel), and mindist vs. residue number (bottom panel) were expressed as a 3-panel plot of each replicate. The first two panels were used to determine the kinetics of protein binding in terms of the time event of protein attachment, or lipid binding time, to each lipid type. The last plot provides important information about the mindist of the nearest-neighbor lipids or water surrounding the protein upon forming a membrane-bound state.
These mindist analysis plots and the VMD analyses provide the basis for determining the lipid kinetics and binding regions of K18 oligomers and for establishing the membrane domain target of the membrane-bound oligomer in the model raft.
Similarly, the time-averaged and replicate-averaged nonbonded interaction energies between the oligomer and each lipid group, or protein-lipid binding energy were collected using the tool, g_energy, from GROMACS [23] for each oligomer/raft complex. Both nonbonded Van Der Waals or Lennard-Jones and electrostatic or Coulomb energies were separately calculated from CG and AA simulations of our protein/raft systems. Details of the use of the above analysis tools for studying protein/raft binding kinetics, behaviors and energetics are described elsewhere [16,17].

Secondary Structure Determination
The secondary structure of each amino acid of the protein at every time-step was calculated using a tool, do_dssp, from GROMACS [23] based on the Define Secondary Structure of Proteins (DSSP) algorithm [33]. DSSP is the standard method of identifying secondary structures, and it is based on electrostatic energies among hydrogen donors and acceptors of the protein. Eight types of secondary structure are assigned, including hydrogen-bond based structures (helical, beta sheet, and turn) as well as non-hydrogenbond based structures (bend and coil). These structures are color coded and represented in a 3D map with the x-axis indicating the simulation time and the y-axis the residue location. Note that VMD was used to visualize the 3D image of the protein secondary structure. However, the protein structure calculation from VMD was based on the STRIDE algorithm [34] which is slightly different from DSSP. The assessment of membrane-induced protein folding, or disordered-to-order transition, of K18 was based on the DSSP only.

CG-Tau Oligomers in Solution
As described in Materials and Methods, CG WT-K18 and CG MBD-K18 oligomers were created by a self-aggregation process by bringing multiple (two or four) monomers near each other in a simulation box to create the initial structure. Oligomers were then created after performing CG MD simulations of the monomers in 0.1 M NaCl, at 310 K and 1 atmospheric pressure. Figure 1F demonstrates the highly dynamic nature of the aggregation process during the first 3 µs of the self-aggregation process starting from four monomers. Similar observation was found for the formation of the dimer as demonstrated in Supplementary Materials Figure SC1. When compared with tetramer, a dimer was readily formed in less than 1 µs ( Figure SC2A), but it took around 2 µs to establish an aggregated tetramer. Similar self-aggregation rates were found for the MBD-K18 in solution. After the protein self-aggregation process, the binding kinetics, domain preference and energetics of these pre-formed oligomers on different lipid rafts were systematically investigated.
Using a statistical and molecular interaction analysis tool, CONAN (see Materials and Methods), Supplementary Materials Figure SB1 shows the time-average and standard deviation protein-residue contact maps of monomeric and tetrameric WT-K18. It is clear that both monomer and tetramer exhibit large residue-contact regions as shown in the average contact maps and large fluctuation of residue-contacts as reflected in the standard deviation of contacts among the residues in the protein.

CG-Tau Oligomers Binding to Lipid Rafts
Three different lipid rafts, CO-raft ( Figure SC2A), GM-raft ( Figure SC2B) and PSraft ( Figure SC2C), in 0.1 M NaCl were used as model membranes. In the absence of an externally added protein, structurally heterogeneous, and highly dynamic lipid domains, DPPC-rich Lo (light green), DLPC-rich Ld (light orange), and mixed DPPC-DLPC Lod (dark green or dark orange) were observed. Figure SC2A demonstrates the domain structures of these Lo, Ld, Lod on both leaflets of the bilayer in both lateral and transverse views of the CO-raft. Figure SC2B also shows the GM1-clusters on one leaflet of the bilayer. Similarly, the PS-clusters which were found on Lo and Lod domains but not in the Ld domains on one side of the leaflet are demonstrated in Figure SC2C. The cholesterol was mainly found in the Lo domains in both lateral and transverse views for all rafts.
Using the g-select tool of GROMACS as mentioned in Methods, the lipid compositions of lipids in the Lo, Ld and Lod domains were determined, and the results are summarized in Figure SC5. For the cholesterol, about~60% of CHOL were distributed to the Lo domains vs.~10% and 30% to the Ld and Lod domains, respectively, for the CO-raft and GMraft. For the PS-raft, the distributions of CHOL were~50%, 15% and 35% in the Lo, Ld and Lod, respectively. Therefore, the presence of the GM-clusters had no effect on the CHOL distributions in Lo, Ld and Lod domains. However, the presence of PS-clusters slightly reduced the CHOL partitioning to the Lo domains and concomitantly increased its partitioning to the Ld and Lod domains. The DPPC % in the Lo domain and DLPC % in the Ld domain reflect the relative sizes of the Lo and Ld domains in each raft. Here, Lo-DPPC% were~20% in all rafts, and Ld-DLPC~10%. Furthermore, the total Lod % which is the sum of Ld-DPPC% and Lo-DLPC% accounted for~70%. Therefore, our raft systems provided the dynamically stable Lo, Ld and Lod domains, as well as the asymmetrically distributed GM-clusters and PS-clusters, to examine the lipid domain-binding preference of tau oligomers on raft surfaces. Protein-residue contact maps ( Figure SB2) of monomeric and tetrameric WT-K18 bound to CO-raft revealed smaller residue-contact regions and fluctuation of residue-contacts were evident when compared to the protein in solution.
To investigate the binding behaviors of tau oligomers on raft surfaces, the pre-formed, self-aggregated oligomers were placed at different positions (see Materials and Methods) in solution to create the starting or 0 µs structures of the CG oligomer and raft simulation replicates. Figure 2 demonstrates the binding behavior of tetramers starting from the solution states and above the surfaces of the CO-raft (Figure 2A), GM-raft ( Figure 2B) and PS-raft ( Figure 2C) to the final equilibrated, membrane-attached or -bound states after 15 µs as shown in Figure 2D-F, accordingly. Figures SC2 and SC3 demonstrate the similar tau dimer binding behaviors to the CO-raft, GM-raft and PS-raft from the solution states to the membrane-bound states.
DPPC% were ~20% in all rafts, and Ld-DLPC ~10%. Furthermore, the total Lod % which is the sum of Ld-DPPC% and Lo-DLPC% accounted for ~70%. Therefore, our raft systems provided the dynamically stable Lo, Ld and Lod domains, as well as the asymmetrically distributed GM-clusters and PS-clusters, to examine the lipid domain-binding preference of tau oligomers on raft surfaces. Protein-residue contact maps ( Figure SB2) of monomeric and tetrameric WT-K18 bound to CO-raft revealed smaller residue-contact regions and fluctuation of residue-contacts were evident when compared to the protein in solution.
To investigate the binding behaviors of tau oligomers on raft surfaces, the preformed, self-aggregated oligomers were placed at different positions (see Materials and Methods) in solution to create the starting or 0 μs structures of the CG oligomer and raft simulation replicates. Figure 2 demonstrates the binding behavior of tetramers starting from the solution states and above the surfaces of the CO-raft (Figure 2A), GM-raft ( Figure  2B) and PS-raft ( Figure 2C) to the final equilibrated, membrane-attached or -bound states after 15 μs as shown in Figure 2D-F, accordingly. Figures SC2 and SC3 demonstrate the similar tau dimer binding behaviors to the CO-raft, GM-raft and PS-raft from the solution states to the membrane-bound states.   The tau oligomer, both WT-K18 and MBD-K18, bound to the Lod domain in the COraft, and the asymmetrically distributed GM-cluster or PS-cluster on one leaflet of the bilayer. For the PS-raft binding, all three independent replicates of tetrameric WT-K18 and MBD-K18 bound exclusively to the PS-clusters within 6 µs. However, one replicate of monomeric WT-K18, one replicate of dimeric WT-K18 and one replicate of monomeric MBD-K18 bound to the leaflet in the absence of PS. In other words, the protein bound to the opposite leaflet, where there were no PS lipids present. In those systems, an extra replicate, or fourth replicate, with a slightly different initial position from the raft surface, was created, and successful PS-cluster binding of the fourth replicate was achieved. This observation suggests that the binding affinity to PS-clusters and Lod domains may be similar for the small tau oligomers. In contrast, all oligomers bound to GM-clusters and no non-specific binding to the "opposite leaflet" were observed in both WT-K18 and MBD-K18. The results further indicate that binding affinity to GM-clusters may be stronger than the binding affinity to Lod and PS-cluster domains, for tau oligomers of all sizes.
No significant disruptions of the domain organization were detected in all our oligomer binding studies. Figure SC5 shows the time-and replicate-averaged domain compositions over the last 5 µs and across all simulation replicates in the presence of membrane-bound monomeric (n = 1), dimeric (n = 2) and tetrameric (n = 4) of CG WT-K18 and CG MBD-K18. No significant changes of the domain compositions were evident before and after the protein attachment for both CG WT-K18 and CG MBD-K18. These results indicate that the lateral organization of the raft domains was unaffected by membrane-bound oligomers.
At the end of the CG simulations, a CG-to-AA resolution transformation step was carried out to convert the final CG structures to the starting AA structures, as demonstrated by overlaying the CG and AA structures in Figure 2D-F. Each converted AA structure then underwent 200 ns of AA simulation, and the final 200 ns AA structures of the oligomer/raft complexes are shown in Figure 2G-I. All tau oligomers remained firmly attached to the raft surfaces during the AA simulations.

Lipid Binding Analysis of K18 Oligomers
A minimum distance (mindist) analysis (see Materials and Methods) based on a 3-panel mindist plot format was performed to investigate the binding kinetics and residue-specific binding domain of K18 oligomers on all three rafts. Figure 3 shows the mindist plots for the representative replicates of CG monomers and dimers of WT-K18 and MBD-K18 on all lipid rafts, and Figure SD1 in Supplementary Materials shows the plots for the representative replicates of CG tetramers and AA tetramers of WT-K18 and WT-K18 for all lipid rafts.
To study the binding kinetics, the mindist vs. time (upper panel of the 3-panel plot) was used to quantify the lipid-binding time of CG K18 oligomers to raft surfaces. Here, the lipid binding time is defined as the time at which the protein undergoes a transition from the solution state with a large mindist value to a stable, membrane-bound state with a small and stabilized mindist value. For example, abrupt declines or transitions in mindist were detected at 1.72 µs and 7.63 µs for the WT-K18 monomer and tetramer in the CO-raft, respectively, as demonstrated in Figure 3A,D. Table 1 summarizes the lipid-binding time for all 54 independently simulated CG replicates. In general, for the CO-raft, the binding times varied over a wide range from 0.1 to 20 µs. All K18 monomers in the CO-raft established stable membrane binding within 10 µs; however, the K18 dimers, particularly the WT-K18 dimers, experienced longer binding times, ranging from 14 to 20 µs. On the other hand, the binding times of K18 to PS-raft were within 6 µs, and those to GM-raft were even shorter, i.e., within 3 µs. These results suggest stronger or more efficient binding of tau oligomers to the anionic GM-and PS-rafts than to the CO-rafts.   The number of contacts vs. time plots (mid-panel of the 3-panel plot) was used to investigate the lipid environment surrounding the membrane-bound oligomers. As expected, the number of contacts between protein and lipid atoms increased with the size of the oligomer, i.e., from monomer to tetramer, for each lipid raft, as demonstrated in Figures 3 and 4. Additionally, the number of contacts between the protein and GM1 lipids in the GM-raft ( Figure 3B,E) was about three times more than those between the protein and other non-GM1-lipids in the same GM-raft, as well as in other rafts, i.e., CO-raft ( Figure 3A,D) and PS-raft ( Figure 3C,F). The observation suggests that GM1-lipids provide a more stable binding environment and a larger number of protein binding sites to K18 than other non-GM1 lipids in all three lipid rafts.
Other than the number of contacts data, Table 1 shows the time-and replicate-averaged lipid compositions and the total number of lipids in the 0.5 nm annular lipid (AL) shell that provide a quantitative evaluation of the lipid environment of the membrane bound-protein.
Overall, the total number of lipids within the 0.5 nm AL shell increased steadily as the size of the oligomer increased. Additionally, the oligomers of all sizes bound to the Lod domain in CO-raft, GM1-clusters in the GM-raft, and PS-clusters in the PS-raft. For the CG-system, the number of lipids in the AL shells for the WT-K18 oligomers was greater than those for the MBD-K18 oligomers, suggesting a stronger affinity of WT-K18 to GM1 when compared to MBD-K18. Yet, no significance difference was found in the AA-system.  To investigate the residue-resolved lipid-binding conformation of oligomers, the timeaveraged mindist over the last 5 µs of simulation vs. protein residue number, or mindist spectrum, of each simulation system was examined (SI-D). Due to the presence of multiple chains, accumulated residue numbers, i.e., 1-130, 131-260, 260-390 and 391-520, are used to label the residue locations A-D, accordingly. The locations of mutations are marked with dashed red lines in the plots, specifically at 45*, 66* and 76* for chain A, 175*, 196* and 206* for chain B, 305*, 362* and 336* for chain C, and 435*, 492* and 466* for chain D. Note that any dips or minima in the mindist vs. residue number plots reflect strong affinity for the protein at those residues where the dips are found. For reference, only dips with mindist less than 1.0 nm were considered.
A difference spectrum (WT-K18-MBD-K18) represents an intuitive way to look at how WT-K18 and MBD-K18 differ in preference of residue-lipid contacts. It is more helpful to consider replicate-averaged mindist than chain-averaged mindist, because the conformational geometry of larger oligomers (i.e., tetramers) prevents all chains from equally accessing the membrane. Even without chain-averaged results, the mindist spectra of higher order oligomers are too convoluted to track meaningful protein-residue binding preferences which are otherwise obvious in monomers and dimers.
Our results reveal that the point mutations significantly decrease the K18's binding affinity and stability in the CO-raft across all sizes of oligomers. For example, the difference spectrum plot of the K18 monomer in the CO-raft as shown in Supplementary Materials Figure SD2 demonstrates a clear dip between the second mutation at residue 66* and the third mutation at residue 76*. The mindist plots of K-18 monomers ( Figure 3A,D) reveal that the WT-K18 binds closely to the membrane (about 0.5 nm) at the second point mutation, residue 66*, whereas an opposite effect is evident with the MBD-K18. The only region where the MBD-K18 oligomers exhibit clean and consistent binding is at the N-terminus of each chain, a prominent binding site across all oligomer types and sizes.
For the GM-raft, most residues in both WT-K18 and MBD-K18 interacted strongly with the GM1 lipids. Figure 3 and Supplementary Materials Figure SD3 show that the minimum distance between the WT-K18 and the GM lipid maintains a relatively flat baseline of about 0.5 nm, with small peaks around residues 73 and 90. Conversely, the MBD-K18 monomer exhibits a less stable baseline, also at about 0.5 nm, but with several small peaks from residues 20-30, 66*, 73, and 90-100. Overall, the point mutations appear to have little effect on the protein binding pattern in the presence of GM1 lipids; however, the distances between protein and membrane components do indicate that the WT-K18 binds more closely to GM1-clusters than does the MBD-K18.
For the PS-raft, both the replicate-averaged WT-K18 monomer and dimer appear to bind more closely with the PS membrane than the MBD-K18, as shown in Figure 3. The mindist difference plot of the monomer shows a clear dip at residue 66*, as does the chain B of the dimer at residue 196* (Supplementary Materials Figure SD4). Notably, the minimum distances between the protein and water molecules also exhibit significant peaks at these points, up to about +0.15 nm. These results suggest that (a) the non-mutated, hydrophobic residues 66* and 196* represent the major binding sites of WT-K18, and (b) the mutated, negatively charged residues of MBD-K18 bind poorly to the anionic PS-lipid domains relative to the other rafts.
Due to the limited phase sampling behaviors of AA simulations, extensive differences in the protein-lipid binding patterns in AA vs. CG are not expected. However, in this work we observed subtle differences in protein-binding patterns in asymmetric rafts. In general, we observed a stronger protein-lipid interactions in AA simulations vs. CG simulations at various regions, e.g., residue 130 of MBD-K18 on GM-raft (Supplementary Materials Figure SD6), residues 76* and 435* of WT-K18 on PS-raft, and residues 42, 492* of MBD-K18 on PS-raft (Supplementary Materials Figure SD7). Certainly, much longer AA simulation times are needed to explore the detailed protein-lipid interactions of K18 on raft surfaces.

Protein-Lipid Interactions of Tau Oligomers
Protein-lipid binding energies were calculated to evaluate the binding affinity of tau oligomers with each lipid type, and the results are shown in Figure 4 for oligomers of different sizes at CG and AA resolutions. All binding energies represent both the timeaverage and replicate-average over the last 5 µs for CG or 50 ns for AA, and across all three independent replicates. Here, nonbonded Lennard-Jones and Coulomb (electrostatic) binding energies and their sum (total binding energy) were determined. For the CGoligomers of all sizes, the Coulomb energy was much less than 10% of the total protein-lipid binding energy. In contrast, for the AA-oligomers, the Coulomb energy was made up of >50% and~20-30% of the total protein-lipid binding energy for diacyl lipids (Figure 4Q-S) and cholesterol ( Figure 4P), respectively. This difference is attributed to the difference in the fine and extensive partial charge distribution of the atomistic structure of the protein and lipid molecules in the AA force fields vs. the CG force fields.
The protein-cholesterol binding energy of WT-K18 was generally stronger than that of MBD-K18 for CG-oligomers of all sizes in all three lipid rafts ( Figure 4A,F,K), with the largest energy being present in the PS-raft. No significant difference was found between WT-K18 and MBD-K18 for AA-tetramers ( Figure 4P).
The protein-DPPC binding energy of WT-K18 was also stronger than that of MBD-K18 for the CG-oligomers of all sizes in all rafts ( Figure 4B,G,L). For the AA-tetramer, the protein-DPPC binding energy of MBD-K18 was slightly higher than that of WT-K18 in CO-raft and GM-raft, but the binding energy of WT-K18 was almost twice as large as that of MBD-K18 in PS-raft ( Figure 4Q), indicating a significantly stronger binding affinity of WT-K18 to DPPC when compared with MBD-K18 at both CG and AA resolutions. In addition, in the AA resolution, the Coulomb energy between protein and DPPC for the WT-K18 was twice as large as that for the MBD-K18, indicating a strong electrostatic interaction between DPPC and WT-K18 protein in the PS-raft.
For DLPC, the protein-lipid binding energy was overall smaller than that for other diacyl lipids for both CG-and AA-oligomers. Interestingly, the protein-DLPC binding energy of MBD-K18 was twice as large as that of WT-K18 for the AA-oligomers in GM-raft as shown in Figure 4R.
For both the GM1 and POPS lipids, WT-K18 consistently exhibited greater proteinlipid binding energies than MBD-K18 across all dimers and tetramers at both CG and AA resolutions, indicating that the WT-K18 oligomers were more strongly attracted to the anionic GM1 and PS domains than other lipid types. In addition, the binding energy between the protein and GM1 was almost 6-10 times stronger than the that between the protein and non-GM1 lipids, e.g., the protein-GM1 binding energy of~−8000 kJ/mol vs. the protein-DPPC binding energy of~−1200 kJ/mol for the AA-tetramer on the same GM-raft surface. Similarly, the protein-POPS binding energy (~−3000 kJ/mol) was about 2 times stronger than that of the protein-DPPC (~−1200 kJ/mol) for the AA-tetramer on the same PS-raft surface.
Since the protein-lipid binding energy depends strongly on the number of lipids surrounding the protein, we have also calculated the normalized protein-lipid binding energy by dividing the total binding energies by the number of annular lipids with a cutoff-distance of 1.2 nm (Supplementary Materials Figure SE1). The choice of 1.2 nm reflects the cut-off energy calculation threshold, as described in previous studies [16,17]. The resulting normalized protein-lipid binding energy plots are shown in Supplementary Materials Figure SE2. The normalized protein-lipid binding energies for CHOL, DPPC and DLPC were within~10 to~40 kJ/mol, and no significant differences were found between protein-lipid binding energies of WT-K18 and MBD-K18 at both CG and AA resolutions, with the exception that the protein-DPPC binding of WT-K18 tetramer was much greater than that of the MBD-K18 tetramer in the PS-raft at the AA resolution (Supplementary Materials Figure SE2N), similar to the trend in total binding energy ( Figure 4Q). Our normalized protein-lipid binding results further reveal the ranking order of the lipid binding affinity of tau-oligomers according to the lipid type. The order is GM1 >> PS > DPPC > DLPC > CHOL.

Structure of Annular Lipids Surrounding Tau Oligomers
The time-and replicate-averaged lipid orientational order parameters vs. carbon number of 0.5 nm AL lipids surrounding the AA WT-K18 and AA MBD-K18 tetramers over the last 50 ns and across three independent simulation replicates were calculated ( Figure 5). Here, the 0.5 nm-AL lipids are selected according to the minimum distances of the atoms of the proteins and lipids that were within 0.5 nm (see Materials and Methods). Figure 5A,B illustrate the selections of 0.5 nm AL lipids of different types surrounding the representative replicates of membrane-bound AA WT-K18 and AA MBD-K18 on CO-raft. Note that only the lipids in one leaflet of the lipid bilayer were selected in both cases. In general, the acyl chain order parameter of the acyl chains of any diacyl chain lipid, DPPC, DLPC, GM1, or POPS, increased progressively, reached a peak between carbon numbers 4 and 8, and declined steadily as the lipid carbon number increased. Note that the chain carbon number starts from 1 (near the polar lipid headgroup) and ends at the terminal methyl carbon at the end of the chain. To improve the statistics of order parameter calculations, the order parameters from both sn-1 and sn-2 chains of the diacyl chain lipids were combined.
For the CO-raft, the acyl chain orders of DPPC in the AL lipids surrounding the AA WT-K18 tetramer were significantly lower than those surrounding the MBD-K18 tetramer for carbon numbers 1-10 ( Figure 5C). Yet, no significant difference in the order parameters of the AL lipids for WT-K18 vs. MBD-K18 was observed for DLPC across all carbon numbers ( Figure 5D). As controls, the acyl chain orders of either DPPC or DLPC in the non-AL lipids, nAL, or lipids outside the AL lipid shells, for the AA WT-K18 were identical to that for the AA MBD-K18. The results suggest that the AA WT-K18 perturbed the DPPC and the DLPC chains' order more than those of the AA MBD-K18 at the top half of the acyl chain region for the CO-raft system. The results indicate that WT-K18, compared to MBD-K18, significantly perturbed the saturated lipids, DPPC, upon binding to the CO-raft.
For the GM-raft, the DPPC acyl chain order of the AL lipids surrounding the AA WT-K18 was slightly lower than those surrounding the AA MBD-K18 for carbon numbers 1-6 ( Figure 5E). However, for DLPC, a significantly lower order parameter of the AL lipids surrounding AA WT-K18 than those surrounding AA-MBD-K18 at low carbon numbers 0-4 was detected. However, the trend was reversed at high carbon numbers 14-16 near the end of the acyl chain. The GM1 acyl chain order of the AL lipids surrounding the AA MBD-K18 was slightly lower than those surrounding the AA WT-K18 at large carbon numbers 14-16. For the nAL lipids, the order parameters of DPPC and DLPC were identical for both AA WT-K18 or AA MBD-K18. However, in the presence of AA WT-K18, the GM1 chain orders of the nAL lipids were very similar to those of the AL lipids. In contrast, in the presence of AA MBD-K18, the GM1 chain orders in the nAL lipids were higher than those of the AL and nAL lipids in the presence of AA WT-K18. The above observations indicate that AA-WT-K18 disrupted the DPPC and DLPC chain orders of the AL lipids more than AA-MBD-K18. The GM1 chain order disruptions induced by AA WT-K18 and AA MBD-K18 are similar in the AL lipids. However, the disruptive effect by AA WT-K18 extended to the nAL lipids. In other words, the AA WT-K18 disrupted the GM1 chain orders to the entire GM1 lipid population (AL and nAL), AA MBD-K18 exerted a localized effect at the AL lipid region only. The results indicate that WT-K18 on the GM-raft significantly disordered the acyl chains of nearest DLPC lipids and all the GM-lipids upon binding to the GM1-cluster.
For the PS-raft, no significant difference in the DPPC order parameters of the AL or nAL lipids surrounding the AA WT-K18 and those surrounding the AA MBD-K18 was found. However, the DLPC order parameters of the AL lipids surrounding the AA WT-K18 were slightly lower than those surrounding the AA MBD-K18 at carbon numbers 10-12. The trend was reversed for the POPS order parameters. Again, for the nAL lipids, the order parameters of all nAL lipids surrounding the AA WT-K18 were identical to those surrounding the AA MBD-K18. The lack of significant perturbation of lipid chain order in POPS by both WT-K18 and MBD-K18 suggests that the interactions of protein with the lipids in PS-raft mainly affect the headgroup region with minimal effects at the acyl chains below the headgroups of the lipids. For the GM-raft, the DPPC acyl chain order of the AL lipids surrounding the AA WT-K18 was slightly lower than those surrounding the AA MBD-K18 for carbon numbers 1-6 ( Figure 5E). However, for DLPC, a significantly lower order parameter of the AL lipids surrounding AA WT-K18 than those surrounding AA-MBD-K18 at low carbon numbers 0-4 was detected. However, the trend was reversed at high carbon numbers 14-16 near

Protein Folding on Raft Surfaces
The protein secondary structures of membrane-bound AA WT-K18 and AA MBD-K18 tetramers on all three different raft surfaces have been calculated using the DSSP (see Materials and Methods) as shown in Supplementary Materials. Within the 200 ns all-atom MD simulations, evidence of transitions from non-hydrogen bonded, disordered structures-such as turn or coil-to hydrogen-bonded, ordered structures-such as bend, helix, or beta-were found for certain regions of the protein. Most of the appearances of the ordered structures, particularly beta-sheet or beta-bridge, are transient. Interestingly, certain replicates, e.g., one replicate of WT-K18 and one replicate of MBD-K18, showed significant beta-sheet formation as shown in Figure 6. In these two replicates, three out of four chains participated in beta-sheet formation. More beta-sheet structures (red strips in the DSSP plot) were detected in WT-K18 than in MBD-K18. The structures and locations of these beta-sheets at the end of the 200 ns-simulation were visualized by VMD (see Materials and Methods), which uses a slightly different secondary structure calculation algorithm than DSSP. Here, 10 beta-sheets in WT-K18 but only 6 beta-sheets in MBD-K18 were detected at 200 ns by VMD. Both inter-chain and intra-chain beta-sheets were found in both WT-K18 and MBD-K18. Therefore, our results indicated that some of our K18 oligomers underwent raft surface-induced folding, i.e., from a disordered state in solution to a partially ordered state in the membrane-bound state, for some of our oligomers within 200 ns of AA simulations in this pilot study. Longer simulation times will be needed to detect more stable beta-sheet structures on the raft surfaces.

Discussion
Binding kinetics and energetics of Tau-K18 oligomers on various membrane domain models were systematically investigated. Among all three co-existing domain types, Lo, Ld, and Lod, all Tau-K18 oligomers bound to the Lod as demonstrated by the mindist vs. time and composition of the annular lipids surrounding the membrane-bound proteins. This Lod domain preference of Tau-K18 binding agrees well with our previous studies on

Discussion
Binding kinetics and energetics of Tau-K18 oligomers on various membrane domain models were systematically investigated. Among all three co-existing domain types, Lo, Ld, and Lod, all Tau-K18 oligomers bound to the Lod as demonstrated by the mindist vs. time and composition of the annular lipids surrounding the membrane-bound proteins. This Lod domain preference of Tau-K18 binding agrees well with our previous studies on highly ordered beta-amyloid fibrils Aβ   [16] and the highly disordered full length Aβ 1-42 oligomers [17]. Our results, therefore, suggest that the Lod domains are common membrane nanodomain regions for both disordered amyloidogenic proteins, such as betaamyloid and tau, to bind to. Interestingly, apart from amyloidogenic proteins, another class of membrane-active protein, e.g., HIV-gp42 fusion peptide, has been shown to bind to these Lod domains. The line tension created at the domain boundaries due to the thickness mismatch between the thicker Lo domain and the thinner Ld domain has been proposed to be the major molecular mechanism driving membrane active protein to Lod domain [35,36]. More computational and experimental studies are needed to explore the biophysical mechanisms of protein binding to those self-assembling and dynamic Lod membrane domains.
For the asymmetrically distributed GM-raft containing co-existing GM1-clusters within the Lo, Ld, and Lod domains, or 4 lipid domains on one leaflet, and 3 lipid domains (Lo, Ld, and Lod) on the "opposite" leaflet, tau proteins of all sizes bind exclusively to the GM1-clusters. This exclusive binding property of all our model tau proteins to the GM1clusters indicates that disordered tau proteins of different sizes have the strong propensity to interact with the GM1-lipids, which are found exclusively on the outer surface of the neuronal membranes [37]. It is well established that initial misfolded and self-aggregated tau proteins are in the cytosol of the neurons [4,6] of the brain. Here, our observations of disordered tau proteins interacting strongly with GM1 lipids, a major biomarker [38], of the outside leaflet of the plasma membranes of neurons, lead us to conclude that GM1 could be a major target of disordered tau oligomers once they are able to migrate and escape from the cytosol of the neurons. Note that a recent work suggests that tau proteins do exit the neurons and translocate to other areas of the brain, including the space between the synaptically associated neurons [5]. In addition, our previous studies on ordered Aβ 17-42 fibrils and disordered Aβ 1-42 oligomers also revealed strong binding of amyloid protein aggregates to GM1-clusters [16,17]. It is believed that the presence of highly hydrophilic and disordered carbohydrate groups of the GM1 can provide strong binding sites to the hydrophilic residues of the disordered tau protein. The detailed molecular mechanism of carbohydrate-amino acid residue interaction that stabilizes the tight-binding of GM1 to disordered tau protein remains to be explored computationally and experimentally.
For the asymmetrically distributed PS-raft containing co-existing PS-clusters within the Lo, Ld and Lod domains, or again 4 lipid domains on one leaflet, and 3 lipid domains (Lo, Ld, and Lod) on the "opposite leaflet", tau tetramer bound exclusively to the PS-clusters on one leaflet. By reducing the charge of Tau-K18 from +10e to +7e, and further decreasing the hydrophobicity of tau, i.e., the mutation from WT-K18 to MBD-K18, a significant decrease in the protein-POPS binding energy as well as a significant difference in the protein-residue minimum distance (mindist) at the mutation region I308E were evident. The results suggest that electrostatic interaction between the cationic K18 and negatively charged PS is a major membrane binding mechanism of Tau-K18 with PS-containing membrane domains, exclusively found in the inner leaflet of neuronal membranes. Our observations agree with previous experimental work on the preferred binding of WT-K18 over MBD-K18 on PS containing liposomes, and hence validate our computational approach of investigating disordered tau oligomer interactions with raft membranes.
Other than identifying the membrane domain targets of Tau-K18, we have also attempted to quantitate the extent of membrane structural disruption, or membrane damage, and the associated membrane-induced beta-sheet formation in the presence of the membrane-bound Tau-K18 oligomers. Membrane damage and beta-sheet formation are considered to be the major factors associated with the interaction of tau aggregates with lipid membranes and contribute significantly to the early molecular events of tauopathies that initiate the onset of Alzheimer's [1,4,6]. Substantial disruptions of lipid structural packing, as quantified by a significant decrease in the orientational order parameter of annular lipids surrounding the membrane-bound tau oligomers, were evident. Interestingly, WT-K18 disrupted the annular lipids more than MBD-K18, particularly the saturated DPPC in all lipid rafts. The observation of a stronger order disruption of unsaturated DLPC by WT-K18 than by MBD-K18 in GM-raft suggests that WT-K18 was able to perturb the unsaturated lipid more effectively than MBD-K18. Since the orientational order of lipid acyl chains are important to maintain the physiological functions of membranes, a large disruption in the lipid order will create detrimental effects on the structure-function relationship of membranes, inducing tauopathies associated with membranes [2].
Although our tau protein is a relatively large amyloid aggregate when compared to, say, Aβ 1-42 (with tetramers of 520 residues vs. 168 residues, respectively), significant but localized beta-sheet formations in the membrane-bound tau aggregates were evident in some of our simulation replicates. For the case of the tau tetramer on the CO-raft, small interchain and intrachain beta-sheet formations were found, and WT-K18 appeared to produce and sustain more beta-sheets than MBD-K18. Certainly, much longer simulation time and advanced sampling techniques are required to further explore the role of membrane domains on membrane-bound beta-sheet formations. Our preliminary results in this study indicate that membrane-assisted beta-sheet formation, starting from a disordered solution state to the membrane-bound state, can be detected for relatively large amyloid protein aggregates like Tau-K18. Experimental observations of beta-sheet formation on the lipid membrane surface, particularly anionic membrane surface, of Tau-K18 has recently been reported [12]. Our atomistic simulation work provides computational evidence that interchain and intrachain beta-sheet formations are not only possible, but likely to occur on anionic membrane surfaces.
The formation of small beta-sheets on membrane-bound amyloid aggregates may represent the early events of creating more membrane-disruptive structures associated with the pore, carpeting and detergent-like models [2]. The higher rate of beta sheet formation in WT-K18 over MBD-K18 suggests that electrostatic interactions play a key role in membrane-templated protein folding [12]. Furthermore, these small beta-sheet structures may act as seeds to attract or recruit more amyloid proteins from solution or adjacent membrane-bound amyloid proteins to create even more membrane-disruptive aggregates that further damage the membranes [4]. Therefore, the observations of small beta-sheets in this pilot study provide the base for future experimental or computational work investigating the nucleation processes of tau on membrane domains [4,6].
In this study, we have used multiscale molecular dynamics simulations to investigate early tau binding events to three structurally distinctive model raft membranes on the microsecond time scale and with both coarse-grained (CG) and atomistic (AA) spatial resolutions. CG simulations allow us to explore protein-lipid binding events up to tens of microseconds, while subsequent AA simulations provide protein-folding and detailed nonbonded interactions, such as electrostatic and van der Waals, that CG simulations do not. Due to the large number of atoms, partial charge distributions of molecules, and a lack of phase sampling capability of conventional AA simulations, limited proteinlipid conformational sampling around the pre-established protein-lipid contact regions can only be explored within a few hundred nanoseconds. Our multiscale simulations represent a useful approach to investigate protein-lipid binding events and protein-folding on surfaces upon binding associated with tau-membrane interactions. Note that this multiscale simulation approach has successfully been applied to other systems, such polymer-ion complex and protein in solution [39][40][41].

Conclusions
Using multiscale molecular dynamic simulations (CG-MD and AA-MD), the binding behaviors of WT-K18 and MBD-K18 oligomers onto phase-separated lipid domains, with and without anionic biomarkers of the inner and outer leaflet of the neuron, were explored. By using models mimicking different domains of realistic neuronal membranes and multiscale simulation, we were able to holistically investigate binding events of disordered tau aggregates on lipid rafts, as well as high-resolution protein folding and protein-lipid conformation. Electrostatic interactions between Tau-K18 and lipid domains indicate the preferential binding regions of tau oligomers, prioritizing in the order of GM1 clusters (outer membrane leaflet's biomarker), then PS cluster (inner membrane leaflet's biomarker), and finally, the mixed Lod domain. Additionally, detailed membrane-order analyses in atomistic resolution suggest significant membrane damage caused by both WT-K18 and MBD-K18 to CO-raft and GM-raft, but not PS-raft, suggesting the membrane disruptive of tau oligomers, especially WT-K18 oligomers, once they delocalize from the inside of a neuron. Beta-sheet motifs in membrane templated-protein folding were also evident, where WT-K18 had higher propensity for beta sheet formations than MBD-K18, suggesting the toxicity of disordered tau oligomers and potential seeding effect for further membrane damage. Overall, results indicated that the three point mutations (V287E, I308E, V318E) in MBD-K18 successfully decreased membrane binding and disruption on neuronal surfaces, but still significantly affected membrane lipid orientations, especially when relocated to outer neuronal leaflets containing GM1 clusters. The technique of multiscale MD simulation and analysis provides a new approach to understand the insights of tauopathic membrane-associated mechanisms and damage.