Alterations of Nonconserved Residues Affect ... - ACS Publications

28 Sep 2015 - domain from human tenascin-C (TNfn3) was rationally redesigned for stability. Our study unveiled that charge−charge interactions reduc...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCB

Alterations of Nonconserved Residues Affect Protein Stability and Folding Dynamics through Charge−Charge Interactions Swarnendu Tripathi,†,‡,∥ Angel E. Garcìa,‡,§ and George I. Makhatadze*,†,‡ †

Department of Biology, ‡Center for Biotechnology and Interdisciplinary Studies, and §Department of Physics, Applied Physics, and Astronomy, Rensselaer Polytechnic Institute, Troy, New York 12180, United States ∥ Department of Physics, University of Houston, Houston, Texas 77204, United States S Supporting Information *

ABSTRACT: Charge−charge interactions play an important role in thermal stability of proteins. We employed an all-atom, native-topology-based model with non-native electrostatics to explore the interplay between folding dynamics and stability of TNfn3 (the third fibronectin type III domain from tenascinC). Our study elucidates the role of charge−charge interactions in modulating the folding energy landscape. In particular, we found that incorporation of explicit charge−charge interactions in the WT TNfn3 induces energetic frustration due to the presence of residual structure in the unfolded state. Moreover, optimization of the surface charge−charge interactions by altering the evolutionarily nonconserved residues not only increases the thermal stability (in agreement with previous experimental study) but also reduces the formation of residual structure and hence minimizes the energetic frustration along the folding route. We concluded that charge−charge interaction in the rationally designed TNfn3 plays an important role not only in enhancing the stability but also in assisting folding.

■

INTRODUCTION Many applications in biotechnology require higher thermal stability of proteins to maintain the native configuration.1,2 In the past decade, studies from several groups combining computational modeling and experiments have demonstrated that thermal stability of protein can be increased substantially by optimizing the surface charge−charge interactions without affecting their biological functions.3−8 However, one major drawback of the simplified model used in these studies for calculation of electrostatic interactions is that they ignore the effect of charge−charge interactions in the unfolded state ensemble of the protein. Also, results from several experimental studies have shown existence of residual structure (e.g., local structural order) in the unfolded state due to electrostatic interactions, which may also impact protein folding and stability.9−12 In some instances a destabilizing effect was observed when charged residues formed favorable interactions in the unfolded state.9,13 Nevertheless, experimentally it is challenging to characterize the unfolded state, as protein samples diverse conformations at this state. In this study, we used an all-atom structure-based model with explicit charge− charge interactions and demonstrated how the electrostatic energy of a rationally designed protein3 enhances the thermal stability by modulating the folding energy landscape. The energy landscape theory provides a generalizable framework that describes the ensemble nature of the protein folding process.14−16 Nevertheless, over the past years several studies have focused on hydrophobic interactions17−20 to gain a broader understanding of protein folding dynamics. However, the role of non-native electrostatic interactions on funneled © XXXX American Chemical Society

energy landscape remains elusive mainly because of its complexity,21 as the electrostatic interactions are both shortand long-range and can be either attractive or repulsive in nature. In this study, by incorporating explicit non-native electrostatic interactions in a native-topology-based all-atom model, we provided a qualitative description on the role of charge−charge interactions on protein folding dynamics and stability. The model system is based on the earlier experimental study by Strickler et al.,3 where the third fibronectin type III domain from human tenascin-C (TNfn3) was rationally redesigned for stability. Our study unveiled that charge−charge interactions reduces the formation of residual structure in the unfolded state of the designed TNfn3, which as a result minimizes the energetic frustration along the folding route. We concluded that optimization of the surface charges in the rationally designed TNfn3 not only enhances the thermal stability by making the electrostatic interactions more favorable in the native state than the unfolded state, in agreement with earlier study,3 but also assists the folding dynamics. The lower degree of evolutionary conservation of the charged residues further corroborates the idea that the sequence of TNfn3 perhaps evolved to favor biological functions over stability.

■

RESULTS Evolutionary Analysis of WT TNfn3. TNfn3 contains 92 amino acids with an immunoglobulin-like fold composed of

Received: September 1, 2015 Revised: September 26, 2015

A

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B seven antiparallel β-strands arranged in two sheets22 (Figure 1). Wild-type (WT) TNfn3 has a total of 27 charged amino acids

Table 1. Mutations of TNfn3 Based on the Study by Strickler et al.3 residue no. (position) 7 (strand A) 19 (strand B) 49a (strand C′) 89 (strand G) a

residue WT (e) Gln Leu Asp Thr

(0) (0) (−1) (0)

residue MUT (e) Lys Lys Lys Lys

(+1) (+1) (+1) (+1)

Charge reversal.

analysis based on more than 200 homologous sequences of WT TNfn3 using the evolutionary trace (ET) server as described in ref 23 (see the Methods section). In Figure 2A, residues of the WT TNfn3 are colored according to their evolutionary importance based on the real

Figure 1. Native structures and contact map of TNfn3. Panels A and B show the structures of WT TNfn3 (in green) with surface (in gray) and charges. In both figures charged residues are represented in licorice. Positively charged residues Arg and Lys are shown in blue, and negatively charged residues Asp and Glu are shown in red. In the WT protein there are 4 Arg, 5 Lys, 10 Asp, and 8 Glu residues. In panel B for MUT TNfn3, residues Gln7 (Q7), Leu19 (L19), Asp49 (D49), and Thr89 (T89) are shown in spheres with blue color to denote the mutations to Lys residue, according to Table 1. In each figure, the asterisk denotes the residue D49 that has charge reversal in MUT TNfn3. Panel C shows residue−residue contact map of TNfn3 based on the native atom−atom contacts. Contacts from the different segments of the proteins are denoted by different colors. The secondary structure of TNfn3 is shown at the side and top of the contact map plot. In the lower half of the figure, the native structure of TNfn3 is shown with the color changing from blue (N-terminal end) to red (C-terminal end). All the three-dimensional images of proteins in this figure and other figures of this paper are obtained using the Visual Molecular Dynamics (VMD) package.45

Figure 2. Evolutionary trace (ET) analysis of WT TNfn3 sequences. Panel A: WT TNfn3 structure is colored based on the real value evolutionary trace (rvET) scores from red (conserved) to purple (nonconserved). Evolutionary conserved residues are shown in a van der Waals representation. Nonconserved residues are shown in the cartoon representation only. Panel B: distributions of rvET scores of WT TNfn3 are shown for all the 92 residues (in gray) and the 27 charged residues (in black) separately. Arg, Lys, Asp, and Glu residues are considered as charged residues.

value ET (rvET) score (also see Figure S1). In Figure 2B, we plotted the distribution of the rvET score separately for the charged residues only and compared that with the distribution of all the residues of WT TNfn3. Although, several residues from both the core and surface of TNfn3 structures are found to be conserved (shown in van der Waals representation in Figure 2A); however, the majority of these residues are uncharged (see Figure 2B and Figure S2). For example, out of 24 residues of TNfn3 with rvET score < 8, only 3 residues are charged. Interestingly, the 4 residues of WT TNfn3 (Gln7, Leu19, Asp49, and Thr89) that were mutated to Lys (Table 1) and exhibited significant increase in stability of the protein3 are

(9 positively charged (Arg and Lys) and 18 negatively charged (Asp and Glu) amino acids) (see Figure 1A) with a net charge of −9e. The rationally designed MUT3 of TNfn3 (Figure 1B) has three additional positive charges (Gln7Lys, Leu19Lys, and Thr89Lys) and one charge reversal (Asp49Lys), and the net charge is −4e (see Table 1), under normal physiological conditions (pH ∼ 6−8). To explore the degree of conservation of the charged residues of TNfn3, we performed evolutionary B

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. Global structural characteristics of TNfn3. Panels A−C show the distributions of the all (heavy) atom radius of gyration, Rg, end-to-end distance, and all (heavy)-atom root-mean-square deviation (RMSD) for the WTpf, WT, and MUT systems. In panel B the end-to-end distances are calculated between residues 5 and 88. All the plots are made at the corresponding TF of each system.

way than the Rg and the end-to-end distance for all the three systems. Specifically, we noticed that although the RMSD shows a two-state-like behavior for WTpf system, however, both WT and MUT clearly indicate the presence of an intermediate (I) state at ∼1.5 nm. Altogether, these analyses revealed that introduction of the charge−charge interactions into the structure based model have similar effect on the global structural properties of both WT and MUT systems of TNfn3; i.e., the RMSD analysis reveals the presence of an intermediate state, I. As it will be shown below, the presence of the I state is identifiable in the WTpf as well, when different reaction coordinate is used. Thermodynamic Properties: MUT TNfn3 Is the Most Stable Protein. To explore the role of charge−charge interactions on the thermodynamic properties of TNfn3, we compared the heat capacity, Cv, for all the three systems in Figure 4A. We observed that relative to the WTpf system there is a significant shift in the Cv peak for both WT and MUT systems. The folding transition temperatures (TF) from the Cv peaks for the WTpf, WT, and MUT systems are found to be at

evolutionarily not conserved (rvET score > 10). On an average, we found that the degree of conservation of the uncharged residues of WT TNfn3 is higher (average rvET score ∼ 10.4) than the charged residues (average rVET score of positively and negatively charged residues are ∼ 14.2 and 12.8, respectively). Effect of Charge−Charge Interactions on Global Properties. We studied the folding dynamics of TNfn3 using a pure-funneled (Go̅-like) structure based model and allatom representation combined with the Debye−Hückel (DH) equation to incorporate interactions between charged residues. Specifically, we studied the three following systems: (i) WTpf (pure-funneled model for wild type) TNfn3, (ii) WT (purefunneled model with charge−charge interactions for the wild type) TNfn3, and (iii) MUT (pure-funneled model with charge−charge interactions for the rationally designed mutant3) TNfn3. From our previous study,24 we found that the regular molecular dynamics (MD) simulations for the WTpf TNfn3 show slow convergence, and the system often gets trapped in an intermediate state. Therefore, for all the systems in this study, we performed extensive replica exchange molecular dynamics (REMD) simulations to overcome this problem (see the Methods section). Structural Properties: Presence of an Intermediate State. We first analyzed the role of the charge−charge interactions on the overall structural properties of the different systems of TNfn3. Our results showed that the (all heavyatom) radius of gyration (Rg) and the end-to-end distance of all the three systems exhibit a two-state-like folding behavior (Figures 3A and 3B). Additionally, we found that the charge− charge interactions have a significant effect on the configuration of the unfolded (U) state ensemble, as shown in Figure 3. Specifically, both Rg and end-to-end distance show shift in the population of the U state ensemble toward the lower values for both WT and MUT systems compared to WTpf. For the U state Rg of WT and MUT indicates almost 30% decrease in the peak value compared to WTpf TNfn3 (from ∼2.5 nm for WTpf to ∼1.75 nm for both WT and MUT), whereas for the native (N) state the peaks are found at ∼1.37 nm for all the three systems (Figure 3A). Similarly, at the U state the end-to-end distance also shows nearly 33% decrease for both WT and MUT compared to WTpf TNfn3 (from ∼6 nm for WTpf to ∼4 nm for both WT and MUT), while for the N state the peak is at ∼1.5 nm for all the three systems. These results clearly indicate the effect of charge−charge interactions on the overall increase in the compactness of the WT and MUT TNfn3 in the U state ensemble, probably due to the presence of residual structure. Interestingly, we found that the distribution of the (all heavyatom) root-mean-square deviation (RMSD) acts in a different

Figure 4. Folding thermodynamics of TNfn3. Panel A: heat capacities at constant volume (Cv) for the WTpf, WT, and MUT systems at different temperatures, T* ≡ kBT/εC (in the reduced unit). Panel B: free energy curves F(Q) are plotted as a function of the global reaction coordinate, Q, for the WTpf, WT, and MUT systems at T ∼ TF of WT. C

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. Electrostatic energy of TNfn3. Panel A shows the plots of the distribution of electrostatic energy, VDH, in the reduced unit, εC, for the WT and MUT TNfn3, respectively at their corresponding T ∼ TF. Panels B and C show two-dimensional map of VDH vs the reaction coordinate Q, where the color represents the joint probability distribution P(VDH,Q). In the color bar, red and blue indicate highest and lowest probability, respectively.

Figure 6. Folding of the different segments of TNfn3. Panels A−C illustrate the formation of the partitioned fraction of native contacts (Qp) for different segments of the protein compared to the average Q for the WTpf, WT, and MUT, respectively, at their corresponding T ∼ TF. Qp for different segments of the protein is indicated by different colors according to the contact map in Figure 1. In each plot, the diagonal line is drawn as a reference for native contacts that ideally follow the average Q and small arrows indicate the significant change in the segments BE and FG. The shaded regions in the plots represent the folded state ensemble (0.6 ≤ Q ≤ 1.0), where the changes are minimal.

0.938, 0.944, and 0.952, respectively, in the reduced unit T* ≡ kBT/εC. In terms of GROMACS simulation temperature this indicates 3 K increases in TF for WT relative to WTpf and an additional 4 K increase in TF for MUT relative to WT. This observation clearly validates that for TNfn3 the charge−charge interactions play a crucial role on the thermal stability. Specifically, the designed MUT shows significant increase in the stability than WT based on the optimization of surface charge−charge interactions, which is in good agreement with the experimental study by Strickler et al.3 Additionally, from the Cv plots (Figure 4A) we observed that the folding of both WT and MUT is significantly less cooperative than WTpf. This is clearly noticeable from the decrease in the height of the Cv curves for both WT and MUT systems than WTpf TNfn3 (Figure 4A). To characterize the change in stability of TNfn3 due to the charge−charge interaction, we plotted the free energy curves F(Q) for all the three systems as a function of the reaction coordinate Q (fraction of total native contacts as described in the Methods section) at T ∼ TF (determined from the peak of the Cv vs T plots in Figure 4A) of WTpf in Figure 4B. We found that for these three systems the N state is located at Q ∼ 0.8. However, for the U state there is a shift in Q value. The Q for WTpf is at ∼0, and due to the electrostatic interactions the Q value increases to ∼0.05 and 0.1 for WT and MUT, respectively. Increase in Q value of WT and MUT is consistent with our structural analyses in the previous section, which show that charge−charge interactions in the U state may induce formation of residual structure and therefore causes compactness in WT and MUT compared to WTpf. Additionally, this analyses revealed that for all the three systems there is an I state at Q ∼ 0.4. Interestingly, from our analyses in the previous

section we did not observe any I state for WTpf in terms of the overall conformation by calculating the Rg, end-to-end distance, or RMSD. Nevertheless, from F(Q) (Figure 4B) we found that the I states in both WT and MUT are significantly more populated than WTpf due to charge−charge interactions. Moreover, increase in the population of the I state in both WT and MUT makes these systems less cooperative than WTpf, which is noticeable from the Cv curves in Figure 4B, as discussed above. From the free energy difference ΔΔFU−N = ΔFU − ΔFN between the U and N states in Figure 4B, we estimated that for WT there is ∼4ϵC increase in the stability relative to WTpf. Similarly, we found ∼5ϵC increase in the stability of MUT relative to WT. This result again confirms that optimization of the surface charge−charge interactions enhances the thermal stability of MUT TNfn3. Electrostatic Energy and Thermal Stability. Next we explored the effect of electrostatic energy, VDH, on the thermal stability of WT and MUT TNfn3. In Figure 5A, we plotted the distribution of V DH for both these systems at their corresponding TF. The distributions clearly indicate an overall shift in VDH for MUT system toward more negative value than WT. By comparing the peak values (indicated by the dashed lines in Figure 5A) of the two distributions, we found ∼−5ϵC change in VDH of MUT TNfn3 compared to WT. To further understand how the charge−charge interactions affect the different states of TNfn3, in Figures 5B and 5C we plotted the joint probability distribution P(VDH,Q) for WT and MUT TNfn3, respectively. We found that the broad distribution of WT in Figure 5A corresponds to the U state ensemble (see Figure 5B), for which VDH is at ∼−5ϵC. On the other hand, VDH at ∼−10ϵC corresponds to the I and N states of WT (Figures 5A and 5B). In addition, we observed that D

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 7. Dynamical cross-correlation map (DCCM) of TNfn3. Panels A−C illustrate the DCCM calculated (using eq 7; see Methods section) over the unfolded (U) state ensemble of the WTpf, WT, and MUT systems, respectively, only for the Cα atoms. Similarly, panels D−F and G−I show the DCCM for the intermediate (I) and native (N) state ensemble, respectively, for the WTpf, WT, and MUT systems. In the color bars red indicates anticorrelated motions with negative values, blue indicates correlated motions with positive values, and greens indicate uncorrelated motions with zero values.

analysis, we first divided the native contacts of the protein in different segments as shown in Figure 1C. Note that we have the same native contacts for all the three systems (see the Methods). In TNfn3 the contacts in the segments AB, BE, CC′, CF, and FG form the β-sheets, whereas the contacts in AF and BG enclose the protein to provide a compact native structure. In Figure 6, we plotted the partitioned fraction of native contacts, Qp(Q), for each segment of the protein as a function of the total fraction of native contacts, Q, for all the three systems. We are particularly interested to see how different parts of the protein fold along the global reaction coordinate Q. As indicated in Figures 6A−C, we mainly focused on the range from 0 ≤ Q ≤ 0.6 for each system, where the changes along the folding route are noteworthy. First, we noted that for WTpf there is a significant backtracking in the BE segment. The backtracking in BE can be seen from the transient increase in Qp at Q ∼ 0.2 and then decrease at Q ∼ 0.4 (see Figure 6A). As we mentioned in our earlier study,24 the backtracking in BE occurs through breaking some of the native contacts between strands B and E mainly to

populations of both the I and N states are much higher with P(VDH,Q) ≥ 0.6 than the U state for which P(VDH,Q) ≤ 0.4 (Figure 5B). This further indicates that charge−charge interactions are more favorable in both the I and N states of WT TNfn3 compared to its U state. Interestingly, for MUT TNfn3 we noted that not only the VDH becomes more negative for all the three states (U, I, and N) compared to the VDH of WT (Figures 5A and 5C), but the N state is the most populated state with P(VDH,Q) ≥ 0.6 (Figure 5B). On the contrary, for both the U and I states P(VDH,Q) is ≤0.4 for MUT (Figure 5C). This result implies that charge−charge interaction is highly favorable in the N state of the rationally designed MUT TNfn3. Energetic Frustrations and Backtracking. We also aimed to investigate the role of charge−charge interactions on the overall folding dynamics of TNfn3. In our earlier study detailed characterization of the folding mechanism of WTpf TNfn3 in the absence of electrostatics was performed.24 Therefore, in this study we especially focused on the folding route analysis of WT and MUT systems. To perform this E

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

indicates that formation residual structure in the U state of these systems (WT and MUT) is perturbed by the charge− charge interactions. On the other hand, the segment AB and the CC′ loop show positive correlations between them (enclosed by small dashed circles in Figures 7A−C); however, the magnitude of cij is low (∼ 0.2), especially for WT and MUT systems because of the charge−charge interactions. At the I state, both positive and negative correlations are observed in different parts of the protein for all the three systems (Figures 7D−F). In particular, we found strong positive correlation within the CC′, CF, and FG segments, which is in fact associated with the formation of the native contacts in these segments (see Figure 1C). By comparing the cij in Figures 7D−F, we noted that the positive correlation in FG is higher (cij = 0.8) in WTpf than both WT and MUT systems for which cij is ∼0.6 for most of the residue pairs in this segment (as indicated by the dashed circles). This is particularly true for WT, for which backtracking is more prominent in the FG segment than MUT (Figures 6B and 6C), and as a result, the correlation is weaker. Additionally, our results showed strong negative correlation (cij ∼ −0.6 to −0.8) between the segment CC′EFG and the segment AB for all the three systems in the I state. This observation is particularly in agreement with the fact that at the I state the native contacts in the segments AB and BE are not formed (see Figure 6). For all the three systems in the N state, we first noted that segments AB and BE also show positive correlation (cij ∼ 0.4; Figures 7G−I) due to the formation of β sheets in these segments (Figure 1C). Additionally, we observed positive correlation between segment CC′ and strand G in WT and MUT systems that is more stronger than WTpf TNfn3, even though there are no native contacts exist (Figure 1C) in these segments. Stronger (negative) correlation (cij ∼ −0.4) is also observed for the segment BC, strand E, and segment FG with strand A in WT and MUT compared to WTpf TNfn3. On the other hand, the positive correlation in segment AB is found to be weaker in WT and MUT than WTpf TNfn3. These analyses revealed that charge−charge interactions induce higher degree of correlation in both the N state of WT and MUT TNfn3 compared to the WTpf, although the magnitude of correlation is weaker in the N state than the U and I states for the three systems.

facilitate the formation of the contacts in CF because of the topological frustration in this protein. For both WT and MUT we observed that the segment BE has significantly more native contacts (Qp of BE is ∼0.3 in Figures 6B and 6C, respectively) than WTpf (Qp ∼ 0.2 for BE in Figure 6A) in the range from 0 ≤ Q ≤ 0.2 before backtracking, which is particularly true for the WT system. This further suggests that charge−charge interactions in TNfn3 have a significant effect to induce formation of native-like residual structure between strands B and E in the U state. Additionally, we found that the segment BE shows less backtracking in both WT and MUT systems than WTpf system (indicated by an arrow in Figures 6A−C). For WTpf the change in Qp value of the segment BE decreases from ∼0.2 to 0.05, whereas for WT and MUT the change in Qp is from ∼0.3 to 0.2. As a result, at the range Q ∼ 0.4−0.5, Qp of BE increases significantly from ∼0.05 to 0.35 for WTpf TNfn3. On the other hand, for the same range of Q, the change in Qp of BE for WT and MUT TNfn3 is small (from ∼0.2 to 0.3). We further observed that other than the segment BE residual structure also forms in the segment CF due to the charge− charge interactions in the U state (indicated by arrows in Figures 6B and 6C at 0.0 ≤ Q ≤ 0.3 for WT and MUT, respectively), specifically for WT. As a result, we also noticed backtracking for the contacts in FG segment at 0.0 ≤ Q ≤ 0.3 for both WT and MUT (indicated by arrows in Figures 6B and 6C). Additionally, we found that backtracking in FG is more pronounced in WT than MUT. On the other hand, backtracking is not observed for this segment in WTpf (see the arrow in Figure 6A). As a result, at Q ∼ 0.2, Qp of FG for WTpf is much higher (∼0.5) than WT and MUT for which Qp of FG is ∼0.3 and ∼0.4, respectively (Figures 6A−C). This indicates that charge−charge interactions cause energetic frustration in the FG segment for both WT and MUT. Nonetheless, the optimization of the surface charges of MUT minimizes the energetic frustration in FG compared to WT (indicated by arrows in Figures 6B,C) by reducing the formation of residual structure in this segment. Because of the backtracking in FG, there is an early initiation of the native contacts in CF for WT and MUT, along the folding route, which is also understood from the crossing of the CF and FG curves in Figures 6B and 6C. Finally, we noted that the charge− charge interactions also affect the contacts in the BG segment. For WT and MUT systems Qp is ∼0.2 for BG at Q ∼ 0.5 (Figures 6B,C), whereas for WTpf TNfn3 the corresponding value of Qp is ∼0.3 (see Figure 6A). Effect of Charge−Charge Interactions on Folding Dynamics. In the earlier section we illustrated the effect of charge−charge interactions on the folding mechanism of the three different systems of TNfn3 using order parameter, Qi, which is based on the native interactions. However, it is also interesting to explore the effect of electrostatics on both the native and non-native interactions together. Therefore, we performed a cross-correlation analysis of the fluctuations at the different states (U, I, and N states) of the three systems (WTpf, WT, and MUT) of TNfn3 by considering the Cα atom pairs from the protein backbone. First, we noted that for all the three systems in the U state there is a negative correlation for the fluctuations between the segments ABCC′ and EFG (Figures 7A−C). This correlation is more prominent between the CC′ and FG segments for which the cross-correlation coefficient cij is ∼−0.8. Additionally, we found that, overall, the negative correlation in the U state is stronger in WTpf than WT and MUT (shown by the dashed circles in Figures 7A−C). This

■

DISCUSSION Experimental studies by Clarke and co-workers25 have shown that the 90-residue long TNfn3 (PDB ID: 1TEN) protein is reasonably stable (∼ 5 kcal/mol at 20 °C). Later, they confirmed that the thermal stability of TNfn3 could be increased by ∼3 kcal/mol with two amino acid (Gly-Lys) extension at the C-terminus.26 Interestingly, a combined computational and experimental study3 unveiled that the stability of the 92-residue long TNfn3 (which is denoted as WT TNfn3 in this study) can be enhanced further (with an increase of ∼1.3 kcal/mol at 25 °C) by rationally modifying the surface charge−charge interactions. In the study by Strickler et al.,3 first the sequence of WT TNfn3 was computationally redesigned using the genetic algorithm.27 The particular sequence that exhibited largest increase in energy from favorable charge−charge interactions due to the fewest number of substitutions was opted for experimental validation. For the selected redesigned protein (which is denoted as MUT TNfn3 in this study) there were total four substitutions compared to the WT TNfn3, including F

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B a charge reversal (Table 1). The study by Strickler et al.3 had two major outcomes: First, they confirmed the importance of surface charge−charge interactions on protein stability, and second, they showed that protein stability could be increased by rational optimization of charge−charge interactions. Nevertheless, the effect of residual charge−charge interactions in the unfolded state of the protein was not considered for the calculations of the energies in their simplified model. In this study, our evolutionary analysis of WT TNfn3 sequence confirmed that the four residues (Gln7, Leu19, Asp49, and Thr89) that were mutated to Lys in the designed MUT TNfn3 to enhance the thermal stability are not conserved. Additionally, evolutionary analysis revealed that the uncharged residues have lower degree of conservation compared to the charged residues, indicating relatively late evolution for the latter. Interestingly, the residues Leu19 and Asp49 of WT TNfn3 that were mutated to Lys in MUT TNfn3 show presence of Lys in sequence of TNfn3 from other species in the evolutionary history (Figure S1). Recently, several groups have shown interest in characterizing the role of charge−charge interactions in protein folding dynamics and stability using theoretical and computational approaches.12,28−33 Azia and Levy29 used a two-bead SBM of proteins with DH term to represent the non-native charge− charge interactions and showed that by placing the charge of the side-chain beads to the farthest heavy atom revealed better correlation between calculated and measured thermal stability from experiment. Naganathan et al.31 used an Ising-like statistical mechanical model to study the effects of two stabilizing charge mutations on the stability of barstar. A recent study by Tzul et al.30 using the Cα-based SBM for the protein with DH term for the electrostatics have shown that this simple model is able to predict the experimental thermodynamics and kinetics of proteins qualitatively. In this study, we implemented an all (heavy)-atom SBM with DH term to represent a more realistic protein model where charges are placed on the side chain (see Methods section). From an experimental study the 92-residue-long WT TNfn3 was predicted to be a two-state folder.26 However, using an allatom structure-based, pure-funnel model (WTpf without charges) in our previous study24 and including the DH term for the charge−charge interactions with pure-funnel model in the present study, we showed that the 92-residue-long WT TNfn3 is a three-state folder. It is important to note that the two-state or three-state behavior of TNfn3 is highly dependent on the choice of the order parameter. As our analyses illustrated that the Rg and end-to-end distance indicate two-state behavior for all the three systems, whereas the RMSD shows two-state behavior only for WTpf and three-state behavior for both WT and MUT TNfn3. On the other hand, the PMF F(Q) revealed three states for all the systems of TNfn3. In this study we further showed that the charge−charge interactions do not affect the global structural properties of WT and MUT TNfn3 at the N state but makes the U state and the I state more compact than WTpf TNfn3 due to the formation of residual structure through charge−charge interactions. Differences in the probability of native contact formation also demonstrates that in the U state charge−charge interactions induce formation of residual structure in the WT TNfn3, specifically in the segments BE and CF compared to the WTpf (Figure S2A). On the other hand, formation of residual structure is found particularly in the segments CC′ and FG in the MUT TNfn3 compared to the WT. Additionally, we noted that the charge−

charge interactions also affect the native contacts at the I and N states of WT and MUT TNfn3 (Figure S2B,C for WT; Figure S2E,F for MUT). However, based on the magnitude of the change in probability of native contact formation, the effect of charge−charge interactions in the N state is found to be less significant compared to the U and I states in both WT and MUT systems. Nonetheless, we found that the electrostatic energy in the N state of the rationally designed MUT TNfn3 is highly favorable than WT TNfn3 (Figures 5B and 5C, respectively). As a result, we noticed further increase in the thermal stability of MUT TNfn3 compared to WT (Figure 4) in good qualitative agreement with the experimental results.3 Frustration in the protein energy landscape has two components: topological and energetic.34−37 In structurebased (pure-funnel) model, where the native interactions are energetically favorable, only topological frustration exists because of the geometric constraints of protein chain. However, incorporation of electrostatics in a structure-based model may induce energetic frustration due to the non-native charge− charge interactions. Comparing the folding route analysis of WTpf and WT TNfn3 in the present study, we showed that charge−charge interactions not only cause energetic frustration through backtracking, but also result in small but detectable formation of native-like residual structure in the U state. Interestingly, we further observed that optimization of the surface charges of TNfn3 significantly reduces the energetic frustration in the MUT TNfn3 compared to the WT. This reduction of energetic frustration upon optimization of charge− charge interactions was previously demonstrated by us using Cα-structure-based model30 and is also in agreement with recent computer simulations of Naganathan et al.31 We also explored the interrelationship between charge− charge interactions and folding mechanism from the correlated fluctuations of the Cα atoms at different states for all the three systems of TNfn3. The correlation map captures both the native and non-native interactions from local fluctuations of Cα positions of each residue. Our study revealed that at the U state (Figure 7) there is a strong anticorrelated motions between the N-terminal and the C-terminal parts for all the TNfn3 systems (Figures 7A−C) and formation of residual structure due to electrostatics reduces the anticorrelated motion in the WT and MU TNfn3 systems. Additionally, we noted that at the I-state strong (positive) correlations are mostly associated with the formation of native contacts. It should be also noted that at this state the native contacts in the segments AB and BE were not formed. We conclude that significant backtracking in the segment BE impedes the formation of native contacts between strands B and E in the I state. Besides, the BE segment charge− charge interactions also impede formation of native contacts in the FG segment through backtracking. Nonetheless, in the case of MUT TNfn3 the optimization of the surface charge reduces the energetic frustration compared to the WT. For the N states of all the three systems some (positive) correlation is also observed between the strands G and C that does not have any native interactions (see Figure 1C), specifically for WT and MUT systems.

■

CONCLUSIONS In the present study, we used an all-atom SBM with electrostatics to illustrate the effect of charge−charge interactions on both the folding dynamics and thermal stability of TNfn3. Our simulation results demonstrate that in the case of WT TNfn3 the incorporation of electrostatics increases the G

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

make the native interactions not too dominating over the electrostatic interactions in our structure-based model. To incorporate the electrostatic effects we used the Debye− Hückel (DH) model following previous studies by other groups12,29 to define the second term in eq 1. The DH potential (VDH) takes into account both the dielectric effect of the solvent and the screening due to the charge−charge interactions by mobile ions in the solvent. VDH is defined as

thermal stability but causes backtracking due to the energetic frustration in the unfolded state. The energetic frustration evolves through the formation of native-like residual structure in the unfolded state that impedes the folding process. Additionally, our study revealed that for a redesigned TNfn3 protein (MUT TNfn3) the optimization of surface charges not only increases the thermal stability, in agreement with previous experiments,3 but also minimizes the energetic frustration by reducing the formation of residual structure, which also agrees with experimental observation that MUT TNfn3 folds faster than the WT TNfn3.30

VDH = kelecB(λD) ∑

■

qiqje−rij / λD

i,j

METHODS Protein Systems. The structure of the WT protein was modeled from the known structure of TNfn3 (PDB ID 1TEN) using Modeler version 7.7.38 The residue number 802-891 from the original PDB file (1TEN) of TNfn3 was renumbered from 1 to 90. Following the work by Hamill et al.,26 a 92-residue protein that was extended at the C-terminal end by the two residues Gly91 and Leu92 was used as a WT TNfn3. Additionally, the incomplete residue Arg1 was modeled in the WT TNfn3. In this study we kept all the interactions of the WTpf, WT, and MUT TNfn3 identical besides the electrostatics, and only the charges in the side chain of the mutated residues were added in the MUT system (Gln7Lys, Leu19Lys, Asp49Lys, and Thr89Lys; see Table 1) without performing the actual mutations. All-Atom Structure-Based Model. To study the effect of charge−charge interactions on the biophysical characteristics of proteins, we used an all-atom structure based model approach. In this model, the force field has a bias toward the native topology. The potential for the all-heavy atoms (nonhydrogen) in our model is defined as follows: V = VAA + VDH

⎛ ε k T ⎞1/2 r B ⎟ λD = ⎜ 2 ⎝ 8πkelecNAe I ⎠

(5)

where kB = 8.314 kJ mol−1 K−1, T = 300 K, NA = 6.022 × 1023 mol−1, and I = 0.5∑iqi2ci is the ionic strength. Here, ci is the concentration of the solution in mol/L. By considering a dilute solution of ionic strength I = 100 mM, eq 5 gives λD = 0.974 nm. In the above equation B(λD) is a salt-dependent coefficient and given as B(λD) =

(2)

ear / λD 1 + ar /λD

(6)

where ar ≈ 1.4 Å is the radius of ion. In case of a dilute solution (λD ≫ ar), B(λD) ≈ 1 and therefore can be neglected in the DH equation. Simulation Details. All simulations in the present paper for all the WT and MUT TNfn3 systems were performed using the Gromacs 4.0.740 without modifying the source code. We used Langevin dynamics with a time step τ = 0.004 ps and friction coefficient of 1 ps−1. To explore the folding thermodynamics, we used the Replica Exchange Molecular Dynamics (REMD) method,41 as implemented in Gromacs,40 for all the systems studied in this paper. We chose 20 replicas close to their corresponding folding temperature, TF, for the WT, WT, and MUT systems. In the reduced unit the temperature range 0.922 ≤ T* ≤ 0.96 with a uniform temperature step ΔT* = 0.002 was used for the consecutive replicas for all the systems, where T* ≡ kBT/εC is in the reduced temperature units. The range of temperatures of REMD simulations for all the systems was decided based on the observation that the proteins remain always folded at the lower range and unfolded at the higher end. For each replica we ran 3 × 108 steps of simulation, which correspond to a total time of 28.8 μs in Gromacs time scale. The data from the last 2.5 × 108 steps of each replica were used for all the analysis. During the simulation replica exchange was

where the first term defines the bonded interactions that is obtained using the united-atom force field ffg43a1 (GROMOS96). A more detailed description of Vbonded can be found in our previous work.24 The second term in eq 2 represents the nonbonded interactions in the protein given as Vnonbonded

−2

where kelec = 1/4πε0 = 138.935 kJ mol nm e , qi (qj) is the sign of the charge for charged atom i (j) (i.e., either +1, −1, or 0), rij is the distance between charged atoms i and j in nm, εr is the dielectric constant which is set to 80, and λD is the screening length (or Debye length) in nm. The positively charged residues Lys and Arg were given charge +1, and the negatively charged residues Glu and Asp were given −1 charge. The charge on the residues was placed at their corresponding central atom of the side-chains charged group as follows: Lys, Nζ; Arg, Cζ; Glu, Cδ; and Asp, Cγ. No charges were incorporated at the N- and C-termini of the protein. Linearization of the Poisson−Boltzmann equation yields the following relations for the Debye length

The first term in the above equation defines the all (heavy)atom potential similar to the approach by Clementi et al.39 that enforces a funneled energy landscape given by

⎡⎛ ⎞12 ⎛ σij ⎞6 ⎤ σij ⎢ ⎜ ⎟ = ∑ εC⎢⎜ ⎟ − 2⎜⎜ ⎟⎟ ⎥⎥ + r ⎝ rij ⎠ ⎦ contacts ⎣⎝ ij ⎠

(4) −1

(1)

VAA = Vbonded + Vnonbonded

εr rij

⎛ σ ⎞12 ∑ εNC⎜⎜ NC ⎟⎟ ⎝ rij ⎠ noncontacts

(3) −1

In eq 3 εC = 0.5 kcal mol is used as the strength of a native contact defined by a cutoff distance 4 Å between the atoms i and j separated by the corresponding amino acid residues |i′ − j′| > 3. rij is the distance between the atoms i and j in a particular conformation, and σij is the distance between the same atoms in the native state structure. Interactions between the atoms which are not in the native-contact are given a pure repulsive potential of strength, σNC = 0.01 kcal mol−1 and σNC = 3.5 Å. It should be noted that the value of εC that we use in the present study is 0.5 times lower than the value used in our previous folding study of TNfn324 (without the effects of electrostatics). The reason behind lowering the value of εC is to H

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

■

ACKNOWLEDGMENTS This work was supported by grants from the National Science Foundation (MCB-0110396, MCB-1330249 to G.I.M and MCB-1050966 to A.E.G). S.T. is partly supported by a training fellowship from the Keck Center of the Gulf Coast Consortia, on the Computational Cancer Biology Training Program, Cancer Prevention and Research Institute of Texas (CPRIT) RP140113.

attempted at every 1000 steps, and data were also written at every 1000 steps. Analyses. To analyze the degree of conservation for each residue of WT TNfn3, we followed the procedure by Mihalek et al.23 using the Evolutionary Trace (ET) server (http:// mammoth.bcm.tmc.edu/ETserver.html). ET server uses hybrid methods that combine evolutionary and entropic information from multiple sequence alignments (MSA). In this method, first the protein sequences homologous to TNfn3 were obtained with Basic Local Alignment Search Tool (BLAST) against the National Center for Biotechnology Information (NCBI) Entrez nonredundant protein sequence database, with the Expect (E) value smaller than 0.05. Second, the sequences identical to human TNfn3, or the short sequences with a length less than 0.8 times the length of human TNfn3, were removed from the sequence set. Third, the MSA was performed, and the results were subject to a phylogenetic tree construction. The basic assumption of the current method is that the more important the residue, the sooner it becomes fixed in the evolutionary branches. Finally, the real-value Evolutionary Trace (rvET) scores for each residue in TNfn3 was calculated based on the phylogenetic tree and were used to indicate the degree of conservation of the residue. In the present study we used Q as the global reaction coordinate to characterize all the systems of TNfn3. Q is the fraction of natively interacting residues that are in contact. We considered two residues in contact if any two atoms from those residues are within the 1.5 times their native distance σij, as defined in the eq 3. All thermodynamic calculations in this work were obtained from the weighted histogram analysis method (WHAM).42,43 The correlation between residue fluctuations is calculated as the dynamic cross-correlation map (DCCM), which is defined by the coefficient as cij =

■

(7)

where Δri and Δrj are the fluctuations of the position vectors of the corresponding C atoms of residues i and j, respectively. The position vectors are obtained by aligning each conformation to the native structure. The ⟨...⟩ in the above equation indicates ensemble average over the different states as following: unfolded (U, Q ∼ 0.0−0.3), intermediate (I, Q ∼ 0.3−0.6), and native states (N, Q ∼ 0.6−1.0). The coefficients cij in this study are calculated using the Bio3D package44 in R (https:// www.R-project.org/).

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b08527. Evolutionary trace (ET) analysis of WT TNfn3 and probability of native contact formation of the three TNfn3 systems (PDF)

■

REFERENCES

(1) Fu, K.; Klibanov, A. M.; Langer, R. Protein stability in controlledrelease systems. Nat. Biotechnol. 2000, 18, 24−25. (2) Manning, M. C.; Chou, D. K.; Murphy, B. M.; Payne, R. W.; Katayama, D. S. Stability of protein pharmaceuticals: an update. Pharm. Res. 2010, 27, 544−575. (3) Strickler, S. S.; Gribenko, A. V.; Gribenko, A. V.; Keiffer, T. R.; Tomlinson, J.; Reihle, T.; Loladze, V. V.; Makhatadze, G. I. Protein stability and surface electrostatics: a charged relationship. Biochemistry 2006, 45, 2761−2766. (4) Gribenko, A. V.; Patel, M. M.; Liu, J.; McCallum, S. A.; Wang, C.; Makhatadze, G. I. Rational stabilization of enzymes by computational redesign of surface charge-charge interactions. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 2601−2606. (5) Grimsley, G. R.; Shaw, K. L.; Fee, L. R.; Alston, R. W.; HuyghuesDespointes, B. M.; Thurlkill, R. L.; Scholtz, J. M.; Pace, C. N. Increasing protein stability by altering long-range coulombic interactions. Protein Sci. 1999, 8, 1843−1849. (6) Marshall, S. A.; Morgan, C. S.; Mayo, S. L. Electrostatics significantly affect the stability of designed homeodomain variants. J. Mol. Biol. 2002, 316, 189−199. (7) Gribenko, A. V.; Makhatadze, G. I. Role of the charge-charge interactions in defining stability and halophilicity of the CspB proteins. J. Mol. Biol. 2007, 366, 842−856. (8) Perl, D.; Mueller, U.; Heinemann, U.; Schmid, F. X. Two exposed amino acid residues confer thermostability on a cold shock protein. Nat. Struct. Biol. 2000, 7, 380−383. (9) Pace, C. N.; Alston, R. W.; Shaw, K. L. Charge-charge interactions influence the denatured state ensemble and contribute to protein stability. Protein Sci. 2000, 9, 1395−1398. (10) Cho, J. H.; Sato, S.; Horng, J. C.; Anil, B.; Raleigh, D. P. Electrostatic interactions in the denatured state ensemble: their effect upon protein folding and protein stability. Arch. Biochem. Biophys. 2008, 469, 20−28. (11) Zhou, H. X. Residual electrostatic effects in the unfolded state of the N-terminal domain of L9 can be attributed to nonspecific nonlocal charge-charge interactions. Biochemistry 2002, 41, 6533−6538. (12) Weinkam, P.; Pletneva, E. V.; Gray, H. B.; Winkler, J. R.; Wolynes, P. G. Electrostatic effects on funneled landscapes and structural diversity in denatured protein ensembles. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 1796−1801. (13) Sali, D.; Bycroft, M.; Fersht, A. R. Surface electrostatic interactions contribute little of stability of barnase. J. Mol. Biol. 1991, 220, 779−788. (14) Onuchic, J. N.; Luthey-Schulten, Z.; Wolynes, P. G. Theory of protein folding: the energy landscape perspective. Annu. Rev. Phys. Chem. 1997, 48, 545−600. (15) Bryngelson, J. D.; Onuchic, J. N.; Socci, N. D.; Wolynes, P. G. Funnels, pathways, and the energy landscape of protein folding: a synthesis. Proteins: Struct., Funct., Genet. 1995, 21, 167−195. (16) Dill, K. A.; Chan, H. S. From Levinthal to pathways to funnels. Nat. Struct. Biol. 1997, 4, 10−19. (17) Clementi, C.; Nymeyer, H.; Onuchic, J. N. Topological and energetic factors: what determines the structural details of the transition state ensemble and “en-route” intermediates for protein folding? An investigation for small globular proteins. J. Mol. Biol. 2000, 298, 937−953. (18) Cheung, M. S.; Garcia, A. E.; Onuchic, J. N. Protein folding mediated by solvation: water expulsion and formation of the

⟨ΔriΔrj⟩ ⟨Δri 2⟩1/2 ⟨Δrj 2⟩1/2

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected]; Ph 518-276-4417; Fax 518-276-2955 (G.I.M.). Notes

The authors declare no competing financial interest. I

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B hydrophobic core occur after the structural collapse. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 685−690. (19) Zhou, R.; Huang, X.; Margulis, C. J.; Berne, B. J. Hydrophobic collapse in multidomain protein folding. Science 2004, 305, 1605− 1609. (20) Dill, K. A.; Fiebig, K. M.; Chan, H. S. Cooperativity in proteinfolding kinetics. Proc. Natl. Acad. Sci. U. S. A. 1993, 90, 1942−1946. (21) Sanchez-Ruiz, J. M.; Makhatadze, G. I. To charge or not to charge? Trends Biotechnol. 2001, 19, 132−135. (22) Leahy, D. J.; Hendrickson, W. A.; Aukhil, I.; Erickson, H. P. Structure of a fibronectin type III domain from tenascin phased by MAD analysis of the selenomethionyl protein. Science 1992, 258, 987− 991. (23) Mihalek, I.; Res, I.; Lichtarge, O. A family of evolution-entropy hybrid methods for ranking protein residues by importance. J. Mol. Biol. 2004, 336, 1265−1282. (24) Tripathi, S.; Makhatadze, G. I.; Garcia, A. E. Backtracking due to residual structure in the unfolded state changes the folding of the third fibronectin type III domain from tenascin-C. J. Phys. Chem. B 2013, 117, 800−810. (25) Clarke, J.; Hamill, S. J.; Johnson, C. M. Folding and stability of a fibronectin type III domain of human tenascin. J. Mol. Biol. 1997, 270, 771−778. (26) Hamill, S. J.; Meekhof, A. E.; Clarke, J. The effect of boundary selection on the stability and folding of the third fibronectin type III domain from human tenascin. Biochemistry 1998, 37, 8071−8079. (27) Ibarra-Molero, B.; Sanchez-Ruiz, J. M. Genetic algorithm to design stabilizing surface-charge distributions in proteins. J. Phys. Chem. B 2002, 106, 6609−6613. (28) Tsai, M. Y.; Zheng, W.; Balamurugan, D.; Schafer, N. P.; Kim, B. L.; Cheung, M. S.; Wolynes, P. G. Electrostatics, structure prediction, and the energy landscapes for protein folding and binding. Protein Sci. 2015, DOI: 10.1002/pro.2751. (29) Azia, A.; Levy, Y. Nonnative electrostatic interactions can modulate protein folding: molecular dynamics with a grain of salt. J. Mol. Biol. 2009, 393, 527−542. (30) Tzul, F. O.; Schweiker, K. L.; Makhatadze, G. I. Modulation of folding energy landscape by charge-charge interactions: linking experiments with computational modeling. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, E259−266. (31) Naganathan, A. N.; Sanchez-Ruiz, J. M.; Munshi, S.; Suresh, S. Are protein folding intermediates the evolutionary consequence of functional constraints? J. Phys. Chem. B 2015, 119, 1323−1333. (32) Zarrine-Afsar, A.; Zhang, Z.; Schweiker, K. L.; Makhatadze, G. I.; Davidson, A. R.; Chan, H. S. Kinetic consequences of native state optimization of surface-exposed electrostatic interactions in the Fyn SH3 domain. Proteins: Struct., Funct., Genet. 2012, 80, 858−870. (33) Reddy, G.; Thirumalai, D. Dissecting Ubiquitin Folding Using the Self-Organized Polymer Model. J. Phys. Chem. B 2015, 119, 11358−11370. (34) Nymeyer, H.; Garcia, A. E.; Onuchic, J. N. Folding funnels and frustration in off-lattice minimalist protein landscapes. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 5921−5928. (35) Shea, J. E.; Onuchic, J. N.; Brooks, C. L., III Exploring the origins of topological frustration: design of a minimally frustrated model of fragment B of protein A. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 12512−12517. (36) Ferreiro, D. U.; Hegler, J. A.; Komives, E. A.; Wolynes, P. G. Localizing frustration in native proteins and protein assemblies. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 19819−19824. (37) Shea, J. E.; Onuchic, J. N.; Brooks, C. L. Energetic frustration and the nature of the transition state in protein folding. J. Chem. Phys. 2000, 113, 7663−7671. (38) Sali, A.; Blundell, T. L. Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 1993, 234, 779−815. (39) Clementi, C.; Garcia, A. E.; Onuchic, J. N. Interplay among tertiary contacts, secondary structure formation and side-chain packing in the protein folding mechanism: all-atom representation study of protein L. J. Mol. Biol. 2003, 326, 933−954.

(40) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (41) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314, 141−151. (42) Ferrenberg, A. M.; Swendsen, R. H. New Monte Carlo technique for studying phase transitions. Phys. Rev. Lett. 1988, 61, 2635−2638. (43) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte Carlo data analysis. Phys. Rev. Lett. 1989, 63, 1195−1198. (44) Grant, B. J.; Rodrigues, A. P.; ElSawy, K. M.; McCammon, J. A.; Caves, L. S. Bio3d: an R package for the comparative analysis of protein structures. Bioinformatics 2006, 22, 2695−2696. (45) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38.

J

DOI: 10.1021/acs.jpcb.5b08527 J. Phys. Chem. B XXXX, XXX, XXX−XXX