Explainable artificial intelligence for precision medicine in acute myeloid leukemia

Artificial intelligence (AI) can unveil novel personalized treatments based on drug screening and whole-exome sequencing experiments (WES). However, the concept of “black box” in AI limits the potential of this approach to be translated into the clinical practice. In contrast, explainable AI (XAI) focuses on making AI results understandable to humans. Here, we present a novel XAI method -called multi-dimensional module optimization (MOM)- that associates drug screening with genetic events, while guaranteeing that predictions are interpretable and robust. We applied MOM to an acute myeloid leukemia (AML) cohort of 319 ex-vivo tumor samples with 122 screened drugs and WES. MOM returned a therapeutic strategy based on the FLT3, CBFβ-MYH11, and NRAS status, which predicted AML patient response to Quizartinib, Trametinib, Selumetinib, and Crizotinib. We successfully validated the results in three different large-scale screening experiments. We believe that XAI will help healthcare providers and drug regulators better understand AI medical decisions.

In the supplementary material, we have included an index of figures, supplementary methods, and supplementary results.The supplementary methods include a description of the missing values in the data as well as a thorough description of the optimization procedure, and several in-silico validations.
On the other hand, the supplementary results include a description of the standard classification of the patients, an analysis of the biomarkers (prevalence of each of them with the aid of the maftools R package), an analysis of the known drug biomarker associations, and a profile of enriched functions of the four subgroups according to the predicted treatment.This heatmap shows the response for each patient to all possible treatments.Treatments are contained in the rows and patients in the columns.Columns are labelled according to subgroup.The content is the sensitivity score in blue represents that the patient is sensitive to that treatment and in red that the patient is resistant to the treatment or the patient is not having the biomarker associated or it has been treated before.______ 10 Figure S5: Results 10-fold Cross-validation.In the Y-axis is plotted the number of folds in which each subgroup appeared.In the X-axis all the different subgroups that appeared in the cross-validation.___ 11 Figure S6: Results Sensitivity.We plotted the sorted ranks of the drug predicted by MOM for each patient.This plot shows that the suggested treatment was the best one in 3% of the cases, within the top 10% in 30% of the cases, within the first quartile in 46% of cases.The statistical significance for each of the thresholds can be stated using a Bernoulli distribution.We also included the p-values according to this distribution using thresholds for 1% (p-value=0.005),5%(p-value=4.58e-11),10%(p-value=2.24e-23) and 25%(p-value= 4.32e-17) ____________________________________________________________ Figure S7: Validation in cell lines using DEMETER Score.In red the cell lines with the biomarker associated with the treatment; in blue the cell lines without the biomarker.On the left, FLT3 Mut  This plot shows the different types of genetic variants, the lateral barplot shows the sum of all the genetic alterations in all patients and colours the type of variant.In the horizontal axis we have the individual patients' information, showing some patients having up to 70 co-occurring mutations.Of the 608 patients, 538 had at least one genetic variant (88.16%).____________________________________________________________________ Figure S14: Genetic Variant Type Summary.Variant classification plot represents the different genetic variant types, in terms of functionality and we see that the most common is missense variant, the Variant type plot represents the type of structural variant, whether they are Single Nucleotide Variants (SNVs) or Indels with a clear predominance of SNVs, SNV classification plot shoed the type of mutational signature that is predominant with the signature C>T that is quite frequent in malignant cancer types followed by C>A which is associated with environmental exposure 10 .The Variants per Sample Plot, tells that there is a median of 8 variants per patient.The variant classification summary plot summarises all the prior plots and, finally, the Top 10 mutated genes plot, shows for each of the top 10 genes the type of mutation.__________________________________________________________________________

Missing Data Processing
We used KNN Impute 1 to impute the missing values.No values were imputed for patients with NA % higher than 80% and drugs whose proportion of missing values was above 70%.These thresholds seem to be high, however, the distribution of the NAs shows that the vast majority of either the patients or the drugs have less than 20 % of missing values (Figure S1).Since the proportion of missing values was reasonably low we did not compare different algorithms to impute them.We used KNNImpute since it has shown to be reliable and widely used.

MILP optimization model
The core of MOM is an integer programming optimization model that predicts the combination of drugs and biomarkers that optimize patient response to treatment (i.e., IC50*).Let us define a treatment as a combination of a drug and a companion biomarker.The solution to the optimization problem consists of a set of treatments that will be applied sequentially to patients in a defined number of steps (one treatment per step).
Let S, P, and T be the total number of possible steps, patients, and treatments included in the study respectively.Let OM be an input    matrix of essentiality, whose elements  , contains the normalized sensitivity score (IC50*) of ex-vivo experiments for a patient p and a treatment t, which fulfills  , ≤ 0. IC50* values are all negative.
Let K be a    a binary matrix whose element   denotes whether a patient p is eligible for the treatment t, i.e., the treatment's companion biomarker is present in the patient, as follows: Let  be a binary      array whose element   states whether a patient p is treated with treatment t in step s, as follows: Let  be a    binary matrix whose element   represents whether a treatment t is used in step s, as follows: Given these variables, the MOM algorithm was built as a MILP optimization problem defined by the following equations. ( The objective function of the MILP problem, equation (1), is to minimize drug sensitivity scores (IC50 * ) for all patients.The sensitivity score will be considered if it is included in the problem solution (  ).As drug sensitivity scores are negative, the MILP solution will intrinsically maximize the number of treated patients, as each included patient adds a negative term to the objective solution.
The proposed MILP problem has four sets of restrictions, namely equation (2) to equation (5).Equation ( 2) is a set of S restrictions stating that each step consists of one treatment.Equation (3) is a set of P restrictions stating that at most one treatment must be used to treat each patient.Equation ( 4) is a set of      restrictions stating that the treatment t can be applied to patient p in step s, only if (i) the patient p is eligible for treatment t based on his/her biomarkers (  = 1), and (ii) the treatment t is used in step s (  = 1).Finally, equation ( 5) is set of      restrictions that impose that the treatments included in the solution must be selected hierarchically, i.e., if we have a patient that would be eligible for two treatments, only the first treatment must be considered in the optimal solution.
The optimization model needs several input variables.The most crucial input variable is the matrix containing all the pre-processed essentiality information regarding each of the individual associations, called IC50* p,t .Figure S2 contains a pre-visualization of this OM matrix, located in the rows all the different individual associations called treatments; and in the columns the patients.This matrix contains the sensitivity score for the patients having the biomarkers associated with each treatment in the row.Also, groups labels shown in the heatmap include relevant biomarker status.The score in blue represents sensitivity while in red refers to resistance or absence of the biomarker.In the Y-axis is plotted the number of folds in which each subgroup appeared.In the X-axis all the different subgroups that appeared in the cross-validation.

Performance of MOM
The application of typical performance measures of machine learning (specificity, accuracy, sensitivity, ROC and PR curves, etc.) to this specific problem is not straightforward.Since MOM suggests a single drug for each patient, the potential contingency matrix will be very unbalanced: for each patient, only the drug suggested by MOM is a positive and all the other treatments are negatives.If the suggested drug is the one with the lowest IC50*, the prediction will be a true positive.Otherwise, it will be a false positive.On the other hand, all the drugs that were not selected are true negatives (except the one with the lowest IC50*).The drug with the lowest IC50*, if not selected by MOM, will be a false negative.This approach to construct the contingency matrix has several drawbacks.First, the classes are strongly unbalanced -for each patient there is only one positive and hundreds of negatives for our data.Furthermore, this approach makes no difference between predicting the second-best drug or predicting the worst drug.
Instead of this approach, we plotted (Figure S6) the sorted ranks of the drug predicted by MOM for each patient.This plot shows that the suggested treatment was the best one in 2% of the cases, within the top 10% in 30% of the cases and so on.The statistical significance for each of the thresholds can be stated using a Bernoulli distribution.We computed the p-values according to this distribution using thresholds for 1%, 5%, 10% (0.005, 4.58e-11, and 2.24e-23 respectively).A "prediction" algorithm that prescribes a drug by chance would show a curve close to the diagonal in this graph.

Patients
Ranking Percentage

Standard Classification of Patients
Current patient stratification guides divide AML patients into three subgroups according to their prognosis, namely favorable-, intermediate-and adverse-risk.Each subgroup is defined by a combination of genetic biomarkers that can be either chromosomal rearrangements, genetic mutations, or allele deletions.Thus, the favorable risk subgroup -a 5-year overall survival (OS) of 45% to 80%-includes 45% of AML patients and is diagnosed mainly through the biomarkers NPM1 Mut , chromosome 16 inversion (inv( 16)), and CEBPA Mut .The intermediate-risk subgroup -5 years OS of 30%-comprises 25% of AML cases and is associated with the internal tandem duplications in the FLT3 gene (FLT3-ITD), and NPM1 WT .Finally, the adverse risk subgroup -5 years OS of 10%-represents 30% of AML cases and has scattered deletions and complex karyotypes as biomarkers 6 .
Despite having disparate prognoses, the treatment is similar: all of them are treated with standard induction cytotoxic therapy ("3+7") 6 with different dosages and aggressiveness depending on the severity.Recently, FLT3 inhibitors have been incorporated as a treatment directed to FLT3-ITD patients, but effective treatments for patients who do not have this alteration remain an unmet clinical challenge 7 .The lack of treatments directed to FLT3 WT patients -70% of all AML patients 7 -demands a new classification guide based on the response to therapy based on the presence or absence of specific genetic biomarkers.

BeatAML Treatment Stages
From the 319 patients contained in the study coming from the BeatAML cohort, the chemotherapy phase veries considerabily from the 319, 52 patients had not treatment stage assigned from which 32 were de Novo AML.The latter Treatment Stage is shown in the following graph (Figure S12a).There are patients who have gone thorugh up to seven different therapy stages.The distribution of the patients accounting for the different chemotherapy phasesexcluding the patients with missing info or who did not take any chemotherapy-is the following (Figure S12b).
Even though, there could be different agreements mentioning whether it could be more interesting to treat the disease according to the clinical stage: de Novo, Relapsed, Refractory, … Our inclusion criterion focused on the disease itself, as the cohort is an adult cohort with median of age 61, the results of the manuscript reflect an approach to treat the adult AML disease, independently from the chemotherapy phase.

Biomarker Analysis
We performed an extensive biomarker analysis to characterize the WES from Beat AML 8 cohort using the maftools 9 R package (version 2.10.05).We intended to understand the different processes that are regulating the biomolecular characterization of this cohort.To do it, we plotted several figures that contain the relevant information concerning the genetic profile of the patients in the cohort.Figure S13 shows that 88.16% of patients have at least one genetic variation, the majority of these variants are missense, and DNMT3A, NPM1 and NRAS are the most commonly mutated genes in this cohort.Moreover, from these variants, the majority correspond to single Nucleotide Variants (SNVs) with the signature C>T that is quite frequent in malignant cancer types followed by C>A which is associated with environmental exposure 10 (Figure S14).We appreciated that the median of genetic variants per patient is 8 variants and that only FLT3 and SRSF2 had insertions i.e.FLT3-ITD.
We tried to understand more in-depth the SNVs changes and classified them into transitions (two-ring purines or one-ring purines changes) and transversions (changes of purines for pyrimidines) we discovered that in this cohort the transitions are more frequent with the most common transition being C>T, followed by a transversion C>A (Figure S15).
Finally, we addressed the biological consequences of this mutational landscape by interrogating the alteration of the most common oncogenic pathways (Figure S18).We saw that RTK-RAS is the most affected pathway, having an alteration in 31 out of 85 genes and it is present in 237 out of 608 patients.Remarkable alterations include NOTCH, WNT, MYC, TP53 and TGF-β pathways.We included a summary by the patient showing the complete pathway alterations (Figure S19).or Indels with a clear predominance of SNVs, SNV classification plot shoed the type of mutational signature that is predominant with the signature C>T that is quite frequent in malignant cancer types followed by C>A which is associated with environmental exposure 10 .The Variants per Sample Plot, tells that there is a median of 8 variants per patient.The variant classification summary plot summarises all the prior plots and, finally, the Top 10 mutated genes plot, shows for each of the top 10 genes the type of mutation.In the Y-axis all genes are included in the pathway, in blue the oncogenes and red tumour suppressor genes.In the X-axis all the samples with RTK-RAS altered and the red marks show the pathway genes altered for each sample.

Individual Drug Biomarker Associations
Drugs provided by the ex-vivo experiments in the Beat AML cohort respond to the following diseases categories (Figure S20).Once we performed the pre-processing in the gene variants data and drug screenings, we obtained the individual associations as mentioned in the main manuscript.Individual associations show how biomarkers can influence sensitivity or resistance for a certain drug.We discovered a total number of 156 significant associations (Table S1), all of them, candidates for patient stratification.

Individual Analysis of the Predicted Subgroups
We tried to understand more in-depth each of the subgroups whose treatment according to MOM is different.The methodology is described in the Methods section of the main manuscript.This section includes the results from the enrichment analysis based on gene expression including all the functions that appeared as statistically enriched from the two conditions up and downregulated.The color in the plots refers to the statistical significance the redder it is the more significant, the bluer it is, the less significant it is.Figures S22-S29, describe all these conditions showing a bar plot, dot plot, upset plot, and emap plot to better understand the relevance and interaction of the obtained figures.

Figure S1 :
Figure S1: NAs distribution in IC50* data.A) shows the percentage of NAs by patients.Patients with more than 80 % of missing data (dotted blue line) were excluded from the analysis.B) This graph shows the percentage of NAs by drug, drugs with more than 70% (dotted red line) were excluded from the analysis.____________________________________________________________________________ 5 Figure S2: Heatmap showing matrix OM p,t .Located in the rows are all the different individual associations called treatments, and in the columns are the patients.This matrix contains the sensitivity score for the patients having the biomarkers associated with each treatment in the row.Also, groups labels shown in the heatmap include relevant biomarker status.The score in blue represents sensitivity while in red refers to resistance or absence of the biomarker.__________________________________ 8 Figure S3: New Patient Stratification showing the score for each of the patients on each of the subgroups.__________________________________________________________________________ 9 Figure S4: Heatmap of all possible treatments labelled by selected treatment.This heatmap shows the response for each patient to all possible treatments.Treatments are contained in the rows and patients in the columns.Columns are labelled according to subgroup.The content is the sensitivity score in blue represents that the patient is sensitive to that treatment and in red that the patient is resistant to the treatment or the patient is not having the biomarker associated or it has been treated before.______ 10 FigureS5: Results 10-fold Cross-validation.In the Y-axis is plotted the number of folds in which each subgroup appeared.In the X-axis all the different subgroups that appeared in the cross-validation.___ 11 FigureS6: Results Sensitivity.We plotted the sorted ranks of the drug predicted by MOM for each patient.This plot shows that the suggested treatment was the best one in 3% of the cases, within the top 10% in 30% of the cases, within the first quartile in 46% of cases.The statistical significance for each of subgroup, on the right NRAS Mut subgroup._________________________________________________________ Figure S8: Overall biomarker effect in the DEMETER cohort.All cell lines with the biomarker for the different subgroups are shown in red.In blue are shown all the cell lines that were not having the associated biomarker.P.value was computed using Wilcoxon's Test.The smaller the DEMETER Score is, the more sensitive the cells are._________________________________________________________ Figure S9: Overall biomarker effect in the DEMETER2 cohort.All cell lines with the biomarker for the different subgroups are shown in red.In blue are shown all the cell lines that were not having the associated biomarker.P.value was computed using Wilcoxon's Test.The smaller the DEMETER2 Score is, the more sensitive the cells are._________________________________________________________ Figure S10: Overall biomarker effect in the CERES cohort.All cell lines with the biomarker for the different subgroups are shown in red.In blue are shown all the cell lines that do not have the associated biomarker.P-value was computed using Wilcoxon's Test.The smaller the CERES Score is, the more sensitive the cells are._________________________________________________________________ Figure S11: Overall biomarker effect in the GDSC cohort.All cell lines with the biomarker for the different subgroups are shown in red.In blue are shown all the cell lines that were not having the associated biomarker.P.value was computed using Wilcoxon's Test.The smaller the IC 50 is, the more sensitive the cells are._________________________________________________________________ Figure S12: Treatment Stages in the BeatAML cohort.A) From the 319 patients contained in the study coming from the BeatAML cohort, the chemotherapy phase veries considerabily.B) Distribution of the patients accounting for the different chemotherapy phases excluding the patients with missing info or who did not take any chemotherapy._____________________________________________________ Figure S13: Mutational Status of Beat AML cohort.

Figure S15 :
Transitions(Ti) and Transversions (Tv) landscape in Beat AML cohort.Ti vs. Tv plot shows the number of the transitions and transversions showing that even when transversions seem most probable to occur, transitions are more present in this cohort.Boxplot showing overall distribution of six different conversions and as a stacked barplot showing the fraction of conversions in each sample.The most common transition is C>T, followed by a transversion of C>A._____________________________Figure S16: Somatic Interactions.

Figure S18 :
Figure S18: Oncogenic Signalling Pathways altered in Beat AML cohort.The barplot on the left represents the proportion of genes that are altered in the pathway, whereas the barplot on the right, shows the proportion of samples that are having an alteration in that pathways._________________ FigureS19: RTK-RAS pathway alterations.In the Y-axis all genes are included in the pathway, in blue the oncogenes and red tumour suppressor genes.In the X-axis all the samples with RTK-RAS altered and the red marks show the pathway genes altered for each sample._________________________________ Figure S20: Beat AML Drugs distribution.Drugs studied in the BeatAML cohort cover a wide variety of different cancers and diseases, from which 24% are designed, experimental, or already prescribed for AML, 16% for other leukemias, 10% are shown concerning multiple myeloma, and 4% for different types of lymphomas.This means that 54% of the drugs have been studied for hematological malignancies._ Figure S21.Relevant Individual Drug-Biomarker Associations.In blue the patients with the biomarker in orange the patients without the biomarker.P-value association is corrected according to the main manuscript.The score shows the differential IC 50 , more negative means more sensitive.____________ Figure S22: FLT3 Mut -Quizartinib Upregulated GO enrichment Functions _______________________ Figure S24: FLT3 WT & Inv(16) -Trametinib GO upregulated enrichment Functions ________________ Figure S25: FLT3 WT & Inv(16) -Trametinib GO downregulated enrichment Functions _____________ Figure S26: FLT3 WT & no Inv(16) & NRAS Mut -Selumetinib GO upregulated enrichment Functions ___ Figure S27: FLT3 WT & no Inv(16) & NRAS Mut -Selumetinib GO downregulated enrichment Functions _ Figure S28: Rest -Crizotinib GO upregulated enrichment Functions ___________________________ Figure S29: Rest -Crizotinib GO downregulated enrichment Functions _________________________

Figure S1 :
Figure S1: NAs distribution in IC50* data.A) shows the percentage of NAs by patients.Patients with more than 80 % of missing data (dotted blue line) were excluded from the analysis.B) This graph shows the percentage of NAs by drug, drugs with more than 70% (dotted red line) were excluded from the analysis.

Figure S2 :
Figure S2: Heatmap showing matrix OM p,t .Located in the rows are all the different individual associations called treatments, and in the columns are the patients.This matrix contains the sensitivity score for the patients having the biomarkers associated with each treatment in the row.Also, groups labels shown in the heatmap include relevant biomarker status.The score in blue represents sensitivity while in red refers to resistance or absence of the biomarker.

Figure S3 :
Figure S3: New Patient Stratification showing the score for each of the patients on each of the subgroups.

Figure S4 :
Figure S4: Heatmap of all possible treatments labelled by selected treatment.This heatmap shows the response for each patient to all possible treatments.Treatments are contained in the rows and patients in the columns.Columns are labelled according to subgroup.The content is the sensitivity score in blue represents that the patient is sensitive to that treatment and in red that the patient is resistant to the treatment or the patient is not having the biomarker associated or it has been treated before.

Figure S5 :
Figure S5: Results 10-fold Cross-validation.In the Y-axis is plotted the number of folds in which each subgroup appeared.In the X-axis all the different subgroups that appeared in the cross-validation.

Figure S7 :
Figure S7: Validation in cell lines using DEMETER Score.In red the cell lines with the biomarker associated with the treatment; in blue the cell lines without the biomarker.On the left, FLT3 Mut subgroup, on the right NRAS Mut subgroup.

Figure S8 :
Figure S8: Overall biomarker effect in the DEMETER cohort.All cell lines with the biomarker for the different subgroups are shown in red.In blue are shown all the cell lines that were not having the associated biomarker.P.value was computed using Wilcoxon's Test.The smaller the DEMETER Score is, the more sensitive the cells are.

Figure S9 :
Figure S9: Overall biomarker effect in the DEMETER2 cohort.All cell lines with the biomarker for the different subgroups are shown in red.In blue are shown all the cell lines that were not having the associated biomarker.P.value was computed using Wilcoxon's Test.The smaller the DEMETER2 Score is, the more sensitive the cells are.

Figure S10 :Figure S11 :
Figure S10: Overall biomarker effect in the CERES cohort.All cell lines with the biomarker for the different subgroups are shown in red.In blue are shown all the cell lines that do not have the associated biomarker.P-value was computed using Wilcoxon's Test.The smaller the CERES Score is, the more sensitive the cells are.

Figure S12 :
Figure S12: Treatment Stages in the BeatAML cohort.A)From the 319 patients contained in the study coming from the BeatAML cohort, the chemotherapy phase veries considerabily.B) Distribution of the patients accounting for the different chemotherapy phases excluding the patients with missing info or who did not take any chemotherapy.

Figure S13 :
Figure S13: Mutational Status of Beat AML cohort.This plot shows the different types of genetic variants, the lateral barplot shows the sum of all the genetic alterations in all patients and colours the type of variant.In the horizontal axis we have the individual patients' information, showing some patients having up to 70 co-occurring mutations.Of the 608 patients, 538 had at least one genetic variant (88.16%).

Figure S14 :
Figure S14: Genetic Variant Type Summary.Variant classification plot represents the different genetic variant types, in terms of functionality and we see that the most common is missense variant, the Variant type plot represents the type of structural variant, whether they are Single Nucleotide Variants (SNVs) or Indels with a clear predominance of SNVs, SNV classification plot shoed the type of mutational signature that is predominant with the signature C>T that is quite frequent in malignant cancer types followed by C>A which is associated with environmental exposure10 .The Variants per Sample Plot, tells that there is a median of 8 variants per patient.The variant classification summary plot summarises all the prior plots and, finally, the Top 10 mutated genes plot, shows for each of the top 10 genes the type of mutation.

Figure S15 :
Figure S15: Transitions(Ti) and Transversions (Tv) landscape in Beat AML cohort.Ti vs. Tv plot shows the number of the transitions and transversions showing that even when transversions seem most probable to occur, transitions are more present in this cohort.Boxplot showing overall distribution of six different conversions and as a stacked barplot showing the fraction of conversions in each sample.The most common transition is C>T, followed by a transversion of C>A.

Figure S16 :
Figure S16: Somatic Interactions.Mutually exclusive or co-occurring set of genes calculated using pair-wise Fisher's Exact test.Associations plotted in green represent Co-occurrence while brown is a sign of mutual exclusivity.Stars are assigned to associations with P< 0.05.We appreciated that FLT3 and NPM1 variants are co-occurrent, and FLT3 and TP53 and NRAS and IDH2 are mutually exclusive respectively.

Figure S19 :
Figure S19: RTK-RAS pathway alterations.In the Y-axis all genes are included in the pathway, in blue the oncogenes and red tumour suppressor genes.In the X-axis all the samples with RTK-RAS altered and the red marks show the pathway genes altered for each sample.

Figure S20 :
Figure S20: Beat AML Drugs distribution.Drugs studied in the BeatAML cohort cover a wide variety of different cancers and diseases, from which 24% are designed, experimental, or already prescribed for AML, 16% for other leukemias, 10% are shown concerning multiple myeloma, and 4% for different types of lymphomas.This means that 54% of the drugs have been studied for hematological malignancies.

Figure S21 .
Figure S21.Relevant Individual Drug-Biomarker Associations.In blue the patients with the biomarker in orange the patients without the biomarker.P-value association is corrected according to the main manuscript.The score shows the differential IC 50 , more negative means more sensitive.

Figure S17: Translocations and SNVs. All translocations
the Whole Exome Sequencing (WES) data provided in the Beat AML cohort we analysed co-occurrence and mutual exclusivity of genetic variants at the gene level.We appreciated that FLT3 and NPM1 variants are co-occurrent (p-value <0.05), and FLT3 and TP53 (p-value <0.05) and NRAS and IDH2 are mutually exclusive respectively (p-value <0.05) (FigureS16).