A recursive framework for predicting the time-course of drug sensitivity

The biological processes involved in a drug’s mechanisms of action are oftentimes dynamic, complex and difficult to discern. Time-course gene expression data is a rich source of information that can be used to unravel these complex processes, identify biomarkers of drug sensitivity and predict the response to a drug. However, the majority of previous work has not fully utilized this temporal dimension. In these studies, the gene expression data is either considered at one time-point (before the administration of the drug) or two time-points (before and after the administration of the drug). This is clearly inadequate in modeling dynamic gene–drug interactions, especially for applications such as long-term drug therapy. In this work, we present a novel REcursive Prediction (REP) framework for drug response prediction by taking advantage of time-course gene expression data. Our goal is to predict drug response values at every stage of a long-term treatment, given the expression levels of genes collected in the previous time-points. To this end, REP employs a built-in recursive structure that exploits the intrinsic time-course nature of the data and integrates past values of drug responses for subsequent predictions. It also incorporates tensor completion that can not only alleviate the impact of noise and missing data, but also predict unseen gene expression levels (GEXs). These advantages enable REP to estimate drug response at any stage of a given treatment from some GEXs measured in the beginning of the treatment. Extensive experiments on two datasets corresponding to multiple sclerosis patients treated with interferon are included to showcase the effectiveness of REP.

Prediction of drug response based on patients' clinical and molecular features is a major challenge in personalized medicine. A computational model that can make accurate predictions can be used to identify the best course of treatment for patients, or to identify therapeutic targets that can overcome drug resistance 1 . Considerable efforts have been made 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 molecular 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][8][9] . 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, Costello et al. 2 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., GI50-the concentration for 50% of maximal inhibition of cell proliferation from a single time-point. In practice, however, for many diseases (e.g. cancers) the response to a drug changes over time due to various reasons such as the development of drug resistance or changes in the progress of the disease. To capture such changes at a molecular level, a collection of temporal gene expression profiles of samples over a series of time-points during the course of a biological process is necessary to provide more insights than a single (or two) time-point(s) 10 . 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 levels (GEXs) 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. For example, time-course gene expression data can be used to not only identify longitudinal phenotypic markers 11,12 , but also assess the association between the time course molecular data and cytokine production in this HIV trial 13 and predict drug response during a treatment 8,14 . In 15 , 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 Methods Figure 1 sketches the idea behind the proposed REP algorithm, where the Fig. 1a-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. 1a. In the following, we explain them in more detail.
Pre-processing. 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 been previously suggested for handling missing values, such as median-imputation 21 and nearest neighbor imputation 22,23 . Instead, 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 GEX data 24 . We note that our low-rank tensor model suggests three factors that uniquely determine the values of GEXs, 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 GEX over time from a set of initial GEX measurements; these estimated values are then used to predict the time-course of drug response.
Towards this goal, we first assume (For high-enough but finite F, any patient-gene-time dataset can be expressed this way. See 25 for a tutorial overview of tensor rank decomposition.) that each GEX is represented as a summation of F triple products from the latent factors of patient, gene and time, respectively. Let us denote g ijk as the jth GEX of patient i recorded at time k. Based on our assumption, we have where a if , b jf and c kf are the latent factors of patient, gene and time, respectively. Suppose that there are J genes measured over K time points. By varying the indices j and k in (1), the expression of the genes in all the timepoints in patient i can be represented as represents a diagonal matrix holding the ith row of A as the main diagonal, which is a latent representation of the ith patient. We use a if to represent the (i; f)-entry of A , b jf to represent the (j; f)-entry of B and c kf to represent the (k; f)-entry of C. Assume that there are I patients in the training set. After collecting {G 1 , . . . , G I } , we stack them in parallel along the patient-axis, which results in a GEX tensor that takes the form of where • is the outer product and a f is the fth column of A , and likewise for b f and c f . Here, we assume that G is the noiseless GEX data and X is the corresponding noisy data with missing values. The relationship between G and X is described as where N is the noise in the data, is the index set of the observed GEXs in X , and P 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., B and C ) are identical for different patients, and the variability among patients is captured by A . In other words, given B and C , each row of the patient factor matrix A uniquely determines the GEXs of the corresponding patient. As we will see later, our model is able to predict unseen GEXs, which also enables to prescreen the drug response for different stages of a treatment.
Assuming non-negative GEXs, (Due to some preprocessing steps such as z-score normalization, the GEX values can be negative. To facilitate our method, we undo these preprocessing steps or use the raw dataset.) we can use non-negative tensor factorization to compute missing GEX values: where many sophisticated algorithms are applicable to optimize (5), e.g., block coordinate descent 25,26 . Intuitively, (5) seeks to identify the lowest rank solution G that best matches the observations X . The regularization terms are added to further encourage low rank and prevent over-fitting. When (5) is solved, we complete the GEX data through where c contains the indices of missing values in X. www.nature.com/scientificreports/ Training. The effects of drugs are usually cumulative over time 27 , 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. Figure 1c shows an overview of REP's pipeline, where drug responses {y(0), . . . , y(t − 1)} in the previous time stages are integrated with the gene expression information z t for predicting the current response y(t). Here, we accumulate the historical responses by concatenating them into a new vector as which is then fed back as an input feature for subsequent drug response prediction. Therefore, at time t, the output of the predictor depends not only on the GEX at that time point, but also the previously observed drug responses. We always insert the drug response from the most recent time point into the first element of ỹ(t) , so that the model can capture a sense of time, and learn from recent/emerging trends in drug response. For the ith patient at time t, we concatenate ỹ i,t and the corresponding GEX vector z i,t together, where ỹ i,t denotes the historical responses of patient i at time t. We then pass z i,t and y i,t to a predictor f (·) to predict drug response, i.e., where f (·) 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, r(·) is a regularizer that imposes a certain structure on θ , and ≥ 0 is a regularization parameter. In the literature, popular regularizers include r(θ) = �θ � 2 2 , r(θ) = �θ � 0 , r(θ) = �θ � 1 and r(θ) = 1 + (θ) , i.e., the indicator function of the non-negative orthant.
Our main idea is to feed back the historical drug responses and then combine them with GEX values to predict the drug response in the future. This is the major difference between our method and the state-of-theart algorithms: prior art ignored the previous drug response outputs. Therefore, the training set for our method is created in a slightly different way. Recall that at each time point, we stack the historic drug responses into a vector. For any patient in the training set, we can further concatenate all such prior response vectors for K time points together, which yields a feedback matrix for patient i as It is also important to mention that our method can predict either binary or non-binary drug responses, e.g., continuous values. When the drug response is binary, the predictor f (z i,t ,ỹ i,t ) will typically be a classifier. When the drug response is continuous, f (z i,t ,ỹ i,t ) will be a regression algorithm.
It is important to note that our approach is really a framework that is applicable no matter what is the choice of the final classification or regression algorithm. Nevertheless, for the purposes of exemplifying and illustrating the merits of our proposed framework, we are particularly interested in support vector machines (SVM and SVR, for classification and regression, respectively), which have shown promising performance in this type of task (cite a few past papers using SVM for this task here). We set θ = [u T , v T ] T and ℓ(·) to be the hinge loss, resulting in where b is the intercept, C denotes a convex set such as ℓ 1 -ball for feature selection and ρ represents 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. In SVMs, we penalize the two-norm of the linear weights equally-the implicit assumption being that features have similar powers. In our context, however, the GEX values are much larger than the drug response labels which are either 1 or −1 . As a result, the GEX values are likely to end up playing a more significant role in the prediction-simply because we cannot scale the labels up to any meaningful level, due to the regularization y 1,1 y 1,2 · · · y 1,K y 2,1 y 2,2 · · · y 2,K . . . www.nature.com/scientificreports/ term. To compensate for this imbalance, in the above formulation, we introduce fixed weight ρ in the cost function. In practice, we recommend to choose a relatively large ρ.
Drug response prediction. Our method can predict the drug response values for a new patient at any time point. Specifically, given the GEXs of a new patient at time t, i.e., x(t) , we first check if there are missing values. If so, we employ the factors B and C to complete x(t) . Let us denote ¯ and ¯ c as the sets of indices of the observed and missing elements in x(t) . According to our model in (1), x(t) can be uniquely determined by B , C and an unknown vector a -a latent representation of this new patient. Thus, for the expression level of the jth gene at time t, we have where n j is the additive noise which is assumed as Gaussian distributed, ⊙ is the Khatri-Rao (column-wise Kronecker) product, and B(j, :) and C(t, :) denote the tth row of B and C , respectively.
Since B and C are known, the problem of estimating a can be formulated as which is a non-negative 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 x(t) should be ≥ F . The GEX vector of the patient is then estimated as which leads to a completed GEX vector The vector z(t) together with the cumulated historical drug response ỹ(t) , are the input data for our predictor f (·) . We estimate the drug response of this patient at time t via It is crucial to mention that in some cases, there might be missing labels in the testing set, such that ỹ(t) cannot be constructed. To handle this scenario, we can use the predicted labels instead of the missing ones to construct ỹ(t) . More specifically, we start from t = 0 and predict y(0), which is used to construct ỹ(1) . Then we use ỹ(1) to predict the response at t = 1 , so on and so forth.
Predicting unseen GEXs. 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 GEXs of all time points up to the one of interest a priori, which is impossible in practice. In this subsection, we provide an efficient solution that allows to predict the unseen GEXs.
Recall that in our model, the GEX of a patient is determined by three factors, i.e., the latent representation of patient-a , the time evolution factor-B and the gene factor matrix-C , where a is different for patients, and needs to be estimated for the new patient. On the other hand, B and C 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 a from the initial GEXs of the new patient. We can simply substitute t = 1 in (14) to find â . Finally, the GEXs for the remaining time points are estimated as Now we have estimated the unseen GEXs for t ≥ 2 , which allows us to predict drug response values for the whole duration of the treatment. We start from x(1) and estimate the drug response for t = 1 as where ỹ(1) = 0 . When ŷ(1) is available, we set ỹ(2) =ŷ(1) . With the GEX estimate x(2) from (18), we can predict ŷ(2) = f (x(2),ỹ(2)) , and so forth for the other time points.

Remark 1
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. ).

Experiments
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 examine two tasks: classification on binary drug response and regression on continuous drug response.
Dataset. We consider two datasets to evaluate the performance of our method.
• The first dataset used is the interferon (IFN)-β time-course dataset which is available in the supplementary of 15 . The dataset was collected from 53 Multiple Sclerosis (MS) patients who received IFN-β treatment for 2 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 seven stages (i.e., timepoints) of the treatment, where the time lag between two adjacent time points was 3 months in the first year, and 6 months in the second year. The responses to the therapy were measured at each time point using the expanded disability status scale (EDSS) which is a method of quantifying disability in multiple sclerosis and monitoring changes in the level of disability over time 28 . Note that EDSS values in the dataset are between 0 and 7, where the higher EDSS values reflect more severe disability. Except for the EDSS at the initial time point, the others were measured after the IFN-β injection at each time point. Therefore, we focus on the prediction of EDSS after t = 1 . In addition to EDSS, whether a patient had good or poor response to each treatment was also recorded, for each patient-and this is the indicator that we seek to predict in our classification experiments. The percentage of good patient responses to individual treatments was 58.5%; the remaining 41.5% responses were poor, on average, across the patient population considered. There are also missing values in this dataset, where most missing values were caused by the absence of patients at some stages. Only 27 patients had records for all stages, while the other 26 patients missed at least one stage, which resulted in the entire GEXs as well as the drug response at that stage being missed. In the following experiments, unless specified otherwise, we employ the 27 full records to evaluate the performance of algorithms, where the final GEX data is of size 27 × 7 × 76 and the response data is of size 27 × 7. • The second dataset is from a Gene Expression Omnibus (GEO) record GSE24427 29 also corresponding to MS.
In the dataset, there are 16 female patients and 9 male patients who received IFN-β therapy for 24 months. During the treatment, the RNA expression values were measured five times: at baseline (before first IFN-β injection), 2 days (before second IFN-β injection), and 1 month (before month-1 IFN-β injection), 12 months (before month-12 IFN-β injection), and 24 months (before month-24 IFN-β injection), respectively. The EDSS values were measured four times: at baseline, after 1 year, 2 years and 5 years of the initial injection, respectively. We use (1) the RNA expressions measured before month-12 to predict the EDSS measured after 1 year treatment, (2) the gene expressions measured before month-24 to predict the EDSS measured after 2 years treatment. There are 47,522 gene probes in this dataset. We employed the python package mygene (mygene: https ://mygen e.info/) to map the probes to gene names, which yielded 19,565 gene names.
Unlike the first dataset, we do not have binary drug response in the second dataset. Therefore, we focus on the prediction of binary drug response on the first dataset (whether or not a patient will have good or poor response), while the prediction of EDSS for both datasets is viewed as a regression task, because EDSS is an ordinal variable (predicting a 6 as a 7 is better than predicting a 6 as a 3; thus mean absolute error and room mean squared error make sense as performance metrics).

Methods for comparison.
We examine the predictive ability on the prediction of binary drug response and ordinal EDSS response. For the binary case, we apply a number of classifiers including two linear models (EN-LR 8 and SVM), one nonlinear model (K-nearest neighbors (KNN) 30 ), and a probabilistic graphical model (discriminative loop hidden Markov model (dl-HMM) 14 ) 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 timecourse gene expression data, while SVM and KNN are widely used classification methods. For ordinal prediction, we implement Elastic Net, Support Vector Regression (SVR) with radial basis function (rbf) krenel, Random Forest and KNN on the two datasets. All methods are implemented via the Python sklearn package with version 0.0. We use default settings for Elastic Net and SVR algorithm. For Random Forest, we set the number of trees in the forest to 20. For KNN, we set the number of neighbors to 10. For each dataset, we first create two versions of training and testing sets, where one is with the drug response feedback described in Fig. 1 and the other one is without feedback. We use REP-ElasticNet, REP-SVR, REP-RandomForest and REP-KNN to denote the respective algorithms with drug response feedback.  where M is the number of samples in the testing set, N is the number of Monte-Carlo tests, y m denotes the ground-truth drug response of the mth testing sample and ŷ m is its estimate with y m (t) and ŷ m (t) being their values at time t, respectively.
We report performance for all methods by using the same training, validation and testing sets. Specifically, we employ leave-one-out cross validation (LOO) for testing, where at each fold, we split 27 patients into a training set with 26 patients and a testing set with one patient. We then hold the testing set and randomly split the training set into two parts, where the first part has 25 patients and the second part has one patient, i.e., the validation set. We train models on the first part and tune hyper-parameters on the second part. Note that for each algorithm, we select its best hyper-parameters-those that yield the highest prediction accuracy on the validation set. Finally, we apply the selected model to the testing set. For a fair comparison, in all experiments, we apply the same missing value imputation method to all algorithms.
For REP-SVM, its hyper-parameters include and ρ in (12)

Results. Parameter selection for tensor completion.
We first study how the hyper-parameters F and µ in (5) affect the prediction performance. The percentage of missing values in GEXs is fixed at 5%. We vary F from 2 to 5 and µ from 0.01 to 100, and report ACC of REP-SVM on the classification task, i.e., predicting good or poor responder. The ACC is calculated using leave-one-out (LOO) cross-validation, where in each fold, we select one patient's record as a testing set that contains a 7 × 76 GEX matrix and a response vector with length 7, while the remaining 26 records are assigned to the training set. It can be seen in Fig. 2 that F ≤ 3 produces better results than F ≥ 4 in general, especially when µ ≥ 1.
Performance evaluation on binary drug response. We now compare the performance of the five algorithms in terms of prediction accuracy and AUC. In the raw data, the missing values in X are only 0.23%. For all the methods, the missing GEXs were completed using non-negative tensor completion in Section II-A, where we set F = 3 and µ = 1 . It can be seen in Table 1 (i.e., the rows with % miss as 0.23) that REP-SVM achieves a higher prediction accuracy compared to other methods and its performance is followed by the SVM and EN-LR algorithms. The KNN and dl-HMM algorithms have relatively low accuracy. We note that REP-SVM takes a similar formulation as the SVM. However, REP-SVM is 2.5% more accurate than SVM, which implies that the recursive structure in REP-SVM is helpful in improving the prediction accuracy.
We sought to determine the effect of missing values on the performance of these methods. For this purpose, we randomly sampled the GEX data and hid the selected entries. As the percentage of missing values increases, all methods suffer performance loss, but REP-SVM's ACC and AUC remain the highest in all cases (see Table 1). We highlight that when the percentage of missing values is 20%, REP-SVM has ACC close to 0.872 and AUC greater than 0.941. EN-LR outperforms the classical SVM method in many cases. When the percentage of missing values increases, the performance of EN-LR and SVM drop significantly, while that of REP-SVM still remains at a high level. For example, when the percentage of missing increases to 15%, the ACC of EN-LR drops to 0.844 and that of SVM drops to 0.837, but ours is 0.887, which indicates that REP-SVM is more robust against missing values. We mention that in this experiment, the ratio of the positive and negative classes is 19/8. So the percentage of the positive class is about 70.4%. We found that most of the predicted labels of KNN were positive, meaning that it cannot distinguish the negative class. The scores produced by KNN were not good enough to separate the two classes. This is why KNN yields seemingly reasonable accuracy but low AUC.
In this example, we evaluate the performance comparison on all patients in Dataset1. As we have mentioned before, there are 26 patients that do not have seven time points records, so they cannot be used in the training step of REP-SVM. Therefore, we use them for testing, and the training is based on the 27 patients with full temporal records, where we set F = 3 and µ = 1 . Note that in the testing data, approximately 18.1% of the GEX values and 22% drug response labels are missing. Here, the missing GEXs are completed through our non-negative tensor completion. We calculate the ACC and AUC based on the known drug response labels. The results are shown in Table 2. We see that REP-SVM outperforms the other competitors in both accuracy and AUC, where it achieves Scientific Reports | (2020) 10:17682 | https://doi.org/10.1038/s41598-020-74725-2 www.nature.com/scientificreports/ about 0.790 in ACC and 0.884 in AUC. EN-LR has the second best performance and is followed by SVM in terms of accuracy, but EN-LR has higher AUC than SVM. In this case, the distributions of missing values in the training and testing sets are very different, where the percentage of missing values in training set is about 0.23% but that in the testing set is about 18.1%. Recall that in Table 1, where the missing values were randomly assigned through a uniform distribution, REP-SVM, EN-LR and SVM have higher ACC and AUC than the results in Table 2, even though the percentage of missing values reaches 20%. This indicates that the distribution of missing values in training and testing sets may affect the performance of drug response predictors. We have shown that under the same completion algorithm, REP-SVM has better performance than the methods, but one may wonder if the same conclusion holds for other types of completion. To answer this question, we further compare REP-SVM, EN-LR, and SVM with mean, median, and KNN imputation. The results are shown  www.nature.com/scientificreports/ in Table 3. We see that all predictors with tensor completion achieve the highest ACC and AUC compared to the standard imputation methods. We also note that REP-SVM continues to have the best performance even when used with less sophisticated imputation methods. Figure 3 shows the top 20 genes selected by REP-SVM. Note that we run REP-SVM ten times on the 27 patients with full temporal records, average the weights corresponding to the genes and then rank the weights to generate the gene ranking. The genes IRF3, IRF4, IRF6 and IRF8 belong to the interferon regulatory transcription factor (IRF) family, which is critical for induction of type I (IFN-α/β ) and type III (IFN-β ) IFN expression [31][32][33] .
Performance evaluation on ordinal drug response. Next, we evaluate the performance of the proposed method on the prediction of EDSS. Since the number of genes is much larger than the number of patients, it is challenging to perform regression directly on such an underdetermined dataset. To handle this issue, we select the genes which are also in the first dataset. This results in 42 overlapped genes (see Supplementary Table S1). Therefore, in the following experiments, we only used the selected 42 genes. Note that all methods use the same genes for regression, so reducing the number of genes will not affect the fairness. The results are shown in Table 4, where the percentage of missing values of GEX in Dataset1 is 0.23% and that in Dataset2 is 0. In both datasets, all methods with drug response feedback outperform their respective version without feedback in terms of MAE and MSE. Table 4 shows the overall performance comparison by averaging all time points. One may also be interested in the performance of REP at each time point. In Supplementary Tables S2-S5,   The above experiments are all based on the actual GEX values with only a small portion of them imputed. In the last experiment, we assume that in the testing set, only the GEX values at the initial time point have been observed. The training set and the way of imputing the missing values therein are kept the same as in the previous experiment in Table 4. Our task is to predict future drug responses for the remaining time points at an initial time point. In this case, the GEX values for the future time points are entirely missing. To handle this issue, we first implement the tensor completion method to learn the latent factor B for GEX and another latent factor C corresponding to the T time points from the training set. Then according to (1), given the GEX vector of the ith patient at the initial time point, we can estimate his/her latent representation denoted by a i , by solving a nonnegative last squares problem, i.e., min a i ≥0 �g 0 − Bdiag(c 0 )a i � 2 2 , where g 0 is the observed GEX values at t = 0 and c 0 is the first row of the estimated temporal factor matrix C corresponding to t = 0 . This problem is convex and it can be optimally solved. Then we can substitute a i into (2) to predict the entire GEX matrix for the ith patient and we use the second to the last column of G for prediction, where the tth column in G i represents the predicted GEX values at time t. Furthermore, we do not have any information about the drug responses for the future except for the one at t = 0 , and thus unable to feed back the actual drug response at t ≥ 1 for subsequent prediction. We have mentioned previously that for such cases, we can feed back the predicted drug responses for future prediction if a previous drug response is unavailable. Therefore, in the testing step, we start from t = 1 where we can use the GEX values and EDSS value at t = 0 to predict t = 1 . Then we use this estimated drug response and the estimated GEX at t = 1 to predict t = 2 , so on and so forth. Table 5 shows the results, where the column 'Using estimated GEX' stands for using estimated GEXs and feeding back estimated drug responses for subsequent prediction, while 'Using actual GEX' stands for using actual GEXs and feeding back actual drug responses. It can be seen in Table 5 that predictors using the predicted GEX with the feedback of predicted drug responses have slightly worse performance than using the actual GEX

Conclusion
We studied the problem of drug response prediction for time-course gene expression data and presented a computational framework (REP) that: (1) has a recursive structure that integrates past drug response records for subsequent predictions, (2) offers higher prediction accuracy than several classical algorithms such as SVM and LR, (3) exploits the tensor structure of the data for missing GEX completion and unseen GEX prediction, (4) can predict drug response of different stages of a treatment from some initial GEX 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.