One of the fundamental challenges in modern medicine is that over time, many diseases develop resistance to drugs. For infectious diseases, this often happens as the microbial species evolve, and as a result, drug resistant species replace those that are not resistant [3]. For genetic disorders, such as cancer, drug resistance arises naturally as each cancer is a result of unique mutations for individuals, while drugs target the condition of an average patient belonging to a population of cancer patients of varying mutation types [25]. New mutations yield cancer cells, which may not fit into the average cancer profile, leading to drug resistant cancer cells [12].
The e↵ectiveness of anti-cancer drugs is often measured by their e↵ects on the growth and death of cancer cells relative to control samples. Two standard metrics for measuring drug potency and e cacy are the half-maximal inhibitory concentration (IC50) and the fraction of viable cells at the highest drug concentration (Emax) [11]. One way to combat drug resistance is to use combinations of drugs [2,13,17]. Doing so aims to increase the e cacy of the drugs by taking advantage of drug-drug interactions, and thus yields an e↵ect on the targeted disorder more than the e↵ect of the individual drugs if they were used independently. This improvement in drug e cacy happens for several reasons. First, the drug combination reduces the chance of developing resistance, for when cancer cells are exposed to multiple drugs, it is more di cult for them to develop resistance to all of them at the same time. This is because, a synergic drug combination accesses context specific multi-target mechanisms, and makes it harder for a cancer cell to develop resistance to multiple mechanisms of action [14]. Secondly, it reduces the risk of toxicity side e↵ects as synergistic combination often reduces the dosage needed for each drug [14,16,22].
There are several databases which maintain the interactions among drugs. Among them, GDSC 2 hosts the e cacy of approximately 2000 clinically relevant drug combinations. These combinations encompass 125 characterized breast, colorectal, and pancreatic cancer cell lines [13]. The synergy between drugs is rare and highly context-dependent. However, when the targeted drugs are combined, they are more likely to be synergistic [13]. At the heart of the drug-drug interaction identification problem lies the fact that trying all possible drug-drug combinations on various cell lines and finding the right dosage of the drug-drug combination su cient to eliminate a targeted percentage of given cancer cell lines is not feasible. This is because there are many drug-drug combinations, let alone the amounts of cell lines to be tested for each drug. This challenge necessitates developing computational strategies for estimating the right drug-drug combinations to be e↵ective on the targeted cell lines. Many prediction methods for single drug response are developed with e↵orts of National Cancer Institute (NCI) and AstraZeneca-Sanger Drug Combination Prediction Dialog for Reverse Engineering Assessments and Methods (DREAM) Challenge [7]. Some of them are based on matrix factorization [26], recommender system based on cell line similarities [23], drug/cell line network [27] or manifold learning [1]. As we include drug-drug combinations, the problem size increases, bringing further challenges. Existing methods for predicting drug-drug interactions consider this data as a set of independent samples [9,28]. The data available for drug-drug interactions however varies both in terms of quantity and distribution across di↵erent drugs as well as cell lines. This inherent nature of variability of drugdrug interaction data introduces substantial disadvantages to existing methods in terms of their accuracy and e ciency.
Our contributions. In this paper, we consider the problem of predicting drugdrug interactions. We develop a novel algorithm named %0AC80; 814AB to eciently and accurately predict the response of a given cell line to a given pair of drugs. More specifically, assume that we are given an incomplete three dimensional dataset, where the first dimension corresponds to a drug (called library drug), the second dimension also corresponds to a drug called anchor drug (i.e., a drug with a specific dosage,) and the third dimension corresponds to cell lines. The value at each entry of this tensor denotes IC50 value of the library drug for that cell line when that library drug is combined with the corresponding anchor drug. We aim to estimate the missing values in this dataset accurately and e ciently. Our method takes advantage of the topological properties of the three dimensional drug-drug interaction data. Fiber structure (i.e., the vector consisting of all values in the given 3D matrix, tensor, when the coordinates of two of the dimensions are fixed) is a common topological pattern in the missing data [5,19,24]. Drug-drug interaction data, however, does not have a uniform fiber structure as some values may be known and others may be missing in each fiber. We call such patterns partial fibers. Our method allows us to learn from multiple dimensions of the given data. Given 3D drug-drug interaction data, at each iteration, our algorithm randomly picks one of the three dimensions and learns from the partial fibers along that dimension. Our method provides a simultaneous update for all the samples along the same partial fiber. Our results demonstrate that %0AC80; 814AB e ciently and accurately predicts the drug responses for the cell lines and it shows superior performance as compared to state-of-the-art methods in terms of both accuracy and running time performance. It benefits from the intrinsic structure of the data for leading to an e cient per-iteration update performing a faster and more accurate imputation of drug-drug interactions.
In this section, we describe our %0AC80; 814AB method for e cient and accurate imputation of the missing drug-drug interactions data. In Section 2.1, we introduce our notations and the preliminary concepts of tensor completion needed to explain our algorithm. In Section 2.2, we describe our %0AC80; 814AB algorithm in detail.
Here, we introduce the notation that we use for describing the core concepts required to understand our algorithm.
Baseline tensor completion. The three-dimensional tensor completion problem aims to estimate the missing values in a given tensor, such that a nonempty strict subset of the entries of this tensor is known apriori. The objective here is to find three-factor matrices such that their inner product produces a tensor of the same size, and the values at the tensor locations corresponding to the apriori known values are as close as possible to the given known values. Once such a tensor is produced, this technique estimates the values generated in the positions corresponding to the missing values as the generated values at those positions of the new tensor.
Let us denote a three-dimensional tensor with rows, columns and depth with ^2 R ⇥ ⇥ , and the value of the entry in this tensor at index (8, 9, :) with G 8, 9,: . We denote the rank of this tensor as ' (' min{ , , } is a positive integer), and the three factor matrices as G 2 R ⇥' , H 2 R ⇥' , I 2 R ⇥' . We denote the 8th row of the factor matrix G with a 8 , the 9th row of the factor matrix H with b 9 and the :th row of the factor matrix I with c : . Similarly, we represent the value at index (8, A) of the factor matrix G with a 8A , the value at index ( 9, A) of the factor matrix H with b 9A and the value at index (:, A) of the factor matrix I with c :A . A nonempty subset of the entries in ^are known. We denote the set of indices (8, 9, :) for which the values of G 8, 9,: are available with ⌦. Using the notation above, the traditional CANDECOMP/PARAFAC (CP) formulation of tensor completion [6,10]
(1)
In ( 1), we present the optimization problem to find the factor matrices G, H, I that are going to be used to complete the missing entries in ^. Classic CP tensor completion algorithm iterates over the available tensor entries (⌦) one by one (8(8, 9, :) 2 ⌦) and updates the corresponding entries of the three factor matrices G, H, I [4]. Once the objective (1) converges, the complete data matrix can be constructed using the learned factor matrices.
Tensor completion with fibers. The distribution of the missing entries over the data tensor ^is not always uniform across the three dimensions. This provides opportunities to benefit from the structure of the missing data. A common structure for the missing data is fiber. A fiber is the vector of all index values (8, 9, :) when the values of any two of these three variables are fixed. For instance, given the values of 8 and 9, fiber is the vector of all the values along the :th dimension [15]. Thus, a fiber in a three dimensional tensor may consist of a varying index in any of the three dimensions. Recent research shows that the optimization function ( 1) can be computed more e ciently when the missing values of the tensor are along the same dimension (i.e., they make up a fiber) as multiple entries of the factor matrices can be computed in one iteration [5].
In this paper, we formulate the drug-drug interaction prediction problem, as the tensor completion problem. The fundamental challenge here is the tensor corresponding to this data does not follow the fiber structure we explain above. To solve this problem, we define a multi-dimensional notation which has more expressive power than the classic tensor notation. We use the "*" symbol to represent the dimension at which the fiber extends. For instance, we denote the fiber which has fixed values 9 and : along the last two dimensions and varies in the first dimension with G ⇤ 9 : = [G 1, 9,: , G 2, 9,: , . . . , G , 9,: ] ) . Similarly, we denote the fiber which has fixed values 8 and : along the first and third dimensions and varies in the second dimension with G 8⇤: = [G 8,1,: , G 8,2,: , . . . , G 8, ,: ] ) . Finally, we denote the fiber which has fixed values 8 and 9 along the first and second dimensions and varies in the third dimension with G 8 9⇤ = [G 8, 9,1 , G 8, 9,2 , . . . , G 8, 9, ] ) . Assume that, regardless of the direction, for all the fibers, within a fiber, either all entries are missing or all entries are available. Let us represent the Hadamard (elementwise) product with ⌦. Let us represent the set of all available fibers along the first dimension with 1 . Similarly, let us represent the set of all available fibers along the second and third dimensions with 2 and 3 , respectively.
Using this notation, depending on the direction of the fiber, we reiterate the original fiber update formulation and rewrite (1) as follows. If the fiber is along the first dimension, we write minimize G,b 9 ,c :
(
Similarly, we rewrite the optimization function based on the second dimension and the third dimension to update the factor matrices, H and I. We update the factor matrices iteratively until the factor matrices represent the actual data with a small loss, i.e., until convergence.
Updating the factor matrices in fiber provides faster convergence than the classic CP tensor completion algorithm because the entire fiber can be updated simultaneously. However, this strategy works only if all the values in a given fiber are missing. Thus, it cannot be applied to solve the drug-drug interaction prediction problem considered in this paper, as this data does not have this characteristic.
We first introduce the drug-drug interaction dataset and highlight the properties of the dataset. Then, we discuss our %0AC80; 814AB algorithm designed for drugdrug interaction dataset.
We use the drug-drug interaction data in the GDSC 2 dataset [13] for pairs of drugs. We refer to the first drug in a pair of drugs in this dataset as library drug. Unlike the library drug, this dataset associates the second drug with a specific dosage of that drug. We refer to the second drug at that dosage as anchor drug.
For the given dosage of the anchor drug, this dataset provides the IC50 value of the library drug when the two drugs are applied on a given cell line. We represent this data using a three-dimensional tensor with rows, columns and depth, ^2 R ⇥ ⇥ . Row 8, represents the 8th library drug. Column 9, represents the 9th anchor drug along with its dosage. Depth :, represents the :th cell line. In summary, each entry G 8, 9,: , denotes the dosage of the 8th library drug required Algorithm 1 %0AC80; 814AB, SGD for tensor completion using partial fibers Input: ^(the 3D incomplete data)
1: initialize G, H, and I 2: repeat 3:
Randomly select a dimension i, j, or k. (Assume that the dimension is i) 4:
Randomly draw an i-fiber, (⇤, 9, :) from ⇤, 9,: 5:
Initialize J 6:
Calculate the gradients 7:
to eliminate 50% of the cell line :, when it is combined with the pre-defined dosage of anchor drug 9. Figure 1 illustrates a representation of the drug-drug interaction data. Black and gray entries represent available and missing data, respectively.
We categorize fibers in two groups: complete fiber and partial fiber. We say that a fiber is complete if all the entries of that fiber are available or they are all missing. Otherwise, we call it a partial fiber. Orthogonal to this grouping, we further categorize fibers into 3 groups based on their directions; 8-fibers, 9-fibers, and :-fibers. -fibers represent all library drugs for a fixed anchor drug/amount, 9, and for a fixed cell line. For instance, in Figure 1(b), the 8-fiber represents all the library drugs for the anchor drug/amount in index 4, and the cell line in index 0. Similarly, the 9-fiber represents all the anchor drugs/amounts for the library drug in index 3, and the cell line in index 0. The :-fiber represents all the cell lines for the library drug in index 3, and the anchor drug/amount in index 4. In Figure 1(b), all three fibers are partial fibers since some values are available and some are missing in each. As the figure demonstrates, each fiber in the 3D representation may have a di↵erent number of missing entries. Also, the missing and present entries can be in any position as well as the (positive) number of positions.
We solve the optimization problem (1) using partial fiber updates. Algorithm 1 presents our %0AC80; 814AB method. The algorithm takes the incomplete tensor ^as input. We start by randomly initializing the factor matrices G, H and I (line 1). At each iteration, we randomly select a dimension among the three dimensions with equal probability (line 3). Upon selecting the dimension, we randomly choose one fiber among all the fibers in the selected direction with a uniform probability of being selected (line 4). In order to keep Algorithm 1 short, we write the subsequent steps assuming that we select the first dimension. However, it is important to note that any of the three dimensions can be chosen at each iteration. The subsequent steps still work the same by using the correct dimensions as indices in the formulation. Selecting the first dimen-sion implies that the current fiber will be selected among all available 8-fibers. Thus, a fiber along the first dimension for a randomly selected 9, : is the vector [G 1 9 : , G 2 9 : , . . . , G 9: ] ) . Next, we construct an ⇥ diagonal matrix J as follows. For all indices 8 2 1, 2, . . . , , if the value G 8 9: in the given tensor is known (i.e., (8, 9, :) 2 ⌦), then we set the diagonal entry J [8, 8] = 1. We set all the remaining entries in J to 0 (line 5 of Algorithm 1). The purpose of this diagonal matrix is to mask the unavailable entries within a fiber and to focus on only the available entries. We represent the missing data with the underscore " ". As an example, let us assume that = 4, and the first and third indices along dimension 8 are missing. We represent the ( 9,:)th fiber with G ⇤ 9 : = [ , G 2, 9,: , , G 4, 9,: ]. The corresponding diagonal matrix J 2 ' 4⇥4 for this example is, J =
2 6 6 6 6 6 6 4 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 3 7 7 7 7 7 7 5
.The next step is to perform parameter update, i.e. we need to update the indices of each factor matrix by using the information we learned from the available data. We use traditional stochastic gradient descent (SGD) as a baseline stochastic optimization algorithm [20]. Traditional SGD requires using one sample at a time. We redesign the traditional SGD to perform an e cient update by simultaneously learning from all the selected fiber entries. We design our SGD algorithm following from the observation that samples come in the form of x ⇤ 9 : fibers, which involves the entire matrix G and the variables b 9 , c : . We represent the objective function for the ( 9, :)th fiber with 5 We represent the step size in SGD with W. Then, the SGD algorithm takes the form as in line 10 of Algorithm 1. We repeat iterations (lines 2 to 11) until convergence, i.e. until the factor matrices are representative of the available data. The same convergence guarantees for the SGD apply for the partial fibers too, hence after convergence, we construct the complete data matrix using the learned factor matrices. Notice that it is also possible to adapt (2.2) to other stochastic gradient descent algorithms, such as Adam [18], Adagrad [8], and SPPA [4] through similar algebraic manipulations.
Advantages. There are three main advantages of partial fibers over the complete fiber setting:
-Data flexibility. If a specific entry in the tensor (8, 9, :) is missing, the complete fiber approach requires at least one of the fibers x ⇤ 9 : , x 8⇤: or x 8 9⇤ to be complete. This is an extremely restrictive constraint which is often not observed in drug-drug interaction datasets. %0AC80; 814AB provides a flexible condition which works for any distribution of missing data. We allow the data to have a di↵erent number of available entries, the entries can be in any order along the same dimension, we benefit by grouping the missing items together and perform the partial fiber update e ciently.
-Algorithmic flexibility. Complete fiber approach requires all fibers to be along the same dimension. %0AC80; 814AB, however, does not have this requirement. It selects the fiber dimension randomly with the same probability to benefit from every dimension equally.
-Running time. %0AC80; 814AB provides further improvement in running time as it performs a parameter update using all the samples within the same fiber.
It allows these samples to be non-sequential and of any size. The more number of available items would mean more simultaneous updates, which results in further time improvement. In the worst case, there will be only one such training sample which is equivalent to considering each sample one by one.
(Fully stochastic, non-fiber state).
Dataset. We use the ⌧⇡(⇠ 2 dataset which contains anchor drugs, library drugs, and their impact on a set of cell lines along with their IC50 values [13]. The dataset contains 160 anchor drugs, 65 library drugs, and 125 cell lines. Approximately 77% of the data is missing.
Experimental setup. We separate the available drug-drug interaction data into two disjoint sets: training and testing, for a given percentage ? 2 (0 : 100) as ?% and (100 ?)% of the available drug-drug interaction data, respectively.
In our experiments, we select the ? to be 85, and this designates 15% of the available data to be used for testing and the remaining for training. However, we do not use all of the remaining data for training. In order to test the potential of our method for di↵erent scenarios of data availability, we further sample from the 85% of the available data to cover 50%, 60%, 70%, 80% of the entire available data. For instance, for 50%, we only use a randomly selected subset of only 60% of the data (although more data is available, and hide the rest of the data from %0AC80; 814AB). To ensure the consistency in our experiments, as we grow the training dataset with more samples, we made sure that every anchor drug, library drug, and cell line pair in 50% is contained within 60%, every sample in 60% is contained within 70%, and every sample in 70% is contained within 80%.
Methods compared. We compare the prediction performance of %0AC80; 814AB with two state-of-the-art methods. These are non-fiber, well-known tensor completion algorithm [6,10] and Medians of Means (MoM) [21] approaches.
Other details. We implement all the algorithms in Python and perform the experiments on a Mac machine with 8-Core Intel Core i9 processors running at 2.3 GHz, and 16GBs RAM.
We start by evaluating how well %0AC80; 814AB predicts the missing values. We apply our %0AC80; 814AB algorithm, non-fiber and MoM algorithms separately to compute the missing values in drug-drug interactions dataset. Then, we calculate the normalized mean squared prediction error over all test samples with di↵erent percentages of data availability. Figure 2 presents the normalized mean squared loss values in the presence of di↵erent data availabilities. We observe that regardless of how much data is available, %0AC80; 814AB performs better than non-fiber [10] and MoM [21] algorithms. %0AC80; 814AB reaches to a smaller loss value than both methods in less than 400 seconds. Increasing the available data from 50% to 70%, we see that the normalized MSE loss value decreases the most for %0AC80; 814AB among the three methods. This implies that as partial fibers get denser (i.e., has more data points) our algorithm converges faster. We observe that the three compared algorithms benefit from the increasing data availabilities in di↵erent amounts. The converged loss value for %0AC80; 814AB, Non Fiber, and MoM [21] decreases by 24.2%, 15.80% and 1.86% respectively as we increase the available data percentage from 50% to 70%. We observe the smallest change in MoM. This implies that the newly added data do not have too much of an impact on the overall prediction error of MoM. We conjecture that this is because MoM selectively uses only subsets of the new data depending on the tensor entry to be predicted.
We hypothesize that the more data we have, the more simultaneous updates %0AC80; 814AB can provide. This implies that we expect that %0AC80; 814AB will converge faster with more data available. Figure 3 shows the convergence of our method when 50% and 80% of data is available. We observe that for both data availabilities, %0AC80; 814AB successfully achieves a small normalized mean squared loss value. This is because, as we provide more data, the number of indices that are updated simultaneously is expected to increase. This allows us to get a smaller loss value in less time. Our results support our hypothesis: Our method for 80% data availability is always below that for 50% in Figure 3. This means that we can reach the same loss value in less time by using more data. (a) Anchor Drugs -Partial Fiber (b) Library Drugs-Partial Fiber Fig. 4. The loss values for each anchor drug and each library drug after convergence and the sorted loss values with a focus on the top 20 drugs with the highest loss values.
Here, we present our results that focus on the drugs which have the most/least accurate predictions. We calculate the individual loss values for each of the anchor drugs and library drugs, both with partial fiber setting and non-fiber setting. We sort the drugs based on their loss values and share the name of the drugs that have the highest loss in each setting. We share the detailed results in Figure 4. We observe that more than half (⇠ 20) of the anchor drugs, and more than half of the library drugs (⇠ 40) reach a very small loss value, in both fiber and non-fiber settings. This implies that both of the methods are e ciently predicting the drug response for most of the drugs. Our results demonstrate that loss values are smaller for the %0AC80; 814AB for both anchor and library drugs. When we focus on 20 drugs which have the highest loss, for anchor drugs we see that 12 out of 20 drugs are common for both %0AC80; 814AB and non-fiber approaches, and if we narrow down and look at the top 10 drugs we see that 7 of them are same for both of the algorithms. The order might be di↵erent but for the same set of drugs we observe poor prediction performance, and we see that the loss is always less with %0AC80; 814AB. The common anchor drugs are as follows; Crizotinib, Trametinib, Oxaliplatin, SB505124, LGK974, AZD8055, Vinorelbine, AZD7762, 5-Fluorouracil, Navitoclax, Olaparib, Sapitinib. Each drug in this set targets di↵erent type of cancer. It is important to note the weaknesses before we can make a trustworthy data completion.
In this paper, we presented a novel algorithm called %0AC80; 814AB for predicting the drug dosage in drug-drug combination setting. Given partial information on the dosage of drug-drug combinations, %0AC80; 814AB e↵ectively predicts the IC50 values for the missing entries by specifically benefiting from the structure of the data without further assumptions. We showed that %0AC80; 814AB is a time e cient, flexible, scalable computational method which can provide accurate predictions.