Probing α/β Balances in Modified Amber Force Fields from a Molecular Dynamics Study on a ββα Model Protein (1FSD)

Abstract

1FSD is a 28-residue designed protein with a ${\beta}{\beta}{\alpha}$ motif. Since this protein displays most essential features of protein structures in such a small size, this model protein can be an outstanding system for evaluating the balance in the propensity of the secondary structures and the quality of all-atom protein force fields. Particularly, this protein would be difficult to fold to its correct native structure without establishing proper balances between the secondary structure elements in all-atom energy functions. In this work, a series of the recently optimized five amber protein force fields [$ff03^*$, $f99sb^*$-ildn, ff99sb-${\phi}^{\prime}$-ildn, ff99sb-nmr1-ildn, ff99sb-${\Phi}{\Psi}$(G24, CS)-ildn] were investigated for the simulations of 1FSD using a conventional molecular dynamics (MD) and a biased-exchange meta-dynamics (BEMD) methods. Among those tested force fields, we found that ff99sb-nmr1-ildn and ff99sb-${\Phi}{\Psi}$(G24, CS)-ildn are promising in that both force fields can locate the native state of 1FSD with a high accuracy (backbone rmsd ${\leq}1.7{\AA}$) in the global free energy minimum basin with a reasonable energetics conforming to a previous circular dichroism (CD) experiment. Furthermore, both force fields led to a common set of two distinct folding pathways with a heterogeneous nature of the transition state to the folding. We anticipate that these force fields are reasonably well balanced, thereby transferable to many other protein folds.

Introduction

Molecular dynamics (MD) simulation has been a powerful computational tool for investigating protein folding charac-teristics. Recently, notable improvements in computer hard-ware and novel computational algorithms for MD simulation methods have enabled fully atomistic MD simulations to reach the time scale of folding events of small proteins. Since empirical force fields directly influence the outcome of MD simulations, accuracy of force fields is critical for unbiased interpretation of in silico protein folding analysis. For this reason, systematic efforts to improve pre-existing all-atom force fields (charmm or amber) have been made by optimizing protein energy parameters of backbone or side chains using various NMR experimental data. One of the immediate concerns in such improved all-atom protein force fields (energy function) is whether the energy functions are balanced to properly define the propensity of secondary structures in mixed folds. Although these new force fields have been refined in many ways compared to their previous counterparts, it remains to be seen whether any bias towards α-helix and β-sheet structures can be alleviated through additional optimization procedures. Unfortunately, in the development of all-atom force fields, balancing the propen-sity of the secondary structures is not a trivial task, because even a small imbalance or error in protein backbone para-meters can facilitate a substantial bias towards α or β struc-tures, thereby producing unrealistic folding free energy land-scapes of α/β mixtures. Therefore, it is particularly important to probe α/β balances in newly optimized protein force fields. For evaluating α/β balances in such new protein force fields, we chose 1FSD1 as a benchmark model system. 1FSD is a 28-residue ββα model protein designed for the zinc-finger domain. This protein consists of the most essential features of protein structure: N-terminal β-hairpin and C-terminal α-helix along with a well-defined hydrophobic core. Thermodynamically, this protein is marginally stable at ambient temperatures, as shown by an earlier circular di-chroism (CD) experiment.1 Particularly, its β-strand is relatively weak due to the lack of backbone hydrogen bonds.2 Despite its computationally tractable size, it is still difficult to fold 1FSD to its correct native structure even with much advanced simulation methods unless empirical potential energy func-tions retain the subtle α/β balance required for 1FSD.

For 1FSD or its double-mutant 1FME3 (Q1E and I7Y), several all-atom simulations with implicit or explicit water solvents using advanced MD simulation methods have been reported.2,4-10 Previous all-atom simulations of 1FSD with a generalized Born (GB) implicit solvation successfully predicted the native structure of 1FSD.5-9 In some of the im-plicit water simulation studies, the backbone energy term in the amber ff99 force field were re-parameterized to improve the propensity of the secondary structures.5,6 Recently, it was found that even without re-optimization, a combination of standard amber ff96 with GB implicit solvation models produce a reliable free energy landscape for 1FSD folding.7 In earlier simulations with explicit water solvents, Wang et al. first applied amber ff0311/TIP3P12 to a replica exchange molecular dynamics (REMD) simulation of 1FSD.4 Unfor-tunately, the REMD simulation time was too short (18 ns) to justify the validity of this force field in the case of 1FSD; this led us to run more test simulations using ff03/TIP3P. Our extensive MD simulation of 1FSD with ff03/TIP3P at 300 K indicated that the β-hairpin segment of 1FSD was prone to unfold and eventually failed to preserve the correct native structure of 1FSD (data not shown). Thus, ff03/TIP3P may not be a suitable choice for the simulation of 1FSD. One of the notable all-atom simulation with explicit solvation is a direct folding study by Shaw et al.13 using a modified charmm force field14 (charmm22*/TIP3P).

Using their specialized high performance computer,15 they applied conventional MD simulations at elevated temper-atures to 12 proteins including 1FME and successfully observed reversible folding events for all those proteins. The results of their simulation of 1FME indicated that its native state was captured within a Cα rmsd value of 1.6 Å at 325 K.13 Despite such remarkable progresses in ab-initio protein folding MD simulations using the charmm22* force field, the melting point of 1FME in this force field was about 0°C, indicating that charmm22* may not yield a thermodynamically stable native state for 1FME at ambient temperatures. Since the double mutation of 1FME would induce more compact hydrophobic core packing than 1FSD, the simulated melting point of 1FSD is expected to be even lower than that of 1FME, which is in contrast with a previously reported experi-mental melting temperature of 42°C. This discrepancy may open up the possibility for more refinements of the charmm22* force field.

Recently, in an effort to increase the accuracy of all-atom force fields, several variants of highly optimized amber force fields have been developed by fitting to quantum mechanical data,16,17 overall helical propensities,18-22 or various NMR data23-25 of several trial peptide and proteins. The main goal of the present work is to provide some insights into this new generation of amber force fields by assessing how well these new energy functions can maintain the thermodynamic stability of α/β mixed folds at ambient temperatures. Here, we report the results of test simulation of 1FSD using a series of recently modified amber force fields in conjunction with explicit water solvation. The following five amber force fields were employed: ff03*,18 ff99sb*-ildn,17,18 ff99sb-φ'- ildn,23 ff99sb-nmr1-ildn,24 and ff99sb-ΦΨ(G24,CS)-ildn.25 The parent force field for these force fields is either ff0311 or ff99sb.16 As a first step to test these force fields for 1FSD, starting from the native structure of 1FSD, we performed a total of 1.0 μs MD simulations at 300 K using five modified amber force fields with the TIP3P solvation model. Among the five force fields that were tested, we found that the two energy functions of ff99sb-nmr1-ildn (nmr1) and ff99sb-ΦΨ(G24, CS)-ildn (cmap) were able to steadily maintain the native NMR structure of 1FSD. The 1 μs MD simulation, however, may not be sufficient to decisively address the stability of the native structure under the thermodynamic principle. Subsequently, these two promising force fields (nmr1 and cmap) were used for directly computing the fold-ing free energy landscape of 1FSD at 300 K. For well-con-verged and effective free energy mapping at that temperature, we employed a biased exchange meta-dynamics (BEMD) simulation method.26 The BEMD is a replica exchange method27 in combination with a meta-dynamics scheme28 and has been proven to be a viable tool for effectively ex-ploring complex conformational spaces of biomolecules. In this work, we found that the two force fields (nmr1 and cmap) fulfill the thermodynamic requirement for the native state by predicting the native structure of 1FSD in the lowest free energy minimum region with the native state stability conforming to an earlier circular-dichoism (CD) experiment.1,2 Similarities and differences in resulting free energy maps and folding pathways were compared in the Result and Discussion section.

 

Methods

The NMR structure of 1FSD was taken from the protein database (PDB entry: 1FSD). Its sequence is QQYTA KIKGR TFRNE KELRD FIEKF KGR. This protein was solvated by 3614 TIP3P water molecules in a 54.3 Å dodecahedron box. Five Cl− ions were added for charge neutrality of the simu-lated system. Then, the system including the protein and water molecules was minimized by the steepest descent method. All the chemical bonds were fixed with the LINCS algorithm.29 The side chain groups of –NH2, –CH3, and –NH3+ were treated as virtual interaction sites.30 With both LINCS and virtual site scheme included, a time step of 4 fs was applied for all the simulations. The minimized system was subject to an equilibration for 200 ps using an isothermal-isobaric MD method at 300 K and 1 atm. For this equilib-ration run, the simulation temperature was maintained at 300 K using the modified Berendsen method of Bussi et al.31 with a coupling constant of 0.1 ps. The system pressure was also fixed with the Berendsen barostat.32 The Lennard-Jones potential with a dispersion correction was computed using a cutoff distance of 10 Å. The electrostatic potential energy was calculated with the particle-mesh Ewald method33 with a cutoff distance of 10 Å. The equilibration run was follow-ed by a constant temperature MD simulation for 1.0 μs using each of the aforementioned five force fields (ff03*, ff99sb*-ildn, ff99sb-φ'-ildn, ff99sb-nmr1-ildn, and ff99sb-ΦΨ(G24, CS)-ildn) in conjunction with the TIP3P water solvation. The ff03* force field18 is a modified version of the amber ff03 force field. In the ff03* force field, the ψ backbone di-hedral potential was corrected. Although this energy func-tion was optimized primarily for helix-coil transitions of polypeptides, the MD simulation result with ff03*/TIP3P appears to be quite promising as the correct folding of both Villin HP35 (three stranded α-helix) and the pin WW domain (three stranded β-sheet) was observed.34 Though the ff99sb* force field18 was optimized in the same way as ff03*, its correction was derived from the ff99sb force field as reported by Simmerling et al.16 The performance of both ff03* and ff99sb* is virtually comparable.35 The ff99sb*-ildn force field is an extended version with side chain corre-ction to the amber ff99sb*. Here, the suffix “-ildn” in the force field name denotes an inclusion of the newly optimiz-ed χ1 torsional potential of Shaw et al.17 The benefit of using the ildn side chain correction is a substantial improvement of PDB rotamer distributions.17 The ff99sb-φ'-ildn energy func-tion23 is augmented with the φ'-backbone term (C-N-Cα-Cβ atoms) in ff99sb. This feature of the φ'-backbone correction resulted in more accurate intrinsic conformational prefer-ences. The ff99sb-nmr1-ildn24 was obtained in such a way as to reproduce NMR chemical shifts of all carbon nuclei Cα, Cβ, and C'. Similarly, the ff99sb-ΦΨ(G24,CS)-ildn25 was optimized against NMR chemical shift data with the back-bone Φ/Ψ cross terms (CMAP scheme) added to the ff99sb-ildn force field. Both ff99sb-nmr1 and ff99sb-ΦΨ(G24,CS)-ildn force fields displayed a substantial improvement over their parent force field (ff99sb).

For the free energy calculation, the BEMD method employ-ing the well-tempered scheme36 was used as an enhanced sampling strategy. This method uses several biased replicas evolving under a time-dependent potential that is construct-ed as the sum of Gaussian type functions deposited along the trajectory in pre-defined collective variable (CV) space. For each biased replica, the Gaussian potential is added in real time and the simulation at a constant temperature was em-ployed to explore the conformational space of each CV. For the deposition of the Gaussian potential, an initial hill size of w = 0.2 kJ/mol was used with a hill deposition time of τG = 2 ps. The canonical distribution of the unbiased replica was reproduced in the limit of w/τG = 0. In practice, however, a finite value of w/τG = 0.1 kJ/mol/ps was employed. For the biased replicas, the well-tempered scheme of Parrinello et al.36 was used in conjunction with a bias factor of (T + ΔT)/T = 10, where T and ΔT are a system and a biased temperature, respectively. As a result of the well-tempered scheme, initially w/τG was a finite value of 0.1 kJ/mol/ps, but this factor gradually decreased during the simulation, since the w value decayed in time depending on the bias factor. Here, a total of four replicas were used in the BEMD simulation at T = 300 K. One of those replicas was a regular (unbiased) replica for a normal MD with an original potential energy function and other three were the biased replicas evolving in time with the meta-dynamics using biased potential energy functions in CV spaces.

In this work, a total of three CVs were used with a choice of their Gaussian widths of σCV1 = 0.1, σCV2 = 0.1, and σCV3 = 0.5 Å. Among them, two CVs represent the amounts of α-helix blocks (CV1) and anti-parallel β-blocks (CV2) used to measure α- and β-secondary structure contents of 1FSD.37 The last CV is the radius of gyration (CV3) of the hydro-phobic core residues consisting of Ala5, Ile7, Phe12, Leu18, Phe21, Ile22, and Phe25. Therefore, the aforementioned three CVs were aimed for efficient sampling of the second-ary structures and tertiary hydrophobic contacts pertaining to 1FSD. The exchange event between any random pair of neighboring replicas was attempted every 20 ps. All BEMD simulations were initiated from the native NMR structure. We performed a total of 2.0 μs BEMD simulation for each replica using ff99sb-nmr1-ildn (nmr1) and ff99sb-ΦΨ(G24, CS)-ildn (cmap). The simulations were performed using the Gromacs 4.5.3 program38 in conjunction with the Plumed 1.3 program.39

The free energy surfaces for these force fields were com-puted using the simulated trajectory from the unbiased replica. For the free energy representation, the root mean square deviations of the backbone (RMSD_bb), the N-terminal β-blocks from residues 3-13 (RMSD_β), and the C-terminal α-helix blocks from residues 14-26 (RMSD_α) were used as reaction coordinates. For these coordinates, the NMR native structure1 was taken as a reference. In addition, the radius of gyration of the hydrophobic core residues (Rg_core) was employed.

 

Results and Discussion

Evaluation of Force Fields with 1-μs MD Simulation at Ambient Temperature. As a preliminary test for evaluating the native state stability and α/β balances of 1FSD, we carried out 1 μs MD simulation at 300 K for each of the aforementioned five amber force fields. Figure 1 shows the time profile of the backbone root mean square deviation (RMSD_bb) value for each force field. Furthermore, the averaged RMSD_bb values for all the tested energy func-tions are given in Figure 2. As shown in Figures 1 and 2, the two force fields of nmr1 and cmap had the capacity for preserving the native NMR structure of 1FSD with minimal RMSD fluctuations during the entire 1 μs MD run. The energy functions of other facilitated immediate (ff03* and ff99sb*-ildn) or gradual (ff99sb-φ'-ildn) unfolding of the native NMR geometry via disruptions of the secondary structure elements. In terms of preserving the native state, the ff99sb-φ'-ildn force field seems better than ff03* and ff99sb*-ildn. The ff99sb-φ'-ildn maintained the β-hairpin segment, but it had the tendency to unfold the α-helical region of 1FSD. Such underestimated helical propensity is in agreement with a previous result from an extensive force field evaluation.35 More detailed information on the secondary structural elements can be seen in Figure 3, where the propensity of secondary structures was traced in time for each force field using the DSSP program.40 This result clearly indicated that when nmr1 or cmap force field is applied, the secondary structures are properly maintained. Compared to the α-helix part, the β-hairpin segment, however, undergoes somewhat enhanced structural fluctuations.

Figure 1.The root mean square deviations of the protein backbone (RMSD_bb) in time profiles for each force field; (a) ff03*, (b) ff99sb*-ildn, (c) ff99sb-nmr1-ildn, (d) ff99sb-φ'-ildn, and (e) ff99sb-ΦΨ(G24,CS)-ildn. This result was obtained from 1 μs MD simulation at 300 K, starting from the native NMR structure of 1FSD.

Figure 2Averaged backbone RMSD values for the five force fields.

Figure 3.Time profiles of the secondary structure propensities for the five force fields. The structure in the right side is the final representative at 1.0 μs. The DSSP program40 was used for the assignment. The ‘N’ letter indicates the N-terminal of protein.

Thermodynamic Stability of the Native State. From the results of our MD analysis, it is evident that the use of either nmr1 or cmap force field may be a better choice for in silico folding simulations targeting 1FSD. For quantitative inter-pretations of the thermodynamic stability of this protein, free energy computations were performed using the nmr1 and cmap force fields. For the purpose of generating folding free energy landscapes, we utilized a total of 2.0 μs trajectory from the unbiased replica in the BEMD. The convergence of the present simulation was monitored with the block averag-ed fold population (Pfold) in each time block from the unbiased replica (Figure 4). As shown in Figure 4, the time profile of the Pfold value indicates unfolding/refolding events and thereby a reasonably fair convergence appears to be reached for the nmr1 and the cmap force fields, respectively. Accordingly, the first 200 ns (nmr1) and 400 ns (cmap) trajectories were discarded in the corresponding free energy mapping. Moreover, after 400 ns, the three biased replicas, which run under a well-tempered scheme, appeared to reach slowly varying states in time, as seen from the time course of average hill size of the depositing potential for each CV space (Figure S1). Having demonstrated the convergence of our simulations, the unbiased trajectory was projected into a proper set of reaction coordinates for the free energy repre-sentation

At first, the RMSD_bb and Rg_core values were used to construct the two-dimensional (2D) free energy surface. In comparison, two free energy maps resulting from both nmr1 and cmap force fields are shown in Figure 5. It is encourag-ing that both force fields at 300 K can locate the native state of 1FSD in the lowest free energy minimum basin centered at (RMSD_bb, Rg_core) = (1.7 Å, 5 Å). From the earlier CD experimental result, using a two-state folding model, one can obtain a free energy difference between the native and unfold basins ( ) at 300 K. Thus, the quality of our free energy surfaces can be evaluated in comparison with the experimental value. Consistent with the experimental value ( = 0.5 kcal/mol at 300 K), the predicted values for nmr1 and cmap were 0.8 kcal/mol and 0.6 kcal/mol, respectively.

Figure 4.200 ns block averaged folded fractions from the trajectory of the unbiased replica for (a) ff99sb-nmr1-ildn and (b) ff99sb-ΦΨ(G24, CS)-ildn. For the definition of the fold state, RMSD_bb < 3.0 Å was used.

Figure 5.Free energy surface with the backbone RMSD (RMSD_bb) and the radius of hydration of the hydrophobic core residues (Rg_core) for (a) ff99sb-nmr1-ildn (b) ff99sb-ΦΨ(G24, CS)-ildn. The hydrophobic core residues consist of Ala5, Ile7, Phe12, Leu18, Phe21, Ile22, and Phe25. The contour spacing is 0.5 kcal/mol.

Comparisons of Free Energy Landscapes and Folding Pathways. As illustrated in Figure 5, the nmr1 and the cmap predicted the native state in the identical basin of the free energy surface. A notable difference between the two free energy maps is the overall topography of the unfold state. Thus, based on the progress of hydrophobic core formation, the nmr1 and the cmap display somewhat different folding pathways. In the case of the nmr1, the native structure was formed gradually with the aggregation of the hydrophobic core residues. On the other hand, in the case of cmap, the collapse of the hydrophobic core initiated the folding of 1FSD. Although the free energy map with (RMSD_bb, Rg_core) gives an overall picture of folding energy landscape of 1FSD, detailed information regarding the secondary structure for-mation is missing. In order to provide a clear view of the folding order of the α- and β-formation, another set of reac-tion coordinates (RMSD_α and RMSD_β) were introduced for subsequent free energy mapping. For such representa-tions in combination with RMSD_bb, the corresponding free energy surfaces with (RMSD_bb, RMSD_α) and (RMSD_bb, RMSD_β) are shown in Figure 6. The free energy map with (RMSD_bb, RMSD_β) shows a gradual formation of the β-hairpin part. Interestingly, the free energy representation with (RMSD_bb, RMSD_α) allowed at least two different folding routes (Paths I and II) to be well resolved. As a consequence, two different transition states (TS), T1 and T2, were associated with the two folding pathways (see Figure 7). Interestingly, both nmr1 and cmap force fields predicted Path I to be the major folding pathway. The observed dual folding pathways are illustrated for each force field in Figure 7. In the first folding pathway (Path I), α-helix is formed first followed by the concurrent formation of β-hairpin and tertiary structure (U1→T1→N). Interestingly, this pathway is consistent with a previously suggested pathway from all-atom simulations with a GB implicit solvation model.5-7 The other folding pathway (Path II), which is rather new, sug-gested the simultaneous formation of the secondary struc-tural elements (α and β segments) and the hydrophobic tertiary core (U2→T2→N).

Figure 6.Free energy surface resulting from (a) ff99sb-nmr1-ildn and (b) ff99sb-ΦΨ(G24, CS)-ildn with the backbone RMSD and the α-helix RMSD (residues 14-26) (top) and with the backbone RMSD and the β-hairpin RMSD (residues 3-13) (bottom).

Duan and his colleagues5 also proposed dual folding path-ways of 1FSD using a balanced amber force field with a GB implicit solvation model. In their study, consistent with Path I described in this work, the major folding pathway facilitat-ed the initial formation of α-helix followed by β-hairpin formation and the completion of hydrophobic core packing. In addition, they claimed that in the minor folding pathway β-hairpin formation initiated the folding of 1FSD. In our simulation, however, the minor pathway (Path II) shows that both α-helix and β-strand propagate simultaneously near the pre-fold loop region intervening α-helix and β-hairpin regions.

Figure 7.Folding pathways resulting from (a) ff99sb-nmr1-ildn and (b) ff99sb-ΦΨ(G24, CS)-ildn.

 

Conclusions

1FSD is a 28-residue designed protein with a ββα motif, which is considered as a highly valuable model system representing simplest possible α/β mixed folds. In general, the native state of 1FSD is managed by a subtle balance between α-helix and β-strand. Therefore, even a small bias in the propensity of the secondary structures can lead to an easy disruption of the native structure at ambient temper-ature. For this reason, 1FSD can be an excellent model system for assessing the bias toward α- or β-strands in all-atom energy functions. In order to provide a practical insight into force field selections in targeting mixed folds, we applied a recently optimized set of five amber force fields to the simulation of 1FSD. We found that the nmr1 and cmap force fields were outstanding in preserving the correct native structure of 1FSD. As a result, these two force fields were further applied to free energy computations using an enhanc-ed sampling method. Our free energy simulations indicated that both nmr1 and cmap force fields guarantee that the native state of 1FSD unequivocally matches the lowest free energy minimum state at 300 K. Moreover, the free energy differences between native and unfold states for both force fields were in accord with the CD experimental value, indicating that a reasonable quality of the folding energetics could also be obtained. Further investigations of the present simulation data led to the identification of dual pathways for the folding of 1FSD. Although these predicted pathways need to be confirmed by experiments, these pathways can provide a useful guideline for understanding multiple fold-ing behaviors of 1FSD.

The nmr1 and cmap force fields were derived from the same parent force field (ff99sb). During the optimization of nmr1 and cmap, the protein backbone dihedral angle para-meters were adjusted, so that MD-generated trajectories reproduced NMR chemical shift data of the tested proteins. In the former, the backbone parameters were re-optimized and in the latter, the backbone parameters defined by the parent force field were retained, but the cross backbone terms similar to the CMAP scheme in charmm41 were includ-ed for the optimization. Since both force fields were NMR-based potentials obtained in a similar fashion, the accuracy of these force fields are expected be at least comparable.

In the previous force field evaluation,35 the nmr1 force field had the highest accuracy among many different force fields, but it still appeared to be slightly biased to helical propensity. Nevertheless, as shown in this work, the nmr1 shows a satisfactory performance in describing one of the difficult α/β model proteins. Although more critical assess-ments on these force fields are needed for future protein folding studies, due to their reasonable balance between α-helix and β-sheet, the nmr1 as well as the cmap are expected to be physically reliable choices for simulating much diverse protein folds.