In this study, we set out to construct cell type‐ and patient‐specific models from the drug response data obtained using the MPS technology (overview of the pipeline in Fig 1 A). We took advantage of our tool CellNOptR (Terfve et al , 2012 ), to train a general network of the underlying intrinsic and extrinsic apoptosis pathways from data obtained for two cell lines and biopsies from four pancreatic tumor patients at different stages. We obtained cell line‐ and patient‐specific continuous logic models based on ordinary differential equations (ODEs). We found these models to be a useful tool to understand specific pathway deregulations and to predict new patient‐specific therapies.

Application of this strategy has been limited so far to in vitro contexts, as the experimental technologies to generate perturbation data require large amounts of material, which are unavailable from most primary tissues such as solid tumors. With recently developed organoid technologies, it became possible to generate large amounts of material ex vivo , enabling such screens in principle. However, they would be associated with large costs and, while recapitulating some of the features of the tumor physiology, the cells unavoidably diverge from the primary tumor as they are grown ex vivo (Letai, 2017 ). We have recently developed a novel strategy based on microfluidics that enables testing apoptosis induction upon a good number of conditions (56 with the current settings, with at least 20 replicates each) starting from as little as one million viable cells. Cells are encapsulated in 0.5 μl plugs together with an apoptosis assay and single or combined drugs. Using valves to control individual fluid inlets allows the automatic generation of plugs with different composition. These Microfluidics Perturbation Screenings (MPS) are suitable to collect such drug response datasets even with the very limited number of cells available from tumor resection biopsies (Eduati et al , 2018 ).

Charting the dynamic wiring of signaling networks is of paramount importance to understand how cells respond to their environment. Identifying the differences in this wiring between normal and cancerous cells can shed light on the pathophysiology of tumors and pave the way for novel therapies (Werner et al , 2014 ; Saez‐Rodriguez et al , 2015 ; Zañudo et al , 2018 ). A powerful tool to gain insight into these processes is to monitor the response of cells to multiple perturbations. When combined with mathematical modeling, such data can be used to determine cell type‐specific wiring phenomena, predict efficacy of drug treatments, and understand resistance mechanisms (Saez‐Rodriguez et al , 2011 ; Klinger et al , 2013 ; Merkle et al , 2016 ; Eduati et al , 2017 ; Hill et al , 2017 ).

Results

Data and modeling of apoptosis pathways Experimental data were generated using a novel Microfluidics Perturbation Screening (MPS) platform, which allowed performing combinatorial drug screening of biopsies from human tumors (Eduati et al, 2018; see Materials and Methods). Data represent caspase‐3 (Cas3 in Fig 1B, marked in blue) activation after perturbation with 10 different compounds including seven kinase inhibitors (targeting IKKs, MEK, JAK, PI3K, EGFR, AKT, PDPK1—inhibited nodes are depicted in red in Fig 1B), one cytokine (TNF, stimulated node, in green), and two chemotherapeutic drugs (gemcitabine and Oxaliplatin). All 10 compounds were tested alone and in all 45 possible pairwise combinations (Fig 1C) on two pancreatic cancer cell lines (AsCP1 and BxPC3) and biopsies from four patients with pancreatic tumors at different stages (one intraepithelial neoplasia, two primary tumors, and one liver metastasis). Measurements were performed 16 h after perturbation, when Cas3 activation was shown to reach a plateau (Eduati et al, 2018). To investigate the signaling mechanisms behind the differential drug responses of our cell lines and patients, we derived a general logic model of apoptosis pathways involved in the regulation of Cas3 (our measurement), which is considered as effector node and indicator of apoptosis. Models were then trained using the patient‐specific experimental data from Eduati et al (2018) to obtain personalized models. The general model (Fig 1B) was built integrating information derived from literature and from public repositories (details in Materials and Methods section). The model describes both intrinsic (mediated by the mitochondria, named Mito in the model) and extrinsic (mediated by tumor necrosis factor receptors, TNFRs) apoptotic signals, including nodes encoding for both anti‐ and pro‐apoptotic effects. We incorporated in the model all nodes perturbed by specific compounds in our screening such as targeted drugs (kinase‐specific inhibitors) and the cytokine TNF. The effect of chemotherapeutic DNA damaging drugs was not included in the model since they inhibit DNA replication rather than acting directly on specific signaling nodes. However, nodes such as p53, which are activated by DNA damaging drugs, are included in the model since they are key elements of different pathways. Since our screening included two AKT inhibitors (i.e., MK‐2206 and PHT‐427) with different mechanisms of action (allosteric and PH domain inhibitors, respectively), they were modeled as acting on two different nodes (AktM and AktP, respectively), both needed for the activation of AKT. The logic model includes AND gates (circles in Fig 1B) when all upstream regulators are needed to activate a node, while cases with multiple independent regulators are considered as OR gates. The logic model is interpreted using the logic‐based ordinary differential equation formalism (logic ODEs; Wittmann et al, 2009) as implemented in CellNOptR (Terfve et al, 2012). This formalism allows to maintain the simple causal structure of logic models, while considering also the dynamic nature of the interactions and the continuous scale for the activation of the nodes, by using ODEs. Thus, they are not limited to capture only binary events as is the case for Boolean models. As previously described (Eduati et al, 2017), we consider one parameter for each edge j → i in the network, which characterize the strength of the regulation of species i dependent on species j and one parameter for each node i, which represents the responsiveness of the node (see Materials and Methods).

Calibration of the apoptosis model for cell lines The parameters of the generic model were fitted separately to the data of each cell line, resulting in specific models tailored to the experimental data for each cell line (more details in Materials and Methods section). Parameter fitting was repeated 10 times, and performances were assessed using different metrics to compare model simulations with the experimental data, showing good and quite robust performances (average metrics for AsPC1 and BxPC3, respectively: Pearson correlation 0.72, 0.74; mean squared error 0.03, 0.02; coefficient of determination 0.5, 0.5; Appendix Fig S1A). Model simulations for the best specific models were compared with the corresponding measured experimental data, showing a very good agreement (Pearson correlation equal to 0.89 and 0.83 for AsPC1 and BxPC3, respectively; Appendix Fig S1B and C). The calibrated models for these two cell lines were then used to uncover potential differentially regulated mechanisms which are behind the different drug responses of the cell lines. Due to the limited number of data and the complex nature of the signaling pathways involved in the activation of apoptosis, not all model parameters can be estimated with the same confidence. In order to estimate the variability of the optimized parameter values, we derived a bootstrapped distribution for each parameter for each cell line, by repeating the optimization 500 times while randomly resampling the data with replacement (Dataset EV1). Results from bootstrap iterations were used to assess whether any specific experiment was essential for constraining the model (Appendix Fig S2) and for making predictions (Appendix Fig S3), by considering the left out experiments (due to resampling with replacement) as validation set. This analysis confirmed that even if a specific condition is missing from the training set, this does not significantly affect the model and the resulting predictions. This is probably due to the fact that individual drugs are used in multiple experiments (all the combinations); therefore, even if one specific condition is missing from the training set, the model is still constrained by all the other conditions. These bootstrap distributions were then used to compare the two cell lines using statistical tests to highlight significant differences (Wilcoxon sum rank test, adjusted P‐value < 0.01, effect size > 0.2, see Materials and Methods) as represented in Fig 2A. The comparison revealed some regulatory mechanisms which are upregulated in either AsPC1 or BxPC3. Additionally, the dynamic of Cas3 activation appears to be faster in BxPC3 (node border in purple). Main differences involve the PI3K‐Akt pathway. The main linear pathway is more active in BxPC3, whereby the negative feedback loop, from p53 to PIP3 mediated by PTEN, is stronger in AsPC1. These differences in the model parameters cause changes in the dynamic behavior of the system (Fig 2B) and are behind the differential activation of Cas3 in response to drugs. Figure 2.Models, predictions, and validation for AsPC1 and BxPC3 cell lines Differentially regulated mechanisms in AsPC1 or BxPC3, highlighted in green or purple depending on whether the corresponding estimated parameters are higher in AsPC1 or BxPC3, respectively. Time course simulation of PI3K‐AKT pathway and related inhibitors in AsPC1 and BxPC3. Lines represent median values, and error bars represent standard deviation from the bootstrapped simulations. Assessment of model predictions by cross‐validation, removing for each repetition all experiments involving one of the eight drugs (the corresponding drug targets are reported in the legend) from the training set and using it as test set. Bad predictions (Pearson Correlation < 0.6—which corresponds to removal of drug targeting PDPK1) are marked with empty dots. New drug combinations predicted to be highly specific for each cell line, limited to targettable nodes. Combinations predicted to be specific for AsPC1 are marked in green and those specific for BxPC3 in purple. Combinations predicted to have no effect in both cell lines are marked with a dot. We then investigated if these differences in dynamic behavior could be derived from their genetic makeup (Garcia‐Alonso et al, 2018). From the proteins in our model, only KRAS is functionally mutated in AsPC1 and TP53 in BxPC3. Furthermore, no known direct regulators of the nodes in our model—from those in Omni‐Path (Türei et al, 2016), a compendium of pathway resources—were mutated. Interestingly, KRAS and TP53 are indeed involved in the pathways that we found to be differentially activated between the two cell lines, but on their own, they cannot explain the differences in pathway structure in terms of strength of regulations. Therefore, information on mutations alone would not be sufficient to describe the dynamics of the pathways that mediate apoptosis upon drug treatment. The same holds when looking at basal transcriptomics (Iorio et al, 2016; Appendix Fig S4), supporting the observation that static data are not sufficient to investigate the dynamics of a complex system. Differences in signaling pathways (Fig 2A) are driven by the functional status of the nodes that cannot be inferred solely from differences in gene expression (Appendix Fig S4).

Model predictions and validation We decided to test the predictive power of our cell line‐specific mathematical models in two ways: (i) using cross‐validation on the existing dataset and (ii) predicting the effect of new drug combinations that can be experimentally tested. First, optimization was repeated randomly selecting eight conditions as validation set and using the remaining 36 for training (bootstrapping 100 times). This procedure was repeated for 20 randomized test sets, and results show a good correlation between the predictions and the measurements in the validation set (average Pearson Corr = 0.7). Even when the cross‐validation was repeated removing all the experiments involving a specific drug each time (instead of random ones; Fig 2C), predictions are still very good (Pearson Corr range 0.66–0.97) for all drugs except PHT‐427 (Akt and PDPK1 inhibitor, Pearson Corr = −0.29). This implies that experiments with PHT‐427 are essential to define the models. Mathematical models were then used to simulate the effect of targeting all nodes of the network in pairwise combinations, therefore simulating the effect of potential new drugs acting on pathways that were not previously tested using MPS. For each cell line, we simulated the effect of 186 new perturbations (12 on individual nodes and 174 on node pairs), by inhibiting the corresponding node in the model. Varying confidence of model parameter estimation from the available data is expected to affect the ability to predict certain conditions. By using the family of models optimized using bootstrap, we obtained a distribution of the predicted activation of Cas3 in response to the different simulated conditions, therefore retaining information on the confidence we have for each prediction. In particular, we identified the predictions which were significantly different between the two cell lines (Wilcoxon sum rank test, adjusted P‐value < 0.01, effect size > 0.3, see Materials and Methods; Appendix Fig S5), focusing on those that involve the inhibition of pairwise combinations of the 12 targettable nodes (JAK, MEK, RAS, NFkB, JNK, TNF, AktM, Mdm2, EGFR, PI3K, BclX, IKKs; Fig 2D). This results in 66 drug combinations, 34 of which are specific for AsPC1 and 13 for BxPC3 and 10 have no effect on both cell lines. We then tested experimentally three of the top combinatorial therapies predicted to be highly specific for one of the cell lines based on our mathematical models (Fig 3A–C). For all three cases, we had concordance between model prediction and experimental validation. Interestingly, all three combinations are targeting the PI3K/Akt pathway, which is the one where most differences were highlighted between the two cell lines. As the KRAS oncogene is mutated in > 95% of pancreatic adenocarcinomas (Jones et al, 2008), inhibition of signaling pathway downstream of KRAS is considered to be an attractive therapeutic approach. Against this background, we tested the combination of the Akt inhibitor MK‐2206 with a MEK inhibitor (trametinib). As predicted by our model, the validation experiments showed that the combination of trametinib and MK‐2206 is more efficacious and synergistic (based on Bliss independence model) in BxPC3 than in AsPC1 cells (Fig 3D, Appendix Fig S6). In line with this, a recent report on preclinical models of pancreatic cancer displayed enhanced efficacy of gemcitabine plus nabPaclitaxel when combined with MK‐2206 and trametinib (Awasthi et al, 2019). KRAS‐mutant colorectal cancer was found to be extremely sensitive to combined inhibition of Bcl‐2/Bcl‐xL and mTORC1/2 (Faber et al, 2014). Interestingly, the combination of navitoclax (a Bcl‐2/Bcl‐xL/Bcl‐W antagonist) with PHT‐427 (targeting PDPK1 and Akt, the bona fide downstream effector of mTORC2) was predicted by our model to be specific for BxPC3 cells. Additionally, loss of function in PTEN was shown to be important for synergistic interaction between MEK and mTOR inhibitors (Milella et al, 2017). While both AsPC1 and BxPC3 are PTEN wild type, our model identified a weaker negative feedback loop mediated by PTEN in the BxPC3 cell line that could justify the observed synergy. This prediction was also confirmed by the experimental data showing that the combination of navitoclax and PHT‐427 is more efficacious and synergistic in BxPC3 than in AsPC1 cells (Fig 3E, Appendix Fig S7). Of note, the combination of Bcl‐2/Bcl‐xL/Bcl‐W antagonists and drugs targeting PDPK1/Akt has thus far not been discussed as a treatment option for pancreatic cancer, further demonstrating the applicability of our model to identify innovative combinatorial approaches. Finally, we tested the combination of a PI3K inhibitor (taselisib) with the Bcl‐2/Bcl‐xL/Bcl‐W antagonist navitoclax. Agents targeting the PI3K pathway in combination with a Bcl‐2 family inhibitor have been previously suggested to be relevant in the context of pancreatic cancer (Tan et al, 2013). Hence, being able to predict the efficacy of this combination for specific patients (or cell lines in this case) would be highly desirable. Our validation experiments showed that the combination of taselisib and navitoclax is more efficacious and synergistic in AsPC1 than in BxPC3 cells (Fig 3F, Appendix Fig S8), confirming our model‐based predictions. Figure 3.In vitro and in vivo experimental validation of model predictions A–C. Model simulations when inhibiting (A) MEK and AktM nodes, (B) BclX, AktP and PDPK1 nodes, (C) BclX and PI3K nodes. Data are shown using notched boxplots: the middle line represents the median, the box limits correspond to the interquartile range and the whiskers extend to the most extreme data point, which is no more than 1.5 times the length of the box away from the box (outliers are represented as dots).

D–F. In vitro experimental validation of the combination of (D) trametinib (MEK inhibitor, anchor drug at 1 μM) and MK‐2206 (Akt inhibitor, 8‐points 1:3 dilution series), (E) navitoclax (BclX inhibitor, anchor drug at 10 μM) and PHT‐427 (AktP and PDPK1 inhibitor, 8‐points 1:2 dilution series), (F) navitoclax (BclX inhibitor, anchor drug at 2.5 μM) and taselisib (PI3K inhibitor, 8 time points 1:3 dilution series). Data shown are for three biological replicates with three technical replicates each (error bars represent standard error of the technical replicates). Corresponding boxplots show the resulting synergy scores (Bliss model) computed for each biological replicate considering all concentrations of the anchor drug and the highest two concentrations of the combined drug. Summary statistics are represented using a horizontal line for the median and a box for the interquartile range. The whiskers extend to the most extreme data point, which is no more than 1.5 times the length of the box away from the box.

G. In vivo validation on cell line derived xenograft mouse models comparing the effect on the mouse models derived from the two cell lines. Data shown are for four mice (error bars represent the standard error—full data are provided as Dataset EV2 and EV3). P ‐values were derived using linear mixed‐effect models to compare longitudinal data (corrected for multiple comparisons using Holm method).

H. In vivo comparison of the combination of trametinib and MK‐2206 with the control condition (vehicle alone) and with the standard of care gemcitabine (error bars represent the standard error—full data are provided as Dataset EV2 and EV3). Encouraged by these in vitro results, we performed in vivo validation on xenograft mouse models (Dataset EV2 and EV3). As predicted by our mathematical models, when treated with a combination of trametinib and MK‐2206, both mouse models showed a significantly different response (P‐value = 0.04, see Materials and Methods for details) with BxPC3‐derived mice showing a stronger response (Fig 3G). For the other treatments, no significant difference was found (P‐value = 0.27 for navitoclax + taselisib, P‐value = 0.22 for navitoclax + PHT‐427). Additionally, we tested the effect of gemcitabine, the standard of care in pancreatic cancer, which showed a similar effect in both out mouse models (P‐value = 0.75). Interestingly, the combination of trametinib and MK‐2206 provides significantly better treatment with respect to gemcitabine in the BxPC3 mouse model (P‐value = 0.04), while this is not the case for the AsPC1 mouse model (P‐value = 0.77; Fig 3H). Overall, these results suggest that our mathematical models have a potential for predicting personalized treatment that can specifically improve in vivo response with respect to the standard of care.

Personalized apoptosis models for patients’ tumors The same fitting pipeline previously described for the cell lines was applied to the data from the four pancreatic tumor biopsies (intraepithelial neoplasia, two primary tumors and one liver metastasis) to obtain personalized models. Patient‐specific parameter distributions were used to investigate patient heterogeneity at the level of mechanisms involved in apoptosis signaling pathways. Results are summarized in Fig 4A, showing the 16 (out of 93) parameters which are different in at least one patient (Kruskal–Wallis rank sum test, adjusted P‐value < 0.01, effect size > 0.2, see Materials and Methods). For these parameters, we also performed post hoc pairwise statistical tests to directly compare all patients (Wilcoxon sum rank test, adjusted P‐value < 0.01, effect size > 0.2, see Materials and Methods). For instance, for the parameter representing the EGFR → JAK regulation, the null hypothesis of equal distribution is not rejected when comparing the two primary tumors (lower two boxes in gold) between themself and with respect to intraepithelial neoplasia (top‐left box, half gold next to the corresponding interaction in Fig 3A); however, it is rejected when comparing each of them with liver metastasis (top‐right box, in cyan). Also, the comparison of intraepithelial neoplasia and liver metastasis suggests that the two samples do not come from different distributions (top boxes, both cyan). Figure 4.Patient‐specific models of signaling pathway Mechanisms which are differentially regulated among the patients are highlighted with thick black lines. Colored squares represent the same distribution (same color) or differential distribution (different colors) across the four patient samples. Patient‐specific predictions of new drug combinations. Predictions are ranked for each patient and color coded to compare the efficacy with the other patients. Overall, the most different sample is the liver metastasis (different from all the others in 7 of the 16 heterogeneous parameters (44%)), especially in the extrinsic apoptosis pathway mediated by complex I, cIAPs, and Cas8. This larger dissimilarity could be justified by the difference both in stage and in tissue, since all other samples were resected from the pancreas. Also, the intraepithelial neoplasia shows a quite high level of dissimilarity (37.5%), localized in particular in the IKKs‐NFkB pathway, which could reflect the more advanced stage of the disease. Interestingly, the two primary tumors are the most similar to each other (similar in 11 out of the 16 parameters, corresponding to 69%). However, significant differences were found especially in the PI3K‐AKT pathway, similarly to what we observed for the two pancreatic cancer cell lines. Importantly, these similarities, which reflect the different tumor stages, were not evident directly from the data, where primary tumor #1 clusters closer to the liver metastasis and primary tumor #2 closer to the intraepithelial neoplasia (Appendix Fig S9).