General Purpose Water Model Can Improve Atomistic Simulations of

Subscriber access provided by La Trobe University Library

Biomolecular Systems

A general purpose water model can improve atomistic simulations of intrinsically disordered proteins Parviz S. Shabane, Saeed Izadi, and Alexey V Onufriev J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b01123 • Publication Date (Web): 13 Mar 2019 Downloaded from http://pubs.acs.org on March 16, 2019

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 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 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.

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 49 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

Computed free energy landscape of an intrinsically disordered protein in OPC water model. Protein conformations with substantially different degrees of compactness are within 2 kT of each other. Conformations with very different secondary structure content co-exist within 1 kT of the global free energy minimum. 790x719mm (72 x 72 DPI)

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

A general purpose water model can improve atomistic simulations of intrinsically disordered proteins Parviz Seifpanahi Shabane,† Saeed Izadi,‡ and Alexey V. Onufriev∗,¶ †Department of Physics, Virginia Tech, Blacksburg, Virginia 24060, United States ‡Department of Biomedical Engineering and Mechanics, Virginia Tech, Blacksburg, Virginia 24060, United States ¶Departments of Computer Science, Virginia Tech, Blacksburg, Virginia 24060, United States §Center for Soft Matter and Biological Physics, Virginia Tech, Blacksburg, VA 24061, USA E-mail: [email protected] Abstract Unconstrained atomistic simulations of intrinsically disordered proteins and peptides (IDP) remain a challenge: widely used, “general purpose” water models tend to favor overly compact structures relative to experiment. Here we have performed a total of 93 µs of unrestrained MD simulations to explore, in the context of IDPs, a recently developed “general-purpose” 4-point rigid water model OPC, which describes liquid state of water close to experiment. We demonstrate that OPC, together with a popular AMBER force field ff99SB, offers a noticeable improvement over TIP3P in producing more realistic structural ensembles of three common IDPs benchmarks: 55residue apo N-terminal zinc-binding domain of HIV-1 integrase (“protein IN”), amyloid β-peptide (Aβ42 ) (residues 1-42), and 26-reside H4 histone tail. As a negative control,

1


Page 2 of 49

Page 3 of 49 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


computed folding profile of a regular globular mini-protein (CLN025) in OPC water is in appreciably better agreement with experiment than that obtained in TIP3P, which tends to over-stabilize the compact native state relative to the extended conformations. We employed Aβ42 peptide to investigate possible influence of the solvent box size on simulation outcomes. We advocate a cautious approach for simulations of IDPs: we suggest that the solvent box size should be at least four times the radius of gyration of the random coil corresponding to the IDP. The computed free energy landscape of protein IN in OPC resembles a shallow “tub” – conformations with substantially different degrees of compactness are within 2 kB T of each other. Conformations with very different secondary structure content coexist within 1 kB T of the global free energy minimum. States with higher free energy tend to have less secondary structure. Computed low helical content of the protein has virtually no correlation with its degree of compactness, which calls into question the possibility of using the helicity as a metric for assessing performance of water models for IDPs, when the helicity is low. Predicted radius of gyration (Rg ) of H4 histone tail in OPC water falls in-between that of a typical globular protein and a fully denatured protein of the same size; the predicted Rg is consistent with two independent predictions. In contrast, H4 tail in TIP3P water is as compact as the corresponding globular protein. The computed free energy landscape of H4 tail in OPC is relatively flat over a significant range of compactness, which, we argue, is consistent with its biological function as facilitator of inter-nucleosome interactions.

1

INTRODUCTION

Intrinsically disordered proteins (IDPs) and intrinsically disordered regions (IDRs) are characterized as proteins or parts (regions) of protein that lack stable secondary and tertiary structures under specific physiological conditions. IDPs are abundant in nature; disordered parts happen in 33.0% of eukaryotic proteins 1 . Most IDPs are involved in various biological 2



functions 2 such as molecular recognition 3,4 , molecular assembly (and amyloidogenesis), cell cycle regulation 4 , signal transduction, and transcription 1,5 . As such, IDPs are implicated in the development of several pathological conditions, including cancer, cardiovascular diseases, neurodegenerative diseases such as Parkinson’s, Huntingtons and Alzheimers disease, and diabetes 6 ; IDPs have been shown to be a promising target for structure-based drug development. However, progress in exploring IDPs and their structure-function relationships is hampered by at least two major obstacles. On one hand, experimental methods such as X-ray crystallography and nuclear magnetic resonance (NMR), which routinely deliver fully atomistic structures of regular globular proteins, can only provide ensemble averages over a large number of highly diverse and disordered states of IDPs 7,8 . In this respect, atomistic molecular dynamics (MD) simulations could offer a unique advantage as they can provide spatial and temporal resolution unavailable from experiment. These simulations can rationalize and guide experiment. In particular, they can generate free energy landscapes, which are particularly useful for IDPs, represented by multiple inter-converting conformations 5,9–13 . Yet, despite the strong promise, faithful unconstrained atomistic simulations of IDPs remain a significant challenge 14 . In part, this is because current gas-phase force fields were historically parameterized to describe native states of regular globular (folded) proteins, and much less attention was given to IDPs in the past. Poor performance for IDPs of force fields originally developed for structured proteins, in various combinations with widely used, “standard” water models is now well-documented 10,14–21 . The problem is that, in general, compared to globular proteins in their native states, atomistic modeling of IDPs is inherently more demanding: these proteins are represented by multiple inter converting conformations with substantial populations. Thus, while a simulation that focuses on the unique native state of a globular protein can be fairly robust to errors in the force-field that over-stabilize the native state, the same errors of just 1 or 2 kB T to reality may lead to a completely wrong relative abundance of substantially different conformations representing the IDP. The

3


Page 4 of 49

Page 5 of 49 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


nature of the problem is, in this respect, similar to the long standing problem with allatom simulations of single-stranded RNA, where notable recent advances 22,23 required careful balancing of the gas-phase and, importantly, solvent parts of the potential. That the water model is critical for accurate atomistic simulations of IDPs was conclusively demonstrated recently 14,15 : popular water models such as TIP3P 24 , TIP4P 24 , and TIP4P-EW 25 , used in combination with several widely accepted force-fields, all lead to overly compact protein conformations, while a special-purpose, “more cohesive”, water model designed to address this very issue 15 offered a noteworthy improvement. This finding is consistent with an earlier report 19 that mimicking a better solvent (in the polymer physics sense 26 ) by strengthening of protein-water interactions was sufficient to recover correct dimensions of several intrinsically disordered proteins. There is no doubt that availability of viable strategies 15,19 towards achieving the goal of accurate unconstrained atomistic simulations of IDPs is significant. For IDPs, recently developed 16,17 special-purpose atomistic force-fields offer noteworthy improvements. However, while noticeably better agreement with experiment is seen for some IDPs, it is lacking for others 17 , thus offering only a partial solution. Deficiencies of water models can be mitigated to various degrees by error cancellation between solute-solute, solute-solvent and solvent-solvent interactions. However, there are fundamental limitations to the accuracy of the resulting total energy function 27 .

One possible route to improving

atomistic simulations of IDPs is the development of specialized water models, e.g. designed 15 to be more “cohesive” than existing models, thus favoring less compact protein conformations. However, the “special purpose water model” approach has at least two potential drawbacks. First, as is the case with any specialized model, its performance outside of the necessary and very limited training/test set 28 may deteriorate even for the types of structures the model is intended for. For example, a recent attempt to re-parametrize TIP3P 17 let to inconsistent results with respect to correcting the over compactness of two IDPs caused by the original model; the need for a better water model was emphasized. Perhaps more trou-

4



blesome is the likelihood that a specialized model may fail outside of the type of structures for which the model is designed for. For example, a water model that works well for several IDPs may result in under stabilized globular proteins 15 that may fail to fold 29 . Being able to simultaneously describe both the compact and extended conformations is important – many proteins of interest have both the disordered and compact regions. Finally, in many cases one is interested in a combination of structure types, e.g. proteins and DNA simultaneously. Generally speaking, a water model designed to work well for a certain class of biomolecules is likely to have its liquid state properties deviate significantly from experiment compared to a model designed to describe the liquid state as best as possible. Unintended consequences of such deviations in biomolecular simulations are hard to predict or control. In light of the above caveats, here we explore whether a general purpose water model can offer an improved atomistic simulations description of IDPs by addressing the same problem discussed above – the over stabilization and the compact protein conformations by common water models. In particular, a number of highly accurate general purpose n-point rigid models were developed recently, see e.g. Ref. 30 for a review. The question is whether the higher accuracy relative to older popular models mentioned above can automatically deliver a better performance for IDPs. Specifically, here we test OPC 31 – a 4-point model that represents a global accuracy optimum in the space of the lowest order electrostatic multipole moments of water molecule. For 11 key liquid state properties against which water models are often benchmarked 25,32,33 , OPC is, on average, within 0.76% of the experiment. The model also shows improved performance in several very different types of atomistic simulations 22,34–36 , including those where a delicate relative balance between RNA conformations was significantly improved 22,23 . Encouraged by the improvements that OPC offers for a common IDPs benchmark protein and peptide (apo N-terminal zinc-binding domain of HIV-1 integrase and amyloid β-peptide (Aβ42 )), we have used the model to simulate conformational ensembles of an important example of IDP – H4 N-terminal histone tail. This (26-residue ) polypeptide has an im-

5


Page 6 of 49

Page 7 of 49 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


portant role in chromatin compaction. Pioneering computational studies 37–42 have already demonstrated strong potential of atomistic simulations for investigating structure-function relationships in histone tails, including the role of post-translational modifications (PTM), such as lysine acetylation, in modulating conformational properties of the tails. Subtle changes in these conformational properties can affect compaction of nucleosome arrays 43 . Atomistic simulations can be instrumental in identifying general trends in the role of histone tails in epigenetic mechanisms. For example, a recent work 38 suggested that the effect of combined acetylations of H4-histone tail are more analogues to a rheostat rather than to a ”binary code”: how many of the sites are acetylated may be more important than which specific ones. The use of OPC water model may offer improved description of simulated conformational ensembles of H4 tail, compared to those based on TIP3P used in most previous studies. Here we explore over-all compaction and free energy landscape of H4-histone, with an eye towards making connection to its biological function.

2

Methods

2.1

Test Cases and Structure Preparation

Three intrinsically disordered proteins and a regular globular protein were used. IDPs: (1) apo N-terminal zinc-binding domain of HIV-1 integrase (called “protein IN” here) (2) amyloid β-peptide (Aβ42 ) and (3) human H4-histone tail. The CLN025 a globular protein is a fast-folding mini protein. Protein IN is one of the few completely disordered proteins that have been studied extensively experimentally 44,45 and by atomistic dynamics simulations 15,17 . To check convergence of the simulations we chose the largest IDP from the three tested here, which is proten IN. For this protein, two different starting conformations were selected: a fully extended state and a folded state. The folded state was selected from one of the two monomers of the first structure in the original PDB file (1WJB 46 ). The fully extended state of IN was produced 6



by heating up the folded state of the protein to 400 K using the implicit solvent protocol described below. Aamyloid β protein is a well-studied IDP 47–54 ; its aggregation and accumulation in the brain tissue is associated with Alzheimer’s disease. Here the initial structures of Aβ42 and H4-histone tail were generated using sequence command in AMBER tools. The following sequence 55 was used for Aβ42 : DAEFRHDSGYEVHHQKLVFFAEDVGSNKGAIIGLMVGGVVIA. For H4-histone tail the sequence used was: SGRGKGGKGLGKGGAKRHRKVLRDNI. The structure created by sequence command is fully extended (essentially a straight line), using this structure for building an explicit solvent water box around it would create an unnecessary large box, too expensive to simulate. To circumvent the problem, the initial structure of Aβ42 and H4 were first simulated for 1 µs in implicit solvent at 300 K to obtain more realistic, compact structures, which were subsequently used as starting points for the explicit solvent simulations. The folded protein studied here is CLN025 56 (PDB entry 2RVD), which is a fast-folding 166 atom mini-protein frequently used in computational studies of protein folding 27,57–59 .

2.2

Simulation Protocols

The initial structures of protein IN, Aβ42 and CLN025 were protonated at pH of 6.5 and H4-histone tail at pH 7.5 using H++ server 60 . AMBER ff99SB force field was used for both implicit and explicit solvent simulations, unless otherwise noted. (We have also performed relatively short test simulations of protein IN using ff12SB+OPC. We did not explore this combination further because the results did not appear as promising compared with ff99SB+OPC combination.)

The GPU implementation of Amber15 61 was used for all the

MD simulations, unless otherwise stated. The simulation speed for various structures and water models is reported in the Supplemental information.

7


Page 8 of 49

Page 9 of 49 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


2.2.1

Explicit Solvent

For explicit solvent simulations, the initial structures of protein IN were solvated in a cubic and rectangular box of water with a minimum solute to wall distance of 4 ˚ A. The box size for the folded and extended structures were 76 ˚ A × 76 ˚ A × 76 ˚ A, and 92 ˚ A × 70 ˚ A × 70 ˚ A, respectively. Two water models were used for the explicit solvent simulations: the recently developed 4-point OPC and the commonly used TIP3P. The initial structure of Aβ42 was solvated in two cubic boxes of OPC water with a minimum solute to wall distance of 15 ˚ A and 24 ˚ A respectively. The second box contained close to double the amount of solvent compared to the first, it was used to assess the effect of solvent box size on the simulated conformational ensemble. The corresponding box sizes 3

are: 67.5 ˚ A × 67.5 ˚ A × 67.5 ˚ A (volume = 307277 ˚ A ) with 32698 total number of atoms, and, ˚3 ) with 69318 total number of for the larger box: 85 ˚ A × 85 ˚ A × 85 ˚ A (volume = 616953 A atoms. All the systems were first neutralized by adding the appropriate amount of counter ions (Na+ or Cl-), then 0.1 M salt (NaCl) was added to each system except for H4-histone tail and Aβ42 where 0.15 M (NaCl) was added. The particle mesh Ewald (PME) method 62 with periodic boundary condition was used for all explicit solvent simulations. The structure was first relaxed with 2000 steps of conjugate-gradient energy minimization, 2

with solute atoms restrained to the initial structure by force constant of 5 kcal/mol/˚ A. The IDP systems were heated to 300 K over 600 ps, with a restraint force constant of 1 2

kcal/mol/˚ A . CLN025 was heated to its experimental melting temperature of 340 K. The 2

A, systems were then equilibrated for 2 ns with a restraint force constant of 0.1 kcal/mol/˚ 2 followed by another 2 ns of equilibration with a restraint force constant of 0.01 kcal/mol/˚ A.

All restraints were removed for the production stage described below. A direct space cutoff of 9 ˚ A was used for all stages of the simulations, Van Der Waals interactions outside the cutoff distance were approximated via a continuum model: vdwmeth = 1 option in AMBER (which can be mimicked by DispCorr = EnerPres in GROMACS). 8



Adding this long-range correction to the LJ potential is the default in AMBER, used in parametrizing and testing of OPC family water models 31,63 . Langevin dynamics 64 was used for temperature regulation with a collision frequency of γ = 1 ps−1 . SHAKE algorithm 65 was used to constrain covalently bound hydrogen atoms. Hydrogen Mass Repartitioning (HMR) 66 was applied only to the most time consuming simulations of IDPs (IN, Aβ42 and H4-histone tail). In HMR, the mass of each hydrogen is changed to 3.024 daltons which is the default value in AMBER. The mass of the other atom bonded to the hydrogen is modified so that the total mass of the pair remains the same. The HMR allows the use of a longer time step (4 fs) for the production step of the IDP systems and thus speeding up the simulation. The same 4 fs time step was used for water (but without HMR), see Table 12. Integration step size of 2 fs was used for the CLN025 simulations. The production simulations in explicit solvent were performed for 15 µs for each trajectory of protein IN; a previous study 15 showed that 10 µs simulation is sufficient to obtain reasonably converged estimates of the radius of gyration of IDPs, which is our main metric of water model performance here. For H4-histone tail and Aβ42 , each trajectory is 10 µs long, and for CLN025 we conducted two simulations, each 4 µs long. For the analyses presented in the Results and Discussion section, snapshots from the MD trajectory were saved every 10 ps. Default values were used for all other simulation parameters. 2.2.2

Implicit Solvent

The implicit solvent simulations were conducted using the GB Neck model (igb = 8 in AMBER). The structure was first relaxed with 1000 steps of conjugate-gradient energy minimization. The systems were heated to 300 K in H4-histone tail case and to 400 K in the protein IN extend case over 20 ps. After that, the production step without cutoff (cut = 999) was conducted. SHAKE algorithm (ntc = 2) was used to constrain covalently bound hydrogen atoms. Langevin dynamics (ntt = 3) was used for temperature regulation with a collision frequency of γln = 1 ps−1 . The time step was 2 fs for all stage of the implicit solvent

9


Page 10 of 49

Page 11 of 49 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


simulation.

2.3

Trajectory Analysis

CPPTRAJ 67 module of Amber Tools 15 was used to calculate root-mean-square deviations (RMSDs) and radius of gyration (Rg ) for all the structures. The number of contacts (NOC) between residues was calculated using an in-house script. Two residues with no covalent bonds between them are assumed to be in contact if the shortest distance between them (minimum distance between two atoms belonging to different residues) is less than 2.8 ˚ A. Since the two simulations (extended and folded) of protein IN appear reasonably converged (Figure 1), the 2D free energy landscapes were calculated using the trajectory from the simulations based on the extended structures as the starting point. The first five microseconds of the simulations were excluded from the analysis. The Rg and NOC were calculated for all of the remaining 250,000 frames of the trajectory. The free energy landscapes were cali ), where culated from bin populations in the space of Rg and NOC using ∆G = −kB T ln( NNtot

kB is Boltzmann’s constant, T is the temperature, Ni is the population of bin i and Ntot is the total number of collected samples of the trajectory. The most populated bins contained more than 10,000 data points. In Figure 5 (left), the length of bins in the Rg coordinate was 0.58 ˚ A, and in the NOC coordinate was 4. In Figure 5 (right), the number of residues in the β-sheet were calculated using the secstruct command in CPPTRAJ. The length of bins in the Rg coordinate was 0.23 ˚ A and in the NOC coordinate was 1. The same method was used for making the free energy landscape of Aβ42 (Figure(3) in the Supplemental information). The percentage of residues in conformation corresponding to each type of secondary structure, Table 8, in zones 1, 2 and 3 specified in Figure 5 were estimated as follows: first, all of the frames of the trajectory in those zones were selected. Then, by using the secstruct command in CPPTRAJ, the secondary structure for those frames of the trajectory was calculated. After that, the number of residues in each frame that had the specific secondary structure conformation was determined, and divided by the number of frames in that zone 10



× 55 (the number of residues in protein IN). For the percentage of secondary structure type in Table II of the Supplemental information was calculated as follows. First by using the secstruct command in the CPPTRAJ, the secondary structure for each trajectory (all frames) was calculated. After that, the number of residues in each frame that had the specific secondary structure was determined, and divided by the total number of frames in the corresponding trajectory × 42 (the number of residues in Aβ42 .) As an example, in a hypothetical situation where all the residues in all the frames were in β-sheet conformation, a 100 % would be reported. The error bars in Table 5 are the standard error of the mean. It was calculated as follows. The entire trajectory was divided into 10 equal intervals 1 µs each. The averages in each interval were calculated, and the standard deviation estimated from the resulting 10 values. 2.3.1

Free Energy Landscape of CLN025

The free energy landscape in Figure 4 was calculated as follows: first, the RMSD relative to the experimental folded structure was calculated for each snapshot. The calculated RMSD values were partitioned into bins of 0.2 ˚ A. The free energy for the conformations represented i ), where kB is the Boltzmann by each of these bins was then calculated as ∆G = −kB T ln( NNtot

constant, T is the temperature, Ni is the population of bin i, and Ntot is the total number of the all samples. Here we followed Ref 57 and references therein. We defined the folded and unfolded states as structures with backbone RMSD values of < 1.5 ˚ A and > 4.5 ˚ A, respectively. this is consistent with the location of the folded and unfolded basins of the free energy landscape. Note that the basin between these two states represents a mis-folded state 68 not seen experimentally 69 . ˚ The free energy change in Table 6 was calculated 57 as ∆G = −kB T ln p(RM SD>4.5A) , p(RM SD 0). Energy is defined as the difference in mean-field solvation energy between a residue within the maximally-collapsed chain and in free solution. The volume fraction (φ) and the radius of gyration (Rg ) are related to each other by φ =

(RgF CN )3 , Rg3

where RgF CN is the radius of the gyration of the fully-compact native 1/2

state (with φF CN = 1, by definition). Radius of gyration of an ideal chain,hRg2 i0 , is fixed by demanding that:

12



φ0 =

(RgF CN )3 3/2 hRg2 i0

=

p

19/27(N − 1)−1/2 ,

Page 14 of 49

(3)

Equation 3 is obtained from Landau’s theory of phase transitions 71 . In Figure 6, for the “polymer theory” estimate (green curve), we assumed that RgF CN = 8.6 ˚ A. This is the minimum value of Rg seen in all of our MD simulations of H4-histone tail in OPC and TIP3P. The number is also in agreement with the experimental data in Table 1. In the absence of an experimental estimate of ǫ for protein IN, we used the available IDP, the cold-shock protein from Thermotoga maritima, CspTm, considered in Ref 15 to benchmark a water model for IDPs. Specifically, the magnitude of ǫ in equation 2 was set to ǫ=1.4, from Figure 1 of reference 72 , (2b, purple line) at zero concentration of the denaturant. Among all the proteins considered in Ref. 72 , the chosen ǫ is the lowest. For the polymer theory fully denatured state (blue curve, Figure 6), we again considered CspTm for consistency (2b, purple line in Figure 1 in reference 72 ), now at 5[M] concentration of the denaturant, which gives ǫ = 0.75. The magnitude of Rg for a random coil as the same size as protein IN and Aβ42 in Table4 and in Table5 were calculated from Equation 3. For protein IN, we assumed that RgF CN = 10.3 ˚ A, and for Aβ42 we assumed that RgF CN = 9.3 ˚ A which are the minimum value of Rg seen in all of our MD simulations of protein IN and Aβ42 in OPC and TIP3P. The number is also in agreement with the experimental data in Table 2 and Table3 . The magnitude of Rg for the fully denatured protein and urea denatured at the same size of H4-histone tail and protein IN in Table 10 and Table 4 was calculated from probability distribution Equation 1. For H4-histone tail, N = 26 and for protein IN, N = 55 and as explained above, ǫ = 0.75.

13


Page 15 of 49 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


2.5

Radius of Gyration of Globular Proteins with the Same Number of Residues as in H4-histone Tail, protein IN and Aβ42 .

To obtain estimates of the Rg s of globular proteins of the same lengths as H4-histone Tail, protein IN, and Aβ42 , we randomly selected five proteins from PDB data bank and calculated their Rg , see Table 1, Table 2 and Table 3. Table 1: Radius of gyration of globular proteins of the same length (26 amino-acids) as H4-histone tail. PDB Entry 5kwo 5kwz 2ndd 5ua8 5u9y Average

Method NMR NMR NMR NMR NMR

Rg ( ˚ A) 8.2 8.6 9.4 12.1 12.2 10.1

Table 2: Radius of gyration of globular proteins of the same length (55 amino-acids) as protein IN. PDB Entry 1kma 1shp 2jot 2mhv 2mz0 Average


Rg ( ˚ A) 11.6 10.6 11.3 11.4 12.2 11.4

Table 3: Radius of gyration of globular proteins of the same length (42 amino-acids) as Aβ42 . PDB Entry 1h5o 1px9 1qdp 1wqk 2mub Average


14

Rg ( ˚ A) 10.5 9.7 10.5 9.6 9.6 10



3

RESULTS AND DISCUSSION

3.1

Protein IN Ensembles in OPC vs. TIP3P Water.

Radius of gyration (Rg ) of a protein is a straightforward metric that directly characterizes compactness of the molecule: the metric has been used widely to evaluate performance of force-fields for IDPs 15,17,73,74 . As we argue below, Rg is a better metric for evaluation of over-all performance of a solute model in the context of IDPs than the amount of secondary structure content when the latter is low, as is the case with protein IN. Here, we investigated ensembles of protein IN simulated in OPC and TIP3P, and compared the resulting Rg values along the trajectory to experiment, Figure 1. We note that the two trajectories starting from the folded and extended states are reasonably converged. With respect to TIP3P, our results are in agreement with previous findings 15 , in that simulated protein IN is unrealistically compact. In fact, the protein stayed compact along the entire trajectory in Figure 1. In contrast, the size of protein IN simulated in OPC fluctuated significantly, visiting both compact and extended states, which is a more realistic picture of an IDP. In fact, a number of structures seen along the OPC-based trajectory reached the value (green line in Figure 1 ) estimated 44 from a FRET experiment, in contrast to the TIP3P-based simulation, where Rg was always considerably lower. We therefore argue here that OPC yields a much more realistic ensemble for protein IN. At the same time, the average of the RG for last 5 µs of simulation of protein IN in OPC is Rg = 15.2 ˚ A, which it is still smaller than the FRET-based estimate of Rg = 23 ˚ A; we offer two possible explanations. First, it is possible that the combination of ff99SB+OPC does not quite reach the level of accuracy needed for a quantitative agreement with experiment, resulting in slight overcompaction of the simulated structures. That observation is consistent with the ∼ 1 kB T bias towards compact states of a small protein, discussed below. Note that for IDPs, with their very shallow free energy landscapes, a 1 kB T bias might account for the deviation from

15


Page 16 of 49

Page 17 of 49 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 1: Time-evolution of the radius of gyration, Rg , of protein IN obtained from all-atom simulations at 300 K, in OPC (left panel) and TIP3P (right panel) water. Each simulation starts from two different conformations: folded (red) and extended (black); these starting structures are shown at t = 0. For the simulation in OPC, representative structures seen along the trajectory are also shown. The blue line indicates an Rg of a globular protein of the same number of residue (55) as protein IN (see “Methods”). The pink line indicates Rg of a random coil corresponding to a protein of the same number of residue (55) as protein IN (see “Methods”). The green line corresponds to an Rg estimate based on FRET experiment for protein IN.

16



Page 18 of 49

FRET seen in Figure 1. At the same time, we do not rule out the possibility that the estimate of Rg = 23 ˚ A, derived 44 from end-to-end distance measured by FRET and assuming a Gaussian coil statistics, might be an overestimation of the true Rg for protein IN. It is surprising that FRET-based Rg estimate is higher than the Rg value of the corresponding random coil R = 21.3 ˚ A, and is almost equal to the Rg of a thermally denatured protein of the same size of protein IN, Table 4. Several theoretical 75–79 and experimental 80–82 studies argue that IDPs are more compact than true random coils: for protein IN, the net charge per residue, f+ − f− = 0.073 is far below the threshold (0.4 < (f+ − f− ) < 0.65) expected 83 for the transition from globule to true Gaussian coil. Thus, one may expect some degree of globular character in protein IN. Table 4: Average Rg for various states of compaction of protein IN. The value in each explicit water model was extracted from the last 5 µs of the corresponding MD simulation. Rg for the globular form is from Table 2; see ”Methods” for the other estimates. State of compaction IDP IDP Globular protein IDP Random coil Thermally denatured 84 Urea denatured

3.2

Method MD simulation With OPC MD simulation With TIP3P Experiment FRET Theory Theory/Experiment Theory

Rg 15.2 11.5 11.4 23 21.3 22.4 24.5

Aβ42 Peptide Ensembles in OPC Water.

Here we have investigated ensembles of Aβ42 peptide simulated in two boxes of OPC water of different volume, Figure 2. In contrast to protein IN, the experimental Rg of Aβ42 lies between the random coil and globular protein values, consistent with expectation. In the simulations, the size of the peptide fluctuated significantly on the microsecond time-scale,

17


Page 19 of 49 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


visiting both compact and highly extended states multiple times, suggesting that our Rg estimates are reasonably converged. . For both solvent boxes the calculated average radius of gyration of Aβ42 is within statistical error margin of the experimental value of ∼ 11 ˚ A; the fluctuations are also mainly within the experimentally observed range 52 , Table 5.

Figure 2: Time-evolution of the radius of gyration, Rg , of protein Aβ42 obtained from allatom simulations at 300K, using OPC water model. The solvent box in the right panel contains twice as many water molecules as the box in the left panel, which is “standard” size, see “Methods”. The blue line indicates Rg of a globular protein with the same number of residues (42) as Aβ42 . The pink line indicates Rg of a random coil corresponding to a protein of the same number of residues (42) as Aβ42 . The green line corresponds to the average of Rg estimate based on Ion Mobility Measurements experiment (IMM) for Aβ42 , the shaded area shows the range of Rg from the IMM experiment 52 . See “Methods” for details. Representative snapshots along the trajectory (shown) are available in PDB format as Supplementary information. Experimentally 53 , Aβ42 is known to have a relatively high percentage of β-sheet, turn and bend conformations. The computed distribution of secondary structure populations is in over-all qualitative agreement with experiment in this respect, Figure 3. The simulation also predicts a noticeable population of structures with 3 or 4 β-sheet fragments, Figure 2, known to play an important role in the aggregation and toxicity 49 of Aβ42 . We did 18



Page 20 of 49

Figure 3: (a) Secondary structure composition of Aβ42 peptide estimated from a collection of different experiments 53 . (b) Predicted secondary structure of Aβ42 in the standard box of OPC water. Trajectory average. (c) Predicted secondary structure of Aβ42 in the double volume box of OPC water. Trajectory average.

Table 5: Average Rg for various states of compaction of protein Aβ42 . For OPC, the averages correspond to the last 5 µs of the corresponding MD simulation. Rg for the globular protein is from Table3; see ”Methods” for the other estimates. State of compaction IDP IDP IDP Globular protein Random coil

0 Method simulation with OPC simulation with OPC (double volume) Experiment (IMM) 52 Experiment (Table3) Theory

19


Method Rg (˚ A) 11±0.4 11.9±0.6 11 to 15 10 18.3

Page 21 of 49 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


not pursue a detailed quantitative analysis of the secondary structure populations of Aβ42 because fully converged estimates of secondary structure of an IDP of this size (42 residues) likely require longer simulation times 15 than those reported here, see also Figures 3 and 4 in the Supplemental information. We have not observed a significant dependence of the average radius of gyration of Aβ42 on the size of the solvent box. A small increase, 0.9 ˚ A , of the average Rg upon doubling of the solvent volume is at the borderline of statistical significance, Table 5. The fluctuation range seen on the 10 microsecond time-scale is similar for both solvent box sizes used here, Figure 2, within uncertainties of a stochastic trajectory (the single short-lived excursion of Rg above the random coil reference seen at the end of the trajectory in the large volume box, while curious, is not statistically ). The predicted trajectory-average β-sheet population of Aβ42 increases by ∼ 30 % upon doubling of the solvent volume, see Table II of the Supplemental information. However, as mentioned above, it is difficult to put a statistical error bar on the secondary structure population average estimated from our 10 microsecond long simulations.

Since confinement is known to affect folding thermodynamics 85,86 of regular proteins, it is instructive to estimate the expected magnitude of the effect for Aβ42 , assuming that the same basic physics is at play for IDPs. Specifically, the relative shift in folding temperature Tf relative to the the value in bulk solvent was estimated 85 as (Tf − Tfbulk )/Tfbulk ∼ (R/RgU )−γ , where R is the radius of the confining potential, and RgU is the radius of gyration of the unfolded state of the protein. For example, for (R/RgU ) = 2, the relative change in Tf was found 85 to be ∼ 5% or less, and it dropped down to 1 % or less for (R/RgU ) = 4. It A , Table 5. is reasonable to take RgU to be the random coil RgRC ; for Aβ42 RgRC = 18.3 ˚ In the standard periodic box simulation set-up, the protein interacts with its replicas in the neighboring box images. Thus, an upper bound on the confining potential radius can be estimated as R = B − RgRC , where B is the solvent box size, which yields R/RgU = B/RgRC − 1. In our simulations of Aβ42 , see “Methods”, B/RgRC = 3.6 for the “standard”

20



box and 4.6 for the “double volume” box, which means that a change of about 5 % or less in the thermodynamics (Tf ) of Aβ42 relative to the bulk may be expected upon confining the peptide in the standard solvent box, while the effect should decrease to about 1 % or less for the “double volume” box. These estimates, along with our simulation results, suggest that the “double volume” solvent box of OPC of the size B/RgRC = 4.6 used here is a reasonably “safe” choice with respect to possible influence of the confinement effects, while B/RgRC = 3.6 may still represent an acceptable compromise between accuracy and speed, certainly with respect to over-all protein size, and for qualitative assessment of its secondary structure composition. A recent study 51 investigated the effects of confined aqueous volume on the dynamic behavior and structural properties of Aβ42 in aqueous media, using TIP3P water model, and a somewhat different simulation protocol. Significant differences in secondary structure were reported in parts of Aβ42 , depending on the confining solvent volume, which is not inconsistent with our findings in Figure3, with the caveat of the statistical uncertainty mentioned above. We note that solvent box sizes used in Ref. 51 and in this work are similar. We also note that structural properties of IDPs are likely to be sensitive to water model properties, just like protein folding landscapes are 27 . Water model properties can be affected by details of computational protocols used 30 . Taking together, our results and those of Ref. 51 call for further investigation of the effect of confining volume on IDPs. Certainly, one can not assume that a bare minimum solvent box and any water model will deliver accurate results. As a concrete rule of thumb, we suggest that the solvent box size should be at least four times the radius of gyration of the random coil corresponding to the IDP.

3.3

Negative Control: Folding of a Small Globular Protein

Thus far we have established that the combination of OPC and ff99SB is a promising choice for intrinsically disordered proteins. Intending to move beyond special-purpose water models, we need to be certain that the specific force-field/water model combination works for regular 21


Page 22 of 49

Page 23 of 49 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


proteins too. The ability of a water model to work well for regular globular proteins is also important for partially disordered proteins – proteins whose parts are disordered (IDRs). This class may be the majority among proteins that are not completely structured. One of the well-known examples of partially disordered proteins is the p53 tumor suppressor protein in which nearly 30% of its structures are disordered 87,88 . To establish the ability of OPC to describe thermodynamics of a regular protein, we have carried out a series of simulations on a small protein CLN025 at its melting temperature – a common benchmark used in computational protein folding studies 57–59,69,89–93 . Figure 4 shows the folding free energy landscapes of CLN025 computed in TIP3P and OPC water. There is little difference between the two landscapes for compact states (RM SD < 3˚ A with the X-ray reference), including the native state. The two landscapes begin to diverge from A; the extended state is markedly less populated in TIP3P compared to OPC. RM SD > 3˚ In fact, the folding free energy of this protein is significantly over predicted by TIP3P, while it is almost within kB T of the experimental reference for OPC, Table 6.

Figure 4: Free energy profile extracted from two 4 µs long all-atom simulations (∼ 30 folding/unfolding events) of mini-protein CLN025 at its experimental melting temperature (340 K), solvated in boxes of TIP3P and OPC water, see “Methods” for details. Representative folded and unfolded structures are shown as insets. 22



Table 6: Computed folding free energy (∆G = GExtended − GF olded ) for the protein CLN025 at the experimental melting temperature (340 K), solvated in boxes of OPC and TIP3P. Water Model OPC TIP3P Experiment

∆G(KB T ) 1.3 4.0 0

The fact that, relative to experiment, the balance between the compact and extended states of this globular protein in TIP3P is strongly shifted towards the compact state, and OPC corrects most of this error, is consistent with the behavior of the two water models seen for IDPs. The remaining 1.3 kB T of the imbalance towards the folded state of CLN025 in OPC is consistent with the underestimation of the radius of gyration of protein IN discussed above.

3.4

Common Features of Water Models That Show Promise for IDPs

As shown previously 15 and discussed above, widely used general purpose water models seem to be poor choices for IDPs, since the resulting structures are too compact compared to experiment. At the same time, a specialized water model TIP4P-D 15 and a general purpose model OPC (which is also a 4-point rigid model) show noticeable improvement in that respect. Here we examine the commonalities between these two water models that show promise for IDPs, and what sets them apart form the widely used general purpose models that does not work well for IDPs. We also examine differences between TIP4P-D and OPC. Key parameters and liquid state properties of several water models are presented in Table 7. The obvious common feature of OPC and TIP4P-D that sets both models apart from TIP3P and TIP4P models is considerably larger values of the LJ parameters. A convincing argument for why stronger dispersion interactions are needed for IDPs was presented previously 15 . Stronger water dispersion increases the relative weight of solute-solvent and

23


Page 24 of 49

Page 25 of 49 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


solvent-solvent terms in the fine balance between solute-solute, solute-solvent and solventsolvent interactions that ultimately determines properties of the solute. As a result, both TIP4P-D and OPC favor less compact solute states compared to the common water models in Table7. In addition, values of the dipole moment for OPC and TIP4P-D is larger, and closer to experiment, than that of any of the models listed in Table 7. Higher dipole moment is suggestive of stronger water-water electrostatic interactions, which further contributes to the shifting of the energy balance mentioned above. At the same time, parameters of TIP4P-D and OPC are not identical, e.g. the dispersion interaction in OPC is somewhat weaker than that of TIP4P-D. Since properties of the liquid state are extremely sensitive to parameters of the water model, it is not surprising that key properties of the two models are also quite different, Table 7. The differences between parameters of TIP4P-D and OPC models are just as important for correct reproduction of water properties in the liquid state of water, most relevant to biomolecular simulations. Table 7: Parameters and key liquid state properties of several n-point rigid water models, including two (in bold) that show noticeable improvement of simulation outcomes for IDPs. Shown are: LJ parameters, the dipole moment, heat of vaporization, self-diffusion coefficient, and the dielectric constant at 300 K. Water Model OPC TIP4P-D TIP3P TIP4P-EW TIP4P/2005 Experiment

12

ALJ (kcal/mol ∗ ˚ A ) 865108 904657 582000 656138 731380

6

BLJ (kcal/mol ∗ ˚ A ) 858 900 595 653 736

µ(D) 2.48 2.43 2.35 2.32 2.30 > 2.6

∆Hvap (kcal/mol) 10.57 11.3 10.26 10.58 10.89 10.52

D(10−5 cm2 s−1 ) 2.3 2.1 5.8 2.6 2.2 2.3

ǫ 78.4 68 96 63 56 78.4

In summary, water models with relatively large values of LJ parameters and of the dipole moment, TIP4P-D and OPC, show promise for IDPs because they are more “cohesive”, favoring less compact solute states. For TIP4P-D the move in this direction was by design. In contrast, OPC was designed as “globally optimal” model to represent liquid state of water only – the improvement with respect to IDPs came out automatically. The observation that OPC is more cohesive compared to commonly used models is consistent with the fact that OPC is also known to improve simulation outcomes where widely used 24



water models result in too compact solute states, e.g. in simulations of RNA structures 22 .

3.5

Free Energy Landscape of Protein IN

Since the introduction of the “folding funnel” concept 94 , free energy landscapes (FEL) of proteins have been used extensively to characterize them. Compared to a single native structure, a FEL along with representative structures on the landscape, provides a much richer characterization of a protein, including its thermodynamics 95 . These landscapes help investigate and rationalize multiple characteristics of the protein, including its stability, mechanisms of folding and molecular identification, and potential for mis-folding and aggregation. For IDPs, which are represented by multiple inter-converting structures, FELs are particularly useful. For example, a recent study 9 reported that free energy landscape of a 40-residue IDP (Aβ40 peptide) resembles an “inverted funnel” – in contrast to regular proteins, this IDP becomes more structured (secondary structure content increases) away from the free energy minimum. In general, all-atom, unrestrained simulations can generate highly detailed, unbiased FELs, but applications to IDPs have so far been limited by force-field deficiencies, that of the available water models discussed in the Introduction. Below we explore the landscape of protein IN in OPC water; based on the results above we have confidence that the generated landscape and its main feature are more realistic than what could have been produced with a commonly used water model such as TIP3P. The computed free energy landscape of protein IN is shown in Figure 5; we employ the commonly used 2D sub-space of { radius of gyration } – { number of residue-residue contacts }, which characterizes protein compaction. The landscape resembles a “shallow tub” – the region of the 2D sub-space where the structures are within 2 kB T of each other, purple and blue in Figure 5 (left panel), is relatively wide in contrast to the sharply pointed “bottom” of the classic folding funnel of globular proteins 96 . Also, in contrast to globular proteins, and in some disagreement with the findings 9 for Aβ40 peptide, the global free energy minimum of protein IN is populated with quite different 25


Page 26 of 49

Page 27 of 49 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


structures. While about half of the conformations in this region are unstructured, conformations with very different β-sheet content can also be found within 1 kB T of the global minimum in free energy, Figure 5 (right panel). The protein has very little helical content, which also changes little with the compaction (see below and the SI for a more detailed characterization of the protein secondary structure).

Figure 5: Free energy landscape of protein IN, extracted from the last 10 µs of 15 µs all-atom simulation at 300 K in OPC water. Free energy, ∆G, is calculated as − log(p(occupancy)). White squares show the selected compactness zones with the arrows pointing to representative structures in each zone. In contrast to the “inverted funnel” suggested 9 for Aβ40, for protein IN we find that higher free energy on the compaction FEL implies less secondary structure, see Table 8. Specifically, as one moves away from the bottom of the “tub”, zone 1 to zone 2, which is 2-3 kB T higher in energy than the global minimum in Figure 5 (left panel), the average β-sheet population decreases markedly, while the helical content remains almost constant and very low. The trend of decreasing secondary structure content for states with higher free energy remains in place for the transition from zone 2 to zone 3. It remains to be seen whether the trend becomes more nuanced within 1-2 kB T of the global free energy minimum. Another conclusion that can be drawn from Table 8 is that the low helical content of 26



Page 28 of 49

protein IN has virtually no correlation with its compactness – fairly compact structures in zone 1 have almost the same % of residues in a helical conformation as much less compact structures in zone 3. The observation has implications for testing of water models for IDPs: experimental helical content, when low, is a poor metric of the ability of the model to predict correct over-all compactness of the IDP. Table 8: Percentage of residues in β-sheet and helical structure conformations in different compactness zones in Figure 5. The error bars indicate standard error of the mean, see “Methods”. Zones Zone1 Zone2 Zone3

3.6

% of residues in β-sheet conformation 13.6 3.8 1.4

% of residues in a helical conformation 3.8 3.0 4.6

Application to H4-histone Tail

As an application and further testing of our OPC-based simulation protocol for IDPs, we investigate conformational ensembles of histone N-terminal tail H4 37,38,97–99 of the nucleosome. Histone tails are flexible and highly positively charged disordered peptides with multiple functions, including facilitation of inter-nucleosome interactions that promote compaction of the genomic DNA into higher order chromatin structures. X-ray structures of the nucleosome 100,101 and Magic-angle spinning NMR experiments 102 indicate high conformational flexibility of the histone tails, consistent with the specific amino acid content of tails. For example, only ∼ 15% of amino-acids in H4 are hydrophobic, and the value of the net charge per residue is high. The latter increases the disorder because of the electrostatic repulsion between like charged residues. We begin with a consistency check: comparison of the predicted average secondary structure content of H4 histone tail with experiment and previous simulations. Both OPC and TIP3P predict low helical percentage consistent with experiment, Table 9; one can probably not differentiate experimentally 54 between the low numbers seen it Table 9. However, as we 27


Page 29 of 49 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


saw from the protein IN example, low helical content tells us little about protein compaction, which is critically important for H4 tail function. We therefore focus on available estimates of Rg for H4 from various sources. Table 9: Average helical propensity of H4 histone tail. Method MD simulation With TIP3P 38 MD simulation With OPC Circular dichroism (CD) 97 Chemical Shift (δ2D method) 103

Helical Propensity, % 3 1.6 ∼0 6.5

The simulated distribution of Rg obtained with OPC water model is compared to TIP3Pbased results in Figure 6: consistent with what we have seen for protein IN, H4 conformations are much more compact in TIP3P compared to OPC. In fact, the TIP3P-based ensemble is almost as compact as globular proteins with the same number of amino-acids, Table 10. In contrast, the much more extended H4 tail seen in OPC is consistent with a theoretical Table 10: Average of Rg for H4 histone tail from different methods. The average Rg for globular protein is from Table 1. The average Rg value for the simulation in OPC was computed from the last 7 µs of the 10 µs simulation. The value for TIP3P was taken from reference 38 . For the “polymer theory” (IDPs) and “polymer theory (fully denatured)” estimates see ”Methods”. Protein Globular Protein (Same Length as H4) H4 Histone Tail H4 Histone Tail H4 Histone Tail Fully Denatured Protein (Same Length as H4)

Method Experiment (Table 1) MD simulation With TIP3P 38 MD simulation With OPC Polymer Theory (IDPs) Polymer Theory (Fully Denatured)

Rg ( ˚ A) 10.1 9.7 14.3 14.8 17.9

estimate based on an established polymer theory, green curve in Figure 6. The average Rg of H4 seen in OPC water is also consistent, and TIP3P-based Rg is inconsistent, with an estimate based on a combination of experimental and simulation results 83 that relate the degree of protein compaction to its net charge per residue, Table 11. While the ensemble of conformations representing H4 tail is clearly much less compact than a globular protein 28



Page 30 of 49

of the same number of amino-acids, the ensemble is not as extended as the corresponding completely denatured protein, (blue curve in Figure 6). The observation is consistent with the notion that IDPs are not pure random coil polymers. Table 11: Rg of H4 histone tail normalized to RgRC (Random coil Rg ). Method MD simulation With TIP3P 38 MD simulation With OPC Prediction Based On The Net Charge Per Residue 83

Rg /RgRC 0.66 0.96 0.9

Figure 6: Probability distributions of Rg of H4-histone tail. The probability distribution obtained in TIP3P (black curve) was taken from reference 38 . The OPC-based distribution (red curve) was generated from the last 7 µs of the 10 µs simulation reported here, see “Methods”. The green and the blue curves correspond, respectively, to the mean-field theory predictions for an IDP and a fully denatured protein of the same size as H4, see “Methods”. Representative snapshots form the MD simulations are colored according to the electrostatic potential on the peptide surface 104 . Consistent with the above, free energy landscapes of H4 tail obtained in TIP3P and OPC are quite different, (Figure 7). We argue below that the OPC-based landscape, Rg vs number of contacts (NOC), is also more consistent with the biological function of this IDP. The presumed function of histone H4 tail is facilitating condensation of multiple nucleosomes 29


Page 31 of 49 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


into higher-order structures (arrays) by enabling attractive interactions between neighboring nucleosomes 105 . Presumably, H4 tails of one nucleosome attaches to the acidic patch area on the outside of the nearby nucleosomes 106 ; this action implies that H4 can extend and retract easily. In this respect, the OPC-based landscape is considerably more attractive than the TIP3P-based one. The shows a transition from a relatively compact to a rather extend structure, indicated by white arrows in Figure 7, costs only about 1 kB T in free energy in OPC, while the same transition can be prohibitively expensive, ∼ 5 kB T in TIP3P.

Figure 7: Free energy landscape of the H4-histone tail, extracted from the last 7 µs allatom simulations at 300 K in OPC (left) and TIP3P (right) water. White arrows exemplify possible transitions from compact to extended conformations, e.g. from Rg ≈ 10.5 ˚ A and ˚ ˚ NOC = 60 to Rg ≈16 A and NOC = 60 or Rg ≈ 18.5 A and NOC = 35. NOC = number of contacts between protein residues.

3.7

Robustness of Computed Water Properties to Hydrogen Mass Repartitioning (HMR)

We assessed the impact of HMR, as well as of doubling of the integration time step, on simulation accuracy by examining how these changes from the standard protocol affect four properties of water known to be sensitive to water model parameters. Several simulations were conducted, Table 12 shows the results. The density and heat of vaporization showed no 30



Page 32 of 49

change compared to the reference simulation, but the dielectric constant and diffusion showed small changes. These changes are well within the range of differences between computed properties of popular widely used water models. We therefore consider the impact of HMR and of switching to 4 fs time-step for water only as acceptable. We attribute the robustness to doubling of the integration time step to robustness of the analytical algorithm used by AMBER to enforce rigid geometry of the water molecule 107 . Table 12: Water properties of OPC water model with and without the use of Hydrogen Mass Repartitioning (HMR). The impact of doubling the integration time step dt from the default value of 2 fs is also examined. Simulation Condition No HMR, dt=2 fs No HMR, dt=4 fs HMR, dt=4 fs Experiment

4

Density(g/cm3 ) 0.997 0.996 0.996 0.997

Self-Diffusion(105 cm2 /s) 2.3 2.0 2.0 2.3

Dielectric Constant 78.4 77.5 77.5 78.4

Heat of Vaporization(kcal/mol) 10.63 10.63 10.63 10.52

CONCLUSION

Unconstrained atomistic simulations of intrinsically disordered proteins (IDPs) remain a challenge despite recent advances in force-fields; a significant part of the problem is that widely used “general purpose” water models predict structures of IDPs that are too compact compared to experiment. Here we have investigated whether a newly developed “general purpose” water model OPC can help – previously, the model showed significant improvement of the simulated ensembles of RNA structures, and its liquid state properties are in close agreement with experiment. Our main conclusions are of two types: methodological conclusions regarding potential benefits of OPC in the simulations of IDPs, and conclusions with respect to ensembles and free energy landscapes of three specific IDPs examined here: a common benchmark “protein IN” (55-residue apo N-terminal zinc-binding domain of HIV-1 integrase ), H4-histone tail (26 residues), and a 42-residue Aβ42 peptide. The main methodological conclusion is that OPC, in combination with a commonly 31


Page 33 of 49 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


used AMBER ff99SB force-field, produces much more realistic ensembles of the IDPs tested here, compared to TIP3P. In TIP3P, “protein IN” and H4-histone tail stay highly compact throughout the simulation, their radii of gyration (Rg ) remain close to the values expected of globular proteins of the corresponding size. In contrast, significant and multiple excursions away from the compact state are observed in OPC. For protein IN, for which a FRET-based estimate of Rg is available, the average value over the trajectory in OPC is significantly closer to the FRET estimate. Possible reasons for the remaining deviation are discussed. The compactness of H4-histone is close to an independent prediction based on an accepted polymer model. The computed Rg of Aβ42 peptide is in agreement with experiment. The computed Rg is nearly independent of the box size in this range, within the statistical error margin of our simulations. A brief analysis of the protein thermodynamics is consistent with this conclusion, but the story may be more nuanced when one looks at secondary structure populations, see below. Importantly, OPC passes the “negative control” – folding of a small globular protein (CLN025). In fact, the predicted folding free energy of CLN025 in OPC is substantially closer to experiment than in TIP3P. However, even in OPC, a small bias towards compact state remains. The predicted helical content of both protein IN and H4-histone tail in OPC and TIP3P is low, consistent with experiment. Based on the analysis of protein IN trajectories we have argued that, when low, protein helical content does not correlate well with the over-all compactness of the protein. The predicted β-sheet populations of Aβ42 are in qualitative agreement with experiment, but may depend on the solvent box size, within the range used in this work; longer trajectories are likely needed to make a quantitative assessment. In the meantime, based on our own analysis, and published results of others, we argue for a cautious approach to choosing solvent box size for simulations of IDPs. We suggest that the solvent box size should be at least four times the radius of gyration of the random coil corresponding to the IDP. All of the long MD simulations in explicit solvent reported here were performed with the

32



help of Hydrogen Mass Repartitioning (HMR), which permits doubling of the integration time-step. We have checked that key water properties either do not change, or change only slightly with the use of HMR compared to no HMR and the default time-step, further supporting the use of the method in atomistic simulations. Folding landscapes of IDPs are expected to be different from the folding funnel of regular proteins. We have investigated the free energy landscape of protein IN in some detail: in the (radius-of-gyration) – (number of residue contacts) space, the landscape resembles a shallow “tub”, in that conformations with substantially different degrees of compactness are within 2 kB T of each other. Conformations with very different secondary structure content co-exist within 1 kB T of the global free energy minimum. Given that ff99SB+OPC combination still has some residual bias towards compact states, the actual “folding tub” of protein IN may be even more shallow. However, similarly to globular proteins, conformations further away from the minimum tend to be less structured. The “shallow tub” feature of free energy landscapes of IDPs may be general, and relevant to their functions. For example, the computed free energy landscape of H4 histone tail in OPC is relatively flat over a significant range of compactness – a 60 % change in the radius of gyration of the tail can occur within 1kB T . We argue that this high predicted flexibility is important for the biological function of H4 histone tail as a facilitator of inter-nucleosome interactions. In summary, the general purpose water model OPC shows significant promise in atomistic simulations of intrinsically disordered proteins. The fact that this water model was not designed specifically for IDPs is suggestive of its transferability and versatility: the model may perform well in simulations of IDPs along with other biomolecules. While the number of test structures used here is necessarily small due to the high computational cost of the required simulations, the over-all conclusion we have reached with respect to the promise of OPC water model for atomistic simulations of IDPs is unlikely to be a statistical oddity. This is because physical reasons behind the improved performance of OPC are similar to those

33


Page 34 of 49

Page 35 of 49 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


behind TIP4P-D water model specifically designed for IDPs, and tested on a larger number of IDPs by now. Still, for a definitive conclusion, further improvements, additional testing of OPC, in combination with other gas-phase force-fields, will likely be needed. Judging by recent published examples of similar methodological works, such comprehensive investigations are now within reach, although will likely require significant, and possibly unique, computational resources.

ACKNOWLEDGMENTS Parviz Seifpanahi thanks Dr. Ramu Anandakrishnan for valuable discussions regarding the folding energy landscape and other suggestions, which greatly improved the manuscript. Parviz Seifpanahi thanks Dr. Igor Tolokh for assistance with several technical issues, and for valuable comments and suggestions, which greatly improved the manuscript. Parviz Seifpanahi thanks Dr. Daniel R. Roe for helping him resolve several issues associated with the use of Amber molecular dynamics package. The authors acknowledge Advanced Research Computing at Virginia Tech for providing computational resources and technical support that have contributed to the results reported within this paper. URL: http://www.arc.vt.edu

Supporting Information Additional details, tables and figures, referenced in this article are contained in the file supplement.pdf. Representative snapshots from the MD simulation of Aβ42 peptide are contained in the file snapshots.zip. This information is available free of charge via the Internet at http://pubs.acs.org.

34



Figure 8: OPC, the general purpose water model, can improve MD simulations of intrinsically disordered proteins. Shown is a free energy landscape of an intrinsically disordered protein computed using OPC water model. Protein conformations with substantially different degrees of compactness are within 2 kT of each other. Conformations with very different secondary structure content co-exist within 1 kT of the global free energy minimum.

TOC References (1) Ward, J.; Sodhi, J.; McGuffin, L.; Buxton, B.; Jones, D. Prediction and Functional Analysis of Native Disorder in Proteins from the Three Kingdoms of Life. Journal of Molecular Biology 2004, 337, 635 – 645. (2) Dunker, A. K.; Silman, I.; Uversky, V. N.; Sussman, J. L. Function and structure of inherently disordered proteins. Current Opinion in Structural Biology 2008, 18, 756 – 764. (3) Dunker, A. K.; Cortese, M. S.; Romero, P.; Iakoucheva, L. M.; Uversky, V. N. Flexible nets. FEBS Journal 2005, 272, 5129–5148. (4) Uversky, V. N.; Oldfield, C. J.; Dunker, A. K. Showing your ID: intrinsic disorder as

35


Page 36 of 49

Page 37 of 49 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


an ID for recognition, regulation and cell signaling. Journal of Molecular Recognition 2005, 18, 343–384. (5) Oldfield, C. J.; Cheng, Y.; Cortese, M. S.; Brown, C. J.; Uversky, V. N.; Dunker, A. K. Comparing and Combining Predictors of Mostly Disordered Proteins. Biochemistry 2005, 44, 1989–2000. (6) Uversky, V. N.; Oldfield, C. J.; Dunker, A. K. Intrinsically Disordered Proteins in Human Diseases: Introducing the D2 Concept. Annual Review of Biophysics 2008, 37, 215–246. (7) Brucale, M.; Schuler, B.; Samor, B. Single-Molecule Studies of Intrinsically Disordered Proteins. Chemical Reviews 2014, 114, 3281–3317. (8) Jensen, M. R.; Zweckstetter, M.; Huang, J.-r.; Blackledge, M. Exploring Free-Energy Landscapes of Intrinsically Disordered Proteins at Atomic Resolution Using NMR Spectroscopy. Chemical Reviews 2014, 114, 6632–6660. (9) Granata, D.; Baftizadeh, F.; Habchi, J.; Galvagnion, C.; De Simone, A.; Camilloni, C.; Laio, A.; Vendruscolo, M. The inverted free energy landscape of an intrinsically disordered peptide by simulations and experiments. Scientific Reports 2015, 5, 15449 EP –. (10) Best, R. B. Computational and theoretical advances in studies of intrinsically disordered proteins. Current Opinion in Structural Biology 2017, 42, 147 – 154. (11) Dunker, A.; Lawson, J.; Brown, C. J.; Williams, R. M.; Romero, P.; Oh, J. S.; Oldfield, C. J.; Campen, A. M.; Ratliff, C. M.; Hipps, K. W.; Ausio, J.; Nissen, M. S.; Reeves, R.; Kang, C.; Kissinger, C. R.; Bailey, R. W.; Griswold, M. D.; Chiu, W.; Garner, E. C.; Obradovic, Z. Intrinsically disordered protein. Journal of Molecular Graphics and Modelling 2001, 19, 26 – 59.

36



(12) Burger, V.; Gurry, T.; Stultz, C. Intrinsically Disordered Proteins: Where Computation Meets Experiment. Polymers 2014, 6, 2684–2719. (13) Duong, V. T.; Chen, Z.; Thapa, M. T.; Luo, R. Computational Studies of Intrinsically Disordered Proteins. The Journal of Physical Chemistry B 2018, 122, 10455–10469. (14) Henriques, J. a.; Cragnell, C.; Skepö, M. Molecular Dynamics Simulations of Intrinsically Disordered Proteins: Force Field Evaluation and Comparison with Experiment. J. Chem. Theory Comput. 2015, 11, 3420–3431. (15) Piana, S.; Donchev, A. G.; Robustelli, P.; Shaw, D. E. Water Dispersion Interactions Strongly Influence Simulated Structural Properties of Disordered Protein States. The Journal of Physical Chemistry B 2015, 119, 5113–5123. (16) Song, D.; Wang, W.; Ye, W.; Ji, D.; Luo, R.; Chen, H.-F. ff14IDPs force field improving the conformation sampling of intrinsically disordered proteins. Chemical Biology & Drug Design 2017, 89, 5–15. (17) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubmuller, H.; MacKerell Jr, A. D. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat Meth 2017, 14, 71–73. (18) Ye, W.; Ji, D.; Wang, W.; Luo, R.; Chen, H.-F. Test and Evaluation of ff99IDPs Force Field for Intrinsically Disordered Proteins. Journal of Chemical Information and Modeling 2015, 55, 1021–1029. (19) Best, R. B.; Zheng, W.; Mittal, J. Balanced Protein-Water Interactions Improve Properties of Disordered Proteins and Non-Specific Protein Association. J Chem Theory Comput 2014, 10, 5113–5124. (20) Huang, J.; MacKerell, A. D. Force field development and simulations of intrinsically disordered proteins. Current Opinion in Structural Biology 2018, 48, 40 – 48. 37


Page 38 of 49

Page 39 of 49 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


(21) Stanley, N.; Esteban-Martn, S.; Fabritiis, G. D. Progress in studying intrinsically disordered proteins with atomistic simulations. Progress in Biophysics and Molecular Biology 2015, 119, 47 – 52. (22) Bergonzo, C.; Cheatham III, T. Improved Force Field Parameters Lead to a Better Description of RNA Structure. J. Chem. Theory. Comput. 2015, 11, 3969–3972. (23) Yang, C.; Lim, M.; Kim, E.; Pak, Y. Predicting RNA Structures via a Simple van der Waals Correction to an All-Atom Force Field. J. Chem. Theory Comput. 2017, 13, 395–399. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics 1983, 79, 926–935. (25) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665–9678. (26) Grosberg, A. I.; Khokhlov, A. R. Statistical physics of macromolecules; New York : American Institute of Physics, 1994. (27) Anandakrishnan, R.; Izadi, S.; Onufriev, A. V. Why Computed Protein Folding Landscapes Are Sensitive to the Water Model. Journal of Chemical Theory and Computation 2018, 15, 625–636. (28) Robustelli, P.; Piana, S.; Shaw, D. E. Developing a molecular dynamics force field for both folded and disordered protein states. Proceedings of the National Academy of Sciences 2018, 115, E4758–E4766. (29) Wu, H.-N.; Jiang, F.; Wu, Y.-D. Significantly Improved Protein Folding Thermody-

38



namics Using a Dispersion-Corrected Water Model and a New Residue-Specific Force Field. The Journal of Physical Chemistry Letters 2017, 8, 3199–3205. (30) Onufriev, A. V.; Izadi, S. Water models for biomolecular simulations. Wiley Interdisciplinary Reviews: Computational Molecular Science 2018, 8, e1347. (31) Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014, 5, 3863–3871. (32) Vega, C.; Abascal, J. L. F. Simulating water with rigid non-polarizable models: a general perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663–19688. (33) Vega, C.; Abascal, J. L. F.; Conde, M. M.; Aragones, J. L. What ice can teach us about water interactions: a critical comparison of the performance of different water models. Faraday Discuss. 2009, 141, 251–276. (34) Gao, K.; Yin, J.; Henriksen, N. M.; Fenley, A. T.; Gilson, M. K. Binding Enthalpy Calculations for a Neutral HostGuest Pair Yield Widely Divergent Salt Effects across Water Models. J. Chem. Theory. Comput. 2015, 11, 4555–4564. (35) Javanainen, M.; Lamberg, A.; Cwiklik, L.; Vattulainen, I.; Ollila, O. H. S. Atomistic Model for Nearly Quantitative Simulations of Langmuir Monolayers. Langmuir : the ACS journal of surfaces and colloids 2018, 34, 2565–2572. (36) Gabrieli, A.; Sant, M.; Izadi, S.; Shabane, P. S.; Onufriev, A. V.; Suffritti, G. B. Hightemperature dynamic behavior in bulk liquid water: A molecular dynamics simulation study using the OPC and TIP4P-Ew potentials. Frontiers of Physics 2017, 13, 138203. (37) Potoyan, D. A.; Papoian, G. A. Energy Landscape Analyses of Disordered Histone Tails Reveal Special Organization of Their Conformational Dynamics. Journal of the American Chemical Society 2011, 133, 7405–7415.

39


Page 40 of 49

Page 41 of 49 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


(38) Winogradoff, D.; Echeverria, I.; Potoyan, D. A.; Papoian, G. A. The Acetylation Landscape of the H4 Histone Tail: Disentangling the Interplay between the Specific and Cumulative Effects. Journal of the American Chemical Society 2015, 137, 6245– 6253. (39) Liu, H.; Duan, Y. Effects of Posttranslational Modifications on the Structure and Dynamics of Histone H3 N-Terminal Peptide. Biophysical Journal 2008, 94, 4579– 4585. (40) Erler, J.; Zhang, R.; Petridis, L.; Cheng, X.; Smith, J.; Langowski, J. The Role of Histone Tails in the Nucleosome: A Computational Study. Biophysical Journal 2014, 107, 2911 – 2922. (41) Yang, D.; Arya, G. Structure and binding of the H4 histone tail and the effects of lysine 16 acetylation. Physical chemistry chemical physics : PCCP 2011, 13, 2911–2921. (42) Zhang, R.; Erler, J.; Langowski, J. Histone Acetylation Regulates Chromatin Accessibility: Role of H4K16 in Inter-nucleosome Interaction. Biophysical Journal 2017, 112, 450–459. (43) Collepardo-Guevara, R.; Portella, G.; Vendruscolo, M.; Frenkel, D.; Schlick, T.; Orozco, M. Chromatin Unfolding by Epigenetic Modifications Explained by Dramatic Impairment of Internucleosome Interactions: A Multiscale Computational Study. J. Am. Chem. Soc. 2015, 137, 10205–10215. (44) Mller-Spth, S.; Soranno, A.; Hirschfeld, V.; Hofmann, H.; Regger, S.; Reymond, L.; Nettels, D.; Schuler, B. Charge interactions can dominate the dimensions of intrinsically disordered proteins. Proceedings of the National Academy of Sciences 2010, 107, 14609–14614. (45) Wuttke, R.; Hofmann, H.; Nettels, D.; Borgia, M. B.; Mittal, J.; Best, R. B.;

40



Schuler, B. Temperature-dependent solvation modulates the dimensions of disordered proteins. Proceedings of the National Academy of Sciences 2014, 111, 5213–5218. (46) Cai, M.; Zheng, R.; Caffrey, M.; Craigie, R.; Clore, G. M.; Gronenborn, A. M. Solution structure of the N-terminal zinc binding domain of HIV-1 integrase. Nat Struct Mol Biol 1997, 4, 567–577. (47) Coskuner, O.; Uversky, V. N. Tyrosine Regulates -Sheet Structure Formation in Amyloid-42: A New Clustering Algorithm for Disordered Proteins. Journal of Chemical Information and Modeling 2017, 57, 1342–1358. (48) L¨ uhrs, T.; Ritter, C.; Adrian, M.; Riek-Loher, D.; Bohrmann, B.; Döbeli, H.; Schubert, D.; Riek, R. 3D structure of Alzheimer’s amyloid-(1–42) fibrils. Proceedings of the National Academy of Sciences 2005, 102, 17342–17347. (49) Morimoto, A.; Irie, K.; Murakami, K.; Masuda, Y.; Ohigashi, H.; Nagao, M.; Fukuda, H.; Shimizu, T.; Shirasawa, T. Analysis of the Secondary Structure of Amyloid (A42) Fibrils by Systematic Proline Replacement. Journal of Biological Chemistry 2004, 279, 52781–52788. (50) Sgourakis, N. G.; Yan, Y.; McCallum, S. A.; Wang, C.; Garcia, A. E. The Alzheimers Peptides A40 and 42 Adopt Distinct Conformations in Water: A Combined MD / NMR Study. Journal of Molecular Biology 2007, 368, 1448 – 1457. (51) Weber, O. C.; Uversky, V. N. How accurate are your simulations? Effects of confined aqueous volume and AMBER FF99SB and CHARMM22/CMAP force field parameters on structural ensembles of intrinsically disordered proteins: Amyloid-b(42) in water. Intrinsically Disordered Proteins 2017, 5, e1377813–e1377813. (52) Baumketner, A.; Bernstein, S. L.; Wyttenbach, T.; Bitan, G.; Teplow, D. B.; Bowers, M. T.; Shea, J.-E. Amyloid -protein monomer structure: A computational and experimental study. Protein Science 2006, 15, 420–428. 41


Page 42 of 49

Page 43 of 49 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


(53) Serpell, L. C. Alzheimers amyloid fibrils: structure and assembly. Biochimica et Biophysica Acta (BBA) - Molecular Basis of Disease 2000, 1502, 16 – 30. (54) Roche, J.; Shen, Y.; Lee, J. H.; Ying, J.; Bax, A. Monomeric A140 and A142 Peptides in Solution Adopt Very Similar Ramachandran Map Distributions That Closely Resemble Random Coil. Biochemistry 2016, 55, 762–775. (55) Brown, A. M.; Bevan, D. R. Influence of sequence and lipid type on membrane perturbation by human and rat amyloid -peptide (142). Archives of Biochemistry and Biophysics 2017, 614, 1 – 13. (56) Honda, S.; Akiba, T.; Kato, Y. S.; Sawada, Y.; Sekijima, M.; Ishimura, M.; Ooishi, A.; Watanabe, H.; Odahara, T.; Harata, K. Crystal Structure of a Ten-Amino Acid Protein. J. Am. Chem. Soc. 2008, 130, 15327–15331. (57) Anandakrishnan, R.; Drozdetski, A.; Walker, R. C.; Onufriev, A. V. Speed of Conformational Change: Comparing Explicit and Implicit Solvent Molecular Dynamics Simulations. Biophys. J. 2015, 108, 1153–1164. (58) Hatfield, M. P.; Murphy, R. F.; Lovas, S. Molecular dynamics analysis of the conformations of a beta-hairpin miniprotein. The journal of physical chemistry. B 2010, 114, 3028–3037. (59) Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How Fast-Folding Proteins Fold. Science 2011, 334, 517–520. (60) Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. H++ 3.0: automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations. Nucleic Acids Research 2012, 40, W537–W541. (61) Case, D.; Berryman, J.; Betz, R.; Cerutti, D.; Cheatham III, T.; Darden, T.; Duke, R.; Giese, T.; Gohlke, H.; Goetz, A.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Ko42



valenko, T., A. amd Lee; LeGrand, S.; Li, P.; Luchko, T.; Luo, R.; Madej, B.; Merz, K.; Monard, G.; Needham, P.; Nguyen, H.; Nguyen, H.; Omelyan, I.; Onufriev, A.; Roe, D.; Roitberg, A.; Salomon-Ferrer, R.; Simmerling, C.; Smith, W.; Swails, J.; Walker, R.; Wang, J.; Wolf, R.; Wu, X.; York, D.; Kollman, P. AMBER 2015 ; University of California, San Francisco, 2015. (62) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems. The Journal of Chemical Physics 1993, 98, 10089– 10092. (63) Izadi, S.; Onufriev, A. V. Accuracy limit of rigid 3-point water models. The Journal of chemical physics 2016, 145 . (64) Pastor, R. W.; Brooks, B. R.; Szabo, A. An analysis of the accuracy of Langevin and molecular dynamics algorithms. Molecular Physics 1988, 65, 1409–1419. (65) paul Ryckaert, J.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys 1977, 327–341. (66) Hopkins, C. W.; Le Grand, S.; Walker, R. C.; Roitberg, A. E. Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning. Journal of Chemical Theory and Computation 2015, 11, 1864–1874. (67) Roe, D. R.; Cheatham, T. E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. Journal of Chemical Theory and Computation 2013, 9, 3084–3095. (68) K¨ uhrová, P.; De Simone, A.; Otyepka, M.; Best, R. B. Force-Field Dependence of Chignolin Folding and Misfolding: Comparison with Experiment and Redesign. Biophys. J. 2012, 102, 1897–1906.

43


Page 44 of 49

Page 45 of 49 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


(69) Honda, S.; Akiba, T.; Kato, Y. S.; Sawada, Y.; Sekijima, M.; Ishimura, M.; Ooishi, A.; Watanabe, H.; Odahara, T.; Harata, K. Crystal Structure of a Ten-Amino Acid Protein. J. Am. Chem. Soc. 2008, 130, 15327–15331. (70) Flory, P. J. Principles of Polymer Chemistry, 1st ed.; 1; Cornell University Press: Ithaca, NY 14850, 1953; Vol. 1. (71) Sanchez, I. C.; Lacombe, R. H. Statistical Thermodynamics of Polymer Solutions. Macromolecules 1978, 11, 1145–1156. (72) Ziv, G.; Haran, G. Protein Folding, Protein Collapse, and Tanfords Transfer Model: Lessons from Single-Molecule FRET. Journal of the American Chemical Society 2009, 131, 2942–2947. (73) M¨ uller-Späth, S.; Soranno, A.; Hirschfeld, V.; Hofmann, H.; R¨ uegger, S.; Reymond, L.; Nettels, D.; Schuler, B. Charge interactions can dominate the dimensions of intrinsically disordered proteins. Proceedings of the National Academy of Sciences 2010, 107, 14609–14614. (74) Zerze, G. H.; Best, R. B.; Mittal, J. Sequence- and Temperature-Dependent Properties of Unfolded and Disordered Proteins from Atomistic Simulations. The Journal of Physical Chemistry B 2015, 119, 14622–14630. (75) Ganguly, D.; Chen, J. Atomistic Details of the Disordered States of KID and pKID. Implications in Coupled Binding and Folding. Journal of the American Chemical Society 2009, 131, 5214–5223. (76) Tran, H. T.; Mao, A.; Pappu, R. V. Role of BackboneSolvent Interactions in Determining Conformational Equilibria of Intrinsically Disordered Proteins. Journal of the American Chemical Society 2008, 130, 7380–7392.

44



(77) Vitalis, A.; Wang, X.; Pappu, R. V. Atomistic Simulations of the Effects of Polyglutamine Chain Length and Solvent Quality on Conformational Equilibria and Spontaneous Homodimerization. Journal of Molecular Biology 2008, 384, 279 – 297. (78) Vitalis, A.; Wang, X.; Pappu, R. V. Quantitative Characterization of Intrinsic Disorder in Polyglutamine: Insights from Analysis Based on Polymer Theories. Biophysical Journal 2007, 93, 1923 – 1937. (79) Tran, H. T.; Pappu, R. V. Toward an Accurate Theoretical Framework for Describing Ensembles for Proteins under Strongly Denaturing Conditions. Biophysical Journal 2006, 91, 1868 – 1886. (80) Soranno, A.; Longhi, R.; Bellini, T.; Buscaglia, M. Kinetics of contact formation and end-to-end distance distributions of swollen disordered peptides. Biophys. J. 2009, 96, 1515–1528. (81) Mukhopadhyay, S.; Krishnan, R.; Lemke, E. A.; Lindquist, S.; Deniz, A. A. A natively unfolded yeast prion monomer adopts an ensemble of collapsed and rapidly fluctuating structures. Proceedings of the National Academy of Sciences 2007, 104, 2649–2654. (82) Crick, S. L.; Jayaraman, M.; Frieden, C.; Wetzel, R.; Pappu, R. V. Fluorescence correlation spectroscopy shows that monomeric polyglutamine molecules form collapsed structures in aqueous solutions. Proceedings of the National Academy of Sciences 2006, 103, 16764–16769. (83) Mao, A. H.; Crick, S. L.; Vitalis, A.; Chicoine, C. L.; Pappu, R. V. Net charge per residue modulates conformational ensembles of intrinsically disordered proteins. Proceedings of the National Academy of Sciences 2010, 107, 8183–8188. (84) Ding, F.; Jha, R. K.; Dokholyan, N. V. Scaling Behavior and Structure of Denatured Proteins. Structure 2005, 13, 1047 – 1054.

45


Page 46 of 49

Page 47 of 49 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


(85) Mittal, J.; Best, R. B. Thermodynamics and kinetics of protein folding under confinement. Proceedings of the National Academy of Sciences 2008, 105, 20233–20238. (86) Lucent, D.; Vishal, V.; Pande, V. S. Protein folding under confinement: A role for solvent. Proceedings of the National Academy of Sciences 2007, 104, 10430–10434. (87) Oldfield, C. J.; Meng, J.; Yang, J. Y.; Yang, M. Q.; Uversky, V. N.; Dunker, A. K. Flexible nets: disorder and induced fit in the associations of p53 and 14-3-3 with their partners. BMC Genomics 2008, 9, S1. (88) Uversky, V. N.; Dav, V.; Iakoucheva, L. M.; Malaney, P.; Metallo, S. J.; Pathak, R. R.; Joerger, A. C. Pathological Unfoldomics of Uncontrolled Chaos: Intrinsically Disordered Proteins and Human Diseases. Chemical Reviews 2014, 114, 6844–6879. (89) Honda, S.; Yamasaki, K.; Sawada, Y.; Morii, H. 10 residue folded peptide designed by segment statistics. Structure 2004, 12, 1507–1518. (90) K¨ uhrová, P.; De Simone, A.; Otyepka, M.; Best, R. B. Force-field dependence of chignolin folding and misfolding: comparison with experiment and redesign. Biophys. J. 2012, 102, 1897–1906. (91) Nguyen, H.; Maier, J.; Huang, H.; Perrone, V.; Simmerling, C. Folding simulations for proteins with diverse topologies are accessible in days with a single physics-based force field and implicit solvent. J. Am. Chem. Soc. 2014, (92) Pang, Y.-P. FF12MC: A revised AMBER forcefield and new protein simulation protocol. Proteins: Structure, Function, and Bioinformatics 2016, 84, 1490–1516. (93) McKiernan, K. A.; Husic, B. E.; Pande, V. S. Modeling the mechanism of CLN025 beta-hairpin formation. The Journal of chemical physics 2017, 147 . (94) Wolynes, P. G. Folding funnels and energy landscapes of larger proteins within the capillarity approximation. Proc Natl Acad Sci U S A 1997, 94, 6170–6175. 46



(95) Ovchinnikov, V.; Cecchini, M.; Karplus, M. A Simplified Confinement Method for Calculating Absolute Free Energies and Free Energy and Entropy Differences. The Journal of Physical Chemistry B 2013, 117, 750–762. (96) Onuchic, J. N.; Socci, N. D.; Luthey-Schulten, Z.; Wolynes, P. G. Protein folding funnels: the nature of the transition state ensemble. Folding and Design 1996, 1, 441 – 450. (97) Wang, X.; Moore, S. C.; Laszckzak, M.; Ausi, J. Acetylation Increases the -Helical Content of the Histone Tails of the Nucleosome. Journal of Biological Chemistry 2000, 275, 35013–35020. (98) Allahverdi, A.; Yang, R.; Korolev, N.; Fan, Y.; Davey, C. A.; Liu, C.-F.; Nordenskild, L. The effects of histone H4 tail acetylations on cation-induced chromatin folding and self-association. Nucleic Acids Research 2011, 39, 1680–1691. (99) Gao, M.; Nadaud, P. S.; Bernier, M. W.; North, J. A.; Hammel, P. C.; Poirier, M. G.; Jaroniec, C. P. Histone H3 and H4 N-Terminal Tails in Nucleosome Arrays at Cellular Concentrations Probed by Magic Angle Spinning NMR Spectroscopy. Journal of the American Chemical Society 2013, 135, 15278–15281. (100) Luger, K.; Mader, A. W.; Richmond, R. K.; Sargent, D. F.; Richmond, T. J. Crystal structure of the nucleosome core particle at 2.8 ˚ A resolution. Nature 1997, 389, 251– 260. (101) Luger, K.; Richmond, T. J. The histone tails of the nucleosome. Current Opinion in Genetics & Development 1998, 8, 140 – 146. (102) Gao, M.; Nadaud, P. S.; Bernier, M. W.; North, J. A.; Hammel, P. C.; Poirier, M. G.; Jaroniec, C. P. Histone H3 and H4 N-Terminal Tails in Nucleosome Arrays at Cellular Concentrations Probed by Magic Angle Spinning NMR Spectroscopy. Journal of the American Chemical Society 2013, 135, 15278–15281. 47


Page 48 of 49

Page 49 of 49 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


(103) Camilloni, C.; De Simone, A.; Vranken, W. F.; Vendruscolo, M. Determination of Secondary Structure Populations in Disordered States of Proteins Using Nuclear Magnetic Resonance Chemical Shifts. Biochemistry 2012, 51, 2224–2231. (104) Gordon, J. C.; Fenley, A. T.; Onufriev, A. An analytical approach to computing biomolecular electrostatic potential. II. Validation and applications. The Journal of Chemical Physics 2008, 129, 075102. (105) Saurabh, S.; Glaser, M. A.; Lansac, Y.; Maiti, P. K. Atomistic Simulation of Stacked Nucleosome Core Particles: Tail Bridging, the H4 Tail, and Effect of Hydrophobic Forces. The Journal of Physical Chemistry B 2016, 120, 3048–3060. (106) Kalashnikova, A. A.; Porter-Goff, M. E.; Muthurajan, U. M.; Luger, K.; Hansen, J. C. The role of the nucleosome acidic patch in modulating higher order chromatin structure. Journal of The Royal Society Interface 2013, 10 . (107) Miyamoto, S.; Kollman, P. A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952–962.

48


General Purpose Water Model Can Improve Atomistic Simulations of

Recommend Documents