Patient-specific Boolean models of signalling networks guide personalised treatments

Abstract

Prostate cancer is the second most occurring cancer in men worldwide. To better understand the mechanisms of tumorigenesis and possible treatment responses, we developed a mathematical model of prostate cancer which considers the major signalling pathways known to be deregulated. We personalised this Boolean model to molecular data to reflect the heterogeneity and specific response to perturbations of cancer patients. A total of 488 prostate samples were used to build patient-specific models and compared to available clinical data. Additionally, eight prostate cell line-specific models were built to validate our approach with dose-response data of several drugs. The effects of single and combined drugs were tested in these models under different growth conditions. We identified 15 actionable points of interventions in one cell line-specific model whose inactivation hinders tumorigenesis. To validate these results, we tested nine small molecule inhibitors of five of those putative targets and found a dose-dependent effect on four of them, notably those targeting HSP90 and PI3K. These results highlight the predictive power of our personalised Boolean models and illustrate how they can be used for precision oncology.

Keywords: computational biology; drug combinations; human; logical modelling; personalised drug; personalised medicine; prostate cancer; simulations; systems biology.

Figure 1.. Workflow to build patient-specific Boolean models and to uncover personalised drug treatments from present work.

We gathered data from Fumiã and Martins, 2013 Boolean model, Omnipath (Türei et al., 2021) and pathways identified with ROMA (Martignetti et al., 2016) on the TCGA data to build a prostate-specific prior knowledge network. This network was manually converted into a prostate Boolean model that could be stochastically simulated using MaBoSS (Stoll et al., 2017) and tailored to different TCGA and GDSC datasets using our PROFILE tool to have personalised Boolean models. Then, we studied all the possible single and double mutants on these tailored models using our logical pipeline of tools (Montagud et al., 2019). Using these personalised models and our PROFILE_v2 tool presented in this work, we obtained tailored drug simulations and drug treatments for 488 TCGA patients and eight prostate cell lines. Lastly, we performed drug-dose experiments on a shortlist of candidate drugs that were particularly interesting in the LNCaP prostate cell line. Created with BioRender.com.

Figure 2.. Prostate Boolean model used in present work.

Nodes (ellipses) represent biological entities, and arcs are positive (green) or negative (red) influences of one entity on another one. Orange rectangles correspond to inputs (from left to right: Epithelial Growth Factor (EGF), Fibroblast Growth Factor (FGF), Transforming Growth Factor beta (TGFbeta), Nutrients, Hypoxia, Acidosis, Androgen, fused_event, Tumour Necrosis Factor alpha (TNFalpha), SPOP, Carcinogen) and dark blue rectangles to outputs that represent biological phenotypes (from left to right: Proliferation, Migration, Invasion, Metastasis, Apoptosis, DNA_repair), the read-outs of the model. This network is available to be inspected as a Cytoscape file in the Supplementary file 1.

Figure 3.. Prostate Boolean model MaBoSS simulations.

(A) The model was simulated with all initial inputs set to 0 and all other variables random. All phenotypes are 0 at the end of the simulations, which should be understood as a quiescent state, where neither proliferation nor apoptosis is active. (B) The model was simulated with growth factors (EGF and FGF), Nutrients and Androgen ON. (C) The model was simulated with Carcinogen, Androgen, TNFalpha, Acidosis, and Hypoxia ON.

Figure 4.. Associations between simulations and Gleason grades (GG).

(A) Centroids of the Principal Component Analysis of the samples according to their Gleason grades (GG). The personalisation recipe used was mutations and copy number alterations (CNA) as discrete data and RNAseq as continuous data. Density plots of Proliferation (B) and Apoptosis (C) scores according to GG; each vignette corresponds to a specific sub-cohort with a given GG. Kruskal-Wallis rank sum test across GG is significant for Proliferation (p-value = 0.00207) and Apoptosis (p-value = 2.83E-6).

Figure 5.. Phenotype score variations and synergy upon combined ERK and MYC_MAX (A and C) and HSPs and PI3K (B and D) inhibition under EGF growth condition.

Proliferation score variation (A) and Bliss Independence synergy score (C) with increased node activation of nodes ERK and MYC_MAX. Proliferation score variation (B) and Bliss Independence synergy score (D) with increased node activation of nodes HSPs and PI3K. Bliss Independence synergy score <1 is characteristic of drug synergy, grey colour means one of the drugs is absent and thus no synergy score is available.

Figure 6.. Model-targeting drugs’ sensitivities across prostate cell lines.

GDSC z-score was obtained for all the drugs targeting genes included in the model for all the prostate cell lines in GDSC. Negative values mean that the cell line is more sensitive to the drug. Drugs included in Table 1 were highlighted. ‘Other targets’ are drugs targeting model-related genes that are not part of Table 1.

Figure 7.. Cell viability assay determined by the fluorescent resazurin after a 48 hours incubation showed a dose-dependent response to different inhibitors.

(A) Cell viability assay of LNCaP cell line response to 17-DMAG HSP90 inhibitor. (B) Cell viability assay of LNCaP cell line response to PI-103 PI3K/AKT pathway inhibitor. (C) Cell viability assay of LNCaP cell line response to NMS-E973 HSP90 inhibitor. (D) Cell viability assay of LNCaP cell line response to Pictilisib PI3K/AKT pathway inhibitor. Concentrations of drugs were selected to capture their drug-dose response curves. The concentrations for the NMS-E973 are different from the rest as this drug is more potent than the rest (see Materials and methods).

Figure 8.. Hsp90 inhibitors resulted in dose-dependent changes in the LNCaP cell line growth.

(A) Real-time cell electronic sensing (RT-CES) cytotoxicity assay of Hsp90 inhibitor, 17-DMAG, that uses the Cell Index as a measurement of the cell growth rate (see the Materials and methods section). The yellow dotted line represents the 17-DMAG addition. (B) RT-CES cytotoxicity assay of Hsp90 inhibitor, NMS-E973. The yellow dotted line represents the NMS-E973 addition.

Figure 9.. PI3K/AKT pathway inhibition with different PI3K/AKT inhibitors shows the dose-dependent response in LNCaP cell line growth.

(A) Real-time cell electronic sensing (RT-CES) cytotoxicity assay of PI3K/AKT pathway inhibitor, PI-103, that uses the Cell Index as a measurement of the cell growth rate (see the Materials and methods section). The yellow dotted line represents the PI-103 addition. (B) RT-CES cytotoxicity assay of PI3K/AKT pathway inhibitor, Pictilisib. The yellow dotted line represents the Pictilisib addition.

Appendix 1—figure 1.. Mean activities by subgroups for gene modules defined from pathways described in ACSN.

(A) And in Hallmarks' gene sets (B) and that are significantly overdispersed over all samples. Blue indicates low pathway activity, red indicates high pathway activity.

Appendix 1—figure 2.. Signed directed interactions between HSP90AA1 and nodes already taken into account in the model.

Appendix 1—figure 3.. shortest paths found between ERG and TMPRSS2 or NKX3-1 by Pypath: no direct interaction is found.

Appendix 1—figure 4.. Boolean toy model to showcase different examples of Boolean formulas.

Appendix 1—figure 5.. Mean probabilities of the nodes characterising the cyclins and proliferation, with nutrients and growth factors as inputs.

We choose initial states for the nodes involved in the cell cycle that correspond to quiescence (cyclins OFF, cell cycle inhibitors Rb and p21 ON), in order to visualise the order of activation of the cyclins: first Cyclin D, then Cyclin B. The mean probabilities reach asymptotic levels because of the desynchronisation of stochastic trajectories in the population.

Appendix 1—figure 6.. Mean probabilities in simulations of mutated models.

(A) Loss-of-function mutation of FOXA1. (B) Loss-of-function mutation of TP53.

Appendix 1—figure 7.. Mean probabilities in simulations of the model with a multiple simulation: the gene fusion TMPRSS2:ERG and a loss-of-function of NKX3-1.

Appendix 1—figure 8.. Data integration in Boolean models to have personalised Boolean models.

Appendix 1—figure 9.. Distribution of 488 TCGA prostate cancer patients’ samples per Gleason grade.

Appendix 1—figure 10.. Associations between simulations and Gleason grades (GG).

Distribution histograms of Apoptosis scores according to GG in three groups (A) and five groups (B). Columns correspond to different personalisation recipes (see Béal et al., 2019 for more details). We found that across 3-stage GG Kruskal-Wallis rank sum test is significant for Apoptosis under the ‘Mut, CNA, and RNA’ recipe (p-value = 2.83E-6) and significant across 5-stage GG (p-value = 1.88E-5). Additionally, we used Dunn’s test to identify which pairs of groups are statistically different focusing on the 3-stage GG and found that grade High is statistically different from grades Low (Bonferroni’s adjusted p-value = 3.3E-3) and Intermediate (Bonferroni’s adjusted p-value = 9.47E-6).

Appendix 1—figure 11.. Associations between simulations and Gleason grades (GG).

Distribution histograms of DNA_repair scores according to GG in three groups (A) and five groups (B). Columns correspond to different personalisation recipes (see Béal et al., 2019 for more details). Kruskal-Wallis rank sum test across 3-stage GG is neither significant for DNA_Repair under the ‘Mut, CNA and RNA’ recipe (p-value = 0.217) nor across 5-stage GG (p-value = 0.0995).

Appendix 1—figure 12.. Associations between simulations and Gleason grades (GG).

Distribution histograms of Invasion scores according to GG in three groups (A) and five groups (B). Columns correspond to different personalisation recipes (see Béal et al., 2019 for more details). Kruskal-Wallis rank sum test across 3-stage GG is significant for Invasion under the ‘Mut, CNA, and RNA’ recipe (p-value = 0.0358), but not significant across 5-stage GG (p-value = 0.134). Using Dunn’s test on the 3-stage GG, we found that grade High is statistically different from grade Intermediate (Bonferroni’s adjusted p-value = 0.037).

Appendix 1—figure 13.. Associations between simulations and Gleason grades (GG).

Distribution histograms of Migration scores according to GG in three groups (A) and five groups (B). Columns correspond to different personalisation recipes (see Béal et al., 2019 for more details). Kruskal-Wallis rank sum test across 3-stage GG is neither significant for Migration under the ‘Mut, CNA, and RNA’ recipe (p-value = 0.173) nor across 5-stage GG (p-value = 0.275).

Appendix 1—figure 14.. Associations between simulations and Gleason Grades (GG).

Distribution histograms of Proliferation scores according to GG in three groups (A) and five groups (B). Columns correspond to different personalisation recipes (see Béal et al., 2019 for more details). Kruskal-Wallis rank sum test across 3-stage GG is significant for Proliferation under the ‘Mut, CNA, and RNA’ recipe (p-value = 0.00207) and across 5-stage GG (p-value = 0.013). Using Dunn’s test on the 3-stage GG, we found that grade High is statistically different from grade Intermediate (Bonferroni’s adjusted p-value = 0.0023).

Appendix 1—figure 15.. Principal Component Analysis of all 488 TCGA patients in 3 Gleason Grades using the vectors of all five phenotypes from the model.

Appendix 1—figure 16.. Principal Component Analysis of all 488 TCGA patients in 5 Gleason Grades using the vectors of all five phenotypes from the model.

Appendix 1—figure 17.. Principal Component Analysis’ barycenters of all 488 TCGA patients grouped in 3 Gleason Grades using the vectors of all five phenotypes from the model.

This is the same figure as Appendix 1—figure 3 in the main text.

Appendix 1—figure 18.. Principal Component Analysis’ barycenters of all 488 TCGA patients grouped in 5 Gleason Grades using the vectors of all five phenotypes from the model.

Appendix 1—figure 19.. Nodes in the Boolean model that have a Proliferation value of at least 30% less the wild type value upon inactivation.

(A) Nodes from aggregating all patient-specific results; (B) Nodes from patients from Gleason Grade 1; (C) Nodes from patients from Gleason Grade 2; (D) Nodes from patients from Gleason Grade 3; (E) Nodes from patients from Gleason Grade 4; (F) Nodes from patients from Gleason Grade 5.

Appendix 1—figure 20.. Nodes in the Boolean model that promote Apoptosis at least 30% more than the wild type value upon inactivation.

(A) Nodes from aggregating all patient-specific results; (B) Nodes from patients from Gleason Grade 1; (C) Nodes from patients from Gleason Grade 2; (D) Nodes from patients from Gleason Grade 3; (E) Nodes from patients from Gleason Grade 4; (F) Nodes from patients from Gleason Grade 5.

Appendix 1—figure 21.. Phenotype simulation results across GDSC prostate cell line-specific Boolean models’ simulation with random initial conditions.

WT stands for wild type model, the original prostate model with no personalisation.

Appendix 1—figure 22.. Phenotype simulation results across GDSC prostate cell line-specific Boolean models’ simulation with different initial conditions.

WT stands for wild type model, the original prostate model with no personalisation.

Appendix 1—figure 23.. Phenotype simulation results across GDSC prostate cell line-specific Boolean models’ simulation with random initial conditions under different personalisation recipes.

Mutations and CNA are always considered as discrete data and RNA expression is always considered as continuous data.

Appendix 1—figure 24.. Principal Component Analysis (PCA) of the RNA dataset used to tailor the prostate cell lines.

Appendix 1—figure 25.. Results of the ssGSEA performed on the RNA dataset used to tailor the prostate cell lines.

Appendix 1—figure 26.. Phenotype probabilities of LNCaP model under random initial conditions.

Appendix 1—figure 27.. Wild type and LNCaP-specific model phenotype probability variations under four different growth conditions.

Androgen stands for androgen presence, EGF for EGF presence, and 00 for lack of androgen and EGF.

Appendix 1—figure 28.. PCA of the 32,258 single and double LNCaP model mutants with combinations of the most probable phenotypes.

Appendix 1—figure 29.. PCA of the 32,258 single and double LNCaP model mutants with the decomposition in single phenotypes.

Appendix 1—figure 30.. Invasion-Migration-Proliferation phenotype probability distribution across all mutants for logical gates.

Bin where wild type value is found has been marked with dark red colour. (A) Phenotype probability using level one single perturbations; (B) Phenotype probability using level two double perturbations.

Appendix 1—figure 31.. Distribution of perturbations on nodes’ logical gates that reduce Invasion-Migration-Proliferation phenotype probability to zero.

(A) Counts of level one single perturbations; (B) Counts of level two double perturbations.

Appendix 1—figure 32.. Drug sensitivity of the seven prostate cell lines.

Rank normalised drug sensitivity (0: most sensitive; 1: most resistant, based on GDSC AUC drug sensitivity metric) for each GDSC drug across prostate cancer cell lines. Drugs are grouped to be predicted effective drugs based on the LNCaP Boolean model (orange) and predicted ineffective drugs (blue). Mann-Whitney U p-values for differences between the rank normalised drug sensitivity between predicted effective and ineffective drugs: (A) LNCaP, P = 0.00041 (more sensitive to LNCaP model-predicted drugs); (B) 22RV1, P = 0.0033 (more sensitive to LNCaP model-predicted drugs); (C) BPH-1, P = 0.31; (D) DU-145, P = 0.0026 (more resistant to LNCaP model-predicted drugs); (E) PC-3, P = 0.15; (F) PWR-1E, P = 0.075; (G) VCaP P = 0.38.

Appendix 1—figure 33.. Model-targeting drugs’ sensitivities across prostate cell lines.

GDSC Z-score was obtained for all the drugs targeting genes included in the model for all the prostate cell lines in GDSC. LNCaP is highlighted in red, the other seven prostate cell lines in blue and the rest of the GDSC cell lines are coloured in grey.

Appendix 1—figure 34.. Phenotype score variations of the LNCaP model upon nodes’ inhibition under EGF growth condition.

Values of the scores are depicted with a colour gradient.

Appendix 1—figure 35.. Phenotype score variations of the LNCaP model upon nodes inhibition under AR, EGF, 00 and AR_EGF growth conditions.

Values of the scores are depicted with a colour gradient.

Appendix 1—figure 36.. Proliferation phenotype score variations of the LNCaP model upon combined nodes inhibition under EGF growth condition.

Appendix 1—figure 4A is a closer look at ERK and MYC_MAX combination and Appendix 1—figure 4B at HSPs and PI3K combination.

Appendix 1—figure 37.. Apoptosis phenotype score variations of the LNCaP model upon combined nodes inhibition under EGF growth condition.

Appendix 1—figure 38.. Bliss Independence synergies scores variations in Proliferation phenotype of the LNCaP model upon combined nodes inhibition under EGF growth conditions.

Bliss Independence synergy score <1 is characteristic of drug synergy. Appendix 1—figure 4C is a closer look at ERK and MYC_MAX combination and Appendix 1—figure 4D at HSPs and PI3K combination, grey colour means one of the drugs is absent and thus no synergy score is available.

Appendix 1—figure 39.. Bliss Independence synergies scores variations in Apoptosis phenotypes of the LNCaP model upon combined nodes inhibition under EGF growth conditions.

Bliss Independence synergy score <1 is characteristic of drug synergy, grey colour means one of the drugs is absent and thus no synergy score is available.

Appendix 1—figure 40.. Hsp90 inhibitors resulted in dose-dependent changes in the LNCaP cell line growth.

(A) Real-time cell electronic sensing (RT-CES) cytotoxicity assay of Hsp90 inhibitor, 17-DMAG, that uses the Cell Index as a measurement of the cell growth rate (see the Material and Methods section). The yellow dotted line represents 17-DMAG addition. The brown dotted lines are indicative of the cytotoxicity assay results at 24 hours (B), 48 hours (C) and 72 hours (D) after 17-DMAG addition. (E) RT-CES cytotoxicity assay of Hsp90 inhibitor, NMS-E973. The yellow dotted line represents NMS-E973 addition. The brown dotted lines are indicative of the cytotoxicity assay results at 24 hours (F), 48 hours (G) and 72 hours (H) after NMS-E973 addition.

Appendix 1—figure 41.. PI3K/AKT pathway inhibition with different PI3K/AKT inhibitors shows dose-dependent response in LNCaP cell line growth.

(A) Real-time cell electronic sensing (RT-CES) cytotoxicity assay of PI3K/AKT pathway inhibitor, PI-103, that uses the Cell Index as a measurement of the cell growth rate (see the Material and Methods section). The yellow dotted line represents PI-103 addition. The brown dotted lines are indicative of the cytotoxicity assay results at 24 hours (B), 48 hours (C) and 72 hours (D) after PI-103 addition. (E) RT-CES cytotoxicity assay of PI3K/AKT pathway inhibitor, Pictilisib. The yellow dotted line represents Pictilisib addition. The brown dotted lines are indicative of the cytotoxicity assay results at 24 hours (F), 48 hours (G) and 72 hours (H) after Pictilisib addition.