How Well Can Implicit Solvent Simulations Explore Folding Pathways

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX-XXX

pubs.acs.org/JCTC

How Well Can Implicit Solvent Simulations Explore Folding Pathways? A Quantitative Analysis of α‑Helix Bundle Proteins Qiang Shao*,†,‡ and Weiliang Zhu†,‡ †

Drug Discovery and Design Center, CAS Key Laboratory of Receptor Research, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, 555 Zuchongzhi Road, Shanghai 201203, China ‡ University of Chinese Academy of Sciences, Beijing 100049, China S Supporting Information *

ABSTRACT: Protein folding has been posing challenges for molecular simulation for decades. Implicit solvent models are sought as routes to increase the capability of simulation, with trade-offs between computational speed and accuracy. Here, we systematically investigate the folding of a variety of α-helix bundle proteins ranging in size from 46 to 102 amino acids using a state-of-the-art force field and an implicit solvent model. The accurate all-atom simulated folding is enabled for six proteins, including for the first time a successful folding of protein with >100 amino acids in implicit solvent. The detailed free-energy landscape analysis sheds light on a set of general principles underlying the folding of α-helix bundle proteins, suggesting a hybrid framework/nucleation-condensation mechanism favorably adopted in implicit solvent condition. The similarities and discrepancies of the folding pathways measured among the present implicit solvent simulations and previously reported experiments and explicit solvent simulations are deeply analyzed, providing quantitative assessment for the availability and limitation of implicit solvent simulation in exploring the folding transition of large-size proteins.

1. INTRODUCTION

Implicit solvent models have been of particular interest which replace explicitly simulated solvent molecules with effective intraprotein interactions, with the representative being the generalized Born (GB) models for the implicit treatment of water solvent8−10 as well as bilipid membrane.11−13 Accordingly, while the full protein resolution remains, the solvent degrees of freedom are removed, and the simulation demands are dramatically reduced. The use of implicit solvent models effectively speeds up the simulation process14,15 and thus can technically expand folding simulation to proteins with diverse sizes and structural topologies.16−19 It has been reported in multiple anecdotal cases that the folding of the conformational states in implicit solvent matches the counterparts measured in explicit solvent simulations or experiments.20−24 On the other hand, an extensive collection of simulation studies revealed that the implicit GB models exacerbate the secondary structure bias of force fields23,25 and produce folding thermodynamic quantities deviated from those measured by explicit solvent simulations.26−29 It is noteworthy that the test proteins in this literature are mainly short peptides (≤35 amino acids (AA)).20−29 Recently, Nguyen et al. reported that the combined use of the AMBER all-atom FF14SBonlysc force field and the GB-Neck2 implicit model (denoted igb8 according to the AMBER input options) can dramatically fold proteins in

Protein folding, a molecular self-assembly process in which a disordered polypeptide chain collapses to form a well-defined three-dimensional native structure, is one of the most perplexing problems in molecular biology.1 As highlighted in the famed Levinthal paradox,2,3 instead of the random search in protein conformational space that is practically impossible considering the vast amount of conformational degrees of freedoms, proteins are supposed to search and converge quickly to their native states through well-defined folding pathways, sometimes in the time-scale of microseconds. The underlying principles that govern the folding process along specific pathway(s) are at the heart of the protein folding problem.1,4,5 Computational simulation provides an access to explore protein folding at the atomistic level. The difficulty in sampling the naturally rugged configuration of free-energy landscape (induced by the complex interactions from amino acid residues with varied physicochemical properties), however, makes allatom molecular dynamics (MD) simulation extraordinarily time-consuming and computationally expensive. Despite the significant advances in specialized computer hardware that extend MD simulation time-scale to hundreds of microseconds and even milliseconds,6,7 the strategies aiming to economize simulation time and computational resources (especially on general-purpose computers) are still desirable, with trade-offs between resolution, accuracy, and computational speed. © XXXX American Chemical Society

Received: July 7, 2017 Published: November 9, 2017 A

DOI: 10.1021/acs.jctc.7b00726 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

previously reported results from various experiments and explicit solvent simulations to quantitatively assess the availability and limitation of implicit solvent models in modeling protein folding.

a large range of sizes (10−92 AA), leaving their folding mechanisms not deeply investigated.16 Therefore, how well the state-of-the-art protein force fields developed independently of implicit solvent models can together explore the detailed folding pathways of large-size proteins that attract much attention from biophysicists still remains unclear. The question also arises as to how far implicit solvent simulation can go in the theoretical study of protein folding. In this work, we systematically investigate the folding behavior of a variety of α-helix bundle proteins ranging in size from 46 to 102 AA to (1) assess the combined use of stateof-the-art force field and implicit solvent model in exploring the folding transitions of large-size proteins and (2) reveal a set of general principles underlying the folding of α-helical proteins in implicit solvent. Proteins that were simulated by all-atom MD simulations include the B domain of protein A (BdpA), engrailed homeodomain protein (EnHD) and its mutant (EnHD_K52A), α3W, α3D, and S824. The protein sequences are given in Table S1 in the Supporting Information (SI). We chose this particular set of test proteins based on three criteria. First, they all exhibit well-defined purely α-helical folds that are more easily accessible than similarly sized α/β proteins by implicit solvent simulation, providing enough test cases for an in silico study (only quite a few α/β proteins have been reported to be folded in implicit solvent conditions which are mainly small proteins and peptides17,18,30). Second, the presence of a great amount of preceding experimental and explicit solvent simulation studies on these test proteins (particularly BdpA, EnHD, EnHD_K52A, and α3D) allows for the comparison in various data (e.g., the melting temperature, the structural features of unfolded state, the folding models, the stability of individual helix segments, and their interplay along the folding pathway, etc.).6,31−51 Third, the detailed analysis of the similarities and discrepancies of the folding transitions of proteins having different sequences but similar topologies facilitates the understanding of the inherent sequence-dependent folding mechanism and the effects from solvation condition.45,46,52 The recently developed AMBER FF14SBonlysc force field and GB-Neck2 implicit model, arguably the best combination for folding diverse proteins in implicit solvent condition,16 were used to model protein atoms and solvation effects, respectively. Also notable is that an integrated-tempering-sampling (ITS) enhanced sampling algorithm 53,54 was utilized in MD simulations (hereinafter referred to as ITS-MD simulations) to make efficient sampling of the configuration space and accurate weighting of the equilibrium ensembles in endurable simulation time that provide an effective access to explore the folding transition pathway(s) on the free-energy surface. Such a method employs the temperatures integrated in single MD trajectory to adjust the sampling distribution over the desired potential energy range so as to enhance the conformational sampling in simulation,53,54 similar to the enveloping distribution sampling method.55,56 The advantage of the sampling technique of ITS-MD has been extensively displayed in various biomolecular simulations including protein folding17,57−61 and conformational changes.62−64 For the six proteins including the one with >100 AA (S824) that is folded for the first time in implicit solvent simulation, the detailed structural characterization of the important transient states (e.g., unfolded and intermediate) that govern the folding pathways is analyzed. The measured folding pathways and folding mechanisms are comprehensively compared to

2. COMPUTATIONAL METHODS 2.1. System Preparation for Implicit Solvent ITS-MD Simulations. The all-atom MD simulations implemented with ITS enhanced sampling algorithm (ITS-MD) were performed by AMBER14 software65 on CPUs. Protein atoms were modeled by the AMBER FF14SBonlysc all-atom force field.16 Meanwhile, the solvent effects were mimicked by the recently improved implicit GB model, GB-Neck2 (igb8):66 (1) in each simulation step, the electrostatic (polar) term of protein solvation free-energy was calculated by such a GB method, and (2) the nonpolar term of solvation free-energy was approximated by ΔGnonpol = γ × SASA, where γ is the surface tension coefficient and SASA is the total solvent accessible surface area of protein. For each protein, a fully expanded structure is first generated by the tleap program in AMBER14, with constant protonation states for ionizable residues corresponding to pH = 7 solvent condition. To obtain different initial structures with varied atom coordinates and velocities for ITS-MD simulations, the expanded structure was subjected to the preprocesses performed by conventional MD including the minimization, heating-up, and short-time equilibrium at 300 K temperature. In these preprocesses, the SHAKE algorithm67 with a relative geometric tolerance of 10−5 was used to constrain all chemical bonds. No other constraints are artificially applied on protein atoms. After 50000 steps of minimization, the temperature of the simulation system was increased to 300 K within 120 ps in a total of 6 stages as suggested by an online tutorial of AMBER for the case study of folding Trp-cage peptide.68 Then the system was equilibrated at 300 K for ∼5 ns using Langevin dynamics with the collision frequency of 3 ps−1. Multiple snapshots selected from the short-time equilibrium process were then used as the initial structures for individual independent long-time trajectories (production runs) of ITSMD simulations. 2.2. Simulation Details of ITS-MD. ITS-MD is an integrated tempering approach that utilizes the temperatures integrated in a single MD trajectory to (1) expand the potential energy sampled in simulation and (2) adjust the sampling distribution in the expanded potential energy range, so as to make a thorough sampling on the potential energy surface.53,54 In such a method, the modified potential energy (V′) under simulation can be written as a summation over the integration of a designed temperature range {βi} N

V′ = −

1 ln ∑ nie−βiV β0 i = 1

(i = 1, ..., N ) (1)

where V is the original potential energy, βi = th

1 kBTi

(kB is the

Boltzmann constant, and Ti is the i point in the temperature range), β0 is the inverse temperature of simulated system, N is the number of temperature points in the integration, and {ni} are the crucial parameters that determine whether a thorough sampling can be achieved in the entire potential energy range in ITS-MD simulation. The potential energy distribution function is derived from eq N 1: p(V ′) = e−βV ′ = ∑i = 1 nie−βiV . Ideally, each term in the B


Article


where Ni’s are the contacts of residue i, dij is the Cα−Cα distance (in unit of Å) between residues i and j at any frame extracted from trajectory, and d0ij is the same distance in the native state. The native contacts are equivalent to the sum of native Cα−Cα contacts between residues that are separated by ≥7 residues in sequence and whose distances are shorter than 10 Å for more than 80% of the time in the folded state. The value of Q is essentially in the range of 0 to 1, with a smaller value indicating a less native-like contact profile. The distance from the native structure is assessed by the “nativeness” of a residue i6

distribution function contributes equally to generate an even distribution in the overall potential energy range. To do so, the optimal values of {ni} parameters are determined in a preliminary iteration process as described in detail in the Supporting Information of ref 54. The long-time production runs are eventually performed using the optimal {ni} for each trajectory of ITS-MD simulation. The biased force in the simulation is N

F′ = −

∑i = 1 βi nie−βiV ∂V ′ = F N ∂r β0 ∑i = 1 nie−βiV

(2)

di = 1 − e−0.5 * MSDi

where F is the original molecular force in the simulation system. In each trajectory of ITS-MD, the step size was 0.002 ps, and the temperature was controlled around 300 K using Langevin dynamics with the collision frequency of 3 ps−1. The SHAKE algorithm67 with a relative geometric tolerance of 10−5 was used to constrain all chemical bonds. The cutoff of 30 Å was used to truncate nonbonded pairs (on an atom-by-atom basis). The surface tension coefficient was set as 0.005 kcal/mol/Å2 (default), and the salt concentration was 0.0 M (default) for implicit solvent mimicking. The calculation data were deposited every 2.0 ps. The data from all production runs were collected together for statistical analysis. 2.3. Trajectory Analyses. 2.3.1. Statistical Analysis Constructing Free-Energy Landscape. While the converged calculation of ITS-MD simulation yields a biased nonBoltzmann distribution, the unbiased distribution function needs to be recovered by multiplying a reweighting factor:

where MSDi is the mean-square deviation (in unit of Å2) from the native structure for a combination of 5 amino-acid residues in sequence centered on the ith residue. The ideally zero value of di implies that the ith residue stays at its native state, whereas the value is close to unity if the residue is deviated from the native structure. We calculated the di value for each snapshot in simulation and then reweighted the ensemble average of di for the unfolded (U) state by ⟨di⟩U = ⟨diW⟩U/⟨W⟩U. 2.3.3. Measurement of Secondary Structure Content, the Percentage of Solvent Accessible Surface Area, and Detailed Side-Chain Native Contacts. These parameters were used to display detailed structure changes in protein folding. The secondary structure of each residue could be evaluated by the ptraj program in AMBER14 following the DSSP method by Kabsch and Sander.70 The secondary structure content of each region of protein was then calculated by averaging over all residues involved. The percentage of solvent accessible surface area (PSASA), which can be used to evaluate the compactness of protein structure (with a large value indicating a compact structure), was calculated relative to the difference between the maximal and minimal values in t he simulation ( SASA − SASA max ). The SASA of proteins were measured using

N

W = e β(V ′ (r) − V (r)) = e−βV (r)/p(V ′) = e−βV (r)/∑ nie−βiV i=1

(3)

The true reweighted ensemble average of a quantity A is calculated by ⟨A⟩ = ⟨AW⟩/⟨W⟩, and the free-energy landscape along any n-dimensional reaction coordinate(s) X = (x1,...,xn) is achieved by ΔG(X ) = −

P(X ) 1 ln Pmax β

SASA min − SASA max

“Naccess” program.71 The minimal (SASAmin) and maximal (SASAmax) values for each protein were organized in Table S2. Of the detailed side-chain native contacts that were defined from experimental native structure of protein, the side-chain contacts among different helices could be used to estimate the construction of the global tertiary structure (the native contacts within individual helices were excluded). The side-chain contacts were considered to be formed as the minimal distance of heavy atoms between the pair of residue side chains is less than 5.0 Å. 2.3.4. Measurement of Melting Temperature (Tm). The folding/unfolding thermodynamic quantities at any desired temperature can be measured from simulation. For instance, the profile of specific heat capacity as a function of temperature can be calculated by

(4)

where P(X) is the unbiased probability distribution along X that is obtained through a histogram analysis of the ITS-MD snapshots taking into account their real visiting probabilities recovered by reweighting factor W; Pmax is the maximum of the distribution to be subtracted to ensure that ΔG(X) = 0 for the lowest free-energy minimum.69 To show the convergence of the ITS-MD simulations, the error bars for the folding free-energy profiles along specific reaction coordinates were calculated using the data from different trajectories: first, for each trajectory, the free-energy profile was calculated, respectively; then the error bars were determined by calculating the standard deviation of free-energy values from all trajectories. 2.3.2. Measurement of Fraction of Native Contacts (Q) and Distance from the Native Structure in the Unfolded State. As defined by Lindorff-Larsen et al.,6 Q is the fraction of native contacts formed at each frame that can be used as a reaction coordinate of protein folding N

Q=

N

∑i =res1 ∑ j =i 1

Cp = (⟨(V 2 × W )⟩/⟨W ⟩ − (⟨V × W ⟩/⟨W ⟩)2 )/kT 2 (7)

The melting temperature of protein can be then extrapolated from the pronounced peak position in the profile.

3. RESULTS 3.1. Repeated Folding and Unfolding Simulated in Implicit Solvent. For each protein in this work, we performed 4−8 simulations, each between 0.4 and 1.6 μs long. The trajectories were started from different initial structures with

1 0

1 + e10(dij − (dij + 1)) Nres ∑i = 1 Ni

(6)

(5) C


Article


sampled in a large range of RMSD and Rg. Accordingly, all proteins can be folded into native-like structures with low RMSDs (Figure 1). More importantly, the folded structures are visited repeatedly, collecting plenty of folding events (Table 1) that allow for exploring folding pathways in a quantitative manner. The free-energy landscape of protein folding can be analyzed in a large temperature range via the statistical analysis using a temperature adjustable reweighting factor (eq 3).17,57−59,61 Figure 2 indicates free-energy profiles spanned along a single folding reaction coordinate (RMSD) for the 6 proteins at various temperatures. It is noteworthy that for each protein, the free-energy profiles calculated from individual trajectories share the similar global features, and the calculated error bars are quite small compared to the detailed free-energy values, suggesting the convergence of the ITS-MD simulation. As shown in Figure 2, the folded (F) states are evidently weighted in the free-energy profiles in the temperature range of 280−360 K. Therefore, for suitable proteins, implicit solvent simulations can not only capture the folded structures from the vast conformational space but also assess their weights on the freeenergy surfaces. Besides the folded states, there are also multiple other distinct states, making the free-energy profiles rugged. The increase in temperature decreases the free-energy

varied atom coordinates and velocities. The simulation parameters of ITS-MD are shown in Table 1, and all relevant Table 1. Simulation Parameters for ITS-MD Simulations under Study protein BdpA EnHD EnHD_K52A α3W α3D S824

PDB code

temp range (K)a

no. of temp (N)b

ttraj (μs)c

NFd

1BDD 1ENH

280−430 280−430 280−430 280−430 280−430 280−430

200 350 350 300 400 350

4.9 5.8 4.7 10.4 9.0 6.8

24 27 25 18 14 11

1LQ7 2A3D 1P68

a

The temperature range {βi} integrated in eq 1. bNumber of temperature points in the desired range (N in eq 1). cThe accumulated simulation time. dNumber of folding (F) events observed in all trajectories.

trajectories represented by the time series of root-mean-square deviation (RMSD) and radius of gyration (Rg) are depicted in Figures S1−S6. The calculations of RMSD (with respect to protein native structures) and Rg are on backbone heavy atoms. We can see that the protein configurations are extensively

Figure 1. Comparison of the structures with lowest RMSDs (blue) in ITS-MD simulations initiated with extended conformations to the experimental native structures (red). Under each structure is shown the protein name, chain length, and backbone RMSD value. D


Article


Figure 2. One-dimensional free-energy profiles as the function of backbone RMSD with respect to the native structures for (A) BdpA, (B) EnHD, (C) EnHD_K52A, (D) α3W, (E) α3D, and (F) S824 at different temperatures (black: 280 K, red: 300 K, blue: 320 K, green: 340 K, purple: 360 K). The representative conformations of folded (F) states indicated in free-energy profiles are presented (colors code the secondary structural elements as defined in VMD software, and side chains involved in native hydrophobic contacts are shown with licorice representation).

Figure 3. One-dimensional free-energy profiles as the function of the fraction of native contacts (Q) for (A) BdpA, (B) EnHD, (C) EnHD_K52A, (D) α3W, (E) α3D, and (F) S824 at 300 K, along with the representative conformations of the distinct unfolded (U) states indicated in free-energy profiles.

process6,32 is used to depict one-dimensional free-energy profiles. The free-energy and structural characterization analyses of folding pathways hereinafter (Figures 3−8) are made at room-temperature (300 K). As shown in Figure 3, most of the proteins are folded via two-state transitions from the unfolded (U) to folded (F) states except for α3D and S824 whose folding pathways consist of intermediate (I) states as well (the representative conformations of the intermediates are depicted in Figure S8). For specific proteins (e.g., BdpA, α3D, and S824), the F states are not at the lowest-free-energy level, implying that these global native structures might be not thermodynamically favorable by the force field and implicit solvent model used here.

difference between the folded state and other states and smooths the free-energy surface in each case. We also calculated the heat capacity profiles (Cp ∝ T) for all proteins using eq 7 and extrapolated their melting temperatures (see the detailed description in the “Computational Methods” section). As shown in Figure S7, the calculated Tm values are not highly deviated from the experimental data, suggesting that the combination of the presently used force field and implicit solvent model could be a reasonable choice for the folding of these α-helix bundle proteins. 3.2. Compact and Native-like Topologies of Unfolded States in Implicit Solvent. To suitably demonstrate the folding pathway, the fraction of native contacts Q that is a better reaction coordinate than RMSD in describing the folding E


Article


proteins rich in α-helix structures.84 In contrast, a variety of small-angle X-ray scattering (SAXS) experiments either could not observe the early structural collapse after the initiation of refolding for proteins including protein L85,86 or observed unchanged Rg in a large range of denaturant concentration for reduced r-RNase A,87,88 implying that the unfolded states could adopt expanded conformations for specific proteins. Despite the controversy between experiments,76 on the other hand, it has been recognized that the unfolded states produced by MD simulations are usually more compact than the experimentally suggested ones for proteins larger than 20 amino acids.77,89−92 Such discrepancies between experiments and simulations are speculated to be induced by the general problem of either the functional form or the force field parametrization.93 Next, we tried to compare the unfolded states measured here to the counterparts simulated by the recent extensive explicit solvent MD of Lindorff-Larsen et al. using the CHARMM* force field.6 For EnHD and α3D, the two proteins simulated by both kinds of simulations, hierarchical clustering analyses were performed using the MMTSB toolset94 to determine the structure ensembles of their unfolded states, of which the 5 most populated clusters are shown in Figure 4. One can see that the unfolded states adopt various compact conformations. In addition, for each protein, some unfolded clusters share similar structural features among the simulations despite the use of different force fields (see the clusters connected by dashed lines in Figure 4A and B). Therefore, the present implicit solvent simulations can quantitatively detect unfolded states of proteins to some extent comparable to explicit solvent simulations.6 To further compare the implicit and explicit solvent simulations, we calculated the α-helix contents of individual helix segments along the folding pathway spanned with Q for EnHD and α3D, in which the helices are defined as Helix 1 to Helix 3 from the N- to C-terminus. As shown in Figure 4C, the formation of the individual helices is sequential in EnHD which follows the same order (Helix 1 ≈ Helix 2 > Helix 3) in implicit and explicit solvent simulations. The sequence of α-helix formation and their structural stability also follow the order of Helix 1 ≈ Helix 2 > Helix 3 for α3D modeled in implicit solvent and such a difference among the three α-helices is diminished in explicit solvent simulation (Figure 4D). The main difference between implicit and explicit solvent simulations is that all helices are largely folded even at Q ≈ 0 in implicit solvent whereas the helices are formed gradually along the folding pathway in explicit solvent, indicating the overestimate of secondary structure content in implicit solvent simulations that occurs at the very beginning of folding. The presence of native-like structural elements in the unfolded state is closely related to the folding mechanism. Figure 5 shows the formation of native-like local structures of individual residues in the unfolded state, using the method by Lindorff-Larsen et al. 6 (see the definition of residue “nativeness” in the “Computational Methods” section). In comparison to the extensive explicit solvent MD simulations by Lindorff-Larsen et al.6 and REMD simulations by Jiang and Wu,31 our results give more native-like structures. A general tendency can be seen that the regions corresponding to αhelices in the native structure are always more native-like than loop regions, consistent with the observation in the REMD simulations.31 This observation supports the speculation that the initiation sites in protein folding are preferentially located in regions which have a high propensity to form the native

Two-state folding behavior (U → F) has been reported for many single-domain proteins including those studied here (e.g., BdpA, EnHD, and EnHD_K52A).6,32−34 The structural properties of unfolded state thus attract much attention because of its fundamental importance in determining the mechanism of protein folding.37,72−76 Our simulations give a highly compact unfolded state containing considerable native αhelix secondary structure content for each protein (see the representative (most populated) conformations of the unfolded states in Figure 3). The compactness of the unfolded structure is also illuminated by quite similar Rg between the folded and unfolded states, and the considerable secondary structure content is indicated by the high percentage of residues forming an α-helix in the unfolded state (Table 2). Table 2. Summary of the Structural Properties of the Folded (F) and Unfolded (U) States of Proteins at 300 K protein

Rg,F (Å)a

Rg,U (Å)a

αU (%)b

BdpA EnHD EnHD_K52A α3W α3D S824

10.13 11.54 11.00 11.76 12.40 13.82

10.52 12.36 11.41 11.86 12.68 13.78

31.80 59.90 51.25 72.52 60.25 62.34

a

Average radius of gyration of F and U states. bAverage percentage of residues forming an α-helix in the U state.

Compact and native-like unfolded states of BdpA and EnHD, the two popular models often used in an in vitro protein folding investigation, have been experimentally suggested,35−37 but no experimental data is available for the unfolded states of other proteins studied here. For instance, the time-resolved fluorescence resonance energy transfer (FRET) experiment by Huang et al. suggested a compact unfolded (denatured) state of BdpA.35 The nuclear magnetic resonance (NMR) and solution X-ray scattering experiment by Mayor et al. suggested that the unfolded state of EnHD protein has an extensive native secondary structure and is significantly compact under native conditions.36,37 Specifically, the parameter measuring the compactness of protein, Rg, is 12.8 ± 0.2 Å for the native state and 18.8 ± 0.5 Å for the denatured state of EnHD in mildly denaturing environments (low pH and ionic strength).36 The Rg of EnHD measured here is quite similar for the folded state (Rg = 11.54 Å), but the data of the simulated unfolded state (Rg = 12.36 Å) is moderately smaller than the experimentally measured one. It is also noteworthy that the measured Rg values of protein unfolded states are at a very similar level as those measured by recent explicit solvent REMD simulations of Jiang and Wu (see the Rg values in Table 2 here and Table 1 in ref 31). The compactness of unfolded state is an active topic in protein folding. It has been speculated that under native conditions, e.g., right after the dilution of chemical denaturant in an in vitro refolding experiment, the unfolded protein chain tends to collapse rapidly to a more compact conformation.77,78 This assumption has been supported by the observations of single-molecule FRET experiments on multiple proteins (IgGbinding domain of protein L, the cold shock protein (Csp) from Thermotoga maritime, and RNase H):74,79−83 the Rg values of the unfolded states of these proteins are decreased with decreasing denaturant concentration. Samanta et al. further predicted that β-sheet proteins are far more collapsible than F


Article


Figure 4. Comparison between the present ITS-MD simulations and previous conventional MD simulations by Lindorff-Larsen et al.6 (A, B) The five most populated clusters of unfolded states of EnHD and α3D. In each part, the unfolded clusters measured from the ITS-MD and conventional MD are shown at the top and bottom, respectively. Under the representative structure of each cluster is shown the cluster order and the average population in the cluster ensemble. The similar clusters observed in both simulations are connected by dashed lines. (C, D) The α-helix contents of individual helix segments in proteins along Q for EnHD and α3D. In each part, the data from the ITS-MD and conventional MD is given at the top and bottom, respectively.

Figure 5. Average distance from the native structure in the unfolded state at 300 K (with a smaller value indicating more native-like for a residue) for (A) BdpA, (B) EnHD, (C) EnHD_K52A, (D) α3W, (E) α3D, and (F) S824. The bottom of each plot shows the secondary structures in the native state.

G


Article


Figure 6. (A−B) Average secondary structure content, the percentage of SASA (PSASA), and the fraction of native contacts (Q) of proteins. (C−D) The secondary structure content and the percentage of SASA as a function of the progress coordinate Q. All profiles are measured at 300 K.

structure in the unfolded state, that is, individual helix regions in α-helix bundle proteins.95 3.3. Unified Picture of the Folding of α-Helical Proteins in Implicit Solvent. A folding process includes the formation of secondary structures and the global structural collapse from a fully expanded chain structure to a compact native one. We investigated the interplay between the two events in the folding process that could reveal the detailed mechanism of protein folding. In this respect, the average secondary structure content (α-helix) and the percentage of solvent accessible surface area (PSASA) were calculated, respectively. The former was averaged over all α-helical residues, and the latter was for the entire proteins. Overall, all the proteins except BdpA have high average αhelix contents (Figure 6A). Accordingly, along the folding pathway spanned with Q, most proteins start with high α-helix contents, reflecting the easy α-helix formation in these proteins (Figure 6C). The relatively low average α-helix content of BdpA could be attributed to the difficult formation of such secondary structure at the early folding stage (Q < 0.2). On the other hand, all proteins have high average percentages of SASA (Figure 6B), most of which (except EnHD and α3W) have large SASA percentages at the very beginning of the folding process already (Figure 6D). Similarly, the values of another parameter measuring protein compactness, Rg, also drop quickly for all proteins except EnHD and α3W (Figure S9), suggesting the easy structural collapse in the most simulation cases. The consequence of the structural collapse is the formation of long-range side-chain contacts including native and nonnative ones. The fraction of the native side-chain contacts among individual helices that can evaluate the construction of the correct tertiary structure was calculated along Q for all proteins. As shown in Figure 7, from the beginning of folding (Q = 0), α3W, α3D, and S824 rapidly and significantly increase interhelix native contacts. Conversely, the remaining proteins only possess a small fraction of native interhelix contacts in their respective unfolded states (Q ≈ 0.2), and their native interhelix contacts are largely increased after Q > 0.3. Interestingly, it is α3W, α3D, and S824 that have the highest α-helix contents in the early stage of the folding process (Figure 6C). Therefore, it seems that the α-helix contents formed in the early folding stage could be somehow correlated to interhelix

Figure 7. Fraction of the native side-chain contacts among individual helices as a function of the progress coordinate Q at 300 K.

native contacts. These observations are to some extent consistent with a mechanism for protein folding: the formation of a subset of crucial long-range native contacts in the early stage of folding could largely contribute to the establishment of a native-like topology and the stabilization of the native secondary structure elements in the unfolded state.96 Combining Figures 3−7, a general scenario of the folding of α-helical single-domain proteins can be deduced: starting from a fully expanded structure, the structural collapse occurs concomitantly with the formation of local α-helix secondary structures to reach a compact unfolded state consisting of a portion of native structure content. The further folding involves barrier-crossing transition(s) that are essentially the local structure rearrangement to form more secondary structures and native tertiary contacts. For detailed proteins under study, the difference in folding pathways is mainly the level of forming secondary structures and native contacts in the early folding stage. Two prevailing models have been proposed to explain the folding mechanism of single-domain proteins, in terms of the interplay between secondary and tertiary structure formation in folding pathways. The nucleation-condensation mechanism97,98 suggests a concerted consolidation of secondary and tertiary interactions as a global protein collapses into a “nucleus” which serves as a template for the further condensation toward protein native structure. On the other hand, the framework model99,100 and related diffusion-collision mechanism101−104 postulate that secondary structures are formed incipiently and H


Article


Figure 8. α-Helix contents of individual helices along Q for (A) BdpA, (B) EnHD_K52A, (C) α3W, and (D) S826 at 300 K, respectively.

of BdpA, corresponding to the experimental observation that the isolated Helix 3 and Helix 2 fragments display higher αhelix propensities than the Helix 1 fragment in aqueous solution.49 In addition, the MC simulations by Yang et al. also observed that the helix fraction in the transition state ensemble (TSE) of BdpA follows the order of Helix 3 > Helix 2 > Helix 1.45 In the K52A mutant of EnHD (EnHD_K52A), the α-helix contents of Helix 1 and Helix 2 is roughly maintained, but the α-helix content of Helix 3 is significantly increased (Figure 8B) in comparison to wild-type EnHD (Figure 4C), explaining the experimental report that K52A mutation stabilizes wild-type EnHD.50 For all of the 6 proteins, a general tendency can be seen that it is the helix at the terminus that is most difficult to be formed (has the lowest α-helix content), although details of the helix at which terminus (N- or C-) may vary. A rather similar observation in the MC simulations on BdpA, Villin, and EnHD is that the TSEs of these three-α-helix bundle proteins always consist of two well-formed helices connecting with each other and a third less folded helix that is isolated from the helixturn-helix motif.45 3.4. Implicit Solvent Simulations Weakly Displaying the Relationship of α-Helix Propensity and Folding Mechanism. The balance of secondary and tertiary interactions depends crucially on the inherent propensities for secondary and tertiary structures. The preceding studies have proposed that the folding mechanism of a protein could slide from nucleation−condensation to framework model as the propensity for forming secondary structure increases.46,48 The helix propensity scale of individual amino acids105 can be used to assess the conformational preference of α-helix in protein. This parameter is accumulated for all residues in a protein, with a small value indicating a high α-helix forming tendency of the entire protein. As shown in Figure 9, the α-helix propensities of proteins follow the order of EnHD ≈ EnHD_K52A > BdpA > α3W ≈ α3D > S824, corresponding to experimentally suggested folding mechanisms shifting from the framework (EnHD)46−48 to the nucleation-condensation model (BdpA and α3D).38,39 No experimental data is available to explain the folding mechanism of α3W and S824. On the other hand, in the present study, the α-helix contents averaged over the entire simulations (Figure 6A) and the α-helix contents formed along

independently, and the formed secondary structure elements dock to form the tertiary structure. As shown in Figures 6 and 7, the percentage of α-helix content is generally larger than the percentage of interhelix native contacts, particularly in the early stages of the folding pathways. In addition, α3W, α3D, and S824 form considerable α-helix contents (>0.6, the highest values among the 6 proteins) and a large portion of native contacts (>0.2) at the very beginning of folding. The α-helix contents (0.5−0.6) and native contacts (≤0.1) of EnHD and EnHD_K52A in their early folding stage are at the middle and the smallest among the proteins, respectively. We can regard the folding of these proteins as a hybrid framework/ nucleation-condensation mechanism: the formation of both secondary and tertiary structure elements occurs at the early folding stage, and the formation of the secondary structure is more ready than the tertiary structure. On the other hand, BdpA has only a small fraction of α-helix content and native contacts in its early folding stage, and both structural elements increase synergetically along the folding pathway, suggesting a pure nucleation-condensation model for this protein. It is noteworthy that previous IR temperature jump and ultrafast fluorescence mixing experiment38 and explicit solvent MD simulation6 reported a nucleation-condensation folding mechanism for α3D; the experimental Φ analysis39 and multiple simulations including MD and Monte Carlo (MC)40−45 suggested a nucleation-condensation mechanism for BdpA. In contrast, the EnHD protein that contains a higher helical propensity has been proposed to slide from nucleationcondensation to frame mechanism.46−48 Therefore, while the present implicit solvent simulation technique can reveal an experimentally supported nucleation-condensation folding mechanism for BdpA, it might be difficult to totally distinguish between the nucleation-condensation and framework mechanisms for other cases such as α3D and EnHD. In comparison to the in vitro folding experiment, the in silico simulation implemented with the implicit solvent model is more likely to fold proteins via a hybrid framework/nucleation-condensation mechanism. The interplay of individual helix segments along the folding pathway was analyzed. As shown in Figure 8A, the formation of Helices 2 and 3 is apparently easier than Helix 1 in the folding I


Article


structure of the former protein whose stabilizing effects are amplified in implicit solvent condition.

4. DISCUSSION AND CONCLUSIONS As a representative of implicit solvent models, the generalized Born (GB) model has been widely adopted for molecular dynamics simulations, with trade-offs between resolution, accuracy, and computational speed. The weakness of GB models in structural modeling23,25 and folding thermodynamics evaluation26−29 has been pointed out by a great variety of simulation studies mainly on small peptide systems. A systematic investigation becomes very necessary to illuminate whether molecular simulations implemented with the current implicit solvent models are qualified to explore the folding pathway of large-size proteins, which represents a prevalent application area of implicit solvent simulations. In this work, the folding transitions of a variety of α-helical proteins with diverse sizes (46−102 AA) are investigated at the atomistic level using the recently developed AMBER FF14SBonlysc force field and GB-Neck2 implicit solvent model that have been reported to capture the native structures for structurally diverse proteins.16 The sophisticated ITS algorithm whose improved conformational sampling ability has been displayed in various biomolecular systems17,57−64 is utilized to make a sufficient sampling on protein configuration space. The detailed folding pathway is subsequently explored on the free-energy surface through statistical analysis of the weights of the conformations sampled in simulation. Six of nine test proteins can be folded by the present implicit solvent simulations, suggesting a possible tendency that proteins with high α-helix propensities are more likely to be folded in implicit solvent condition. The exceptional case is S824 which has relatively low α-helix propensity. Its successful folding could be probably induced by plenty of salt-bridges in its native structure which can be overstabilized by the implicit solvent model. To the best of our knowledge, it is the first report for the successful folding of protein with more than 100 amino acids in implicit solvent simulation. The detailed analyses on the 6 folded proteins provide energy landscape views for their folding pathways. An indispensable unfolded state is quantitatively identified for each protein, having compact native-like topology with partial formation of the secondary structure and tertiary contacts. Such structural features of unfolded states have been also suggested in the experiments for BdpA35 and EnHD36,37 and explicit solvent simulations for EnHD and α3D.6,31 A unified picture for the folding of α-helix bundle proteins is thus provided: along a single dominant route, the structural collapse and secondary structure formation occur concomitantly to reach a distinct unfolded state that serves as a “nucleus” for the further folding toward the protein native structure.6,110 Along such a folding pathway, the proportion and degree of formation of secondary and tertiary interactions are varied depending on the detailed protein systems, most likely producing a mixed framework/ nucleation-condensation mechanism. We notice that discrepancies certainly exist among implicit and explicit solvent simulations as well as experiments. First, although the compactness of the unfolded states produced by the present implicit solvent simulations is very similar to that in previous explicit solvent simulations6,31 (see Rg values in Table 2 here and Table 1 in ref 31 for EnHD and α3D, and the comparison of the unfolded clusters between the present simulations and ref 6 in Figure 5), the simulated unfolded states

Figure 9. Distribution of proteins in the plot of the helix propensity scale (that are accumulated for all amino acid residues within a protein).

the folding pathway (Figure 6C) are in the order of α3W > α3D ≈ S824 > EnHD ≈ EnHD_K52A > BdpA, corresponding to the simulated folding mechanisms shifting from the mixed framework/nucleation-condensation model of EnHD, EnHD_K52A, α3W, α3D, and S824 to the nucleation-condensation model of BdpA. Therefore, the α-helix propensities reflected by the present implicit solvent simulations are deviated from the real ones determined by amino acid sequences, leading to the difficulty in distinguishing between the nucleation-condensation and framework mechanisms for specific proteins. As a comparative test, we also tried to simulate the folding of three other α-helix bundle proteins (λ-repressor, Fap-NRα, and S836) using the same ITS-MD simulation protocols to further indicate the feasibility of using the AMBER FF14SBonlysc force field coupled with the GB-Neck2 implicit model to fold proteins with varied α-helix propensities. As shown in Figure 9, these three proteins have relatively low α-helix propensities. It is interesting to observe that the native structures of the three proteins cannot be folded by implicit solvent simulations (see the trajectories in Figures S10−S12 and the structures with lowest RMSDs sampled in simulations in Figure S13). Similarly, the recent REMD simulations using the same force field and implicit solvent model indicated that the incorrectly arranged structure with large RMSD but not the native structure of λ-repressor is much favored in the implicit solvent environment.16 Therefore, it is reasonable to speculate that the currently available implicit solvent simulation method is more suitable for the folding of proteins with high α-helix propensities with the exceptional case of S824. To figure out why S824 can be folded in implicit solvent but S836 is not although they have high sequence identity and share similar four-α-helix bundle native structures, we checked their sequences and the detailed intraprotein interactions in native structures. As shown in Figure S14, the sequence difference between S824 and S836 is mainly focused on 4 positions where S824 contains more charged amino acid residues. Accordingly, the native structure of S824 possesses more salt-bridges than S836. Specifically, an ensemble of saltbridges is formed in the motif of Helix 1-turn-Helix 2 that is contributed by residues of HIE20, Lys27, Asp28, HIE31, Asp32, and Asp34 (Figure S14). Previous studies indicated that GB implicit models overstabilize salt-bridges.106,107 In addition, it has been suggested that the α-helix stability is partially dependent on the relative population of native and/or nonnative salt-bridges.108,109 Therefore, the selective folding of S824 but not S836 by implicit solvent simulation should be attributed to the more salt-bridges contained in the native J


Article

Journal of Chemical Theory and Computation may be more compact than the experimentally suggested ones. For instance, the Rg value of the unfolded state of EnHD (12.36 Å) is smaller than the experimental data (18.8 ± 0.5 Å).36 The preceding studies have indicated a general tendency that the unfolded states of proteins larger than 20 amino acids produced by MD simulations are usually more compact than the experimentally measured ones,77,89,90 probably induced by the general problem with either the functional form or the force field parametrization.93 Second, the present implicit solvent simulations give more native-like structures in comparison to explicit solvent simulations (see the average percentage of residues forming α-helix in the unfolded states in Table 2 here and Table 1 in ref 31 and Table S2 in ref 6). The extensive explicit solvent MD simulations by Piana et al. revealed that the compactness of the unfolded state is supposed to correlate with the α-helix content within the state, with a more compact unfolded state containing a larger fraction of αhelix structures.93 On the other hand, a series of the comparative studies of implicit and explicit solvent MD simulations indicated the overabundance of α-helix secondary structures in the native states of various short peptides, which could be attributed to the errors in the electrostatic component of solvation free-energy (ΔGpol).23,25,106 Meanwhile, Chen and Brooks have displayed that besides the dominant electrostatic solvation free-energy, the approximate evaluation of the nonpolar solvation free-energy (ΔGnonpol) also has a severe limitation that overestimates the strength of pairwise nonpolar interactions.111 As shown in Figure 6D, the barrierless structural collapse occurs rapidly in the very early stage of folding to reach a global compact unfolded structure for each protein. Meanwhile, the α-helix secondary structure is also readily formed in the very early stage of folding. As a result, the overestimate of secondary structure content in the unfolded state in implicit solvent could be attributed to the errors in both the nonpolar and electrostatic components of solvation freeenergy. In summary, implicit solvent models offer an attractive alternative approach for studying protein folding in silico, which, however, may have high risk for providing convincing evidence to reveal the folding mechanism of protein. The present study shows a rush in both structural collapse and secondary structure formation in the early stage of folding in implicit solvent. As a result, the rate-limiting steps of protein folding become mainly the accomplishment of secondary structure elements and native contacts in compact protein structure. Therefore, although implicit solvent simulations can produce transient states to some extent similar to the counterparts by explicit solvent simulations, they are deviated from experimentally suggested ones. The compactness of transient states is supposed to be induced by the general problem of force field parametrization,93 and the overabundance of native-like structure elements along the folding pathway is mainly induced by the overestimate of secondary structure by the implicit solvent model.23,25,106 These results suggest that future efforts might benefit from codevelopment of implicit solvent models with force fields.

■

■

and average values of SASA (SASAavg) for proteins; simulation trajectories of 6 proteins; heat capacity profiles (Cp ∝ T); representative structures of the intermediate states of α3D and S824; average radius of gyration (Rg) along the fraction of native contacts Q; simulation trajectories of λ-repressor, Fap-NRα, and S836 proteins; comparison of the structures with lowest RMSDs in ITS-MD simulations to the experimental native structures of λ-repressor, Fap-NRα, and S836 proteins; comparison of amino acid sequence of S824 and S836 and salt-bridges formed by unconserved residues in the native structures (PDF)

AUTHOR INFORMATION

Corresponding Author

*Phone: +86 21 50806600-1304. E-mail: [email protected]. cn. ORCID

Qiang Shao: 0000-0001-6460-3095 Weiliang Zhu: 0000-0001-6699-5299 Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 21373258), National Basic Research Program (Grant No. 2014CB910400), and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) under Grant No. U1501501. The simulations were run at TianHe 1 supercomputer in Tianjin and TianHe II supercomputer in Guangzhou, China. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank D.E. Shaw Research for graciously sharing MD simulation data.

■

REFERENCES

(1) Dill, K. A.; MacCallum, J. L. The protein-folding problem, 50 years on. Science 2012, 338, 1042−1046. (2) Karplus, M. The levinthal paradox: Yesterday and today. Folding Des. 1997, 2, S69−S75. (3) Zwanzig, R.; Szabo, A.; Bagchi, B. Levinthals paradox. Proc. Natl. Acad. Sci. U. S. A. 1992, 89, 20−22. (4) Daggett, V.; Fersht, A. The present view of the mechanism of protein folding. Nat. Rev. Mol. Cell Biol. 2003, 4, 497−502. (5) Englander, S. W.; Mayne, L. The nature of protein folding pathways. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 15873−15880. (6) Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How fastfolding proteins fold. Science 2011, 334, 517−520. (7) Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H. F.; Shaw, D. E. Biomolecular simulation: A computational microscope for molecular biology. Annu. Rev. Biophys. 2012, 41, 429−452. (8) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990, 112, 6127−6129. (9) Edinger, S. R.; Cortis, C.; Shenkin, P. S.; Friesner, R. A. Solvation free energies of peptides: Comparison of approximate continuum solvation models with accurate solution of the poisson-boltzmann equation. J. Phys. Chem. B 1997, 101, 1190−1197. (10) Cramer, C. J.; Truhlar, D. G. Implicit solvation models: Equilibria, structure, spectra, and dynamics. Chem. Rev. 1999, 99, 2161−2200.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00726. Amino acid sequence of proteins; minimal (SASAmin) and maximal (SASAmax) values used for measuring PSASA K


Article

Journal of Chemical Theory and Computation (11) Im, W.; Feig, M.; Brooks, C. L. An implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteins. Biophys. J. 2003, 85, 2900−2918. (12) Spassov, V. Z.; Yan, L.; Szalma, S. Introducing an implicit membrane in generalized Born/solvent accessibility continuum solvent models. J. Phys. Chem. B 2002, 106, 8726−8738. (13) Feig, M.; Brooks, C. L. Recent advances in the development and application of implicit solvent models in biomolecule simulations. Curr. Opin. Struct. Biol. 2004, 14, 217−224. (14) Tanner, D. E.; Chan, K. Y.; Phillips, J. C.; Schulten, K. Parallel generalized Born implicit solvent calculations with NAMD. J. Chem. Theory Comput. 2011, 7, 3635−3642. (15) 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. (16) Nguyen, H.; Maier, J.; Huang, H.; Perrone, V.; Simmerling, C. Folding simulations for proteins with diverse topologies are accessible in days with a physics-based force field and implicit solvent. J. Am. Chem. Soc. 2014, 136, 13959−13962. (17) Shao, Q.; Shi, J.; Zhu, W. Determining protein folding pathway and associated energetics through partitioned integrated-temperingsampling simulation. J. Chem. Theory Comput. 2017, 13, 1229−1243. (18) Voelz, V. A.; Bowman, G. R.; Beauchamp, K.; Pande, V. S. Molecular simulation of ab initio protein folding for a millisecond folder NTL9(1−39). J. Am. Chem. Soc. 2010, 132, 1526−1527. (19) Zagrovic, B.; Snow, C. D.; Shirts, M. R.; Pande, V. S. Simulation of folding of a small α-helical protein in atomistic detail using worldwide-distributed computing. J. Mol. Biol. 2002, 323, 927−937. (20) Bursulaya, B. D.; Brooks, C. L. Comparative study of the folding free energy landscape of a three-stranded β-sheet protein with explicit and implicit solvent models. J. Phys. Chem. B 2000, 104, 12378−12383. (21) Bottaro, S.; Lindorff-Larsen, K.; Best, R. B. Variational optimization of an all-atom implicit solvent force field to match explicit solvent simulation data. J. Chem. Theory Comput. 2013, 9, 5641−5652. (22) Chen, J.; Im, W. P.; Brooks, C. L. Balancing solvation and intramolecular interactions: Toward a consistent generalized Born force field. J. Am. Chem. Soc. 2006, 128, 3728−3736. (23) Robinson, M. K.; Monroe, J. I.; Shell, M. S. Are AMBER force fields and implicit solvation models additive? A folding study with a balanced peptide test set. J. Chem. Theory Comput. 2016, 12, 5631− 5642. (24) Maffucci, I.; Contini, A. An updated test of AMBER force fields and implicit solvent models in predicting the secondary structure of helical, β-hairpin, and intrinsically disordered peptides. J. Chem. Theory Comput. 2016, 12, 714−727. (25) Roe, D. R.; Okur, A.; Wickstrom, L.; Hornak, V.; Simmerling, C. Secondary structure bias in generalized Born solvent models: Comparison of conformational ensembles and free energy of solvent polarization from explicit and implicit solvation. J. Phys. Chem. B 2007, 111, 1846−1857. (26) Zhou, R. H.; Berne, B. J. Can a continuum solvent model reproduce the free energy landscape of a β-hairpin folding in water? Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12777−12782. (27) Nymeyer, H.; Garcia, A. E. Simulation of the folding equilibrium of α-helical peptides: A comparison of the generalized Born approximation with explicit solvent. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 13934−13939. (28) Tan, C. H.; Yang, L. J.; Luo, R. How well does poissonboltzmann implicit solvent agree with explicit solvent? A quantitative analysis. J. Phys. Chem. B 2006, 110, 18680−18687. (29) Cumberworth, A.; Bui, J. M.; Gsponer, J. Free energies of solvation in the context of protein folding: Implications for implicit and explicit solvent models. J. Comput. Chem. 2016, 37, 629−640. (30) Jang, S.; Kim, E.; Pak, Y. Free energy surfaces of miniproteins with a ββα motif: Replica exchange molecular dynamics simulation with an implicit solvation model. Proteins: Struct., Funct., Genet. 2006, 62, 663−671.

(31) Jiang, F.; Wu, Y. D. Folding of fourteen small proteins with a residue-specific force field and replica-exchange molecular dynamics. J. Am. Chem. Soc. 2014, 136, 9536−9539. (32) Best, R. B.; Hummer, G.; Eaton, W. A. Native contacts determine protein folding mechanisms in atomistic simulations. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 17874−17879. (33) Munoz, V.; Campos, L. A.; Sadqi, M. Limited cooperativity in protein folding. Curr. Opin. Struct. Biol. 2016, 36, 58−66. (34) Baxa, M. C.; Freed, K. F.; Sosnick, T. R. Quantifying the structural requirements of the folding transition state of protein A and other systems. J. Mol. Biol. 2008, 381, 1362−1381. (35) Huang, F.; Lerner, E.; Sato, S.; Amir, D.; Haas, E.; Fersht, A. R. Time-resolved fluorescence resonance energy transfer study shows a compact denatured state of the B domain of protein A. Biochemistry 2009, 48, 3468−3476. (36) Mayor, U.; Grossmann, J. G.; Foster, N. W.; Freund, S. M. V.; Fersht, A. R. The denatured state of engrailed homeodomain under denaturing and native conditions. J. Mol. Biol. 2003, 333, 977−991. (37) Religa, T. L.; Markson, J. S.; Mayor, U.; Freund, S. M. V.; Fersht, A. R. Solution structure of a protein denatured state and folding intermediate. Nature 2005, 437, 1053−1056. (38) Zhu, Y.; Alonso, D. O. V.; Maki, K.; Huang, C. Y.; Lahr, S. J.; Daggett, V.; Roder, H.; DeGrado, W. F.; Gai, F. Ultrafast folding of α3D: A de novo designed three-helix bundle protein. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 15486−15491. (39) Sato, S.; Religa, T. L.; Daggett, V.; Fersht, A. R. Testing proteinfolding simulations by experiment: B domain of protein A. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 6952−6956. (40) Berriz, G. F.; Shakhnovich, E. I. Characterization of the folding kinetics of a three-helix bundle protein via a minimalist langevin model. J. Mol. Biol. 2001, 310, 673−685. (41) Garcia, A. E.; Onuchic, J. N. Folding a protein in a computer: An atomic description of the folding/unfolding of protein A. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 13898−13903. (42) Alonso, D. O. V.; Daggett, V. Staphylococcal protein A: Unfolding pathways, unfolded states, and differences between the B and E domains. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 133−138. (43) Ghosh, A.; Elber, R.; Scheraga, H. A. An atomically detailed study of the folding pathways of protein A with the stochastic difference equation. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 10394− 10398. (44) Favrin, G.; Irback, A.; Wallin, S. Folding of a small helical protein using hydrogen bonds and hydrophobicity forces. Proteins: Struct., Funct., Genet. 2002, 47, 99−105. (45) Yang, J. S.; Wallin, S.; Shakhnovich, E. I. Universality and diversity of folding mechanics for three-helix bundle proteins. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 895−900. (46) Gianni, S.; Guydosh, N. R.; Khan, F.; Caldas, T. D.; Mayor, U.; White, G. W. N.; DeMarco, M. L.; Daggett, V.; Fersht, A. R. Unifying features in protein folding mechanisms. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 13286−13291. (47) Mayor, U.; Guydosh, N. R.; Johnson, C. M.; Grossmann, J. G.; Sato, S.; Jas, G. S.; Freund, S. M. V.; Alonso, D. O. V.; Daggett, V.; Fersht, A. R. The complete folding pathway of a protein from nanoseconds to microseconds. Nature 2003, 421, 863−867. (48) Daggett, V.; Fersht, A. R. Is there a unifying mechanism for protein folding? Trends Biochem. Sci. 2003, 28, 18−25. (49) Bai, Y. W.; Karimi, A.; Dyson, H. J.; Wright, P. E. Absence of a stable intermediate on the folding pathway of protein A. Protein Sci. 1997, 6, 1449−1457. (50) Stollar, E. J.; Mayor, U.; Lovell, S. C.; Federici, L.; Freund, S. M. V.; Fersht, A. R.; Luisi, B. F. Crystal structures of engrailed homeodomain mutants - implications for stability and dynamics. J. Biol. Chem. 2003, 278, 43699−43708. (51) Zhang, C.; Ma, J. Folding helical proteins in explicit solvent using dihedral-biased tempering. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 8139−8144. L


Article

Journal of Chemical Theory and Computation (52) Nickson, A. A.; Wensley, B. G.; Clarke, J. Take home lessons from studies of related proteins. Curr. Opin. Struct. Biol. 2013, 23, 66− 74. (53) Gao, Y. Q. An integrate-over-temperature approach for enhanced sampling. J. Chem. Phys. 2008, 128, 064105. (54) Yang, L.; Liu, C. W.; Shao, Q.; Zhang, J.; Gao, Y. Q. From thermodynamics to kinetics: Enhanced sampling of rare events. Acc. Chem. Res. 2015, 48, 947−955. (55) Christ, C. D.; van Gunsteren, W. F. Enveloping distribution sampling: A method to calculate free energy differences from a single simulation. J. Chem. Phys. 2007, 126, 184110. (56) Lin, Z.; van Gunsteren, W. F. Enhanced conformational sampling using enveloping distribution sampling. J. Chem. Phys. 2013, 139, 144105. (57) Shao, Q.; Shi, J.; Zhu, W. Enhanced sampling molecular dynamics simulation captures experimentally suggested intermediate and unfolded states in the folding pathway of trp-cage miniprotein. J. Chem. Phys. 2012, 137, 125103. (58) Shao, Q.; Gao, Y. Q. Temperature dependence of hydrogenbond stability in β-hairpin structures. J. Chem. Theory Comput. 2010, 6, 3750−3760. (59) Shao, Q.; Gao, Y. Q. The relative helix and hydrogen bond stability in the B domain of protein A as revealed by integrated tempering sampling molecular dynamics simulation. J. Chem. Phys. 2011, 135, 135102. (60) Shao, Q.; Wang, J.; Shi, J.; Zhu, W. The universality of β-hairpin misfolding indicated by molecular dynamics simulations. J. Chem. Phys. 2013, 139, 165103. (61) Shao, Q.; Zhu, W.; Gao, Y. Q. Robustness in protein folding revealed by thermodynamics calculations. J. Phys. Chem. B 2012, 116, 13848−13856. (62) Shao, Q.; Xu, Z.; Wang, J.; Shi, J.; Zhu, W. Energetics and structural characterization of the ″DFG-flip’’ conformational transition of B-raf kinase: A SITS molecular dynamics study. Phys. Chem. Chem. Phys. 2017, 19, 1257−1267. (63) Shao, Q. Enhanced conformational sampling technique provides an energy landscape view of large-scale protein conformational transitions. Phys. Chem. Chem. Phys. 2016, 18, 29170−29182. (64) Yin, Y. D.; Yang, L.; Zheng, G. Q.; Gu, C.; Yi, C. Q.; He, C.; Gao, Y. Q.; Zhao, X. S. Dynamics of spontaneous flipping of a mismatched base in DNA duplex. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 8043−8048. (65) Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Gohlke, H.; et al. AMBER 14; University of California: San Francisco, CA, 2014. (66) Nguyen, H.; Roe, D. R.; Simmerling, C. Improved generalized Born solvent model parameters for protein simulations. J. Chem. Theory Comput. 2013, 9, 2020−2034. (67) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numericalintegration of cartesian equations of motion of a system with constraints - molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341. (68) Walker, R. Case study: All atom structure prediction and folding simulations of a stable protein (folding Trp-cage peptide). http:// ambermd.org/tutorials/basic/tutorial3/section4.htm (accessed Nov 15, 2017). (69) Garcia, A. E.; Sanbonmatsu, K. Y. Exploring the energy landscape of a β-hairpin in explicit solvent. Proteins: Struct., Funct., Genet. 2001, 42, 345−354. (70) Kabsch, W.; Sander, C. Dictionary of protein secondary structure: Pattern recognition of hydrogen-bonded and geometrical features. Biopolymers 1983, 22, 2577−2637. (71) Hubbard, S. J.; Thornton, J. M. Naccess, computer program; Department of Biochemistry and Molecular Biology, University College London: 1993. (72) Schuler, B.; Eaton, W. A. Protein folding studied by singlemolecule FRET. Curr. Opin. Struct. Biol. 2008, 18, 16−26.

(73) Mok, K. H.; Kuhn, L. T.; Goez, M.; Day, I. J.; Lin, J. C.; Andersen, N. H.; Hore, P. J. A pre-existing hydrophobic collapse in the unfolded state of an ultrafast folding protein. Nature 2007, 447, 106− 109. (74) Merchant, K. A.; Best, R. B.; Louis, J. M.; Gopich, I. V.; Eaton, W. A. Characterizing the unfolded states of proteins using singlemolecule FRET spectroscopy and molecular simulations. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 1528−1533. (75) Borgia, A.; Williams, P. M.; Clarke, J. Single-molecule studies of protein folding. Annu. Rev. Biochem. 2008, 77, 101−125. (76) Sosnick, T. R.; Barrick, D. The folding of single domain proteins - have we reached a consensus? Curr. Opin. Struct. Biol. 2011, 21, 12− 24. (77) Kohn, J. E.; Millett, I. S.; Jacob, J.; Zagrovic, B.; Dillon, T. M.; Cingel, N.; Dothager, R. S.; Seifert, S.; Thiyagarajan, P.; Sosnick, T. R.; et al. Random-coil behavior and the dimensions of chemically unfolded proteins. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 12491−12496. (78) Chan, H. S.; Dill, K. A. Polymer principles in protein-structure and stability. Annu. Rev. Biophys. Biophys. Chem. 1991, 20, 447−490. (79) Sherman, E.; Haran, G. Coil-globule transition in the denatured state of a small protein. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 11539−11543. (80) Nettels, D.; Gopich, I. V.; Hoffmann, A.; Schuler, B. Ultrafast dynamics of protein collapse from single-molecule photon statistics. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 2655−2660. (81) Nettels, D.; Muller-Spath, S.; Kuster, F.; Hofmann, H.; Haenni, D.; Ruegger, S.; Reymond, L.; Hoffmann, A.; Kubelka, J.; Heinz, B.; et al. Single-molecule spectroscopy of the temperature-induced collapse of unfolded proteins. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 20740−20745. (82) Ziv, G.; Haran, G. Protein folding, protein collapse, and tanford’s transfer model: Lessons from single-molecule FRET. J. Am. Chem. Soc. 2009, 131, 2942−2947. (83) Kuzmenkina, E. V.; Heyes, C. D.; Nienhaus, G. U. Singlemolecule FRET study of denaturant induced unfolding of RNase H. J. Mol. Biol. 2006, 357, 313−324. (84) Samanta, H. S.; Zhuravlev, P. I.; Hinczewski, M.; Hori, N.; Chakrabarti, S.; Thirumalai, D. Protein collapse is encoded in the folded state architecture. Soft Matter 2017, 13, 3622−3638. (85) Plaxco, K. W.; Millett, I. S.; Segel, D. J.; Doniach, S.; Baker, D. Chain collapse can occur concomitantly with the rate-limiting step in protein folding. Nat. Struct. Biol. 1999, 6, 554−556. (86) Jacob, J.; Krantz, B.; Dothager, R. S.; Thiyagarajan, P.; Sosnick, T. R. Early collapse is not an obligate step in protein folding. J. Mol. Biol. 2004, 338, 369−382. (87) Jacob, J.; Dothager, R. S.; Thiyagarajan, P.; Sosnick, T. R. Fully reduced ribonuclease A does not expand at high denaturant concentration or temperature. J. Mol. Biol. 2007, 367, 609−615. (88) Wang, Y.; Trewhelia, J.; Goldenberg, D. P. Small-angle X-ray scattering of reduced ribonuclease A: Effects of solution conditions and comparisons with a computational model of unfolded proteins. J. Mol. Biol. 2008, 377, 1576−1592. (89) Piana, S.; Klepeis, J. L.; Shaw, D. E. Assessing the accuracy of physical models used in protein-folding simulations: Quantitative evidence from long molecular dynamics simulations. Curr. Opin. Struct. Biol. 2014, 24, 98−105. (90) Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. How robust are protein folding simulations with respect to force field parameterization? Biophys. J. 2011, 100, L47−L49. (91) Best, R. B.; Zheng, W. W.; Mittal, J. Balanced protein-water interactions improve properties of disordered proteins and nonspecific protein association. J. Chem. Theory Comput. 2014, 10, 5113− 5124. (92) Piana, S.; Donchev, A. G.; Robustelli, P.; Shaw, D. E. Water dispersion interactions strongly influence simulated structural properties of disordered protein states. J. Phys. Chem. B 2015, 119, 5113− 5123. (93) Piana, S.; Lindorff-Larsen, K.; Dirks, R. M.; Salmon, J. K.; Dror, R. O.; Shaw, D. E. Evaluating the effects of cutoffs and treatment of M


Article

Journal of Chemical Theory and Computation long-range electrostatics in protein folding simulations. PLoS One 2012, 7, e39918. (94) Feig, M.; Karanicolas, J.; Brooks, C. L. MMTSB tool set: Enhanced sampling and multiscale modeling methods for applications in structural biology. J. Mol. Graphics Modell. 2004, 22, 377−395. (95) Modig, K.; Jurgensen, V. W.; Lindorff-Larsen, K.; Fieber, W.; Bohr, H. G.; Poulsen, F. M. Detection of initiation sites in protein folding of the four helix bundle ACBP by chemical shift analysis. FEBS Lett. 2007, 581, 4965−4971. (96) Lindorff-Larsen, K.; Rogen, P.; Paci, E.; Vendruscolo, M.; Dobson, C. M. Protein folding and the organization of the protein topology universe. Trends Biochem. Sci. 2005, 30, 13−19. (97) Fersht, A. R. Optimization of rates of protein folding: The nucleation-condensation mechanism and its implications. Proc. Natl. Acad. Sci. U. S. A. 1995, 92, 10869−10873. (98) Nolting, B.; Agard, D. A. How general is the nucleationcondensation mechanism? Proteins: Struct., Funct., Genet. 2008, 73, 754−764. (99) Kim, P. S.; Baldwin, R. L. Specific intermediates in the folding reactions of small proteins and the mechanism of protein folding. Annu. Rev. Biochem. 1982, 51, 459−489. (100) Kim, P. S.; Baldwin, R. L. Intermediates in the folding reactions of small proteins. Annu. Rev. Biochem. 1990, 59, 631−660. (101) Karplus, M.; Weaver, D. L. Protein folding dynamics. Nature 1976, 260, 404−406. (102) Karplus, M.; Weaver, D. L. Diffusion-collision model for protein folding. Biopolymers 1979, 18, 1421−1437. (103) Karplus, M.; Weaver, D. L. Protein-folding dynamics: The diffusion-collision model and experimental data. Protein Sci. 1994, 3, 650−668. (104) Zhou, Y. Q.; Karplus, M. Folding of a model three-helix bundle protein: A thermodynamic and kinetic analysis. J. Mol. Biol. 1999, 293, 917−951. (105) Pace, C. N.; Scholtz, J. M. A helix propensity scale based on experimental studies of peptides and proteins. Biophys. J. 1998, 75, 422−427. (106) Shell, M. S.; Ritterson, R.; Dill, K. A. A test on peptide stability of AMBER force fields with implicit solvation. J. Phys. Chem. B 2008, 112, 6878−6886. (107) Geney, R.; Layten, M.; Gomperts, R.; Hornak, V.; Simmerling, C. Investigation of salt bridge stability in a generalized born solvent model. J. Chem. Theory Comput. 2006, 2, 115−127. (108) Bierzynski, A.; Kim, P. S.; Baldwin, R. L. A salt bridge stabilizes the helix formed by isolated C-peptide of RNase A. Proc. Natl. Acad. Sci. U. S. A. 1982, 79, 2470−2474. (109) Shoemaker, K. R.; Kim, P. S.; Brems, D. N.; Marqusee, S.; York, E. J.; Chaiken, I. M.; Stewart, J. M.; Baldwin, R. L. Nature of the charged-group effect on the stability of the C-peptide helix. Proc. Natl. Acad. Sci. U. S. A. 1985, 82, 2349−2353. (110) Onufriev, A.; Case, D. A.; Bashford, D. Structural details, pathways, and energetics of unfolding apomyoglobin. J. Mol. Biol. 2003, 325, 555−567. (111) Chen, J.; Brooks, C. L. Implicit modeling of nonpolar solvation for simulating protein folding and conformational transitions. Phys. Chem. Chem. Phys. 2008, 10, 471−481.

N


How Well Can Implicit Solvent Simulations Explore Folding Pathways

Recommend Documents