University of Southern California Dissertations and Theses

Molecular dynamics studies of protein aggregation in unbounded and confined media

(USC Thesis Other)

Molecular dynamics studies of protein aggregation in unbounded and confined media

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Molecular Dynamics Studies of Protein Aggregation in Unbounded and Confined Media by Size Zheng A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirement for the Degree DOCTOR OF PHILOSOPHY (CHEMICAL ENGINEERING) December 2018 Copyright 2018 Size Zheng i DEDICATION To my son, Leon Zheng, and my wife, Nian Liu. ii ACKNOWLEDGMENT Here I would like to thank all the people who made this dissertation possible. First and foremost, I would like to thank my two advisors, Professor Muhammad Sahimi and Professor Kathrine Shing. It is their great mentorships that support me to go through the whole PhD journey. When I got trapped, they always pull me up, encouraging me to think differently and guiding me back to the right path; when we had disagreements, they always listen patiently, providing me with broad resources and knowledge and helping me to improve; when we obtained achievements, they always carefully push me further, inspiring me to read more and learn more and reminding me we could keep doing better. I could never go thus far without their guidance. Meanwhile, I would like to thank my qualifying committee members: Professor Aiichiro Nakano for his great advice and help on the optimization of code and the improvement of performance; Professor Rajiv Kalia for his positive feedback and inspiration on the program parallelization; Professor Ted Lee for his probing questions and the contribution of the experimental confirmation of our proposed hypothesis. I also would like to express a special gratitude to Professor Priya Vashishta. It was on his class that I was first time introduced molecular dynamics and high-performance computing. His teaching inspired me a great interest on molecular dynamics researches. My sincere thanks also go to Professor Tao Wei for his help on my initial learning of molecular simulation back to master degree and his encouragement and support on AIChE conferences; Dr. Leili Javidpour for her patiently replying my emails full of questions on the program framework; my officemates, Congyue Wang and Meng Wang, for their help and support on the projects and daily life; the staff, Tina Silva and Shokry Bastorous, for their care and help on my teaching assistant job in the lab; Professor Hung Nguyen and Professor Feng Ding, for offering help on DMD development. Last but not least, I would like to give my deepest gratitude to my family. It is their love and support that encourage me to keep going forward and overcome all the obstacles. Thank you very much! iii TABLE OF CONTENTS DEDICATION ................................................................................................................................................. i ACKNOWLEDGMENT ....................................................................................................................................ii TABLE OF CONTENTS................................................................................................................................... iii LIST OF TABLES ............................................................................................................................................ vi LIST OF FIGURES ......................................................................................................................................... vii ABSTRACT ..................................................................................................................................................... x 1. INTRODUCTION .................................................................................................................................... 1 1.1 Overview of protein aggregation ...................................................................................................... 1 1.2 Aggregate structure and aggregation mechanism ............................................................................ 2 1.3 Molecular simulation and protein model .......................................................................................... 4 1.4 Summary of molecular simulation studies on protein aggregation .................................................. 6 1.5 Protein aggregation in cellular environment .................................................................................... 9 2. DISCONTINUOUS MOLECULAR DYNAMICS AND MOLECULAR MODEL ............................................. 11 2.1 Discontinuous molecular dynamics ................................................................................................. 11 2.2 Molecular model .............................................................................................................................. 14 2.2.1 Intermediate-resolution coarse-grained model ....................................................................... 14 2.2.2 High-resolution all-atom model ................................................................................................ 19 2.2.3 Confined environment model ................................................................................................... 22 3. DYNAMICS OF PROTEIN AGGREGATION. I. UNIVERSAL SCALING IN UNBOUNDED MEDIA .............. 24 3.1 Introduction ..................................................................................................................................... 24 3.2 Configuration of DMD simulation.................................................................................................... 24 3.2.1 Models of peptides.................................................................................................................... 24 iv 3.2.2 Molecular simulation ................................................................................................................ 25 3.3 Results .............................................................................................................................................. 26 3.3.1 Formation of β-sheets and fibrils .............................................................................................. 27 3.3.2 Number and size of the aggregates .......................................................................................... 32 3.3.3 Number of β-sheet forming peptides ....................................................................................... 40 3.3.4 Number of hydrogen bonds ...................................................................................................... 42 3.3.5 Radius of gyration of the aggregates ........................................................................................ 44 3.3.6 Diffusivity ................................................................................................................................... 45 3.4 Discussion: the catalytic effects ...................................................................................................... 46 3.5 Summary .......................................................................................................................................... 47 4. DYNAMICS OF PROTEIN AGGREGATION. II. DYNAMIC SCALING IN CONFINED MEDIA ..................... 49 4.1 Introduction ..................................................................................................................................... 49 4.2 Configuration of DMD simulation.................................................................................................... 50 4.2.1 Models of peptides.................................................................................................................... 50 4.2.2 Models of confined media ........................................................................................................ 51 4.2.3 Molecular simulation ................................................................................................................ 51 4.3 Results .............................................................................................................................................. 53 4.3.1 Dependence of the protein secondary structure on the geometry of confined media ........... 54 4.3.2 Dynamic scaling of the aggregation process ............................................................................. 57 4.3.3 Number of peptides that form β-sheets ................................................................................... 63 4.3.4 Distribution of various types of hydrogen bonds ..................................................................... 66 4.3.5 Effect of the peptide’s size ........................................................................................................ 68 4.3.6 Effect of the confined medium’s size ........................................................................................ 70 4.4 Dynamic scaling and fractal protein aggregates ............................................................................. 73 v 4.5 Differences between aggregations in confined and unbounded media ......................................... 75 4.6 Summary .......................................................................................................................................... 76 5. MOLECULAR DYNAMICS STUDY OF THE STRUCTURE OF POLY-GLYCINE-ALANINE IN AGGREGATION AND FOLDING............................................................................................................................................ 78 5.1 Introduction ..................................................................................................................................... 78 5.2 Configuration of DMD simulation.................................................................................................... 79 5.2.1 Models of proteins .................................................................................................................... 79 5.2.2 Molecular simulation ................................................................................................................ 80 5.3 Results .............................................................................................................................................. 81 5.3.1 System temperature ................................................................................................................. 81 5.3.2 Formation of β-sheets ............................................................................................................... 82 5.3.3 Unique helical structures .......................................................................................................... 85 5.4 Discussion and summary ................................................................................................................. 91 6. APPENDIX ........................................................................................................................................... 93 REFERENCE ................................................................................................................................................ 99 vi LIST OF TABLES Table 2-1 Parameter values of PRIME model ........................................................................................... 15 Table 2-2 Parameter values for hydrogen bond ....................................................................................... 18 Table 4-1 Summary of the exponents of dynamic scaling ........................................................................ 60 vii LIST OF FIGURES Figure 1-1 Cartoon representation of the quaternary structure of an Aβ 1-40 fibril, viewed parallel and perpendicular to z axis. ............................................................................................................................... 3 Figure 2-1 Potential profiles in DMD ........................................................................................................ 12 Figure 2-2 Flow chart of simulation process ............................................................................................. 13 Figure 2-3 Schematic diagram of PRIME protein model. .......................................................................... 14 Figure 2-4 Diagram of hydrogen bonding interaction in PRIME model. .................................................. 18 Figure 2-5 Schematic diagram of all-atom protein model. ....................................................................... 20 Figure 2-6 Potential profiles in all-atom model. ....................................................................................... 20 Figure 2-7 Diagram of hydrogen bonding interaction in all-atom model. ............................................... 21 Figure 2-8 Schematic diagram of wall interaction in a cylindrical tube or a spherical cavity. ................. 23 Figure 3-1 Determination of the folding temperature of A 10 from a plot of the potential energy U versus dimensionless temperature T. .................................................................................................................. 27 Figure 3-2 Examples of spherical and nearly spherical particles that are formed at a very early stage of the aggregation. ........................................................................................................................................ 28 Figure 3-3 Four stages of the aggregation process of 16 A 30 peptides. ................................................... 29 Figure 3-4 The structures of the β-sheet (left) and fibrils (right). ............................................................ 30 Figure 3-5 The Ramachandran plot in terms of <U>. ............................................................................... 31 Figure 3-6 Dynamic evolution of the total number of aggregates N(t) formed. ...................................... 32 n is the number of peptides initially in the system. ................................................................................. 32 Figure 3-7 Rescaled plot of the data in Figure 3-6 with X=n/t 0.77 and Y=t/N 0.60 . ...................................... 34 Figure 3-8 Dynamic evolution of the total number of aggregates N(t) formed for the 48-peptide system. ................................................................................................................................................................... 34 Figure 3-9 Dynamic evolution of the mean aggregate size S(t). ............................................................... 36 Figure 3-10 Dynamic evolution of the mean aggregate size S(t) in the 48-peptide system. ................... 36 Figure 3-11 The large aggregate formed by 16 A 30 at 360 K. ................................................................... 38 Figure 3-12 Same system as of Figure 3-11, but at 420 K......................................................................... 39 viii Figure 3-13 Dynamic evolution of N β (t), the number of peptides that are involved in the β-sheet formation. ................................................................................................................................................. 41 Figure 3-14 Rescaled plot of the data of Figure 3-13 with X=nt 0.75 and Y=tN β . ........................................ 42 Figure 3-15 Dynamic evolution of the numbers of various types of hydrogen bonds formed in the systems initially containing (a) 8 (b) 16 (c) 32 and (d) 48 A 10 , and (e) 16 A 30 . ........................................................ 43 Figure 3-16 Dynamic evolution of radius of gyration R g of aggregates of various sizes in the 48-peptide system. ...................................................................................................................................................... 44 Figure 3-17 Dependence of the aggregates’ diffusivity D on their mass M. ............................................ 45 Figure 3-18 Logarithmic plot of the dependence of the diffusivity D on the aggregates mass M. .......... 46 Figure 4-1 Four stages of the aggregation of 16 A 30 peptides in a spherical cavity, whose diameter is 1.5 times of the peptide length when they completely stretch out. ............................................................. 54 Figure 4-2 Comparison of the aggregates of 10 A 30 peptides in a small spherical cavity whose diameter is the same as the peptides’ length in their most stretched form. .............................................................. 55 Figure 4-3 Four stages of the aggregation of 16 A 30 peptides in Tube3 (cylindrical tube with aspect ratio L/D=3). ....................................................................................................................................................... 56 Figure 4-4 Dynamic evolution of the total number of aggregates, N(t). .................................................. 57 Figure 4-5 Comparison of dynamic evolution of the total number aggregate N(t) of the A 30 peptides in various geometries. .................................................................................................................................. 58 Figure 4-6 The universal collapsed curve for the total numbers of aggregates N(t), with X=n/t z’ and Y=t/N ω . ................................................................................................................................................................... 59 Figure 4-7 Power-law scaling of the total number of aggregates N(t) for the aggregation of 48 A 10 peptides. ................................................................................................................................................................... 61 Figure 4-8 Power-law scaling of the mean aggregate size S(t) during aggregation of 48 A 10 peptides. .. 63 Figure 4-9 Comparison of the dynamic evolution of the number of peptides N β involved in the formation of β-sheets during aggregation of 16 A 30 peptides in various environments. ......................................... 64 Figure 4-10 The universal collapsed curve for N β (t), the number of peptides that form β-sheets, with Y=tN β ω versus X=nt z ................................................................................................................................... 65 Figure 4-11 Definition of the dihedral angles 𝜑 and 𝜓 . ........................................................................... 66 ix Figure 4-12 Dynamic evolution of the number of various types of hydrogen bonds formed during aggregation of the A 10 peptides in various geometries. ........................................................................... 67 Figure 4-13. Same as in Figure 4-12, but for the hydrogen bonds formed during the aggregation of A 30 peptides. ................................................................................................................................................... 68 Figure 4-14 Comparison of the effect of peptide size ℓ on the number of peptides N β that are involved in the formation of β-sheets. ........................................................................................................................ 69 Figure 4-15 Comparison of the effect of peptide size ℓ on the total number of aggregates N(t). .......... 69 Figure 4-16 Effect of the diameter D of Tube3 on the total number of aggregates N(t) and the number of peptides N β (t) that are involved in the formation of β-sheets. ................................................................ 71 Figure 4-17 Same as in Figure 4-16, but in a spherical geometry. ........................................................... 72 Figure 4-18 Logarithmic plots of the direct correlation functions C(r) for the aggregates formed by A 10 and A 30 peptides ....................................................................................................................................... 74 Figure 5-1 Folding temperatures of (GA) 25 and (GA) 50 ............................................................................. 80 Figure 5-2 Aggregation of (GA) 25 and (GA) 50 at different temperatures .................................................. 82 Figure 5-3 Four stages of the aggregation of 10 (GA) 25 ............................................................................ 83 Figure 5-4 Detail view of an ordered β-sheet aggregate. ......................................................................... 84 Figure 5-5 Dynamic evolution of the total number of aggregates N(t) formed ....................................... 85 Figure 5-6 Different secondary structures and the corresponding Ramachandran plots ........................ 86 Figure 5-7 Diagram of dihedral angels 𝜑 and 𝜓 in a peptide. .................................................................. 87 Figure 5-8 HB contact map including different secondary structures ...................................................... 88 Figure 5-9 Structure details of different secondary structures. ............................................................... 89 Figure 5-10 Potential energy and hydrogen bond number of each type of structures ........................... 90 Figure 5-11 Identical structures of β-helix (left) and double-helix (right) formed in poly-GR system. .... 92 Figure 6-1 Sample system of a collision of a pair of atoms (2D) .............................................................. 93 Figure 6-2 Diagram of the two roots for the solution to the quadratic formula used to calculate the future collision time. ............................................................................................................................................ 95 Figure 6-3 Diagram of the different collision situations. .......................................................................... 95 Figure 6-4 Flowsheet of soft-sphere collision interactions in DMD. 183 .................................................... 98 x ABSTRACT This dissertation presents the studies of dynamics of protein aggregation in (i) unbounded media, i.e. bulk solutions, and (ii) confined environments, e.g. cylindrical tubes and spherical cavities. Computer simulation is used to investigate the aggregation problem. An intermediate-resolution coarse-grained model of proteins as well as a high-resolution all-atom model of proteins are employed together with discontinuous molecular dynamics (DMD). This combination is approved to be fast and accurate to simulate the aggregation process of multiple-protein systems. We first investigate the aggregation of multiple polyalanine peptides in unbounded media. Different numbers of identical polyalanine peptides are simulated to mimic β-amyloid aggregation in bulk solutions. In all cases we observe fibril-like aggregates, which have been implicated in many neuro- degenerative diseases, such as Alzheimer's, Parkinson's, amyotrophic lateral sclerosis, and prion. Various properties of the aggregates have been computed, including dynamic evolution of the total number of aggregates, mean aggregate size, number of the peptides that contribute to the formation of β-sheets, number of various types of hydrogen bonds formed in the systems, radius of gyration of the aggregates, and diffusivity of the aggregates. These results show that the aggregation process follows a nucleation scenario, which is also observed in experiments. Many of the quantities follow dynamic scaling, similar to those found in the aggregation of colloidal clusters. In particular, at long times the mean aggregate size 𝑆 (𝑡 ) grows with time as 𝑆 (𝑡 )~𝑡 𝑧 , where 𝑧 is the dynamic exponent. The existence of such dynamic scaling was recently confirmed by experiments. Next, we investigate the aggregation of the same polyalanine peptides but in confined environments. The medium, as a model of a crowded cellular environment, is represented by a spherical cavity, as well as cylindrical tubes with two aspect ratios. The aggregation process leads to the formation of β-sheets and eventually fibrils. Several important properties of the aggregation process, including dynamic evolution of the total number of aggregates, mean aggregate size, and number of peptides that contribute to the formation of β-sheets, have been computed. We show, similar to the unconfined media, xi that the computed properties follow dynamic scaling, characterized by power laws. The exponents that characterize the power-law in spherical cavities are shown to agree with those in unbounded media in the same protein density, while the exponents for the aggregation in cylindrical tubes exhibit sensitivity to the geometry of the system. The effects on aggregation from the protein size, i.e. the number of amino acids in a protein, as well as the size of the confined media, have also been studied. Similarities and differences between the aggregations in confined and unconfined media are described. Last, we use DMD combined with an all-atom model of proteins to investigate the aggregation of poly- glycine-alanine in unbounded media. The simulation shows poly-glycine-alanine is an aggregation-prone protein and will eventually form β-sheet-rich aggregates, which matches the experiment observations. Dynamic evolution of the total number of aggregates has been calculated. In the simulations, we observe two unique helical structures, β-helix and double-helix, frequently form in different configurations. We investigate their structures in detail and the analysis results show that such helical structures are very stable and share the characteristics of both α-helix and β-sheet. Identical phenomena are also observed in poly-glycine-arginine systems. We propose that proteins having the sequence of (GX) n , where X could be any non-glycine amino acid and n is the repeat length, would share the same folding structures of β- helix and double-helix, and it is the glycine in the repeat that contributes to this characteristic. 1 1. INTRODUCTION 1.1 OVERVIEW OF PROTEIN AGGREGATION Proteins are among the most important biomolecules in lives. They perform a variety of essential functions: they are not only the basic building materials of our tissues, like collagen, they also help to transport molecules from cells to cells, like hemoglobin; when our bodies have attacks, the immunoglobulins in our immune systems will protect us from bacteria and virus invasions; in our daily metabolism, enzymes accelerate thousands of biochemical reactions and have our body work much more efficiently. We could list tons of more examples here. It is never exaggerated to say that there will be no life without proteins. A protein will perform its function once it properly folds into a native state — a unique, three- dimensional (3D) equilibrium structure. However, in some cases, proteins would not fold correctly. Depending on their environments, even properly folded proteins can spontaneously change their conformations and take on misfolded states. Such states give rise to an important biological phenomenon, protein aggregation. Protein aggregates may take on a variety of shapes and structures. They would form large structures in the form of protein plaques, which are often toxic and their deposition in tissues is believed to be a major contributing factor to many neuro-degenerative diseases, collectively known as amyloidosis, which includes Alzheimer's, Parkinson's, amyotrophic lateral sclerosis (ALS, or the Lou Gehrig's disease), and prion (transmissible spongiform encephalopathies) disease. There are more than 40 such diseases, many of which are fatal. 1-4 In addition, protein aggregation has taken on technological importance. Fabrication of new nanoscale materials using proteins is of much current interest. 5-6 Protein aggregation can also limit the yields of biologically active proteins at most, if not all stages of the manufacturing process. 7-9 During storage, temperature fluctuations, as well as the existence of air and solid interfaces in the container, may cause the proteins to be misfolded, a state prone to aggregate. The thereafter aggregates will reduce the efficiency of formulations and lead to enhanced immunogenicity. 9 During purification, for example, if E. coli is the chosen host for recombinant protein expression, the overexpressed proteins 2 often form intracellular aggregates that are called inclusion bodies (IBs). Although the IBs are relatively pure, they are not biologically active until refolded. But, as they refold, they compete with the aggregation process that can severely limit the refolding efficiency. 9-10 Despite the obvious importance of protein aggregation in both pathology and biotechnology, however, the understanding of its mechanism and physical basis are far from complete. 1.2 AGGREGATE STRUCTURE AND AGGREGATION MECHANISM In the past several decades, numerous works has been done intending to unveil the mysteries of protein aggregation. Two general categories can be summarized to describe the morphological features of protein aggregates: amorphous and fibrillar. Since amorphous aggregates lack of long-range order, they usually only have coarsely descriptions of the molecular organizations. They were originally thought to only contain unfolded proteins that held together by random associations between hydrophobic residues. In contrast, some aggregates retain native activity, and the proteins inside could even be recovered from inclusion bodies in mild denaturing conditions. 11 In addition, staining experiments show that some proteins directed to inclusion bodies contain at least some ordered structural elements common to fibrillar aggregates. 12-13 These suggest that amorphous aggregates are composed of a combination of native- and fibrillar-like structures. Fibrillar aggregates, however, contain long-range order and, as a result, have been described in greater detail than amorphous forms. 14-17 These aggregates have been shown to be composed of substantial amounts of β-sheet structures and lacking in helical content. For example, high-resolution X-ray diffraction and solid-state NMR data reveal that the fibrils formed by 40-residue β-amyloid peptides (Aβ 1- 40 ), which associated with Alzheimer’s disease, have a cross-β structural organization that has not been observed in globular proteins 18 , see Figure 1-1. This cross-β organization has actually been observed in the aggregates of several unrelated amyloid peptides, suggesting this is a regular structural element in fibrillary aggregates. 19-21 3 Figure 1-1 Cartoon representation of the quaternary structure of an Aβ 1-40 fibril, viewed parallel and perpendicular to z axis. Each peptide has two β-strands (red and blue) separated by a turn (green). The neighboring peptides form hydrogen bonds with each other along the strands (red to red and blue to blue), creating a β-sheet folded back on itself. The two folded sheets are stacked together at their C terminus to form the fibril 18 . Four mechanisms have been proposed to describe the conformational transformation and assembly of normally soluble proteins into ordered aggregates 22 . The first is the templated assembly mechanism (TA), 23-24 in which a soluble random-coil peptide binds to a preassembled β-sheet-rich nucleus and then undergoes a rate-determining structural change. The second is the monomer-directed conversion mechanism (MDC), 25 in which a misfolded peptide, which adopts a conformation that is identical to that found in a fibril, binds to another soluble monomer in a rate-determining step, converts that monomer to the same misfolded conformation and then dissociates, with both peptides rapidly add to the end of a growing fibril. The third is the nucleated polymerization mechanism (NP), 26-28 which is characterized by the rate-limiting formation of a nucleus followed by the rapid addition of monomers to the growing end of the nucleus. The fourth is the nucleated-conformational conversion mechanism (NCC), 29 in which the formation of small amorphous aggregates precedes the rate-limiting formation of a critical nucleus and subsequent rapid growth of large fibrils through the addition of globular multimers to the fibril ends. All the above-mentioned mechanisms require that the system goes through a nucleation event before fibrils can be formed: the nucleus in MDC is a monomeric peptide, whereas the nuclei in TA, NP and NCC are 4 oligomers. Nguyen and Hall’s MD simulation study 30 shows that the combination of three of the four mechanisms, TA, NP and NCC, could be used to describe the simulation kinetics. 1.3 MOLECULAR SIMULATION AND PROTEIN MODEL Undoubtedly, experiments 1, 13, 31-37 have been playing a vital role in gaining deeper understanding of protein aggregation and their deposition on the surface of tissues. However, considering that protein aggregation is a rapid process, during which the peptides or proteins undergo rapid conformational changes; when coupled with the solvent effects, one faces significant challenges in understanding the properties of the aggregates’ structures in experimental studies. Thanks to the revolutionary development in computing power of computers, molecular simulation study of the aggregation process has become highly valuable, not mention the fact that including or deleting a particular effect is much straightforward in computer simulations. For this reason, computational studies, especially molecular simulations, provide another important tool to study protein aggregation phenomenon and gain insights that may be too difficult to obtain by experiments alone. In molecular simulations, there are two major techniques used: Monte Carlo (MC) and molecular dynamics (MD). MC methods are quite useful for simulating systems with many coupled degrees of freedom. During simulations, atoms are moved random distances in random directions, and the potential energy of each sampled configuration is calculated. If these motions lead a decrease in potential energy, then they are accepted. Otherwise, a random number will be generated between zero and one and compared with the probability given by the Boltzmann factor, 𝑒 −∆𝜀 /𝑘𝑇 , where ∆𝘀 is the potential energy change before and after the motions, 𝑘 is the Boltzmann constant and 𝑇 is the temperature. If the random number is smaller, then the motions are accepted, otherwise rejected. This process will be repeated until the simulation finishes. The configurational energy is record for further analysis. MC methods can be very fast; however, a limitation is that motions are difficult to simulate for dense systems or molecular structures that are highly complex or large. 5 Classic MD, instead, implements Newton’s equations of motion to study the physical movements of atoms. At each time step during simulation, the pairwise forces are summed to determine the net force on each atom. The acceleration vector is then determined. Based on the current position and velocity of the atom, new position and velocity at the next time step can be calculated by an algorithm, like Verlet or leapfrog algorithm. Repeat this process until the required time scale is reached. Atoms’ coordinates and velocities are recorded based on the dump rate for further analysis. However, to avoid data overflow, the time step has to be very small, usually 1~2 fs; also, the computation of the pairwise forces and potential energy on each atom is expensive: MD could be very time consuming to simulate big systems. Here we introduce an alternative of the classic MD, discontinuous molecular dynamics (DMD), which is an event-driven MD method. It was first introduced in 1959 by Alder and Wainwright 38 for simulations of hard spheres. Later it was used by Rapaport 39-41 for simulation of polymer chains. Nowadays it was adopted for simulations of protein-like polymers. It is extremely fast and suitable for the simulations of large systems on long time scales. In our studies, we combine DMD method with different protein models to perform the simulations. We will come back and talk about DMD in detail in Chapter 2. Different molecular models have been developed to represent proteins in molecular simulations, ranging from lattice one-bead coarse-grained model to off-lattice all-atom model. We can divide them into three categories: high-resolution models, low-resolution models and intermediate-resolution models. High resolution models are based on a realistic representation of protein geometry and energetics. They represent the proteins with all-atom details and are widely used in many popular simulation packages containing well-tested force fields such as AMBER, 42-43 CHARMM, 44 GROMACS, 45-46 and DISCOVER. 47 They usually explicitly represent every solvent molecules as well. High resolution, of course, provides high details, but it also makes the simulation extremely expensive, which precludes its application to problems involving large conformational changes or longtime scales. Low-resolution models, instead, are based on a coarse-grained representation of protein geometry and energetics. They usually simplify atoms into a functional group and use a single bead to reflect the composite properties of the group. Solvent is often implicitly represented in order to enhance computational efficiency. The potential functions are also low-resolution, including the following three categories: (1) Go-type potentials 48-50 whose parameters are chosen to favor the protein’s known native 6 state, (2) potentials based on the relative hydrophobicity of the side chains, and (3) knowledge-based potentials in which statistical data from the Protein Data Bank on residue-residue contacts are used to infer side-chain/side-chain potentials. Despite the simplicity, low-resolution models have provided valuable insights into the basic principles of protein folding and aggregation due to their ability to monitor large conformational changes for longtime scales. Their weakness is, of course, their inability to make definitive statements about the folding/aggregation of specific proteins. Obviously high-resolution and low-resolution models represent two ends of a spectrum of models with varying levels of details. Which type of model to use depends on what situations are being investigated. Trade-offs will be must: the larger the number of degrees of freedom to be simulated, the simpler the model must be. To find the balance between the high-resolution and low-resolution models, people gradually add more beads to the low-resolution models, leading to the third category, intermediate- resolution models. Inside this category, various levels of details of protein representations were designed, including two-bead approaches (one backbone bead and one side-chain bead 51-55 ), three-bead approaches (one backbone bead and two side-chain beads 56 ), four-bead approaches, 57 five-bead approaches, 58-64 and seven-bead approaches. 65-66 Energy functions for the intermediate-resolution models include not only those available in low-resolution models as mentioned above, but also hydrogen-bonding potentials, multibody terms, burial terms (in which the strength of the hydrophobic interaction depends on the extent of burial), and special potentials for disulfide bonds and proline. 1.4 SUMMARY OF MOLECULAR SIMULATION STUDIES ON PROTEIN AGGREGATION At early time, most simulation studies on protein aggregations employed low-resolution models and MC methods. Gupta and Hall 67 performed dynamic MC simulations on lattice multi-chain systems containing 20 to 40 of protein chains having the sequence of (PHPPHPPHHPPHHPPHPPHP) (H: hydrophobic residue, P: polar residue). They found that aggregation arises from association among partially folded intermediates, instead of random coil states. 7 Nguyen and Hall 68 then conducted the same simulations on 50 of the peptides as Gupta and Hall and investigated how changing the rate of chemical or thermal renaturation affects the folding and aggregation behavior of the peptides. Their predictions regarding different renaturation protocols are supported by experimental results on the lysozyme system. 69-71 Ding et al. 72 applied an off-lattice model combined with two-bead approach and Go potentials to study the aggregation of a system containing 8 Src SH3 domain proteins. They observed a fibrillar double β- sheet structure formed by packing the flexible segments (called RT-loops) from different proteins. The structure had inter-β-strand spacing distances that are similar to those observed in experiments. The researches above all used low-resolution models. Instead, most computer simulations using high- resolution protein models have been devoted to the study of systems containing, however, already- formed amyloid fibrils. Li et al. 73 tested the stability of a fibril structure based on synchrotron X-ray studies. MD simulations in explicit solvent were performed using AMBER 4.1 package with the Cornell et al. 74 force field on a 16- mer untwisted β-sheet stack for 2.5 nanoseconds, and on a 48-mer protofilament for 600 picoseconds. They found that the 16-mer structure remains stable in the original untwisted conformation, while the 48-mer structure lessens its original 15°-twist over time. George et al. 75 conducted MC simulations β-amyloid peptide using the CHARMM package to determine the difference between inter-peptide interactions in Aβ(1-40) and Aβ(1-43). They found well-defined enhanced side-chain interactions between residues 41, 42, and 43 and residues on adjacent molecules. They also examined the interaction between the surface of an Aβ(1–43) fibril and a proposed aggregation-inhibitor, the molecule IDOX. Liotta and co-workers 76-77 performed MD simulations using the AMBER package to study the stability of amyloid fibrils containing 8 β-sheets, each with 14 Aβ(10–35) fragments. They found large spatial and temporal fluctuations about a central core whose “micelle” or “molten-globule-like” properties served to lower the entropic penalty for fibril formation, and thus to enhance the overall stability of the fibril structure. 8 High-resolution simulation studies of the formation of fibrils in systems containing random coils have been conducted but they are generally limited to early events in the assembly of a few peptides. Gsponer et al. 78 performed MD simulations on a system containing 3 heptapeptides, GNNQQNY, from the yeast prion Sup35 using the CHARMM program and observed the formation of a β-sheet. Mager and co- workers 79 conducted MD simulations on a system containing 2 to 4 Aβ(1–42) peptides using the AMBER program and observed the formation of a β-sheet. Klimov and Thirumalai 80 conducted long MD simulations on systems containing three Aβ(16–22) peptides using the MOIL program and they observed the formation of an anti-parallel β-sheet through an obligatory α-helical intermediate. Although the abovementioned studies offer considerable insight into the properties of fibrils and the mechanisms by which they are formed, the systems considered do not contain enough peptides to mimic the nucleus that stabilizes the large fibrils observed in experiments. Given current computational capabilities, simpler models are required to simulate spontaneous fibril formation in multiprotein systems. Intermediate-resolution protein models, which are essentially a compromise between the simplified chain models described above and the detailed all-atom models like used in AMBER, CHARMM, and GROMACS, have been used in recent years to simulate the aggregation of many-protein systems. Hall and co-workers used an intermediate-resolution protein model developed by Smith and Hall, 81-84 named PRIME, to study the aggregation of large systems by applying DMD simulations. 30, 85-88 They conducted equilibrium simulations on systems containing many polyalanine peptides and examined the thermodynamics of fibril formation. Using the same system, they then did further investigation on the kinetics of fibril formation. They simulated the whole process of fibril formation from the random coil state to the fibril state. They found that fibril formation of polyalanines incorporate features that combine the characteristics of three experimental models: the templated assembly, nucleated polymerization, and nucleated conformational conversion models. 30 They also used their improved PRIME model to simulate the aggregations of different poly-peptides consisted of alanine, glycine and valine. They showed that subtle changes in sequence will result in slightly different fibril structures. They also found that the ability to form fibrils will increase as the chain length and the length of the stretch of hydrophobic residues increase. 9 Ding et al. 89 and Urbanc et al. 90-91 used the similar intermediate-resolution protein model combined with DMD to simulate protein structure transition, folding and aggregation. Their model is also four-beads coarse-grained model (three on backbone and one on side-chain), but has slightly different values on bead diameters, bond lengths and energy parameters. Ding et al. 92 also developed a more detailed model and added additional effective side-chain atoms into the four-bead model. After that, they further improved their model to all-atom resolution. 93 Although the all-atom model limits the simulation time scale, the combination of DMD and the all-atom model allowed them to achieve the full folding process of different proteins from ab initio states to native or near-native states. This really highlights the potential power of DMD to simulate large systems for relatively longtime scale. 1.5 PROTEIN AGGREGATION IN CELLULAR ENVIRONMENT All the studies mentioned above are about the protein aggregation in bulk solutions. However, as far as the role of protein aggregation in the development of the aforementioned diseases is concerned, the phenomenon occurs in a heterogeneous and highly crowded cellular environment, which, more generally, we refer to as a confined medium. The arising excluded volume interactions either due to macromolecular crowding, or localization of the protein in a small volume, or the extra interactions between protein and the confining boundary can have significant effects on protein folding and aggregation. Thus, it is well-known and well-understood that the properties of many phenomena that occur in a confined medium are vastly different from those in an unbounded medium or bulk. Experimental and theoretical studies 94-97 have reported that confinement will enhance the stability of the native state of protein instead of the non-native conformations. Simulation researches 98-103 also showed the same phenomenon that confinement of denatured proteins in the limited space accelerates folding when compared with that in bulk. Other studies about protein aggregation in confined environments, however, reported that in some circumstances, the “crowding” environment will accelerate the formation of large aggregates for both native 104-108 and non-native 96, 109 aggregation. The competition between the accelerations of folding and aggregating in confined media will be great interesting to be investigated. The effect on the non-native aggregation, especially amyloids, deserves more attentions. 10 The structure of this dissertation will be as follows: in the next chapter, we will talk about DMD method in detail as well as two protein models, one is an intermediate-resolution coarse-grained, called PRIME, and the other is a high-resolution all-atom. Chapter 3 110 will present the studies of dynamics of protein aggregation in unbounded media, i.e. in bulk solutions. Chapter 4 111 will present the studies of dynamics of protein aggregation in confined media, e.g. cylindrical tubes and spherical cavities. Chapter 5 will present the molecular dynamics study of the structure of poly-glycine-alanine in aggregation and folding. The last part will be an appendix describes the derivation of the algorithms used in DMD method. 11 2. DISCONTINUOUS MOLECULAR DYNAMICS AND MOLECULAR MODEL 2.1 DISCONTINUOUS MOLECULAR DYNAMICS Discontinuous molecular dynamics (DMD) is an extremely fast alternative to the traditional molecular dynamics. It was first introduced in 1959 by Alder and Wainwright 38 for simulations of hard spheres. Later it was used by Rapaport 39-41 for simulations of polymer chains. Nowadays it was adopted for simulations of protein-like polymers. It is extremely fast and suitable for the simulations of large systems on longtime scales. During DMD simulation, potentials applied on particles (coarse-grained beads or atoms) are approximated by discontinuous step-functions of inter-particle distance 𝑟 (see Figure 2-1), instead of continuous potential functions applied in the traditional MD. Thus, no forces will be exerted on particles until their pairwise distance becomes equal to the point of a discontinuity on the potential profile. When particles encounter a potential discontinuity, we call that they have an “event”. Between events, particles move at constant velocities. Because the energy change during each event is known, the post- event velocities can be calculated by solving the conservation of momentum and the conservation of energy simultaneously. Thus, the trajectory of a particle can be simulated discontinuously between events. The detailed derivations of the equations of motion can be found in Appendix. 12 Figure 2-1 Potential profiles in DMD (a) hard-sphere collision (b) attractive square well interaction (c) repulsive soft sphere interaction (d) covalent bond and auxiliary bond interaction The simulation code for DMD is significantly different from that for the traditional MD (see Figure 2-2) due to the discontinuous movement of particles from one event to the next at known velocities, and the algorithm could become much more complex due to the event linking, time scheduling and especially the explicit hydrogen bonding mechanism. r U(r) hard sphere r U(r) square well r U(r) shoulder soft sphere r U(r) covalent bond auxiliary bond 13 Figure 2-2 Flow chart of simulation process The program always keeps a set of raw data that stores the basic information of each particle, including the static information, like sequence, type, mass etc., and the dynamic information, like coordinate, velocity etc. Nothing will be modified until an event is committed. In each step, the program will make a copy of the data of the target atom(s) involving in the nearest event and do the calculations only on this copy, i.e. the raw data will not be touched. Since DMD is an event-driven simulation, the results of the previous event would have voided the prediction of the current. Thus, any hazards that would have affected the current event will be evaluated. If there exist any hazards, the prediction of the current event will be voided. In this case, the program will re-predict the event for the target atom(s). Otherwise, the prediction for the current event will be committed and then the program will predict a new future event for the target atom(s). After that, the raw data will be updated and the time tree (stores the event sequence) will be re-scheduled. The program will keep looping this process until finishing the assigned simulation time. 14 2.2 MOLECULAR MODEL Different protein models can be used in DMD ranging from low-resolution coarse-grained to high- resolution all-atom. 2.2.1 INTERMEDIATE-RESOLUTION COARSE-GRAINED MODEL In the early stage of our studies, we use an intermediate-resolution coarse-grained model, named PRIME, developed by Smith and Hall 81-84 . It correctly reproduces the backbone geometry and has been highly successful in reproducing several important and experimentally determined features of the proteins under bulk conditions. In PRIME, every amino acid is represented by four beads: N, C α , C, and R (see Figure 2-3). The N bead represents the amide nitrogen and its hydrogen; the C α bead represents the α- carbon and its hydrogen; the C bead represents the carbonyl carbon and oxygen; the R bead represents the side chain, which are assumed to have the same diameter of CH 3 (the size chain of alanine). The peptide bond is assumed to be in the trans configuration. All bonds’ lengths and angles, as well as the distance between consecutive C α beads, are fixed at their experimentally measured values. Table 2-1 presents all the relevant parameters of the model. Solvent molecules are not explicitly included in the model. The effect of solvent is factored into the energy function as a potential of mean force. Figure 2-3 Schematic diagram of PRIME protein model. Narrow solid lines represent covalent bonds and wide dashed lines represent auxiliary bonds. The latter helps to maintain the correct backbone geometry. 30 15 Table 2-1 Parameter values of PRIME model Beads Diameter (Å) N 3.300 C α 3.700 C 4.000 R 4.408 Bonds Length (Å) N i -C α,i 1.460 C α,i -N i+1 1.510 C i -N i+1 1.330 C α,i -R i 1.531 Auxiliary Bonds Length (Å) N i -C i 2.450 C α,i -N i+1 2.410 C i -C α,i+1 2.450 C α,i -C a,i+1 3.800 N i -R i 2.440 C i -R i 2.490 Bond Angles Angle (deg) ∠N i -C α,i -C i 111.0 ∠C α,i -C i -N i+1 116.0 ∠C i -N i+1 -C α,i+1 122.0 ∠R i -C α,i -N i 109.6 ∠R i -C α,i -C i 110.1 Beads in the protein are subject to five different types of forces during events: (1) infinite repulsion due to excluded volume effect (hard-sphere collision), (2) infinite attraction between the beads connected with covalent bonds and auxiliary bonds, (3) attraction between pairs of backbone beads during hydrogen bond formation, (4) attraction between hydrophobic side chains, and (5) repulsion from the neighboring beads of hydrogen bond (soft-sphere collision). All these forces can be represented by the following step-functions, 16 𝑈 𝑖𝑗 (𝑟 ) = { ∞ 𝑟 ≤ 𝜎 𝑖𝑗 0 𝑟 > 𝜎 𝑖𝑗 (Eqn 2-1) This is for a hard-sphere potential. 𝑟 is the distance between beads 𝑖 and 𝑗 ; 𝜎 𝑖𝑗 is the radius sum of beads 𝑖 and 𝑗 . 𝑈 𝑖𝑗 (𝑟 ) = { ∞ 𝑟 ≤ 𝜎 𝑖𝑗 −𝘀 𝜎 𝑖𝑗 < 𝑟 ≤ 𝜆𝜎 𝑖𝑗 0 𝑟 > 𝜆𝜎 𝑖𝑗 (Eqn 2-2) This is for a square-well or shoulder potential depending on the sign of 𝘀 (positive for well depth, attractive; negative for shoulder height, repulsive). 𝜆𝜎 𝑖𝑗 is the well/shoulder width. 𝑈 𝑖𝑗 (𝑟 ) = { ∞ 𝑟 ≤ 𝑙 (1 − 𝛿 ) 0 𝑙 (1 − 𝛿 ) < 𝑟 < 𝑙 (1 + 𝛿 ) ∞ 𝑟 ≥ 𝑙 (1 + 𝛿 ) (Eqn 2-3) This is for a bond interaction. 𝑙 is the bond length and 𝛿 is the bond vibration tolerance. Note that, as explained in the work by Takada et al. 57 , interaction involving a pair of beads that are very close on the chain are more faithfully represented by an interaction between the atoms themselves, not the united atoms or the coarse-grained beads. In other words, for example, neighboring carbonyl carbon atoms experience interactions based on the diameter of a carbon atom, not the diameter of C bead. Consequently, for interactions between pairs of beads that separated by three or fewer bonds, the beads are allowed to overlap each other by up to 25% of their diameters; for those separated by four bonds, they are allowed to overlap up to 15%. Backbone beads and C α and R beads are connected by covalent bonds. Bonds are allowed to move freely over a small range between 𝑙 (1 − 𝛿 ) and 𝑙 (1 + 𝛿 ). The choice of 𝛿 defines the acceptable range of bond vibration. In the simulation, 𝛿 is chosen to be 0.02375. Auxiliary bonds are used between next- neighboring beads along the backbone chain to hold backbone angles fixed. For example, an auxiliary bond is placed between N i and C i in Figure (2-3) to maintain N i -C a,i -C i angle to be near its ideal value. Auxiliary bonds are also placed between consecutive pairs of C α beads to maintain their distances near the experimentally determined value. The combination of C a,i -N i+1 , C i -C a,i+1 , and C a,i -C a,i+1 auxiliary bonds and the C a,i -C i , C i -N i+1 , and N i+1 -C a,i+1 covalent bonds holds the inter-amino acid group (C a,i -C i -N i+1 -C a,i+1 ) in 17 the trans configuration. Auxiliary bonds are also placed between R bead on side-chain and N and C beads on backbone. This is to hold the side-chain bead fixed relative to the backbone so that the residue is a L- isomer. Hydrogen bond (HB) can be formed between N and C beads on backbones. Four criteria have to be met (see Figure 2-4): (1) the distance between N i and C j beads is 4.2 Å; (2) the neighboring beads, C i-1 and C a,i , C a,j and N j+1 , are at appropriate positions; (3) neither N i nor C j bead is already involved in a hydrogen bond; (4) N i and C j beads are separated by at least three intervening residues or they are in different peptides. To ensure criteria (2) is satisfied, it requires that the four bead-pairs N i -C a,j , N i -N j+1 , C j -C a,i , and C j -C i-1 shown in Figure (2-4) are separated by a distance greater than 𝑑 𝑖𝑗 , whose values (given in Table 2- 2) are chosen to maintain the hydrogen bond in its ideal orientation. When two beads N i and C j approach each other at the distance 𝑑 𝑖𝑗 = 4.2Å, the program will evaluate the total potential energy change by checking the five auxiliary interactions. The potential energy change can be −𝘀 HB , 0, 𝘀 HB , 2𝘀 HB and 3𝘀 HB , depending on the orientations of N i , C j , and their neighboring beads. If the kinetic energy is enough to overcome the total potential energy change, the forward reaction will happen. N i and C j beads will form a hydrogen bond and change their types into N i ’ and C j ’, respectively, and interact with each other according to the parameters of their new types (the interactions between N i /C j bead and the other beads will not be affected). Otherwise, the two beads N i and C j will not change their types and instead may undergo a hard-sphere collision in the near future. The thermal fluctuation distorts the orientation of the hydrogen bond and large fluctuations may break the hydrogen bond. Once the two atoms N i ’ and C j ’ come again to the exact distance 𝑑 𝑖𝑗 = 4.2Å, a reverse reaction may happen. We will evaluate the potential energy change by using the similar way during HB formation, and decide if the HB would break such that N i ’ and C j ’ will change back to their original types, N i and C j , respectively. 18 Figure 2-4 Diagram of hydrogen bonding interaction in PRIME model. Solid lines indicate the covalent and auxiliary bonds; solid arrows indicate the hydrogen bond and its orientation; dashed lines indicate the auxiliary bonds helping to maintain the orientation. 112 Table 2-2 Parameter values for hydrogen bond Pairs 𝒅 𝒊 𝒋 (Å) N i -C α,j 5.00 N i -N j+1 4.74 C j -C α,i 4.86 C j -C i-1 4.83 The DMD simulations performed in this dissertation all follow canonical ensemble, i.e. NVT ensemble (constant system size, volume and temperature). To maintain the temperature around the target value, we implement Andersen thermostat. 113 Within this procedure, all beads in the system are subjected to random, infrequent collisions with ghost particles having the same mass as themselves. The post-event velocity of a bead colliding with a ghost particle is chosen randomly from a Maxwell-Boltzmann 19 distribution centered at the simulation temperature. The collision frequency is controlled at 1% of the total number of events. In DMD simulations, the solvent molecules are not explicitly included. They are modeled implicitly by means of a square well attraction (potential of mean force) between two hydrophobic side-chain beads. The side-chain beads, R beads, are assigned as either “H” or “P”: “H” means it is hydrophobic and such R beads will attract with each other when their distance is equal to 1.5𝜎 𝑖 𝑗 as long as they are separated by at least three intervening residues in the same peptide or they are in the different peptides; “P” means the side-chain bead is polar. The interactions between two polar beads or between a polar and a hydrophobic bead will be a regular hard-sphere collision. 2.2.2 HIGH-RESOLUTION ALL-ATOM MODEL In order to apply DMD to more complex systems, we upgraded the molecular model and potential library from coarse-grained to all-atom based on the work of Ding et el 93 . Here we would like to give our special thanks to Dr. Ding for providing the latest force field data (the parameters on the paper were outdated). The basic concepts and algorithm used in the all-atom model are similar to those in the coarse-grained model: it still uses auxiliary bonds to maintain the configuration of molecules, and it still uses the step functions to represent the forces. But after the updates, the model can represent all the atoms explicitly (except for the light hydrogens), see Figure 2-5; the step functions are also more realistic and closer to the continuous profile, see Figure 2-6 (a), which can handle both attractive and repulsive forces simultaneously, unlike the coarse-grained model that only has either one of them. The update also adds a new type of force: the dihedral constraint. Although it is also maintained by an auxiliary bond, it has more complex profile than a regular bond, see Figure 2-6 (b). 20 Figure 2-5 Schematic diagram of all-atom protein model. The solid lines represent the covalent bonds and the dashed lines represents the auxiliary bonds. 93 Figure 2-6 Potential profiles in all-atom model. (a) Non-bonded interactions. The continuous dashed curve corresponds to the VDW and solvation interactions between a pair of atoms. The solid step lines represent the discretized potential. (b) Potential profile of a dihedral constraint. 93 21 Since the all-atom model represents the oxygen and hydrogen atoms explicitly, the structure of hydrogen bond is also optimized. The new hydrogen bond structure is shown in Figure 2-7. Figure 2-7 Diagram of hydrogen bonding interaction in all-atom model. H i is the hydrogen atom; A j is the acceptor and could be different types of atoms; D i is the donor and X j is the heavy atom directly bonded to A j . The thick solid lines indicate the covalent bonds between atoms; the thin solid lines indicate the hydrogen bond; the dashed lines indicate the auxiliary bonds for helping maintain the hydrogen bond orientation. 93 By adjusting the types of A j atoms, the hydrogen bond can be formed between backbone & backbone, backbone & side chain, and side chain & side chain. The donors include backbone amide nitrogen and side chain nitrogen and oxygen in NH and OH groups; the acceptors include backbone carbonyl oxygen and side chain nitrogen and oxygen in CN and CO groups. The mechanisms of hydrogen bond forming and breaking are similar to the coarse-grained model: once the hydrogen, H i , and the acceptor, A j , reach the interaction distance, the distances between H i X j and D i A j , which define the orientation of hydrogen bond, will be evaluated; the total potential energy change, ΔE, between H i /A j and their neighboring atoms is also evaluated before and after the putative hydrogen bond formation. If these distances satisfy the preset range, and the total kinetic energy is enough to overcome the potential energy change, ΔE, the hydrogen bond will be allowed to form, and forbad otherwise. After the formation of hydrogen bond, the atom H i and A j will become H i ’ and A j ’, respectively, and they will interact with the other atoms by using the new pairs of potential. 22 2.2.3 CONFINED ENVIRONMENT MODEL To mimic the crowding environment in the real cell, we introduce three types of confined environments: (a) parallel pores, (b) cylindrical tubes and (c) spherical cavities. The confined wall is simulated as a smooth surface, which means that there are no explicit atoms on the walls. The wall surface is like a multi-layer shell, and each layer represents a potential step. To simulate parallel pore is very straightforward: once the position of a wall is set, the event time in the nearest future when an atom interact with the wall can be easily determined, and the after-interaction velocity can be computed directly: it is just a one-dimension problem and only the velocity perpendicular to the wall’s surface will be changed. However, for a cylindrical tube or a spherical cavity, the situation will become a little bit more complex. Due to the curvature of the surface and the implicit atoms, it would be difficult to determine the exact direction and position on the wall at which an atom will interact in the nearest event, both of which are required to compute the interaction potential and the after-interaction velocities. To deal with this issue, our DMD program uses an indirect way, see Figure 2-8: it is easy, instead, to compute the instant distance of an atom from the center of either a cylinder tube or a spherical cavity, L. Then the distance between the atom and the wall will be D = R – L, where R is the diameter of the tube or cavity. Thus, the potential region in the step function can be determined by D. Next, instead of doing a lot of trigonometric calculations or trying to find a pseudo-atom on the wall to model the interaction, we can consider this as an interaction between the target atom and a huge pseudo-atom at the center of the tube or cavity. In other words, any interaction between an atom and the wall will be like the atom inside another huge atom interacts internally with the shell of the huge atom. The interaction types, however, must be switched. For example, a well-capture event between the atom and the pseudo-atom will become a well-escape event between the atom and the wall; similarly, a hard-sphere collision event will become a well-bounce event, etc. Thus, the DMD program can use the same algorithm to compute both the interactions (a) between a pair of atoms and (b) between an atom and the wall. 23 Figure 2-8 Schematic diagram of wall interaction in a cylindrical tube or a spherical cavity. The big solid circle is the wall’s boundary, which is infinitely repulsive; the big dashed circle is the potential shell, either attractive or repulsive (the wall may have multiple layers of shells). The small black dashed circle is the original position of the target atom; the small black solid circle is the future position of the target atom after time t. The small red circle is the pseudo-atom at the center. The black line is the diameter of the wall, R; the red line is the distance between the target atom and the pseudo-atom, L; the blue line is the distance between the target atom and the wall, D = R – L. We develop the whole simulation package from scratch by C language, including the simulation program itself as well as the corresponding analysis programs. The package is published under MIT license and is free to download from https://github.com/alexzheng42/sDMD.git. 24 3. DYNAMICS OF PROTEIN AGGREGATION. I. UNIVERSAL SCALING IN UNBOUNDED MEDIA 3.1 INTRODUCTION In this chapter, we will report our simulation results for protein aggregation in an unbounded domain 110 . The motivation for doing so is twofold. One is to develop a benchmark for the aggregation phenomenon, so that the effect of confinement can be understood better. The second reason is that many of the quantities that we have computed have not been calculated before, but they shed light on the dynamics of the aggregation. The structure of this chapter is as follows. In Section 2 we describe the models of peptides as well as the simulation configurations we use in DMD. The results are presented in Section 3. In Section 4 we briefly describe an important phenomenon that contributes to the formation of neurotoxic oligomers after the formation of fibrils but is absent in our molecular simulations. The last section summarizes the chapter. 3.2 CONFIGURATION OF DMD SIMULATION 3.2.1 MODELS OF PEPTIDES Similar to our previous work on the effect of confinement on the stability and folding of proteins, 103, 114 as well as on their diffusion, 115 in the present work we also carry out molecular modeling and simulation of de novo designed α-family of proteins. 57-58, 61, 116-117 The model that we use is a mesoscopic one with a ten-residue sequence and one type of amino acid, each consisting of four beads. Thus, the amino acid sequence of the model is (AAAAAAAAAA), or A 10 . The Aβ, which is believed to be a major component of neurotic plaques in Alzheimer’s disease, typically contains 39-42 amino acids. Thus, we also simulate a much larger protein with a thirty-residue sequence of the same amino acid, or A 30 , in order to not only have a model that is closer to the known experimental information, but also test further the dynamic scaling theory that is proposed below for the aggregation of proteins. 25 Depending on their locations, the amino acids in the sequence may have different polarities, which is hydrophobic (H) or polar (P). The polarity sequence of the ten-residue model is PHHHHHHHHP, or (PH 8 P), while that of the thirty-residue model is (PH 28 P). Each peptide with ten residues folds 103 onto an α-helix with 6 hydrogen bonds, while those with 30 residues fold with 26 of such bonds. With all the amino acids being similar, they mimic the molecular model of an alanine. The configuration that we utilize here is different from what we used in our previous simulations of single proteins in confined media, 103, 114-115 and is intended to mimic the hydrophobic part of a β-amyloid (Aβ). The protein model that we utilize in the simulations is PRIME, 30, 84-85, 88, 112, 118-120 which has been described in detail in Chapter 2. 3.2.2 MOLECULAR SIMULATION In the simulations, energy is measured in units of 𝘀 HB , the strength of the HB interaction, and temperature is expressed in units of a temperature scale 𝑇 𝑠 such that, 𝒩 𝑘 𝐵 𝑇 𝑠 = 𝘀 HB , where 𝑘 𝐵 is the Boltzmann’s constant, and 𝒩 is the Avogadro’s number. In principle, 𝘀 HB depends on many parameters, including the type of the neighboring atoms and the structure of the molecule. In polyalanine, for example, 𝘀 HB varies between 3.5 and 8.6 kcal/mol. 121 MD simulation indicates that 𝘀 HB also depends on whether water is present in the system, and that 𝘀 HB for α-helices is between 1.9 and 5.6 kcal/mol. 122 Thus, we assumed 𝘀 HB = 6 kcal/mol (25 kJ/mol). This implies that 𝑇 𝑠 = 3000 K. The dimensionless temperature 𝑇 ∗ is then defined by 𝑇 ∗ = 𝑇 /𝑇 𝑠 . For convenience, we delete the * sign. All the simulations are carried out at 𝑇 = 1.1𝑇 𝑓 , which is above the folding temperature 𝑇 𝑓 of the peptides (see below). We note that many of the previous DMD simulations of protein dynamics were carried out at temperatures even higher than what we consider in our study. At low temperatures the peptides fold onto either α- helices or hair-pins. If they fully fold onto α-helices, the number of hydrogen bonds in them will be equal to the number of amino acids in a single peptide minus four (hence, 6 for A 10 and 26 for A 30 ). If, however, they fold onto hair-pins, the number of hydrogen bonds will depend on the degree of symmetry of the molecular structure. In the simulation, the effect of the solvent is taken into account implicitly through adjusting the value of 𝘀 HP for hydrophobic interaction. Here we set it to 𝘀 HP = 0.1𝘀 HB . 26 Four systems that contained 8, 16, 32, and 48 peptides of A 10 are simulated. Also simulated is a system with 16 A 30 . The sizes of simulation boxes for the first four systems are, respectively, 11.00 3 , 13.87 3 , 17.47 3 , and 19.97 3 nm 3 , while the system’s size for the case of A 30 is 19.97 3 nm 3 , all of which correspond to a same protein density. Initially, the proteins are inserted at randomly-selected positions within the simulation box. Their initial velocities are selected from the Maxwell-Boltzmann distribution at 𝑇 = 0.083 (250 K). The system is then thermalized for 100 ps, after which it is heated up to 𝑇 = 0.16 (500 K) and thermalized again for 200 ps. After that, the system is quenched by 50 K/100 ps until the temperature drops to 𝑇 = 1.1 𝑇 𝑓 (360 K for A 10 and 420 K for A 30 ), under which temperature the systems are run for up to 90 ns. During the simulation, if two proteins form a hydrogen bond or have hydrophobic interaction with each other for at least 500 simulation steps, we consider them as an aggregate or cluster, in which case they are moved together. Five independent series of simulations are performed for each case, in order to assess the effect of the initial conditions. Thus, the results represent averages of the five MD runs. 3.3 RESULTS We first determine the folding temperature 𝑇 𝑓 of the proteins by the standard method of computing the temperature-dependence of the average potential energy of a single protein. Figure 3-1 presents the result for A 10 indicating that its folding temperature is about 𝑇 𝑓 ≈ 325 K. For A 30 we obtain 𝑇 𝑓 ≈ 380 K. We then carry out long DMD simulations to study the dynamics of protein aggregation. Several important properties of the phenomenon are computed, most of which were not computed and reported in the previous molecular simulations of protein aggregation. In what follows we present the results and describe their implications. 27 Figure 3-1 Determination of the folding temperature of A 10 from a plot of the potential energy U versus dimensionless temperature T. Arrows indicate the locations of the folding temperature and the corresponding potential energy. 3.3.1 FORMATION OF β-SHEETS AND FIBRILS An important aspect of protein aggregation is the formation of β-sheets, sometimes referred to as β- pleated sheets. The sheets consist of β-strands that are connected laterally, where a β-strand is a stretch of polypeptide chain that typically contains 3-10 amino acids and is in an extended conformation in the backbone. To check whether a β-sheet has actually formed, we check the number of hydrogen bonds formed between the backbones of two β-strands. If the number is larger than half the number of amino acids of the participated peptide, we consider they have formed a β-sheet. β-sheets are, generally speaking, twisted and pleated. When large β-sheets overlay each other, they form fibrils that are believed to play an important role in the development of many neurodegenerative diseases, such as Alzheimer’s. Experimental observations 123-125 indicate that the formation of fibrils is preceded at the very early stage of aggregation by the emergence of small aggregates or micelles that are more or less spherical. When they join, they form β-sheets or proto-fibrils. The micelles represent the initial assembly stage and essentially set the path to the eventual formation of “mature” fibrils. PRIME model that we use in this 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 T -7 -6 -5 -4 -3 -2 -1 0 U 28 work does produce such spherical or almost spherical small particles at a very early stage of the aggregation process. Figure 3-2 presents examples of such particles that are obtained by simulating aggregation of 16 A 30 . They are produced after about 8 ns of simulations. At the similar early stage of the aggregation, our simulations also produce small disordered aggregates that are elongated. The reason could be the structure of our model, which has polar amino acids on both terminals of the proteins. In a less coarse 126 or an all-atom model, the proteins can have polar and non-polar terminals, which help producing more symmetrical spherical micelles. Figure 3-2 Examples of spherical and nearly spherical particles that are formed at a very early stage of the aggregation. Figure 3-3 presents four stages of the formation of β-sheets and fibrils. Initially, all the peptides are in a random coil configuration. After some time, small aggregates/spherical micelles are formed, shown in Figure 3-3(b). As the aggregation process continues, the conformation of the spherical aggregates changes as they convert into β-sheets, as shown in Figure 3-3(c), while they grow simultaneously with the monomers surrounding them. Finally, the β-sheets overlay each other and form fibrils, a snapshot of which is shown in Figure 3-3(d). 29 Figure 3-3 Four stages of the aggregation process of 16 A 30 peptides. (a) The initial configuration. All the peptides were random coils. (b) Small spherical micelles are formed. (c) Large β-sheets formed. (d) β-sheets overlay each other and formed fibrils. The typical structure of the β-sheets is shown in Figure 3-4. Hydrogen bonds have been formed between the N beads and C beads on the backbones, contributing to the formation of β-sheets. We also show in Figure 3-4 the details of a fibril structure which formed according to PRIME model. The non-terminal R beads are hydrophobic based on the simulation settings. Such beads attract with each other if they are 30 close. After forming the β-sheets, the R beads are on the surface, hence providing the opportunity for the β-sheets to overlay and form fibrils. Figure 3-4 The structures of the β-sheet (left) and fibrils (right). An effective way of visualizing energetically allowed regions of the backbone dihedral angles, 𝜓 and 𝜙 (see Figure 3-5), of amino acid residues in proteins is the Ramachandran plot. 127 In a polypeptide, the N– C α and C α –C bonds in backbone are relatively free to rotate. The rotations are represented, respectively, by the torsion angles 𝜓 and 𝜙 . Figure 3-5 presents a diagram for the dihedral angles as well as a Ramachandran plot for the alanine model in PRIME calculated by simulating a system containing 16 A 10 . The plot is constructed based on the normalized average potential energy <U> (including intra- and inter- hydrophobic forces and hydrogen bonding). The energy map is divided into a 72 x 72 matrix, with the value in each 5° bin being the average of sum potential energy of all the amino acids that share the same range of the (𝜓 , 𝜙 ) values during the entire simulation. Note that because we use DMD simulations, the potential energy is discontinuous. Thus, the only places in which there are changes in the potential energy are those where hydrophobic and hydrogen bonding interactions occur. The large red spot in the top left corner of the plot represents the β-sheet conformations, while the one to the left of the plot’s center is due to the α-helices. 31 Figure 3-5 The Ramachandran plot in terms of <U>. The large red spot at the top left corner is due to the β-sheets, while the one in the middle left corresponds to α-helices structure. 𝜑 and 𝜓 are the dihedral angles (see the top). The α-helix area of our plot agrees almost exactly with both the experiment 128-129 and quantum- mechanical (QM) calculations. 130-131 For the β-sheet area of the plot, our results are less precise than the experimental data and the QM computations: ours is larger and has only one “hotspot” - the one at the top left corner in Figure 3-5 - whereas the QM computations and the experimental plots have two smaller 32 spots in the same region. An all-atom MD simulation 132 using AMBER force field also produced the similar results and our plots match them very well. The reason that our computed Ramachandran plot is not as accurate is the mesoscale nature of PRIME. Recall that PRIME combines the H and N atoms into N beads, α-C and α-H atoms into C α beads, and the β-C and O atoms into C beads. Thus, the dihedral angles between the backbone beads are a little different from those in an all-atom model. But despite its relative simplicity, DMD simulation and PRIME are capable of reproducing important qualitative aspects of proteins and, therefore, are reliable enough for molecular simulation of their aggregation. Figure 3-6 Dynamic evolution of the total number of aggregates N(t) formed. n is the number of peptides initially in the system. 3.3.2 NUMBER AND SIZE OF THE AGGREGATES Figure 3-6 presents the evolution of the number of aggregates in the four systems that initially contained 8, 16, 32, and 48 A 10 (the results for 16 A 30 are also included in this and the following figures, but they are discussed separately below). It is clear, as mentioned in the introduction, that the aggregation process is very rapid. Note that the inter-peptide interactions are of hydrophobic and hydrogen bond types. In the DMD simulations the hydrophobic forces are represented by a square-well potential with a 0 10 20 30 40 50 60 70 80 90 t (ns) 0 5 10 15 20 25 30 35 40 45 50 N(t) 8 16 32 48 16(A 30 ) n 33 large size, ranging from 0.441 to 0.661 nm, whereas hydrogen bonds can form only at a distance 0.42 nm. On the other hand, the potential well depth of a hydrogen bond is large - 6 kcal in our simulations - ten times larger 30, 103 than that of the hydrophobic force. As the peptides diffuse in the system, the hydrophobic R beads (the side-chain residues) attract with other similar beads. At temperatures around 𝑇 = 1.1𝑇 𝑓 , the attractive forces dominate the behavior of the system, keeping the peptides close to one another and giving rise to larger probability of the N and C beads forming hydrogen bonds. At such temperatures the hydrogen bonds are stable and difficult to break. Figure 3-6 suggests that all the numerical results may be collapsed onto a universal curve. Here we have three variables, 𝑁 (𝑡 ), the total number of aggregates at time 𝑡 , 𝑛 , the number of initial peptides in the solution, and the time 𝑡 . Thus, we introduce two scaled variables, 𝑋 = 𝑛 /𝑡 𝑧 ′ and 𝑌 = 𝑡 /𝑁 𝜔 . We then determine the most accurate values of 𝜔 and 𝑧 ′ that collapse all the data in Figure 3-6 onto a single curve. We find that 𝜔 ≈ 0.60 ± 0.03, 𝑧 ′ ≈ 0.77 ± 0.04 (Eqn 3-1) Figure 3-7 presents the data collapse. Thus, one has 𝑁 −𝜔 = 𝑡 −1 𝑓 (𝑛 /𝑡 𝑧 ′ ) (Eqn 3-2) where 𝑓 (𝑥 ) is the scaling function shown in the figure. Other combinations of the three variables were also attempted, but none resulted in the collapse of the numerical data. Eqn 3-2 is similar to the dynamic scaling for aggregation of colloidal clusters, 133-134 studied extensively in the past. At long times, 𝑁 (𝑡 ) is expected to decay as the aggregates join, while the density of the system is constant. If the system contains a large enough number of initial peptides, then at long times we expect, 𝑁 (𝑡 )~𝑡 −𝜉 (Eqn 3-3) 34 Figure 3-7 Rescaled plot of the data in Figure 3-6 with X=n/t 0.77 and Y=t/N 0.60 . Figure 3-8 Dynamic evolution of the total number of aggregates N(t) formed for the 48-peptide system. The slope of the straight line is about 0.84. Figure 3-8 presents 𝑁 (𝑡 ) for the case in which we initially had 48 peptides in the system. It shows that at long times 𝑁 (𝑡 ) roughly follows the power law Eqn 3-3 with 𝜉 ≈ 0.84 ± 0.05. We shall return to this point shortly. 0.00 0.50 1.00 1.50 2.00 Log (t) 0.8 1.0 1.2 1.4 1.6 Log [N(t)] 35 We now define a mean aggregate size 𝑆 (𝑡 ) by the standard formula used in percolation 135-136 and aggregation of clusters, 133-134 𝑆 (𝑡 ) = ∑𝑠 2 𝑛 𝑠 (𝑡 ) 𝑠 ∑𝑠 𝑛 𝑠 (𝑡 ) 𝑠 (Eqn 3-4) where 𝑛 𝑠 is the number of aggregates (including single peptides) that contain 𝑠 beads. Figure 3-9 presents the results. At long times, 𝑆 (𝑡 ) grows with time 𝑡 as a power law, 𝑆 (𝑡 )~𝑡 𝑧 (Eqn 3-5) For small numbers of initial peptides, the exponent 𝑧 depends on their number. For large 𝑛 , however, namely, 𝑛 = 32 and 48, 𝑧 approaches a constant 𝑧 ′ ≈ 0.77 ± 0.04. This is shown in Figure 3-10 that presents the plot of 𝑆 (𝑡 ) versus 𝑡 for 𝑛 = 48. According to the dynamic scaling for aggregation of clusters, the exponent 𝜉 , 𝑧 and 𝑧 ′ are in fact identical, and the closeness of our estimates of the three exponents confirms this. In the scaling theory of aggregation of clusters 133-134 , 𝑧 is referred to as the dynamic exponent. The similarity between the aggregation of proteins and that of random clusters, which the latter is a well-studied problem in colloid and materials science, is important and, to our knowledge, is pointed for the first time in our study. We will return to this shortly. 36 Figure 3-9 Dynamic evolution of the mean aggregate size S(t). n is the number of the initial peptides. Figure 3-10 Dynamic evolution of the mean aggregate size S(t) in the 48-peptide system. The line representing fit of S(t) at long times has a slope of about 0.77. 0 10 20 30 40 50 60 70 80 90 t (ns) 0 100 200 300 400 500 600 700 S(t) 8 16 32 48 n 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Log (t) 1.50 1.75 2.00 2.25 2.50 2.75 3.00 Log [S(t)] 37 In the study of aggregation of colloidal clusters in three-dimensional systems, the dynamic exponent 𝑧 was found to depend on the diffusivity 𝐷 of the aggregates. In particular, it is found that if one assumes 137 , 𝐷 ~𝐷 0 𝑀 𝛾 (Eqn 3-6) where 𝑀 is the mass of a cluster and 𝐷 0 is a constant, then, the dynamic exponent 𝑧 depends on 𝛾 . For the particular case of 𝛾 = −1, i.e., when the diffusivity of a cluster (or an aggregate) is inversely proportional to its size, Meakin et al. 137 found that 𝑧 ≈ 0.80 ± 0.05. This is consistent with our estimate of 𝑧 . As we will show below, in our DMD simulations the effective diffusivity is indeed roughly proportional to the inverse of aggregates’ size, in agreement with what Meakin et al. 137 had reported, even though their simulations had been carried out by a Monte Carlo method and no molecular interactions had been taken into account. The apparent similarity between the aggregation of colloidal clusters and that of proteins deserves further discussion. One may argue that the interactions between colloids of single particles are isotropic, resulting in the clusters being able to grow in different directions during the formation. On the other hand, the interactions of short alanine-based peptides, which are chained particles, are non-isotropic. Therefore, such peptides form β-sheets that eventually aggregate together into ordered fibril structures. The difference is accentuated when the size of the simulated systems grows, hence allowing the formation of an extended and asymmetric fibril structure. In this case, each β-sheet is infinitely long. Would it be that we get the similarity only due to our simulated systems are relatively small, containing just up to 48 short peptides, so each fibril structure is relatively spherical and similar to a colloidal cluster? To answer this question, the size effect needs to be addressed. As mentioned earlier, we carried out simulations with 16 A 30 , the peptide three times larger than A 10 . Their folding temperature rises, of course. Figure 3-11 presents the very large final aggregate of A 30 at 360 K. Since the temperature is below their 𝑇 𝑓 , the proteins preserved their α-helix structure. So, this might refer to as a “native” aggregation. In this case, aggregation is due to hydrophobic forces, rather than hydrogen bonds, and is fast, resulting quickly in a stabilized aggregate. The structure looks like a twisted rope. 38 Figure 3-11 The large aggregate formed by 16 A 30 at 360 K. Different colors represent different proteins. Once we raise the temperature to 420 K, other aggregate structures emerge and are shown in Figure 3- 12. The configurations include fibrils that consist of two crossing layers of β-sheets, with each sheet consisting of stretched peptides, shown in Figure 3-12(a). One also has fibrils that consist of two layers of β-sheets that are parallel with each other, with some of the peptides forming intra-peptide β-sheets or hairpins, which are shown in Figure 3-12(b). Finally, one also has bended β-sheets because the stretched peptides are very long and flexible, thus, bend in the middle. Such aggregation process is slower than that at 360 K, since in this case it is hydrogen bonding that contributes the aggregation. It will take longer time for the peptides to find the optimal positions in order to form stable inter-molecular hydrogen bonds. 39 Figure 3-12 Same system as of Figure 3-11, but at 420 K. (a) - (c) indicate the various types of aggregate structures. 40 In fact, despite significant differences between the size of A 10 and A 30 as well as their aggregate structures, we find that the latter also follow the dynamic scaling that we presented above for the shorter peptides. In Figure 3-6 and Figure 3-7, as well as Figure 3-13 and Figure 3-14 below, we have included the results for the A 30 proteins. They all follow the dynamic scaling and the collapse of the data. We conclude, therefore, that the dynamic scaling is a general feature of our model of protein aggregation. On the other hand, our proposed dynamic scaling was developed based on an intermediate-resolution molecular model of proteins. It remains to be confirmed by precise measurements for actual proteins. The significance of the apparent similarity between the aggregation of colloidal clusters and that of proteins is at least threefold. One is that one may use the vast knowledge that already exists for the aggregation of colloidal clusters to gain deeper understanding of protein aggregation. The second is that one may be able to utilize the theoretical approaches to colloidal aggregation to study protein aggregation (see also the discussions in Section 4 below). Finally, despite the complexities of molecular structures of various proteins, as well as the interactions between them, the existence of the dynamic scaling indicates that certain features of protein aggregation are universal, independent of the molecular structures and, therefore, simpler models of proteins and their aggregation can also shed light on this important set of phenomena. 3.3.3 NUMBER OF β-SHEET FORMING PEPTIDES As already discussed, the formation of β-sheets is an important aspect of protein aggregation. Thus, it is important to study the number of peptides that participate in the formation of the sheets. Indeed, as time passes, large numbers of peptides participate in the formation process. This is demonstrated in Figure 3-13, where we plot the time-dependence of 𝑁 𝛽 (𝑡 ), the number of peptides participating in the formation of β-sheets at time 𝑡 , for different numbers of the initial peptides in the system. 𝑁 𝛽 (𝑡 ) increases with 𝑡 , but because the total number of peptides in the system is constant, it eventually approaches a constant value. 41 Figure 3-13 Dynamic evolution of N β (t), the number of peptides that are involved in the β-sheet formation. n is the number of initial peptides. Similar to 𝑁 (𝑡 ), the total number of aggregates at time 𝑡 , 𝑁 β (𝑡 ) also appears to follow a universal dynamic scaling. Indeed, if we plot 𝑡 𝑁 β versus 𝑛 𝑡 𝑧 , the results in Figure 3-14 collapse almost completely onto a single curve, which is nearly a straight line, implying that, 𝑁 β ~𝑡 −1 ℎ(𝑛 𝑡 𝑧 ) (Eqn 3-7) where ℎ(𝑥 ) is another scaling function, which is roughly linear. 0 10 20 30 40 50 60 70 80 90 t (ns) 0 5 10 15 20 25 30 35 40 45 50 N b (t) 8 16 32 48 16(A 30 ) n 42 Figure 3-14 Rescaled plot of the data of Figure 3-13 with X=nt 0.75 and Y=tN β . 3.3.4 NUMBER OF HYDROGEN BONDS The formation of hydrogen bonds plays an important role in the aggregation of proteins and, thus, it is important to study not only the number of such bonds that are formed over time, but also their types, namely α, β, and any other types. We define a hydrogen bond as an α-helix type by checking whether the two beads forming the bond are in the same peptide and their sequence interval is 4. A hydrogen bond is considered as a β-sheet type is based on the values of the dihedral angles, 𝜙 and 𝜓 , of the corresponding beads. If 160° < 𝜙 < 60° and, at the same time, 90° < 𝜓 < 180°, the hydrogen bond is of β- sheet type. Other types of hydrogen bonds include β-turn, 3(10)- and π-helix hydrogen bonds. As mentioned earlier, we have used the configuration (PHHHHHHHHP) for A 10 , where P represents a polar side chain and H represent a hydrophobic one. Such configurations are more amenable to interacting with each other, if they are even relatively close. Therefore, the amino acids in the middle parts of the peptides have higher opportunities to form inter-peptide hydrogen bonds. They also form more stable hydrogen bonds comparing with the N and C terminals, because the former has two neighboring beads to help stabilizing the newly formed hydrogen bonds while the latter only has one. 43 Figure 3-15 Dynamic evolution of the numbers of various types of hydrogen bonds formed in the systems initially containing (a) 8 (b) 16 (c) 32 and (d) 48 A 10 , and (e) 16 A 30 . 0 10 20 30 40 50 60 70 80 90 t (ns) 0 10 20 30 40 50 60 70 80 90 100 110 120 N HB a b Other Total type (b) 0 10 20 30 40 50 60 70 80 90 t (ns) 0 20 40 60 80 100 120 140 160 180 200 220 N HB a b Other Total type (c) 0 10 20 30 40 50 60 70 80 90 t (ns) 0 50 100 150 200 250 300 350 N HB a b Other Total type (d) 0 10 20 30 40 50 60 70 80 90 t (ns) 0 5 10 15 20 25 30 35 40 45 50 55 N HB a b Other Total type (a) 0 10 20 30 40 50 60 70 80 90 100 t (ns) 0 50 100 150 200 250 300 N HB α α Other Total type (e) 44 Figure 3-15 presents the results. As long as a peptide is isolated, the α-helix hydrogen bonds are formed spontaneously, and their number increases. But as the peptides diffuse in the solution and come close to each other, aggregates and β-sheets form and thus, the number of β type as well as other types of hydrogen bonds increases, accompanying by the decrease of the number of α-helix hydrogen bonds. 3.3.5 RADIUS OF GYRATION OF THE AGGREGATES Another interesting property of the aggregates is their center-of-mass radius of gyration, 𝑅 𝑔 , and its evolution with time. The increase with time of 𝑅 𝑔 is, of course, a measure of the growth of protein aggregates. Figure 3-16 presents 𝑅 𝑔 of aggregates of various sizes in the system with 48 A 10 . For a single A 10 with an α-helix structure, 𝑅 𝑔 ≈ 0.7 nm, whereas for one with a stretched coil structure, 𝑅 𝑔 ≈ 0.95 nm. Figure 3-16 Dynamic evolution of radius of gyration R g of aggregates of various sizes in the 48-peptide system. 0 10 20 30 40 50 60 70 80 90 t (ns) 1 2 3 4 5 6 7 R g 2 (nm 2 ) size=02 0 10 20 30 40 50 60 70 80 90 t (ns) 1 2 3 4 5 6 7 R g 2 (nm 2 ) size=04 0 10 20 30 40 50 60 70 80 90 t (ns) 1 2 3 4 5 6 7 R g 2 (nm 2 ) size=08 0 10 20 30 40 50 60 70 80 90 t (ns) 1 2 3 4 5 6 7 R g 2 (nm 2 ) size=10 0 10 20 30 40 50 60 70 80 90 t (ns) 1 2 3 4 5 6 7 R g 2 (nm 2 ) size=18 0 10 20 30 40 50 60 70 80 90 t (ns) 1 2 3 4 5 6 7 R g 2 (nm 2 ) size=26 45 The fact that for the small aggregates, namely, those made of 2 and 4 peptides, 𝑅 𝑔 remains constant means that some small aggregates remain isolated from the rest and survive for a long time. The largest fluctuations in 𝑅 𝑔 belong to the intermediate-size aggregates, which are keeping breaking up and then re-aggregating again. Moreover, the fact that 𝑅 𝑔 of the largest aggregate remains constant implies that it is stable and does not break up, at least over the period of time that the DMD simulations are carried out. In addition, the fluctuating 𝑅 𝑔 for the large aggregates may represent the conversion of single-layer β-sheets into multiple-layer fibrils as well as the reversion. 3.3.6 DIFFUSIVITY Figure 3-17 Dependence of the aggregates’ diffusivity D on their mass M. n is the number of initial peptides. Although the effect of the solvent is not explicitly taken into account, we also estimated the “diffusivity” 𝐷 of the aggregates of various sizes by calculating the mean-square displacements 〈𝑅 2 (𝑡 )〉 of the center of mass, and using the relation 〈𝑅 2 (𝑡 )〉 = 6𝐷𝑡 at long times when Fickian diffusion is reached. The results are presented in Figure 3-17. They follow the expected trends, namely, that the larger an aggregate the slower it diffuses. More importantly, the logarithmic plot of the diffusivity 𝐷 versus the 02 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 aggregate mass M 0.5 1.0 1.5 2.0 2.5 3.0 D 8 16 32 48 n 46 aggregate mass 𝑀 shown in Figure 3-18 indicates that 𝐷 follows Eqn 3-6 with 𝛾 ≈ −1, completely consistent with the expected value of the dynamic exponent 𝑧 , when the diffusivity is inversely proportional to the size of a cluster. Figure 3-18 Logarithmic plot of the dependence of the diffusivity D on the aggregates mass M. Solid line corresponds to the 32-peptide system and dashed line corresponds to the 48-peptide system. 3.4 DISCUSSION: THE CATALYTIC EFFECTS As described in the introduction, neurodegenerative diseases, such as Alzheimer’s disease (AD), have been linked to misfolding of certain types of proteins/peptides that aggregate together and form amyloid fibrils with a common β-sheet structure, and our simulations do reproduce such molecular structure. For example, a most important molecular structure associated with the AD is 138 the 42-residue β-amyloid, or Aβ42. Several studies over the past decade have indicated 139-142 that such molecule as well as other prefibrillar oligomeric aggregates of amyloidogenic peptides is the major toxic molecule that causes those infamous neurodegenerative diseases. So, the question is, what is the molecular mechanism that gives rise to the formation of such toxic aggregates? 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 Log (M) -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 Log (D) 32 48 n 47 It has been demonstrated that 143 the surface of high molecular-weight fibrillar aggregates, as well as monomers, acts as very effective catalysts that facilitate the formation of Aβ42 aggregates and similar toxic materials that contribute to the AD and similar diseases. In this mechanism, both fibrils and monomers participate in a catalytic nucleation reaction 144-145 that significantly increases 143, 146-148 the rate of aggregation of Aβ42 and similar toxic molecules. Our DMD simulations do not include the catalytic effect. In order to carry out molecular simulations that allow reactions, one would need a reactive force field 149 that allows the catalytic reaction to occur. To our knowledge, such a reactive force field does not currently exist. Even if it does exist, one would need to carry out an all-atom simulation of the aggregation process and the sideway catalytic reaction, which is, however, way beyond the current computational power. The question, then, is whether the dynamic scaling that we have discovered for the aggregation of peptides and formation of fibrils will still hold for the aggregation of Aβ42 and similar fibrils, if the catalytic effects are included? While the question cannot be surely addressed by molecular simulation, one may develop other theoretical and computational approaches to study the problem. Since our simulations indicate that the aggregation of peptides has qualitative similarities with the coagulation of colloidal aggregates, one possible approach is through the Smoluchowski equation 150 by allowing a reaction term in the equation. Another possible approach could be based on the nonlinear master equation. 145 At the same time, however, we recall that the dynamic scaling for the aggregation of clusters is very robust. Thus, we conjecture that the general form of the dynamic scaling may survive after allowing the catalytic reactions to occur, although the numerical values of the exponents 𝑧 and 𝜉 may change. This remains to be seen, however, and we hope to report the results of such a study in the near future. 3.5 SUMMARY In this chapter, PRIME model is utilized together with extensive DMD simulation to study the aggregation of proteins in an unbounded fluid. Various properties of the aggregates as well as the process are computed, including dynamic evolution of the total number of aggregates, mean aggregate size, number of peptides that contribute to the formation of β-sheets, number of various types of hydrogen bonds, 48 radius of gyration of the aggregates, and the aggregates’ diffusivity. We show that many of such quantities follow dynamic scaling, similar to those for the aggregation of colloidal clusters. In particular, at long times the mean aggregate size 𝑆 (𝑡 ) grows with time as 𝑆 (𝑡 )~𝑡 𝑧 , where 𝑧 is the dynamic exponent. To our knowledge, this is the first time that the qualitative similarity between the aggregations of proteins and colloidal clusters has been pointed out. Given that the dynamic scaling is obtained with an intermediate-resolution model of proteins, it remains to be seen whether the scaling can be confirmed by experiments. In the next chapter, we will report the aggregation of proteins in confined media, as models of the “crowded” cellular environments. We will also discuss the effects from the confinements on the dynamics of protein aggregation. We will show that there are significant differences between protein aggregation in a confined medium and in an unbounded fluid that we present in this chapter. ACKNOWLEDGMENTS We thank an anonymous referee whose constructive comments and criticisms helped us improve the quality of the results and provided the impetus for carrying out the simulations with much larger proteins. 49 4. DYNAMICS OF PROTEIN AGGREGATION. II. DYNAMIC SCALING IN CONFINED MEDIA 4.1 INTRODUCTION In Chapter 3 110 , we have discussed the protein aggregation phenomenon in unbounded domains (bulk conditions). We demonstrated that many of the reported properties follow a universal dynamic scaling (see below). In particular, we showed that at long times the mean aggregate size 𝑆 (𝑡 ) grows with time as, 𝑆 (𝑡 )~𝑡 𝑧 (Eqn 4-1) where 𝑧 is called the dynamic exponent. This proposed dynamic scaling has been confirmed recently by experiments on insulin aggregation at low pH (about 1.6) 151 . In that study, the structure of the prefibrillar intermediates is characterized through a combination of dynamic light scattering, small-angle neutron scattering, atomic force microscopy, and circular dichroism. However, it is well-known and well-understood that the properties of many phenomena that occur in a confined medium, in which the folding and aggregation of proteins occur in vivo, are vastly different from those happen in an unbounded medium or bulk. Thus, studying protein aggregation in unbounded fluids only does not suffice for understanding of the phenomenon in biological environments. Now the question is, what are the major differences, as well as similarities, between aggregations of proteins in confined media and bulk conditions? In particular, we wish to understand whether the dynamic scaling for protein aggregation that discussed in Chapter 3 still holds in confined media, or whether it should be modified, or it does not hold at all. In this chapter, we will report the results of extensive molecular simulations of protein aggregation in confined media as models of the “crowded” cellular environments 111 . One important goal of the present study is to gain an understanding of the confinement effect on the dynamics of aggregation. We show that there are significant differences between protein aggregations in a confined medium and in an unbounded fluid. 50 The structure of this chapter is as follows. In Section 2 we describe the models of peptides and confined media as well as the simulation configurations we use in DMD. The results are presented in Section 3. Section 4 gives a discussion of the root cause of the dynamic scaling and the associated power laws in protein aggregation, while Section 5 describes the differences between protein aggregations in confined and unconfined media. The last section summarizes the chapter. 4.2 CONFIGURATION OF DMD SIMULATION 4.2.1 MODELS OF PEPTIDES Similar to our work reported in Chapter 3, 110 as well as our earlier DMD simulations that studied the effect of confinement on the stability and folding of proteins 103, 114 and on their diffusion, 115 here we model and simulate de novo-designed α-family of proteins. 58, 61, 116-117, 152 We utilize a mesoscopic model of proteins with a ten-residue sequence and one type of amino acid. The amino acid is to mimic alanine (A) and each consists of four beads. Its sequence in the model protein therefore is (AAAAAAAAAA), or A 10 . The β-amyloid, or Aβ, which is generally believed to be an important part of neurological plaques in Alzheimer’s disease, typically contains 39-42 amino acids. Therefore, we also report the results of DMD simulations of a much larger protein, A 30 , with a thirty-residue sequence of the same amino acid, which is much closer to the current known experimental information on the structure of Aβ. More limited simulations are also carried out with the A 15 , A 20 and A 25 proteins, as described below. The amino acids in the sequence may be hydrophobic (H) or polar (P) depending on their locations, with the polarity sequence of A 10 being PHHHHHHHHP, or (PH 8 P), and that of A 30 being (PH 28 P), and similarly for the other peptides. The model is intended to mimic the hydrophobic part of Aβ. The mesoscopic model of proteins that we utilize is PRIME, 30, 81, 84-85, 88, 103, 110, 112, 114-115, 118-120 which has been described in Chapter 2 of this dissertation. 51 4.2.2 MODELS OF CONFINED MEDIA Two models of confined media are utilized. One is a spherical cavity of diameter 𝐷 . Four of such cavities are used to simulate the aggregation of 8, 16, 32, and 48 A 10 peptides with the corresponding diameters of the cavities being, respectively, 8.0, 10.0, 12.7 and 14.5 nm. The cavities’ diameters are varied in order to keep the mass density of proteins the same in every case, and equal to that in the bulk, so that comparisons with the bulk conditions are meaningful. A fifth cavity of diameter 14.5 nm is used to simulate the aggregation of 16 A 30 peptides. The second model of the confined media is a circular tube of diameter 𝐷 and length 𝐿 whose ends are closed. Two of such tubes are used in the simulations. One has an aspect ratio 𝐿 /𝐷 = 2. The tube’s diameters for simulating 8, 16, 32, and 48 A 10 peptides are, respectively, 5.6, 7.0, 8.8, and 10.0 nm. A fifth tube of the same type with 𝐷 = 10 nm is used to simulate the aggregation of 16 A 30 peptides. We refer to this system as Tube2. The second tube, referred to as Tube3 and used for simulating the aggregation of the same numbers of A 10 peptides, has an aspect ratio of 𝐿 /𝐷 = 3 with 𝐷 = 4.8, 6.0, 7.7, and 8.8 nm. A fifth tube3 with 𝐷 = 8.8 nm is also used to simulate the aggregation of 16 A 30 peptides. We also carry out simulations of aggregation for the same numbers and types of peptides under bulk conditions and adjust the dimensions of the cubic simulation boxes so that for each case simulated in a confined medium, we can simulate the same in the bulk with the same mass density of proteins. It is assumed in all cases that the walls of the confined media are made of carbon, due to the ubiquity of carbon in biological cells. 4.2.3 MOLECULAR SIMULATION As already pointed out, it is currently not feasible to carry out all-atom MD simulation of protein aggregation in confined media, particularly if we are to simulate the phenomenon over a relatively long period of time for a large system. The same is true about a mesoscopic model, such as PRIME that we utilize, if complex force fields are used. Thus, we employ DMD simulation. 38, 40, 153-156 The following will describe the details of the simulation configurations in the present work. 52 Since the interaction potentials are simplified while still being kept realistic, we are able to carry out simulations on the order of several microseconds in real times. In particular, the interactions between all the beads, except for hydrogen bonding (HB) or hydrophobic attraction, are hard-core repulsions. The potential for HB and hydrophobic interactions is of square-well type. The N and C beads on the same peptide separated by at least three intervening residues, or located on different peptides, may have HB interactions with a strength of 𝘀 HB and at the range of about 4.2 Å, under the premise of no bead has contributed to more than one HB. An α-HB is one between C i and N i+4 , on the same peptide. We refer to other HBs as non-α. Temperature is expressed in units of a temperature scale 𝑇 𝑠 such that, 𝒩 𝑘 𝐵 𝑇 𝑠 = 𝘀 HB , with 𝑘 𝐵 being the Boltzmann’s constant and 𝒩 the Avogadro’s number. Energy is measured in units of 𝘀 HB , which itself depends on the types of neighboring atoms, the structures of molecule, 121, 157 and whether water is explicitly presented in the system. In polyalanine one has, 121 3.5 kcal/mol < 𝘀 HB < 8.6 kcal/mol, whereas MD simulation 122 indicated that 𝘀 HB for α-helices is, 1.9 kcal/mol < 𝘀 HB < 5.6 kcal/mol. Based on such considerations, we assumed that 𝘀 HB = 6 kcal/mol (25 kJ/mol), implying that 𝑇 s = 3000 K 110 . The dimensionless temperature 𝑇 ∗ is then defined by, 𝑇 ∗ = 𝑇 /𝑇 𝑠 , but for the sake of convenience we delete the superscript ∗. Since proteins of various sizes would have different folding temperatures, the simulations for the proteins A p are carried out at 𝑇 = 1.1𝑇 𝑓 (𝑝 ) , where 𝑇 𝑓 (𝑝 ) is the folding temperature of A p with 𝑝 residues. The relatively-high selected temperatures keep the proteins in some sort of unfolded state during the simulations, hence accelerating the aggregation. Many of the previous DMD simulations of protein dynamics were carried out at temperatures even higher than what we consider here. The temperature is held constant by using Andersen thermostat 113 . Solvent is implicitly presented and its effect is taken into account by adjusting the value of energy parameter 𝘀 HP for hydrophobic interactions. In our simulations, we set it as 𝘀 HP = 0.1𝘀 HB . To commence the simulation, the proteins are initially randomly inserted into a simulation box or a confined medium with initial velocities that are selected from the Maxwell-Boltzmann distribution at low temperatures. We begin the simulation at 𝑇 = 0.083 (250 K) for 100 ps. This usually leads to all the peptides forming α-helices, fully or partially. All the bonds and angles are also relaxed to their equilibrium 53 states. We then increase the temperature to 𝑇 = 0.167 (500 K) and thermalize the system for 2 ns, which leads to the peptides being in random-coil state. The system is then cooled down to the target temperature, 𝑇 = 1.1𝑇 𝑓 (𝑝 ) , with steps of 30 K per 100 ps. Five independent series of simulations are carried out for each case, in order to assess the effect of the initial conditions. Thus, the results represent the averages of all five MD runs. In Chapter 3 we reported the aggregation of proteins under bulk conditions, but the results were for a low-density system. In the present work the mass density of proteins in the confined media is high. Thus, to make meaningful comparison with the bulk conditions, we also perform new simulations under bulk conditions in which the mass density of proteins is the same as in the confined media. 4.3 RESULTS The first step of the simulation is to determine the folding temperature 𝑇 𝑓 of the different sizes of proteins. To do so, we use the standard method based on computing the temperature-dependence of the average potential energy of a single protein. The results are 𝑇 𝑓 ≈ 0.108, 0.118, 0.121, 0.124, and 0.127 for, respectively, A 10 , A 15 , A 20 , A 25 , and A 30 peptides. In order to simulate the aggregations of all the peptides of various sizes at an equal distance from their folding temperatures, we set, as mentioned earlier, the target simulation temperature 𝑇 = 1.1 𝑇 𝑓 (𝑝 ) . In what follows we present the results and describe their implications. 54 Figure 4-1 Four stages of the aggregation of 16 A 30 peptides in a spherical cavity, whose diameter is 1.5 times of the peptide length when they completely stretch out. Initially, all the peptides are in random-coil configuration. At the last stage, all the peptides have joined to form a single large fibril. 4.3.1 DEPENDENCE OF THE PROTEIN SECONDARY STRUCTURE ON THE GEOMETRY OF CONFINED MEDIA Figure 4-1 presents four stages of the aggregation of 16 A 30 peptides in a spherical cavity of diameter 𝐷 = 14.5 ≈ 1.5ℓ nm, where ℓ is the length of the fully extended peptides. All the peptides are initially in random-coil configuration. After some time, the peptides begin aggregating into amorphous/ spherical micelles, in agreement with experimental observations under bulk conditions, 123-125 see Figure 4-1(6 ns). At longer times, shown in Figure 4-1(13 ns), the β-sheets begin growing, during which hydrogen bonds are formed between the N and C beads inter-molecularly. In the end, as Figure 4-1(50 ns) indicates, t = 0 ns 6 ns 13 ns 50 ns 55 large β-sheets overlay each other and eventually form fibrils. According to PRIME model of proteins, the formation of fibrils comes about because some of the side-chain beads are hydrophobic and attract with each other if they are close enough. After forming the β-sheets, the side-chain beads are on the surface, hence providing the opportunity for the sheets to form fibrils, which are believed to play an important role in the development of many diseases, such as Alzheimer’s. If, however, we study the same process in a smaller spherical cavity, the shapes of the aggregates could be strikingly different from those in Figure 4-1, manifesting better the effect of confinement. Shown in Figure 4-2 are the final configurations of two large aggregates in a cavity with 𝐷 /ℓ = 1. Figure 4-2(left) indicates that if the β-sheets are eventually formed, their structure would be twisted due to the restriction imposed by the curvature of the cavity’s wall. Figure 4-2(right) shows that it is possible that no β-sheets and, thus, no fibrils, form at all due to the severely restricted structure of a confined medium. These are striking different from what happens in the aggregation of proteins under bulk conditions. Figure 4-2 Comparison of the aggregates of 10 A 30 peptides in a small spherical cavity whose diameter is the same as the peptides’ length in their most stretched form. On the left, a β-sheet has formed, but its shape is twisted due to the wall curvature, whereas on the right, no sheet has formed, and one has only an amorphous aggregate. t = 50 ns t = 50 ns 56 Similar phenomena are seen in cylindrical tubes. Shown in Figure 4-3 are four stages of the aggregation of 16 A 30 peptides in Tube3. In this case the peptides aggregate into amorphous/spherical cluster(s), after which the β-sheets begin to grow. The time at which the growth of sheets begins varies. Similar to the spherical cavity, due to the wall effect and confinement, the stretched β-sheets are more likely to fold/roll onto themselves, an example of which is also shown at the bottom in Figure 4-3. Figure 4-3 Four stages of the aggregation of 16 A 30 peptides in Tube3 (cylindrical tube with aspect ratio L/D=3). Initially, all the peptides are in random-coil configuration. Eventually, two stable β-sheets have formed at 91 ns. Due to the confinement, the formed β-sheets may fold/roll, and finally form an aggregate whose configuration is shown at the bottom. Note, however, that the initial distribution of peptides in the tube is important. t = 0 ns 46 ns 12 ns 91 ns 57 Figure 4-4 Dynamic evolution of the total number of aggregates, N(t). n is the initial number of peptides. 4.3.2 DYNAMIC SCALING OF THE AGGREGATION PROCESS First, Figure 4-4 presents the dynamic evolution of total number of aggregates 𝑁 (𝑡 ) with various initial numbers of A 10 in the spherical cavity, Tube2 and Tube3, and compares them with the results in unbounded fluid with the same mass density of proteins. The results indicate that the aggregation is very rapid. As pointed out in Chapter 3, 110 the inter-peptide interactions consist of hydrophobic forces and hydrogen bonds. The former is represented in the simulation by a large square-well potential ranging from 0.44 to 0.66 nm. The latter, however, can form only at a distance of 0.42 nm, but their potential well depth is large - 6 kcal/mol in our simulations - which is ten times larger 103, 118 than that of the 0 10 20 30 40 50 t (ns) 0 5 10 15 20 25 30 35 40 45 50 N(t) 8 16 32 48 16(A 30 ) n Tube2 0 10 20 30 40 50 t (ns) 0 5 10 15 20 25 30 35 40 45 50 N(t) 8 16 32 48 16(A 30 ) n Tube3 0 10 20 30 40 50 t (ns) 0 5 10 15 20 25 30 35 40 45 50 N(t) 8 16 32 48 16(A 30 ) n Bulk 0 10 20 30 40 50 t (ns) 0 5 10 15 20 25 30 35 40 45 50 N(t) 8 16 32 48 16(A 30 ) n Sphere 58 hydrophobic force. When the proteins diffuse in the solution, the hydrophobic side-chains represented by R beads attract with the similar beads. At low enough temperatures the attractive forces are dominant, keeping the peptides close to one another and making it more likely that the N and C beads form HBs that are stable and difficult to break. Next, Figure 4-5 presents the comparison of 𝑁 (𝑡 ) for the aggregation of A 30 peptides in various geometries. Initially, the peptides aggregate and begin to form β-sheets in the spherical cavity faster than in Tube2, followed by the bulk condition, while the aggregation in Tube3 is the slowest. After some time, however, due to the relatively large aspect ratio of the tubes, the rate of aggregation in the bulk system exceeds those of the two tubes. At the end of the simulations, all the peptides in bulk, spherical cavity and Tube2 form a single large aggregate, whereas they form two smaller aggregates in Tube3. This is likely due to the fact that Tube3 is narrower (also 15% narrower than Tube2) under the same mass density of proteins, and its geometry gives negative effects on the aggregation. We will discuss this in detail below. Figure 4-5 Comparison of dynamic evolution of the total number aggregate N(t) of the A 30 peptides in various geometries. 0 20 40 60 80 100 t (ns) 0 2 4 6 8 10 12 14 16 N(t) Bulk Sphere Tube2 Tube3 59 Figure 4-6 The universal collapsed curve for the total numbers of aggregates N(t), with X=n/t z’ and Y=t/N ω . n is the initial number of aggregates in the system. As first proposed by us in Chapter 3, 110 the numerical data of Figure 4-4 may be collapsed onto a universal curve. We introduce two scaled variables, 𝑋 = 𝑛 /𝑡 𝑧 ′ and 𝑌 = 𝑡 /𝑁 𝜔 , and determine the most accurate values of 𝜔 and 𝑧 ′ that collapse all the data in Figure 4-4 onto a single curve. The results are shown in Figure 4-6 and indicate that the data collapse is good. More specifically, the excellent collapse of the data for shorter peptides indicates that, under the conditions studied, the aggregation process does follow universal scaling laws, and that the finite-size effects are insignificant. The small deviations of the results from the universal curve in the case of 16 A 30 peptides in the two tubes for intermediate values of the variable X are due to the large size of the peptides, but even in these two cases the data collapse is reasonable. It is likely that any deviations from the scaling laws for larger peptides arise from the 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 X 0 5 10 15 20 25 Y 8 16 32 48 16(A 30 ) n Tube2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 X 0 5 10 15 20 25 Y 8 16 32 48 16(A 30 ) n Tube3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 X 0 5 10 15 20 25 Y 8 16 32 48 16(A 30 ) n Bulk 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 X 0 5 10 15 20 25 Y 8 16 32 48 16(A 30 ) n Sphere 60 expected larger finite-size effect, rather than from failure to reach the scaling regime. This is supported by the fact that the proposed dynamic scaling has been verified in the recent experiments 151 . Thus, one may write, 𝑁 −𝜔 = 𝑡 −1 𝑓 (𝑛 /𝑡 𝑧 ′ ) (Eqn 4-2) where 𝑓 (𝑥 ) is the scaling function shown in Figure 4-6 with the property that 𝑓 (𝑥 ) → 0 for 𝑥 → 0 (i.e. at long times). The computed 𝜔 and 𝑧 ′ for all the cases are listed in Table 4-1. The exponents for the confined media are in agreement with those in an unbounded domain. Table 4-1 Summary of the exponents of dynamic scaling Exponent Bulk Sphere Tube2 Tube3 𝜔 0.30 ± 0.07 0.31 ± 0.08 0.34 ± 0.06 0.32 ± 0.07 𝑧 1.07 ± 0.03 0.98 ± 0.03 0.46 ± 0.02 0.61 ± 0.05 𝑧 ′ 1.05 ± 0.05 1.06 ± 0.05 1.06 ± 0.04 1.02 ± 0.05 𝘁 1.00 ± 0.03 1.07 ± 0.03 0.89 ± 0.05 1.04 ± 0.06 (0.42 ± 0.03) It is clear that as the aggregation proceeds, 𝑁 (𝑡 ) decays with time, while the total mass density of proteins remains a constant. Thus, given a large enough initial number of peptides, we expect at long times to have, 𝑁 (𝑡 )~𝑡 −𝜁 (Eqn 4-3) Figure 4-7 presents a plot of log 𝑁 (𝑡 ) versus log 𝑡 for the case with 48 peptides. The values of the exponent 𝘁 are listed in Table 4-1. Note that to within the estimated errors, the numerical values of 𝘁 are the same for the spherical cavity and bulk. In the case of the two tubes, however, the problem is a bit more complex. While 𝘁 for the aggregation in Tube2 is close to those in the spherical cavity and bulk, 61 it is not so in the case of Tube3. This might seem surprising because, intuitively, one might expect the opposite. However, considering that to hold the volume constant, compared with Tube3, Tube2 is stretched out without significantly increasing its lateral space. As a result, the peptides in Tube2 could diffuse better axially than Tube3 but at the same time need to move further to meet with each other compared with that in the spherical cavity and bulk, hence resulting in a slower rate of aggregation and rate of decrease in 𝑁 (𝑡 ), i.e. smaller 𝘁 . For Tube3, the geometry effect is more significant. Indeed, as Figure 4-7(Tube3) indicates, there is a crossover from the early time regime for which 𝘁 ≈ 1.04 ± 0.06, to the regime at longer times for which 𝘁 ≈ 0.42 ± 0.03. The exponent for the earlier time, when the peptides have still not “felt” the longer length of the tube, is in agreement with those in the other cases. Figure 4-7 Power-law scaling of the total number of aggregates N(t) for the aggregation of 48 A 10 peptides. The Tube3 geometry indicates two distinct scaling regimes. 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Log (t) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 Log [N(t)] Tube2 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Log(t) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 Log[N(t)] Tube3 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Log (t) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 Log [N(t)] Bulk 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Log (t) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 Log [N(t)] Sphere 62 From a physical point of view, the observed crossover in the narrow tubes is important, manifesting the strong effect of confinement on the aggregation process. Note that the crossover occurs at about log 𝑁 (𝑡 ) ≈ 0.6, i.e. 𝑁 (𝑡 ) ≈ 4. This indicates that when the total number of aggregates is about 4, under which situation the average size of the individual aggregate has been quite large, the aggregates are distributed along the tubes’ axis; thus, to decrease 𝑁 (𝑡 ) further, these large aggregates must join together, which, compared with that in a spherical cavity, is quite difficult in a narrow tube. In fact, the system no longer acts as a 3D structure, but rather close to a 1D medium. The dynamic scaling provides other predictions for the aggregation process. In particular, if 𝑛 𝑠 (𝑡 ) is the number of aggregates (including single peptides) at time 𝑡 that contain 𝑠 beads, then, as suggested in Chapter 3, 110 a mean aggregate size 𝑆 (𝑡 ) is defined in a manner similar to those used in the studies of percolation processes 135-136 and aggregation of colloidal clusters, 133, 158 namely, 𝑆 (𝑡 ) = ∑𝑠 2 𝑛 𝑠 (𝑡 ) 𝑠 ∑𝑠 𝑛 𝑠 (𝑡 ) 𝑠 (Eqn 4-4) At long times 𝑆 (𝑡 ) grows with time 𝑡 as a power law given by Eqn 4-1, making it possible to estimate the exponent 𝑧 . Figure 4-8 presents the results for the system with 48 A 10 peptides. The mean size of aggregates increases with time, and as Figure 4-8 indicates, there is clearly a scaling regime in which Eqn 4-1 is satisfied. The estimates for the exponent 𝑧 are all listed in Table 4-1. We also obtained similar results for the larger proteins A p , with 𝑝 = 15, 20, 25 and 30. According to the dynamic scaling, the exponents 𝘁 , 𝑧 and 𝑧 ′ should all be equal. In the case of 𝑧 ′ and 𝘁 , they are indeed equal for all cases to within their estimated errors. The exponent 𝑧 is also equal to 𝘁 and 𝑧 ′ for the bulk and spherical cavity cases, but not so for the two tubes. The reason is similar to the behavior of 𝑁 (𝑡 ) explained earlier: the large aspect ratio of the two tubes means that they are not really 3D, but quasi-1D. It is, indeed, known that the dynamic exponent z is more sensitive to the geometry of the system than the other exponents. 63 Figure 4-8 Power-law scaling of the mean aggregate size S(t) during aggregation of 48 A 10 peptides. The Tube3 geometry indicates two distinct scaling regimes. 4.3.3 NUMBER OF PEPTIDES THAT FORM β-SHEETS Of particular interest is time-dependence of 𝑁 β (𝑡 ), the number of peptides that participate in the formation of β-sheets. 𝑁 β (𝑡 ) increases with 𝑡 , but because the total number of peptides in the system is constant, it eventually approaches a constant value. For example, we show in Figure 4-9 the results for the aggregation of 16 A 30 peptides, which has the similar results for the smaller proteins. In the early stages of the aggregation, the peptides in the spherical cavity aggregate faster and begin to form β- sheets earlier than in Tube2 and bulk, with Tube3 being the slowest. Later on, due to the large aspect ratio of the tubes, the bulk system exceeds with the tubes. 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 Log(t) 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 Log[S(t)] Bulk 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 Log(t) 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 Log[S(t)] Sphere 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 Log(t) 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 Log[S(t)] Tube2 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 Log(t) 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 Log[S(t)] Tube3 64 Figure 4-9 Comparison of the dynamic evolution of the number of peptides N β involved in the formation of β-sheets during aggregation of 16 A 30 peptides in various environments. In this particular case, only the peptides in the bulk system eventually form a single large β-sheet/fibril. Confinement, as a model of the crowded cellular environment, plays a fundamental role in the formation of β-sheets. If confinement is severe, then the sheets that may eventually form will have a twisted structure because they must adapt themselves to the curvature of the confining medium’s walls. This is clearly demonstrated in Figure 4-2 and Figure 4-3. Another possibility, which we have observed in our simulations, is that no β-sheet is formed at all, at least not within any reasonable computational (and, perhaps, experimental) time. For example, in the case for which the results are shown in Figure 4-9, at the end of the simulations only the peptides in the unbounded domain form a single large β-sheet/fibril; others do not, which is due to the confinement. Similar to 𝑁 (𝑡 ), 𝑁 β (𝑡 ) also appears to follow dynamic scaling. Indeed, if we plot 𝑌 = 𝑡 𝑁 𝛽 𝜔 versus 𝑋 = 𝑛 𝑡 𝑧 , where 𝑧 is the same dynamic exponent as before, the data for the 48-A 10 system collapse onto a single curve with 𝜔 ≈ 0.95 ± 0.05. Moreover, the universal curve is nearly a straight line. This is shown in Figure 4-10, which implies that 𝑁 β ~𝑡 −1 ℎ(𝑛 𝑡 𝑧 ) (Eqn 4-5) 0 20 40 60 80 100 t (ns) 0 2 4 6 8 10 12 14 16 N b (t) Bulk Sphere Tube2 Tube3 65 where ℎ(𝑥 ) is another scaling function, which is roughly linear. Two important points should be mentioned here. One is that the results shown in Figure 4-10 represent only those cases in which the β- sheets are eventually formed. We discard the cases in which even after long times the sheets have not formed. The second point is that the collapse of the data for the largest peptides, A 30 , is delayed. The rescaled data for A 30 join the universal curve for larger values of 𝑋 = 𝑛 𝑡 𝑧 , i.e. longer times, which is again an indication of the effect of confinement on the aggregation of large peptides. Figure 4-10 The universal collapsed curve for N β (t), the number of peptides that form β-sheets, with Y=tN β ω versus X=nt z n is the initial number of aggregates in the system. Note that the results for A30 peptides deviate from the universal curve at short to intermediate times, but eventually join the collapsed curve. 0 500 1000 1500 2000 2500 3000 X 0 250 500 750 1000 1250 1500 1750 2000 Y 8 16 32 48 16(A 30 ) n Tube2 0 500 1000 1500 2000 2500 3000 X 0 250 500 750 1000 1250 1500 1750 2000 Y 8 16 32 48 16(A 30 ) n Tube3 0 500 1000 1500 2000 2500 3000 X 0 250 500 750 1000 1250 1500 1750 2000 Y 8 16 32 48 16(A 30 ) n Sphere 0 500 1000 1500 2000 2500 3000 X 0 250 500 750 1000 1250 1500 1750 2000 Y 8 16 32 48 16(A 30 ) n Bulk 66 4.3.4 DISTRIBUTION OF VARIOUS TYPES OF HYDROGEN BONDS An important aspect of protein aggregation is the formation of hydrogen bonds (HBs). A HB is of α-helix type if the two beads connected by the bond are in the same peptide and their sequence interval is four. Whether a HB is of β-sheet type, or simply β type, depends on the dihedral angles 𝜙 and 𝜓 shown in Figure 4-11. If −160° < 𝜙 < −60° and 90° < 𝜙 < 180°, the HB is of β type. Other types of HBs include the β-turns, π-helix, and others. So long as a peptide is isolated, the α-HBs are formed spontaneously with their number increasing over time, but as the proteins diffuse in the solution and come close to each other, the aggregates and β-sheets are formed and, thus, the number of β and other types of HBs increases as well, while the number of α-HBs decreases. Note that A 10 peptides that we simulate, which have the configuration (PHHHHHHHHP) (where P represents a polar side-chain and H represents a hydrophobic side-chain), are more amenable to interacting with each other if they are even relatively close. The amino acids in the middle form more stable HBs when compared with the N and C terminals, because the middle amino acids have two neighboring beads to help stabilizing the newly-formed HBs, while the terminals only have one. Figure 4-11 Definition of the dihedral angles 𝜑 and 𝜓 . 67 Confinement and the size of peptides play important roles in the formation of HBs. Figure 4-12 presents the time-dependence of 𝑁 HB , the number of HBs formed during the aggregation of A 10 peptides. As the results indicate, in the cases of the bulk and spherical cavity systems, the α and β HBs have identical profiles, and those numbers in the two tubes are slightly smaller, which explains the slower rate of the formation of β-sheets in the two tubes. Figure 4-12 Dynamic evolution of the number of various types of hydrogen bonds formed during aggregation of the A 10 peptides in various geometries. Figure 4-13 presents the corresponding results for A 30 peptides, manifesting significant differences with those shown in Figure 4-12. In Figure 4-12, HBs form much slower. Besides, we can observe that, at the same period of time when the steady state has been reached under the bulk conditions, the A 30 peptides 0 10 20 30 40 50 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Tube2 0 10 20 30 40 50 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Tube3 0 10 20 30 40 50 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Bulk 0 10 20 30 40 50 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Sphere 68 in the spherical cavity, however, are still forming β and other types of HBs other than α, and similarly in the two tubes. Figure 4-13. Same as in Figure 4-12, but for the hydrogen bonds formed during the aggregation of A 30 peptides. Note the significant differences with Figure 4-12, which is due to the size of peptides. 4.3.5 EFFECT OF THE PEPTIDE ’S SIZE We emphasize that all the results that have been presented so far for the A 10 and A 30 peptides also hold true for the peptides of other lengths. Figure 4-14 depicts the number of peptides having various lengths ℓ that are involved in the formation of β-sheets in the spherical cavity and in Tube3. All the cases share the same mass density of proteins. The reason we choose these two geometries is because the 0 20 40 60 80 100 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Tube2 0 20 40 60 80 100 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Tube3 0 20 40 60 80 100 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Bulk 0 20 40 60 80 100 t (ns) 0 50 100 150 200 250 300 350 400 N HB a b Other Total type Sphere 69 aggregations in the spherical cavity and Tube3 have, respectively, the least and the most fluctuations. The trends of the plots are also clear: the longer the peptides, the slower the aggregation. Figure 4-14 Comparison of the effect of peptide size ℓ on the number of peptides N β that are involved in the formation of β-sheets. Figure 4-15 Comparison of the effect of peptide size ℓ on the total number of aggregates N(t). Figure 4-14 indicates that the peptide length has a significant effect on the formation of β-sheets. However, comparing with that shown in Figure 4-15, it seems not important to the aggregation. Based on the results presented so far (also see below), they show that the mass density of proteins, instead, is 70 a primary factor that influences the aggregation, whereas the length, at least for the peptides that we have simulated, does not have a significant effect on initiating aggregation. 4.3.6 EFFECT OF THE CONFINED MEDIUM ’S SIZE Clearly, the size of a confined medium relative to the peptide’s length in its most stretched configuration plays a very important role in the aggregation, the formation of β-sheets and eventually fibrils. Thus, we run a series of simulations in the spherical cavities and tubes of various diameters. The results for 𝑁 (𝑡 ) and 𝑁 β (𝑡 ) during the aggregation of A 10 and A 30 in Tube3 and the spherical cavity are shown in Figure 4- 16 and Figure 4-17. In the case of A 10 peptides, due to their small size and relatively rigid structure — larger peptides of A 30 type are more flexible — they are able to diffuse relatively freely even in the tubes, hence indicating that the effect of confinement on their formation of β-sheets is weak, while the mass density of proteins plays a more important role, i.e. higher density accelerates initiating the aggregation. In the case of A 30 peptides, the effect of the confinement size on the aggregation, although similar to that of A 10 , is more significant in Tube3 than that in the spherical cavity. This has been discussed before in the geometry effect on the aggregation. However, it strongly influences the formation of β-sheets. Our simulations indicate that in roughly half of the cases, large peptides do not form any β-sheet even after 100 ns of simulations; only amorphous aggregates are formed (again, the average data shown here only counts those cases eventually formed β-sheet during the simulations). However, once the β-sheet nucleus forms, it grows very fast and is stabilized by the confined environment. Similar results are obtained in both Tube3 and the spherical cavity. These observations are strikingly different from those in the bulk conditions. 71 Figure 4-16 Effect of the diameter D of Tube3 on the total number of aggregates N(t) and the number of peptides N β (t) that are involved in the formation of β-sheets. Left and right figures show, respectively, the results of A 10 and A 30 peptides. Note the significant differences between the two cases in β-sheet formation, which is due to the size of the peptides. 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N b (t) 4.2 5.0 5.9 6.7 7.5 D Tube3-A 10 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N b (t) 6.0 7.3 8.5 9.7 10.9 D Tube3-A 30 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N(t) 4.2 5.0 5.9 6.7 7.5 D Tube3-A 10 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N(t) 6.0 7.3 8.5 9.7 10.9 D Tube3-A 30 72 Figure 4-17 Same as in Figure 4-16, but in a spherical geometry. In a previous paper 114 we studied the folding of α-helix proteins of various lengths and demonstrated for the first time, to our knowledge, that nanopores with purely repulsive walls destabilize smaller α-helices, reducing their folding temperatures comparing with those in bulk solutions. The size of nanopores, however, does not affect the structure of folded state. In addition, we demonstrated that the destabilization is accompanied unexpectedly by entropic stabilization of the misfolded β-strand like states. The entropic effect was shown to be stronger for more severe confining conditions. In the present work the aggregated states are the results of misfolded or partially misfolded proteins joining together and behaving like a large misfolded cluster. Thus, the same type of conclusion is also applicable to them. In particular, if there exists some state of the aggregated proteins that has a geometry similar to the confined media, such as β-strands in the slit pores, then, such aggregates will 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N b (t) 10 12 14 16 18 D Sphere-A 30 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N(t) 6.9 8.3 9.7 11.1 12.5 D Sphere-A 10 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N(t) 10 12 14 16 18 D Sphere-A 30 0 10 20 30 40 50 t (ns) 0 2 4 6 8 10 N b (t) 6.9 8.3 9.7 11.1 12.5 D Sphere-A 10 73 have larger entropy in the confinement comparing with a case in which all the proteins fold. In such a case the aggregates compete strongly with the folded states, which is in contrast to that in an unbounded fluid, similar to what we studied previously. 114 Moreover, in general, the aggregated states may have smaller volumes than the completely folded ones, which implies again that they lose smaller entropy in a confined medium than in bulk. In the other words, the severely confined media are beneficial to the aggregates’ stability. 4.4 DYNAMIC SCALING AND FRACTAL PROTEIN AGGREGATES The dynamic scaling in the aggregation of proteins that we have described is characterized by power laws, which are characteristics of fractal structures. Indeed, as it is well known, 158 the aggregation of particles leads to the formation of self-similar fractal structures. Therefore, the question that arises is whether the protein aggregates represent fractal structures. To address this question, we computed the direct correlation function 𝐶 (𝑟 ), defined by 𝐶 (𝑟 ) = 1 𝑉 ∑ 𝑠 (𝑟 ′ )𝑠 (𝑟 + 𝑟 ′ ) 𝑟 ′ (Eqn 4-6) where 𝑉 is the volume of system, and 𝑆 (𝑟 ) = 1 if a point at r belongs to the aggregate, and 𝑆 (𝑟 ) = 0 otherwise. For a self-similar fractal aggregate 𝐶 (𝑟 ) decays as a power law 𝐶 (𝑟 )~𝑟 𝐷 𝑓 −𝑑 (Eqn 4-7) where 𝑟 = |𝑟 |, and 𝐷 𝑓 and 𝑑 are, respectively, the fractal dimension of the aggregates and the Euclidean dimension of the space in which the aggregates form. 𝐶 (𝑟 ) is computed for two cases, namely, the largest aggregates that are formed as the result of the aggregations of 48 A 10 and 16 A 30 proteins, both in a spherical cavity. The results are presented in Figure 4-18. From the slopes of straight lines in the logarithmic plots of 𝐶 (𝑟 ) versus 𝑟 we obtain with 𝑑 = 3, 𝐷 𝑓 ≈ 1.48 ± 0.03 for both cases. One may also estimate 𝐷 𝑓 at various time intervals as the aggregation proceeds, which provides insight into how the structure of aggregates evolves. Thus, we also estimated 𝐷 𝑓 after several time intervals. We found that 74 the aggregates are compact at the early stages of aggregations but become sparser and take on a fractal structure at later stages. Figure 4-18 Logarithmic plots of the direct correlation functions C(r) for the aggregates formed by A 10 and A 30 peptides Fodera et al. 159 proposed a theoretical model of an aggregate formed as a result of the aggregation of charged proteins. They expressed the free energy of the aggregate as a function of two key parameters, the effective number of charges on a single protein, 𝑛 𝑐 , and the fractal dimension 𝐷 𝑓 . Due to the coarse nature of our model of proteins and use of DMD computations, the aggregation process that we simulate corresponds to the limit 𝑛 𝑐 = 0. We do, however, find that our results are in qualitative agreement with the prediction of their model in that particular case, as we find that in the initial stages of aggregation the aggregate is compact, but as more proteins join the growing cluster, the growing aggregate takes on a sparser structure with 𝐷 𝑓 < 3. This is also in agreement with the results that we reported in Chapter 3 for the aggregation in unbounded fluids, which indicated that one obtains compact, sphere-like clusters at the early stages of aggregation (see Figure 3-2) and later they convert to β-sheets and fibrils, which are sparser. Fodera et al. 159 also demonstrated that the transition between fibrils, spherulites and amorphous aggregates can be described by a multifractal model, by which they meant that 𝐷 𝑓 changes as the aggregation proceeds. This is also in qualitative agreement with our simulations. 75 The aggregation process that we study in this work is of reversible type, because the bonds formed between the peptides can be broken and formed again. Zaccone et al. 160 studied the aggregation and gelation of colloidal particles as a function of attraction energy between particles. Their model included the effect of thermal breakup and, therefore, the process was reversible and similar to what we study here. The model predicted that the growth of compact aggregates, as opposed to the more open structures that we obtain at long times, depends on the depth of attraction between peptides. Since PRIME, the model that we use for the proteins, however, is designed for proteins with specific atoms and beads, the depth of attraction is fixed. This related work could be done in the future by using a more generic model, which would not discuss in detail here. 4.5 DIFFERENCES BETWEEN AGGREGATIONS IN CONFINED AND UNBOUNDED MEDIA The universality of the dynamic scaling should not, however, mask the important differences between protein aggregations in confined media and in unbounded fluids. Aside from the fact that large proteins diffuse in confined media much slower than in the bulk, there are significant differences between the two cases, a summary of which is as follows: (1) Provided that the peptides are large enough, which in our simulations means A 20 or larger, there are significant variations of time for the nucleation of β-sheets. This agrees with the experimental observations 161-162 in crowded cellular environments. This does not happen in the simulations in unbounded fluids. (2) While the aggregates in unbounded domains are usually long single-layer β-sheets, or at most double-layer fibrils, the small peptides, such as A 10 , form short fibrils with more than two layers in confined media, whereas much larger peptides form rolled/folded aggregates. (3) If the confinement is severe, the formation of β-sheets and fibrils may be postponed for a long time, or may not even occur at all, although the latter case remains to be tested by experiments. This is not the case in unbounded fluids. (4) The geometry of a confined medium has a strong influence on the rate of aggregation. Due to the perfect symmetry of spherical cavities, the fluctuations of the local mass density are very 76 small, as a result of which the peptides aggregate faster than in circular tubes, which the latter cases contain much larger fluctuations that delay the formation of aggregates. (5) Although both spherical and cylindrical systems exhibit significant variations in the time that it takes to initiate the nucleation of β-sheets, once the nucleus is formed, the peptides in spherical cavities, and perhaps in all symmetric confined media, grow much faster into β-sheets/fibrils than in those geometries as cylindrical tubes with large aspect ratios. (6) Although the results can be represented by the dynamic scaling in Eqn 4-1 and Eqn 4-3, and even though the scaling exponents 𝑧 and 𝘁 for the spherical cavity agree with those in an unbounded fluid, the same is not true comparing with the case in cylindrical tube. The exponents exhibit sensitivity to the geometry of confined media. (7) As the biological systems in which protein aggregation occurs are unlikely to have perfect symmetry, the results indicate the significance of the geometry of confined media, and the invalidity of the common assumption that the most important features of protein aggregation can be understood by studying the phenomenon in an unbounded medium alone. 4.6 SUMMARY In this chapter, the second in a series devoted to molecular modeling of protein aggregation, we utilize PRIME, a mesoscale model of proteins, together with extensive discontinuous molecular dynamics simulation to study the phenomenon in confined media as models of the crowded cellular environments. The confined media are represented by spherical cavities as well as cylindrical tubes with two aspect ratios. We computed several important properties of the aggregation process, including dynamic evolution of the total number of aggregates, mean aggregate size, number of peptides that contribute to the formation of β-sheets, and number of various types of hydrogen bonds. We show, similar to the study in bulk conditions in Chapter 3, that the computed properties in confined media also follow dynamic scaling. The scaling exponents 𝑧 and 𝘁 are defined in Eqn 4-1 and Eqn 4-3, respectively, and their values are listed in Table 4-1. Moreover, the estimated exponents in the spherical cavity, a perfectly symmetrical structure, are shown to agree with those under the bulk conditions. For the cases in tubes, however, the exponents appear to be sensitive to the details of the confinement 77 geometry. As the biological systems in which protein aggregation occurs are unlikely to have perfect symmetrical structures, the results indicate the significance of the confined media’s geometry on the aggregation process. The existence of dynamic scaling in the aggregation of proteins in confined media, which also exists in the aggregation of colloidal particles, does indicate that many features of protein aggregation are universal, independent of the molecular structure. Therefore, simpler models of proteins and their aggregation can as well shed much light on this important set of phenomena. 159-160 ACKNOWLEDGMENT We thank Dr. Ted Lee for advance information on the experiments by his group on the dynamic scaling of aggregation of proteins, and Dr. Leili Javidpour for useful discussions. We thank the two anonymous reviewers whose constructive comments helped us improve the quality of the paper. The computations were carried out using the computer cluster of the High-Performance Computing Center at the University of Southern California. 78 5. MOLECULAR DYNAMICS STUDY OF THE STRUCTURE OF POLY- GLYCINE-ALANINE IN AGGREGATION AND FOLDING 5.1 INTRODUCTION Expansion of a hexanucleotide GGGGCC (G 4 C 2 ) repeat in the gene chromosome 9 open reading frame 72 (C9ORF72) is widely believed the most common genetic cause of amyotrophic lateral sclerosis (ALS) and frontotemporal dementia (FTD) 163-167 . There has been three distinct, but non-exclusive hypotheses described the mechanisms of this repeat that causes the disease: (1) loss of function due to decreased expression of the C9ORF72 protein, (2) gain of function due to the accumulation of toxic RNA foci 168 , and (3) production of unnatural dipeptide repeat (DPR) proteins and their aggregates via repeat-associated non-ATG (RAN) translation 164, 166 . The third is the most recent and was proposed by Zu et al 169 , who found that although the G 4 C 2 repeat is in the non-coding region of the C9ORF72 gene, there is an unconventional mechanism of translation — repeat-associated non-ATG (RAN) translation — which can bidirectionally translate the repeat in all reading frames into five distinct DPRs: poly-glycine-alanine (poly-GA), poly-glycine-arginine (poly-GR), poly-proline-alanine (poly-PA), poly-proline-arginine (poly-PR) and poly-glycine-proline (poly-GP) 170-174 . Experiments show that these DPRs accumulate and cluster in the characteristic p62-positive TDP-43 negative inclusions found in the affected brain regions, defining C9ORF72 pathology 170-171 . Further, recent research finds that the majority of these characteristic inclusions contain poly-GA proteins, which is highly aggregation-prone compared to the other DPRs and believed substantially contributes to the neurodegeneration in ALS and FTD 172, 175-176 . Zhang et al. 175, 177 , May et al. 178 , Ash et al. 174 and Lin et al. 179 show that poly-GA proteins in cultured cells and primary neurons form insoluble high molecular weight aggregates, and demonstrate that its toxicity is associated with reduced dendritic branching, impairment of the ubiquitin-proteasome system (UPS) and induction of endoplasmic reticulum (ER) stress. Despite such facts, the underlying mechanisms are still poorly understood due to the limited structural information that is available on poly-GA aggregates. Some experiments suggested that poly-GA proteins would form flat, ribbon-type fibrils that contain a characteristic cross β-sheet structure 180 . A more recent work claimed that the aggregates consist of 79 densely-packed twisted ribbons 166 . These structures are similar to the aggregate of β-amyloid peptides observed in Alzheimer’s disease, which are believed to share a lot of similar properties with respect to both structure and mechanism 163-164, 180-181 . In the previous chapters we have used DMD combined with a coarse-grained protein model, named PRIME 81-84 , to simulate the aggregation of polyalanines, whose configuration is intended to mimic the hydrophobic part of β-amyloid, and successfully obtained fibril aggregates 110-111 . In this chapter, we use our upgraded DMD program combined with an all-atom protein model to simulate the aggregation of poly-GA proteins, with the goal of exploring the aggregate’s structure and the underlying mechanism at detailed molecular level. The structure of this chapter is as follows. In Section 2 we describe the models of proteins as well as the simulation configurations used in DMD. Section 3 presents the results. Section 4 gives a discussion and summarizes the chapter. 5.2 CONFIGURATION OF DMD SIMULATION 5.2.1 MODELS OF PROTEINS Because we use an all-atom model to perform the simulations, we can directly use a protein builder to create the protein models. Here, we use the structure builder integrated in Chimera, developed by RBVI and UCSF, to construct the poly-GA protein with 25 repeating units, i.e. (GA) 25 . Zheng et al. 175 suggests that one should use longer repeat lengths, such as 50 repeating dipeptides, to evaluate poly-GA aggregation and toxicity to better mimic C9ORF72 associated ALS/FTD. Thus, we also construct the poly- GA protein with 50 repeats, i.e. (GA) 50 . We then use GROMACS to convert the model from PDB to GRO format and remove the light hydrogens on the protein, after which the model is ready to be applied into the DMD program. The details of the all-atom model have been described in Chapter 2. The details of the potential library can be found in https://github.com/alexzheng42/sDMD. 80 5.2.2 MOLECULAR SIMULATION Two systems containing 10 (GA) 25 and 10 (GA) 50 are simulated in the simulation boxes of dimensions 12 x 12 x 12 nm 3 (10 mmol/L) and 15 x 15 x 15 nm 3 (5 mmol/L), respectively. They, however, share identical mass density, 304.9 g/L. Each system is simulated at three different temperatures, namely, 262.5 K, 275 K and 287.5 K. In fact, the temperature used in the simulation is expressed in units of a temperature scale 𝑇 𝑠 , which based on the unit system of the program it gives 𝑇 𝑠 ≈ 500 K. Thus, the corresponding dimensionless temperatures are 𝑇 ∗ = 0.525, 0.55 and 0.575, respectively. They are below, around and higher than the folding temperatures of (GA) 25 and (GA) 50 , see Figure 5-1. Figure 5-1 Folding temperatures of (GA) 25 and (GA) 50 REMD of single (GA) 25 and (GA) 50 produces their Cv* (reduced heat capacity) profiles, which indicate the folding temperatures are 278 K ( 𝑇 ∗ = 0.556, solid curve) and 272 K ( 𝑇 ∗ = 0.544, dashed curve), respectively, for our protein models. 200 250 300 350 400 T (K) 0 500 1000 1500 2000 Cv * (GA) 25 200 250 300 350 400 0 1000 2000 3000 4000 Cv * (GA) 50 81 Initially, the poly-GA proteins are randomly positioned into the simulation box. The initial velocity of each atom is given from the Maxwell-Boltzmann distribution at 250 K (𝑇 ∗ = 0.5). The system is stabilized under this low temperature for 5000 timesteps (250 ps). After that, it is heated up to 500 K (𝑇 ∗ = 1.0) and quenched by 50 K every 5000 timesteps, until it drops to the target temperature. The system then keeps running for up to 800000 timesteps (40 ns). We performed ten independent simulations for each system at each of the three temperatures in order to assess the effect of the initial conditions. The results represent the averages. 5.3 RESULTS 5.3.1 SYSTEM TEMPERATURE Temperature is one of the most important factors that affects the protein aggregation. By using replica- exchange molecular dynamics (REMD) we calculate the folding temperatures of (GA) 25 and (GA) 50 , which turns out to be 0.556 and 0.544, respectively (see Figure 5-1). However, due to the nature of DMD — the potential profile is discontinuous — the computed folding temperature may not be as precise as we would like it to be. Thus, we perform simulations at three different temperatures around the folding temperature. The snapshots of (GA) 25 and (GA) 50 systems at each temperature are shown in Figure 5-2. As we can see in Figure 5-2, for both (GA) 25 and (GA) 50 at 𝑇 ∗ = 0.525, almost all proteins have either folded by themselves or joined the aggregates. Similarly, at 𝑇 ∗ = 0.55, the majority of the proteins have folded and/or aggregated, with the rest are still in aggregation process and may need more time to form stable structures. At 𝑇 ∗ = 0.575, however, most of the proteins are still in random-coil state, and the folded and aggregated are not stable either: the intra- and inter-molecular hydrogen bonds (HBs) keep forming and breaking and the VDW attraction is week. For this reason, we would only focus on the results from the simulations at 𝑇 ∗ = 0.525 and 0.55. 82 Figure 5-2 Aggregation of (GA) 25 and (GA) 50 at different temperatures Aggregation of (GA) 25 at 𝑇 ∗ = (a) 0.525, (b) 0.55 and (c) 0.575, and (GA) 50 at 𝑇 ∗ = (d) 0.525, (e) 0.55 and (f) 0.575. The secondary structures are calculated by STRIDE algorithm in VMD. The color codes are: white = coil, green = turn, purple = α-helix, yellow = β-sheet, and red = π-helix. Red arrow points inter-molecular β-sheet, black arrow points β-helix, and blue arrow points double-helix. 5.3.2 FORMATION OF β-SHEETS Figure 5-3 shows the four stages of the aggregation of 10 (GA) 25 at 𝑇 ∗ = 0.55. At the beginning, all the proteins are in random-coil configurations. After a short time, they have begun to fold and aggregate, forming intra- and inter-molecular β-sheets. As the aggregation process continues, more and more monomers join the aggregate and keep forming more β-sheets. Since the protein is long, some residues also fold into α-helices. At the end, seven proteins have aggregated into an ordered β-sheet-rich structure. The α-helices are also forced to break and change to β-strands. The results indicate that poly- 83 GA is indeed an aggregation-prone protein and its aggregate is rich in β-sheets. Similar results are obtained in the other simulation runs, see Figure 5-2. The aggregates are stable; however, we can see that the β-strands inside usually do not align well: they are bended and twisted and may be also mixed with turns and coils in between. This could be due to the limited simulation time. In some individual simulation, we do obtain more ordered structures, such as that shown in Figure 5-4. Figure 5-3 Four stages of the aggregation of 10 (GA) 25 Initially, all the proteins are in random-coil configuration. At the last stage, seven out of ten proteins have aggregated into an ordered β-sheet-rich structure. 84 Figure 5-4 Detail view of an ordered β-sheet aggregate. The dynamic evolution of the total number of aggregates is shown in Figure 5-5. At the higher temperature 𝑇 ∗ = 0.55, the systems have obviously larger fluctuations, which is due to the continuous formations and destructions of the inter-molecular interactions, especially the HBs. Thus, the systems aggregate slower. On the other hand, although the systems of (GA) 25 and (GA) 50 have different concentrations, 10 mmol/L and 5 mmol/L, respectively, they aggregate with the same trend. This is consistent with our observation in Chapter 4 that the mass density of proteins, which here is identical, is the primary factor that influences the aggregation, whereas the protein’s length does not have a significant effect on initiating aggregation. 85 Figure 5-5 Dynamic evolution of the total number of aggregates N(t) formed (a) at 𝑇 ∗ = 0.525, (b) at 𝑇 ∗ = 0.55 5.3.3 UNIQUE HELICAL STRUCTURES In all the aggregations we notice some unique structures that are worth paying more attentions, see Figure 5-6. Unlike the regular secondary structures of folded and aggregated proteins, e.g. α-helix and β-sheet (Figure 5-6 (a), (b)), we observe a structure of helix but constructed by β-sheet type HBs instead of α-helix type HBs in Figure 5-6 (c), and a structure of β-hairpin but is twisted into a double-helix shape in Figure 5-6 (d). Figure 5-6 shows the renderings of different secondary structures along with the Ramachandran plots of their glycine and alanine residues. Ramachandran plot, developed by Ramachandran et al. 127 , is used to demonstrate the energetically allowed regions for dihedral angels 𝜙 and 𝜓 (see Figure 5-7) of amino acid backbone. It also provides a good way to visualize the distribution of dihedral angels in a protein structure. Different types of secondary structures and HBs have their specific favored regions in Ramachandran plot. 86 Figure 5-6 Different secondary structures and the corresponding Ramachandran plots 87 Figure 5-7 Diagram of dihedral angels 𝜑 and 𝜓 in a peptide. In a polypeptide, the N–Cα and Cα–C bonds in backbone are relatively free to rotate. The rotations are represented, respectively, by the torsion angles 𝜓 and 𝜙 . Here the hydrogen and carbonyl oxygen are omitted, and the red bead represent the side-chain. In α-helix, the dihedral angle coordinates (𝜙 , 𝜓 ) of both alanine and glycine cover the area (−89 < 𝜙 < −39°, −66 < 𝜓 < −16°), which is located at lower-left of the center in a Ramachandran plot. In β-sheet, however, it becomes a little more complex: for non-glycine amino acids, for both parallel and anti- parallel types, the coordinates would cover the area (−130 < 𝜙 < −111°, +119 < 𝜓 < +145°) at the top-left of a Ramachandran plot; while for glycine, except the area at the top-left corner, it also covers the area at the bottom-right (+135 < 𝜙 < +180°, −180 < 𝜓 < −135°), which the latter can be considered as a continuation of the β-sheet region 182 . (Note that the angle ranges shown here are the most energetically favored, but not necessary the only allowed) For the β-helix and double-helix in Figure 5-6 (c) and (d), the dihedral angles of glycine drop at the bottom-right corner, while those of alanine drop at the top-left corner. This could be used to characterize the formations of the two helices, and the distinction between them is that β-helix forms parallel β- sheets, whereas double-helix forms anti-parallel β-sheet. We can see more details of the structures in the HB contact map in Figure 5-8. For an α-helix, the HBs form between the amino acids at 𝑖 and 𝑖 ± 4, see Figure 5-9 (a). Therefore, its pattern in the contact map is a pair of parallel strips directly adjacent to the diagonal, 𝑥 = 𝑦 . For a β-helix, on the other hand, unlike α-helix that forms HBs between 𝑖 and 𝑖 ± 4, 88 and π-helix between 𝑖 and 𝑖 ± 5 (very unstable), each glycine 𝑖 forms two HBs with the alanines at 𝑖 + 5 and 𝑖 + 7 simultaneously (similar to HBs in a parallel β-sheet, which are between 𝑖 and 𝑗 − 1, 𝑗 + 1), see Figure 5-9 (c). Its pattern, therefore, is two pairs of parallel strips adjacent to the diagonal, but separated wider than α-helix. For a double-helix, each amino acid 𝑖 forms two HBs with another amino acid 𝑗 , which follows the anti-parallel β-sheet pattern in the contact map as cross-diagonal strips (parallel to the diagonal, 𝑥 = −𝑦 ). But it is also different from the regular β-hairpin. Figure 5-8 shows that for β-hairpin, there is only one cross-diagonal strip on each side of the map, while there are two for double-helix. This is due to the helical structure that allows adjacent amino acids, 𝑖 − 1, 𝑖 and 𝑖 + 1 (not those near the terminals), to form HBs with the amino acids at 𝑗 , 𝑘 and 𝑗 − 2, respectively, where 𝑗 and 𝑗 − 2 located in the upper turn and 𝑘 in the lower turn of the helix (𝑘 < 𝑗 − 2), see Figure 5-9 (d). Therefore, every amino acid can form HBs, instead of only half of them can form HBs in a regular β-hairpin 89 , see Figure 5-9 (b); similar in β-helix. Figure 5-8 HB contact map including different secondary structures 89 Figure 5-9 Structure details of different secondary structures. (a) α-helix. Each amino acid 𝑖 forms two HBs, one with amino acid 𝑖 − 4 and the other with 𝑖 + 4. Glycine forms HBs with glycine and alanine forms HBs with alanine. (b) β-hairpin. Anti-parallel β-sheet. Amino acid 𝑖 forms two HBs with another amino acid 𝑗 . But only half of the amino acids form HBs. (c) β-helix. Parallel β-sheet. Each glycine 𝑖 forms two HBs, one with the alanine at 𝑖 + 5 and the other with the alanine at 𝑖 + 7. (d) double-helix. Anti-parallel β-sheet. Each amino acid 𝑖 forms two HBs with another amino acid 𝑗 . The white numbers are the lengths of hydrogen bonds, units of Å. Formation of more HBs indicates that such helical structure is energetically preferred. As stated above, for a protein with 𝑛 amino acids, when it folds into a β-hairpin, at most 𝑛 /2 of the amino acids can form HBs, while it folds into an α-helix, it could form 𝑛 − 4 HBs, with the latter being much more stable (𝑛 > 90 8). The transition from β-sheet to α-helix is coupled with a transition to random-coil 89, 182 , which is a very slow process at low temperature. However, benefitted from the helical structure, a β-helix and double- helix can have the other 𝑛 /2 amino acids also form HBs without transitions and become as stable as an α-helix. This can be seen in Figure 5-10. Figure 5-10 Potential energy and hydrogen bond number of each type of structures A (GA) 10 unit of each structure at T* = 0.55 is picked and its potential energy and HB number are calculated during the last 1.5 ns of simulations. 91 In Figure 5-10, a (GA) 10 unit of each structure is used to calculate the intra-molecular potential energy and HB number. The results of helical structures are very similar: the potential energy keeps around -75 kcal/mol, and the average HB number is around 19, which indicates that just about all the amino acids in the repeats have formed HBs (only the HBs from amino groups on backbone are counted). However, the data from β-hairpin are much different: its potential is much higher with an average value of about -35 kcal/mol and has significant less HBs formed — only 7 intra-molecular HBs. In other words, β-hairpin is less stable comparing with the helical structures. These results all indicate that β-helix and double- helix can share the characteristics of both α-helix and β-sheet. 5.4 DISCUSSION AND SUMMARY The two special helical structures — β-helix and double-helix — frequently form in the simulations at 𝑇 ∗ = 0.525, 0.55 and even 0.575. It is important to understand why these types of structures would form in ploy-GA proteins. It is presumably due to the glycine in the repeat. Figure 5-6 shows that in β-helix and double-helix the dihedral angles of glycine drop at the bottom-right corner instead of the regular top- left, which region only glycine is able to reach in all 20 types of amino acids. To support this hypothesis, we also simulate a DPR protein consists of 25 repeats of glycine-arginine, i.e. (GR) 25 , which is another toxic DPR protein associated with ALS and FTD. Interestingly, we obtain the same structures, see Figure 5-11. Thus, we propose that the repeats of (GX) n , where X could be any non-glycine amino acid and n is the repeat length, would be the root of these unique helical structures, and it is the glycine which has exceptional wide ranges of dihedral angels that makes the structures possible. 92 Figure 5-11 Identical structures of β-helix (left) and double-helix (right) formed in poly-GR system. In summary, in this chapter we use the upgraded DMD program combined with an all-atom molecular model to simulate the poly-GA protein, which is widely believed the main toxic DPR protein associated with ALS and FTD, and investigate its structure in aggregation and folding. We successfully obtain β- sheet-rich aggregates and at the same time observe two special helical structures, β-helix and double- helix, during the aggregation process. These two structures frequently form in the simulations of poly- GA and poly-GR systems at different temperatures. The analysis shows the details of these helical structures and indicates they are very stable. We propose that proteins having the sequence that contains repeat (GX) n , where X could be any non-glycine amino acid and n is the repeat length, would be capable of forming these two structures of β-helix and double-helix, and it is the glycine in the repeat that contributes to this characteristic. 93 6. APPENDIX This appendix will describe how DMD works. Assume we have two particles, 𝑖 and 𝑗 , each with diameter σ. The coordinates of the particles are the vector 𝑟 ⃗ 𝑖 = [𝑥 𝑖 , 𝑦 𝑖 , 𝑧 𝑖 ] and 𝑟 ⃗ 𝑗 = [𝑥 𝑗 , 𝑦 𝑗 , 𝑧 𝑗 ] and they have velocity vector 𝑣 ⃗ 𝑖 and 𝑣 ⃗ 𝑗 (see Figure 6-1). The collision event and collision time can be determined by using the relative position and relative velocity, 𝑟 ⃗ 𝑖𝑗 and 𝑣 ⃗ 𝑖𝑗 . Figure 6-1 Sample system of a collision of a pair of atoms (2D) A collision event is dependent of whether the distance between particles 𝑖 and 𝑗 will become equal to the average molecular diameter 𝜎 𝑖𝑗 = (𝜎 𝑖 + 𝜎 𝑗 )/2 sometime in the future. Using superscript “o” to represent an initial condition, the relative position vector at any time in the future is found by the following equation, 𝑟 ⃗ 𝑖𝑗 = 𝑟 ⃗ 𝑖𝑗 𝑜 + 𝑣 ⃗ 𝑖𝑗 𝑜 (𝑡 − 𝑡 𝑜 ) (Eqn 6-1) So, if 𝑟 ⃗ 𝑖𝑗 and 𝑣 ⃗ 𝑖𝑗 have opposite sense, the particles will be approaching and a collision is possible. We may check to see if the dot product 𝑟 ⃗ 𝑖𝑗 ∙ 𝑟 ⃗ 𝑖𝑗 becomes equal to 𝜎 𝑖𝑗 2 in the future at some time increment (𝑡 − 𝑡 𝑜 ), 94 𝜎 𝑖𝑗 2 = 𝑟 ⃗ 𝑖𝑗 ∙ 𝑟 ⃗ 𝑖𝑗 (Eqn 6-2) 𝜎 𝑖𝑗 2 = [𝑟 ⃗ 𝑖𝑗 𝑜 + 𝑣 ⃗ 𝑖𝑗 𝑜 (𝑡 − 𝑡 𝑜 )] ∙ [𝑟 ⃗ 𝑖𝑗 𝑜 + 𝑣 ⃗ 𝑖𝑗 𝑜 (𝑡 − 𝑡 𝑜 )] (Eqn 6-3) 𝜎 𝑖𝑗 2 = 𝑟 ⃗ 𝑖𝑗 𝑜 ∙ 𝑟 ⃗ 𝑖𝑗 𝑜 + 2𝑟 ⃗ 𝑖𝑗 𝑜 𝑣 ⃗ 𝑖𝑗 𝑜 (𝑡 − 𝑡 𝑜 ) + 𝑣 ⃗ 𝑖𝑗 𝑜 𝑣 ⃗ 𝑖𝑗 𝑜 (𝑡 − 𝑡 𝑜 ) 2 (Eqn 6-4) Here we define a variable 𝑏 𝑖𝑗 as a dot product, 𝑏 𝑖𝑗 ≡ 𝑟 ⃗ 𝑖𝑗 ∙ 𝑣 ⃗ 𝑖𝑗 = |𝑟 ⃗ 𝑖𝑗 ||𝑣 ⃗ 𝑖𝑗 | cos 𝛽 (Eqn 6-5) If 𝑏 𝑖𝑗 > 0, the relative position |𝑟 ⃗ 𝑖𝑗 | will be increasing with time, which is to say, the two particles will move further away with each other. If 𝑏 𝑖𝑗 = 0, the relative velocity and relative position vectors are perpendicular. If 𝑏 𝑖𝑗 < 0, the relative position |𝑟 𝑖𝑗 | will be decreasing and there may be a collision. So, after removing the superscript “o” on 𝑟 ⃗ 𝑖𝑗 and 𝑣 ⃗ 𝑖𝑗 , Eqn 6-4 can be rewritten as follows, 𝜎 𝑖𝑗 2 = 𝑟 𝑖𝑗 2 + 2𝑏 𝑖𝑗 (𝑡 − 𝑡 𝑜 ) + 𝑣 𝑖𝑗 2 (𝑡 − 𝑡 𝑜 ) 2 (Eqn 6-6) or 𝑣 𝑖𝑗 2 (𝑡 − 𝑡 𝑜 ) 2 + 2𝑏 𝑖𝑗 (𝑡 − 𝑡 𝑜 ) + (𝑟 𝑖𝑗 2 − 𝜎 𝑖𝑗 2 ) = 0 (Eqn 6-7) Solving for time increment by using the quadratic formula, we can have, 𝑡 𝑐 = 𝑡 − 𝑡 𝑜 = −𝑏 𝑖𝑗 + 𝑠 time √𝑏 𝑖𝑗 2 − 𝑣 𝑖𝑗 2 (𝑟 𝑖𝑗 2 − 𝜎 𝑖𝑗 2 ) 𝑣 𝑖𝑗 2 (Eqn 6-8) where 𝑠 time = ±1. The discriminant, 𝑏 𝑖𝑗 2 − 𝑣 𝑖𝑗 2 (𝑟 𝑖𝑗 2 − 𝜎 𝑖𝑗 2 ), in Eqn 6-8 has to be non-negative for a collision to occur. If it is negative, then the particles will miss each other. The smaller real root with 𝑠 time = −1 is used for the solution to the quadratic as shown in Figure 6-2. The smaller real root occurs when the particles collide at 𝑟 ⃗ 𝑗 𝑠 − . The location labeled 𝑟 ⃗ 𝑗 𝑠 + is the root that will not occur physically. 95 Figure 6-2 Diagram of the two roots for the solution to the quadratic formula used to calculate the future collision time. Only the smaller real root with 𝑠 time = −1 is used. The other root would not occur physically. Figure 6-3 Diagram of the different collision situations. (a) Approaching particles that will experience a potential change event at the shell boundary; (b) particles in the shell that will have a hard-sphere collision; (c) particles in the shell that will miss a hard-sphere collision, but experience a potential change event after they pass; (d) particles moving away that will experience a potential change event at the shell boundary. If the particle has a potential shell (see Figure 6-3), the formula to calculate the future event time can be modified to, 𝑡 𝑐 = 𝑡 − 𝑡 𝑜 = −𝑏 𝑖𝑗 + 𝑠 time √𝑏 𝑖𝑗 2 − 𝑣 𝑖𝑗 2 (𝑟 𝑖𝑗 2 − 𝑑 𝑖𝑗 2 ) 𝑣 𝑖𝑗 2 (Eqn 6-9) where 𝑑 𝑖𝑗 is the actual distance of centers for the event. 𝑟 ⃗ 𝑗 𝑠 + 𝑣 ⃗ 𝑖𝑗 j i 𝑟 ⃗ 𝑗 𝑠 − (a) (b) (c) (d) 96 Calculation of the velocity changes upon an event will be a little more complex. We use prime (‘) to denote the state after the event. The change in potential energy for the event is indicated by ∆𝘀 = 𝘀 ′ − 𝘀 , which may be zero, positive or negative depending on whether it is a hard-sphere collision, or particles enter or escape from a potential shell (in potential well situation), see Figure 6-3. Total energy is conserved upon an event between particles 𝑖 and 𝑗 , 1 2 𝑚 𝑖 𝑣 ⃗ 𝑖 ∙ 𝑣 ⃗ 𝑖 + 1 2 𝑚 𝑗 𝑣 ⃗ 𝑗 ∙ 𝑣 ⃗ 𝑗 + 𝘀 = 1 2 𝑚 𝑖 𝑣 ⃗ 𝑖 ′ ∙ 𝑣 ⃗ 𝑖 ′ + 1 2 𝑚 𝑗 𝑣 ⃗ 𝑗 ′ ∙ 𝑣 ⃗ 𝑗 ′ + 𝘀 ′ (Eqn 6-10) which can be rearranged, 1 2 𝑚 𝑖 (𝑣 ⃗ 𝑖 − 𝑣 ⃗ 𝑖 ′ ) ∙ (𝑣 ⃗ 𝑖 + 𝑣 ⃗ 𝑖 ′ ) = − 1 2 𝑚 𝑗 (𝑣 ⃗ 𝑗 − 𝑣 ⃗ 𝑗 ′ ) ∙ (𝑣 ⃗ 𝑗 + 𝑣 ⃗ 𝑗 ′ ) + ∆𝘀 (Eqn 6-11) When the two particles contact, the relative position will be, 𝑟 ⃗ 𝑖𝑗 𝑐 = 𝑟 ⃗ 𝑖 𝑐 − 𝑟 ⃗ 𝑗 𝑐 (Eqn 6-12) The change in momentum for an event is proportional to the collision vector 𝑟 ⃗ 𝑖𝑗 𝑐 , 𝑚 𝑖 (𝑣 ⃗ 𝑖 ′ − 𝑣 ⃗ 𝑖 ) = −𝑚 𝑗 (𝑣 ⃗ 𝑗 ′ − 𝑣 ⃗ 𝑗 ) ∝ 𝑟 ⃗ 𝑖𝑗 𝑐 (Eqn 6-13) or (𝑣 ⃗ 𝑖 ′ − 𝑣 ⃗ 𝑖 ) 𝑚 𝑗 = − (𝑣 ⃗ 𝑗 ′ − 𝑣 ⃗ 𝑗 ) 𝑚 𝑖 = 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 (Eqn 6-14) where 𝜙 is the parameter we will try to determine. Insert Eqn 6-14 into Eqn 6-11, 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 ∙ (𝑣 ⃗ 𝑖 + 𝑣 ⃗ 𝑖 ′ ) = 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 ∙ (𝑣 ⃗ 𝑗 + 𝑣 ⃗ 𝑗 ′ ) − 2∆𝘀 (𝑚 𝑖 𝑚 𝑗 ) (Eqn 6-15) Rearrange Eqn 6-14 also leads to, 𝑣 ⃗ 𝑖 ′ = 𝑣 ⃗ 𝑖 + 𝑚 𝑗 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 (Eqn 6-16) and 97 𝑣 ⃗ 𝑗 ′ = 𝑣 ⃗ 𝑗 − 𝑚 𝑖 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 (Eqn 6-17) Plugging Eqn 6-16 and Eqn 6-17 into Eqn 6-15 to eliminate the primed variables, 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 ∙ (2𝑣 ⃗ 𝑖 + 𝑚 𝑗 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 ) = 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 ∙ (2𝑣 ⃗ 𝑗 − 𝑚 𝑖 𝜙 𝑟 ⃗ 𝑖𝑗 𝑐 ) − 2∆𝘀 (𝑚 𝑖 𝑚 𝑗 ) (Eqn 6-18) Define 𝑑 𝑖𝑗 2 ≡ 𝑟 ⃗ 𝑖𝑗 𝑐 ∙ 𝑟 ⃗ 𝑖𝑗 𝑐 , 𝑣 ⃗ 𝑖𝑗 ≡ 𝑣 ⃗ 𝑖 − 𝑣 ⃗ 𝑗 and 𝑏 𝑖𝑗 𝑐 ≡ 𝑟 ⃗ 𝑖𝑗 𝑐 ∙ 𝑣 ⃗ 𝑖𝑗 , and rearrange Eqn 6-18, (𝑚 𝑖 + 𝑚 𝑗 )𝑑 𝑖𝑗 2 𝜙 2 + 2𝑏 𝑖𝑗 𝑐 𝜙 + 2∆𝘀 (𝑚 𝑖 𝑚 𝑗 ) = 0 (Eqn 6-19) Applying the quadratic formula, 𝜙 = −2𝑏 𝑖𝑗 𝑐 ± √ 4(𝑏 𝑖𝑗 𝑐 ) 2 − 8(𝑚 𝑖 + 𝑚 𝑗 )𝑑 𝑖𝑗 2 ∆𝘀 /(𝑚 𝑖 𝑚 𝑗 ) 2(𝑚 𝑖 + 𝑚 𝑗 )𝑑 𝑖𝑗 2 (Eqn 6-20) Define the reduced mass, 𝑚 𝑟 = 𝑚 𝑖 𝑚 𝑗 𝑚 𝑖 +𝑚 𝑗 , then, 𝜙 = −𝑏 𝑖𝑗 𝑐 + 𝑠 vel √ (𝑏 𝑖𝑗 𝑐 ) 2 − 2𝑑 𝑖𝑗 2 ∆𝘀 /𝑚 𝑟 (𝑚 𝑖 + 𝑚 𝑗 )𝑑 𝑖𝑗 2 (Eqn 6-21) where the sign of 𝑠 vel depends on the context of the event (see Figure 6-4). For an event where ∆𝘀 = 0, 𝑠 vel = −1 is the only reasonable solution. For a hard-sphere collision, the distance is the inter-particle diameter, i.e. 𝑑 𝑖𝑗 = 𝜎 𝑖𝑗 , and the relation becomes, 𝜙 = −2𝑏 𝑖𝑗 𝑐 (𝑚 𝑖 + 𝑚 𝑗 )𝜎 𝑖𝑗 2 (Eqn 6-22) Based on the context of the event, DMD calculation will follow the flowsheet below, 98 Figure 6-4 Flowsheet of soft-sphere collision interactions in DMD. 183 99 REFERENCE 1. Dobson, C. M., Protein misfolding, evolution and disease. Trends Biochem Sci 1999, 24 (9), 329- 32. 2. Sipe, J. D., Amyloidosis. Annu Rev Biochem 1992, 61, 947-75. 3. Ross, C. A.; Poirier, M. A., Protein aggregation and neurodegenerative disease. Nat Med 2004, 10 Suppl, S10-7. 4. Stefani, M.; Dobson, C. M., Protein aggregation and aggregate toxicity: new insights into protein folding, misfolding diseases and biological evolution. J Mol Med (Berl) 2003, 81 (11), 678-99. 5. Zhang, S., Fabrication of novel biomaterials through molecular self-assembly. Nat Biotechnol 2003, 21 (10), 1171-8. 6. Koga, T.; Higuchi, M.; Kinoshita, T.; Higashi, N., Controlled self-assembly of amphiphilic oligopeptides into shape-specific nanoarchitectures. Chemistry 2006, 12 (5), 1360-7. 7. Frokjaer, S.; Otzen, D. E., Protein drug stability: a formulation challenge. Nat Rev Drug Discov 2005, 4 (4), 298-306. 8. Chi, E. Y.; Krishnan, S.; Randolph, T. W.; Carpenter, J. F., Physical stability of proteins in aqueous solution: mechanism and driving forces in nonnative protein aggregation. Pharm Res 2003, 20 (9), 1325- 36. 9. Cellmer, T.; Bratko, D.; Prausnitz, J. M.; Blanch, H. W., Protein aggregation in silico. Trends Biotechnol 2007, 25 (6), 254-61. 10. Clark, E. D., Protein refolding for industrial processes. Curr Opin Biotechnol 2001, 12 (2), 202-7. 11. Patra, A. K.; Mukhopadhyay, R.; Mukhija, R.; Krishnan, A.; Garg, L. C.; Panda, A. K., Optimization of inclusion body solubilization and renaturation of recombinant human growth hormone from Escherichia coli. Protein Expr Purif 2000, 18 (2), 182-92. 12. Azakami, H.; Mukai, A.; Kato, A., Role of amyloid type cross beta-structure in the formation of soluble aggregate and gel in heat-induced ovalbumin. J Agric Food Chem 2005, 53 (4), 1254-7. 13. Sunde, M.; Blake, C., The structure of amyloid fibrils by electron microscopy and X-ray diffraction. Adv Protein Chem 1997, 50, 123-59. 14. Antzutkin, O. N.; Balbach, J. J.; Leapman, R. D.; Rizzo, N. W.; Reed, J.; Tycko, R., Multiple quantum solid-state NMR indicates a parallel, not antiparallel, organization of beta-sheets in Alzheimer's beta- amyloid fibrils. Proc Natl Acad Sci U S A 2000, 97 (24), 13045-50. 100 15. Petkova, A. T.; Leapman, R. D.; Guo, Z.; Yau, W. M.; Mattson, M. P.; Tycko, R., Self-propagating, molecular-level polymorphism in Alzheimer's beta-amyloid fibrils. Science 2005, 307 (5707), 262-5. 16. Paravastu, A. K.; Leapman, R. D.; Yau, W. M.; Tycko, R., Molecular structural basis for polymorphism in Alzheimer's beta-amyloid fibrils. Proc Natl Acad Sci U S A 2008, 105 (47), 18349-54. 17. Luhrs, T.; Ritter, C.; Adrian, M.; Riek-Loher, D.; Bohrmann, B.; Dobeli, H., . . . Riek, R., 3D structure of Alzheimer's amyloid-beta(1-42) fibrils. Proc Natl Acad Sci U S A 2005, 102 (48), 17342-7. 18. Petkova, A. T.; Yau, W. M.; Tycko, R., Experimental constraints on quaternary structure in Alzheimer's beta-amyloid fibrils. Biochemistry 2006, 45 (2), 498-512. 19. Wiltzius, J. J.; Sievers, S. A.; Sawaya, M. R.; Cascio, D.; Popov, D.; Riekel, C.; Eisenberg, D., Atomic structure of the cross-beta spine of islet amyloid polypeptide (amylin). Protein Sci 2008, 17 (9), 1467-74. 20. Luca, S.; Yau, W. M.; Leapman, R.; Tycko, R., Peptide conformation and supramolecular organization in amylin fibrils: constraints from solid-state NMR. Biochemistry 2007, 46 (47), 13505-22. 21. Kajava, A. V.; Aebi, U.; Steven, A. C., The parallel superpleated beta-structure as a model for amyloid fibrils of human amylin. J Mol Biol 2005, 348 (2), 247-52. 22. Kelly, J. W., Mechanisms of amyloidogenesis. Nat Struct Biol 2000, 7 (10), 824-6. 23. Griffith, J. S., Self-replication and scrapie. Nature 1967, 215 (5105), 1043-4. 24. Uratani, Y.; Asakura, S.; Imahori, K., A circular dichroism study of Salmonella flagellin: evidence for conformational change on polymerization. J Mol Biol 1972, 67 (1), 85-98. 25. Prusiner, S. B., Novel proteinaceous infectious particles cause scrapie. Science 1982, 216 (4542), 136-44. 26. Beaven, G. H.; Gratzer, W. B.; Davies, H. G., Formation and structure of gels and fibrils from glucagon. Eur J Biochem 1969, 11 (1), 37-42. 27. Hofrichter, J.; Ross, P. D.; Eaton, W. A., Kinetics and mechanism of deoxyhemoglobin S gelation: a new approach to understanding sickle cell disease. Proc Natl Acad Sci U S A 1974, 71 (12), 4864-8. 28. Jarrett, J. T.; Lansbury, P. T., Jr., Seeding "one-dimensional crystallization" of amyloid: a pathogenic mechanism in Alzheimer's disease and scrapie? Cell 1993, 73 (6), 1055-8. 29. Serio, T. R.; Cashikar, A. G.; Kowal, A. S.; Sawicki, G. J.; Moslehi, J. J.; Serpell, L., . . . Lindquist, S. L., Nucleated conformational conversion and the replication of conformational information by a prion determinant. Science 2000, 289 (5483), 1317-21. 30. Nguyen, H. D.; Hall, C. K., Kinetics of fibril formation by polyalanine peptides. J Biol Chem 2005, 280 (10), 9074-82. 101 31. Chiti, F.; Webster, P.; Taddei, N.; Clark, A.; Stefani, M.; Ramponi, G.; Dobson, C. M., Designing conditions for in vitro formation of amyloid protofilaments and fibrils. Proc Natl Acad Sci U S A 1999, 96 (7), 3590-4. 32. Sunde, M.; Serpell, L. C.; Bartlam, M.; Fraser, P. E.; Pepys, M. B.; Blake, C. C., Common core structure of amyloid fibrils by synchrotron X-ray diffraction. J Mol Biol 1997, 273 (3), 729-39. 33. Kelly, J. W., The alternative conformations of amyloidogenic proteins and their multi-step assembly pathways. Curr Opin Struct Biol 1998, 8 (1), 101-6. 34. Serpell, L. C., Alzheimer's amyloid fibrils: structure and assembly. Biochim Biophys Acta 2000, 1502 (1), 16-30. 35. Serpell, L. C.; Smith, J. M., Direct visualisation of the beta-sheet structure of synthetic Alzheimer's amyloid. J Mol Biol 2000, 299 (1), 225-31. 36. Caughey, B.; Lansbury, P. T., Protofibrils, pores, fibrils, and neurodegeneration: separating the responsible protein aggregates from the innocent bystanders. Annu Rev Neurosci 2003, 26, 267-98. 37. Hou, L.; Zagorski, M. G., Sorting out the driving forces for parallel and antiparallel alignment in the abeta peptide fibril structure. Biophys J 2004, 86 (1 Pt 1), 1-2. 38. Alder, B. J.; Wainwright, T. E., Studies in Molecular Dynamics. I. General Method. The Journal of Chemical Physics 1959, 31 (2), 459-466. 39. Rapaport, D. C., Molecular dynamics simulation of polymer chains with excluded volume. Journal of Physics A: Mathematical and General 1978, 11 (8), L213. 40. Rapaport, D. C., Molecular dynamics study of a polymer chain in solution. The Journal of Chemical Physics 1979, 71 (8), 3299-3303. 41. Rapaport, D., The event scheduling problem in molecular dynamic simulation. Journal of Computational Physics 1980, 34 (2), 184-201. 42. Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G., . . . Weiner, P., A new force field for molecular mechanical simulation of nucleic acids and proteins. Journal of the American Chemical Society 1984, 106 (3), 765-784. 43. Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A., An all atom force field for simulations of proteins and nucleic acids. Journal of Computational Chemistry 1986, 7 (2), 230-252. 44. Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M., CHARMM: a program for macromolecular energy, minimization, and dynamics calculations. Journal of computational chemistry 1983, 4 (2), 187-217. 45. Levitt, M., Molecular dynamics of native protein. I. Computer simulation of trajectories. J Mol Biol 1983, 168 (3), 595-617. 102 46. Levitt, M.; Hirshberg, M.; Sharon, R.; Daggett, V., Potential energy function and parameters for simulations of the molecular dynamics of proteins and nucleic acids in solution. Computer physics communications 1995, 91 (1), 215-231. 47. Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T., Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase- trimethoprim, a drug-receptor system. Proteins 1988, 4 (1), 31-47. 48. Taketomi, H.; Ueda, Y.; Go, N., Studies on protein folding, unfolding and fluctuations by computer simulation. I. The effect of specific amino acid sequence represented by specific inter-unit interactions. Int J Pept Protein Res 1975, 7 (6), 445-59. 49. Go, N.; Taketomi, H., Respective roles of short- and long-range interactions in protein folding. Proc Natl Acad Sci U S A 1978, 75 (2), 559-63. 50. Go, N.; Taketomi, H., Studies on protein folding, unfolding and fluctuations by computer simulation. III. Effect of short-range interactions. Int J Pept Protein Res 1979, 13 (3), 235-52. 51. Liwo, A.; Pincus, M. R.; Wawak, R. J.; Rackovsky, S.; Scheraga, H. A., Calculation of protein backbone geometry from alpha-carbon coordinates based on peptide-group dipole alignment. Protein Sci 1993, 2 (10), 1697-714. 52. Liwo, A.; Pincus, M. R.; Wawak, R. J.; Rackovsky, S.; Scheraga, H. A., Prediction of protein conformation on the basis of a search for compact structures: test on avian pancreatic polypeptide. Protein Sci 1993, 2 (10), 1715-31. 53. Liwo, A.; Oldziej, S.; Kazmierkiewicz, R.; Groth, M.; Czaplewski, C., Design of a knowledge-based force field for off-lattice simulations of protein structure. Acta Biochim Pol 1997, 44 (3), 527-47. 54. Hardin, C.; Luthey-Schulten, Z.; Wolynes, P. G., Backbone dynamics, fast folding, and secondary structure formation in helical proteins and peptides. Proteins 1999, 34 (3), 281-94. 55. Hardin, C.; Eastwood, M. P.; Prentiss, M.; Luthey-Schulten, Z.; Wolynes, P. G., Folding funnels: the key to robust protein structure prediction. J Comput Chem 2002, 23 (1), 138-46. 56. Wallqvist, A.; Ullner, M., A simplified amino acid potential for use in structure predictions of proteins. Proteins 1994, 18 (3), 267-80. 57. S. Takada, Z. L.-S., P. G. Wolynes, Folding dynamics with nonadditive forces: A simulation study of a designed helical protein and a random heteropolymer. J. Chem. Phys. 1999, 110 (23), 11616-11629. 58. Sun, S.; Thomas, P. D.; Dill, K. A., A simple protein folding algorithm using a binary code and secondary structure constraints. Protein Eng 1995, 8 (8), 769-78. 59. Sun, S.; Luo, N.; Ornstein, R. L.; Rein, R., Protein structure prediction based on statistical potential. Biophys J 1992, 62 (1), 104-6. 103 60. Sun, S., Reduced representation approach to protein tertiary structure prediction: statistical potential and simulated annealing. J Theor Biol 1995, 172 (1), 13-32. 61. Sun, S., Reduced representation model of protein structure prediction: statistical potential and genetic algorithms. Protein Sci 1993, 2 (5), 762-85. 62. Irback, A.; Sjunnesson, F.; Wallin, S., Hydrogen bonds, hydrophobicity forces and the character of the collapse transition. J Biol Phys 2001, 27 (2-3), 169-79. 63. Irback, A.; Sjunnesson, F.; Wallin, S., Three-helix-bundle protein in a Ramachandran model. Proc Natl Acad Sci U S A 2000, 97 (25), 13614-8. 64. Favrin, G.; Irback, A.; Wallin, S., Folding of a small helical protein using hydrogen bonds and hydrophobicity forces. Proteins 2002, 47 (2), 99-105. 65. Forcellino, F.; Derreumaux, P., Computer simulations aimed at structure prediction of supersecondary motifs in proteins. Proteins 2001, 45 (2), 159-66. 66. Koretke, K. K.; Luthey-Schulten, Z.; Wolynes, P. G., Self-consistently optimized energy functions for protein structure prediction by molecular dynamics. Proc Natl Acad Sci U S A 1998, 95 (6), 2932-7. 67. Gupta, P.; Hall, C. K.; Voegler, A. C., Effect of denaturant and protein concentrations upon protein refolding and aggregation: a simple lattice model. Protein Sci 1998, 7 (12), 2642-52. 68. Nguyen, H. D.; Hall, C. K., Effect of rate of chemical or thermal renaturation on refolding and aggregation of a simple lattice protein. Biotechnol Bioeng 2002, 80 (7), 823-34. 69. Yoshii, H.; Furuta, T.; Yonehara, T.; Ito, D.; Linko, Y.-Y.; Linko, P., Refolding of denatured/reduced lysozyme at high concentration with diafiltration. Bioscience, biotechnology, and biochemistry 2000, 64 (6), 1159-1165. 70. Maeda, Y.; Koga, H.; Yamada, H.; Ueda, T.; Imoto, T., Effective renaturation of reduced lysozyme by gentle removal of urea. Protein Eng 1995, 8 (2), 201-5. 71. Hevehan, D. L.; De Bernardez Clark, E., Oxidative renaturation of lysozyme at high concentrations. Biotechnology and bioengineering 1997, 54 (3), 221-230. 72. Ding, F.; Dokholyan, N. V.; Buldyrev, S. V.; Stanley, H. E.; Shakhnovich, E. I., Molecular dynamics simulation of the SH3 domain aggregation suggests a generic amyloidogenesis mechanism. J Mol Biol 2002, 324 (4), 851-7. 73. Li, L.; Darden, T. A.; Bartolotti, L.; Kominos, D.; Pedersen, L. G., An atomic model for the pleated beta-sheet structure of Abeta amyloid protofilaments. Biophys J 1999, 76 (6), 2871-8. 74. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M., . . . Kollman, P. A., A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemical Society 1995, 117 (19), 5179-5197. 104 75. George, A. R.; Howlett, D. R., Computationally derived structural models of the beta-amyloid found in Alzheimer's disease plaques and the interaction with possible aggregation inhibitors. Biopolymers 1999, 50 (7), 733-41. 76. Lakdawala, A. S.; Morgan, D. M.; Liotta, D. C.; Lynn, D. G.; Snyder, J. P., Dynamics and fluidity of amyloid fibrils: a model of fibrous protein aggregates. J Am Chem Soc 2002, 124 (51), 15150-1. 77. Morgan, D. M.; Lynn, D. G.; Lakdawala, A. S.; Snyder, J. P.; Liotta, D. C., Amyloid structure: Models and theoretical considerations in fibrous aggregates. Journal of the Chinese Chemical Society 2002, 49 (4), 459-466. 78. Gsponer, J.; Haberthur, U.; Caflisch, A., The role of side-chain interactions in the early steps of aggregation: Molecular dynamics simulations of an amyloid-forming peptide from the yeast prion Sup35. Proc Natl Acad Sci U S A 2003, 100 (9), 5154-9. 79. Mager, P. P., Molecular simulation of the primary and secondary structures of the Abeta(1-42)- peptide of Alzheimer's disease. Med Res Rev 1998, 18 (6), 403-30. 80. Klimov, D. K.; Thirumalai, D., Dissecting the assembly of Abeta16-22 amyloid peptides into antiparallel beta sheets. Structure 2003, 11 (3), 295-307. 81. Voegler Smith, A.; Hall, C. K., alpha-helix formation: discontinuous molecular dynamics on an intermediate-resolution protein model. Proteins 2001, 44 (3), 344-60. 82. Voegler Smith, A.; Hall, C. K., Bridging the gap between homopolymer and protein models: A discontinuous molecular dynamics study. The Journal of Chemical Physics 2000, 113 (20), 9331-9342. 83. Smith, A. V.; Hall, C. K., Assembly of a tetrameric alpha-helical bundle: computer simulations on an intermediate-resolution protein model. Proteins 2001, 44 (3), 376-91. 84. Smith, A. V.; Hall, C. K., Protein refolding versus aggregation: computer simulations on an intermediate-resolution protein model. J Mol Biol 2001, 312 (1), 187-202. 85. Nguyen, H. D.; Hall, C. K., Phase diagrams describing fibrillization by polyalanine peptides. Biophys J 2004, 87 (6), 4122-34. 86. Nguyen, H. D.; Hall, C. K., Spontaneous fibril formation by polyalanines; discontinuous molecular dynamics simulations. J Am Chem Soc 2006, 128 (6), 1890-901. 87. Marchut, A. J.; Hall, C. K., Effects of chain length on the aggregation of model polyglutamine peptides: molecular dynamics simulations. Proteins 2007, 66 (1), 96-109. 88. Wagoner, V. A.; Cheon, M.; Chang, I.; Hall, C. K., Computer simulation study of amyloid fibril formation by palindromic sequences in prion peptides. Proteins 2011, 79 (7), 2132-45. 89. Ding, F.; Borreguero, J. M.; Buldyrey, S. V.; Stanley, H. E.; Dokholyan, N. V., Mechanism for the alpha-helix to beta-hairpin transition. Proteins 2003, 53 (2), 220-8. 105 90. Urbanc, B.; Cruz, L.; Yun, S.; Buldyrev, S. V.; Bitan, G.; Teplow, D. B.; Stanley, H. E., In silico study of amyloid beta-protein folding and oligomerization. Proc Natl Acad Sci U S A 2004, 101 (50), 17345-50. 91. Urbanc, B.; Borreguero, J. M.; Cruz, L.; Stanley, H. E., Ab initio discrete molecular dynamics approach to protein folding and aggregation. Methods Enzymol 2006, 412, 314-38. 92. Ding, F.; Buldyrev, S. V.; Dokholyan, N. V., Folding Trp-cage to NMR resolution native structure using a coarse-grained protein model. Biophys J 2005, 88 (1), 147-55. 93. Ding, F.; Tsao, D.; Nie, H.; Dokholyan, N. V., Ab initio folding of proteins with all-atom discrete molecular dynamics. Structure 2008, 16 (7), 1010-8. 94. Brinker, A.; Pfeifer, G.; Kerner, M. J.; Naylor, D. J.; Hartl, F. U.; Hayer-Hartl, M., Dual function of protein confinement in chaperonin-assisted protein folding. Cell 2001, 107 (2), 223-33. 95. Zhou, H. X.; Dill, K. A., Stabilization of proteins in confined spaces. Biochemistry 2001, 40 (38), 11289-93. 96. Minton, A. P., Implications of macromolecular crowding for protein assembly. Curr Opin Struct Biol 2000, 10 (1), 34-9. 97. Eggers, D. K.; Valentine, J. S., Molecular confinement influences protein structure and enhances thermal protein stability. Protein Sci 2001, 10 (2), 250-61. 98. Friedel, M.; Sheeler, D. J.; Shea, J.-E., Effects of confinement and crowding on the thermodynamics and kinetics of folding of a minimalist β-barrel protein. The Journal of Chemical Physics 2003, 118 (17), 8106-8113. 99. Baumketner, A.; Jewett, A.; Shea, J. E., Effects of confinement in chaperonin assisted protein folding: rate enhancement by decreasing the roughness of the folding energy landscape. J Mol Biol 2003, 332 (3), 701-13. 100. Mittal, J.; Best, R. B., Thermodynamics and kinetics of protein folding under confinement. Proc Natl Acad Sci U S A 2008, 105 (51), 20233-8. 101. Takagi, F.; Koga, N.; Takada, S., How protein thermodynamics and folding mechanisms are altered by the chaperonin cage: molecular simulations. Proc Natl Acad Sci U S A 2003, 100 (20), 11367-72. 102. Klimov, D. K.; Newfield, D.; Thirumalai, D., Simulations of beta-hairpin folding confined to spherical pores using distributed computing. Proc Natl Acad Sci U S A 2002, 99 (12), 8019-24. 103. Javidpour, L.; Tabar, M. R.; Sahimi, M., Molecular simulation of protein dynamics in nanopores. I. Stability and folding. J Chem Phys 2008, 128 (11), 115105. 104. Behe, M. J.; Englander, S. W., Sickle hemoglobin gelation. Reaction order and critical nucleus size. Biophysical Journal 1978, 23 (1), 129-145. 106 105. Drenckhahn, D.; Pollard, T. D., Elongation of actin filaments is a diffusion-limited reaction at the barbed end and is accelerated by inert macromolecules. J Biol Chem 1986, 261 (27), 12754-8. 106. Lindner, R. A.; Ralston, G. B., Macromolecular crowding: effects on actin polymerisation. Biophys Chem 1997, 66 (1), 57-66. 107. Herzog, W.; Weber, K., Microtubule formation by pure brain tubulin in vitro. The influence of dextran and poly(ethylene glycol). Eur J Biochem 1978, 91 (1), 249-54. 108. Wilf, J.; Gladner, J. A.; Minton, A. P., Acceleration of fibrin gel formation by unrelated proteins. Thromb Res 1985, 37 (6), 681-8. 109. van den Berg, B.; Ellis, R. J.; Dobson, C. M., Effects of macromolecular crowding on protein folding and aggregation. EMBO J 1999, 18 (24), 6927-33. 110. Zheng, S.; Javidpour, L.; Shing, K. S.; Sahimi, M., Dynamics of proteins aggregation. I. Universal scaling in unbounded media. J Chem Phys 2016, 145 (13), 134306. 111. Zheng, S.; Shing, K. S.; Sahimi, M., Dynamics of proteins aggregation. II. Dynamic scaling in confined media. J Chem Phys 2018, 148 (10), 104305. 112. Nguyen, H. D.; Marchut, A. J.; Hall, C. K., Solvent effects on the conformational transition of a model polyalanine peptide. Protein Sci 2004, 13 (11), 2909-24. 113. Andersen, H. C., Molecular dynamics simulations at constant pressure and/or temperature. The Journal of Chemical Physics 1980, 72 (4), 2384-2393. 114. Javidpour, L.; Sahimi, M., Confinement in nanopores can destabilize alpha-helix folding proteins and stabilize the beta structures. J Chem Phys 2011, 135 (12), 125101. 115. Javidpour, L.; Tabar, M. R.; Sahimi, M., Molecular simulation of protein dynamics in nanopores. II. Diffusion. J Chem Phys 2009, 130 (8), 085105. 116. Regan, L.; DeGrado, W. F., Characterization of a helical protein designed from first principles. Science 1988, 241 (4868), 976-8. 117. Guo, Z.; Thirumalai, D., Kinetics and thermodynamics of folding of a de novo designed four-helix bundle protein. J Mol Biol 1996, 263 (2), 323-43. 118. Nguyen, H. D.; Hall, C. K., Molecular dynamics simulations of spontaneous fibril formation by random-coil peptides. Proc Natl Acad Sci U S A 2004, 101 (46), 16180-5. 119. Marchut, A. J.; Hall, C. K., Side-chain interactions determine amyloid formation by model polyglutamine peptides in molecular dynamics simulations. Biophys J 2006, 90 (12), 4574-84. 120. Guo, J.-t.; Hall, C. K.; Xu, Y.; Wetzel, R., Modeling Protein Aggregate Assembly and Structure. In Computational Methods for Protein Structure Prediction and Modeling: Volume 1: Basic Characterization, Xu, Y.; Xu, D.; Liang, J., Eds. Springer New York: New York, NY, 2007; pp 279-317. 107 121. Ireta, J.; Neugebauer, J.; Scheffler, M.; Rojo, A.; Galván, M., Density Functional Theory Study of the Cooperativity of Hydrogen Bonds in Finite and Infinite α-Helices. The Journal of Physical Chemistry B 2003, 107 (6), 1432-1437. 122. Sheu, S. Y.; Yang, D. Y.; Selzle, H. L.; Schlag, E. W., Energetics of hydrogen bonds in peptides. Proc Natl Acad Sci U S A 2003, 100 (22), 12683-7. 123. Knowles, T. P.; Fitzpatrick, A. W.; Meehan, S.; Mott, H. R.; Vendruscolo, M.; Dobson, C. M.; Welland, M. E., Role of intermolecular forces in defining material properties of protein nanofibrils. Science 2007, 318 (5858), 1900-3. 124. Yong, W.; Lomakin, A.; Kirkitadze, M. D.; Teplow, D. B.; Chen, S. H.; Benedek, G. B., Structure determination of micelle-like intermediates in amyloid beta -protein fibril assembly by using small angle neutron scattering. Proc Natl Acad Sci U S A 2002, 99 (1), 150-4. 125. Sabate, R.; Estelrich, J., Evidence of the existence of micelles in the fibrillogenesis of beta-amyloid peptide. J Phys Chem B 2005, 109 (21), 11027-32. 126. Yu, X.; Wang, Q.; Zheng, J., Structural determination of Abeta25-35 micelles by molecular dynamics simulations. Biophys J 2010, 99 (2), 666-74. 127. Ramachandran, G. N.; Ramakrishnan, C.; Sasisekharan, V., Stereochemistry of polypeptide chain configurations. J Mol Biol 1963, 7, 95-9. 128. Hovmoller, S.; Zhou, T.; Ohlson, T., Conformations of amino acids in proteins. Acta Crystallographica Section D 2002, 58 (5), 768-776. 129. Parchansky, V.; Kapitan, J.; Kaminsky, J.; Sebestik, J.; Bour, P., Ramachandran Plot for Alanine Dipeptide as Determined from Raman Optical Activity. J Phys Chem Lett 2013, 4 (16), 2763-8. 130. Scott, K. A.; Alonso, D. O.; Sato, S.; Fersht, A. R.; Daggett, V., Conformational entropy of alanine versus glycine in protein denatured states. Proc Natl Acad Sci U S A 2007, 104 (8), 2661-6. 131. Tsai, M. I.; Xu, Y.; Dannenberg, J. J., Ramachandran revisited. DFT energy surfaces of diastereomeric trialanine peptides in the gas phase and aqueous solution. J Phys Chem B 2009, 113 (1), 309-18. 132. V. Hornak; R. Abel; A. Okur; B. Strockbine; A. Roitberg; Simmerling, C., Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Structure, Function, and Bioinformatics 2006, 65 (3), 712-725. 133. Vicsek, T.; Family, F., Dynamic Scaling for Aggregation of Clusters. Physical Review Letters 1984, 52 (19), 1669-1672. 134. Meakin, P., Fractals, Scaling and Growth far from Equilibrium. Cambridge University Press: London, 1998. 108 135. Stauffer, D.; Aharony, A., Introduction to percolation theory. 2nd ed.; Taylor & Francis: London ; Washington, DC, 1992; p xiii, 181 p. 136. Sahimi, M., Applications of percolation theory. Taylor & Francis: London ; Bristol, PA, 1994; p xiv, 258 p. 137. Meakin, P.; Vicsek, T.; Family, F., Dynamic cluster-size distribution in cluster-cluster aggregation: Effects of cluster diffusivity. Physical Review B 1985, 31 (1), 564-569. 138. Selkoe, D. J., Folding proteins in fatal ways. Nature 2003, 426 (6968), 900-4. 139. Bucciantini, M.; Giannoni, E.; Chiti, F.; Baroni, F.; Formigli, L.; Zurdo, J., . . . Stefani, M., Inherent toxicity of aggregates implies a common mechanism for protein misfolding diseases. Nature 2002, 416 (6880), 507-11. 140. Walsh, D. M.; Klyubin, I.; Fadeeva, J. V.; Cullen, W. K.; Anwyl, R.; Wolfe, M. S., . . . Selkoe, D. J., Naturally secreted oligomers of amyloid beta protein potently inhibit hippocampal long-term potentiation in vivo. Nature 2002, 416 (6880), 535-9. 141. Kayed, R.; Head, E.; Thompson, J. L.; McIntire, T. M.; Milton, S. C.; Cotman, C. W.; Glabe, C. G., Common structure of soluble amyloid oligomers implies common mechanism of pathogenesis. Science 2003, 300 (5618), 486-9. 142. Haass, C.; Selkoe, D. J., Soluble protein oligomers in neurodegeneration: lessons from the Alzheimer's amyloid beta-peptide. Nat Rev Mol Cell Biol 2007, 8 (2), 101-12. 143. Cohen, S. I. A.; Linse, S.; Luheshi, L. M.; Hellstrand, E.; White, D. A.; Rajah, L., . . . Knowles, T. P. J., Proliferation of amyloid-β42 aggregates occurs through a secondary nucleation mechanism. Proceedings of the National Academy of Sciences 2013, 110 (24), 9758-9763. 144. Ruschak, A. M.; Miranker, A. D., Fiber-dependent amyloid formation as catalysis of an existing reaction pathway. Proc Natl Acad Sci U S A 2007, 104 (30), 12341-6. 145. Cohen, S. I.; Vendruscolo, M.; Dobson, C. M.; Knowles, T. P., Nucleated polymerization with secondary pathways. II. Determination of self-consistent solutions to growth processes described by non-linear master equations. J Chem Phys 2011, 135 (6), 065106. 146. Jeong, J. S.; Ansaloni, A.; Mezzenga, R.; Lashuel, H. A.; Dietler, G., Novel mechanistic insight into the molecular basis of amyloid polymorphism and secondary nucleation during amyloid formation. J Mol Biol 2013, 425 (10), 1765-81. 147. Arosio, P.; Cukalevski, R.; Frohm, B.; Knowles, T. P.; Linse, S., Quantification of the concentration of Abeta42 propagons during the lag phase by an amyloid chain reaction assay. J Am Chem Soc 2014, 136 (1), 219-25. 148. Cohen, S. I.; Arosio, P.; Presto, J.; Kurudenkandy, F. R.; Biverstal, H.; Dolfe, L., . . . Linse, S., A molecular chaperone breaks the catalytic cycle that generates toxic Abeta oligomers. Nat Struct Mol Biol 2015, 22 (3), 207-13. 109 149. Farah, K.; Muller-Plathe, F.; Bohm, M. C., Classical reactive molecular dynamics implementations: state of the art. Chemphyschem 2012, 13 (5), 1127-51. 150. Smoluchowski, M. V., Drei Vortrage uber Diffusion, Brownsche Bewegung und Koagulation von Kolloidteilchen. Physik. Zeit. 1916, 17, 557-585. 151. Mazwi, K.; Lee, C. T., (personal communication). 152. Takada, S.; Luthey-Schulten, Z.; Wolynes, P. G., Folding dynamics with nonadditive forces: A simulation study of a designed helical protein and a random heteropolymer. The Journal of Chemical Physics 1999, 110 (23), 11616-11629. 153. Bellemans, A.; Orban, J.; Van Belle, D., Molecular dynamics of rigid and non-rigid necklaces of hard discs. Molecular Physics 1980, 39 (3), 781-782. 154. Smith, S. W.; Hall, C. K.; Freeman, B. D., Molecular Dynamics for Polymeric Fluids Using Discontinuous Potentials. Journal of Computational Physics 1997, 134 (1), 16-30. 155. Tsao, D.; Minton, A. P.; Dokholyan, N. V., A didactic model of macromolecular crowding effects on protein folding. PLoS One 2010, 5 (8), e11936. 156. Hall, D.; Minton, A. P., Macromolecular crowding: qualitative and semiquantitative successes, quantitative challenges. Biochim Biophys Acta 2003, 1649 (2), 127-39. 157. Markovitch, O.; Agmon, N., Structure and energetics of the hydronium hydration shells. J Phys Chem A 2007, 111 (12), 2253-6. 158. Rak, M., Fractals, Scaling and Growth Far from Equilibrium, Paul Meakin, Cambridge Non‐linear Science Series 5, Cambridge University Press, Cambridge, 1998, ISBN 0‐521‐45253‐8, pp. 674, £75 (hardback). International Journal of Numerical Modelling: Electronic Networks, Devices and Fields 1999, 12 (6), 493-494. 159. Fodera, V.; Zaccone, A.; Lattuada, M.; Donald, A. M., Electrostatics controls the formation of amyloid superstructures in protein aggregation. Phys Rev Lett 2013, 111 (10), 108105. 160. Zaccone, A.; Crassous, J. J.; Ballauff, M., Colloidal gelation with variable attraction energy. J Chem Phys 2013, 138 (10), 104908. 161. Szavits-Nossan, J.; Eden, K.; Morris, R. J.; MacPhee, C. E.; Evans, M. R.; Allen, R. J., Inherent variability in the kinetics of autocatalytic protein self-assembly. Phys Rev Lett 2014, 113 (9), 098101. 162. Knowles, T. P.; White, D. A.; Abate, A. R.; Agresti, J. J.; Cohen, S. I.; Sperling, R. A., . . . Weitz, D. A., Observation of spatial propagation of amyloid assembly from single nuclei. Proc Natl Acad Sci U S A 2011, 108 (36), 14746-51. 163. Edbauer, D.; Haass, C., An amyloid-like cascade hypothesis for C9orf72 ALS/FTD. Curr Opin Neurobiol 2016, 36, 99-106. 110 164. Freibaum, B. D.; Taylor, J. P., The Role of Dipeptide Repeats in C9ORF72-Related ALS-FTD. Front Mol Neurosci 2017, 10, 35. 165. Renton, A. E.; Majounie, E.; Waite, A.; Simon-Sanchez, J.; Rollinson, S.; Gibbs, J. R., . . . Traynor, B. J., A hexanucleotide repeat expansion in C9ORF72 is the cause of chromosome 9p21-linked ALS-FTD. Neuron 2011, 72 (2), 257-68. 166. Guo, Q.; Lehmer, C.; Martinez-Sanchez, A.; Rudack, T.; Beck, F.; Hartmann, H., . . . Fernandez- Busnadiego, R., In Situ Structure of Neuronal C9orf72 Poly-GA Aggregates Reveals Proteasome Recruitment. Cell 2018, 172 (4), 696-705 e12. 167. DeJesus-Hernandez, M.; Mackenzie, I. R.; Boeve, B. F.; Boxer, A. L.; Baker, M.; Rutherford, N. J., . . . Rademakers, R., Expanded GGGGCC hexanucleotide repeat in noncoding region of C9ORF72 causes chromosome 9p-linked FTD and ALS. Neuron 2011, 72 (2), 245-56. 168. Taylor, J. P.; Brown, R. H., Jr.; Cleveland, D. W., Decoding ALS: from genes to mechanism. Nature 2016, 539 (7628), 197-206. 169. Zu, T.; Gibbens, B.; Doty, N. S.; Gomes-Pereira, M.; Huguet, A.; Stone, M. D., . . . Ranum, L. P., Non-ATG-initiated translation directed by microsatellite expansions. Proc Natl Acad Sci U S A 2011, 108 (1), 260-5. 170. Mori, K.; Arzberger, T.; Grasser, F. A.; Gijselinck, I.; May, S.; Rentzsch, K., . . . Edbauer, D., Bidirectional transcripts of the expanded C9orf72 hexanucleotide repeat are translated into aggregating dipeptide repeat proteins. Acta Neuropathol 2013, 126 (6), 881-93. 171. Zu, T.; Liu, Y.; Banez-Coronel, M.; Reid, T.; Pletnikova, O.; Lewis, J., . . . Ranum, L. P., RAN proteins and RNA foci from antisense transcripts in C9ORF72 ALS and frontotemporal dementia. Proc Natl Acad Sci U S A 2013, 110 (51), E4968-77. 172. Mori, K.; Weng, S. M.; Arzberger, T.; May, S.; Rentzsch, K.; Kremmer, E., . . . Edbauer, D., The C9orf72 GGGGCC repeat is translated into aggregating dipeptide-repeat proteins in FTLD/ALS. Science 2013, 339 (6125), 1335-8. 173. Gendron, T. F.; Bieniek, K. F.; Zhang, Y. J.; Jansen-West, K.; Ash, P. E.; Caulfield, T., . . . Petrucelli, L., Antisense transcripts of the expanded C9ORF72 hexanucleotide repeat form nuclear RNA foci and undergo repeat-associated non-ATG translation in c9FTD/ALS. Acta Neuropathol 2013, 126 (6), 829-44. 174. Ash, P. E.; Bieniek, K. F.; Gendron, T. F.; Caulfield, T.; Lin, W. L.; Dejesus-Hernandez, M., . . . Petrucelli, L., Unconventional translation of C9ORF72 GGGGCC expansion generates insoluble polypeptides specific to c9FTD/ALS. Neuron 2013, 77 (4), 639-46. 175. Zhang, Y. J.; Jansen-West, K.; Xu, Y. F.; Gendron, T. F.; Bieniek, K. F.; Lin, W. L., . . . Petrucelli, L., Aggregation-prone c9FTD/ALS poly(GA) RAN-translated proteins cause neurotoxicity by inducing ER stress. Acta Neuropathol 2014, 128 (4), 505-24. 111 176. Mackenzie, I. R.; Frick, P.; Grasser, F. A.; Gendron, T. F.; Petrucelli, L.; Cashman, N. R., . . . Neumann, M., Quantitative analysis and clinico-pathological correlations of different dipeptide repeat protein pathologies in C9ORF72 mutation carriers. Acta Neuropathol 2015, 130 (6), 845-61. 177. Zhang, Y. J.; Gendron, T. F.; Grima, J. C.; Sasaguri, H.; Jansen-West, K.; Xu, Y. F., . . . Petrucelli, L., C9ORF72 poly(GA) aggregates sequester and impair HR23 and nucleocytoplasmic transport proteins. Nat Neurosci 2016, 19 (5), 668-677. 178. May, S.; Hornburg, D.; Schludi, M. H.; Arzberger, T.; Rentzsch, K.; Schwenk, B. M., . . . Edbauer, D., C9orf72 FTLD/ALS-associated Gly-Ala dipeptide repeat proteins cause neuronal toxicity and Unc119 sequestration. Acta Neuropathol 2014, 128 (4), 485-503. 179. Lin, G.; Mao, D.; Bellen, H. J., Amyotrophic Lateral Sclerosis Pathogenesis Converges on Defects in Protein Homeostasis Associated with TDP-43 Mislocalization and Proteasome-Mediated Degradation Overload. Curr Top Dev Biol 2017, 121, 111-171. 180. Chang, Y. J.; Jeng, U. S.; Chiang, Y. L.; Hwang, I. S.; Chen, Y. R., The Glycine-Alanine Dipeptide Repeat from C9orf72 Hexanucleotide Expansions Forms Toxic Amyloids Possessing Cell-to-Cell Transmission Properties. J Biol Chem 2016, 291 (10), 4903-11. 181. Hipp, M. S.; Park, S. H.; Hartl, F. U., Proteostasis impairment in protein-misfolding and - aggregation diseases. Trends Cell Biol 2014, 24 (9), 506-14. 182. Hovmoller, S.; Zhou, T.; Ohlson, T., Conformations of amino acids in proteins. Acta Crystallogr D Biol Crystallogr 2002, 58 (Pt 5), 768-76. 183. Elliott, J. R.; Lira, C. T., Supplement. In Introductory chemical engineering thermodynamics, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 2012.

Abstract (if available)

Abstract This dissertation presents the studies of dynamics of protein aggregation in (i) unbounded media, i.e. bulk solutions, and (ii) confined environments, e.g. cylindrical tubes and spherical cavities. Computer simulation is used to investigate the aggregation problem. An intermediate‐resolution coarse‐grained model of proteins as well as a high‐resolution all‐atom model of proteins are employed together with discontinuous molecular dynamics (DMD). This combination is approved to be fast and accurate to simulate the aggregation process of multiple‐protein systems. ❧ We first investigate the aggregation of multiple polyalanine peptides in unbounded media. Different numbers of identical polyalanine peptides are simulated to mimic β‐amyloid aggregation in bulk solutions. In all cases we observe fibril‐like aggregates, which have been implicated in many neuro‐degenerative diseases, such as Alzheimer's, Parkinson's, amyotrophic lateral sclerosis, and prion. Various properties of the aggregates have been computed, including dynamic evolution of the total number of aggregates, mean aggregate size, number of the peptides that contribute to the formation of β‐sheets, number of various types of hydrogen bonds formed in the systems, radius of gyration of the aggregates, and diffusivity of the aggregates. These results show that the aggregation process follows a nucleation scenario, which is also observed in experiments. Many of the quantities follow dynamic scaling, similar to those found in the aggregation of colloidal clusters. In particular, at long times the mean aggregate size S(t) grows with time as S(t)∼tᶻ, where z is the dynamic exponent. The existence of such dynamic scaling was recently confirmed by experiments. ❧ Next, we investigate the aggregation of the same polyalanine peptides but in confined environments. The medium, as a model of a crowded cellular environment, is represented by a spherical cavity, as well as cylindrical tubes with two aspect ratios. The aggregation process leads to the formation of β‐sheets and eventually fibrils. Several important properties of the aggregation process, including dynamic evolution of the total number of aggregates, mean aggregate size, and number of peptides that contribute to the formation of β‐sheets, have been computed. We show, similar to the unconfined media, that the computed properties follow dynamic scaling, characterized by power laws. The exponents that characterize the power‐law in spherical cavities are shown to agree with those in unbounded media in the same protein density, while the exponents for the aggregation in cylindrical tubes exhibit sensitivity to the geometry of the system. The effects on aggregation from the protein size, i.e. the number of amino acids in a protein, as well as the size of the confined media, have also been studied. Similarities and differences between the aggregations in confined and unconfined media are described. ❧ Last, we use DMD combined with an all‐atom model of proteins to investigate the aggregation of poly‐glycine‐alanine in unbounded media. The simulation shows poly‐glycine‐alanine is an aggregation‐prone protein and will eventually form β‐sheet‐rich aggregates, which matches the experiment observations. Dynamic evolution of the total number of aggregates has been calculated. In the simulations, we observe two unique helical structures, β‐helix and double‐helix, frequently form in different configurations. We investigate their structures in detail and the analysis results show that such helical structures are very stable and share the characteristics of both α‐helix and β‐sheet. Identical phenomena are also observed in poly‐glycine‐arginine systems. We propose that proteins having the sequence of (GX)ₙ, where X could be any non‐glycine amino acid and n is the repeat length, would share the same folding structures of β‐helix and double‐helix, and it is the glycine in the repeat that contributes to this characteristic.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Control the function, dynamics, aggregation of proteins with light illumination

PDF

Molecular dynamics simulation study of initial protein unfolding induced by the photo-responsive surfactants, azoTAB

PDF

Dynamics of water in nanotubes: liquid below freezing point and ice-like near boiling point

PDF

Stability and folding rate of proteins and identification of their inhibitors

PDF

Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale

PDF

A molecular dynamics study of interactions between the enzyme lysozyme and the photo-responsive surfactant azobenzene trimethylammonium bromide (azoTAB)

PDF

Molecular-scale studies of mechanical phenomena at the interface between two solid surfaces: from high performance friction to superlubricity and flash heating

PDF

Molecular- and continuum-scale simulation of single- and two-phase flow of CO₂ and H₂O in mixed-layer clays and a heterogeneous sandstone

PDF

Experimental study and atomic simulation of protein adsorption

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Exploring properties of silicon-carbide nanotubes and their composites with polymers

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow

PDF

Controlling the form-dynamics-function relationship of proteins with light illumination

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

Structural studies of the IAPP membrane-mediated aggregation pathway

PDF

Understanding dynamics of cyber-physical systems: mathematical models, control algorithms and hardware incarnations

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Dynamic modeling and simulation of flapping-wing micro air vehicles

Asset Metadata

Creator Zheng, Size (author)

Core Title Molecular dynamics studies of protein aggregation in unbounded and confined media

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Chemical Engineering

Publication Date 10/08/2018

Defense Date 08/21/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag confined media,dynamic scaling,molecular dynamics,OAI-PMH Harvest,protein aggregation

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Sahimi, Muhammad (committee chair), Nakano, Aiichiro (committee member), Shing, Katherine (committee member)

Creator Email alexzheng42@gmail.com,sizezhen@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-76032

Unique identifier UC11671404

Identifier etd-ZhengSize-6795.pdf (filename),usctheses-c89-76032 (legacy record id)

Legacy Identifier etd-ZhengSize-6795.pdf

Dmrecord 76032

Document Type Dissertation

Format application/pdf (imt)

Rights Zheng, Size

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA