Fibril–Barrel Transitions in Cylindrin Amyloids - American Chemical

Jul 3, 2017 - Fibril−Barrel Transitions in Cylindrin Amyloids. Huiling Zhang,. †. Wenhui Xi,. †,¶. Ulrich H. E. Hansmann,*,¶ and Yanjie Wei*,â...
0 downloads 0 Views 2MB Size
Subscriber access provided by UNIVERSITY OF CONNECTICUT

Article

Fibril-Barrel Transitions in Cylindrin Amyloids Huiling Zhang, Wenhui Xi, Ulrich H. E. Hansmann, and Yanjie Wei J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00383 • Publication Date (Web): 03 Jul 2017 Downloaded from http://pubs.acs.org on July 6, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Fibril-Barrel Transitions in Cylindrin Amyloids Huiling Zhang,∗,†,§ Wenhui Xi,∗,‡,¶,§ Ulrich H. E. Hansmann,∗,¶ and Yanjie Wei∗,‡ †Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, 518055,China ‡Shenzhen Institutes of Advanced Technology, Chinese Academy of Science, Shenzhen, 518055,China ¶Dept. of Chemistry & Biochemistry, University of Oklahoma, Norman, OK 73019, USA §Contributed equally to this work E-mail: [email protected]; [email protected]; [email protected]; [email protected]

Abstract We introduce Replica-Exchange-with-Tunneling (RET) simulations as a tool for studies of the conversion between polymorph amyloids. For the eleven-residue amyloidforming cylindrin peptide we show that this technique allows for a more efficient sampling of the formation and interconversion between fibril-like and barrel-like assemblies. We describe a protocol for optimized analysis of RET simulations that allows us to propose a mechanism for formation and interconversion between various cylindrin assemblies. Especially, we show that an inter-chain salt bridge between residues K3 and D7 is crucial for formation of the barrel structure.

1

Introduction

A wide spectrum of illnesses, ranging from metabolic to neurodegenerative diseases, is characterized by presence of amyloid fibrils. 1–3 Examples are type 2 diabetes, cataract, Alzheimer’s 1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

or Parkinson’s, to name only a few. 4 Prerequisite for the emergence of mature fibrils is the self-assembly of monomers into a critical nucleus. In this process, smaller oligomers are formed that in the case of Alzheimer’s disease appear to be the primary neurotoxic agents responsible for memory loss. 5,6 Monomers, oligomers and amyloid fibrils exist in an equilibrium of interchanging structures characterized by polymorphism, with the most neurotoxic oligomers likely being off-pathway, 7–9 and having conformational phenotypes with distinct biochemical and cellular properties. 10 Formation of fibrils and oligomers are in such scenario competing processes. As some amyloid fibrils act as an instrument for storing and releasing of hormones, one can speculate that under non-disease conditions, peptides and proteins are either in their functional form (that is folded or naturally disordered), or stored in fibrils. 11,12 On the other hand, under disease conditions, off-pathway oligomers dominate and lead to the onset of the disorders. Such hypotheses are difficult to test as the ensemble of amyloids, and the interconversion between the various structures, can rarely be probed in experiments or on a computer. For instance, simulations suffer from the problem that the assembly of mis-folded proteins into oligomers and fibrils happens on time scales (seconds to days) that are beyond the reach of most all-atom molecular dynamics simulations. Hence, exhaustive explorations of protein aggregate landscapes are rarely possible, disallowing often the calculation of the relative weight of the various structures and the pathways connecting them. 13,14 In principle, this sampling problem can be alleviated by using replica-exchange-molecular dynamics and other generalized-ensemble techniques, 15,16 where the computational costs grow no longer exponentially with protein size but only as a power law. However, the exponent is of order four, and pre-factors can be large, effectively restricting the application of these methods to systems much smaller than desired. We have recently proposed 17 to overcome some of the shortcomings of replica exchange molecular dynamics by integrating ideas from hybrid MC/MD 18 into the replica-exchange protocol. 19–21 When combined with exchange moves between “physical” models and such relying on Go-type force fields that bias toward distinct

2

ACS Paragon Plus Environment

Page 2 of 26

Page 3 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

configurational states, this Replica Exchange with Tunneling (RET) 17 approach is especially suitable for simulating structural conversions. 22 For this reason, we use the RET approach in the present paper to study the association of mis-folded proteins into small soluble oligomer, and their conversion into insoluble fibrils, processes that are difficult to probe experimentally.

(a) fibril

(b) barrel

Figure 1: The αB - crystalline segment in (A) its non-toxic fibril and (B) its toxic barrel form. Part of the experimental challenges is the difficulty of resolving structures for the transient and polymorphic toxic oligomers. Only recently have suitable models become available. One example 9 is an eleven-residue segment of αB - crystalline that can aggregate in two forms, both shown in Figure 1. Besides as amyloid fibril with the characteristic in-register β strands (Figure 1a), the segment can also form in solution a hexamer, where out-of-register antiparallel β-strands arrange as a cylindrical barrel (Figure 1b). The Eisenberg group, which resolved this hexamer structure, has proposed that this barrel motif is shared by all (or at least many) toxic oligomers in amyloid diseases. 9 This hypothesis is supported by recent simulations 23 and by experiments of Pan et al. 24 who observed β-barrels of tetrameric Aβ(1-40) peptides with antiparallel β-turn-β motifs similar to the cylindrin assemblies. Because of the presumed universality we have chosen cylindrin to study with RET simulations 3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the assembly and interconversion between amyloid assemblies. Besides extending RET simulations to amyloids, we also go beyond our earlier work 22 by introducing a protocol for optimized analysis of RET generated data. Using this optimized RET approach we will study the transitions between barrel and fibril assemblies, and comment on the question whether these toxic oligomers are on or off pathway to fibril formation.

2

Materials and Methods

2.1

Replica Exchange with Tunneling (RET)

In its standard implementation, replica-exchange sampling 19–21 aims to enforce a random walk in temperature as a way to escape out of local minima in order to achieve faster convergence at a (low) target temperature. However, the exchange move between neighboring temperatures often leads to a proposal state that is exponentially suppressed, but where if accepted the two replica will evolve rapidly to states with similar weights as before the exchange. In replica-exchange-with-tunneling (RET) 17,22 this problem is tackled by “tunneling” ˆ at a temperature through the unfavorable “transition state”, replacing configuration A by B T1 , and B by Aˆ at the temperature T2 of the neighbor replica, where, for instance, A = (qA , vA ) denotes a state characterized by coordinates qA of all its atoms and the associated velocities vA . This RET move consists of four parts: 1. Starting with a short microcanonical molecular dynamics run, the configuration A evolves at temperature T1 to a configuration A′ = (qA′ , vA′ ) with total energy Etot = EP ot + Ekin = E1 . In a similar way, the other replica, sitting at a temperature T2 , evolves from a configuration B to a new configuration B ′ that has the same energy E2 as configuration B. 2. In the second step, the two replica are exchanged, and at the same time their velocities 4

ACS Paragon Plus Environment

Page 4 of 26

Page 5 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

rescaled such that the total energies at the two temperatures stay the same: Etot (B ′′ ) = E1 and Etot (A′′ ) = E2 . Here, A′′ = (qA′ , vA′′ ) is the configuration now sitting at T2 and B ′′ = (qB′ , vB′′ ) the new configuration at T2 , and the velocities are rescaled by

vA′′ = vA′

s

E2 − Epot (qA′ ) Ekin (vA′ )

and vB′′ = vB′

s

E1 − Epot (qB′ ) . Ekin (vB′ )

(1)

This exchange move generates a “transition state” where the rescaled velocities lead to kinetic energies that compensate for unfavorable potential energies. 3. At temperature T1 the B ′′ evolves now again by microcanonical molecular dynamics converting kinetic energy into potential energy and vice versa; the final configuration ˆ = (ˆ B qB , vˆB ) has a comparable potential energy to the configuration A and a velocity distribution typical for T1 . In a similar fashion evolves the configuration A′′ at temperature T2 to a configuration Aˆ = (ˆ qA , vˆA ) that has a velocity distribution typical for temperature T2 and a similar potential energy as the start configuration A at this temperature. ˆ are compared with the start config4. In the last step, the final configurations Aˆ and B urations A and B, and the exchange move accepted or rejected with probability

exp (−β1 (Epot (ˆ qA ) − Epot (qB )) − β2 (Epot (ˆ qB ) − Epot (qA )))

with β = 1/kB T . (2)

If rejected, the molecular dynamics simulations will continue at temperature T1 with configuration A, and at temperature T2 with configuration B. Note that the above exchange probability is an approximation. Hence, detailed balance is preserved only in the limit of sufficiently large systems and long microcanonical segments ˆ B −→ B ′ and B ′′ −→ B. ˆ For details and the correct (but difficult to A −→ A′ , A′′ −→ A; evaluate) expression for the exchange probability, see Ref. 22 . When simulating conformational transitions in systems with competing attractors one can 5

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

combine the RET move with exchange moves between “physical” models and such relying on Go-type force fields that bias toward the various attractors. One way of implementing this idea is to utilize a potential energy made out of three terms Epot = Ephys +EGo +λEλ . The first term, Ephys , is given by a “physical” energy function (such as Amber 25 or CHARMM 26 ) and a suitable solvent model. The second term EGo describes a structure-based energy function (“Go-model”) designed to bias a simulation toward a pre-selected structure. The third term Eλ describes the similarity between the two models (such as the Hamming distance between the corresponding contact maps) and depends on the system, with the parameter λ tuning how strongly the physical model is biased by the Go-term. The various replica differ in this parameter, starting at a replica where the Go-model biases the ”physical” system strongly to one attractor, followed by a ladder of replica with decreasing bias toward this attractor, until reaching a replica with λ = 0, i.e., where the physical model is not biased toward this attractor. In our case, the ladder is continued to the other side with increasing parameters λ, coupling now the ”physical” model to a second Go-model that steers the system toward the other attractor. The replica-exchange move now leads to a walk in λ-space between replica where the Go-term biases the physical models towards one attractor, and on the other side such where the bias is toward another attractor. In its most straightforward implementation of this approach data will be only analyzed for the replica where λ = 0, i.e., where the physical model is not biased by other terms. However, information from replica with λ 6= 0 can be also utilized by using re-weighting techniques. 27 The weight of a configuration sampled at a replica with λ = λi is proportional to exp(−βλi Eλ ), i.e., by re-weighting with a term proportional to exp(+βλi Eλ ) one obtains its weight at λ = 0. For combining information from all replica one should use multiple histogram re-weighting approaches such as the weighted histogram analysis method (WHAM). 28 “Feeding” of physical systems by Go-models has been previously proposed by us 29–31 and others 32 as a way to avoid the intrinsic bias in Go-model simulations, but the approach is

6

ACS Paragon Plus Environment

Page 6 of 26

Page 7 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

prone to low acceptance rates, a problem that is overcome here by the RET move as shown for various systems in Ref. 22 Note that with the above defined energy function, Eq. 1 for the rescaling of velocities takes the form

vA′′ = vA′

s

E2 − Ephys (qA′ ) − EGo (qA′ ) − λ2 Eλ (qA′ ) Ekin (vA′ )

vB′′ = vB′

s

E1 − Ephys (qB′ ) − EGo (qB′ ) − λ1 Eλ (qB′ ) , Ekin (vB′ )

(3)

and RET moves are accepted with probability 

exp −β1



(1) ∆Ephys

(i)

+

(1) ∆Ego

+

(1) λ1 ∆Eλ



(i)

− β2



(2) ∆Ephys

+

(2) ∆Ego

+

(2) λ2 ∆Eλ



(4)

(i)

qA ) − Ephys (qB ); ∆EGo and ∆Eλ are defined accordingly. where ∆Ephys = Ephys (ˆ

2.2

Implementation and Technical Details

The cylindrin sequence consists of the eleven residues KVKVLGDVIEV, and six of these peptides are needed to form the barrel structure deposited in the Protein Data Bank under identifier 3SGO. The barrel can also be formed by three tandem repeats of this sequence, i.e. by three peptides GKVKVLGDVIEVGGKVKVLGDVIEV, and the corresponding structure was also derived by Laganowsky and coworkers 9 and deposited in the Protein Data Bank under identifier 3SGR. For convenience we choose this trimer instead of the 3SGO hexamer for our investigation, and refer to the barrel structure in the manuscript as C3. A fibril model build out of KVKVLGDVIEV peptides was provided to us kindly by the Eisenberg group. Arranging two of these structures we have built a fibril model for the tandem repeat by connecting them and adding glycine at suitable positions. After minimizing the resulting configurations we then obtained the trimer structure used our investigations and denoted as F3 in this manuscript. We use the pdb2gmx module of the GROMACS package, 33 with the CHARMM27 force 7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 26

field 26 and no water selected, to generate the topology and .gro (coordinate) files for physical models of both C3 and F3; and similarly the SMOG server 34 for the corresponding Go-model files. Master topology and .gro files for the C3 and F3 models were generated by merging the corresponding topology and .gro files for the physical and Go-model. The two resulting systems have each 1689 atoms. In order to generate start configurations we start from linear configurations and heat the systems up for 5 ns to T = 800 K, and then cool them down over 2 ns to T=340 K. Less than five inter-chain hydrogen bonds are found in the resulting configurations between the three chains, each forming a random coil. Our simulations rely on an in-house modification of the GROMACS package, 33 version 4.6.5 (available from the authors), and use a potential energy function made of three terms:

Epot = Ephys + EGo + λEλ

(5)

The first term, Ephys , is given by the physical interactions between the atoms in the protein as described by CHARM27 force field with CMAP corrections 26 and an implicit solvent 35 (calculating the Born radii at every step and setting a cutoff of 1.5 for the calculation of Born radii) to approximate the interaction of the protein with the surrounding solvent. We choose as Go-term EGo the SMOG energy function, 34,36 i.e. we use the SMOG-server at http://smog-server.org to generate a set of parameters for the C3 and F3 models described above. The so-generated Go-model is an all-atom model without hydrogens and atomic charges, where the interactions are given by a van der Waals contact term with the contacts defined by a shadow map. We scale the masses in the Go-model (and the corresponding potential energy terms and forces) such that the temperature scale is the same as in the physical model. The physical model and the Go-model are coupled by a third term proposed in Ref. 37 and tested in Ref., 32 where the parameter λ describes how strongly the physical model is

8

ACS Paragon Plus Environment

Page 9 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

biased by the Go-term. This coupling energy term has the form:  1 2   (∆ (i, j))   2   B A+ S + fmax ∆(i, j) Eλ = ∆ (i, j)    B    A− S − fmax ∆(i, j) ∆ (i, j)

−ds < ∆(i, j) < ds ∆(i, j) > ds

(6)

∆(i, j) < −ds

where ∆ij = δphys (i, j) − δgo (i, j) and δ(i, j) is the distance between Cα -atoms i and j in the respective model. The two factors A and B ensure continuity of Eλ at the preselected switching distances ∆ij = ±ds, fmax is the maximum force as ∆ij → ±∞, and the parameter S controls how quickly this value is approached. In our simulations we calculated the two factors A and B from pre-set S = 1, fmax = 0 and ds = 3˚ A by

A=



fmax − ds S



dsS+1

1 1 1 and B = ( + )ds2 − ( + 1)fmax ds . 2 S S

(7)

We found it easier to implement our approach into GROMACS by using two λ = 0 replica, one with the physical system and a Go-model leading to structure F3 (replica 7), and one with the physical system and a Go-model leading to structure C3 (replica 8). Only the configurations of the physical systems are exchanged between the two λ = 0 replica. Physical and Gomodel do not interact in both λ = 0 replica, while they become increasingly tighter couples with growing values of λ. The maximum λ values were carefully chosen to allow the Go model to pull within a reasonable time frame the physical model into either the C3 or the F3 assembly. The chosen 16 λ - values are : λi = 0.022, 0.017, 0.013, 0.009, 0.006, 0.003, 0.001, 0.0, 0.0, 0.0015, 0.004, 0.008, 0.015, 0.022, 0.03, 0.04, with replica 0-6 exhibiting a bias toward F3, and replica 9-15 toward C3. Note the larger λmax value on the C3 side, reflecting our observation that assembly into the barrel form is more difficult. The microcanonical segments in each RET steps are 1ps long to allow for relaxation (“tunneling”) before and after the exchange move, i.e this segment is chosen such that the velocity distribution before and after the RET move are the same. This is necessary to avoid 9

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

introducing a bias by the RET move. We chose as time between RET moves 50 ps. All our RET simulations run for 175 × 106 steps (350ns). Note that our implementation of the RET algorithm into GROMACS requires that the thermostat temperatures differ by 0.01 K, i.e go from 340 K to 340.16 K. The temperatures are controlled in our simulations by the stochastic v-rescaling method. 38 The cut-off of van der Waals(vdW) and electrostatic interactions are 1.5nm. Bond lengths of hydrogen atoms are constrained for all residues with the LINCS algorithm 39 in the physical model where we capped the N- and C- termini with methyl groups. The integration time step in the leap-frog algorithm is 2 fs. As described above, physical and go-model are independent for λ = 0, i.e., the physical model is for this replica neither biased toward barrel nor toward fibril structures. Hence, all our results are either from measurements at λ = 0 replica, or use re-weighting to correct for the bias introduced to the physical model in a replica with λ 6= 0. At such a replica, configurations are sampled with a weight ∝ e−βλEλ . Hence, for the calculation of thermodynamic averages, measurements of replica with λ 6= 0 are re-weighted according to the Eλ term of the corresponding configuration by a factor ∝ eβλEλ as discussed earlier in the method section. When calculating free energy landscapes we found that the results were more stable when the re-weighting was done for histogram bins and contributions from a replica capped for poorly sampled bins.

3

Results and Discussions

In order to show that RET simulations are suitable tools for probing conversions between amyloids, we demonstrate first that they lead to a faster walk of replica through λ-space, i.e. between replica where the Go-term biases the physical model towards barrel-like configurations (replica 0-6) and on the other side such where the bias is toward fibril-like forms (replica 9-15). The faster this walk is, the more independent configurations will be sampled in the λ = 0 replica (7 and 8) where the physical model is not biased by any Go-model term. In

10

ACS Paragon Plus Environment

Page 10 of 26

Page 11 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(a)

(b)

Figure 2: A typical walk of one of the copies of our system in λ-spacestarting from replica 5, between λ values that bias toward the barrel-like structure of Figure 1a (replica 0-6) and such that bias toward the fibril-like form of Figure 1b (replica 9-15) as seen in simulations that use our RET move (a). For comparison we show also corresponding walks for two copies in a simulation that relies on the standard Hamilton Exchange move. One of the copies started from replica 5 (red line) and walks only between replica (0-6) that bias toward the fibril structure F3; while the other (blue line), starting from replica 14, remains in the region of replica (9-15) where the bias is toward the barrel structure C3. (b) In order to allow for a comparison we show for both methods the same time interval (0-100ns). our case a total of 43 completed walks are observed, with a typical example shown in Figure (2a). For comparison we show in Figure (2b) also the walks of two replica of the same system, using the same energy function and λ-distribution, as obtained from a simulation where the common Hamilton Replica Exchange move is used and the exchange of replica between neighboring λ values accepted with probability exp (β∆λ∆Eλ ). The faster diffusion through λ-space in RET simulations is due to their higher acceptance rates (kept at ≈ 50%) over that regular Hamilton Exchange where the exchange rate (when excluding the two central replica with λ = 0) was on average ≈ 16%. Arriving at comparable exchange rates as RET would require a much larger number of replica, as was seen by us also in previous work. 22 The faster walk through λ-space leads to a better sampling at the physical (λ = 0) replica. We have checked the convergence of the RET simulation by showing in Figure 3 the free energy as function the root-mean-square deviation with either the fibril or barrel assembly for three time intervals, from 200 ns to 300 ns, and from 300 ns to 350 ns. While

11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3: Free energy (in units of RT) as function of the root-mean-square-deviation (RMSD) with either (a) the fibril or (b) the barrel form. Measurements taken in the time interval 100 ns -200 ns (red) of the simulation are compared with that taken in the interval 200 ns -300 ns (blue) and the last 50 ns (300 ns - 350 ns: black). Note that the first 100 ns have been discarded to allow for equilibrization. subject to fluctuations, the curves resemble each other, and we conclude that with our RET approach 100 ns are sufficient to derive at an equilibrium distribution in the λ = 0 replica. Our hope is that the improved sampling efficiency of RET allows one to explore exhaustively the energy landscape of this amyloid system as needed for understanding the thermodynamics of amyloid formation and the transitions between the various polymorphic states. We show in Figure 4 the landscape as collected from the central (λ = 0) replica where the physical model is not biased by a Go-term. The landscape is projected on the root-meansquare-deviation with respect to either fibril (x-axis) or barrel (y-axis) form. Representative structures are shown for the main valleys. For comparison we show also the landscapes of replica 0 where there is a maximal bias for the fibril form, and the landscape of replica 15 which has the strongest bias toward the barrel form. The strong bias is clearly seen in the landscapes of these two replica: replica 0 has the dominant minimum at the fibril structure, and replica 15 at the barrel structure. On the other hand, the λ = 0 landscape covers both regions with fibril and such with barrel-like oligomers. However, the resolution of this landscape is limited, with less distinction between barrels (or fibril) and monomer population than experimentally expected, suggesting that our statistics is not sufficient to obtain an 12

ACS Paragon Plus Environment

Page 12 of 26

Page 13 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 4: Free energy landscape (in units of RT), as obtained from our RET simulation, projected on the root-mean-square-deviation (RMSD) with either the fibril or the barrel form. In (a) we show the landscape as obtained from the central replica where λ = 0, i.e. where the physical model is not biased by either a Go-term favoring fibril or a term favoring barrel forms. Typical structures are shown for the main valleys. In (b) we show the same landscape as derived for replica 0 which is the one where the Go-term leads to maximal bias for the fibril form. In a similar way we show in (c) the landscape for replica 15 which has the maximal bias for the barrel-like form. understanding of the switching mechanism and the nature of the transition states. Much more detail can be seen when data from all 16 replica are included by re-weighting the data of a given replica from its λ-value to the bias-free case λ = 0. The so-obtained landscape is, together with typical configurations for the main minima, shown in the Figure 5a. For illustration, we have overlaid on this landscape a typical walk of one system through the whole λ-space, i.e, we show a walk starting at the replica 0 favoring maximally fibril-like configurations up to replica 15 which is the one with the strongest bias toward barrel-like assemblies. Note that this walk does not necessarily represents a typical pathway as RET relies on an artificial dynamics for sampling. This is why we emphasize the walk in replica space in Figure 5b by the color-coding of the curve that shows how the root-mean-squaredeviation (RMSD) to the fibril-like assembly changes along this pathway (a similar plot could be drawn for the RMSD to the barrel-like assembly). As the system walks through 13

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

λ-space, it changes from a structure similar to the barrel-like assembly to one close to the in-register fibril form. On this way, the inter-chain hydrogen bonds break (see the lower inset of Figure 5b). This leads to an intermediate state where the hydrophobic core is exposed to the solvent. As a consequence, the barrel breaks up and the three tandem-repeat peptides unroll and re-arrange into the flat fibril assembly. Note that in this process the hydrogen bonds between the peptides break and form again, however the ones between the two strands of a tandem repeat do not dissolve (this can be seen by the second inset, in the upper left corner of Figure 5b), i.e., the system switches from fibril to barrel (or barrel to fibril) without that the chains dissolve into unfolded, random configurations. This chain of events suggests that the barrel configurations are on-pathway to fibril formation, which is surprising as previous studies by us 40 and others 9 rather indicated that the barrel-like oligomers are off-pathway intermediates to fibril formation. For instance, Laganowski et al. 9 used targeted molecular dynamics to study the transition from the barrellike cylindrin assembly to antiparallel fibril and found a high free energy barrier between the two structures, suggesting that the barrel-like configurations are off-pathway. While their transition pathway shares some similarity with our observation about the gradually breaking of the barrel and dissolution of out-of-register hydrogen bonds, in the target molecular dynamic runs the out-of-register hydrogen-bonds break together with barrel structure. These differences may be either because their fibril model is made of KVKVLGDVIEV peptides and not of the tandem-repeats, where the connection between the two units may lead to differences in hydrogen bonding and steric constraints, or they may result from biases introduced by the targeted molecular dynamics technique. 41 We can also not exclude the possibility that use of an implicit solvent distorts the equilibrium between folded and unfolded chains. Such imbalances have been suggested in several previously published studies 42–44 but seem to be less important in amyloid simulations that rely on CHARMM force field variants. 45 Another explanation for the discrepancy may be that RET employs an artificial kinetics, and dynamical information can therefore not be extracted directly from this walk. Hence,

14

ACS Paragon Plus Environment

Page 14 of 26

Page 15 of 26

RMSD to barrel / nm

(a)

RMSD to fibril / nm

(b) 1.6

RMSD to fibril / nm

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1.2

0.8

0.4

260

265

270

275

280

285

Time / ns

Figure 5: (a) Free energy landscape (in units of RT), as obtained from our RET simulation, projected on the root-mean-square-deviation (RMSD) with either the fibril or the barrel form. Values are calculated by re-weighting data from all 16 replica to the case of λ = 0, i.e., for the case that the physical model is not coupled to a biasing Go-term. Overlaid on the landscape is a typical pathway of a system crossing the whole λ-space between replica strongly biasing toward fibril forms to such at the other end favoring barrel forms. (b) Rootmean-square-deviation (RMSD) to the fibril-structure as calculated along this walk. The coloring marks the replica visited through this walk. The inset in the lower right quadrant shows the number of inter-chain hydrogen bonds seen in the barrel (blue) and the fibril assembly (red). The number of intra-chain hydrogen bonds along this weight is shown in the second inset, in the upper left corner.

15

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

while the described walk through λ-space is a possible pathway and is consistent with the topology of the landscape, it does not describe necessarily a common transition pathway between the two oligomers. In order to understand the transition mechanism one needs to analyze the landscape of the system in more detail. For this purpose we use here a variant of Markov state models (MSM) 46 and transition path theory(TPT). 47 MSM allows one to classify all states into different discrete groups and to identify all kinetically relevant states, while TPT allows one to calculate a transition probability matrix and therefore the transitions between these states. It has been shown earlier 48,49 that these techniques are applicable in replica-exchange and other generalizedensemble techniques. Our implementation relies on a cluster analysis, using the Ward method for hierarchical clustering in R; 50 with as input a RMSD matrix calculated by the g rms module in GROMACS from the combined configurations of all 16 trajectories, sampled every 80 ps after omitting the first 100 ns (a total of 24,000 configurations). This method allows to classify configurations in layered groups, and we found a division into six clusters optimal. The weight of each cluster was obtained by summing up the weights of all cluster members as obtained by re-weighting toward the physical case of λ = 0. We then calculated the transition rates between these clusters as seen in the trajectories of the 16 replica. Both quantities are shown in Figure (6a) as an overlay on the free energy landscape with a cluster represented by a circle around the central member, with the radius of that circle given by the logarithm of cluster size. The logarithm of the measured transition rates described the transition probability between clusters and is indicated in the figure by the thickness of the arrows. Such an analysis is possible because the equilibrium configurations are sampled after the simulation has converged. Re-weighting therefore allows one to determine the weight of cluster members, and the free energy difference between the various clusters. A similar cluster distribution for a more fine-grained grouping with 150 clusters is shown in Figure (6 b), but now the size of the colored circles indicates the total number of inter-chain hydrogen bonds that members of the cluster share with the fibril structure (red, upper left corner) or

16

ACS Paragon Plus Environment

Page 16 of 26

Page 17 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

the barrel structure (purple, upper right corner). Blue circles mark the number of intra-chain hydrogen bonds shared with both structures (lower left corner), and grey circles the number of inter-chain salt bridges K3-D7 that members of a cluster share with the barrel structure (lower right corner). The hydrogen bond numbers are calculated with the g hbond tool in GROMACS.

Figure 6: (a) Overlay of the cluster distribution on the free energy landscape of Figure (5). The clusters are marked by a circle around the central member, with the radius of that circle given by the logarithm of cluster size. The arrows mark the transition between clusters, with the thickness of the arrows given by the logarithm of the transition rates. b) Four other cluster distributions projected on the same coordinates as the free energy landscape, only now the radius of the circles is given by the number of inter-chain bonds that are also seen in the fibril structure (upper left corner); by the number of inter-chain contacts also seen in the cylindrin barrel (upper right corner); by the number of intra-chain contacts (lower left corner); or by the number of salt chain contacts that are also seen in the cylindrin barrel (lower right corner).

Based on the cluster distributions in Figure (6) we identify five regions in the landscape, marked also in Figure (7a). The first region is that of un-structured unfolded and not associated chains. In the second region are the tandem-repeats loosely folded into their Ushaped structure with the intra-chain hydrogen bonds partially formed. In the third region, the repeats have proper U-shape and all intra-chain hydrogen bonds formed. However, there are still no inter-chain hydrogen bonds formed and the three peptides have not yet 17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

associated. This is the key transition region, and out of it evolve by association of the three peptides either the fibril-like structures in the fourth region of the landscape, or the barrellike structures of the fifth region. Taking into account the transition rates shown in Figure (6a) we propose the following conversion and aggregation mechanism, shown in Figure (7b). As monomers the tandem-repeat peptides prefer disordered, coil states, but can transiently also assume hydrogen-bonded U-shaped configurations. The lifetime of these configurations depends on the number of formed intra-chain hydrogen bonds, and with increasing number of these bonds stable β-hairpins are formed. Such configurations constitute the key transition state on the pathway to both fibril and barrel assemblies. Crucial for the association into these two aggregates is not only the formation of stabilizing inter-chain hydrogen bonds but also that of the salt bridge K3-D7 characteristic for the barrel-like assembly. Formation of this salt bridge seems to guide the chains into the barrel form instead of allowing them to assemble into fibrils.

Figure 7: Proposed mechanism for formation and interconversion of fibril and barrel-like structures of cylindrin assemblies. The key-regions are also marked in the free energy landscape that is shown here again.

18

ACS Paragon Plus Environment

Page 18 of 26

Page 19 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

4

Conclusions

We have shown that replica exchange with tunneling (RET) simulations allow for an efficient simulation of formation and interconversion between fibril-like and barrel-like assemblies of the amyloid-forming cylindrin peptide. Our simulations profit from two factors. First, the RET move leads to faster walk between replica where the system is biased toward fibril assemblies and such where it is biased toward barrel-like aggregates. The net-effect is a more effective sampling of independent configurations at the replica where λ = 0, i.e., where the physical model is not biased. The second advantage is that our implementation allows for inclusion of information from all replica. Hence, while at replica with λ 6= 0 the physical model is biased toward either fibril or barrel structure, this bias is accounted for and corrected through reweighting to the λ = 0 replica. Both effects allowed a detailed exploration of the free energy landscape of cylindrin assemblies, which let us propose the mechanism for formation and interconversion of the various assemblies that is sketched in Figure 7b. Its main element is that the transition between the two polymorphs does not involve unfolding of the chains but only their dissociation and re-association. Crucial for formation of the barrellike oligomer is the salt-bridge between K3-D7 which guides the association of the peptides into this form instead of the energetically more favorable fibril. Future studies will show whether this mechanism is common for cylindrin-like oligomers. We plan now simulations of out-of-register S-shaped Aβ1−42 peptides in either fibril-like or barrel-like oligomers to address this question.

Acknowledgement We thank Erik Alred for his help to this project as it was at an early stage. H.Z and W.X. have equally contributed to this work and should be considered both ”First Author”. U.H. thanks SIAT for kind hospitality during two month-long visits in 2015 and 2016. The simulations in this work were done using Tianhe-2 cluster at Guangzhou Na19

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tional Supercomputing Center and the SCHOONER cluster of the University of Oklahoma. We acknowledge financial support from NSF CHE-1266256, the National Key Research and Development Program of China under Grant No. 2016YFB0201305, the Science Technology and Innovation Committee of Shenzhen under grant No. JCYJ20140901003939036 and JCYJ20160331190123578, Guangdong Provincial Department of Science and Technology under grant No. 2016B090918122.

References (1) Murphy, R. M. Peptide aggregation in neurodegenerative disease. Annu. Rev. of Biomed. Eng. 2002, 4, 155–174. (2) Chiti, F.; Dobson, C. M. Protein misfolding, functional amyloid, and human disease. Annu. Rev. Biochem 2006, 75, 333–366. (3) Jucker, M.; Walker, L. C. Self-propagation of pathogenic protein aggregates in neurodegenerative diseases. Nature 2013, 501, 45–51. (4) Eisenberg, D.; Jucker, M. The Amyloid State of Proteins in Human Diseases. Cell 2012, 148, 1188–1203. (5) Kirkitadze, M. D.; Bitan, G.; Teplow, D. B. Paradigm shifts in Alzheimer’s disease and other neurodegenerative disorders: The emerging role of oligomeric assemblies. J. Neurosci. Res. 2002, 69, 567–577. (6) Lesne, S.; Koh, M. T.; Kotilinek, L.; Kayed, R.; Glabe, C. G.; Yang, A.; Gallagher, M.; Ashe, K. H. A specific amyloid-beta protein assembly in the brain impairs memory. Nature 2006, 440, 352–357. (7) Necula, M.; Kayed, R.; Milton, S.; Glabe, C. G. Small molecule inhibitors of aggregation

20

ACS Paragon Plus Environment

Page 20 of 26

Page 21 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

indicate that amyloid beta oligomerization and fibrillization pathways are independent and distinct. J. Biol. Chem. 2007, 282, 10311–10324. (8) Rangachari, V.; Moore, B. D.; Reed, D. K.; Sonoda, L. K.; Bridges, A. W.; Conboy, E.; Hartigan, D.; Rosenberry, T. L. Amyloid-beta(1-42) rapidly forms protofibrils and oligomers by distinct pathways in low concentrations of sodium dodecylsulfatet. Biochemistry 2007, 46, 12451–12462. (9) Laganowsky, A.; Liu, C.; Sawaya, M. R.; Whitelegge, J. P.; Park, J.; Zhao, M. L.; Pensalfini, A.; Soriaga, A. B.; Landau, M.; Teng, P. K.; Cascio, D.; Glabe, C.; Eisenberg, D. Atomic View of a Toxic Amyloid Small Oligomer. Science 2012, 335, 1228–1231. (10) Petkova, A. T.; Leapman, R. D.; Guo, Z. H.; Yau, W. M.; Mattson, M. P.; Tycko, R. Self-propagating, molecular-level polymorphism in Alzheimer’s beta-amyloid fibrils. Science 2005, 307, 262–265. (11) Maji, S. K.; Perrin, M. H.; Sawaya, M. R.; Jessberger, S.; Vadodaria, K.; Rissman, R. A.; Singru, P. S.; Nilsson, K. P. R.; Simon, R.; Schubert, D.; Eisenberg, D.; Rivier, J.; Sawchenko, P.; Vale, W.; Riek, R. Functional Amyloids As Natural Storage of Peptide Hormones in Pituitary Secretory Granules. Science 2009, 325, 328–332. (12) Maji, S. K.; Schubert, D.; Rivier, C.; Lee, S.; Rivier, J. E.; Riek, R. Amyloid as a depot for the formulation of long-acting drugs. Plos Biology 2008, 6, 240–252. (13) Anand, P.; Hansmann, U. H. Internal and environmental effects on folding and dimerisation of Alzheimer’s -amyloid peptide. Mol. Sim. 2011, 37, 440–448. (14) Berhanu, W. M.; Hansmann, U. H. In Biomolecular Modelling and Simulations; Karabencheva-Christova, T., Ed.; Advances in Protein Chemistry and Structural Biology; Academic Press, 2014; Vol. 96; pp 113 – 141.

21

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(15) Hansmann, U. H.; Okamoto, Y. Prediction of Peptide Conformation by Multicanonical Algorithm: A New Approach to the Multiple-Minima Problem. J. Comput. Chem. 1993, 14, 1333–1338. (16) Yasar, F.; Jiang, P.; Hansmann, U. H. E. Multicanonical molecular dynamics simulations of the N-terminal domain of protein L9. Euro. Physical. Lett. 2014, 105, 30008. (17) Yasar, F.; Bernhardt, N. A.; Hansmann, U. H. E. Replica-exchange-with-tunneling for fast exploration of protein landscapes. J. Chem. Phys. 2015, 143, 224102. (18) Duane, S.; Kennedy, A. D.; Pendleton, B. J.; Roweth, D. Hybrid Monte Carlo. Phys. Lett. B 1987, 195, 216–222. (19) Geyer, G.; Thompson, E. Annealing Markov Chain Monte Carlo with Applications to Ancestral Inference. J. Amer. Statist. Assoc. 1995, 90, 909–920. (20) Hukushima, K.; Nemoto, K. Exchange Monte Carlo method and application to spin glass simulations. J. Phys. Soc. (Japan) 1996, 65, 1604–1608. (21) Hansmann, U. H. E. Parallel tempering algorithm for conformational studies of biological molecules. Chem. Phys. Lett. 1997, 281, 140–150. (22) Bernhardt, N. A.; Xi, W. H.; Wang, W.; Hansmann, U. H. E. Simulating Protein Fold Switching by Replica Exchange with Tunneling. J. Chem. Theory Comput. 2016, 12, 5656–5666. (23) Do, T. D.; LaPointe, N. E.; Nelson, R.; Krotee, P.; Hayden, E. Y.; Ulrich, B.; Quan, S.; Feinstein, S. C.; Teplow, D. B.; Eisenberg, D. Amyloid β-protein C-terminal fragments: formation of cylindrins and -barrels. J. Am. Chem. Soc 2016, 138, 549–557. (24) Pan, J. X.; Han, J.; Borchers, C. H.; Konermann, L. Structure and Dynamics of Small Soluble A beta(1-40) Oligomers Studied by Top-Down Hydrogen Exchange Mass Spectrometry. Biochemistry 2012, 51, 3694–3703. 22

ACS Paragon Plus Environment

Page 22 of 26

Page 23 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(25) Salomon-Ferrer, R.; Case, D. A.; Walker, R. C. An overview of the Amber biomolecular simulation package. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2013, 3, 198–210. (26) Foloppe, N.; MacKerell, A. D. All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data. J. Comput. Chem. 2000, 21, 86–104. (27) Ferrenberg, A. M.; Swendsen, R. H. New Monte-Carlo Technique for Studying PhaseTransitions. Phys. Rev. Lett. 1988, 61, 2635–2638. (28) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules .1. The Method. J. Comput. Chem. 1992, 13, 1011–1021. (29) Meinke, J. H.; Hansmann, U. H. E. Protein simulations combining an all-atom force field with a Go term. J. Phys.:Cond. Mat. 2007, 19, 285215. (30) Berhanu, W. M.; Jiang, P.; Hansmann, U. H. E. Folding and association of a homotetrameric protein complex in an all-atom Go model. Phys. Rev. E 2013, 87, 014701. (31) Jiang, P.; Hansmann, U. H. E. Modeling Structural Flexibility of Proteins with GoModels. J. Chem. Theory Comput. 2012, 8, 2127–2133. (32) Zhang, W.; Chen, J. Accelerate Sampling in Atomistic Energy Landscapes Using Topology-Based Coarse-Grained Models. J. Chem. Theory Comput. 2014, 10, 918– 923. (33) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435–447.

23

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(34) Noel, J. K.; Whitford, P. C.; Sanbonmatsu, K. Y.; Onuchic, J. N. SMOG@ctbp: simplified deployment of structure-based models in GROMACS. Nucl. Acids Res. 2010, 38, W657–W661. (35) Onufriev, A.; Bashford, D.; Case, D. A. Exploring protein native states and largescale conformational changes with a modified generalized born model. Proteins: Struct. Funct. Bioinf. 2004, 55, 383–394. (36) Whitford, P. C.; Noel, J. K.; Gosavi, S.; Schug, A.; Sanbonmatsu, K. Y.; Onuchic, J. N. An all-atom structure-based potential for proteins: Bridging minimal models with allatom empirical forcefields. Proteins: Struct. Funct. Bioinf. 2009, 75, 430–441. (37) Moritsugu, K.; Terada, T.; Kidera, A. Scalable free energy calculation of proteins via multiscale essential sampling. J. Chem. Phys. 2010, 133, 224105. (38) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (39) Hess, B. P-LINCS: A parallel linear constraint solver for molecular simulation. J. Chem. Theory Comput. 2008, 4, 116–122. (40) Berhanu, W. M.; Hansmann, U. H. E. The stability of cylindrin -barrel amyloid oligomer modelsA molecular dynamics study. Proteins: Struct. Funct. Bioinf. 2013, 81, 1542– 1555. (41) Roe, D. R.; Bergonzo, C.; Cheatham, T. E. Evaluation of Enhanced Sampling Provided by Accelerated Molecular Dynamics with Hamiltonian Replica Exchange Methods. J. Phys. Chem. B 2014, 118, 3543–3552. (42) Akin-Ojo, O.; Wang, F. The quest for the best nonpolarizable water model from the adaptive force matching method. J. Comput. Chem. 2011, 32, 453–62.

24

ACS Paragon Plus Environment

Page 24 of 26

Page 25 of 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(43) Ravikumar, K. M.; Huang, W.; Yang, S. Coarse-grained simulations of protein-protein association: an energy landscape perspective. Biophys. J. 2012, 103, 837–45. (44) Yang, S.; Ravikumar, K. M.; Levine, H. Energy evaluation of beta-strand packing in a fibril-forming SH3 domain. J. Phys. Chem. B 2013, 117, 13051–7. (45) Auer, S.; Dobson, C. M.; Vendruscolo, M.; Maritan, A. Self-templated nucleation in peptide and protein aggregation. Phys. Rev. Lett. 2008, 101, 258101. (46) Shukla, D.; Hernndez, C. X.; Weber, J. K.; Pande, V. S. Markov State Models Provide Insights into Dynamic Modulation of Protein Function. Acc. Chem. Res. 2015, 48, 414–422. (47) E, W.; Vanden-Eijnden, E. Transition-path theory and path-finding algorithms for the study of rare events. Ann. Rev. Phys. Chem. 2010, 61, 391–420. (48) Kelley, N. W.; Vishal, V.; Krafft, G. A.; Pande, V. S. Simulating oligomerization at experimental concentrations and long timescales: A Markov state model approach. J. Chem. Phys. 2008, 129, 214707. (49) Wu, H.; Paul, F.; Wehmeyer, C.; Noe, F. Multiensemble Markov models of molecular thermodynamics and kinetics. Proc. Natl. Acad. Sci. (U S A) 2016, 113, E3221–30. (50) R Core Team, R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing: Vienna, Austria, 2014.

25

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry

26

ACS Paragon Plus Environment

Page 26 of 26