p38α Mitogen-Activated Protein Kinase Is a Druggable Target in Pancreatic Adenocarcinoma

p38 mitogen-activated protein kinases are signaling molecules with major involvement in cancer. A detailed mechanistic understanding of how p38 MAPK family members function is urgently warranted for cancer targeted therapy. The conformational dynamics of the most common member of p38 MAPK family, p38α, are crucial for its function but poorly understood. Here we found that, unlike in other cancer types, p38α is significantly activated in pancreatic adenocarcinoma samples, suggesting its potential for anti-pancreatic cancer therapy. Using a state of the art supercomputer, Anton, long-timescale (39 μs) unbiased molecular dynamics simulations of p38α show that apo p38α has high structural flexibility in six regions, and reveal potential catalysis mechanism involving a “butterfly” motion. Moreover, in vitro studies show the low-selectivity of the current p38α inhibitors in both human and mouse pancreatic cancer cell lines, while computational solvent mapping identified 17 novel pockets for drug design. Taken together, our study reveals the conformational dynamics and potentially druggable pockets of p38α, which may potentiate p38α-targeting drug development and benefit pancreatic cancer patients.

INTRODUCTION p38 mitogen-activated protein kinases (MAPKs) play critical roles in cellular responses, proliferation, survival, cell cycle, and migration in cancer. p38 MAPK family includes p38α (MAPK14), p38β (MAPK11), p38γ (MAPK12), and p38δ (MAPK13). The four p38 MAPK family members have different tissue expression patterns, with p38α being ubiquitously expressed at significant levels in most cell types. The p38 MAPKs function in a cell context-dependent manner (1)(2)(3)(4). However, a detailed mechanistic understanding of how p38 MAPK family members function is still not well-understood. A major challenge will be to determine when and how to specifically target p38 MAPK for disease therapy.
To gain insights into its molecular mechanisms and design potential therapeutics, the structure of MAP kinase p38α has been extensively studied over the last two decades (5)(6)(7)(8). X-ray crystallography and nuclear magnetic resonance (NMR) showed that p38α consists of two domains, a 135-residue N-terminal domain mainly composed of β-sheets and a 225-residue C-terminal domain mainly composed of αhelices, in between of which lies the catalytic site, i.e., the ATP-binding pocket (9,10). Despite its seemingly rigid crystal structures, p38α is highly dynamic, which is supported by various evidence. Firstly, in the majority (∼78%) of the p38α crystal structures, e.g., 3L8X (11), 3OCG (12), and 2NPQ (13), the glycine-rich loop and/or the activation loop is invisible in the electron density map, indicating the structural flexibility of these loops. Secondly, hydrogen-exchange mass spectrometry (HX/MS) study showed that phosphorylation of the activation loop induced conformational changes in p38α (14); While X-ray crystallographic study showed that phosphorylation brings the N-terminal domain and C-terminal domain closer (15). Thirdly, NMR experiments of p38α in apo and ligand-bound forms suggested that the ATP-binding pocket is highly flexible even after ligand binding (16,17). These studies reveal the important roles of conformational dynamics in the activation and catalysis of p38α.
Pharmacologically, p38α is considered as a potential target for various diseases such as inflammatory diseases and cancer. The potent small-molecule inhibitor SB203580 competitively binds the ATP-binding pocket and has been widely used to study p38 MAPK functions (5). Nevertheless, given the similarity of the ATP-binding site in different kinases, SB203580 also targets non-p38 protein kinases, usually at higher concentrations (18). Over the years, many p38α inhibitors have been developed, for example BIRB796, which leads to a conformational reorganization that prevents ATP binding and activation (19). The C-terminal domain has also been predicted to have the flexibility for potentially binding differently shaped compounds (20). For p38α inhibitor development, allosteric small-molecule inhibitors that target other regions of the kinase are warranted and they might reduce the off-target effects in drug applications.
Clinically, p38 MAPK inhibitors show significant effects in pre-clinical animal models but repeatedly failed in clinical trials (21). In cancer study, p38α and p38β increase cell proliferation and invasion of colon cancer, follicular lymphoma, ovarian cancer, and more recently, pancreatic cancer (22)(23)(24)(25). It is indicated that targeting p38 MAPKs, especially p38α, should be explored for cancer therapy. However, several phase I clinical studies with p38 MAPK small-molecule direct or indirect inhibitors show hepatic, neurological, gastrointestinal, and cardiovascular toxicities, indicating that the commonly studied inhibitors are not highly selective (26). Currently, one of the major challenges in pancreatic cancer drug development is to overcome drug off-target effects and drug resistance (26,27). It is rational to speculate that next-generation highly selective p38 MAPK inhibitors may exhibit less adverse effects. However, it still requires further investigation in cancer patients.
In this work, by screening pancreatic adenocarcinoma (PDAC) patient samples (n = 40) and TCGA database, we demonstrated that p38 MAPKs, especially p38α are highly expressed and activated. p38α blockades show significant antitumor effect in PDAC cells, but are not selective enough. To pave the way for highly selective inhibitor development, conformational dynamics of p38α was examined. Anton supercomputer long-timescale (39 µs) MD simulations using both AMBER and OPLS force fields show the p38α flexibility in six regions. Mechanistically, p38α MD simulations reveal a "butterfly" motion that might be important for p38α catalytic function. In addition, computational solvent mapping reveals 17 novel pockets that are potentially druggable for cancer therapy. To our knowledge, this is the first comprehensive study of both the conformational dynamics and potentially druggable pockets in p38α. This study provides insights into understanding the molecular mechanism of p38α function and into developing potential drugs with high specificity and selectivity against PDAC.

RESULTS p38 MAPKs Expression Correlates With Poor Prognosis in PDAC Patients
To investigate the role of p38 MAPKs in cancer, we screened a panel of human tumor tissues in The Cancer Genome Atlas (TCGA) that spontaneously express p38 MAPKs. We have found that the majority of the tumor tissues express similar level of p38 MAPKs to its adjacent healthy tissues (black labeled, Figure 1A, Figures S1A-C), while uterine carcinosarcoma (UCS), uterine corpus endometrial carcinoma (UCEC), and chromophobe renal cell carcinoma (KICH) show decreased p38α (MAPK14), p38β (MAPK11), and p38γ (MAPK12) expression compared with control pancreas (green labeled, Figure 1A, Figures S1A,B). Surprisingly, between more than 30 types of cancer, only pancreatic adenocarcinoma (PAAD or PDAC) shows a significant increase of p38α, p38β, and p38γ compared with healthy control pancreas (red labeled, Figure 1A, Figures S1A,B), suggesting that p38 MAPKs may be involved in PDAC development. To further validate these findings, we detected the mRNA level for all four members of the p38 MAPK family. Indeed, p38α, p38β, and p38γ are all increased markedly, while p38δ (MAPK13) shows a trend of increase ( Figure 1B). p38α shows a significant dominance among these genes ( Figure 1C), indicating its importance among four members. We next investigated the relationship of p38 with PDAC clinical outcomes. Surprisingly, p38α showed a trend of increasing with advanced stage, and strongly correlated with poor overall survival in PDAC patients (Figures 1D,E). These results suggest that p38α plays a role in PDAC progression and correlates with poor prognosis in PDAC patients.

p38α Expression Correlates With Adipose Markers in PDAC Tissues
It was reported that proliferating cancer cells may take up exogenous lipids and activating endogenous lipid biosynthesis (28), and tumor implanted in adipose environment show significant lipid metabolic reprogramming (29). Considering that PDAC is one of the tumors that adjacent to the adipose environment, we tested the correlation of p38 MAPKs and lipid droplet marker perilipin (PLIN) family in the PDAC database. Surprisingly, p38α strongly correlated with PLIN 2 and 3, two small lipid droplets markers, but not PLIN1, 4, and 5, which   are the big lipid droplets marker ( Figure 1F). It was plausible that small lipid droplets exist inside pancreatic tumor cells or in mobilized adipocytes in the tumor microenvironment. These findings support that p38α correlates with adipose-rich PDAC and may be involved in cancer lipid metabolism.

p38α Is Activated in PDAC Patient Samples
To validate our findings, we measured p38α levels in human PDAC patients. In this study, the patient PDAC sample and adjacent pancreas were compared. In accordance with previous report (30), PDAC tissues show significant infiltration of inflammatory cells, stromal cellular components, and activated fibroblasts (Figure 2A). Next, we detected p38α in these samples. As expected, p38α is highly expressed in human PDAC tissues (Figures 2B,D; Figure S2). Consistent with this result, western blot and immunohistochemistry staining show highly activated phospho-p38α in human PDAC tissues (Figures 2B,D). Furthermore, pathological analysis shows that both p38α and phospho-p38α are mainly located in epithelial cancer cells ( Figure 2C). Western blot of various cell lines demonstrates that PDAC cancer cells and endothelial primary cells expressed high levels of p38α protein, whereas human THP-1 macrophage-like cell line, and human stromal fibroblasts lacked p38α expression ( Figure 2E). These findings further support our notion that p38α is activated in PDAC cells. These pilot clinical findings validate our TCGA data.

p38α Blockade Enhances Apoptosis of Human and Mouse PDAC Cells
To gain further mechanistic insights of p38α blockade on cancer, human pancreatic adenocarcinoma cells MiaPaca-2, Panc-1, and mouse PDAC cell line Panc02 were projected for p38α blockade. The canonical inhibitor SB203580, which can block both p38α and p38β (31), and ralimetinib, a small molecule inhibitor specifically for p38α under phase 2 clinical trial (32), were tested. As expected, both drugs significantly inhibit PDAC cell viability, with the IC50 around 56.89µM ( Figure 2F). These findings suggest the potential of p38α as a therapeutic target for PDAC.

p38α Is Highly Dynamic Locally and Globally
Despite its potential for PDAC therapy, our understanding of the structural and functional basis of p38α is still poor. To investigate the structural dynamics of p38α, we applied the state of the art supercomputer, Anton for simulation of our potential targets (33). Three simulation runs of apo p38α and one simulation run of ligand-bound p38α ( Table 1) for the AMBERff99SB-ILDN and OPLS-AA/L force fields (abbreviated as AMBER and OPLS) were performed ( Table 2) and the convergence was evaluated ( Figure S3). We plot the root mean square deviations (RMSDs) of the protein backbone from their crystallographic positions for all simulation runs. In AMBER simulations, RMSDs of apo p38α range from 2 to 5 Å ( Figure 3A, left panel, red, orange, and green lines); in OPLS simulations, RMSDs of the apo p38α run go up to 8 Å ( Figure 3A, right panel, red line). In contrast, ligand-bound p38α is below 4 Å in both AMBER and OPLS simulations ( Figure 3A, blue lines). The RMSD discrepancy between AMBER and OPLS simulations is possibly due to their different preference for secondary structure propensities. The RMSD discrepancy between apo and ligand-bound p38α suggests that the conformational flexibility of p38α is dampened upon ligand binding. It should be noted that as the N terminal and C terminal ends of p38α can be quite floppy ( Figure S4), these ends were truncated in MD simulations. These results show that apo p38α exhibits higher structural flexibility than ligandbound p38α.
To show global structural excursions of apo p38α, we visually inspected trajectories of simulation runs and found several representative conformations ( Figure 3B, pink ribbons) that were aligned to the crystal structure ( Figure 3B, light blue ribbons). In AMBER apo3, we see a spike in the RMSD plot at ∼1,790 ns ( Figure 3A, left panel, green line), which attributes to the loss of helicity in αD helix and "closing" motion of the glycine-rich loop ( Figure 3B, AMBER apo3, magenta circles; Figure S5). In OPLS apo1, RMSD rises rapidly to ∼8 Å ( Figure 3A, right panel, red line), which can be mainly attributable to the "closing" motion of the glycine-rich loop that is accompanied by large conformational changes in both the N-terminal and C-terminal domains ( Figure 3B, OPLS apo1 magenta circle). In OPLS apo2, RMSD climbs up to ∼4.5 Å at ∼1,000 ns and descends to 3.5 Å at ∼2,000 ns ( Figure 3A, right panel, orange line), which identify the "closing" and "opening" of the glycine-rich loop ( Figure 3B, OPLS apo2, magenta circle). In OPLS apo3, RMSD rises steadily to ∼4 Å ( Figure 3A, right panel, green line), which classifies the "closing" of the glycinerich loop ( Figure 3B, OPLS apo3, magenta circle). Among these structural excursions, OPLS apo1 seems to have the highest RMSD at 7.0 Å ( Figure 3B, OPLS apo1, magenta circle), probably due to the global conformational change of both the N-terminal and C-terminal domains. It is also intriguing to see that in both OPLS apo1 and OPLS apo2, the glycine-rich loop is deformed significantly compared with the β-hairpin structure in the crystal structure ( Figure 3B, OPLS apo1 and OPLS apo2).
Furthermore, to understand the local structural dynamics of p38α, residue-based root mean square fluctuations (RMSFs) were calculated and mapped onto the p38α structure ( Figure 3C). In general, we found that local structural dynamics is highly dependent on the secondary structure of the local region. In The atom names are consistent with the crystal structure (PDB code: 1A9U).  . Taken together, these data suggested p38α is highly dynamic in the entire kinase structure and in the specific secondary structures of local regions.

"Butterfly" Motions Potentially Contribute to p38α Enzymatic Catalysis
To gain further insights on the catalysis mechanism of p38α, we performed a principal component analysis (PCA) on p38α trajectories. For both AMBER and OPLS simulations, the trajectories of three apo p38α runs and one ligand-bound p38α run were combined to build a structure ensemble for PCA. The two most dominant PCs, PC1 and PC2, account for a large fraction of the correlated motions (S7A). Both AMBER and OPLS force fields produce considerable structural excursions (as big as 100 Å) from the crystal structures (Figures 4A,B). AMBER simulations show a more restricted structural excursion pattern ( Figure 4A) while OPLS simulations have a more expanded pattern ( Figure 4B). These results are in accordance with our previous findings ( Figure 3A). Structural morphing between the two extreme structures along the PC1 and PC2 axis was performed (Movies S1-S4). For PC1 and PC2, we overlaid the two extreme structures onto each other (Figures 4C-F). Note that the pink and blue structures represent the + and -axis. In AMBER simulations, PC1 represents the "butterfly" motion between the N-terminal and C-terminal domains ( Figure 4C; Movie S1), where the Nterminal domain and the C-terminal domain resembling the wings of "butterfly"; PC2 represents the "twist" motion, where the N-terminal domain and the C-terminal domain rotates relative to each other ( Figure 4D; Movie S2). In OPLS simulations, PC1 represents a large conformational change where both the N-terminal domain and the C-terminal domain move toward to each other, completely occluding the ATP-binding pocket ( Figure 4E; Movie S3); PC2 represents a "butterfly" motion ( Figure 4F; Movie S4). Notably, although the similarity of PCs in AMBER and OPLS simulations are limited ( Figure S7B), it appears that both force fields have a reasonable agreement with the experimental crystal structures (Figures S7C,D).

Simulations Agree With Experimental NMR Observables
Chemical shifts have been important indicators of local backbone conformations (35). For verification of our results, we compared simulated and experimental chemical shifts of apo p38α. Linear regression of chemical shifts for five atom types in two representative simulation runs is performed (Figures 5A,B). The correlation coefficient (r 2 ) suggest that there is excellent agreement for atom CB (r 2 : 0.98, 0.99) and atom CA (r 2 : 0.88, 0.89); reasonable agreement for atom C and N (r 2 : 0.42, 0.64); and bad agreement for atom H (r 2 : 0.19, 0.28). The correlation coefficient of the five atoms ranks similarly in other simulation runs ( Figure S8). The top six outlier amino acid residues were highlighted ( Figures 5A,B, red dots). To confirm the PPM results, we also performed chemical shift predictions using SHIFTX2. As expected, similar results were observed, indicating the accuracy of our simulation data (Figures S9-S11). Despite different structural flexibility in AMBER and   The thickness of the tube is normalized within each simulation run with thicker tubes corresponding to higher RMSFs. The color of the tube is normalized using all simulation runs using the "rainbow" gradient with warmer color corresponding to higher RMSFs. Color normalization was done by setting the color range to 0-8.3 Å (the largest RMSF value in all simulations). The N-terminal end (residue 4-13) and C-terminal end (residue 345-354) are not shown and excluded from the color normalization for better visualization. (D) RMSFs of p38α backbone in AMBER and OPLS simulations. The color scheme is same as in RMSD. OPLS simulations, there seems to be no significant difference between the accuracy of the chemical shifts. Furthermore, experimental residual dipolar couplings (RDCs) of apo p38α were applied for accuracy validation ( Table 3). We compare simulated and experimental RDCs for 39 residues (Figure S12A) in apo p38α (36). The probability density function of Pearson's R (between simulated and experimental RDCs) indicated that MD simulations using both force fields agree generally well with RDC experiments (Figure S12B). Linear regressions of simulated and experimental RDCs show the same trend except OPLS apo1 (r 2 is 0.00), possibly due to the surface charge distribution difference for all 85 charged residues (Table 4, Figure S13). The other five simulation runs agree generally well (r 2 is 0.47, 0.34, 0.38, 0.35, 0.19 for AMBER apo1, AMBER apo2, AMBER apo3, OPLS apo2, OPLS apo3; Figures 5C,D). Taken together, these data suggest the high quality and accuracy of our simulations.

Dynamics of p38α Inhibitor in the ATP-Binding Pocket
The commonly used p38α inhibitor, SB203580, was suggested as a potential drug for cancer therapy (37). To study the drugkinase dynamics of p38α, we measured the distances of two major interactions: (1) the hydrogen bond between the backbone of MET109 and SB203580; (2) the π-π stacking interactions between the phenyl sidechain of TYR35 and SB203580. In AMBER simulations, the hydrogen bond is well-maintained; the π-π stacking interaction is frequently disrupted with ring-ring distance up to 13 Å ( Figure 6A). In OPLS simulations, although with minor fluctuations, both the hydrogen bond and the π-π stacking interaction are well-maintained ( Figure 6A). SB203580 is less stable in AMBER simulations than in OPLS simulations, which is consistent with the wider range of RMSD fluctuations in AMBER simulations ( Figure 3A, blue lines), suggesting that the ligand stability might affect the protein flexibility and vice versa.

Current p38α Inhibitor Is Not Highly Selective
Although a variety of p38α inhibitors have been developed with enhanced specificity (38), most of these inhibitors are ATP competitors. Due to the similarity of the ATP-binding site of different kinases, off-target effects remain one of the biggest obstacles for the clinical application of p38α inhibitors. We tested the p38α inhibitors in PDAC cancer cells and various host cell lines. Indeed, both first and second generation of p38α inhibitors impedes cell viability of healthy host fibroblasts, and monocytes in a similar pattern as of cancer cells ( Figure 6B). Surprisingly, in PDAC cells treated with two classical inhibitors, the phosphorylation of other kinases such as AKT, ERK is significantly altered under a relatively modest concentration ( Figure 6C). Interestingly, the total amount of AKT was reduced under 100 µM SB203580 and 1 µM ralimetinib treatment, suggesting that an indirect regulatory pathway might be involved ( Figure 6C). These findings suggest current p38α inhibitor is not selective enough. New approaches to inhibit the p38α needs to be explored for PDAC treatment.

ASP354
Indices are consistent with that in distance matrix plots.

Potential Drug Binding Sites in p38α
To explore potential ligand-binding pockets in p38α, computational solvent mapping was performed using FTMap (39,40). Pockets with consensus cluster strength S ≥ 16 are considered druggable. Using these criteria, 19 druggable pockets were found ( Table 5): αD-L13 pocket, L4-L7 pocket, P-αC pocket, αG-L14 pocket, αE-L16 pocket, αE-L16b pocket, P-L12 pocket, L12 pocket, αE-β7 pocket, αC-αL16 pocket, β5-αL16 pocket, P-L16 pocket, MKI pocket, ATP pocket, αE-αF pocket, αH-MKI pocket, αF-αG pocket, β2-L4 pocket, αE-L4 pocket ( Figure 6D). These pockets have a varied frequency of occurrences in simulations ( Table 6) and are scattered in the entire p38α structure (Figure 6E). For comparison, FTMap analysis performed on crystal structure ensemble consisting of 196 crystal structures of p38α identified 6 druggable pockets: P-L16 pocket, ATP pocket, L4-L7 pocket, αG-L14 pocket, β5-αL16 pocket, MKI pocket, which partially verified our simulation results. Besides the well-studied ATP pocket and MKI pocket, the remaining 17 pockets are novel pockets that may be explored for further cancer drug development. Interestingly, a comparison between apo and ligand-bound simulations indicates that the presence of the ligand seems to prevent the exposure of a few binding pockets in AMBER ( Figure S14A) and OPLS simulations ( Figure S14B). Taken together, we provide potential new drug-binding sites that may pave the way for new generation drug design targeting p38α.

DISCUSSION
The role of p38 MAPKs in cancer is still under intense investigation. What is the role of p38 MAPK family members in cancer? The one-word answer is "context-dependent." Some human tumors, such as HCCs, have lower p38 MAPK activity than non-tumorigenic tissues (41). Several reports showed that p38α is a tumor suppressor. In support of that view, negative regulators of p38α, such as the phosphatases PPM1D and DUSP26, are overexpressed in human tumors (42,43). However, bearing the evidence of p38α signaling in tumor suppression, mutations of p38α have not been consistently identified in human tumors. It is suspected that cancer cells may benefit from the p38α signaling pathway. In support of that view, p38α blockades show significant tumor suppressing effects in vivo in various cancer types (44,45), suggesting the dual roles of p38α signaling in cancer. Here in our work, a thorough screening for more than 30 types of tumors and their counterparts shows that p38α expression is significantly increased in PDAC tissues, which drives our curiosity. We originally hypothesized that p38α may also be activated in PDAC cells and p38α blockades may benefit PDAC patients. Our experimental data support this hypothesis, and SB203580 significantly inhibit human and mouse PDAC cell growth. Moreover, in our recently unpublished data, SB203580 reduces the activity of pancreatic tumor derived-macrophages in in vivo models. Our data suggest that p38α blockades may block various cell components in PDAC tumor microenvironment, and may serve as a potential target for PDAC therapy. Clinically, obesity is associated with cancer risk (46), and the adipose tissue microenvironment supports cancer development, metastasis (47), and drug resistance (29). For tumors predominantly occur in the adipocyte-rich microenvironments such as breast, prostate, ovarian, colon, and pancreatic cancers, the degree of adipose tissue involvement is often correlated with poor prognosis (48). Our data shows that MAPK14 is correlated with PLIN2 and PLIN3, but not PLIN1, 4, and 5. Notably, PLIN2 and 3 is the marker for small lipid droplets. These results may reflect that cancer cells induce lipolysis of surrounding adipocyte, which is previously reported in breast cancer (49).
After establishing the role of p38α in PDAC cells, we turn to computational biology tools for more detailed information on the atomistic level for this interesting target. We performed molecular dynamics study on p38α conformational dynamics  on the state of the art, highly parallel supercomputer, Anton. Supercomputer-powered biomolecular simulation provides us abundant information related to p38α function. Initially, a powerful and direct metrics, all-too-all RMSDs were calculated to examine the convergence of MD simulations. We observed the frequent appearance of low RMSD stripes in both AMBER and OPLS simulations (Figure S3A, dark blue areas), indicating a reasonable convergence. The results show that across different simulation runs, the conformational similarity is quite low, especially for OPLS apo1 (6.3 µs), which is significantly different from all other simulation runs ( Figure S3B, mid panel, lower panel, red areas), possibly due to large conformational change of the glycine-rich loop in OPLS apo1 simulation ( Figure S5). We found 11 p38α crystal structures with glycine-rich loop conformations that at least partially resemble the unfolded conformation found in MD simulations, suggesting that the novel conformation of the glycine-rich loop might be realistic.
The high flexibility we found in local regions of apo p38α is generally in line with previous simulation studies. Kuzmanic (53). These data support that the conformational flexibility of local regions might be a common feature in protein kinases, presumably conferring adaptability in enzyme catalysis and protein-protein interactions.
The "butterfly" motion identified in PCA analysis is a quite interesting protein dynamics motion. Similar motions were seen in simulations of other enzymes as well such as the adenylate kinase (54). Many two-lobe enzymes may show similar type of motions during catalysis. Meanwhile, it is possible that both "conformational selection" and "induced fit" (55)(56)(57) contribute to the ligand-binding because the PC1-PC2 subspace of apo and ligand-bound p38α has both overlapping and non-overlapping regions. Nevertheless, the PCA result is complicated by two factors. Firstly, the initial structure of apo and SB203580-bound p38α simulations are similar (with RMSD of only 0.5 Å). Therefore, it is highly likely to have structural overlaps. Secondly, the apo and SB203580-bound p38α crystal structures come from mouse and human, respectively. The two structures differ at residue 48 and 263 (HIS and ALA for the mouse form; LEU and THR for the human form), which may be a potential source of error.
Despite the different conformational changes, linear regression of simulated and experimental chemical shifts shows similarly good agreement in AMBER and OPLS simulations, indicating that chemical shift is not sensitive to conformational changes. In contrast, linear regression of simulated and experimental RDCs shows a lower correlation coefficient in OPLS apo1 than other simulation runs, suggesting that RDC is sensitive to conformational changes. Large conformational changes of both N-terminal and C-terminal domains in OPLS apo1 may also contribute to the relatively worse RDC agreement. Since accurate RDC predictions are dependent on the precise distribution of the surfaces charges (58), we calculated all-to-all distance matrix for all 85 charged residues (Table 4) in p38α including ARG, LYS, ASP, and GLU for the simulation run AMBER apo1, AMBER apo2, and OPLS apo1. The charged residue distance matrices have a much bigger discrepancy between AMBER apo1 and AMBER apo2 than between AMBER apo1 and OPLS apo1 (Figure S13), suggesting that surface charge distribution can be an important contributing factor in RDC predictions.
Newly discovered potential druggable pockets in this study may benefit millions of PDAC patients. Notably, ATP pocket, MAP kinase insert (MKI) pocket, and αC-αL16 pocket have been experimentally verified (13,59). Several novel pockets are in proximity to the binding site of upstream activators or downstream substrates of p38α and can potentially be explored to modulate its activity. There are at least four pockets that can be potentially used for this purpose: (1) αD-L13 pocket, which can be used to block docking of kinase substrates; (2) αC-αL16 pocket, which can be used to lock p38α in inactivate conformation; (3) αG-L14 pocket, which can be used to block p38α-TAB1 interaction (60); (4) αE-β7 pocket, which can be used to block p38α-MK2 interaction (61) and p38α-MKP5 interaction (62). Interestingly, the binding of the inhibitor seems to veil many potential binding pockets in AMBER and OPLS simulations (Figure S14). Further, virtual screenings and experimental validations are needed to confirm the druggability of these novel pockets.
Taken together, our study provides an interesting potential target for combating PDAC and offers detailed conformational dynamics information on p38α protein. Our findings have also paved potential avenues for developing the new classes of p38α-targeting drug, which may overcome the side effects of current therapy and benefit PDAC patients.

Cell Culture
Murine

Human Patient Samples
All studies related to clinical human samples were approved by the Ethical Review Committee in Shuguang Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai, China. Pancreatic tumor samples and adjacent pancreatic tissues were collected from cancer patients with written informed permission.

RNA Extraction and Quantitative Real-Time PCR
Total RNAs were extracted from tumor tissues and cultured cells using RNAsimple Total RNA kit (Cat. No. DP419, TIANGEN, China). Total RNA from each sample was reversely transcribed using an All-in-One cDNA Synthesis SuperMix (Cat. No. B24408, Bimake, China). Reverse transcription was performed at 42 • C for 60 min, followed by 70 • C for 5 min to inactivate the enzyme activity. cDNA samples were stored at −20 • C and subjected to qPCR using a StepOnePlus Real-Time PCR System Simulation System Preparation p38α protein structures were prepared using Maestro 9.0 (Schrödinger, LLC). The crystal structures of apo p38α (PDB code: 1P38) (9) and SB203580-bound (PDB code: 1A9U) p38α (10) were used as initial structures. The protein was capped with acetyl group (ACE) at the N-terminus and N-methyl group (NME) at the C-terminus to prevent artificial chargecharge interactions. AMBERff99SB-ILDN (64) and OPLS-AA/L (65) force fields were used for the protein. A box size of 85.0 × 85.0 × 85.0 Å 3 was used to ensure the enclosure of the entire protein. Apo p38α or SB203580-bound p38α system was solvated in water described by the TIP3P model (66) to mimic physiological condition and neutralized using nine sodium ions (67).

Ligand Parameterization
To derive parameters for the SB203580, we used different approaches. For AMBER simulation, the bonded parameters and non-bonded parameters for SB203580 were parameterized using General Amber Force Field (GAFF) (68). The partial charges were derived using the restrained electrostatic potential (RESP) method on the RED server (69,70). For OPLS simulations, the bonded parameters, non-bonded parameters, and partial charges were derived by analogy method as previously described (71). In both simulations of ligand-bound p38α, dihedral restraints were applied to eight dihedral angles to keep the ligand conformations close to the crystal structures. The force constant for dihedral restraints is 4184 kJ/mol.rad. Each equilibrium angle ( Table 1) was obtained from four SB203580-bound p38α crystal structures [PDB codes: 1A9U (10), 1PME (72), 2EWA (17), and 3GCP (73)].

Molecular Dynamics Simulations
Before the production simulation, energy minimization, and equilibration were performed on our in-house clusters using Desmond 2.4 (74). Energy minimization was initially carried out using the steepest descent algorithm for 1,000 steps and the L-BFGS algorithm for 1,000 steps. The energy-minimized structure was incrementally heated up from 50 to 298 K (50, 100, 150, 200, 250, 298 K) with 100 ps at each temperature in the NVT ensemble using Berendsen thermostat (75). Next, the protein was equilibrated at 298 K for an additional 10 ns in the NVT ensemble using the Nosé-Hoover thermostat (76,77). For van der Waals interactions and short-range electrostatic interactions, a 9 Å cutoff was used; for long-range electrostatic interactions, Particle Mesh Ewald (PME) was used (78). The "bonded, " "near, " and "far" time steps are 1, 1, and 3 fs, respectively.

Anton Supercomputer Molecular Dynamics Simulations
For production simulations, we performed simulations on the Anton supercomputer (Pittsburgh supercomputing center, Pittsburgh, PA) (33). Production simulations were performed in the NVT ensemble with the temperature maintained at 298 K using Nosé-Hoover thermostat (76,77). For van der Waals interactions and short-range electrostatic interactions, the cutoff was automatically determined by the "guess_chem" utility on Anton for optimized performance. Gaussian Split Ewald (GSE) (79) method was used for long-range electrostatic interactions. M-SHAKE algorithm (80) was used to constrain all bonds, enabling a 2 fs time step. Snapshots of the simulations were saved at an interval of 0.1 ns for analysis. Visualizations of simulation snapshots and trajectories were done using both PyMOL 1.8 (PyMOL 2015) and VMD 1.9 (81). For AMBER and OPLS simulations, three independent simulations of apo p38α and one simulation of SB203580bound p38α were performed with detailed parameters shown in Table 2.

Calculations of RMSD and RMSF
To obtain root mean square deviations (RMSDs) in p38α simulations, the trajectory (non-terminal residues, residue 14-344) at an interval of 0.1 ns was fitted to the initial structure on protein backbone; RMSDs were calculated for each snapshot over simulation time. To obtain all-to-all RMSDs, the trajectory at an interval of 1 ns was fitted to the initial structure on the protein backbone before all-to-all RMSD matrix was generated using pytraj (82), a Python package binding to cpptraj program (83). To obtain root mean square fluctuations (RMSFs) in p38α simulations, the trajectory at an interval of 0.1 ns was fitted to the initial structure on protein backbone; RMSFs were then calculated for each residue.

Calculations of Chemical Shifts and Residual Dipolar Couplings
To obtain chemical shifts of the backbone atoms in p38α from MD simulations, we performed ensemble-based chemical shift calculations on MD trajectory at an interval of 1 ns using the PPM program (86). PPM predicts the chemical shifts of CA, CB, C, N, and H based solely on the physical and chemical properties of the protein structure. For comparison, we also performed chemical shift calculations on MD trajectory at an interval of 1 ns using the SHIFTX2 program (87). SHIFTX2 predicts the chemical shifts of CA, CB, C, N, and H using both structure-and sequence-based criteria. Experimental p38α chemical shifts were obtained from Biological Magnetic Resonance data bank (BMRB entry number: 6468) (88); the entry includes 219 1 H chemical shifts, 684 13 C (including CA, CB, and C) chemical shifts, and 219 15 N chemical shifts for apo p38α. For comparison of simulated and experimental chemical shifts of atom CA, CB, C, N, and H, different residues were used depending on the availability of experimental chemical shifts for specific atoms in BMRB entry 6,468. To obtain residual dipolar couplings (RDCs) of each residue in p38α, we performed RDC calculations on structural snapshots at an interval of 1 ns using the PALES program (58,89). Experimental RDC data for 39 residues of apo p38α and the parameters in RDC ( Table 3) were obtained from previous report (36). Specifically, a concentration of bacteriophage Pf1 at 20 mg/mL and a pH of 6.0 was used. The sodium chloride concentration was optimized using a range from 0.005 to 0.5 M and 0.2 M was used for final RDC calculations.

FTMap Analysis
A standalone version of FTMap (39,40) was used to perform calculations on MD snapshots at an interval of 1 ns. FTMap analysis mainly consists of three steps: (1) independent docking of 16 organic probes on the protein using a fast Fourier transform correlation approach; (2) refinement of the probe positions using energy minimizations; (3) clustering and ranking of the resulting poses to identify consensus sites (CSs

Classification of Binding Pockets
Binding pockets are named and classified based on its adjacent secondary structures, which harbor pocket lining residues ( Table 5). To determine the identity of a pocket, the intersecting volume of the bounding sphere for the probe clusters and the pocket lining residues were calculated for each pocket using the equation as described previously (90): for d ≥ r 1 + r 2 π 12d r 1 + r 2 − d 2 d 2 + 2d (r 1 + r 2 ) − 3(r 1 − r 2 ) 2 for |r 1 − r 2 | < d < r 1 + r 2 4 3 π(min {r 1 , r 2 }) 3 for 0 ≤ d ≤ |r 1 − r 2 | V o is the intersecting volume; r 1 , r 2 are the radii of the bounding sphere of the probe cluster and of the pocket lining residues, respectively; d is the distance between the two sphere centers. A specific probe cluster falls in a pocket category if it has the maximum intersecting volume with that pocket. A scaling factor of 0.75 on r 2 was used to avoid the overestimation of the bound site volume. All pockets and subpockets in the ATP-binding site are classified as ATP pocket.

Statistical Analysis
Statistical analyses of biological experiments were performed using the standard two-tailed Student t-test, and P < 0.05 was considered statistically significant.

DATA AVAILABILITY STATEMENT
The datasets analyzed in this study can be found in The Cancer Genome Atlas (https://portal.gdc.cancer.gov/).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics committee, Shuguang Hospital, Shanghai University of Traditional Chinese Medicine. The patients/participants provided their written informed consent to participate in this study.