Predicting cooperative drug effects through the quantitative cellular profiling of response to individual drugs.

Quantitative prediction of cellular responses to drugs and drug combinations is a challenging and valuable topic in pharmaceutical research. In the past decade, microarray technology has become a routine tool for monitoring genome-wide expression changes and has been widely adopted for exploring drug response in the pharmaceutical field. However, how to predict the synergistic effect of drug combinations using microarray data is a challenging task. In this article, we report a simple prediction framework based on the genome-wide and quantitative profiling of cellular responses to individual drugs. By exploring the differential expression profiles, our correlation-based strategy can reveal the synergistic effects of drug combinations. The comparison with gold-standard experimental results demonstrates the strengths and weaknesses in relation to prediction based only on cellular response to individual drugs. Specifically, the prediction strategy may work for a drug combination whose individual drugs show related transcriptomic mechanisms but not for others.

An understanding of the complex biological responses of the human organism to drugs is key to investigating the efficacy and safety of compounds in drug development. Many genomic features, including DNA methylation patterns, messenger RNA levels, and protein expression or metabolite profiles, may be used for monitoring biological responses. Microarray is currently the least expensive high-throughput technology for simultaneously monitoring genome-wide expression profiling of transcriptional response to drugs.

Microarray data have been systematically explored in model organisms to elucidate the drug mechanism of action, and coexpression analysis enables the inference of functional roles for genes that respond coherently to drug perturbations. The National Cancer Institute's NCI-60 project and the Connectivity Map have extended the concept of genome-wide gene expression profiles of drug response to human cell lines.[1,2] The NCI-60 project screened 60 human tumor cell lines against more than 100,000 compounds and constructed a public repository for basal gene expression and drug sensitivity information. The Connectivity Map project generated genome-wide expression profiles both before and after drug treatment for 1,309 compounds and constructed a drug network by comparing ranked lists of up- and downregulated genes.[3] More recently, Iskar et al.[4] identified a large set of drug-induced transcriptional modules from genome-wide microarray data of drug-treated human cell lines and rat liver and examined their conservation across tissue types and organisms. The large-scale Cancer Cell Line Encyclopedia project extended such analysis and generated gene expression, chromosomal copy number, and massively parallel sequencing data from 947 human cancer cell lines. They analyzed the data using pharmacological profiles for 24 anticancer drugs across 479 of the cell lines and identified genetic-, lineage-, and gene expression–based predictors of drug sensitivity.[5] All these efforts have enabled the systematic translation of molecular data into biological knowledge and therapeutic possibilities.

However, it is well known that redundancy and multifunctionality are inherent characteristics of biological systems. Most complex diseases are usually caused by the interweaving of multiple biological processes that are redundant and robust to the perturbation of one single drug.[6,7] Thus, selecting proper drug combinations that can yield a synergistic effect is thought to be an effective way of countering biological buffering and increasing therapeutically relevant selectivity.[8,9,10,11,12] By modulating the activity of distinct proteins through drug combinations, we can gain not only the potential of overcoming the redundancy underlying pathogenic processes but also fewer side effects.[13,14] For example, estrogen, which modulates feeding behavior and energy expenditure, has been suggested as a therapy for obesity and type 2 diabetes; however, it can also cause unwanted effects, such as tumor growth promotion. By chemically linking estrogen to the hormone GLP-1 to selectively target metabolically relevant tissue, Finan et al.[15] obtained a more selective and potent metabolic drug with fewer side effects.

Despite these recent advances, understanding the molecular mechanisms underlying drug combinations is challenging. High-throughput screening is useful to identify possible drug combinations; however, it is impractical to screen all possible drug combinations for all possible indications. We still have no practical experiment strategy to identify possible effective drug combinations. Several computational methods have recently been proposed.[16,17,18,19,20,21,22] Molecular and pharmacological features of drugs are the common basis of these computational methods. By integrating several of these features, such as their targets, indications, or chemical structures, researchers can identify the similarity between each drug pair. The drug pairs with significant high similarity can be considered candidates for possible effective drug combinations. Signaling pathways and human protein–protein interaction networks are two other resources on which many models have been based. Jin et al.[18] proposed an enhanced Petri-Net model to identify the targeted signaling networks of two drugs separately and their combination. By comparing identified drug effects, they could verify whether a drug combination shows synergistic effect.

Only a few studies have been devoted to predicting synergistic effects from gene microarray data.[18,19] Their common strategy is comparing the drug effects derived from the microarray treated by drug combination with those from the corresponding two microarrays treated by two drugs separately. In our study, we investigated whether cellular responses of drug combinations can be modeled quantitatively by the cellular profiling of response to individual drugs. We demonstrated quantitative modeling based on the expression profiles acquired after treatment with 14 individual compounds from the NCI-DREAM Drug Sensitivity Prediction Challenge. Because the challenge is greater in the absence of cellular response profiling of drug combinations, we needed to find common features from their respective response data. Our hypothesis underlying the current method is that the differential expressed genes after administration of one drug can express the action mechanism of this compound. The combination of two drugs with correlated or anticorrelated differential expression profiles may have, respectively, a synergistic or an antagonistic effect.

Result

Overview of the analysis

Recently, the Dialogue on Reverse Engineering Assessment and Methods (DREAM) (http://www.the-dream-project.org/challenges/nci-dream-drug-sensitivity-prediction-challenge) consortium designed a community-based, collaborative competition for the systematic and objective validation of computational methods to predict the cooperative effects of 14 distinct drugs/compounds (M. Bansal et al., personal communication). This project generated gene expression data only for samples treated by individual drugs, not pairwise combination in a human β-cell lymphoma cell line lymphoma cell line (), which is a key difficulty. Predictions from 31 different methods have been experimentally assessed, and the present method (with slight modifications here) has been ranked fourth. Such results confirm that our method is a simple but promising one to predict cooperative effects of drug pairs.

The goal of this study was to rank the activity of all 91 pairwise combinations of 14 distinct drugs/compounds in a human diffuse large B-cell lymphoma (DLBCL) cell line (LY3) from the most synergistic to the most antagonistic (). We used only the baseline expression profile and the expression profiles treated with each individual drug/compound, with a limited number of conditions and replicates to predict synergistic effects of drug/compound combinations. The basic assumption of our method is that the expression properties of differential expressed genes obtained by comparing the expression profile of a cell line before and after treatment with a drug reflect the effect of this drug to that cell line. We further assume that the correlation of expression profiles of differential expressed gene sets for two drugs will reflect the synergistic and antagonistic effects of the two drugs (see and Methods). All results and codes are provided in the Supplementary Information.

The predictions were then evaluated against an experimentally assessed “gold standard” generated by evaluating compound combination activity in vitro (M. Bansal et al., personal communication). A measure of the concordance between the predicted rank list and the gold standard known as the concordance index (c-index) was computed to score the prediction. In essence, the c-index computes the fraction of item pairs whose relative order in the two rank lists is coincident.

Differential expressed genes indeed provide valuable information for predicting cooperative drug effects

In our method, we set a P value threshold of 0.05 to choose the significantly differentially expressed genes as signature genes used for the following correlation analysis. Using this signature gene set, we can get a c-index of 0.599, which is very competitive compared to the top two methods, with c-indexes of 0.613 and 0.605, respectively (M. Bansal et al., personal communication). To show its effectiveness, we also adopted a pathway-based approach for the current task. Specifically, we mapped the differential score of expressed genes using the Human Protein Reference Database (HPRD) and then detected the significant “active” pathway which corresponds to the maximal weighted connected subnetwork by the BioNet method. However, it only shows equivalent performance compared to the correlation-based strategy for this task the in the current setting with weighted c-indexes of 0.608 and 0.599, respectively.

We conducted a random perturbation experiment to test the significance of the selected differentially expressed genes. For each drug combination, we randomly extracted the same number of signature genes as the overlap we obtained using a P value of 0.05. Then we calculated the correlation of expression profiles of these two drugs based on the randomly selected signature genes. The selection and computation were repeated three times, and the one with the maximal absolute value was selected as the synergistic score. Then we used the scoring scripts to calculate one c-index. The entire procedure was repeated 1,000 times (the distribution of these c-index values is plotted in ). We can see that the differential analysis indeed has the power to identify the most relevant signature genes for the predicting task (P value = 3.39 × 10−5; t-test).

We further assessed the robustness of the P-value threshold to see how it affected the results. shows the values of the c-index for each threshold. We obtained the more informative predictions with a P value = 0.05 than with others, especially those with a larger value. We should also note that with a P value <0.001, we could get only a limited number of differential genes that could not be used for the correlation analysis.

The drug combinations with top consistent predictions show similar functional mechanisms

We further inspected our predictions compared with the experimental “gold-standard” rank (). We found that some top predictions that also have relative top ranking by the gold standard show highly similar functional mechanisms. For example, the predicted top 1 drug combination (mitomycin and etoposide) was ranked 3 in the gold-standard list. We analyzed their overlapping signature genes by functional enrichment analysis and extracted connected subnetworks by mapping them onto a protein interaction network. We found that the common signature genes are enriched in DNA repair, DNA metabolic process, DNA replication, cell cycle biological processes, and the p53 signaling pathway (). More specifically, we observed that treatment with these two drugs can cause upregulation of p53 signaling genes (CDKN1A, PPM1D, TP53I3, SERPINE1, DDB2, MDM2, RRM2B, FAS, GADD45B, SESN2, SESN1, and GADD45A) and downregulation of cell cycle genes (E2F3, E2F5, DBF4, SKP2, MCM3, CDC27, MCM4, WEE1, MCM6, CDKN1A, PLK1, BUB1, BUB1B, MDM2, ANAPC7, GADD45B, GADD45A, CCNA2, and TFDP1) in the diffuse large B-cell lymphoma LY3 cell lines by checking their overlapping signature genes.

We further found that existing studies have shown that these two drugs have similar action mechanisms. It has been reported that mitomycin is activated in vivo to a bifunctional and trifunctional alkylating agent which binds to DNA and leads to cross-linking and inhibition of DNA synthesis and function.[23,24,25] Moreover, previous studies have shown that mitomycin is a cell cycle phase-nonspecific agent. Etoposide inhibits DNA synthesis by forming a complex with topoisomerase II and DNA. This complex induces breaks in double-stranded DNA and prevents repair by topoisomerase II binding. Accumulated breaks in DNA prevent entry into the mitotic phase of cell division and lead to cell death. The key point is that etoposide acts primarily in the G2 and S phases of the cell cycle.[26,27] All the analysis demonstrates that the common action mechanism of these two drugs at the transcriptomic level explains why the current strategy can predict it well.

The drug combination monastrol and trichostatin A, ranked 18 in the gold-standard list, has been predicted to be rank 19 among all (). Through functional enrichment analysis, we found that the overlapping signature genes of the two drugs are involved in DNA replication and cell cycle pathway. Previous studies have demonstrated that these two drugs are the protein kinase inhibitors and can induce cell cycle arrest. Monastrol targets kinesin-related motor protein KIF11, which is required for establishing a bipolar spindle and for mitosis. Blocking of KIF11 can prevent centrosome migration and causes cell cycle arrest in mitosis.[28,29,30] Trichostatin A is a classical histone deacetylase inhibitor. It has been demonstrated to induce cell cycle arrest, promote cell apoptosis, and inhibit metastasis.[31] We believe that the similar functional roles of these two drugs, reflected at the transcriptomic level, enable our method to reliably predict their synergistic effect.

However, we found that some drug combinations cannot be predicted well. For example, the top drug combination on the gold-standard list—doxorubicin and H-7—is predicted to be ranked 83 by our method. In enrichment analysis, we found that doxorubicin mainly alters the expression of genes in the p53 signaling pathway, cell cycle, and DNA replication, while H-7 mainly affects steroid biosynthesis, the spliceosome, and the ribosome. Previous studies have also shown that doxorubicin interacts with DNA by intercalating and inhibiting macromolecular biosynthesis, and induces a break in the DNA chain, preventing the DNA double helix from being resealed, thereby stopping the process of replication,[32][33][34] but H-7 (protein kinase C-α antibody) is a protein kinase inhibitor.[35] Protein kinase C has been shown to act as a tonic negative regulator of basal steroidogenesis in Y1 cells by suppressing the expression of mRNA encoding the steroid synthetic enzymes.[36] This analysis may show that the predicting task based only on gene expression in treatment with individual drugs would fail when the effects of individual drugs are different. It will be necessary to incorporate more types of cell-response information to improve the accuracy of prediction.

The drug's different time spans to generate their maximal effects can affect prediction

Drugs work on cells differently, and their maximal effects can occur at different time points. We used a pairwise corresponding comparison between time points, which may limit the predicting ability. To test if the different time spans to generate maximal drug effects affect the prediction, we added six cross-time point comparisons and then chose the best one for the ranking. Surprisingly, based on this modified strategy, we got a better prediction with c-index = 0.628 (), which is even higher than the top two methods with c-indexes of 0.613 and 0.605 in the challenge. We found that this improvement is mainly due to the prediction of the drug combination blebbistatin and H-7. This combination is ranked 7 in the gold-standard list. In our original strategy, this combination was ranked 79, whereas in the modified version it was ranked 3. This analysis shows that considering the time point of maximal drug effects is helpful to improve prediction accuracy.

Discussion and Conclusion

Predicting the synergistic effects of drug combinations using the quantitative gene expression profiles of cellular response to drugs is challenging. The current study shows that prediction based only on the gene expression profiles in treatment with individual drugs, not in combination, is even more challenging. We attempted to define the gene signatures that reflect the effect of all individual drugs using differential expression analysis and calculated the correlations of those gene signatures to predict the synergistic effect of drug/compound combinations. In other words, we assumed that the effects of individual drugs were measured in the expression levels of their regulated genes. The expression correlation of common regulated genes of two individual drugs was used to approximate the synergistic effect.

The method presented here is very simple but shows very competitive performance compared with other methods in the challenge. We have shown that the differential gene signatures are more informative than randomly selected ones, and the method can capture common functional genes for some top-ranked drug pairs. The pathway analysis shows that the cellular response to individual drugs can partially explain how it works. The significant consistency between our predictions and the experimental gold standard demonstrates that this approach can effectively identify synergistic drug combinations via the cellular response to individual drugs. A potential limitation of this method is that it relies on the transcriptional signature of individual drugs and omits other complicated combinatorial transcriptional features of drug combinations in real application. The key hypothesis here is based on the assumption that similar effects should have similar gene signatures no matter what effects the two drugs have. However, this may be true for chemotherapy drugs but not true for other targeted-therapy drugs. Thus, the prediction strategy here may work for a drug combination whose individual drugs show related transcriptomic mechanisms but not for others. Another potential limitation is that the correlation-based method is not totally in sync with the prediction of the synergy. However, the performance compared to other methods demonstrates that the correlation strategy can still capture some underlying connections between drugs.

The simple method for predicting cooperative drug effects confirms the potential predictive ability based only on the quantitative cellular profiling of response to individual drugs. With the development of systems biology, we believe that the current method can be potentially improved in several aspects: (i) generating treated expression profiles of cell lines with drug combinations and building a model to compare them with those of each drug and (ii) increasing the number of samples and time points. Moreover, integrating other types of response data may be helpful to improve the performance. For example, Winter et al.[37] integrated transcriptional data with proteomic data to dissect the synergy between two multikinase inhibitors in chronic myelogenous leukemia cell lines. Finally, we note that integrating pathway and network data to improve the interpretation of synergy prediction will be necessary in future studies.

Methods

Data

We were given the genomic data on the LY3 diffuse large B-cell lymphoma cell line for this study. These data were collected both before and after treatment with isolated compounds. More detail is provided at http://www.the-dream-project.org/challenges/nci-dream-drug-sensitivity-prediction-challenge.

Treated and untreated gene expression profiles were used in this study. The treated data included a set of 150 gene expression profiles corresponding to a set of 14 perturbations (14 compounds in DMSO), at three time points (6, 12, and 24 hours), at one concentration corresponding to the IC20 of the compound at 24 hours, in triplicate. The untreated data included the gene expression profiles of cells exposed to DMSO in eight replicates at all three time points. All gene expression profiles were profiled using an Affymetrix U219 96 array plate (Santa Clara, CA).

We also obtained the IC20 of the compounds and the synergistic activity of drug combinations assessed within the same experiment in triplicate to avoid differences in compound potency or cell line physiology in different experiments. The experimental synergistic activity was used to rank the computational predictions as “gold standard.”

We downloaded a human binary protein interaction network from the HPRD database (http://www.hprd.org/), which consists of 30,047 proteins and 41,327 interactions.

Preprocessing

We found that the gene expression profiles of the three replicates were highly correlated with each other (Pearson correlation coefficient, r > 0.99). We used the average of the three replicates (after normalization) as the unified profile of each drug.

Method

In the following, we describe the detailed procedure of our method ():

The gene expression profile after administration of DMSO vehicle was regarded as the baseline profile for computing the differential expression after treatment with each drug at three different time points.

For each drug, we extracted the genes with an absolute differential expressed z-score > 1.96 (P < 0.05) across all genes for the same time point as signature genes.

For each pair of drugs, we calculated their correlation of expression profiles based on the overlap of their signature genes.

We selected the one with maximal absolute value among three correlation values corresponding to three time points as the synergistic score for each pair of drug and ranked the drug pairs according to this score.

Network-based method

We mapped the differential expression genes onto the HPRD network and then extracted the maximal weighted connected subnetwork using the BioNet package[38] to do the same analysis as in our correlation-based method. The detailed procedure is as follows:

Extract the significantly differentially expressed genes if their corresponding |z-score| > 1.96.

Download a human binary protein interaction network from the HPRD database (http://www.hprd.org/).

Assign the nodes of the significantly differentially expressed genes weights as their absolute z-score. The other nodes were assigned weights based on their corresponding |z-score| – 1.96.

Maximal score connected subnetwork is extracted using the R package BioNet.

Do the same downstream analysis as in our present method and calculate the weighted c-index.

The c-index

Here, a modified version of the c-index called the probabilistic c-index, to account for the probabilistic nature of the gold standard, was used in authoritative scoring scripts. The probabilistic c-index was calculated by comparing the predicted rank list of drug pairs to the experimentally determined ranked list of drug combinations. Given a pair of drug combinations j and k, a prediction rank was assigned a score given by the probability that the noisy gold standard supports this prediction. The probabilistic c-index is the average, for all drug pairs, of the probabilities that the gold standard supports the prediction.

Functional enrichment analysis

We adopted the DAVID online tool (http://david.abcc.ncifcrf.gov/) to do functional enrichment analysis. The pathways and Gene Ontology terms with adjusted P values (Bonferroni correction) <0.01 were selected as significant ones.

Author contributions

S.Z. and J.Z. designed the research. J.Z. and S.Z. performed the research, J.Z., X.-S.Z., and S.Z. analyzed the data. S.Z. and J.Z. wrote the manuscript.

Conflict of Interest

The authors declared no conflict of interest.

Study Highlights

Table 1

The top-ranked drug combinations per the gold-standard method and their predicted ranking by our method