Prediction of drug response based on patients’ clinical and molecular features is a major challenge in personalized medicine. Great effort has been devoted to identify molecular biomarkers of drug sensitivity and to develop computational models to predict drug response based on these biomarkers. Gene expression data is one of the most commonly used data type in these studies, due to their high predictive ability, and numerous methods have been proposed for drug response prediction based on gene expression data [1, 2, 3, 4, 5, 6, 7]. However, many existing methods only use basal gene expression data (i.e., gene expression values before administration of the drug) and hence can only capture the influence of the steady state of the cells on their response to a drug. For example, the authors of  analyzed 44 drug response prediction methods that employed gene expression profiles of breast cancer cell lines taken before treatment to predict dose-response values, e.g., –the concentration that inhibited cell growth by 50% after 72 hours of treatment. Similar work can also be found in [2, 3, 4], which only incorporate gene expression data from a single time-point. A collection of temporal gene expression profiles of samples over a series of time-points during the course of a biological process can provide more insights than a single (or two) time-point(s) . Therefore, developing algorithms that can predict the drug response over time using time-course gene data is of great interest.
With the advancement of gene sequencing technologies, collecting gene expression data over multiple time-points and their matched drug response values is now feasible. In parallel with these technological developments, there has been growing interest in the application of machine learning methods to analyze the time-course gene expression data. In
, the authors proposed an integrated Bayesian inference system to select genes for drug response classification from time-course gene expression data. However, the method only uses the data from the first time-point, and hence does not benefit from the additional temporal information. Lin et al.,
presented a Hidden Markov model (HMM)- based classifier, in which the HMM had fewer states than time points to align different patient response rates. This discriminative HMM classifier enabled distinguishing between good/bad responders. Nevertheless, choosing the number of states for this HMM is a major practical issue. In addition, this method cannot handle missing data and it requires the full knowledge of GELs in all time-pointsa priori. This implies that the HMM may not be able to predict drug response at multiple stages in future time points, since the corresponding GELs are not measurable.
The time-course gene expression data contains the GELs of different patients over a series of time points, which can be indexed as patient-gene-time and represented as a three-dimensional tensor. Motivated by this, several tensor decomposition based algorithms have been proposed. For example, Taguchi  employed tensor decomposition to identify drug target genes using time-course gene expression profiles of human cell lines. Li and Ngom  proposed a higher-order nonnegative matrix factorization (HONMF) tool for classifying good or bad responders from a latent subspace corresponding to patients learned from HONMF. One limitation of this work is that the latent subspace may not have discriminative ability in classifying patients, since it is learned without accounting for the class-label information. Moreover, this method simply discards samples with missing values, causing unnecessary information loss.
Recently, Fukushima et al., 
developed an algorithm for joint gene selection and drug response prediction for time-course data. The method uses Elastic-Net (EN) to select a set of genes that show discrimination of patients’ drug responses throughout the treatment. The selected genes are then passed to a logistic regression (LR) classifier for drug response prediction. But in real applications, due to the existence of noise and missing values in the data, finding discriminative genes for all patients may be difficult. In fact, several studies have shown that it is more viable to find genes that have consistent discrimination in a subset of samples along the time series[13, 14, 15]. Therefore, relying only on discriminative gene selection but without modifying classification algorithms may not achieve satisfactory performance.
In this paper, we take a different approach for time-course drug response prediction. We hypothesize that a patient’s drug response at a given time-point can be inferred from the response at a previous time point. This means that not only the GELs but also the past response results can be integrated to identify the drug response for a subsequent time point. We develop a REcursive Prediction (REP) algorithm to predict the drug response of samples using their time-course gene expression data and their drug response at previous time-points. REP has a built-in recursive structure that exploits the intrinsic time-course nature of the data through integrating past drug responses for subsequent prediction. In other words, in REP, not only the GELs but also the past drug responses are treated as features for drug response prediction. Furthermore, by taking into consideration the intrinsic tensor structure of the time-course gene expression data and leveraging identifiability of low-rank tensors, REP can alleviate the noise corruption in GEL measurements, complete missing GELs and even predict GELs for subsequent time points. These features enable REP to evaluate drug response at any stage of a given treatment from some GELs measured in the beginning of a treatment. Experiments on real data are included to demonstrate the effectiveness of the REP algorithm.
Fig. 1 sketches the idea behind the proposed REP algorithm, where the subfigures 1(a)–1(c) show the pre-processing, model training and prediction of our method, respectively. The tensor structure of time-course gene expression data is shown in Fig. 1(a). In the following, we explain them in more detail.
One major issue in using gene expression data for drug response prediction is the existence of missing values. To overcome this problem, we first impute the missing values during pre-processing. Various methods have bene previously suggested for handling missing values, such as median-imputation and nearest neighbor imputation [17, 18]. We employ a low-rank tensor model to fit the time-course gene expression dataset such that the missing values can be completed. Our supporting hypothesis is that genes never function in an isolated way, but oftentimes groups of genes interact together to maintain the complex biological process, which results in correlation in the GEL data . We note that our low-rank tensor model suggests three factors that uniquely determine the values of GELs, i.e., the factors corresponding to patient, gene and time, respectively (see Fig. 1). As we will see later, our model allows us to estimate the variation of GEL over time from a set of initial GEL measurements; these estimated values are then used to predict the time-course of drug response.
Towards this goal, we first assume111For high-enough but finite , any patient-gene-time dataset can be expressed this way. See  for a tutorial overview of tensor rank decomposition. that each GEL is represented as a summation of triple products from the latent factors of patient, gene and time, respectively. Let us denote as the th GEL of patient recorded at time . Based on our assumption, we have
where , and are the latent factors of patient, gene and time, respectively. Suppose that there are genes measured over time points. By varying the indices and in (1), the expression of the genes in all the time-points in patient can be represented as
where , , . In this equation, represents a diagonal matrix holding the th row of as the main diagonal, which is a latent representation of the th patient. We use to represent the -entry of , to represent the -entry of and to represent the -entry of .
Assume that there are patients in the training set. After collecting , we stack them in parallel along the patient-axis, which results in a GEL tensor that takes the form of
where is the outer product and is the th column of , and likewise for and . Here, we assume that is the noiseless GEL data and is the corresponding noisy data with missing values. The relationship between and is described as
where is the noise in the data, is the index set of the observed GELs in , and is the operator that keeps the entries in and zeros out the others.
The model in (2) indicates that the gene and time factors (i.e., and ) are identical for different patients, and the variability among patients is captured by . In other words, given and , each row of the patient factor matrix uniquely determines the GELs of the corresponding patient. As we will see later, our model is able to predict unseen GELs, which also enables to prescreen the drug response for different stages of a treatment.
Assuming nonnegative GELs222 Due to some preprocessing steps such as z-score normalization, the GEL values can be negative. To facilitate our method, we undo these preprocessing steps or use the raw dataset.
Due to some preprocessing steps such as z-score normalization, the GEL values can be negative. To facilitate our method, we undo these preprocessing steps or use the raw dataset., we can use nonnegative tensor factorization to compute missing GEL values:
where many sophisticated algorithms are applicable to optimize (II-A), e.g., block coordinate descent [21, 20]. Intuitively, (II-A) seeks to identify the lowest rank solution that best matches the observations . The regularization terms are added to further encourage low rank and prevent over-fitting. When (II-A) is solved, we complete the GEL data through
where contains the indices of missing values in .
The effects of drugs are usually cumulative over time , i.e., drug doses taken in the past will affect the current response. This implies that the drug response in the past time-points may help predict the current response. Based on this hypothesis, we propose a recursive prediction algorithm, henceforth referred to as REP for simplicity, which enables to integrate past drug response records with gene expression values for subsequent drug response predictions. Fig. 1(c) shows an overview of REP’s pipeline, where drug responses in the previous time stages are integrated with the gene expression information for predicting the current response . Here, we accumulate the historical response as
which is then fed back as an input feature for subsequent drug response prediction. The intuition behind (7) is that if the responses at adjacent time stages are consistent, then the response at the next stage is more likely to continue the same pattern instead of going in the opposite direction. If the responses at two adjacent time stages are opposite, i.e., one is good and the other one is poor, their effect will be canceled out in (7). In this case, the predictor will focus more on the GELs.
Mathematically, we have
where denotes the accumulated historical responses of patient at time , is the predictor which can be trained by minimizing the following cost function
where contains the parameters of the predictor,
is the loss function of a classifier such as hinge loss and cross-entropy loss,is a regularizer that imposes a certain structure on , and is a regularization parameter. In the literature, popular regularizers include , , and , i.e., the indicator function of the nonnegative orthant.
For the purpose of illustration, we are particularly interested in but not limited to the SVM classifier. We set and as hinge loss, resulting in
where is the intercept, denotes a convex set such as
-ball for feature selection andrepresents the importance of the response at a previous time point on the subsequent one. In our formulation, the drug response feedback plays an important role and it can be viewed as a “must-have” feature. Therefore, the hyper-parameter is critical for the performance of REP. In practice, we recommend to choose a relatively large in an incremental fashion.
The major difference between (II-B) and the standard SVM-based drug response predictors lies in the feedback : prior art ignored the previous drug response outputs. For the initial time point, there is no response feedback, i.e., does not exist. So how do we choose ? One way could be which can be interpreted as uncertainty, i.e., we do not know how the patient responds to a treatment. Then we construct the following feedback matrix:
where . Finally, the features in the training set is formed by concatenating in along the gene-axis as shown in the left-bottom corner of Fig. 1(b), and the training labels are
Ii-C Drug response prediction
Our method can predict the drug response values for a new patient at any time point. Specifically, given the GELs of a new patient at time , i.e., , we first check if there are missing values. If so, we employ the factors and to complete . Let us denote and as the sets of indices of the observed and missing elements in . According to our model in (1), can be uniquely determined by ,
and an unknown vector–a latent representation of this new patient. Thus, for the expression level of the th gene at time , we have
is the additive noise which is assumed as Gaussian distributed,is the Khatri-Rao (column-wise Kronecker) product, and and denote the th row of and , respectively.
Since and are known, the problem of estimating can be formulated as
which is a nonnegative least squares (NLS) problem and can be optimally solved. We note that to obtain a unique estimate , the number of available gene expression entries in should be . The GEL vector of the patient is then estimated as
which leads to a completed GEL vector
The vector together with the cumulated historical drug response , are the input data for our predictor . We estimate the drug response of this patient at time via
Ii-C1 Predicting Unseen GELs
Previously, we have explained how to predict drug response for patients at a certain time point. However, in practice, we are more interested in knowing the drug response of a few time-points in the future from the beginning of a treatment. This requires to know the GELs of all time points up to the time-point of interest a priori, which is impossible in practice. In this subsection, we provide an efficient solution that allows to predict the unseen GELs.
Recall that in our model, the GEL of a patient is determined by three factors, i.e., the latent representation of patient–, the time evolution factor– and the gene factor matrix–, where is different for patients, and needs to be estimated for the new patient. On the other hand, and are common gene and time evolution bases that reflect different types of patients, as determined from historical patient data – the training data. Therefore, the problem boils down to the estimation of from the initial GELs of the new patient. We can simply substitute in (14) to find . Finally, the GELs for the remaining time points are estimated as
Now we have estimated the unseen GELs for , which allows us to predict drug response values for the whole duration of the treatment. We start from and estimate the drug response for as
where . When is available, we set . With the GEL estimate from (18), we can predict , and so forth for the other time points.
Note that here we substitute predicted drug responses for the unseen drug responses. Clearly, when actual drug responses for past time ticks are available, they should be used. We only do the substitution here for a preliminary assessment of how well a patient is likely to respond over time, before the beginning of treatment – which is naturally a more ambitious goal.
In this section, we provide some numerical experiments to showcase the effectiveness of REP for drug response prediction from time-course gene expression data. We apply a number of drug response prediction methods including two linear models (EN-LR  and SVM), one nonlinear model (
-nearest neighbor (KNN)), and a probabilistic graphical model (discriminative loop hidden Markov model (dl-HMM) ) to real-world time-course data. We did not include SVM with nonlinear kernels (e.g. Gaussian), since its performance was inferior compared to the linear kernel. Note that EN-LR and dl-HMM were specifically designed for prediction of drug response values based on time-course gene expression data, while SVM and KNN are widely used classification methods.
The dataset used is the interferon (IFN)- time-course dataset which is available in the supplementary of . The dataset is collected from 53 Multiple Sclerosis (MS) patients who received IFN- treatment for two years. The gene expression data (microarray) was obtained from peripheral blood mononuclear cells of the patients, which contained the expression levels of 76 pre-selected genes over 7 stages (i.e., time-points) of the treatment. However, there are missing values in this dataset, where most missing values were caused by the absence of patients at some stages. Among the patients, only 27 patients have records for all stages, while the other 26 patients missed at least one stage, which resulted in the entire GELs as well as the drug response at that stage being missed. For ease of comparison, we only employ the 27 full records in our experiments, where the final GEL data is with dimension and the response data is with dimension .
Iii-B Evaluation metric
We use prediction accuracy (ACC) and area under receiver operating characteristic (ROC) curve (AUC) to evaluate the performance of REP, where ACC is defined as:
and ROC curve is created by plotting true positive rate (TPR) versus false positive rate (FPR), defined as
In the equations above, TP, FP, FN and TN stand for the number of true positives, false positives, false negatives and true negatives, respectively. The ACC is calculated using 27-fold cross validation, where in each fold, we select one patient’s record as a testing set that contains a GEL matrix and a response vector with length 7, while the remaining 26 records are assigned to the training set.
The hyper-parameters of EN-LR, SVM and KNN are tuned through 5-fold cross validation, while those of REP are tuned slightly differently: we first split the training data into a training set and a validation set consisting of 25 and 1 records, respectively, and then determine the final parameters as those that achieve the highest prediction accuracy in the validation set. It should be noted that the parameters of REP are tuned once but not in a CV way, meaning that REP is tuned not as accurately as the others. For REP, its hyper-parameters include the rank , and in (II-B), which are selected from and . The SVM method solves the following problem:
where its parameter is tuned from ; for EN-LR, we set which is a hyper-parameter balancing the ridge and LASSO regularizations; for KNN, the number of neighbors is selected from . After that, we apply the trained classifier to the testing data to calculate ACC. We implemented REP, EN-LR, SVM and KNN in Python 3.7. Since the authors of dl-HMM have published their MATLAB codes (http://www.cs.cmu.edu/~thlin/tram/), we implemented these two methods in MATLAB. For both HMM methods, we choose the number of hidden states as 4 which is as suggested in .
Fig. 2 shows the prediction accuracy of the five algorithms. In this example, the missing values is only 0.23%. For all the methods, the missing GELs were completed using nonnegative tensor completion in Section II-A. It can be seen that REP achieves a higher prediction accuracy compared to other methods and its performance is followed by the EN-LR and SVM algorithms. The KNN and dl-HMM algorithms do not perform well. We observed that dl-HMM was very sensitive to the data, and it frequently produced not-a-number (NaN) values. Particularly for patients with ID 1117161, 1215133 and 995640, the original codes of dl-HMM always produced NaN of the three patients’ responses. This is the main reason why dl-HMM had low prediction accuracy. We note that REP takes a similar formulation as the SVM. However, REP is 6.2% more accurate than SVM, which implies that the recursive structure in REP is helpful in improving the prediction accuracy. In Fig. 3, the ROC curves are plotted, which shows that the proposed method has the highest AUC value.
Next, we sought to determine the effect of missing values on the performance of these methods. For this purpose, we randomly sampled the GEL data and hid the selected entries. As the percentage of missing values increases, we can see in Table I that all methods suffer performance loss, but REP’s ACC and AUC remain the highest in all cases. We highlight that when the percentage of missing values is smaller than 15%, REP has ACC close to 90% and AUC greater than 95%. EN-LR outperforms the classical SVM method in many cases. When the percentage of missing values increases, the performance of EN-LR and SVM drops significantly, while that of REP still remains at a high level. For example, when percentage of missing increases to 15%, the AUC of EN-LR drops to 87.4%, but ours is 95.4%, which indicates that REP is more robust against missing values.
The above examples were conducted under the full knowledge of GELs for each time point. However, we are more interested in predicting the drug response in the beginning of a treatment, meaning that we only have a few GELs at the first time point while all the subsequent GELs are unavailable. Specifically, we can first learn factors using nonnegative tensor completion, then follow the procedures in Section II-C1 to predict the GELs, and finally, use the estimated GELs for drug response prediction. To show the effectiveness of REP, we only used the GELs of new patients at the first time point, i.e., , while hiding the GELs for other time points through to . We studied the performance of REP where only GELs for the first time-points were available. Fig. 4 shows the prediction accuracy of REP by using estimated GELs, where the percentage of missing stands for the missing values in the training set. We see that REP with estimated GELs achieves about 85% accuracy. Although this performance is not as good as REP with the ground truth GELs (see Fig. 2), we must emphasize that such performance is obtained by only knowing the GELs from the first stage of the treatment, which is more meaningful in practice.
We studied the problem of drug response prediction for time-course gene expression data and presented a computational algorithm (REP) that: (i) has a recursive structure that integrates past drug response records for subsequent predictions, (ii) offers higher prediction accuracy than several classical algorithms such as SVM and LR, (iii) exploits the tensor structure of the data for missing GEL completion and unseen GEL prediction, (iv) can predict drug response of different stages of a treatment from some initial GEL measurements. The achieved performance improvement for real data application suggests that our method serves as a better predictor of drug response using the time-course data.
Conflict of Interest: none declared
-  J. C. Costello, L. M. Heiser, E. Georgii, M. Gönen, M. P. Menden, N. J. Wang, M. Bansal, P. Hintsanen, S. A. Khan, J.-P. Mpindi et al., “A community effort to assess and improve drug sensitivity prediction algorithms,” Nature biotechnology, vol. 32, no. 12, p. 1202, 2014.
-  C. Suphavilai, D. Bertrand, and N. Nagarajan, “Predicting cancer drug response using a recommender system,” Bioinformatics, vol. 34, no. 22, pp. 3907–3914, 2018.
-  C. Qian, N. D. Sidiropoulos, M. Amiridi, and A. Emad, “From gene expression to drug response: A collaborative filtering approach,” in ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019, pp. 7465–7469.
-  F. Zhang, M. Wang, J. Xi, J. Yang, and A. Li, “A novel heterogeneous network-based method for drug response prediction in cancer cell lines,” Scientific reports, vol. 8, no. 1, p. 3355, 2018.
-  U. McDermott, “Cancer cell lines as patient avatars for drug response prediction,” Nature genetics, vol. 50, no. 10, p. 1350, 2018.
-  Z. Bar-Joseph, A. Gitter, and I. Simon, “Studying and modelling dynamic biological processes using time-series gene expression data,” Nature Reviews Genetics, vol. 13, no. 8, p. 552, 2012.
-  A. Fukushima, M. Sugimoto, S. Hiwa, and T. Hiroyasu, “Elastic net-based prediction of ifn- treatment response of patients with multiple sclerosis using time series microarray gene expression profiles,” Scientific reports, vol. 9, no. 1, p. 1822, 2019.
-  Y. Luan and H. Li, “Clustering of time-course gene expression data using a mixed-effects model with b-splines,” Bioinformatics, vol. 19, no. 4, pp. 474–482, 2003.
-  S. E. Baranzini, P. Mousavi, J. Rio, S. J. Caillier, A. Stillman, P. Villoslada, M. M. Wyatt, M. Comabella, L. D. Greller, R. Somogyi et al., “Transcription-based prediction of response to ifn using supervised computational methods,” PLoS biology, vol. 3, no. 1, p. e2, 2004.
-  T.-h. Lin, N. Kaminski, and Z. Bar-Joseph, “Alignment and classification of time series gene expression in clinical studies,” Bioinformatics, vol. 24, no. 13, pp. i147–i155, 2008.
Y.-H. Taguchi, “Identification of candidate drugs using tensor-decomposition-based unsupervised feature extraction in integrated analysis of gene expression between diseases and drugmatrix datasets,”Scientific reports, vol. 7, no. 1, p. 13733, 2017.
-  Y. Li and A. Ngom, “Non-negative matrix and tensor factorization based classification of clinical microarray gene expression data,” in 2010 IEEE international conference on bioinformatics and biomedicine (BIBM). IEEE, 2010, pp. 438–443.
-  L. Zhao and M. J. Zaki, “Tricluster: an effective algorithm for mining coherent clusters in 3d microarray data,” in Proceedings of the 2005 ACM SIGMOD international conference on Management of data. ACM, 2005, pp. 694–705.
-  D. Jiang, J. Pei, M. Ramanathan, C. Tang, and A. Zhang, “Mining coherent gene clusters from gene-sample-time microarray data,” in Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining. ACM, 2004, pp. 430–439.
-  K. Ø. Mikalsen, F. M. Bianchi, C. Soguero-Ruiz, and R. Jenssen, “Time series cluster kernel for learning similarities between multivariate time series with missing data,” Pattern Recognition, vol. 76, pp. 569–581, 2018.
-  K. Moorthy, M. Saberi Mohamad, and S. Deris, “A review on missing value imputation algorithms for microarray gene expression data,” Current Bioinformatics, vol. 9, no. 1, pp. 18–22, 2014.
-  S. Dudoit, J. Fridlyand, and T. P. Speed, “Comparison of discrimination methods for the classification of tumors using gene expression data,” Journal of the American statistical association, vol. 97, no. 457, pp. 77–87, 2002.
-  Y. Li, A. Ngom, and L. Rueda, “Missing value imputation methods for gene-sample-time microarray data analysis,” in 2010 IEEE Symposium on Computational Intelligence in Bioinformatics and Computational Biology. IEEE, 2010, pp. 1–7.
-  A. Kapur, K. Marwah, and G. Alterovitz, “Gene expression prediction using low-rank matrix completion,” BMC bioinformatics, vol. 17, no. 1, p. 243, 2016.
-  N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos, “Tensor decomposition for signal processing and machine learning,” IEEE Transactions on Signal Processing, vol. 65, no. 13, pp. 3551–3582, July 2017.
-  Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” SIAM Journal on imaging sciences, vol. 6, no. 3, pp. 1758–1789, 2013.
-  W. A. Craig, “Pharmacokinetic/pharmacodynamic parameters: rationale for antibacterial dosing of mice and men,” Clinical infectious diseases, vol. 26, no. 1, pp. 1–10, 1998.
-  R. Parry, W. Jones, T. Stokes, J. Phan, R. Moffitt, H. Fang, L. Shi, A. Oberthuer, M. Fischer, W. Tong et al., “k-nearest neighbor models for microarray gene expression analysis and clinical outcome prediction,” The pharmacogenomics journal, vol. 10, no. 4, p. 292, 2010.