Cell-Specific Computational Modeling of the PIM Pathway in Acute Myeloid Leukemia

The potential of personalized medicine is dependent on our ability to translate the molecular context of patients' tumors into interpretable clinical outcomes. Successful steps have been taken to accurately predict tumor progression and response to treatment from molecular disease markers (1, 2). Using tumor cell line–based compound screening, we can provide robust readouts of cellular responses to multiple compounds. This information can be used to systematically train computational models of the molecular signaling pathways contributing to drug sensitivity and resistance in various cancer settings, and to propose novel drug targets and combination approaches. Cell line screens have provided some success in explaining or predicting drug responses by driver gene mutations (2–4); however, in many cases the true mechanism of resistance remains elusive or more complex. Most predictive methods routinely used today use correlative statistics or feature-based learning techniques such as machine learning, while network methods remain scarce despite their potential for extracting mechanistic insights and actionable biomarkers.

The molecular heterogeneity within cancer types further complicates the prediction of tumor cell behavior determining a patient's drug response. Multiple somatic mutations, epigenetic events or otherwise deregulated gene/protein expression may contribute to driving the disease. This is true in acute myeloid leukemia (AML), where patients may harbor somatic mutations in a number of potential oncogenes, including FLT3, MLL, TYK2, FGFR1, PDGFRA, IDH1, DNMT3A, affecting expression of downstream signaling for example through PIM kinases (5–7). FLT3 internal tandem duplications (FLT3-ITD) and PIM overexpression are associated with poor prognosis in AML patients, motivating the development of small molecule inhibitors targeting these proteins (8, 9). Incomplete signaling inhibition or the presence of multiple molecular alterations that reduce a tumors dependency on any one target may result in drug resistance (10, 11). This may be overcome through rational drug combinations; however, optimal approaches are rarely obvious and high-throughput combination screening is complex and expensive with limited success shown.

With an aging population, the incidence of AML is increasing, with the number of new cases per year approaching 20,000 in the United States alone. AML therefore presents a large unmet clinical need, with overall 5-year survival rates remaining at around 25%. Most patients will respond to initial cytoreductive therapy, but a large proportion will relapse with emergence of drug-resistant clones. Given that bone marrow transplantation as the only curative therapy is not an option for many patients, a better understanding of the regulatory pathways causing leukemic transformation and in particular the emergence of resistance will be essential to improve treatment outcomes in AML.

Computational simulations of cancer cell signaling have the potential to overcome both the limitation of cell line diversity and in vitro screening throughput. Computational modeling approaches can be used to capture and integrate knowledge with molecular and phenotypic data to better understand the genetic and signaling dependencies determining a drug's mechanism of action. The models should be unique to the tumor cell context, include key proteins and their interactions while accounting for influential gene mutations, and would ideally extend to other molecules involved in cell signaling. Execution of such models should demonstrate the intracellular signaling activity as it is triggered by different mutations and different therapeutic modalities, resulting in different cell phenotypes.

Computational models based around Boolean networks, pioneered by Kauffman (12) as a model for genetic regulatory networks, have been demonstrated for interpretation of large data sets as well as for drug discovery (13–15). In these models, relationships are represented in a dynamic network with discrete time steps. Signaling molecules represented by nodes in a network can have two states (hence a Boolean network) and edges are directed and may be activating or inhibitory; however, this can oversimplify biological signaling where molecules often exist in multiple states with interactions that are rarely binary. Qualitative Networks (QN) make an extension to Boolean networks to allow variables to range over larger discrete domains by replacing Boolean functions with algebraic functions (further details in Supplementary Methods; ref. 16). Specifically, the graphical tool Bio Model Analyzer (BMA; ref. 17; available at http://biomodelanalyzer.org/) has previously been used to encapsulate chronic myeloid leukemia (CML) cell signaling information from >150 publications in a QN model (18) able to then successfully recapitulate multiple independent experimental results. Another extension to Boolean networks is provided by Quantitative Modeling approaches, allowing variables to range over nondiscrete values and so capturing more complex relationships, but only feasible for much smaller, well-studied systems (19).

In this study, we use QNs to model the protein signaling connecting genetic aberrations in FLT3, TYK2, PDGFRA, or FGFR1 to cell proliferation/apoptosis via the PIM and PI3K kinases for four AML cell lines, accounting for their unique genetic and phenotypic diversity. Construction and analysis of the biological QN model was achieved in BMA (17). By incorporating cell-specific–context switches in the model for four cell lines, we were able to accurately model response and resistance to a pan-PIM kinase inhibitor AZD1208 and the FLT3 inhibitor AC220 and validated experimentally our predictions. The model provides a useful tool for AML research, and the approach offers value to drug discovery and early development.

AZD1208, selumetinib, and AZD5363 were synthesized by AstraZeneca R&D and diluted in dimethyl sulfoxide (Sigma-Aldrich). AC220 and pictilisib were purchased externally.

Cell lines (CMK, EOL1, HL60, KASUMI3, KG1A, MOLM13, MOLM16, MONOMAC6, MV411, NOMO1, OCIAML2, OCIM1, and OCIM2) were purchased from ATCC (http://www.atcc.org/) cell bank and passaged in our laboratory for fewer than 6 months after receipt or resuscitation. ATCC uses morphology, karyotyping, and PCR-based approaches to confirm the identity of human cell lines and to rule out both intra- and interspecies contamination.

All cells were cultured and assayed as previously described in ref. 8 and in Supplementary Methods.

For single agent, GI₅₀ were calculated from the ratio of the 72-hour treatment to 72-hour DMSO control, after subtraction of the day 0 data from each measurement. The dose–response data were fitted using Xlfit (Microsoft Excel). For combination, percent growth inhibition was determined using the Chalice software with values of 0 to 100%, indicating antiproliferation (fewer number of cells than the vehicle control but greater than or equal to the number of cells at the start of treatment) and values of 101% to 200%, indicating cell death (fewer cells than at the start of treatment). Day 0 values were subtracted from the day 3 treatments. The combination Indexes (CI) and Synergy scores were determined using the software program Chalice (Zalicus), and CI determination was made at the ED₅₀ value. Synergy was determined by the Loewe additivity model.

Full methods for gene expression microarray, whole-exome DNA sequencing, and Theranositcs proteomics are in the Supplementary Methods; a brief description follows.

Cell line lysate was generated from logarithmic growing CMK, EOL1, HL60, KASUMI3, KG1A, MOLM13, MOLM16, MONOMAC6, MV411, NOMO1, OCIAML2, OCIM1, and OCIM2 cell lines. Lysate was sent to Expression Analysis (http://www.expressionanalysis.com/) for gene expression analysis on Affymetrix Human Genome U133 Plus 2.0 Array. Expression results were fRMA normalized, log₂ transformed, and expression was averaged by gene symbol across probesets.

Cell lines lysate was generated from logarithmic growing CMK, EOL1, HL60, KASUMI3, KG1A, MOLM13, MOLM16, MONOMAC6, MV411, NOMO1, OCIAML2, OCIM1, and OCIM2. Lysate was sent to Expression Analysis (http://www.expressionanalysis.com/) for whole-exome DNA sequencing and processed with the BCBio pipeline (https://bcbio-nextgen.readthedocs.org). Paired analysis of the parental and resistant cell lines was performed to using FreeBayes (20), MuTect (21), and VarDict (GitHub) to call resistance–specific mutations.

Cells were treated with AZD1208 or AC220 as single agent or in combination for 3 or 24 hours. Lysate was prepared and shipped to Theranostics Health for reverse-phased protein array (RPPA) experiments.

The relative linear log₂ RPPA values were categorized for use in executable network modeling (Supplementary Fig. S1) on a 5 point scale from 0 to 4 fitted to the distribution of the values (Supplementary Table S1).

Cells were treated with AZD1208, AZD5363, or AC220 as single agent, combination, or resistance as experimentally described. Whole-cell extracts were fractionated by SDS-PAGE and transferred to a nitrocellulose membrane in transfer buffer (500 mmol/L glycine, 50 mmol/L Tris–HCl, 0.01% SDS, 10% methanol) buffer at 20 volts for 90 minutes using a semi-dry transfer apparatus according to the manufacturer's protocols (Invitrogen). The membranes are blocked with 10% nonfat milk in TBS-T (10 mmol/L Tris, pH 8.0, 150 mmol/L NaCl, 0.5% Tween 20) for 1 hour and washed three times with TBS-T and exposed to primary antibodies in 5% milk in TBS-T against pPRAS40 (CST 2997), p4EBP1 Ser 65 (CST 9451), pBAD (CST 9296), pp70S6 (CST 9206), pS6 (CST 4858), pERK (CST 9106), pElF4B (CST 8151), pAKT(CST 4058), α-tubulin (CST 2144), or β-actin (CST 4970) at 4°C O/N. Membranes are washed three times for 10 minutes and incubated with a 1:10,000 dilution of horseradish peroxidase–conjugated anti-mouse or anti-rabbit antibodies (CST 7074) for 1 hour at room temperature. After washing the membranes three times for ten minutes, signals were visualized using the ECL system (Thermo Scientific).

We confirmed the robustness of our finding for MOLM16 cells treated with 2 μmol/L AZD1208 for 3 hours (Supplementary Table S2) by applying a LC-MS/MS phosphorylation proteomic approach. Additional details are provided in the Supplementary Note.

We investigated phenotypic and cell signaling responses by RPPA. Because the PIM gene family is often overexpressed (22–24) and FLT3-ITD's are prevalent in AML (7), we treated the cells with the pan-PIM kinase inhibitor AZD1208 and the potent selective FLT3 inhibitor AC220 (Qizartinib) as monotherapy and in combination (25) and compared AML cell lines were treated with DMSO, 1 μmol/L AZD1208, 6 nmol/L AC220, or the respective combination, for 3 or 24 hours. Lysates were generated, protein values were assessed by RPPA and quadrant median normalized (QMN) protein levels calculated (Supplementary Table S3). Statistically significant total protein and phosphorylation changes were determined by log₂ QMN differences greater than or equal to 0.5 and Wilcoxon rank-sum tests P values less than or equal to 0.1 (Supplementary Table S4).

To identify potential genetic alterations associated with sensitivity to the pan-PIM kinase inhibitor AZD1208, we surveyed gene variants by whole-exome DNA sequencing (Supplementary Table S5) and prioritized by AML disease relevance (7). Although cells sensitive to AZD1208 do harbor AML relevant PDGFRA, FGFR1, FLT3, and MLL genetic variants, only a small number of cell lines harbor the same variant, thereby failing to reach statistical significance in association to drug response (Fig. 1A). Basal cell line PIM1 mRNA expression tends to be higher in sensitive lines, as previously shown at the protein level (8), underlying the importance of compound target expression alongside the interplay with genetic alterations for sensitivity. However, cells harboring pathway relevant genetic alterations or overexpressing PIM exhibit varied response to treatment, calling for a deeper examination of the cell signaling relating genotype to phenotype to provide a better understanding of the molecular dependencies underlying PIM inhibitor sensitivity in AML cell lines.

Given the wide variability of response to therapeutic agents across AML cell lines, we explored the differences in phospho-protein signaling downstream of PIM for AML cell lines MOLM16, MV411, EOL1, and KG1A. RPPA measurements taken 24 hours post pan-PIM inhibitor treatment reproduced published findings (8) of reduced BAD phosphorylation in the MOLM16 cell line and reduced S6 pS235/236 in EOL1 (Fig. 1B). To estimate response of cell lines, growth inhibition was quantified according to the number of viable cells after culturing with different concentrations of pan-PIM inhibitor, AZD1208, and FLT3 inhibitor, AC220, in combination (Fig. 1C; Supplementary Fig. S2). Directional de-phosphorylation signaling trends seen in RPPA for BRAF pS445, EIF4B pS406, mTOR pS2481, and global BAD phosphorylation were confirmed by PhosphoScan mass spectrometry in MOLM16 cells after 3-hour treatment with AZD1208 (Supplementary Table S6).

In order to model the observed genotypic and phenotypic differences between the AML cell lines, we proposed a workflow for developing a cell-specific context network model using the BMA tool from cell line molecular information (Fig. 2). We generated an initial generalized model from the manual curation of 68 publications (Supplementary Table S7) for the AML cell line. The initial model contains a canonical set of 64 interactions among 32 interacting proteins connected to 2 cell phenotypes/behaviors of apoptosis and cell proliferation (http://www.bioc.cam.ac.uk/fisher/aml; GeneralModel.json, Supplementary Table S8). All values at nodes range from 0 to 4 to represent the phosphorylation activity from the transformed RPPA data, with 0 representing low to no activity and 4 representing abnormal over activity. The cellular behavior outcome for each disease state is reflected by the two terminal downstream nodes, which model the outcome for cellular abnormal proliferation and apoptosis rates. The generalized model of AML signaling was able to capture only partial abnormal cell behavior for untreated cells, capturing the abnormal low apoptosis levels for both MOLM16 and MV411, and showing an increase in proliferation, yet not capturing the magnitude of the increase. In addition, known perturbations such as simulating inhibition of PIMs in the model, showed the expected trend line of decreased proliferation, yet did not exhibit the expected effect on apoptosis levels.

We incorporated multiple gene mutation switches to construct cell-specific context model (http://www.bioc.cam.ac.uk/fisher/aml; http://www.bioc.cam.ac.uk/fisher/aml/CellSpecificAML.json/view; Table 1). We iteratively refined the target function for each internal node to reflect the levels of phosphorylation activity as measured by the transformed RPPA data for each cell line as well as the qualitative activity reported in the literature in accordance to gene mutations. A cell-specific context in the model is simulated by setting the switches for the driver mutations found in that cell to 1, while all other switches to mutations that were assessed as non-driver are set to 0 (Supplementary Table S9). As a result of different set of mutations “turned on” the protein activity exhibited by the model will differ between cell lines (Supplementary Table S10). Additional details on data processing and model construction are in the Supplementary notes.

The executable QN model (Fig. 3A) was built on the RPPA and growth inhibition of MOLM16 and MV411 cell lines, harboring TYK2 mutation and FLT3-ITD, respectively (Fig. 1B). For each cell line across each treatment (Fig. 3B), the mean square error (MSE) observed between the transformed RPPA values and modeled signaling activity ranged from 0.17 to 0.26 and median of 0.2 (0.3 to 0.57 in the generalized model, median of 0.41), with the lowest seen for untreated MV411 and MOLM16 cell lines and the highest for the MOLM16 cell line treated with AZD1208. Meanwhile, across each protein, the MSE observed between the transformed RPPA values and each protein signaling activity ranged from 0 to 0.55 and median of 0.25 (0.5 to 0.88 in the generalized model, median of 0.58), with the lowest seen for BAD and BCR and the highest seen for p27.

Equally as important, the cell-specific context model performed well in predicting cellular response as measured both by growth inhibition and markers of reduced proliferation and increased apoptosis (Fig. 1B and C). The model accurately predicted (Fig. 3C) the reduction in proliferation as a result of treatment with AZD1208 single agent, AC220 single agent, and drug combinations in MV411 cells. Although under predicting the magnitude of increase in apoptosis for AC220 single agent, the model accurately predicted the directional responses with increases in apoptosis for AZD1208 single agent, AC220 single agent, and combination treatments in MV411 cells.

In addition to predicting differential phenotypic responses in each cell line, the model highlights key signaling events that may underlie the mechanism for each. We validated the robustness of events suggested for MOLM16 using mass spectrometry. Most importantly, the mass spectrometry corroborated the decreased EIF4B pS406 phosphorylation after AZD1208 treatment, contributing to decrease in proliferation, as well as the decrease in BAD pS112 and pS155 after AZD1208 treatment, which increases apoptosis. A key differentiating feature of MOLM16 cell lines is the lack of hyperactivity from the MAPK (Ras Raf MEK ERK) and AKT–mTORC1 pathways post AZD1208 treatment, supported by dephosphorylation at downstream BRAF pS445 and mTOR pS2481 in the mass spectrometry data.

Testing the adaptability of the model to new cell-specific contexts, we "turned on" new genetic alterations FIP1L1-PDGFRA fusion to simulate EOL1 cell line and FGFR1 fusion to simulate KG1A cell line. The apoptosis range was expanded to span the full dynamic range seen in these cell lines, yet no further refinement of the model was performed. The model reflected the cellular signaling changes observed in RPPA data (Fig. 3D) where the MSE ranged from 0.18 to 0.28 with the lowest seen for KG1A cell line treated with single agent AZD1208 and the highest for the KG1A cell line treated with AC220. Across each protein, the MSE observed between each protein signaling activity and the transformed RPPA values ranged from 0 to 0.49, with the lowest seen for AKT and the highest seen for BAD. The model also performed well in predicting cellular response (Fig. 3E). For the proliferation and apoptosis cell behaviors, the model accurately predicted the cellular responses seen in KG1A for AZD1208, AC220, and combination treatments, as well as the cell behaviors for EOL1 with AZD1208 treatment (AC220 was not tested for EOL1).

The model also replicated variations in sensitivity, such as EOL1, reacting with reduced apoptosis to PIM inhibition when compared with MOLM16.

A by-product of refining the QN model to capture cell type–specific signaling is a graphical and descriptive representation of cell-specific signaling dynamics between proteins in the network (Fig. 3A). By simulating the QN model, we were able to test our assumptions regarding the signaling dependencies between proteins, as described by the target functions (Table 1). For instance, despite FLT3-ITD being upstream of PIM1, the effect revealed by the iterative optimization of the model was less than other interacting proteins, also suggested by the RPPA measurements (Fig. 1A), leading to BAD overactivity in MOLM16 but not MV411. The target function of AKT shows that it is dependent on the activity of the FLT3-ITD and FGFR1 fusion, reflecting the accumulation of evidence for AKT/mTOR pathway role in AML (suggested previously by ref. 8). The target function of S6 reflects the dominant over activation of it via AKT–mTOR pathway, additive to the activity of the MAPK pathway, and leading to antiapoptotic cell behavior of MV411 and KG1A. At the same time, the target function of BAD accumulates with activity of the MAPK pathway and of PIM1 direct phosphorylation of all three sites of BAD (26), leading to the antiapoptotic behavior observed in MOLM16.

As a first independent test of the AML cell–specific model, we assessed its ability to replicate in silico a sample of protein and phenotypic cell line responses to drug treatment reported in the literature but not used as part of model construction or refinement. We replicated each in vitro experiment by turning on a respected set of mutations and adding the new examined inhibitor to the model, then observing the predicted protein expression. All eight protein changes were successfully predicted (Supplementary Table S11). The model successfully predicted cell-specific response to compounds including failure of a MEK inhibitor to induce apoptosis in EOL1 (27); insensitivity of KG1A to the combination of AKT, PDK1, and FLT3 inhibitors (28); and the growth inhibition induced on EOL1 by combining PIM and AKT inhibition (11). The decrease in cell proliferation of MV411 in response to the mTORC (29) inhibitor was not recapitulated; however, Willems and colleagues (29) attribute the decrease in proliferation to eIF4E decreased expression, which was accurately replicated by the model.

To assess the potential to prioritize synergistic combinations through in silico hypothesis testing with these models, we assessed the PIM inhibitor AZD1208, the AKT inhibitor AZD5363, MEK inhibitor selumetinib (selumetinib, ARRY-142886), FLT3 inhibitor AC220, and PI3K inhibitor pictilisib across the 4 AML cell lines (Fig. 4A; Supplementary Fig. S3) also summarized in Supplementary Table S12. For each cell line, we simulated inhibition of the drug targets first as single agents and then as combinations with PIM inhibition. We validated each combination in each cell line experimentally across a dose range for each agent (Fig. 4B).

The MOLM16 cell line was correctly predicted to be hypersensitive to the PIM inhibitor, resulting in almost complete cell kill, and no additional effect was predicted in combination with other inhibitors.

In contrast, the MV411 context model, which harbors a FLT3-ITD, correctly predicted a strong synergy between AZD1208 and AC220 combination attributed to apoptotic effect, evident even at lower dosage of combined treatments. Very weak synergy with mild apoptosis was correctly predicted in MV411 in combination with either MEK or PI3K inhibition.

Meanwhile, EOL1 was correctly predicted to gain apoptotic synergic effect with the PIM and AKT inhibitor combination, as well as the PIM and PI3K inhibitor combination. Surprisingly, and the only synergy of the 16 combinations not predicted by the model, EOL1 also exhibited a synergic effect with the AZD1208 and AC220 combination. AC220 efficacy has previously only been reported in FLT3 driven tumors; however, these data suggest efficacy from AC220 in PDGFRA mutated tumors potentially through inhibition of PDGFRA-driven AKT/PI3K and MAPK signaling.

Finally, in the KG1A model, which harbors an activating FGFR1 fusion, we did not see a co-occurrence of high apoptosis and high growth inhibition for any of the combination treatments, validated as well by the in vitro assays. Our model suggests that the persistent insensitivity of KG1A may be derived by the high levels of cMyc, which is not directly targeted by any of the combinations.

Four separate populations of MOLM16 cells were made resistant to PIM inhibition by growth in the presence of increasing doses of the compound over a 4-month period until resulting cell populations were able to maintain growth at 1 μmol/L AZD1208. While the parental MOLM16 cell has a 50 nmol/L GI₅₀ in a 3-day MTS proliferation assay, all four resistant populations had GI₅₀s greater than 9 μmol/L to AZD1208 over the same 3-day growth period (Supplementary Fig. S4A). RPPA measurements were taken for the parental and resistant cell lines.

We predicted candidate genetic causes of resistance by iteratively perturbing all individual and pairs of nodes in the parental MOLM16 model, and choosing those leading to similar signaling activity and phenotype as observed in the resistant populations, quantified by lower MSE (Supplementary Fig. S4B and S4C). This resulted in four different resistant contexts, one for each resistant cell population (Fig. 5A). All contexts show overactivation through RAS/PI3K as well as AKT/mTOR signaling, supported by RPPA (Supplementary Fig. S4B). Interestingly, the different resistance contexts differ in their strength of altered signaling where resistant cell population R1 has a higher activity for both pathways and resistant cell population R3 has lower activity for the AKT pathway. The predicted and observed pathway signaling suggests increased signaling activity through 4EBP1, EIF4B, S6, and BAD contributing to resistance. In particular, it highlights the AKT–S6 pathway as a major cause for the decreased apoptosis compared with MOLM16 parental when treated with AZD1208.

Whole-exome DNA-seq was performed to identify potential protein altering genetic variants that could be driving the AZD1208 resistance. All variant calls with significant differences from the parental line (Supplementary Table S13) were further parsed to highlight genes encoding proteins that have BIOGRID interactions (Supplementary Table S14) to the RAS/PI3K and/or the AKT/mTOR signaling pathways (Fig. 5B).

Using the four resistant MOLM16 context models, we predicted possible treatments to overcome resistance by simulating inhibition at each point through systematic addition of an inhibitor node to the network. In line with signaling changes, introduction of an AKT inhibitor, AZD5363, to the resistant populations was predicted to overcome the AZD1208 resistance by blocking the abnormal PRAS40, 4EBP1, and S6 activity (Fig. 5D). To test this prediction, parental MOLM16- and AZD1208-resistant populations were treated with and without 1 μmol/L AZD5363 for 1 hour. The resistant populations responded to AKT inhibition with AZD5363 by decreased pS235/235 S6 ribosomal protein and pT246 PRAS40 (Fig. 5D), providing strong evidence for inhibition of AKT/MTOR signaling. The decrease in AKT/mTOR signaling was accompanied with an increase in cleaved PARP, indicating increased apoptosis and highlighting the dependency on this signaling pathway during AZD1208 resistance in MOLM16 cells.

Qualitative models provide coarse-grained descriptions useful for systems whose mechanistic underpinnings remain incomplete. The range of qualitative modeling approaches provides two major types of simplifications: Boolean models relax the activity of biological entities to binary (ON or OFF), alternatively, the relation of entities may be relaxed to simple logic operators (AND, OR, NOT). We explored the use of alternative approaches and robustness of findings by building a Boolean model and an AND/OR model via the same pipeline. For single-agent PIM-inhibitor treatment, the Boolean model was able to reasonably predict the proliferation and apoptotic responses in MOLM16 and KG1A, partially predicted proliferation response in EOL1, but poorly predicted responses in MV411 (Supplementary Fig. S5). The MV411 cell line was correctly predicted to response well to FLT3 inhibition. The model was not, however, able to predict treatment combination synergies (Supplementary Fig. S5). Because the Boolean model is simpler and easier to construct than a qualitative model it offers a useful tool for investigating single-agent treatment in larger networks.

The AND/OR gated model recaptured most of the responses to single treatments, as well as synergistic combinations, revalidating the predictions made by our model (Supplementary Fig. S6). The synergistic response of KG1A to the combination of AZD1208 and AC220 was the only response not recaptured. This phenotype is likely derived by S6 additive activity from the MAPK pathway and AKT–mTORC1, which cannot be accurately described using AND/OR gates. AND/OR models may be generated by automated tools (30) and can serve well as an initial model scaffold. However, more complex relationships such as those in our model between BAD, S6, 4EBP1, TSC2, and EIF4B in AML need to be further refined.

The success of personalizing treatments for AML patients by tailoring to respective genetic alterations that characterize cancer subtypes has so far been limited. Moreover, drug responses seen in genetically matched patients or representative cell lines show considerable diversity (7, 10). By integrating both genomic and baseline proteomic data from AML cell lines with known tumor-driving genetic events we generated an AML network model capturing cell-context–specific signaling in the PIM kinase pathway. We developed a workflow methodology for constructing a network model with cell-specific context switches, which focuses on iterative refinement of the target function to reflect literature and experimental evidences. Users may also consider applying automated tools to decipher the target functions, such as CellNOpt-cFL tool developed by Morris and colleagues (30), and follow by manual refinement of the target functions.

The resulting cell-specific model captures cell-specific signaling and response to cancer therapeutics and provides virtual cell line models in which to test hypotheses for tailored therapy in silico. The cell-specific model significantly reduced the prediction error for both the baseline training data and on-treatment changes in protein expression compared with the generalized model. This is unsurprising because a generalized AML model insufficiently explains the heterogeneity in the mutational landscape and protein-signaling dynamics reported across different cell lines, for example, a lack of signaling through AKT unique to cells with mutations in TYK2.

The cell-specific model accurately and directly recapitulated published experimental results for reported changes in expression in all 8 cases, and 9 out of 10 responses in cell behavior. These results are particularly remarkable when considering the potential variability in signaling and phenotypic output over time, and the focus of these models on the cells steady state reflected by model stability.

We progressed to experimentally validate predictions made with the cell-specific model. The MV411 context model captured the signaling impact of the FLT3-ITD to correctly predict induction of apoptosis after treatment with PIM and PI3K inhibitors, and no effect with PI3K inhibitor alone (11). For the cell line KG1A, we identified contribution of high cMyc activity to cell proliferation and correctly predicted insensitivity to inhibition of targets thought to be elevated by the FGFR fusion (28) including AKT, PDK1, and FLT3. The EOL1 context model identified previously unreported combination synergy between PIM and PI3KCA inhibitors, validated through increased tumor growth inhibition. This could lead to patients treated with lower doses of the inhibitors if the same efficacy is achieved by combinations and thereby reducing the risk of toxicity.

Model discrepancies highlight potential gaps in the captured network knowledge and hypotheses that warrant further investigation. For example, our model fails to capture BCR and ERK overexpression following treatment in EOL1 and KG1A cell lines. This cannot be resolved through simple optimization of the current network, suggesting a potential gap in our knowledge of how the MAPK pathway influences these mechanisms (Fig. 3D). We found that Siendones and colleagues (31) had also previously hypothesized the coexistence of transduction signal event, triggering the MAPK pathway independent of the FLT-ITD event, and coupled with poor response to FLT3 inhibitor. Investigating this discrepancy may shed new light on the resistance mechanism of these patients to FLT inhibitors.

Furthermore, using the MOLM16 context model, we were able to systematically explore genetic changes that may render the cell resistant to PIM inhibition. Exome sequencing and subsequent drug combination treatment of MOLM16 cell populations with acquired resistance to AZD1208 confirmed our predicted mechanistic dependency on AKT signaling and AKT inhibition as a second-line therapy to overcome resistance.

By accurately predicting drug responses and combination synergies, and providing the mechanistic insight on the proteins driving the response, we highlight the ability of simulated models and virtual experimentation to prioritize effective therapies accompanied with associated predictive and dynamic biomarkers. Successful drug combinations could significantly augment therapy options for AML patients by overcoming innate and acquired resistance to drugs. Simulated qualitative models potentially offer a virtual platform to screen, discover and prioritize drug combinations in silico, focusing experimental approaches to validation. Comprehensive genetic diagnosis using targeted exome sequencing is already entering the clinic in major teaching hospitals. When coupled with emerging mass cytometry analysis (PMID: 26095251), all the biological information to build patient specific qualitative networks models may soon be available from frontline diagnostics data.

Taken together, the complexity of signaling pathways and the large number of resistance mechanisms mean that executable cellular models that are easily and quickly interpretable, like the ones we have presented here, are key for pinpointing potential combination therapies for different cancer types and subtypes. Furthermore, scaling these executable models to simulate patient-specific cancers paves the way for improved personalized treatments and enhanced precision medicine choices.

S. Grosskurth is a senior scientist/bioinformatics at AstraZeneca and a senior scientist III/bioinformatics at Abbvie. No potential conflicts of interest were disclosed by the other authors.

Conception and design: D. Silverbush, S. Grosskurth, D. Wang, F. Powell, J. Dry, J. Fisher

Development of methodology: D. Silverbush, S. Grosskurth, D. Wang, J. Fisher

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Dry

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Silverbush, S. Grosskurth, J. Dry, J. Fisher

Writing, review, and/or revision of the manuscript: D. Silverbush, S. Grosskurth, D. Wang, F. Powell, B. Gottgens, J. Dry, J. Fisher

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): F. Powell

Study supervision: J. Dry, J. Fisher

We would like to thank Dennis Huszar and Kirsten McEachern for their knowledge surrounding the PIM signaling network, Mika Ahdesmaki for performance of the DNA-seq alignment and variant calling, and Greg O'Connor for generation of the AZD1208-resistant MOLM16 cells. A central element of this study relied on single agents and combination cell line treatment; in particular, we would like to thank Suping Wang and Erica Keaton for AZD1208 combination screen data, Keith Dillman for sample preparation and Western validation for FLT3 combinations and resistant cell lines; all from AstraZeneca. We would also like to thank Bloodwise for supporting B. Gottgens, and the Israeli Ministry of Science, Technology and Space and Edmond J. Safra Center for Bioinformatics at Tel-Aviv University for supporting D. Silverbush.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.