Multifunnel Energy Landscapes for Phosphorylated Translation Repressor 4E-BP2 and Its Mutants.

Upon phosphorylation of specific sites, eukaryotic translation initiation factor 4E (eIF4E) binding protein 2 (4E-BP2) undergoes a fundamental structural transformation from a disordered state to a four-stranded β-sheet, leading to decreased binding affinity for its partner. This change reflects the significant effects of phosphate groups on the underlying energy landscapes of proteins. In this study, we combine high-temperature molecular dynamics simulations and discrete path sampling to construct energy landscapes for a doubly phosphorylated 4E-BP218-62 and two mutants (a single site mutant D33K and a double mutant Y54A/L59A). The potential and free energy landscapes for these three systems are multifunneled with the folded state and several alternative states lying close in energy, suggesting perhaps a multifunneled and multifunctional protein. Hydrogen bonds between phosphate groups and other residues not only stabilize these low-lying conformations to different extents but also play an important role in interstate transitions. From the energy landscape perspective, our results explain some interesting experimental observations, including the low stability of doubly phosphorylated 4E-BP2 and its moderate binding to eIF4E and the inability of phosphorylated Y54A/L59A to fold.

Introduction

Many cellular processes such as signaling pathways,[2−5] membrane transport,[6−8] and gene transcription[9−12] and translation[13−16] are orchestrated by highly dynamic phosphorylation/dephosphorylation of proteins. As the most common type of reversible post-translational modification (PTM) of proteins, phosphorylation can exert a profound influence on protein structures and functions.[17] Upon phosphorylation, proteins can undergo fundamental structural transformations, leading to altered binding affinities for their partners. Many major conformational changes upon phosphorylation have been observed in intrinsically disordered proteins (IDPs)[18,19] and regions (IDRs).[20−22]

One of the most dramatic structural changes induced by phosphorylation reported so far was observed in eukaryotic translation initiation factor 4E (eIF4E) binding protein 2 (4E-BP2).[18] Human 4E-BP2 is a 120-residue intrinsically disordered protein[23] with transient secondary structures.[24] Nonphosphorylated or minimally phosphorylated 4E-BP2s can compete with eukaryotic translation initiation factor 4G (eIF4G) for binding to eIF4E through a conserved YxxxxLΦ binding motif,[25,26] where x is any amino acid and Φ represents a hydrophobic amino acid. In complex with eIF4E, the binding motif adopts α helical conformations (Figure B).[27] Upon multiphosphorylation at five sites (T37, T46, S65, T70, and S83), 4E-BP2 folds into a four-stranded β-sheet (Figure A) and the structural change dissociates it from eIF4E.[18] Among the five sites, T37 and T46 are phosphorylated first[28] and the dual phosphorylation is sufficient to trigger the structural transformation, while phosphorylation at the other sites further stabilizes the folded state.

Figure 1

(A) Four-stranded β-sheet structure of phosphorylated 4E-BP218–62 (PDB ID: 2MX4),[18] with the two phosphorylated threonines (pT37 and pT46) shown as sticks. (B) unphosphorylated 4E-BP247–65 in complex with eIF4E (PDB ID: 3AM7).[27] 4E-BP247–62 is colored with the same color scheme as in (A), and the rest is in gray.

Several computational studies have been conducted to study the folding mechanism of this system and the origin of stabilization by phosphorylation. Gopi et al.[29] studied different phosphorylated states using the statistical mechanical Wako–Saitô–Muñoz–Eaton (WSME) model.[30,31] They found that phosphorylation-induced stabilization in this system originated from strong electrostatic interactions between phosphate groups and nearby Arg residues, as generally observed in other systems.[17] In the same study, a folding mechanism was also proposed starting with pairing of the third and fourth β-strands. Since the phosphorylated residues are located within the two turn motifs, the folding of turn motifs with/without phosphorylation was studied in atomic detail by Bomblies et al. and the folding free energy of each turn motif was estimated to be −2 kJ/mol upon phosphorylation.[32] Zeng et al.[33] identified several intermediates during unfolding via all-atom replica exchange molecular dynamics (REMD) simulations starting from the folded structure, but REMD simulations starting from unfolded conformations failed to find the experimentally determined folded state.

There are still some interesting properties of this system that are not well understood either experimentally or computationally. These include the following: (1) the doubly phosphorylated form (pT37pT46) is only weakly stable and has moderate binding affinity (weaker than that of the unphosphorylated form but stronger than that of fully phosphorylated form) for eIF4E.[18] The moderate binding affinity indicates the existence of binding-competent minor state(s), which is in agreement with NMR results.[18] However, detailed structural characteristics of these minor conformers and interconversions between different conformers cannot be inferred with existing experimental methods. (2) The stability of the fold is very sensitive to mutations. Two mutants with phosphorus-mimicking residues (Asp or Glu), T37D/T46D and T37E/T46E, cannot fold. Mutating glycine located in turn regions to valine (G39V/G48V) also prevents folding.[18] While the destabilizing effect of these mutations can be explained in terms of unfavorable electrostatic interactions or steric clashes within the turn regions, other mutants are hard to understand. For instance, phosphorylated Y54A/L59A [p(Y54A/L59A)], with two residues in the canonical binding motif (Y54xxxxL59Φ) substituted by alanine, cannot fold but retains the two β-turns. Another mutant, p(D33K), was predicted to have a more stable fold compared to the doubly phosphorylated wild-type (pWT) by Gopi et al.[29] Neither of the two mutant stabilities can be understood from pWT’s experimental structure.

To elucidate these properties from the perspective of energy landscapes, we combine high-temperature molecular dynamics (MD) simulations and techniques based on geometry optimization to construct the energy landscapes for pWT and the two interesting mutants. High-temperature MD simulations encourage energy barrier crossings and are popular in folding/unfolding studies of slowly folding proteins.[34] Here, we use high-temperature simulations to identify possible minor states. Several interesting minor conformers of pWT are extracted from high-temperature trajectories and further optimized via basin-hopping;[35] kinetically relevant paths connecting them to the folded state are then constructed via discrete path sampling (DPS).[36,37] Some of the minor states can explain not only the moderate binding affinity but also some nuclear Overhauser effect (NOE) signals observed experimentally, which are not consistent with the major state. Energy landscapes for p(D33K) and p(Y54A/L59A) are also constructed. The potential and free energy landscapes for these three systems are multifunneled and exhibit different relative stabilities for the experimentally determined fold, in line with available experimental data. Folding from a representative temperature-unfolded state to the native state is also analyzed, and key intermediates involved are highlighted and discussed.

Methods

Force Fields and Solvent Models

Parameters for phosphorylated Thr (pThr) were obtained from Case et al.,[38] and some atom types (see the Supporting Information (SI) for details) were modified to make them compatible with the ff14SB force field.[39] For explicit solvent simulations, RSFF2C (ff14SB plus residue-specific CMAP terms) was used with TIP3P water. This combination was previously validated for several folded proteins and IDPs.[40] For basin-hopping and discrete path sampling, an implicit solvent generalized Born (GB) model, GB-OBC[41] (igb = 2 in Amber(42)) was used and CMAP terms in RSFF2C were refitted for this solvent model to reproduce each residue’s intrinsic conformational preferences observed in a coil library, following a previously described procedure.[40]

In addition, we checked the appropriateness of interactions between phosphate groups and Thr, Arg, and Lys in the igb2 model. Inspired by Okur et al.’s work,[43] we used three model peptides containing pThr to check these interactions: Thr–Ala–Ala–pThr (TAApT), Arg–Ala–Ala–pThr (RAApT), Lys–Ala–Ala–pThr (KAApT), and another three well-studied peptides: Ser–Ala–Ala–Glu (SAAE), Arg–Ala–Ala–Glu (RAAE), and Lys–Ala–Ala–Glu (KAAE) for a consistency check. For each peptide, short REMD simulations (50 ns per replica) were conducted with TIP3P water (igb5 for the SAAE peptide) and the igb2 model, respectively. Potential of mean forces (PMFs) were calculated using histogram analysis along a specific coordinate (see the SI for details), following Okur et al.[43] Guided mainly by the PMFs in TIP3P water, we optimized the intrinsic Born radii of three atom types: hydroxyl hydrogen (HO), hydrogen atoms within Arg sidechains (HN+), and oxygen atoms (OX) within the phosphate group.

MD Simulations of pWT in TIP3P Water

For pWT, we conducted MD simulations in TIP3P water at four temperatures (300, 400, 450, and 480 K). The first model of the NMR ensemble determined experimentally[18] was extracted and capped with Ace (N-terminal) and Nme (C-terminal) using Pymol.[44] Then, it was solvated in a truncated octahedral box with 2802 TIP3P water molecules and 3 Na+ ions were added for charge neutralization. After energy minimization and 5 ns NPT equilibration at 300 K, a 2 μs NVT production run was conducted at 300 K. For each high temperature, we performed three independent 1 μs NVT production runs starting from the same structure but with different initial velocities.

Starting Structures and Basin-Hopping

Starting from each of the 20 structures (capped with Ace and Nme) in the NMR ensemble, 1000 steps of basin-hopping[35] were performed to locate lower-energy folded structures. The one with the lowest potential energy was selected to represent the folded state (F). For minor states, the system visited a vast set of minima with diverse structures during the high-temperature MD simulations. However, it would be computationally expensive to construct an energy landscape covering the whole of this conformational space. Instead, we extracted four intermediate states from the high-temperature simulation trajectories. Specifically, since we were interested in minor states that could potentially bind to eIF4E, we extracted all snapshots containing a C-terminal helix within the binding motif from the nine high-temperature MD trajectories. Then, based on the Cα root-mean-square deviation (RMSD) of residues 19–60, we conducted a clustering analysis for these structures using the hierarchical agglomerative approach implemented in Amber’s cpptraj. From the five most populated states (i.e., clusters) (see Figure S3), we selected three states with diverse secondary structure features: a state with four β-strands and a C-terminal helix (H1), another one similar to state H1 but with β1 almost detached from β4 (H2), and another where the detached β1 becomes an N-terminal helix (H3). We also chose a state without any terminal helix and with both β1 and the long loop detached from the other three β-strands (3β) because it was frequently observed as a metastable intermediate during high-temperature unfolding and was also proposed as an important folding intermediate in a previous computational study.[29] Finally, we extracted a representative unfolded state (U) at 480 K. A clustering analysis based on the Cα RMSD of residues 37–46 was conducted for the trajectories after unfolding. From the most populated cluster, the structure with the lowest backbone RMSD to the folded state was extracted. For each of the minor states (H1–3 and 3β) and the unfolded state (U), 1000 steps of basin-hopping were conducted to find locally more stable conformations (further steps should eventually lead to the global minimum).

Construction of Energy Landscapes via DPS

To survey the potential energy landscapes, for each of the five states (F, H1–3, and 3β), the structure with the lowest energy from basin-hopping was used as a starting point. The OPTIM program[36] was used to find initial discrete paths between each pair of starting structures. The initial stationary point databases were further expanded using the SHORTCUT and UNTRAP schemes in PATHSAMPLE.[37] SHORTCUT attempts to shorten the pathways and circumvent high energy barriers. The UNTRAP scheme[45] was used to remove artificial kinetic traps, and this step was conducted after removing some unphysical structures with incorrect hydrogen bonds involving phosphate groups from the databases (see Results section for details) to focus the sampling only on physical structures and pathways. For p(D33K) and p(Y54A/L59A), the starting structures were prepared by mutating the five starting structures of pWT by Pymol. The energy landscapes were constructed using the same procedure as above. For pWT, the folding pathway from the temperature-unfolded state (U) sampled during high-temperature simulations to the folded state was also constructed and optimized following a similar procedure. Convergence of the databases was monitored by checking interconversion rates between low-lying states.

After construction of the stationary point databases, free energy landscapes for three systems were calculated using the harmonic superposition approximation (HSA).[46] Free energy minima and transition states were recursively regrouped[47] using a threshold of 10 kcal/mol. Dijkstra’s shortest path analyses[48] with appropriate edge weights, involving the logarithm of the branching probability, were performed on the regrouped databases to find interconversion pathways making the largest contributions to the steady-state rate constants. The corresponding rate constants were calculated using a graph transformation approach.[49] Finally, the computed energy landscapes were visualized using disconnectivity graphs.[50,51]

Results and Discussion

Intrinsic Born Radii Optimization

In the NMR structure of pWT (PDB ID: 2MX4),[18] phosphate groups form stabilizing hydrogen bonds (HBs) with nearby hydroxyl groups of Thr, including pThr37 OX–Thr41 HO and pThr46 OX–Thr50 HO. Furthermore, because of the negative charge on the phosphate groups, they can potentially form salt bridges with positively charged Arg/Lys sidechains. We thus checked the appropriateness of these interactions in the igb2 model with small model peptides (TAApT, RAApT, and KAApT).

PMFs of the model peptides in implicit solvent igb2 were first calculated and compared to those in TIP3P water, revealing that the interactions between pThr and the three residues (Thr, Arg, and Lys) with the original intrinsic Born radii were too weak (Figure , in red). A simple fix would be to adjust the intrinsic Born radii of these interacting atoms, but we could not determine the origin of these weak interactions because each interaction involves two groups. We further checked the PMFs of three additional model peptides (SAAE, RAAE, and KAAE), and we found that some adjustments were needed to match the reference PMFs. Specifically, the interaction between Arg and Glu in RAAE was slightly stronger than the TIP3P results, and so the intrinsic Born radii of protons within the Arg sidechain (HN+) were decreased from 1.3 to 1.22 Å, while for Glu the radii were not altered because the PMF of KAAE in igb2 agrees well with that for TIP3P. The radius of Hγ within the Ser sidechain (HO) was increased from 1.2 to 1.35 Å to make the PMF of SAAE in igb2 match with that in igb5, which reproduced the PMF in TIP3P very well according to Nguyen et al.[52] Then, we fixed the radii of these two atom types and optimized the radii of phosphate oxygen (OX) in TAApT and RAApT, producing a value of 1.68 Å. Finally, the KAApT peptide was used to verify the radii modifications and the resulting PMF agrees well with that in TIP3P water. The PMFs after radii optimization are illustrated in Figure .

Figure 2

PMFs before (in red) and after (in green) intrinsic Born radii modifications compared to reference PMFs (in black). The references are PMFs in TIP3P water except for SAAE, where the PMF in igb5 was reported[52] to agree with TIP3P result and was used as the reference.

High-Temperature MD Simulations

At 300 K, pWT remained stable during the short 2 μs simulation with a backbone RMSD to NMR structure (PDB ID: 2MX4)[18] less than 1.5 Å most of the time (Figure S1). This result validates the force field/explicit solvent combination we employ. To encourage barrier crossing, high-temperature simulations were conducted at 400, 450, and 480 K to sample minor states. At 400 K (Figure ), complete unfolding was only observed in one trajectory (run B). In run A, the system was stable during the first 800 ns, except that at about 520 ns a C-terminal helix was formed within the binding motif and the angle between β3 and β4 grew (this state was also observed in run B at 400 K and run A at 450 K, see Figure S2 top). Then, β1 detached and formed a small helix. In the end, the β turn stabilized by pT46 was retained, and a small helix involving pT37 was formed. In run C, the system visited several metastable states before folding back to a native-like state. At 450 K (Figure S2 top), we also observed several states containing the C-terminal helix (run A) and some states containing two or three β-strands (run B and run C). At 480 K (Figure S2 bottom), complete unfolding occurred before 500 ns, but states with transient secondary structures and native contacts were formed later.

Figure 3

Fraction of native contacts during three independent explicit solvent MD simulation runs at 400 K. Some intermediate states are shown below.

During the high-temperature simulations, a transient α helix within the canonical binding motif was frequently observed. Based on examination of the complex structure,[27] it seems, with some structural re-arrangement prior to or upon binding, several states containing this C-terminal helix can bind to eIF4E, and some of them contain considerable folded structure. The formation of an α helix within the canonical binding motif prior to binding eIF4E can explain the moderate binding affinity of pWT to eIF4E. Although the existence of binding-competent minor states was suspected earlier, no direct observation of these states experimentally or computationally has been reported before, to the best of our knowledge.

Minor States and NOE Restraints

The unusually high number of NOE violations (about 5% of NOE restraints are completely unsatisfied by the folded state)[18] is a key piece of evidence supporting the existence of minor states. The largest cluster of NOE violations involves residues 42–45 and residues 31–34. In states H1 and H3 (as in Figure ), a large portion of these NOE violations are satisfied because of the extended pairing between β2 and β3. Two representative examples (P31-HA_S44-HA and D33-HB_T45-HN) are shown in Figure . It was previously suspected that alternative states with extended β2−β3 pairing are unlikely to retain β1−β4 pairing.[18] Here, we show that, in state H1, the β2−β3 pairing is longer than that in the folded sate F, and the longer β2−β3 pairing does not break the pairing between β1 and β4, but the β-strands are twisted because of the shortened loop region. In state H3, the detached β1 forms a helix, which is anchored close to pT46 by hydrogen bonds between R20 and pT46. However, no NOE signals between this helix and the turn region were observed in experiments, suggesting that this helix is probably transient and the linker between it and β2 is likely to be flexible.

Figure 5

Potential energy disconnectivity graphs (ΔE = 10 kcal/mol) for pWT (left), p(D33K) (middle), and p(Y54A/L59A) (right). Minima with unphysical HBs are highlighted in red. Representative structures for the principal funnels are shown below.

Figure 4

Examples of NOE restraints satisfied by minor states H1 (A) and H3 (B). Structures with the lowest potential energies from basin-hopping runs are shown. The structures are colored with the same rainbow scheme as in Figure , where β1, β2, β3, and β4 are roughly in blue, green, yellow, and orange, respectively.

Unphysical Hydrogen Bonds Involving Phosphate Groups

During discrete path sampling, we obtained some minima (red branches in Figure ) with an unphysically large number of hydrogen bonds (HBs) between the pThr sidechain and other residues, especially Arg. Similar structures were observed in our explicit solvent MD simulations and previously in REMD folding simulations.[33] Some snapshots were extracted and are shown in Figure S4. In the first two states (Figure S4a,b), pT37 forms HBs with three Arg residues and more complicated HB networks are formed in the latter two states (Figure S4c,d), with Arg sidechains bridging the two pThr residues together. In these states, at least one oxygen atom within the phosphate groups is involved in as many as five HBs. We note that it is not unusual for one pThr to form HBs with three Arg residues according to our survey in the Protein Data Bank (data not shown), but it is certainly unphysical for one oxygen atom in the phosphate group to form more than three HBs because it has at most three lone pairs of electrons. This problem is probably due to improper treatment of HBs in classical force fields, where directional HBs are approximated by nondirectional electrostatic and van der Waals interactions. This approximation seems appropriate for systems without complicated HB networks but more likely to cause problems for systems containing strong HB acceptors (such as phosphate groups) and multiple HB donors. These HBs can not only hinder sampling by trapping the system in unphysical states, which was probably responsible for the previous unsuccessful folding simulations,[33] but will also fail to reproduce the correct HB networks observed in phosphorylated proteins.[17]

To eliminate energetic frustration caused by these unphysical HBs from the energy landscapes, we removed all structures where any of the three terminal oxygen atoms in phosphate groups formed more than three HBs (structures with unphysical HBs are marked in red in Figure ). All further analyses are conducted for databases containing the remaining stationary points. Most of the incorrect HBs exist in the funnels containing states F, H3, and 3β. For p(Y54A/L59A), the funnel containing H3 actually disappeared after removing these incorrect structures. This result indicates how incorrect potentials can perturb the underlying energy landscape. Before removing the incorrect structures, the H3 state is separated from the other states by a high energy barrier, which is hard to cross during MD simulations. For more accurate simulations of phosphorylated systems, further force field improvements are needed.

Potential Energy Landscape for pWT

The potential energy landscape for pWT is multifunneled with the folded state F slightly more favorable than the four minor states (Figure ), resembling that of an intrinsically disordered protein.[1] Since the existence of an α helix within the binding motif may be indicative of the binding affinity to eIF4E, all minima with the binding motif in α helical conformation were highlighted based on results from DSSP calculations (Figure , left). For states H1 and H3, most low-lying minima in these two funnels have the C-terminal helix. A mixture of helical and nonhelical conformers is observed in the branch containing state H2, indicating the helix instability in this state. In the middle graph, minima are colored according to their backbone RMSD to the global minimum. As expected, both states H1 and H2 have relatively low RMSD to state F because they share a similar backbone conformation, while states H3 and 3β are more distant to state F in terms of RMSD because the first β-strand (β1) is dissociated from the other strands. We also compare the number of HBs between sidechains of the two pThr residues and other residues in different states (Figure , right). Although HBs involving pThr stabilize the folded state, more HBs involving pThr are formed in states 3β and H3. These additional stabilizing HBs compensate for the loss of backbone HBs due to the breaking of β1−β4 pairing, making them very close to the folded state in potential energy. The number of HBs involving pThr in states H1 and H2 are similar to the folded state, while fewer HBs are formed in high-lying minima outside the five major funnels. Overall, these results clearly demonstrate the stabilizing effect of HBs involving phosphate groups.

Figure 6

Potential energy disconnectivity graph (ΔE = 10 kcal/mol) for pWT with different color schemes. Left: minima with the binding motif forming an α helix are highlighted in red. Middle: colored by backbone RMSD to the global minimum (the minimum with the lowest energy in state F). Right: colored by the number of HBs formed between sidechains of the two pThrs and other residues.

Potential Energy Landscapes for the Two Mutants

Similar to pWT, the potential energy landscapes for p(D33K) and p(Y54A/L59A) are also multifunneled, but some significant differences exist among the three systems. As shown in Figure , for p(D33K), three states (F, 3β, and H3) have comparable potential energies, so the landscape is more frustrated than for pWT. This difference probably arises from the loss of HBs involving Asp33 on substituting Asp33 by Lys33. In pWT, Asp33 forms two or more HBs with other residues (including Tyr54) in state F, but only one (for states H1, H2, and H3) or none (for 3β) in the minor states (Figure B, left), indicating some contribution of this residue to the relative stability of the folded state. In p(D33K), Lys33 only forms HBs in state H1 (Figure B, middle). p(Y54A/L59A) exhibits a different pattern, and 3β has the lowest potential energy. This result agrees with the experimental observation that this mutant cannot fold but still retains both β-turns;[18] state 3β of this mutant is clearly stabilized by HBs involving phosphate groups (Figure A, right). State H3 has very high potential energy, leaving four principal funnels, rather than five for the other two systems. Furthermore, state F has relatively high potential energy, probably because of the loss of HBs between Tyr54 and other residues, which contribute to the stability of states H3 and F in pWT and p(D33K) (Figure S6).

Figure 7

Potential energy disconnectivity graphs (ΔE = 10 kcal/mol) for pWT (left), p(D33K) (middle), and p(Y54A/L59A) (right). (A) Minima are colored according to the number of hydrogen bonds formed between the two phosphate groups and other residues. (B) Minima are colored according to the number of hydrogen bonds between Asp33 (for pWT and p(Y54A/L59A)) or Lys33 (for p(D33K)) and other residues.

Free Energy Landscapes

For all three systems, the free energy landscapes are multifunneled, with states F, H1, and 3β lying very close in energy and states H2 and H3 lying higher (Figure ). State F has the lowest free energy for pWT, while state H1 is slightly more favorable for p(Y54A/L59A). Similar to the features observed in the potential energy landscapes, among the three systems, p(D33K) has the most frustrated landscape with four states (F, H1, H3, and 3β) having comparable energies, separated by considerable barriers.

Figure 8

Free energy landscapes (ΔE = 5 kcal/mol) for pWT (left), p(D33K) (middle), and p(Y54A/L59A) (right) at 298 K with a regrouping threshold[47] of 10 kcal/mol.

Conversion from Minor State 3β to the Folded State

For pWT, the fastest interconversion rate between the folded state and minor states is observed between state 3β and the folded state F. As shown in Figure , folding from state 3β starts with structural adjustments involving both termini to form an intermediate state where HBs between Tyr34/Tyr54 and backbone bring the long loop over the three β-strands. The second barrier corresponds to condensation of the N-terminus beside strand β4. Extending the β1−β4 pairing from this transition state follows a downhill path and, subsequently, fixes the long loop region. The rate constant calculated for conversion from state 3β to state F is 1.2 × 10–5 s–1 and that for the backward conversion is 1.5 × 10–6 s–1, suggesting that the interstate conversion is slow in the absence of binding partners. This mechanism is different from the pathway previously proposed by Gopi et al.,[29] where the long loop is completely formed before β1−β4 pairing. The difference can be explained by the existence of HBs between Arg20 and pThr46 in state 3β. Since these HBs already bring the N-terminus close to strand β4, initiation of β1−β4 pairing from this terminus is energetically more favorable compared to breaking this HB and winding the loop over the three β-strands. Furthermore, in our mechanism, the final formation of the long loop actually synchronizes with the extension of β1−β4 pairing, and thus seems more plausible, considering the compensation between enthalpy gain and entropy loss. The role of Tyr54 during the conversion is revealed in atomic detail: it forms a hydrogen bond with Gln29 through the hydroxyl group in the sidechain, which helps to fix the loop region. It also stabilizes the folded state through hydrogen bonding with Asp33. This mechanism can explain the inability of mutant p(Y54A/L59A) to fold and the low stability of the corresponding folded state.

Figure 9

Fastest conversion path from state 3β to the folded state F for pWT. Stationary points in the kinetic transition network were regrouped based on a threshold of 10 kcal/mol. Representative states are shown with key residues in sticks.

Folding from a Temperature-Unfolded State to the Folded State

We also studied the folding mechanism from a temperature-unfolded state to the folded state for pWT. The temperature-unfolded state (U) was extracted from high-temperature MD simulations at 480 K, as described in Methods section. The fastest folding pathway is presented in Figure A with selected states shown along the path. The pathway can be roughly divided into two stages. The system undergoes several structural transformations before reaching step 7, where the two β-turns are almost formed. Intermediate states in this stage have high energies and no obvious secondary structure. The highest barrier during folding is observed at the beginning of this stage. From step 7, the system goes through a transition state where the β2 and β3 form partial pairing and then reaches a state with three β-strands following a downhill path. This structure is similar to the minor state 3β described earlier, with HBs between Arg20 and pThr46. Folding from this point to the folded state thus follows a similar pathway, as described in Figure . The free energy difference between the temperature-unfolded state and the folded state is 31.8 kcal/mol, much larger than the experimentally estimated unfolding free energy (2.61 kcal/mol) for the full-length pT37pT46.[18] This difference is understandable because the unfolded state used here was sampled in high-temperature simulations and does not correspond to the experimentally assumed unfolded state (i.e., an intrinsically disordered ensemble) at room temperature. With this effect taken into consideration, the differences between the folding mechanism and that proposed by Gopi et al.[29] can also be explained: here, we find that, among the 4 β-strands, β2−β3 pairing forms first, while Gopi et al. suggested that β3−β4 pairing was energetically more favorable. Apart from the alternative starting unfolded state, another difference may lie in the different potentials used (i.e., atomic model vs coarse-grained model at the residue level).

Figure 10

(A) Fastest folding pathway from a temperature-unfolded state (U) to the folded state (F) for pWT at 293 K (experimental temperature). Stationary points in the kinetic transition network were regrouped[47] based on a threshold of 10 kcal/mol. The number of steps corresponds to the number of free energy minima along the path. Representative structures of selected sates are shown. (B) Number of native HBs (blue) and nonnative HBs (red) formed between the two phosphate groups and other residues along the pathway.

We computed the number of HBs formed between pThr sidechains and other residues in representative structures. HBs observed in the folded state (global minimum in the stationary database) are defined as native HBs, while the others are nonnative. Figure B shows how the number of native and nonnative HBs evolve along the folding process. Many nonnative HBs involving phosphate groups are observed in most states along the path, except for those close to the folded state and TS5. Formation of several states, including steps 1, 5, 8, and 11 (i.e., the folded state), are accompanied by an increase of HBs involving pThr sidechains, while most transition states involve a decrease in nonnative HBs (Figure A). The ubiquitous HBs involving phosphate groups along the folding pathway indicates that they are important in determining the folding of this protein.

Characteristics of the Heat Capacities

From the databases of minima, we also calculated heat capacity features (Figure ) for the three systems in the harmonic approximation and identified the primary minima that contribute to the Cv peaks using the temperature gradient of the occupation probabilities.[53] For pWT, Cv exhibits a small peak at 199 K and a larger peak at 315 K in the temperature range from 50 to 500 K. The first peak corresponds to interconversion of minima within state F, while the second peak is related to the conversion from state F to state H1 (Figure S7). The experimental temperature (293 K) is located within the second peak, where major state F is dominant and the binding-competent minor state H1 is accessible. This result explains both the relative stability of the folded state and the moderate binding affinity of pWT to eIF4E. Only one sharp peak (at 241 K) was observed for p(Y54A/L59A) and it corresponds to interconversion between state 3β and state H1 (Figure S9). This result suggests the existence of a mixture of the less ordered state (3β) and the more ordered state (H1) at the experimental temperature, which is possible considering the partial structural order revealed by NMR experiments.[18] The Cv curve of p(D33K) also has two peaks, but they have similar heights and correspond to interconversion between states. The first (at 135 K) peak corresponds to interconversion of state H3 and states F and 3β, while the second peak (at 216 K), overlapping with the only peak for p(Y54A/L59A), corresponds to interconversion between these three states and H1 (Figure S8). With no experimental data for this mutant available for comparison, we predict p(D33K) to be more structurally heterogeneous at room temperatures than the other two systems based on the Cv features.

Figure 11

Heat capacity curves for pWT (black), p(D33K) (blue), and p(Y54A/L59A) (red).

Conclusions

Doubly phosphorylated 4E-BP2, pT37pT46, adopts a marginally stable four-stranded β-sheet fold with eIF4E-binding-competent minor states and exhibits great structural sensitivity to mutations. This study combined high-temperature MD simulations and discrete path sampling[36,37] to elucidate these properties. High-temperature MD simulations sampled several minor states with the binding motif (YxxxxLΦ) folded into a small α helix. These α helix containing minor states can explain the moderate binding affinity to eIF4E, and among them, two minor states also explain some NOE signals that cannot be assigned to the folded state. Energy landscape explorations further revealed that these minor states lie close to the folded state (F) in energy with stabilizing contributions from hydrogen bonds (HBs) involving the two phosphate groups. The fastest interconversion was observed between the folded state and a state with only three β-strands (state 3β). The interconversion pathway is different from that suggested by Gopi et al.,[29] indicating that alternative folding mechanisms can be adopted by the system depending on the initial state in question.

Two mutants of the doubly phosphorylated 4E-BP2 were also examined in this study. Both p(D33K) and p(Y54A/L59A) have multifunnel energy landscapes, but, among the three variants, p(D33K) has the most frustrated landscape, with several states lying close in energy, while p(Y54A/L59A) prefers the three-stranded β-sheet (3β), in line with experimental observations.[18]

For pWT, we also studied the folding pathway from a temperature-unfolded state extracted from MD simulations at 480 K to the folded state. Although the unfolded state is not representative of the unfolded ensemble for this system at room temperature, the pathway reveals the important role of phosphate groups during folding. They stabilize intermediate states and are also, at least partially, responsible for the energy barriers along the pathway by forming hydrogen bonds with other residues.

During both MD simulations and discrete path sampling, we observed an unphysically large number of HBs between phosphate groups and other residues, especially Arg. In some structures, one single oxygen atom in the phosphate group can form as many as five HBs, indicating an improper treatment of HBs within classical force fields. As a consequence, the energy landscapes were contaminated by artificial frustration (low-lying conformations with incorrect HBs). These incorrect structures can be removed systematically from the stationary point databases, but we cannot see a straightforward way to avoid them during MD simulations unless some modifications are made to the force fields. Specific modifications would be worthwhile considering the abundance and importance of phosphorylation.