1. Biophysics and Structural Biology
Download icon

Changes in the free-energy landscape of p38α MAP kinase through its canonical activation and binding events as studied by enhanced molecular dynamics simulations

  1. Antonija Kuzmanic
  2. Ludovico Sutto
  3. Giorgio Saladino
  4. Angel R Nebreda
  5. Francesco L Gervasio Is a corresponding author
  6. Modesto Orozco Is a corresponding author
  1. The Barcelona Institute of Science and Technology, Spain
  2. University College London, United Kingdom
  3. Catalan Institution for Research and Advanced Studies (ICREA), Spain
  4. Joint BSC-CRG-IRB Program in Computational Biology, Spain
  5. University of Barcelona, Spain
Research Article
Cited
0
Views
807
Comments
0
Cite as: eLife 2017;6:e22175 doi: 10.7554/eLife.22175

Abstract

p38α is a Ser/Thr protein kinase involved in a variety of cellular processes and pathological conditions, which makes it a promising pharmacological target. Although the activity of the enzyme is highly regulated, its molecular mechanism of activation remains largely unexplained, even after decades of research. By using state-of-the-art molecular dynamics simulations, we decipher the key elements of the complex molecular mechanism refined by evolution to allow for a fine tuning of p38α kinase activity. Our study describes for the first time the molecular effects of different regulators of the enzymatic activity, and provides an integrative picture of the activation mechanism that explains the seemingly contradictory X-ray and NMR data.

https://doi.org/10.7554/eLife.22175.001

Introduction

p38 Ser/Thr kinases are mitogen-activated protein kinases (MAPKs) involved in the regulation of multiple cellular processes, including cell proliferation, differentiation, senescence, and death (Cuadrado and Nebreda, 2010). The p38 subgroup of MAPKs comprises four members (α-δ) where only p38α is expressed ubiquitously at high levels (Cuadrado and Nebreda, 2010). p38 signaling is strongly activated by environmental stresses (e.g. osmotic shock, ionizing radiation), and biological stimuli, such as growth factors and inflammatory cytokines (Trempolec et al., 2013). Furthermore, p38α has been implicated in several pathological conditions, for example, chronic inflammatory diseases, cancer, and heart and neurodegenerative diseases (Yong et al., 2009, Denise Martin et al., 2012), which is why elucidation of its activation mechanism is of therapeutical importance (Yong et al., 2009; Hammaker and Firestein, 2010).

Like most kinases, p38α is composed of two lobes: the smaller N-terminal lobe, consisting mostly of β-sheets, and the α-helical C-terminal lobe. The lobes are linked by a flexible hinge that forms the ATP-binding site together with structural elements from both lobes (marked regions in Figure 1). In the canonical activation pathway, MAPK kinases (MAP2Ks) dually phosphorylate the TGY sequence in the activation loop (A-loop; Figure 1a and Video 1) which, according to X-ray crystallography (Zhang et al., 2011a), triggers its large conformational rearrangements and the formation of a characteristic β-sheet motif away from the ATP- and substrate-binding sites (Figure 1b). The new position of the A-loop brings N- and C-lobes closer and facilitates the reorientation of key residues which participate in the stabilization of ATP and the catalytically crucial Mg2+ ions (Taylor and Kornev, 2011) (i.e. the universally conserved K53, E71, and D168). Moreover, p38α and other MAPKs possess two distinct structural elements (Figure 1): (1) the MAPK insert (MKI), which forms a lipid-binding site (Diskin et al., 2008) and (2) the L16 loop that extends from the C-lobe to the N-lobe where it folds into the C-terminal L16 helix. At its C-lobe segment, the L16 loop contains an acidic patch, called the common docking (CD) motif (Tanoue et al., 2000), which together with the hydrophobic groove, formed by the linker between αD and αE helices, and the β78 reverse turn (also termed the ED site) (Chang et al., 2002) (Figure 1a), defines the peptide docking site recognized by the conserved docking motif present in activators, substrates, and regulators of MAPKs. The docking site enhances the specificity of MAPKs’ interactions, as well as its activity, albeit through an unknown molecular mechanism (Tanoue et al., 2000). All the structural elements are detailed in the Video 1.

Structural models of (a) the inactive (PDB ID: 3S3I), and (b) the active, dually phosphorylated, p38α (PDB ID: 3PY3).

Main structural elements are colored as follows: A-loop in purple (inactive) and green (active), αC helix in red, L16 loop in slate, P-loop in cyan, hinge in yellow, and MKI in orange. Key residues are shown as sticks and labeled appropriately. The ATP-binding site is indicated by the grey filled ellipse, while the αE helix, and CD and ED sites by the green dashed ellipse.

https://doi.org/10.7554/eLife.22175.002
Video 1
Detailed presentation of p38α structural elements and binding sites which play an important role in its activation.
https://doi.org/10.7554/eLife.22175.003

The canonical activation mechanism outlined above, derived from the analysis of phosphorylated and unphosphorylated crystal structures of p38α (Zhang et al., 2011a) has recently been challenged by an NMR study of p38α in solution (Tokunaga et al., 2014), which surprisingly showed that the dual phosphorylation of the A-loop induces only small conformational changes to p38α, raising doubts about the reliability of conformational states associated with the canonical activation. To understand the changes that take place upon p38α phosphorylation and to explain the observed discrepancies between the NMR and X-ray results, we employed state-of-the-art simulation techniques. Specifically, we used unbiased molecular dynamics (MD) simulations and parallel tempering metadynamics (PT-metaD) (Bussi et al., 2006; Bonomi and Parrinello, 2010; Sutto et al., 2012), a method that combines an enhanced sampling technique with a multi-replica approach and is thus able to converge very complex conformational free-energy surfaces as a function of a few relevant variables, as has been shown previously for several protein kinases (Berteotti et al., 2009; Lovera et al., 2012; Sutto and Gervasio, 2013; Lovera et al., 2015; Marino et al., 2015) and other proteins (Bonomi et al., 2007; Papaleo et al., 2014; Lambrughi et al., 2016). In total, we obtained a high level of sampling (>115 μs) which allowed us to explore the molecular mechanism of p38α canonical activation in unprecedented detail and reconcile the apparently contradictory X-ray and NMR data.

Results

We performed PT-metaD in the well-tempered ensemble to obtain free-energy surfaces (FESs) for human p38α in four different states: (1) the unphosphorylated, (2) the apo dually phosphorylated, (3) the dually phosphorylated state with bound ATP-Mg2+, and (4) the dually phosphorylated state with bound ATP-Mg2+ and MK2 docking peptide (Figure 2 and Figure 2—figure supplement 1), obtaining total sampling times of 24.25, 24.85, 10.35, and 9.86 μs, respectively. To analyze the systems, we projected the FESs along two collective variables (CVs) - CV1 and CV2, both of which are based on the A-loop contacts that measure the distance from the crystallographically observed inactive (CV1) and active (CV2) conformations (see Materials and methods), and which were successfully used previously to study protein kinases (Sutto et al., 2012; Sutto and Gervasio, 2013; Marino et al., 2015). When CV1 < 0.2 and CV2 > 0.8, then the system assumes an inactive conformation (as in Figure 1a). On the other hand, when CV1 > 0.6 and CV2 < 0.2, the A-loop has formed the characteristic β-sheet motif (as in Figure 1b) and the enzyme has adopted the canonical active state. As the representative structure of each highlighted minimum, we chose the central structure of the most populated cluster present in the minimum (see Materials and methods for further details).

Figure 2 with 6 supplements see all
Free-energy surface for p38α (apo unphosphorylated, apo dually phosphorylated, dually phosphorylated with bound ATP-Mg2+, and dually phosphorylated with bound ATP-Mg2+ and MK2 peptide) as a function of CV1 and CV2 which indicate the distance from the reference inactive and active structures, respectively.

The contour lines are drawn every two kcal mol−1. Representative structures for the selected minima are shown and colored as follows: A-loop in purple, αC helix in red, L16 loop in slate, P-loop in cyan, and hinge in yellow. Key residues are shown as sticks and labeled appropriately. Part of the FES of the phosphorylated p38α that is populated by structures which are likely to bind ATP are highlighted with a grey rectangle (upper right panel).

https://doi.org/10.7554/eLife.22175.004

Unphosphorylated p38α

FES of the unphosphorylated state (upper left panel, Figure 2) and the representative structures of the explored minima (Figure 2 and Figure 2—figure supplement 1, enlarged structures in Figure 2—figure supplement 2) show that the unphosphorylated protein neither explores the active conformation, nor does it populate a state able to accommodate ATP. Specifically, in the more inactive of the two free-energy minima - (CV1, CV2) = (0.2, 0.7) (upper left panel in Figure 2, Figure 2—figure supplement 2a), the A-loop is positioned over the C-lobe forming contacts with residues involved in the substrate binding and phosphotransfer (Table 1). The P-loop is quite flexible and the K53-E71 salt bridge, which stabilizes ATP, is displaced from the C-lobe due to the broken regulatory spine that is otherwise formed with D168 upon ATP binding (Taylor and Kornev, 2011). Furthermore, E71 interacts with the MAPK-conserved R67 from the αC helix that stabilizes pT180 in the active form (Figure 1b, present in >90% of structures in the most populated cluster). On the other hand, the A-loop of the second minimum contains a short helix which places the conserved R173 in direct contact with the K53-E71 salt bridge, as well as D168 of the conserved DFG motif (Figure 2, Figure 2—figure supplement 2b), while the P-loop is partly collapsed into the ATP-binding site. We observe similar autoinhibitory interactions for the other two higher energy minima explored by p38α (Supplementary file 1, Figure 2—figure supplement 2c,d). In both cases, R173 is in contact with D168 which otherwise stabilizes ATP upon binding (Kumar et al., 1995; Bellon et al., 1999). Furthermore, in the more inactive minimum (Figure 2—figure supplement 2c), the P-loop is collapsed into the ATP-binding site; while the other explores a conformation reminiscent of the Src-like inactive conformation (Xu et al., 1999) (Figure 2—figure supplement 2d) in which the A-loop forms a short helix stacked against the αC helix. This configuration has been previously observed also in CDK (De Bondt et al., 1993) and EGFR kinases (Sutto and Gervasio, 2013). These findings are in perfect agreement with experimental studies reporting that the unphosphorylated p38α shows almost no affinity for ATP (Tokunaga et al., 2014; Frantz et al., 1998), confirming the strong regulation of the enzyme by phosphorylation.

Table 1

Interactions of key residues observed in the free-energy minima. Residue pairs with interaction occupancy >75% in the most populated clusters of the selected minima are bolded, while the occupancy for the rest of the pairs is in the range of 60–75%.

https://doi.org/10.7554/eLife.22175.011
SystemMinimumResidue pairs
p38α0–2 kcal mol−1
(CV1, CV2 = 0.2, 0.7)
D176-R189, D177-R186, E178-K152, E178-R173
0–2 kcal mol−1
(CV1, CV2 = 0.2, 0.6)
R173-E71, D168-R173, D177-R186
p38α-pTpY0–2 kcal mol−1E178-R149, E178-R189, D177-R67, pY182-R173
4–6 kcal mol−1
(CV1, CV2 = 0.45, 0.4)
D176-R149, D176-R189, D177-R70, pT180-R67, R173-E328
p38α-pTpY (ATP-Mg2+)0–2 kcal mol−1pT180-R67, pT180-R70, pY182-R186, pY182-R189, E178-R70
p38α-pTpY (ATP-Mg2+ + MK2 peptide)0–2 kcal mol−1
(CV1, CV2 = 0.7, 0.1)
pT180-R67, pT180-R173, D112-ATP
0–2 kcal mol−1
(CV1, CV2 = 0.4, 0.4)
pY182-R67, pY182-R70, D177-R173, D112-ATP

To confirm that the extended conformation of the A-loop attributed to the active state is highly unstable in the unphosphorylated protein, we used unbiased MD simulations that were started from the conformation observed in the X-ray structure of p38α-pTpY (PDB ID: 3PY3) (Figure 1b) (see Materials and methods). We performed ten 500-ns-long independent runs at 380 K, all of which quickly showed a change in the A-loop conformation toward the inactive state and sampled conformations that fall in the regions explored in metadynamics simulations. Furthermore, in one of the independent runs the A-loop managed to break the β-sheet motif in the simulated time, so we extended this run to 2.675 μs to allow the protein to transition to the inactive state (Figure 3a). Indeed, the protein mostly remained in the global minima observed in metadynamics simulations, and the final conformation showed very similar features to the representative inactive structures (Figure 3b). Root-mean-square fluctuations (RMSF) reveal a high flexibility of several structural elements besides the A-loop (i.e. αC, αD, and the MKI helices, Figure 3b) that occurs due to the opening of the lobes and the local unfolding (in the case of αD helix). Both these events have been previously linked to large conformational changes occurring in other kinases (Sutto and Gervasio, 2013; Whitford et al., 2007; Shan et al., 2013) but have been so far unobserved for p38α.

(a) Free-energy surface for the unphosphorylated p38α, obtained from metadynamics simulations, as a function of CV1 and CV2 which indicate the distance from the reference inactive and active structures, respectively.

The contour lines are drawn every two kcal mol−1. The CV values for the unbiased MD simulation of the same system are drawn every 0.5 ns and shown in a grey gradient as a function of time. (b) The conformation assumed by the protein at the end of the unbiased MD simulation. The backbone is colored according to the RMSF values, with purple indicating higher flexibility regions, while the key residues are shown as sticks and labeled appropriately.

https://doi.org/10.7554/eLife.22175.012

Dually phosphorylated p38α

Surprisingly, once p38α is dually phosphorylated, it explores conformations that are not largely different from the ones observed for the unphosphorylated protein (upper right panel in Figure 2 and Figure 2—figure supplement 1, enlarged structures in Figure 2—figure supplement 3), contrary to the reported X-ray data (Zhang et al., 2011a). In the global minimum, the A-loop moves away from the ATP-binding site, mostly due to the electrostatic interactions formed by the phosphorylated residues which subsequently also stabilize the P-loop and prevent it from collapsing into the ATP binding pocket (Table 1). However, the sampled conformation is unlikely to bind ATP due to the formation of a salt bridge between R49 (positioned beneath the P-loop) and D112 (found at the beginning of the αD helix) (Figure 2, Figure 2—figure supplement 3a; distance histogram is shown in Figure 2—figure supplement 4). This interaction, which seems to be specific for the phosphorylated protein, locks the hinge and occludes the ATP-binding site. However, the protein also explores a higher energy minimum that comes closer to the active one (Figure 2, Figure 2—figure supplement 3b): the A-loop extends toward the L16 loop, yet does not form the β-sheet motif associated with the active conformation. Regardless, the higher energy conformation allows the A-loop to form a number of electrostatic interactions with the αC helix and the L16 loop (Table 1), while pT180 forms a salt bridge with the highly conserved R67 preventing the P-loop from collapsing into the ATP-binding site. Furthermore, the αD helix partly unfolds allowing the ATP-binding pocket to assume a more open conformation with the key residues aligned for ATP binding. This observation explains the protein’s low affinity toward ATP in the phosphorylated state (Tokunaga et al., 2014; Frantz et al., 1998) and demonstrates that the dual phosphorylation can help define a receptive ATP-binding site but is not sufficient for the protein to adopt the conformation captured by X-ray crystallography (Figure 1b). In fact, to populate that state, we would need to add 16–18 kcal mol−1 of free energy to the phosphorylated system (minimum located at (CV1, CV2) = (0.6, 0.2) in the upper right panel of Figure 2, Figure 2—figure supplement 3d).

To further confirm the instability of the conformation observed by X-ray crystallography, we again used unbiased MD simulations that were started from the said conformation (PDB ID: 3PY3) (Figure 1b) (see Materials and methods). We performed ten 1-μs-long independent runs at 380 K, most of which showed a change in the A-loop conformation toward the regions explored by the protein in metadynamics simulations (upper right panel in Figure 2). Furthermore, in one of the independent runs the A-loop broke the β-sheet motif, so we extended this run to 2 μs during which the protein settled in the aforementioned higher energy minima that seems likely to bind ATP (Figure 4a). The final conformation shows similar features to the representative structure of this minimum, that is, the A-loop is positioned away from the ATP-binding site and forms several electrostatic interactions with the αc helix and the L16 loop (Figure 4b). RMSF values reveal a high flexibility of several structural elements besides the loops (i.e. αD, and the MKI helices, Figure 4b) that occurs due to local unfolding.

(a) Free-energy surface for the dually phosphorylated p38α, obtained from metadynamics simulations, as a function of CV1 and CV2 which indicate the distance from the reference inactive and active structures, respectively.

The contour lines are drawn every two kcal mol−1. The CV values for the unbiased MD simulation of the same system are drawn every 0.5 ns and shown in a grey gradient as a function of time. (b) The conformation assumed by the protein at the end of the unbiased MD simulation. The backbone is colored according to the RMSF values, with purple indicating higher flexibility regions, while the key residues are shown as sticks and labeled appropriately.

https://doi.org/10.7554/eLife.22175.013

Dually phosphorylated p38α with bound ATP-Mg2+

The FES plot for the phosphorylated system with bound ATP-Mg2+ complex shows a clear shift toward a more active conformation (i.e. a displacement toward right and bottom in the lower left panel of Figure 2 and Figure 2—figure supplement 1, enlarged structures in Figure 2—figure supplement 5). While the representative structure of the global minimum does not form the β-sheet motif associated with the active state (Figure 2—figure supplement 5a), its A-loop extends toward the L16 loop and connects to the αC helix via the coordination of pT180 by R67 and R70 (as observed in the active X-ray structure in Figure 1b). The β-sheet motif does form in the higher energy minimum centered at (CV1, CV2) = (0.6, 0.1); however, its representative structure (lower left panel in Figure 2, Figure 2—figure supplement 5b) shows pT180 forming less contacts with the αC helix and the L16 loop (Supplementary file 1). The lack of these interactions, together with the disorder present in the P-loop, breaks the contact of E71 with Mg2+ ions resulting in a less stable complex with ATP. Overall, albeit the binding of ATP stabilizes an active-like structure, it is not enough for the protein to display the active conformation observed in the crystal structure (Figure 1b). While the A-loop does adopt an extended conformation and forms the β-sheet motif at the higher energy state, it still remains quite mobile and the phosphorylated residues do not necessarily form all the contacts observed in the X-ray structures.

Once again, we used unbiased MD simulations to give further support to our free-energy surfaces. We performed ten 1-μs-long independent runs at 380 K of the dually phosphorylated p38α with bound ATP-Mg2+ starting with the A-loop in the extended conformation (as in Figure 1b) (see Materials and methods). Most of the simulations explored only the local minimum close to the starting structure, while one of the runs showed more conformational diversity, so we extended it to 1.55 μs which allowed the protein to transition to the global minimum (Figure 5a). The final conformation shows the disruption of the β-sheet motif in the A-loop and the calculated RMSFs reveal its high flexibility (Figure 5b). Furthermore, pT180 remains in contact with the conserved arginine residues of the αC helix, while pY182 moves away from H228, similarly to the representative structure of the global minimum in the free-energy surface. Overall, RMSF values show a lower flexibility of the protein upon the addition of the ATP compared to the simulations of the apo systems (Figures 3b and 4b).

(a) Free-energy surface for the dually phosphorylated p38α with bound ATP-Mg2+, obtained from metadynamics simulations, as a function of CV1 and CV2 which indicate the distance from the reference inactive and active structures, respectively.

The contour lines are drawn every two kcal mol−1. The CV values for the unbiased MD simulation of the same system are drawn every 0.5 ns and shown in a grey gradient as a function of time. (b) The conformation assumed by the protein at the end of the unbiased MD simulation. The backbone is colored according to the RMSF values, with purple indicating higher flexibility regions, while the key residues are shown as sticks and labeled appropriately.

https://doi.org/10.7554/eLife.22175.014

The active conformation

While the dual phosphorylation and the ATP-binding drive the enzyme toward an active conformation, there is an obvious need for additional elements to complete the transition. Particularly, we wondered whether the binding of the MK2 docking peptide could provide the necessary stability to form the β-sheet motif. The FES for this system indeed shows a deepening of the minimum at (CV1, CV2) = (0.7, 0.1) which is associated with the formed β-sheet motif (lower right panel in Figure 2, enlarged structures in Figure 2—figure supplement 6). The representative structure of the minimum (Figure 2—figure supplement 6a) confirms this observation and shows pT180 near the three conserved arginines (R67, R70, and R173), while pY182 again displays higher mobility and moves away from H228 (Table 1, Figure 1b). Furthermore, the addition of the MK2 peptide does not eliminate the global minimum at (CV1, CV2) = (0.4, 0.4), observed also in the presence of ATP (two lower panels in Figure 2). However, it flattens the FES which would indicate an easier transition between the two described minima. The representative structure of the minimum (Figure 2—figure supplement 6b) shows a conformation in which the A-loop is extended away from the ATP-binding site, yet without the formed β-sheet motif and with switched positions of phosphorylated residues compared to the X-ray structure (Table 1, Figure 1b).

In addition, we performed unbiased MD simulations of the dually phosphorylated p38α with bound ATP-Mg2+ and the MK2 docking peptide. We generated ten 1-μs-long independent runs at 380 K, obtaining a total accumulated sampling time of 10 μs (see Materials and methods). The protein remains mostly in the regions associated with the extended conformation of the A-loop with formed β-sheet motif, but it also explores the conformational space closer to the other global minimum (also the global minimum in the system without the MK2 peptide) (Figure 6a). Interestingly, the N-terminal part of the MK2 docking peptide binds the ED site and stabilizes the αDE linker hindering the flexibility of the αD helix and allowing the D112 to reorient and stabilize ATP by forming hydrogen bonds with ribose (Figure 2—figure supplement 6, Figure 6b). While we also observe this interaction in the absence of the MK2 peptide, it is not as stable as in its presence which explains the phosphotransfer enhancement upon MK2 peptide docking (Tokunaga et al., 2014).

(a) Free-energy surface for the dually phosphorylated p38α with bound ATP-Mg2+ and MK2 docking peptide, obtained from metadynamics simulations, as a function of CV1 and CV2 which indicate the distance from the reference inactive and active structures, respectively.

The contour lines are drawn every two kcal mol−1. The CV values (calculated every 0.5 ns) for the unbiased MD simulations of the same system are shown as orange contours based on their density. (b) Histograms of the minimal distances between ATP and D112 calculated from the unbiased simulations of the dually phosphorylated p38α with bound ATP-Mg2+ with or without the MK2 docking peptide (colored grey and red, respectively). The average values are shown as dashed lines, while the bin size is 0.5 Å.

https://doi.org/10.7554/eLife.22175.015

Our results suggest that p38α predominantly samples conformations with the fully extended A-loop only upon double phosphorylation and the binding of ATP and the docking peptide. This might seem in contrast to the model obtained from X-ray crystallography (Figure 1b). However, careful analysis of the crystal contacts (which form between the molecules due to their orientation and packing within the crystal, Figure 7) reveals numerous interactions formed between functionally important structural elements, such as the L16 loop. These contacts mimic the stabilizing effects of p38α binders. Particularly, we found an atypically long expression His-tag of a neighboring molecule bound to the hydrophobic docking groove which, upon closer inspection (see below), has the same effect on the protein as the MK2 docking peptide did in our simulations.

Contact map for interactions formed between the residues of the reported X-ray structure (PDB ID: 3PY3, p38α-pTpY) (x-axis) and its symmetry related neighbors (y-axis).

The residues are considered to be in contact if the minimal distance between any of their atoms is under 4 Å. The contacts formed by the docked elements are indicated by the solid ellipses, while the ones formed between functionally important elements are indicated by the dashed ellipses. Secondary structure elements are also shown with α-helices in purple, β-sheets in green, and loops in grey. Missing residues are shown as light grey rectangles.

https://doi.org/10.7554/eLife.22175.016

To prove that this is the case, we simulated a crystal supercell of the apo dually phosphorylated p38α composed of 4 unit cells and 16 protein molecules (see Materials and methods), which allowed us to maintain crystal contacts and obtain 8 μs of sampling for ~170,000 atoms. As shown by Figure 8, the A-loop is very stable in the conformation captured by the X-ray model - the β-sheet motif is unperturbed during the simulation, as well as the contacts made by pT180 with the three conserved arginines (R67, R70, and R173). However, pY182 is more mobile and at times breaks away from H228 and the conserved arginines of the P+1 loop (R186 and R189). This is consistent with the electron density which is not so well defined for pY182 as it is for the rest of the A-loop. Furthermore, we could also qualitatively match the B-factors reported for the X-ray structure (PDB ID: 3PY3) (Figure 8b,d), while the most mobile regions in our simulations (mainly, the N- and C-terminal tails) are missing from the deposited electron density. We would like to stress that it is impossible for our MD simulations to quantitatively match the reported B-factors for the following reasons: (1) While the crystals were grown at 293 K, the data were collected at 100 K. Since MD simulations and their algorithms have not been developed for such low temperatures, we performed our crystal simulations at 293 K which increased the overall atomic motions; (2) The experimental data collected for the crystals in general are time and ensemble averages of 1023 molecules which are in the end represented by a single model, while its B-factors then reflect both the static and the dynamic disorder of the crystal. Finally, two independent studies had shown that the refinement process underestimates the disorder present in the crystal thereby producing lower B-factors (Janowski et al., 2013; Kuzmanic et al., 2014).

(a) Free-energy surface for the dually phosphorylated p38α, obtained from metadynamics simulations, as a function of CV1 and CV2 which indicate the distance from the reference inactive and active structures, respectively.

The contour lines are drawn every two kcal mol−1. The CV values (calculated every 0.5 ns) for the unbiased MD simulations of the 3PY3 crystal are shown as orange contours based on their density. (b) The average conformation of the protein calculated from the unbiased MD simulations. The protein is colored according to the calculated RMSF values after the backbone alignment with respect to the 3PY3 X-ray structure. (c) Backbone RMSD values calculated with respect to the 3PY3 X-ray structure either with or without N- and C-terminal tails. The average values over 16 molecules are shown by the black lines, the extreme values by the red lines, while the error bars depicting one standard deviation are shown every 25 ns. (d) X-ray structure of the dually phosphorylated p38α (PDB ID: 3PY3) colored by refined B-factors.

https://doi.org/10.7554/eLife.22175.017

Thus, our analysis and simulations strongly suggest that the deposited crystal structure does not fully capture the apo conformation that p38α assumes under physiological conditions and it most likely suffers from artifacts due to the crystal contacts.

Discussion

In the clever paraphrase of Tolstoy’s Anna Karenina, Noble et al. (Noble et al., 2004) stated how all active kinases are alike but an inactive kinase is inactive after its own fashion - referring to the attractiveness of the inactive form in the everlasting search for inhibitor selectivity. Kinases, however, seem to have a few surprises left within their folds, even when it comes to the seemingly uniform active state. In the case of p38α, the protein requires dual phosphorylation of the TGY motif for its activation and subsequent ATP binding. Based on the apo X-ray structure of p38α-pTpY (Zhang et al., 2011b), the phosphorylation event has been perceived as the trigger for the conformational change in which the A-loop moves away from the ATP-binding site enabling the N- and C-lobes to reorient and align the residues necessary for the ATP stabilization. Unlike CDK5 (Berteotti et al., 2009), EGFR (Sutto and Gervasio, 2013), and Src (Shukla et al., 2014; Foda et al., 2015), p38α never explores the αC-out conformation (neither in any of the deposited X-ray structures - Figure 9a nor in our simulations), in which the helix rotates by >45° (breaking the K53-E71 salt bridge in the process). This likely occurs due to the extensive range of contacts formed between the αC and the L16 helix, which is unique to MAPKs, indicating that the observed feature could be family-specific. Furthermore, MAP kinases contain two conserved arginine residues (R67 and R70) in the αC helix that readily form salt bridges with numerous negatively charged residues in the A-loop (Figure 10), thereby stabilizing various A-loop conformations which are not necessarily observed by X-ray crystallography. Despite an overwhelming amount of structural data (with 206 structures of human p38α in the PDB), just 47 structures have complete A-loops (i.e. without any missing or zero-occupancy residues), out of which only one structure contains the A-loop devoid of any crystal contacts (PDB ID: 2BAJ). Unfortunately, this structure has been solved in the DFG-out conformation with a bound inhibitor and it cannot be representative of the apo unphosphorylated state. However, the mere lack of information indicates just how flexible the A-loop is and why it is almost impossible to resolve it unless it forms crystal contacts. Our extensive simulations also indicate that the unphosphorylated p38α never explores a conformation that could possibly bind and stabilize ATP - mostly due to A-loop interactions with the DFG-Asp which is crucial for ATP stabilization (Kumar et al., 1995). It is worth noting that related proteins, like ERK2, are able to bind ATP in the unphosphorylated state (Lee et al., 2005; Sours et al., 2008; Canagarajah et al., 1997; Khokhlatchev et al., 1998) most likely due to longer A-loops (Figure 10), which are then easier to stabilize away from the ATP-binding site.

Histograms of the minimal distances between (a) K53-E71, (b) R49-D112, and c) R49-E160 calculated from p38α X-ray structures currently deposited in PDB (March, 2017).

K53-E71 distances were calculated for structures without any mutations, while R49-D112 and R49-E160 for structures without missing R49, D112, and E160 residues. The average values are shown as red dashed lines, while the bin size is 0.5 Å.

https://doi.org/10.7554/eLife.22175.018
Multiple alignments of the αC helix and the activation loop (from the DFG motif to the P+1 loop) sequences of several human kinases.

Selected positively charged residues are colored green, negatively charged ones magenta, while the catalytically important residues in the αC helix are bolded, as well as the phosphorylated residues in the activation loop.

https://doi.org/10.7554/eLife.22175.019

Phosphorylation, however, does change the free-energy surface and enables the A-loop to form salt bridge interactions with the autoinhibitory arginine residues - R67 and R173. The former of the two, as mentioned above, is MAPK-specific and its importance is highlighted further by the R67A mutation that decreases p38α kinase activity (Gum et al., 1998), and the R67I somatic mutation that was found in colon adenocarcinoma (COSMIC: COSM1444062). The newly formed interactions keep the A-loop away from the ATP-binding site which is occluded by the formation of the R49-D112 salt bridge. ATP binding is therefore only possible in less populated conformations of the protein, explaining the experimentally observed low ATP-binding affinity (Tokunaga et al., 2014; Frantz et al., 1998). If one looks solely at X-ray structures where R49 is ~15 Å away from D112 (Figure 9b), the formation of this particular salt bridge might seem like an artifact arising from force field or methodological imperfections. However, R49 is a part of the loop positioned below the very flexible P-loop (whose residues are often unresolved) and we have often observed its interaction with E160 of the ED site in our simulations (both unbiased and biased), regardless of the protein’s phosphorylation state (R49-E160 is ~9 Å in the analyzed X-ray structures - Figure 9c). We believe that the reason for these discrepancies lies in the fact that the loop containing R49 forms crystal contacts in all the 206 deposited X-ray structures which would most likely restrict its mobility and prevent the interactions observed in MD simulations.

While the binding of ATP-Mg2+ helps p38α come closer to the active form, it is insufficient to lead to a complete shift, as proposed by the NMR model (Tokunaga et al., 2014). The full transition, however, does occur in our simulations upon binding of the MK2 docking peptide, which has been previously shown to enhance the catalytic steps (Tokunaga et al., 2014). We hypothesize that the docking peptide prevents the formation of the R49-D112 salt bridge in the phosphorylated protein and brings in the L16 loop. This allows the A-loop to form additional stabilizing contacts, as suggested in the case of ERK2 (Zhou et al., 2006) and causes the observed 20-fold increase in the ATP-binding affinity by the MK2 docking peptide (Tokunaga et al., 2014). The changes that the docking peptide introduces into the protein structure are particularly interesting when it comes to the crystal contact analysis of the dually phosphorylated p38α X-ray structure with the active A-loop conformation. We found an atypically long expression tag of the neighboring molecule bound in the docking site which, together with other crystal contacts, induces the active state stabilized within the crystal. These observations show how important it is to carefully analyze symmetry-related molecules and they call for caution in the interpretation of deposited X-ray structures, as they can be misleading.

Overall, our extensive PT-metaD calculations converge to free-energy surfaces that together with the unbiased MD simulations offer unparalleled insight into the molecular mechanism of p38α canonical activation elucidating a complex conformational behavior that is consistent with experimental observations. Despite the initial discrepancies between previous structural studies, our findings reconcile the deposited models from X-ray crystallography with the NMR data. Considering the importance of p38α for both normal physiology and pathological processes, we hope the knowledge obtained in this study will help to target the protein with more specificity.

Materials and methods

Protein structures and their preparation

The X-ray structure used as the starting conformation for the PT-metaD simulations of inactive unphosphorylated and phosphorylated human p38α is deposited under the Protein Data Bank (PDB) code 3S3I. This structure was also used with the structure of the active phosphorylated mouse p38α (PDB ID: 3PY3) to build the human protein in the latter conformation by using ICM (RRID:SCR_014878) (Abagyan et al., 1994). The coordinates of the ATP molecule and the Mg ions were taken from the active phosphorylated human p38γ (PDB ID: 1CM8), while the coordinates of the MK2 docking peptide (residues 370–393) were taken from the p38α-MK2 complex (PDB ID: 2OKR). The missing atoms in 3S3I and MK2 peptide were built with PDB2PQR (Dolinsky et al., 2004), while the residues T180 and Y182 were phosphorylated using Vienna-PTM server (Margreitter et al., 2013). All proteins consist of the sequence E4-P352 and the standard residue numbering for p38α is applied throughout.

MD simulations

Each simulated system was generated using CHARMM27 force field (MacKerell et al., 1998; Foloppe and MacKerell, Jr., 2000; MacKerell and Banavali, 2000) at pH 7.4. The protonation states of the residues were determined both by H++ server (Anandakrishnan et al., 2012) and PROPKA3.0 (Li et al., 2005) which gave identical results leaving all the residues in their usual charge states, except for His312 which was protonated. Since the pKa2 of phosphoric acid is measured at 6.7 and 7.2 (Peacock and Nickless, 1969; Lide, 2005), the charge of phosphoric residues was set at −1. The systems were solvated with ~20,000 TIP3P water molecules (Jorgensen et al., 1983; Mahoney and Jorgensen, 2000) in a triclinic box with periodic boundary conditions, while Na+ and Cl- ions were added to reach neutrality and the final concentration of 50 mM (the total number of atoms was ~65,000). Energy of the simulated system was initially minimized in two cycles of steepest-descent energy minimization. The initial velocities for the atoms were taken from Maxwell distribution at 100 K, and the system was subsequently heated to 300 K in five steps of 50 K simulated for 200 ps at constant volume each using velocity rescale thermostat (Bussi et al., 2007). In parallel, atomic position restraints for the protein (as well as the ATP molecule and Mg2+ ions if present) were uniformly relaxed (with the restraint spring constant going from 25,000 kJ mol−1 nm−2 to 0 kJ mol−1 nm−2 in steps of 5000 kJ mol−1 nm−2). In the end, the system was equilibrated for additional 200 ps under the final conditions. The magnesium ions were restrained to the ATP molecule, that is its last two phosphate atoms (at 3 and 4 Å) and their adjoining oxygen atoms (two restraints per Mg2+ at 2.4 Å) throughout the equilibration and production runs using a restraint spring constant of 150 kJ mol−1 nm−2. The weak restraints were used to make sure that the ions do not unbind at higher temperature replicas. The production simulations were generated using the GROMACS 4.6.7 biomolecular simulation package (RRID:SCR_014565) (Pronk et al., 2013) with the PLUMED2.1 plugin (Bonomi et al., 2009; Tribello et al., 2014) with a 2-fs integration step. Constant volume conditions were employed throughout with the velocity rescale thermostat (Bussi et al., 2007). Bond lengths were constrained using LINCS (Hess et al., 1997), while van der Waals interactions were treated with a cutoff of 10 Å. Electrostatic interactions were computed using the particle mesh Ewald method (Essmann et al., 1995) with the direct sum cutoff of 10 Å and the Fourier spacing of 1.2 Å.

In the case of unbiased simulations (unphosphorylated p38α with the A-loop in the extended conformation, apo dually phosphorylated p38α with the A-loop in the extended conformation, the dually phosphorylated p38α with bound ATP-Mg2+, and the dually phosphorylated p38α with bound ATP-Mg2+ and MK2 docking peptide), 10 independent runs of each system were equilibrated to a temperature of 380 K by two additional steps of 40 K compared to the aforementioned protocol, while the position restraints for the protein were uniformly relaxed in steps of 2500 kJ mol−1 nm−2. The system was equilibrated for additional 200 ps under constant pressure using Berendsen barostat (Berendsen et al., 1984). The production simulations were generated under the same conditions as above with the addition of constant pressure of one bar using Parrinello-Rahman barostat (Parrinello and Rahman, 1981) and coordinate output of 10 ps. In the case of unphosphorylated system, each independent run was initially 500 ns long. One run that showed a broken β-sheet motif after 500 ns was extended to a total of 2.675 μs. For the apo phosphorylated system, each independent run was 1 μs long. One of the trajectories showed a broken β-sheet motif early on and was extended to 2 μs. Each independent run of the phosphorylated system with bound ATP-Mg2+ was also initially simulated for 1 μs, while one run that exhibited more conformational diversity was extended to 1.55 μs to transition to the global minimum. In the case of the dually phosphorylated system with bound ATP-Mg2+ and MK2 docking peptide, each independent run was 1 μs, amounting to a total sampling time of 10 μs. The magnesium ions were restrained to the ATP molecule throughout the equilibration and production runs in the same way as in the metadynamics runs.

For the simulation of p38α protein crystal, a supercell of the dually phosphorylated mouse p38α (PDB ID: 3PY3) containing four unit cells (2 × 2 × 1), each with four asymmetric units (i.e. 16 protein molecules in total), was built using PyMOL (RRID:SCR_000305) (Schrodinger, 2015). Missing residues of the N- and C-terminal tails were added to each protein using Modeller (RRID:SCR_008395) (Sali and Blundell, 1993). Protein chain whose tail was built buried within the crystal was chosen as the starting structure to avoid steric clashes caused by the lack of periodic boundary conditions during the modelling process. The supercell was rebuilt in PyMOL and used for pKa calculations with PROPKA3.0 (Li et al., 2005) leaving all the residues in the same charged state as in the previous simulations, while among the newly built residues His-10 was protonated. The system was solvated with ~24,000 TIP3P water molecules (Jorgensen et al., 1983; Mahoney and Jorgensen, 2000) in a triclinic box using the dimensions of the original crystal and periodic boundary conditions. Na+ and Cl- ions were added to reach neutrality and the ionic strength of sodium citrate used in the original experiment (with the final concentration of 600 mM). The total number of atoms in the system was ~170,000. The system was minimized and equilibrated as described above to the temperature of 293 K which corresponds to the temperature used for growing crystals. Canonical ensemble at 293 K was maintained for the production run which was 500 ns long. By using a supercell system, crystal contacts were present throughout the simulation, while the interaction of molecules with themselves through periodic boundary conditions was avoided and 8 µs of sampling were obtained for the dually phosphorylated p38α in the crystal environment.

Enhanced sampling

First, a parallel tempering metadynamics (PT-metaD) with 10 replicas per system was performed at increasing temperatures (300, 310, 320, 330, 341, 352, 363, 375, 387, and 400 K) in which only the potential energy was biased to achieve an exchange rate of 30% between the neighboring replicas. The obtained bias was saved and used further in the production runs as a constant bias. The production runs of PT-metaD in the well-tempered ensemble (Bonomi and Parrinello, 2010) were performed for each system consisting of 10 replicas, meaning that a Gaussian potential was deposited in the collective variable space every 2 ps with the height W=W0eV(s,t)/(f1)T, where W0 = 5 kJ mol−1 is the initial height, V(s,t) is the bias potential at time t and CV value s,f = 10 is the biasing factor, and T is the temperature of the replica. The collective variables used in the study are distances in contact map space from the inactive (CV1) and active (CV2) A-loop conformation. CV(R)=1/NγϵΓ(Dγ(R)Dγ(Rref))2, where Dγ(R) is a sigmoidal function that measures the degree of formation of the contact γ in the structure R and is defined as Dγ(R)=Wγ1(rγ/rγ0)n1(rγ/rγ0)m.rγ. is the contact distance in the structure R,rγ0 is the contact distance in either the reference inactive or active conformation depending on which the contact γ is specific for, Wγ is the weight of the contact and is set to one for regular contacts and three for salt bridges, N is a normalization constant, n = 6, and m = 10. The widths of the Gaussians are σ1,2 = 0.05. The set of contacts Γ defining the contact maps was determined as follows: minimized structures of inactive (PDB ID: 3S3I) and active (PDB ID: 3PY3) conformations were used to determine all the pairs of Cα, Cβ or backbone O atoms that fall within a 5 Å distance cutoff. To relax and lose possible crystal contacts, 10 independent equilibration runs were generated for both the inactive and the active phosphorylated systems. Only those atomic pairs that appear in five or more equilibrated structures and that specifically discriminate between the inactive and active conformations were kept, that is, those that appear either in one or the other, thereby discarding the pairs common to both conformations. A set of 97 regular contacts was obtained this way to which seven salt-bridge contacts involving residues of the A-loop that were specific to the active conformation were added and weighed three times more than regular contacts.

Convergence of metaD simulations

Initially, FE surfaces were produced with a bias corresponding to an exchange rate of 30% between neighboring replicas and were subsequently refined by lowering the bias on the potential energy to achieve an exchange rate of 10% and run for additional 150 ns to convergence. The simulations were considered to be converged once the free energy in the bidimensional projections and in the monodimensional projections does not change more than 2 kcal mol−1 in the last 100 ns. This criterion resulted in a 2.425-μs-long simulation for the unphosphorylated, 2.485-μs-long simulation for the phosphorylated system, 1.035-μs-long simulation for the phosphorylated system with bound ATP, and 0.986-μs-long simulation for the phosphorylated system with bound ATP and MK2 docking peptide.

Analysis of structures populating the free-energy minima

Free-energy surfaces (FES) for the two CVs were calculated by integrating the bias deposited during the simulations. For each system, several basins from the FES plot were chosen for further analysis. All the structures that fell within each basin were clustered using the gromos algorithm (Daura et al., 1999) and the g_cluster routine (GROMACS) by using the RMSD of Cα atoms as the distance between the structures and the cutoff value of 2 Å. The central structure of the most populated cluster of each basin (i.e. the structure with the smallest distance to all the other members of the cluster) was chosen as the representative of the basin. Furthermore, central structures of top five most populated clusters (which on average contain 75% of all structures) were inspected to make sure that the main structural features in terms of P-loop and A-loop positions were captured by these clusters. To avoid highlighting interactions that occur only sporadically in the simulations, only those salt-bridge interactions that are maintained in the majority of structures of the most populated cluster, as well as most of structures of the selected basin, were highlighted in the tables throughout the paper.

Analysis of X-ray structures and their symmetry-related neighbors

Symmetry-related molecules of mouse apo p38α-pTpY (PDB ID: 3PY3) that are closer than 4 Å to the original structure were generated by PyMOL (Schrodinger, 2015). The same distance criterion was used to determine which pairs of residues of neighboring molecules come into contact.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
    CRC Handbook of Physics and Chemistry, Internet Version 2005
    1. DR Lide
    (2005)
    Boca Raton, Florida: CRC Press.
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
    The PyMOL Molecular Graphics System
    1. LLC Schrodinger
    (2015)
    The PyMOL Molecular Graphics System, Version 1.8.
  49. 49
  50. 50
  51. 51
  52. 52
    New advances in metadynamics
    1. L Sutto
    2. S Marsili
    3. FL Gervasio
    (2012)
    Wiley Interdisciplinary Reviews: Computational Molecular Science 2:771–779.
    https://doi.org/10.1002/wcms.1103
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64

Decision letter

  1. Yibing Shan
    Reviewing Editor; DE Shaw Research, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Molecular mechanism of canonical activation of p38α MAP kinase" for consideration by eLife. Your article has been favorably evaluated by John Kuriyan (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

This manuscript reports simulation studies aiming to elucidate the activation mechanism of p38α kinase. Using unbiased molecular dynamics simulations and calculation of free energy surfaces of the system, the authors attempt to reconcile the apparent contradictions between the X-ray structures of p38α (showing large structural changes associated with kinase dual-phosphorylation) and the NMR findings (the phosphorylation does not lead to significant chemical shift perturbation).

Essential revisions:

1) The free energy surface (FES) of ATP-bound phosphorylated system shows that the active state is accessible now, but unbiased MD simulations showed that only with substrate binding the system is stable at the active basin. To complete the paper's argument, FES of phosphorylated p38α with ATP and substrate bound should be calculated besides the unbiased simulations at 380 K. Also, showing unbiased MD on the phosphorylated and ATP-bound system without substrate and mapping the trajectory to the two CVs analogous to Figure 4 would be helpful to this argument.

2) p38α should be simulated with its crystal contacts in unbiased simulations to support the notion that crystal contacts maintain apo p38α in its otherwise unstable active conformation in crystal environment.

3) The paper is descriptive of selected conformations representative of conformational ensembles at various free energy basins; but, it is not clear how the discussed conformations are selected to represent the ensemble of each basin. The selections should be justified based on characterization of the ensemble of the entire free-energy basin; this justification is missing as is.

Other important points:

1) For phosphorylated p38α, a histogram/time-series should be shown for R49-D112 salt-bridge (ascribed to occlude ATP binding). The regions in the free energy surface allowing ATP binding should also be shown. Reviewers are concerned that the free energy surfaces of the apo unphosphorylated system are very similar to the apo phosphorylated, but the latter system binds ATP much stronger than the former. (In Tokunaga et al., 2014, the ATP affinity is 430 μM, while Frantz et al., 1998 only reported the Km[ATP] to be 25 μM. Unphosphorylated systems binds ATP much weaker at millimolar affinity.) How is the μM ATP binding reconciled with the R49-D112 salt-bridge? Another concern is that, by the energy surface of the apo phosphorylated system, the ATP-bound X-ray structure is as much as 16~18 kcal/mol higher in free energy than the inactive conformation, which seems unlikely.

2) The manuscript reports certain structural features differ considerably from any of the known crystallographic structures (e.g. 3S3I, 3PY3, 2OKR and 1CM8), which calls to question the structural details extracted from the PT-metaD trajectories (e.g. the activation-loop helix (purple in Figure 3B, and the R49-D112 salt bridge, which requires D113 to move 18 Å from its positions in the X-ray structures). The reliability of these simulations generated conformations need be calibrated in Discussion.

3) Finally, how the simulation results explain the seemingly contradictory X-ray and NMR data needs to be further discussed in more explicit terms. The contradictory is understood to be the lack of chemical shift perturbation (CSP) between unphosphorylated and phosphorylated yet the x-ray structures are the closed and open alternative structures of the A-loop. The NMR study also shows a large CSP for ATP binding to phosphorylated p38α. The minima in the FE surfaces in Figure 2 differ between the 3 states (unphosphorylated, phosphorylated and ATP bound phosphorylated) yet have overlapping regions. The overlap could explain the lack of CSP for unphosphorylated and phosphorylated (as stated in the subsection “Dually phosphorylated p38α”), but how the large CSP for bound ATP is rationalized by the FE surface, since, as stated in the manuscript, "the representative structure of the global minimum [for ATP-bound] is still not too different from the dominant one in the apo phosphorylated state." This should be explained in more explicit terms.

https://doi.org/10.7554/eLife.22175.023

Author response

Essential revisions:

1) The free energy surface (FES) of ATP-bound phosphorylated system shows that the active state is accessible now, but unbiased MD simulations showed that only with substrate binding the system is stable at the active basin. To complete the paper's argument, FES of phosphorylated p38α with ATP and substrate bound should be calculated besides the unbiased simulations at 380 K. Also, showing unbiased MD on the phosphorylated and ATP-bound system without substrate and mapping the trajectory to the two CVs analogous to Figure 4 would be helpful to this argument.

Following reviewers’ suggestions, we have now complemented our manuscript with the FES of the phosphorylated p38α with bound ATP and the MK2 docking peptide (with overall sampling of >9 µs). We have changed Figure 2 and its supplements accordingly, as well as the manuscript text to reflect the addition of new data. We have also performed unbiased MD simulations of the phosphorylated p38α with and without bound ATP which resulted in two new figures (Figure 4 and Figure 5) and >20 µs of sampling of these systems. What is more, we have extended our unbiased simulations of the phosphorylated system with bound ATP and the MK2 docking peptide (reaching 10 µs of sampling) whose results are shown in an updated version of Figure 4 (Figure 6 in the revised version). The additional simulations give further support to the claims we have made in the original submission.

2) p38α should be simulated with its crystal contacts in unbiased simulations to support the notion that crystal contacts maintain apo p38α in its otherwise unstable active conformation in crystal environment.

We thank the reviewers and the editor for the useful advice. The results of the new simulations strengthen the overall conclusions. We have simulated a crystal supercell of the dually phosphorylated mouse p38α (PDB ID: 3PY3) containing 4 unit cells (2x2x1), each with 4 asymmetric units (16 copies of the protein in total). By using a supercell system, we have been able to maintain crystal contacts, avoid the interaction of molecules with themselves through periodic boundary conditions, and achieve 8 µs of sampling (500 ns of the whole supercell x 16 copies of the protein). We have described the modelling process and the simulation setup in the Methods section (subsection “MD simulations”, last two paragraphs), added the results as Figure 8, and made appropriate changes in the text referring to new data (subsection “The active conformation”, fourth paragraph).

As shown by Figure 8, the A-loop is very stable in the conformation captured by the X-ray model – the β-sheet motif is unperturbed during the simulation, as well as the contacts made by pT180 with the three conserved arginines (R67, R70, and R173). However, pY182 is more mobile and at times breaks away from H228 and the conserved arginines (R186 and R189). This is consistent with the electron density which is not so well defined for pY182 as it is for the rest of the A-loop (see Author response image 1). Finally, we show an excellent qualitative agreement with the B-factors reported for the X-ray structure (3PY3), while (as expected) the most mobile regions in our simulations (mainly, the N- and C-terminal tails) are missing from the deposited electron density. While comparing the crystallographic and MD-derived B-factors it is also important to keep in mind that: 1) The data were collected at 100 K. Since MD force-fields and algorithms have not been developed for such low temperatures, we have performed our crystal simulations at 293 K which proportionally increases the overall atomic motions; 2) The experimental data collected for the crystals in general are time and ensemble averages of 1023 molecules which are in the end represented by a single model, while its B-factors then reflect both the static and the dynamic disorder of the crystal. Indeed, two independent studies had shown that the refinement process underestimates the disorder present in the crystal thereby producing lower B-factors (Janowski et al., 2013, JACS, doi: 10.1021/ja401382y; Kuzmanic et al., 2014, Nat. Commun., doi: 10.1038/ncomms4220).

3) The paper is descriptive of selected conformations representative of conformational ensembles at various free energy basins; but, it is not clear how the discussed conformations are selected to represent the ensemble of each basin. The selections should be justified based on characterization of the ensemble of the entire free-energy basin; this justification is missing as is.

We analyzed each selected free energy basin through a clustering method (as described in the Methods section) and we chose the central structures of the most populated clusters as the representative ones we had shown and discussed in the paper. Furthermore, we checked the central structures of top 5 most populated clusters (which on average contain 75% of all structures) to make sure that the main structural features in terms of P-loop and A-loop positions were captured by these clusters. To avoid highlighting interactions that occur only sporadically in the simulations, we also added the tables (Table 1 and Supplementary file 1) with salt-bridge interactions that are maintained in the majority of structures of the most populated cluster, as well as most of structures of the selected basin. In order to avoid confusion among the readers, we now state clearly in the main text how the structures were chosen (Results, first paragraph and subsection “Enhanced sampling”). We have also modified the aforementioned tables that now include a measure of how often the interaction occurs in the most populated cluster of each bin.

Other important points:

1) For phosphorylated p38α, a histogram/time-series should be shown for R49-D112 salt-bridge (ascribed to occlude ATP binding). The regions in the free energy surface allowing ATP binding should also be shown. Reviewers are concerned that the free energy surfaces of the apo unphosphorylated system are very similar to the apo phosphorylated, but the latter system binds ATP much stronger than the former. (In Tokunaga et al., 2014, the ATP affinity is 430 μM, while Frantz et al., 1998 only reported the Km[ATP] to be 25 μM. Unphosphorylated systems binds ATP much weaker at millimolar affinity.) How is the μM ATP binding reconciled with the R49-D112 salt-bridge? Another concern is that, by the energy surface of the apo phosphorylated system, the ATP-bound X-ray structure is as much as 16~18 kcal/mol higher in free energy than the inactive conformation, which seems unlikely.

Following reviewers’ suggestions, we have now included the histogram for the R49-D112 salt bridge as Figure 2—figure supplement 4, which shows that almost all the structures in the most populated cluster of the global minimum have this interaction formed (>16000 structures). We have also highlighted the minimum in the free energy surface of the phosphorylated p38α which is represented by structures that we believe can bind ATP (and which are also explored by the unbiased simulations).

Furthermore, as reviewers have noted, the free energy surfaces of apo unphosphorylated and apo phosphorylated p38α are indeed similar, but they are not the same, nor are their minima populated by the same conformations. In general, free energy surface of a protein reflects its ability to populate certain states defined by collective variables. This means that at a given temperature all the conformations of the system that we have highlighted in the paper will be accessible; however, their population will be different.

Based on our analyses, apo unphosphorylated p38α does not seem to explore conformations that can bind ATP – either the highly flexible P-loop is collapsed into the ATP-binding site, or the R173 of the A-loop forms interactions that are typically viewed as autoinhibitory (e.g., forming a salt bridge with D168 of the DFG motif which stabilizes ATP upon binding). Upon the A-loop phosphorylation, two residues change their properties significantly which allows them to stabilize interactions which are inaccessible to the unphosphorylated system. However, the global minimum of the phosphorylated system is dominated by the conformation with the formed R49-D112 salt bridge which would prevent ATP binding. Nevertheless, this does not mean that the phosphorylated p38α cannot bind ATP, seeing how a higher-energy minimum is populated by conformations that can bind ATP (which are also explored by the unbiased simulations). Basically, once p38α is phosphorylated, it exists as a mixture of all the described conformations that are present at different concentrations, therefore determining the macroscopic property of ATP binding to the phosphorylated system. This agrees with the experimental observations mentioned by the reviewers, as well as the observation that the binding of the MK2 docking peptide has a profound effect on ATP binding, i.e. it causes a 20-fold increase in binding which most likely occurs due to the breaking of the salt bridge and the overall stabilization of the αD helix. However, our results do not exclude the possibility of induced fit mechanism in which the addition of ATP to the phosphorylated system would shift the populations towards structures that are more likely to bind ATP.

Finally, we have shown with our metadynamics simulations (and confirmed additionally with unbiased MD simulations) that the apo phosphorylated structure captured by X-ray crystallography is 16-18 kcal/mol higher in energy and therefore rather inaccessible to the phosphorylated system in the absence of ATP. However, the presence of bound ATP in the system causes a shift towards this conformation which is then at 2-4 kcal/mol compared to the global minimum, not at 16-18 kcal/mol. This observation agrees with the NMR data which show large chemical shift perturbations only upon ATP binding, not the dual phosphorylation.

We have amended our manuscript throughout to discuss the points above more clearly (e.g., subsections “Unphosphorylated p38α” and “Dually phosphorylated p38α with bound ATP-Mg2+”, and Discussion).

2) The manuscript reports certain structural features differ considerably from any of the known crystallographic structures (e.g. 3S3I, 3PY3, 2OKR and 1CM8), which calls to question the structural details extracted from the PT-metaD trajectories (e.g. the activation-loop helix (purple in Figure 3B, and the R49-D112 salt bridge, which requires D113 to move 18 Å from its positions in the X-ray structures). The reliability of these simulations generated conformations need be calibrated in Discussion.

We agree with the reviewers that there is always a possibility that certain structural features arise from either force field or methodological imperfections, we have thus modified the Discussion which now mentions these possible pitfalls (Discussion, first and second paragraphs). However, we also find that crystal contacts introduce artifacts that can lead to misinterpretation of the structural features of p38α (as we have already described above for 3PY3). When it comes to the position of the A-loop in the unphosphorylated state, out of 206 currently deposited X-ray structures of human p38α (with 224 individual chains), only 47 molecules have complete A-loops (i.e., without any missing or zero occupancy residues). From this reduced set, a single deposited structure has the A-loop in a conformation in which the loop does not form any crystal contacts (PDB ID: 2BAJ). Unfortunately, this structure has been solved in the DFG-out conformation with a bound inhibitor and it cannot be representative of the apo unphosphorylated state. However, the mere lack of information indicates just how flexible the A-loop is and why it is almost impossible to resolve it unless it forms crystal contacts. Thus, it is not surprising that the A-loop conformations we observe do not match the ones present in X-ray structures. Furthermore, we have shown that extensive crystal contacts are in fact responsible for the conformation observed in 3PY3 (the only X-ray structure of the dually phosphorylated p38α), while we suspect that the same applies for 1CM8 (dually phosphorylated p38γ with bound ATP). While this structure does not have anything bound in the docking site, its A-loop forms crystal contacts with the N-terminal loops of a symmetry related molecule (see Author response image 2). Unfortunately, the electron density of this molecule is unavailable, so we are unable to inspect the quality of the modelled A-loop and its crystal contacts. We have also included this analysis in the manuscript (Discussion).

Author response image 2
Contact maps for interactions formed between the residues of the reported X-ray structure (PDB ID: 1CM8, p38γ-pTpY) (x-axis) and its symmetry related neighbors (y-axis).

The residues are considered to be in contact if the minimal distance between any of their atoms is under 4 Å. The contacts formed by the A-loop are indicated by the solid ellipse. Secondary structure elements are also shown with α-helices in purple, β-sheets in green, and loops in grey. Missing residues are shown as light grey rectangles.

https://doi.org/10.7554/eLife.22175.022

As for the position of R49 and its interactions, the deposited X-ray structures are again quite misleading. R49 is a part of the loop positioned below the very flexible P-loop (whose residues are often unresolved) and is on average 9 Å away from the C-lobe residues in the deposited X-ray structures. However, in our simulations (both unbiased and biased), we have often observed R49 interacting with E160 of the ED site, regardless of the protein’s phosphorylation state. We believe that the reason for this discrepancy (as well as the formation of the R49-D112 salt bridge) lies in the fact that R49 loop forms crystal contacts in all 206 deposited X-ray structures which would most likely restrict its mobility and prevent the interactions observed in MD simulations. We have included these findings in the revised version of the manuscript (Discussion, second paragraph).

3) Finally, how the simulation results explain the seemingly contradictory X-ray and NMR data needs to be further discussed in more explicit terms. The contradictory is understood to be the lack of chemical shift perturbation (CSP) between unphosphorylated and phosphorylated yet the x-ray structures are the closed and open alternative structures of the A-loop. The NMR study also shows a large CSP for ATP binding to phosphorylated p38α. The minima in the FE surfaces in Figure 2 differ between the 3 states (unphosphorylated, phosphorylated and ATP bound phosphorylated) yet have overlapping regions. The overlap could explain the lack of CSP for unphosphorylated and phosphorylated (as stated in the subsection “Dually phosphorylated p38α”), but how the large CSP for bound ATP is rationalized by the FE surface, since, as stated in the manuscript, "the representative structure of the global minimum [for ATP-bound] is still not too different from the dominant one in the apo phosphorylated state." This should be explained in more explicit terms.

As also discussed above, the free energy surfaces do not only reflect which areas have been explored by p38α in MD simulations, but also how often. This means that even though the free energy surface of the phosphorylated system with bound ATP overlaps with the apo one, these systems predominantly explore different conformations. Upon ATP binding the A-loop can finally explore regions that are almost inaccessible in the absence of ATP and it also forms more stable contacts with the αC helix in the global minimum. These changes, along with the fact that ATP is in direct contact with the binding site residues, would explain the differences observed by NMR. We have modified the manuscript text to clarify these points (e.g., subsection “The active conformation”).

https://doi.org/10.7554/eLife.22175.024

Article and author information

Author details

  1. Antonija Kuzmanic

    1. Institute for Research in Biomedicine (IRB Barcelona), The Barcelona Institute of Science and Technology, Barcelona, Spain
    Contribution
    AK, Conceptualization, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0003-2815-5605
  2. Ludovico Sutto

    1. Department of Chemistry, University College London, London, United Kingdom
    Contribution
    LS, Conceptualization, Investigation, Methodology
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-4084-8562
  3. Giorgio Saladino

    1. Department of Chemistry, University College London, London, United Kingdom
    Contribution
    GS, Formal analysis, Investigation, Methodology
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-3234-5762
  4. Angel R Nebreda

    1. Institute for Research in Biomedicine (IRB Barcelona), The Barcelona Institute of Science and Technology, Barcelona, Spain
    2. Catalan Institution for Research and Advanced Studies (ICREA), Barcelona, Spain
    Contribution
    ARN, Conceptualization, Supervision, Funding acquisition, Project administration, Writing—review and editing
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-7631-4060
  5. Francesco L Gervasio

    1. Department of Chemistry, University College London, London, United Kingdom
    Contribution
    FLG, Conceptualization, Resources, Supervision, Investigation, Methodology, Writing—review and editing
    For correspondence
    1. f.l.gervasio@ucl.ac.uk
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0003-4831-5039
  6. Modesto Orozco

    1. Institute for Research in Biomedicine (IRB Barcelona), The Barcelona Institute of Science and Technology, Barcelona, Spain
    2. Joint BSC-CRG-IRB Program in Computational Biology, Barcelona, Spain
    3. Department of Biochemistry, University of Barcelona, Barcelona, Spain
    Contribution
    MO, Conceptualization, Resources, Supervision, Funding acquisition, Project administration, Writing—review and editing
    For correspondence
    1. modesto.orozco@irbbarcelona.org
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-8608-3278

Funding

Seventh Framework Programme (FP7/Marie Curie Actions/COFUND, no. 600404)

  • Antonija Kuzmanic

Engineering and Physical Sciences Research Council (EP/M013898/1)

  • Giorgio Saladino
  • Francesco L Gervasio

European Research Council (Project ID 294665)

  • Angel R Nebreda

Horizon 2020 (BioExcel (ID: 676559))

  • Modesto Orozco

Horizon 2020 (ELIXIR-Excellerate (ID: 675728))

  • Modesto Orozco

Ministerio de Economía y Competitividad (BIO2015/64802-R)

  • Modesto Orozco

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Dr. R Soliva and Dr. Lucia Diaz Bueno for many helpful discussions. We acknowledge PRACE for awarding us access to resource MareNostrum based in Spain at Barcelona Supercomputing Center. AK was partly supported by the European Union Seventh Framework Programme (FP7/Marie Curie Actions/COFUND, no. 600404). FLG and GS acknowledge financial support from EPSRC (grant EP/M013898/1). ARN was supported by a grant from the European Commission (ERC, project ID: 294665). MO is an ICREA research fellow. MO acknowledges the support from two grants from the European Commission H2020 program (BioExcel project ID: 676559 and ELIXIR-Excellerate project ID: 675728) and the Spanish MINECO (BIO2015/64802-R). We gratefully acknowledge institutional funding from the Spanish Ministry of Economy, Industry and Competitiveness (MINECO) through the Centers of Excellence Severo Ochoa award, and from the CERCA Programme of the Catalan Government.

Reviewing Editor

  1. Yibing Shan, Reviewing Editor, DE Shaw Research, United States

Publication history

  1. Received: October 10, 2016
  2. Accepted: April 6, 2017
  3. Version of Record published: April 26, 2017 (version 1)

Copyright

© 2017, Kuzmanic et al

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 807
    Page views
  • 185
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Scopus, Crossref.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Genomics and Evolutionary Biology
    2. Microbiology and Infectious Disease
    Kevin S Bonham et al.
    Research Article