Myeloma cell lines used in this study were U266 (17), JJN3 (18), CAG (19), INA6 (20), OH2 (21, 22), IH1 (23), KJON, VOLIN (24), URVIN, and FOLE. The FOLE MM cell line (Misund et al., unpublished), has biallelically-inactive TP53, and was used as a positive control for testing of the MDM2-P53 interaction inhibitor nutlin-3A. The MM cell lines were cultured in RPMI 1640 medium (ThermoFisher Scientific, Waltham, MA, USA) supplemented with 2mM L-glutamine, 1% Penicillin-Streptomycin, 1x Sodium Pyruvate (NaPur), and fetal bovine serum (FBS) at 10% (JJN3, CAG, INA6, VOLIN, URVIN) or 15% FBS (U266), or 10% human serum (Sigma-Aldrich, Saint-Louis, MO, USA) (IH1, OH2, KJON). The culture medium was supplemented with 2 ng/mL IL6 (ThermoFisher Scientific) for culturing of INA6, KJON, IH1, OH2, URVIN, and VOLIN. The cells were expanded, aliquoted and cryopreserved until experimental assays were performed. See Table 1 for MM cell line characteristics.

Drugs (n=33) were selected based on mutated targets in the cell lines (see Tables 2 , 3 for drugs and their corresponding mutated targets), and were added to 384-well TC-microplates (Greiner #781080) using an acoustic dispenser (Echo 550, Labcyte Inc., CA, USA). Each drug was tested at five concentrations ranging from 1 to 10,000 nM. Experiments on cell lines were done on freshly thawed cells. MM cells (5000 cells per well in 25 ul volume) were transferred into plates using an automatic dispenser (Certus Flex, Fitz Gyger, Thun, Switzerland), and incubated at 37°C for 72h. Cell viability was assessed by the CellTiter-Glo luminescence assay (Promega, WI, USA) according to the manufacturer’s instructions. Luminescence was recorded with an Envision Xcite plate reader (Perkin Elmer, MA, USA). The raw concentration-response data were analyzed using the KNIME software (AG, Zurich, Switzerland) and Rstudio Team (Boston, MA) (25). Normalization of the response readout was done to the negative (0.1% DMSO) and positive (100 µM benzethonium chloride) controls.

Phospho flow assays were performed on freshly thawed myeloma cells as described previously (26). Antibody-stained samples were run on a BD LSR Fortessa and output data were analyzed in Cytobank (https://cellmass.cytobank.org/cytobank/). Raw data were transformed to an arcsinh ratio relative to the signal of an isotype control, which was set to zero.

Curve-fits of normalized concentration-response data used the function drm from the R package drc (https://www.r-project.org/) with the four parameter log-logistic model, LL.4, or the logistic model, L.4, where LL.4 failed to converge. To quantify drug responses, a modified drug sensitivity score (DSS) was calculated for each drug (27). In this modified function, area under the curve was calculated using a response-window from 100% to 10%, and a concentration-window from the minimum concentration tested to the concentration where the viability reached 10% (threshold %). DSS type 1 was used, without the term for division by the logarithm of the upper limit. The DSS scores were calculated on a scale of 0-100. A high DSS therefore indicates that MM cells are drug-sensitive, while a low DSS indicates drug-resistant MM cells.

Statistical analyses of output data from drug sensitivity screens and phospho flow assays were performed in GraphPad Prism v8 (San Diego, CA, USA). Welch’s t test [(*p < 0.05, ***p < 0.001, ns (not significant)] was used to compare two means, as indicated in the respective figure legends. ClustVis (28) was used for unsupervised clustering of the DSS values and column annotations indicating genes with the presence of point mutations. No scaling is applied to rows. Columns were clustered using Euclidean distance and Ward linkage. The significance of the correlation coefficients between DSS values and phospho flow readouts was assessed by Pearson’s r test for Figures 1B, C and Supplementary Figure 2 , where the p-values prior to multiple-testing correction are shown.

Predictive models for drug response (as defined by the DSS) were constructed using either mutation or phospho flow data as covariates. In each case, a type of multi-response penalized linear regression model known as tree-lasso (29, 30) was used for feature selection and to assess the predictive power of the selected features for drug response. Regression coefficients for all 33 drugs were penalized jointly in a hierarchical fashion, with penalties weighted by the strengths of the correlations among drug responses, as defined by the height parameters in the hierarchical tree structure of the DSS data ( Supplementary Figure 1 ). For optimal results, the tree height parameters were cut-off at 0.7, on a scale where the root of the tree is at a height of 1 and the leaves are at a height of zero.

Both the mutation-dependent and protein-dependent models were trained using leave-one-out cross validation to compute the Mean Squared Error (MSE) as a function of the tuning parameter. The tree-lasso cost function was optimized using the sub-function tree.lasso from the R package IPFStructPenalty, available at https://github.com/zhizuio/IPFStructPenalty. In each model, the tuning parameter Λ was set to a value at which the cross-validated MSE becomes essentially flat and the smoothing parameter for Smoothing Proximal Gradient descent (SPG) optimization was set to μ = 10^-4. For the mutation-dependent model, the tuning parameter was set to Λ = 25 (with a cross-validated MSE of 175), and for the protein-dependent model it was set to Λ = 24 (with a cross-validated MSE of 530). The tree-lasso regression coefficients of each model were then computed using the optimal values of the tuning parameter Λ ( Figures 2A , 1A, B ).

The amount of variance explained by each model was also assessed by the R squared value, which was computed for each tuned model as R² = 1 - RSS/TSS, where the Residual Sum of Squares (RSS) and Total Sum of Squares (TSS) were both defined using all 9 cell lines at once. The R squared value was R² = 0.39 for the mutation-dependent model, and R² = 0.58 for the protein-dependent model.

The covariates for the protein-dependent model were defined as the standardized phospho flow data, i.e., for each protein the phospho flow arcsinh ratios were centralized at their mean value across samples and normalized by their standard deviation prior to the tree-lasso model calculations. As a consistency check, an overall matching between the signs of the regression coefficients and those of the protein-drug Pearson correlation coefficients was observed ( Supplementary Figure 1B ), since the standardization of the phospho flow data allows analogous interpretations for these two types of coefficients.

Targeted DNA sequencing was performed using the SureSelect Human Kinome kit (Agilent Technologies), with capture probes targeting 3.2 Mb of the human genome, including exons and untranslated regions (UTRs) of all known kinases and selected cancer‐related genes (n = 612). Paired-end sequencing reads of 100-bp length were aligned to the human reference genome (hg19) with Novoalign (version 2.08.3), followed by filtering and realignment with GATK tools and Picard. Single point mutations were identified using MuTect (version 1.1.4). Variant consequence annotation was performed with ANNOVAR, using RefSeq as the underlying transcript model. Details of the complete variant calling pipeline that was applied on the MM cell lines have been described previously (31). The mean sequencing coverage of the kinome was 477x, and we obtained a minimum coverage of 150x for 87.6% (VOLIN) and 90.0% (KJON-1) of the kinome targeted regions. Variant calling performed on cancer cell lines without utilizing matched normal samples is bound to generate a mix of germline variants and somatic variants. We therefore set up a set of filtering procedures to both i) exclude known germline variants, and ii) enrich for coding, cancer-associated variants. Specifically, we excluded all variants that overlapped with germline variants found in the 1000 Genomes Project phase 3 (minor allele frequency > 1% in any population), and NHLBI Exome Sequencing Project (minor allele frequency > 0.1% in any population) (32, 33). In addition, we excluded variants present in the single nucleotide polymorphism database (dbSNP) (build 138) (34) that had no clinical associations (as given from ClinVar cross-references). Finally, we restricted the variant set to coding variants (missense, stop-gain/stop-loss, frameshift/non-frameshift, splice site donor/acceptor) in known cancer census genes (COSMIC version 68) (35). All variants were subjected to a functional annotation workflow that included UniProt KB (functional protein properties) and Polymorphism Phenotyping v2 (PolyPhen-2) web server (computational predictions of effect of coding variants) (http://genetics.bwh.harvard.edu/pph2/) (36).

In order to identify functionally relevant and actionable mutations in MM, nine MM cell lines ( Table 1 ) were subjected to targeted high-throughput DNA sequencing. Short sequencing reads were processed with a variant calling pipeline and subsequent variant filtering procedure, from which we detected a total of 136 mutated genes in the cell lines ( Supplementary Table 1 ). Mutated genes included NRAS/KRAS, BRAF, RAF-1, TP53, PIK3CA, PIK3R3, MTOR, FLT1, FLT3, FLT4, EGFR, and SYK, which are known to be frequently mutated in MM and other related blood malignancies (37–39), as well as to have a therapeutic potential in MM (39) ( Tables 2 , 3 ).

To prioritize the identified gene variants for functional impact, all protein-coding alterations were subjected to analysis with PolyPhen2 (see Materials and Methods), a tool that provides computational predictions for the functional impact of amino acid changes. We identified 31 amino acid changes with predicted damaging effects which were druggable (39), including the most oncogenic RAS mutated isoform, G12 and Q61, as well as PIK3CA catalytic subunit mutations, p.H1047R and p.W590C (40, 41) ( Tables 2 , 3 and Supplementary Table 1 ).

Based on the mutation analysis, we designed a drug library consisting of 33 clinically relevant drugs targeting identified druggable gene products ( Table 3 ). Drug sensitivity screens were then performed on the 9 MM cell lines.

An unsupervised clustering of the cell lines by their DSS profiles showed that MM cell lines with PI3K and NRAS/KRAS mutations displayed high sensitivity to MEK and PI3K inhibitors, respectively ( Figure 1A ). Concentration-response effects on the MM cell viability to individual inhibitors with and without associated targets are shown in Figures 1B, C . These findings support our hypothesis that the presence of RAS and PI3K mutations may increase MM cell sensitivity to MEK and PI3K inhibitors.

We observed variability in efficacy among drugs targeting the PI3K signaling pathway. The most effective drugs were temsirolimus (mean DSS ± SD = 51.96 ± 13.66), taselisib (mean DSS ± SD = 41.82 ± 17.10), pictilisib (mean DSS ± SD = 32.50 ± 13.66), and copanlisib (mean DSS = 19.83 ± 18.65), while treatment with idelalisib and GSK2636771 resulted in markedly low DSS (mean DSS = 2.29 ± 2.5, mean DSS = 0.26 ± 0.46, respectively) ( Figure 2A ). Among the MEK inhibitors, trametinib induced the highest response (mean DSS ± SD = 47.20 ± 27.46), followed by binimetinib (mean DSS ± SD = 21.88 ± 21.54), cobimetinib (mean DSS ± SD = 21.11 ± 23.80), and U0126 (mean DSS ± SD = 11.83 ± 13.85) ( Figure 2A ).

Interestingly, the MM cell lines CAG and URVIN, harboring PIK3CA and PIK3R3 mutations, respectively, exhibited the highest sensitivity to taselisib and temsirolimus treatment, while JJN3, harboring mutations in PIK3CA and mTOR genes, showed less sensitivity to taselisib compared to CAG and URVIN cell lines. On the other hand, JJN3 gained sensitivity to temsirolimus that also target mTOR ( Figure 2B ). The MM cell lines IH1, OH2, and KJON, with NRAS/KRAS mutations, were found highly sensitive to trametinib ( Figure 2B ). These results indicate that sensitivity to kinase inhibitors may be associated with corresponding pathway mutations.

To study whether the mutational status could predict drug sensitivity, we compared treatment responses to PI3K and MEK inhibitors between cell lines with and without selected mutations. We found that cell lines with PI3K pathway mutations (i.e. PIK3CA, PIK3R3, and mTOR) or RAS gene mutations (i.e. p.G12V, p.G12D, p.G13D, p.Q61K, and p.Q61H) ( Table 2 ) exhibited significantly higher DSS relative to cell lines with WT forms of the listed genes when comparing aggregated data for each class of drugs tested ( Figures 3A, B , left panels). The three PI3K inhibitors, taselisib, copanlisib, and pictilisib, each showed a trend towards association with PI3K mutational status, but this was not statistically significant ( Figure 3A , right panels). Of the four MEK inhibitors tested, a statistically significant association was found between the response to binimetinib and RAS mutational status whereas trends towards statistical significance were observed for trametinib, cobimetinib, and U0126 ( Figure 3B , right panels).

Interestingly, by applying a tree-lasso regression model as described in Materials and Methods, a few mutations were found to be associated in varying degrees with high sensitivity to several drugs, as indicated by the positive values in the sparse matrix of regression coefficients ( Figure 3C ). In particular, the mutational status of the PIK3CA and NRAS genes were selected as predictors for drug sensitivity to PI3K and MEK inhibitors ( Figure 3C ).

We also noticed that, while the prediction analyses demonstrated that the mutational status may be useful to predict drug sensitivity, the predictive impact of NRAS mutations on the drug sensitivity to MEK inhibitors was higher than the one observed for PIK3CA mutations on PI3K inhibitors, as indicated by the color-intensity of the heatmap in Figure 3C .

The inhibition of the MDM2-TP53 interaction by nutlin-3A is an attractive strategy to stabilize the P53 mediated apoptosis in various WT TP53 tumors, including MM (42–44), and therefore worth investigating. Interestingly, we observed differing sensitivities to the MDM2-TP53 inhibitor nutlin-3A between the TP53 WT IH1 and the biallelically TP53 mutant MM cell line FOLE. ( Figure 3D ). When stratifying nutlin-3A responses on TP53 mutational status in all MM cell lines tested, we observed a significantly higher drug sensitivity in the TP53 WT cell lines, in agreement with earlier reports (14, 15, 45) ( Figure 3E ).

Having demonstrated that gene mutations in MM cell lines are linked to drug sensitivities towards PI3K and MEK inhibitors, we were then interested in studying how expression levels of the pathway (phospho) proteins correlated with their drug sensitivity.

In order to identify relevant signaling proteins whose observed basal expression or phosphorylation levels across different cell lines can explain corresponding DSS values for different drugs, our first approach was a correlation analysis using Pearson correlation coefficients. Correlation coefficients for all possible protein-drug pairs, including the 33 drugs in our library, and the 31 relevant signaling proteins (16) selected for our study, were computed using the phospho flow readouts and DSS values from all the MM cell lines ( Supplementary Figure 1B ). Multiple significance tests of the Pearson correlation coefficients followed by a false discovery rate controlling procedure revealed that this method is not sufficiently robust to capture significant correlations. Our solution was to model all the drugs jointly using the tree-lasso regression approach described in Materials and Methods. By penalizing regression coefficients jointly according to the hierarchical clustering tree of the correlations among DSS values of different drugs, we were able to achieve enough sensitivity to capture protein-drug correlations that had been lost in our previous, naïve approach ( Figure 4A ). We found that PI3K and MEK inhibitors emerged as separate clusters, each with a common set of predictive variables, as indicated by a few non-zero tree-lasso regression coefficients ( Figure 4A ). We note that since the covariates in our model were defined as the standardized phospho flow readouts, the tree-lasso regression coefficients have a similar interpretation to correlation coefficients (i.e., positive or negative values for the coefficients have similar meanings in both cases), and therefore can be regarded as a special type of correlation.

Notably, using the tree-lasso regression model we found that high levels of phospho-TBK1 inversely correlated with low drug sensitivity to MEK inhibitors, including trametinib, binimetinib, and cobimetinib ( Figure 4A ). Although these correlations seem to be in agreement with individual Pearson correlation plots of the TBK1 levels against DSS for trametinib and binimetinib (upper plots in Figures 4B, C ), we emphasize that their true significance is only revealed by an integrated model that explores the data structure of all the drugs to define a new, more robust, type of correlation.

As also observed in our previous study (16), we found that high levels of Bcl-2 and phospho-Bcl-2 predicted low sensitivity to trametinib, binimetinib, and cobimetinib ( Figures 4A-C middle and lower plots). Once again, individual Pearson correlation plots ( Figures 4B, C middle and lower plots) are not enough to identify these relationships, since their statistical significance were only seen after using the more robust tree-lasso regression approach for variable selection. Also, worth mentioning is our finding that phospho-MEK1 (pS298) is positively correlated (in the sense of positive tree-lasso regression coefficients) with sensitivity to the same MEK inhibitors (trametinib, binimetinib, and cobimetinib). Moreover, we observed that high phosphorylation levels of AKT (pS473) correlated with high responses to PI3K inhibitors taselisib, pictilisib and copanlisib ( Figure 4A ). In this case, Pearson correlation p-values are too high even before multiple-testing correction ( Supplementary Figures 2A-D ), however, tree-lasso regression was still sensitive enough to capture a relatively weak correlation.

Taken together, our results indicate that the expression and phosphorylation levels of signaling proteins can inform on drug sensitivity, and therefore, can provide relevant information as part of functional drug testing studies.