I Introduction
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 [1] analyzed 44 drug response prediction methods that employed gene expression profiles of breast cancer cell lines taken before treatment to predict doseresponse 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 timepoint. A collection of temporal gene expression profiles of samples over a series of timepoints during the course of a biological process can provide more insights than a single (or two) timepoint(s) [8]. Therefore, developing algorithms that can predict the drug response over time using timecourse gene data is of great interest.
With the advancement of gene sequencing technologies, collecting gene expression data over multiple timepoints 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 timecourse gene expression data. In
[9], the authors proposed an integrated Bayesian inference system to select genes for drug response classification from timecourse gene expression data. However, the method only uses the data from the first timepoint, and hence does not benefit from the additional temporal information. Lin et al.,
[10]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 timepoints
a 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 timecourse gene expression data contains the GELs of different patients over a series of time points, which can be indexed as patientgenetime and represented as a threedimensional tensor. Motivated by this, several tensor decomposition based algorithms have been proposed. For example, Taguchi [11] employed tensor decomposition to identify drug target genes using timecourse gene expression profiles of human cell lines. Li and Ngom [12] proposed a higherorder 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 classlabel information. Moreover, this method simply discards samples with missing values, causing unnecessary information loss.
Recently, Fukushima et al., [7]
developed an algorithm for joint gene selection and drug response prediction for timecourse data. The method uses ElasticNet (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 timecourse drug response prediction. We hypothesize that a patient’s drug response at a given timepoint 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 timecourse gene expression data and their drug response at previous timepoints. REP has a builtin recursive structure that exploits the intrinsic timecourse 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 timecourse gene expression data and leveraging identifiability of lowrank 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.
Ii Method
Fig. 1 sketches the idea behind the proposed REP algorithm, where the subfigures 1(a)–1(c) show the preprocessing, model training and prediction of our method, respectively. The tensor structure of timecourse gene expression data is shown in Fig. 1(a). In the following, we explain them in more detail.
Iia Preprocessing
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 preprocessing. Various methods have bene previously suggested for handling missing values, such as medianimputation
[16] and nearest neighbor imputation [17, 18]. We employ a lowrank tensor model to fit the timecourse 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 [19]. We note that our lowrank 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 timecourse of drug response.Towards this goal, we first assume^{1}^{1}1For highenough but finite , any patientgenetime dataset can be expressed this way. See [20] 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
(1) 
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 timepoints in patient can be represented as
(2) 
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 patientaxis, which results in a GEL tensor that takes the form of
(3) 
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
(4) 
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 GELs^{2}^{2}2
Due to some preprocessing steps such as zscore 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:(5) 
where many sophisticated algorithms are applicable to optimize (IIA), e.g., block coordinate descent [21, 20]. Intuitively, (IIA) seeks to identify the lowest rank solution that best matches the observations . The regularization terms are added to further encourage low rank and prevent overfitting. When (IIA) is solved, we complete the GEL data through
(6) 
where contains the indices of missing values in .
IiB Training
The effects of drugs are usually cumulative over time [22], i.e., drug doses taken in the past will affect the current response. This implies that the drug response in the past timepoints 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
(7) 
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
(8) 
where denotes the accumulated historical responses of patient at time , is the predictor which can be trained by minimizing the following cost function
(9) 
where contains the parameters of the predictor,
is the loss function of a classifier such as hinge loss and crossentropy 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
(10) 
where is the intercept, denotes a convex set such as
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 “musthave” feature. Therefore, the hyperparameter is critical for the performance of REP. In practice, we recommend to choose a relatively large in an incremental fashion.Remark 1
The major difference between (IIB) and the standard SVMbased 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:
(11) 
where . Finally, the features in the training set is formed by concatenating in along the geneaxis as shown in the leftbottom corner of Fig. 1(b), and the training labels are
(12) 
IiC 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(13) 
where
is the additive noise which is assumed as Gaussian distributed,
is the KhatriRao (columnwise Kronecker) product, and and denote the th row of and , respectively.Since and are known, the problem of estimating can be formulated as
(14) 
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
(15) 
which leads to a completed GEL vector
(16) 
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
(17) 
IiC1 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 timepoints in the future from the beginning of a treatment. This requires to know the GELs of all time points up to the timepoint 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
(18) 
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
(19) 
where . When is available, we set . With the GEL estimate from (18), we can predict , and so forth for the other time points.
Remark 2
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.
Iii Simulations
In this section, we provide some numerical experiments to showcase the effectiveness of REP for drug response prediction from timecourse gene expression data. We apply a number of drug response prediction methods including two linear models (ENLR [7] and SVM), one nonlinear model (
nearest neighbor (KNN)
[23]), and a probabilistic graphical model (discriminative loop hidden Markov model (dlHMM) [10]) to realworld timecourse data. We did not include SVM with nonlinear kernels (e.g. Gaussian), since its performance was inferior compared to the linear kernel. Note that ENLR and dlHMM were specifically designed for prediction of drug response values based on timecourse gene expression data, while SVM and KNN are widely used classification methods.Iiia Dataset
The dataset used is the interferon (IFN) timecourse dataset which is available in the supplementary of [9]. 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 preselected genes over 7 stages (i.e., timepoints) 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 .
IiiB 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 27fold 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 hyperparameters of ENLR, SVM and KNN are tuned through 5fold 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 hyperparameters include the rank , and in (IIB), which are selected from and . The SVM method solves the following problem:
(20) 
where its parameter is tuned from ; for ENLR, we set which is a hyperparameter 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, ENLR, SVM and KNN in Python 3.7. Since the authors of dlHMM 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 [10].
IiiC Results
Metrics  miss  REP  ENLR  SVM  KNN  dlHMM 

ACC  5  94.8  90.5  88.4  70.4  65.1 
10  92.6  88.4  88.7  69.3  59.6  
15  89.4  85.2  81.6  64.3  59.6  
20  88.5  84.7  81.5  66.8  55.6  
AUC  5  96.9  93.0  91.2  50.2  53.9 
10  96.8  90.1  88.6  49.4  54.7  
15  95.4  87.4  84.3  51.9  54.7  
20  93.2  88.1  87.1  52.3  56.4 
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 IIA. It can be seen that REP achieves a higher prediction accuracy compared to other methods and its performance is followed by the ENLR and SVM algorithms. The KNN and dlHMM algorithms do not perform well. We observed that dlHMM was very sensitive to the data, and it frequently produced notanumber (NaN) values. Particularly for patients with ID 1117161, 1215133 and 995640, the original codes of dlHMM always produced NaN of the three patients’ responses. This is the main reason why dlHMM 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%. ENLR outperforms the classical SVM method in many cases. When the percentage of missing values increases, the performance of ENLR 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 ENLR 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 IIC1 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 timepoints 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.
Iv Conclusion
We studied the problem of drug response prediction for timecourse 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 timecourse data.
Conflict of Interest: none declared
References
 [1] 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.
 [2] C. Suphavilai, D. Bertrand, and N. Nagarajan, “Predicting cancer drug response using a recommender system,” Bioinformatics, vol. 34, no. 22, pp. 3907–3914, 2018.
 [3] C. Qian, N. D. Sidiropoulos, M. Amiridi, and A. Emad, “From gene expression to drug response: A collaborative filtering approach,” in ICASSP 20192019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019, pp. 7465–7469.
 [4] F. Zhang, M. Wang, J. Xi, J. Yang, and A. Li, “A novel heterogeneous networkbased method for drug response prediction in cancer cell lines,” Scientific reports, vol. 8, no. 1, p. 3355, 2018.
 [5] U. McDermott, “Cancer cell lines as patient avatars for drug response prediction,” Nature genetics, vol. 50, no. 10, p. 1350, 2018.
 [6] Z. BarJoseph, A. Gitter, and I. Simon, “Studying and modelling dynamic biological processes using timeseries gene expression data,” Nature Reviews Genetics, vol. 13, no. 8, p. 552, 2012.
 [7] A. Fukushima, M. Sugimoto, S. Hiwa, and T. Hiroyasu, “Elastic netbased 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.
 [8] Y. Luan and H. Li, “Clustering of timecourse gene expression data using a mixedeffects model with bsplines,” Bioinformatics, vol. 19, no. 4, pp. 474–482, 2003.
 [9] 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., “Transcriptionbased prediction of response to ifn using supervised computational methods,” PLoS biology, vol. 3, no. 1, p. e2, 2004.
 [10] T.h. Lin, N. Kaminski, and Z. BarJoseph, “Alignment and classification of time series gene expression in clinical studies,” Bioinformatics, vol. 24, no. 13, pp. i147–i155, 2008.

[11]
Y.H. Taguchi, “Identification of candidate drugs using tensordecompositionbased unsupervised feature extraction in integrated analysis of gene expression between diseases and drugmatrix datasets,”
Scientific reports, vol. 7, no. 1, p. 13733, 2017.  [12] Y. Li and A. Ngom, “Nonnegative 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.
 [13] 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.
 [14] D. Jiang, J. Pei, M. Ramanathan, C. Tang, and A. Zhang, “Mining coherent gene clusters from genesampletime microarray data,” in Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining. ACM, 2004, pp. 430–439.
 [15] K. Ø. Mikalsen, F. M. Bianchi, C. SogueroRuiz, 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.
 [16] 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.
 [17] 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.
 [18] Y. Li, A. Ngom, and L. Rueda, “Missing value imputation methods for genesampletime microarray data analysis,” in 2010 IEEE Symposium on Computational Intelligence in Bioinformatics and Computational Biology. IEEE, 2010, pp. 1–7.
 [19] A. Kapur, K. Marwah, and G. Alterovitz, “Gene expression prediction using lowrank matrix completion,” BMC bioinformatics, vol. 17, no. 1, p. 243, 2016.
 [20] 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.
 [21] 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.
 [22] 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.
 [23] R. Parry, W. Jones, T. Stokes, J. Phan, R. Moffitt, H. Fang, L. Shi, A. Oberthuer, M. Fischer, W. Tong et al., “knearest neighbor models for microarray gene expression analysis and clinical outcome prediction,” The pharmacogenomics journal, vol. 10, no. 4, p. 292, 2010.
Comments
There are no comments yet.