A Kernelisation Approach for Multiple d-Hitting Set and Its Application in Optimal Multi-Drug Therapeutic Combinations

Therapies consisting of a combination of agents are an attractive proposition, especially in the context of diseases such as cancer, which can manifest with a variety of tumor types in a single case. However uncovering usable drug combinations is expensive both financially and temporally. By employing computational methods to identify candidate combinations with a greater likelihood of success we can avoid these problems, even when the amount of data is prohibitively large. Hitting Set is a combinatorial problem that has useful application across many fields, however as it is NP-complete it is traditionally considered hard to solve exactly. We introduce a more general version of the problem (α,β,d)-Hitting Set, which allows more precise control over how and what the hitting set targets. Employing the framework of Parameterized Complexity we show that despite being NP-complete, the (α,β,d)-Hitting Set problem is fixed-parameter tractable with a kernel of size O(αdkd) when we parameterize by the size k of the hitting set and the maximum number α of the minimum number of hits, and taking the maximum degree d of the target sets as a constant. We demonstrate the application of this problem to multiple drug selection for cancer therapy, showing the flexibility of the problem in tailoring such drug sets. The fixed-parameter tractability result indicates that for low values of the parameters the problem can be solved quickly using exact methods. We also demonstrate that the problem is indeed practical, with computation times on the order of 5 seconds, as compared to previous Hitting Set applications using the same dataset which exhibited times on the order of 1 day, even with relatively relaxed notions for what constitutes a low value for the parameters. Furthermore the existence of a kernelization for (α,β,d)-Hitting Set indicates that the problem is readily scalable to large datasets.


Introduction
Typically the selection of a drug therapy for a disease is limited to a single drug, however diseases such as cancer may present as a heterogeneous mix of subtypes of the general disease. In cases such as these multi-drug therapies may prove more effective than single drug therapies, and many trials have been conducted to this end [1][2][3]. Furthermore combinations of drugs may allow a more targeted approach for a selection of subtypes of a disease, while minimizing effects on unaffected cells. Unfortunately with the abundance of compounds available for the treatment of many conditions of interest, the time and expense in testing even all two drug combinations may be prohibitive. Therefore a smarter approach is needed. Vazquez [4] introduces the HITTING SET problem for this task in the context of oncological drug therapy. The HITTING SET problem is a combinatorial problem that proves extremely useful in modeling a large variety of problems in many domains including protein network discovery [5], metabolic network analysis [6], diagnostics [7][8][9], gene ontology [10] and gene expression analysis [11,12].

The Hitting Set Problem
HITTING SET is a combinatorial problem that models the problem of selecting a small group of elements to represent or cover a collection of sets. Such a group that covers every set in the collection is called a hitting set. Finding such a set without any constraint is simple, however if we required that the size of the hitting set be relatively small, the problem becomes computationally challenging (NP-complete in a formal sense). This difficulty in obtaining solutions with desirable qualities thus requires more thoughtful approaches.
We now give some technical details and formal definitions of the problems of interest.
HITTING SET is equivalent to the SET COVER problem [13], and when otherwise unrestricted, is equivalent to the RED/BLUE DOMINATING SET [14] problem and is related to the k-FEATURE SET [15] problem.
The decision version of the HITTING SET problem is defined as follows:

HITTING SET
Instance: A set S and a collection C(2 S and an integer k. Question: Is there a set S'(S with DS'Dƒk such that for every c [ C we have c\S'=w? The set S' is called a hitting set for C, or simply a hitting set. For an element s [ S' and an element C [ C if s [ c we say that s hits c. This problem is NP-complete even when the maximum size of each element of C is two (by equivalence with VERTEX COVER [13]) and W ½2-complete for parameter k; Cotta and Moscato [16] give a parameterized proof via k-FEATURE SET and Paz and Moran [17] give a proof which along with the equivalence of HITTING SET and SET COVER leads to the same result, though predates the parameterized complexity framework. However if we restrict the cardinality of the elements of C to d the problem, while remaining NP-complete, becomes fixed-parameter tractable where d is a constant and the parameter is k [18]. In this case the problem is known as the HITTING SET FOR SETS OF SIZE d or d-HITTING SET problem. We note that HITTING SET has several equivalent formulations, in particular we choose to use the bipartite graph representation where S and C form the two partite vertex sets of the graph and an edge sc corresponds to the element s [ S being an element of c [ C. This allows us to employ some simplifying graph theoretic terminology and techniques. We generalize this problem to include the case where we may want the elements of C to be hit more than once. In particular this includes the case where we ask if all the sets of C can be hit a times, but extends to the case where the elements of C can be hit up to a times. We encode this by the use of a hitting function g. Our problem then becomes the a-MULTIPLE d-HITTING SET (or (a,d)-HITTING SET): When g(c)~a for all c [ C, (a,d)-HITTING SET can be (1z ln d)approximated in time O(a : DCD : DSD) [19], but cannot be approximated with a factor of (1{e) ln n for any e [ (0,1) unless NP(DTIME(n log log n ) [20].

Results and Discussion
The Fixed-Parameter Tractability of (a,d)-Hitting Set As we prove in the Materials and Methods section, the (a,d)-HITTING SET problem is fixed-parameter tractable, and indeed a more general variant the (a,b,d)-HITTING SET problem is also fixed parameter tractable when we take the maximum degree d of the class vertices C as a constant and the size k of the hitting set and the maximum desired coverage a as a joint parameter. Though the problem is formally hard -which would normally give the intuition that an exact solution would be too expensive to compute -the fixed-parameter tractability indicates that it is likely that we can obtain an exact solution efficiently. Armed with this knowledge we proceed with the experiments of the following section, where we use the drug response data of the NCI60 anti-tumor drug screening program to determine a sets of drugs that hit cancerous cell lines multiple times. These drug sets are than mathematically supportable candidates for combination chemotherapies. Moreover we are able to tune the nature of the hitting sets via the numbers k, a and b, which allows us to control which cell lines are targetted (and which are specifically not) and how much each cell line is hit in the solution.

A Comparative Application
The NCI60 human tumor anti-cancer drug screen dataset [21] was established in the 1980s as an enabling tool for anti-cancer drug development. Included in this dataset is response data for over 40,000 drugs against the 60 cell lines of the dataset. Vazquez [4] highlights the utility of a hitting set approach in developing multi-drug therapies for heterogeneous malignancies; given the plethora of available compounds, testing multi-drug combinations exhaustively is prohibitive if not impossible. Applying hitting set to efficacy data measured on an individual basis for each compound allows us to determine possible drug combinations that would provide the best chance of efficacy against many cancer types. Using the GI50 response NCI60 dataset (available from the DTP website [22]) Vazquez uncovers a minimum hitting set with three compounds that cumulatively gives a good response with all cell lines in the dataset, where a response is considered good if it is more than two standard deviations above the mean of the ztransformed response data. Vazquez uses first a greedy highestdegree-first approach to give an estimate of the maximum size of a minimum hitting set, followed by either an exhaustive search or simulated annealing, depending on the size of the hitting set. Vazquez reports times for such approaches on the order of one day on a desktop computer.
We revisit Vasquez's experiment, using data reduction (though it is not necessary to employ the more complex rules given in the kernelization proof) with IBM ILOG CPLEX [23] as the kernel solver by framing the problem as a integer programming problem. We use the same threshold for the z-transformation to identify significant response levels. Using this approach we reduce the time to solve the instance to less than 5 seconds, where most of the time is spent loading and reducing the data, with CPLEX solving the integer programming instance in approximately 0:08 milliseconds. Furthermore this approach guarantees optimality in the size of the hitting set.
From here we employ more a more recent version of the NCI60 dataset (2009 as compared to Vazquez's 2006). At the time of writing, the latest NCI60 dataset includes 14 additional cell lines, however we remove these, as there is insufficient response data in the dataset, leading to inflated hitting set sizes. The latest data also includes a further 2281 compounds. We note that employing the new GI50 response data we are able to uncover 3 element hitting sets involving compounds not available in the earlier dataset (an example is given in Table 1 and Figure 1), in particular Everolimus (NSC 733504) a drug now used for the treatment of advanced renal cancer which is also giving positive results in phase II trials for metastatic melanoma [24,25]. However there have recently been some concerns over the provenance of some of the cell lines in the NCI60 dataset. In particular Lorenzi et al. [26] suggested that the MDA-N cell line, nominally a breast cancer cell line is in fact similar the M14 and MDA-MB-435 cell lines, and thus should be is in fact a melanoma cell line. Chambers [27] however suggests that although M14 and MDA-MB-435 are identical cell lines, they may not in fact be melanoma cell lines. We do not attempt to resolve this dispute, however with regard to this, and as a indication of the flexibility of the method we employ we consider both the case where MDA-N is a breast cancer cell line and the the case where MDA-N is a melanoma cell line.
Employing the (a,b,d)-HITTING SET model gives more flexibility in what kind of therapy we would like to pursue. For instance, by choosing g 1~2 for all vertices, we are able to find a hitting set that hits every cell line at least twice (see Table 2). However the size of this hitting set is 6, which is likely to be beyond the point where the trade off between anti-cancer efficacy and side effects is acceptable. Fortunately we can exploit (a,b,d)-HITTING SET more intelligently. For example we may wish to find a hitting set that specifically targets breast cancer cell lines -for which we set all breast cancer cell line vertices to have g 1~1 and all other cell lines to have g 2~0 . This gives a hitting set that hits only breast cancer cell lines, which may be useful in minimizing unwanted peripheral damage to non-breast cancer cells. This gives a hitting set with three elements. In the case where we considered MDA-N to be a breast cancer cell line (see Table 3 and Figure 2) this set includes the compound deoxypodophyllotoxin, which is known to induce apoptosis [28]. If we consider MDA-N as a melanoma cell line we obtain a different hitting set (see Table 4 and Figure 3). If we relax our requirements an allow other cell lines to be hit at most once we can obtain a hitting set that hits the breast cancer cell lines more (Table 5 and Figure 4). The results when we set g 1 to 2 for all breast cancer lines are given in Table 6 and Figure 5 (including MDA-N) and Table 7 and Figure 6 (excluding MDA-N). We note particularly that in the case where MDA-N is included, the optimal hitting set uncovered includes Docetaxel, a well known anti-cancer agent [29] for several cancer types including breast cancer. Interestingly Docetaxel is also currently included in several clinical trials examining its potential as part of a multi-drug therapy [30][31][32][33][34].
In another example, we may wish to target melanoma cell lines exclusively, and furthermore, we may wish to attack each cell line with at least two drugs at once. However in this case (where g 1~2 for melanoma cell lines and g 2~0 for all others) the minimal hitting set size is 6 (or 5 if MDA-N is included as a melanoma cell line - Table 8 and Figures 7 & 8). Considering that a therapeutic cocktail involving 6 compounds may have excessive side effects, we can relax the requirements, and allow g 2~1 for non-melanoma cell lines. In this case we find that the smallest hitting set is of size 3. By altering the focus when solving the kernel by fixing the hitting set size (k) at 3 and maximizing the total degree of the vertices in the hitting set, subject to the g 1 and g 2 constraints, we can obtain the minimal size hitting set that hits our targets as much as possible, within the bounds given by the constraints. This results in the hitting sets in Tables 9 & 10 and Figures 9 & 10. Of note is AZD6244, which is currently involved in 21 anti-cancer drug trials [35] and has been identified as a potent kinase inhibitor [36,37].

Conclusion
Given the size of modern datasets, and the expectation that they will only get larger, it is clear that we require efficient approaches to solving important computational biology problems. The first phase of any such approach is simply defining the problem at hand. Unfortunately once clearly stated, many such problems are NP-hard or worse. However this need not mean that we must resort to inexact or approximate approaches, which could be undesirable in a field such as drug selection. Parameterized Complexity provides a toolkit for dealing with nominally hard problems, and identifying cases where despite super-polynomial running times, we may still expect good performance.
The drug selection problem as examined here is one such problem. It is modeled well by the d-HITTING SET problem, which is fixed-parameter tractable when parameterized by the maximum size of the hitting set. Therefore we can expect that despite being NP-complete, it would be relatively quick to solve when these parameters are small. However we demonstrate that the much more flexible variant (a,b,d)-HITTING SET is also fixed-parameter tractable, with only the addition of a single parameter -the maximum of the minimum number of times any vertex should be hit. With (a,b,d)-HITTING SET we are able to better control the nature of the hitting set uncovered, and thus tailor any such hitting set to a useful set of constraints, such as limits on which cell lines are to be hit, the maximum any of these can be hit and of course the minimum number of times any cell line should be hit. Moreover we can solve this problem quickly, and guarantee optimality -without any notable restrictions on the parameters and constants. This allows the quick generation of possible drug combinations for testing, with guarantees of a certain baseline performance, eliminating the need to exhaustively test all possible combinations, which would be financially and temporally prohibitive.
In brief this paper provides a robust and flexible methodology for multiple drug selection, which can easily be applied to other domains that are modeled by the d-HITTING SET problem, with a sound theoretical background as to why and how the problem can be solved efficiently, despite its NP-completeness. Moreover the existence of a kernelization for (a,b,d)-HITTING SET indicates that even without using a specialized commercial solver such as CPLEX, the problem is readily scalable to large datasets. Given the speed at which we are able to solve instances with on the order of 40,000 vertices, we can expect that much larger datasets are also solvable in a reasonable time.
A future extension that may be of interest would be to somehow encode in the problem the notion that some hitting vertices are incompatible, e.g., two compound may have severe adverse interactions, and thus can never be used together as a therapy, regardless of their individual usefulness.

Dataset and Computational Method
The dataset primarily employed is the NCI60 DTP Human Tumor Cell Line Screen, available from [22]. We use the version released in October 2009, and downloaded in April 2010. The raw dataset is presented as a series of cell line and compound pairs, along with the GI50 response measurement (the method for producing the measurements is also detailed by [22]) for that pair plus concentration information and statistical information. Where there are multiple entries for the same compound-cell line pair, we select the entry resulting from the experiment using the highest concentration of the compound. We extract this data into a matrix cross indexed by the NSC number of the compound and the name of the cell line. Where an entry does not exist for a given compound-cell line pair, we enter ''NA'' for that entry in the matrix.
Once the data is in this matrix format we threshold the data according to the method used by Vazquez [4] whereby the raw data is subject to a z-transformation over a logarithmic scale and then any value above a certain threshold expressed in terms of the standard deviation to 1, and anything below, including ''NA'' values, to 0. In line with Vazquez we choose two standard deviations as our particular threshold for this paper, though this is adjustable.
We then construct a graph for the hitting set instance using the Java Universal Network/Graph Framework (JUNG) [38] with the SetHypergraph class, representing each compound with a vertex and each cell line with a (hyper)edge which carries a weight indicating the number of times that edge is to be hit. This graph is then reduced to remove vertices of zero degree, edges with no incident vertices (which are noted as technically this would indicate a no instance unless that edge does not require hitting) and vertices that are only adjacent to edges that require zero hits. This basic reduction alone typically reduces the number of vertices significantly, bringing the graph within a reasonable size for immediate processing. From a theoretical standpoint the constant d is of importance, for the graph constructed as stated, d~4741 (as we allow the natural value, rather than imposing an external limit). In practice a d value of this magnitude proves perfectly workable, and returning to the theoretical viewpoint indicates that the instance is in a sense already kernelized.
Once the graph is reduced, we construct an integer programming instance equivalent of the problem given the graph, and pass this instance to CPLEX [23] (version 11.200) and search for an optimal solution to one of two objective functions, given the constraints of the number of hits for each cell line (given by the g 1 value). The first objective function simply minimizes the size of the hitting set (k), for the second objective function we fix the size of the hitting set, and maximize the number of hits on vertices where no maximum number of hits has been set (the g 2 value). As part of this search CPLEX may apply some unspecified proprietary reduction process.
The figures were created using yEd Graph Editor [39]. The computer hardware employed is a Dell PowerEdge III Dual Xeon 5550 server with 32Gb of RAM, operating Red Hat Linux 64 bit EL 4 Server.    Table 3. Minimal hitting set targeting only breast cancer.   Parameterized Complexity. A parameterized (decision) problem is a formally defined computational problem consisting of three components; the input, a special part of the input called the parameter, and the question. Following Flum and Grohe's [40] definition we may assume that the parameter is derived from a polynomial time computable mapping from the input to the natural numbers. A parameterized problem P is fixed-parameter tractable if there is an algorithm A such that for every instance (x,k) where x is the input, k is the parameter and DxD~n, A correctly answers YES or NO in time bounded by f (k)p(n) where p is a polynomial and f is a computable function.

Theoretical Background and Kernelization Proof
A polynomial time kernelization (or just kernelization) is a polynomial time mapping that given an instance (x,k) of a parameterized problem produces a new instance (x',k') of the problem such that: 1. x is a YES-instance if and only if x' is a YES-instance, 2. k'ƒk and 3. Dx'Dƒg(k') for some computable function g.
It is easy to see that if a problem has kernelization, then it is fixed-parameter tractable. It is also easy to prove that if a problem is fixed-parameter tractable, then it has a kernelization [41].
Parameterized complexity has a fully developed theory for determining when a problem is unlikely to be fixed-parameter tractable, but as this is not necessary for this work, we refer the reader to the monographs of Flum and Grohe [40] and Downey and Fellows [42] for full discussion, and simply state that if a problem is W ½t-hard or W ½t-complete for any t [ N z , then the problem is not fixed-parameter tractable unless certain complexity theoretic assumptions are false, which seems unlikely.

The Fixed-Parameter Tractability of (a,d)-Hitting Set
Our kernelization for (a,d)-HITTING SET follows the basic format of Abu-Khzam's kernelization for d-HITTING SET [18].
Let (G,k) be an instance of (a,d)-HITTING SET which we assume to have been preprocessed for nonsense input such as vertices c [ C with d(c)wd or d(c)vg(c). Therefore we may assume that for all c [ C we have g(c)ƒd(c)ƒd and that for all vertices s [ S we have d(s) §0.
We first apply Reduction Rules 1 to 3 exhaustively, before applying Rules 4 and 5.: by 1, delete s from G and reduce k by 1. Finally, delete c from G.
Lemma 1 Reduction Rule 1 is sound.
Proof. If such a vertex c exists, then all its neighbors in S must be in the hitting set, and we can remove them from the graph after suitably noting the effect for the vertices of N(N(c)).
Note in particular that this rule effectively allows us to assume that m is at most d{1. This will be used implicitly in Reduction Proof. If two such vertices c and b exist, then any hitting set that hits c at least g(c) times will hit b at least g(b)vg(c) times.
Let B(S be a set of size d{1 vertices such that B is the pairwise intersection of the neighborhoods of a vertex set N(C. Let N i~f n [ NDg(n)~ig.
Reduction Rule 4: Let B(S and N(C be vertex sets as described. For each i [ ½1,a such that DN i Dwk add a vertex c to C with g(c)~i and edges such that N(c)~B and delete N i from G.
Lemma 4 Reduction Rule 4 is sound.
Proof. Let (G,k) be a YES-instance of (a,d)-HITTING SET. Then there is a set A(S with DADƒk that hits each element c of C at least g(c) times. Assume that there are sets B and N as described in the reduction rule and that for some i [ ½1,a we have that DN i Dwk. Let A Ni be the subset of A that hits N i . Assume further that A Ni 6(B, then for each n [ N i there is at least one other vertex in  Table 5. Minimal hitting set targeting breast cancer but allowing other cell lines to be hit.  Table 6. Minimal hitting set hitting breast cancer twice, and no others, with MDA-N.   Table 7. Minimal hitting set hitting breast cancer twice, and no others, without MDA-N.  In this case we allow non-breast cancer cell lines to be hit at most once. By relaxing the restriction on hitting non-breast cancer cell lines, we obtain a hitting set which hits more of the breast cancer cell lines repeatedly. The trade-off being that other cell lines are also affected, increasingly the likelihood that non-cancerous cells are also affected by the treatment, as the compounds are less specific to a particular genetic signature. Only cell lines with at least one adjacent compound are shown. doi:10.1371/journal.pone.0013055.g004 A Ni , but then DAD §DA Ni Dwk, which contradicts the assumption that (G,k) is a YES-instance. Therefore the set N i must be hit by B, so we may restrict our search to the intersection. .   Let W (C be a maximal set of pairwise weakly related vertices. Let B(S be a set of vertices, and denote by W B the set of vertices of W whose neighborhood is a superset of B. Delete W B,t from G.  Table 9. Minimal hitting set targeting melanoma, without MDA-N.

NSC Number
Compound Name  Lemma 7 Reduction Rule 5 is sound. Proof. We defer the proof of the bound on the size of W until the proof of Lemma 8.
Let (G,k) be a YES-instance of (a,d)-HITTING SET. Then there is a set A(S that hits S sufficiently. For sets of size d{1, Reduction Rule 4 proves the soundness of the first iteration of the outer loop.
For each other iteration, assume that the iteration for sets of size j holds, then let B be set of size j{1 where DW B,t Dwk d{jz1 for some t. If DA\BDvt then by the pigeon hole principle there is some vertex v [ A that is in at least k d{j neighborhoods of vertices in W B,t , but then B|fvg is a set that is the intersection of at least k d{j neighborhoods of vertices in some subset of W , contradicting the correctness of the previous iteration. Therefore the entire set of vertices hitting each W B,t vertex is contained within B if DW B,t Dwk d{j , so we may replace W B,t with a single vertex.
Note also that for each element of W there is at most 2 d sets B, so we may iterate through all sets in time O(ad2 d : DW D), so we can perform the replacements in polynomial time.
Lemma 8 If (G,k) is a YES-instance of (a,d)-HITTING SET, reduced under Reduction Rules 1 to 5, then DV (G)Dƒ(dz1)ak d .
Proof. If (G,k) is a YES-instance of (a,d)-HITTING SET, then there is a set A(S such that for every s [ S we have DN(s)\AD §g(s) with DADƒk. Claim 9 C~W . By construction, every vertex in C with degree at most d{1 is in W . Assume there is some c [ C with d(c)~d and c 6 [W , then there must be some vertex c' [ C such that DN(c)\N(c')Dwd{1, but then as the degree of any vertex in C is at most d, N(c)~N(c'), and Reduction Rule 3 would apply. Therefore there are no vertices from C not in W .
Claim 10 DW Dƒak d . As A hits each vertex of W at least once, by Reduction Rule 5 each element of A as a singleton is in the neighborhood of at most ak d{1 vertices from C. Therefore DW DƒDAD : ak d{1 ƒak d . The key difference with the case where we consider MDA-N to be a non-melanoma cell line is that in this case we obtain a hitting set that hits the melanoma cell lines slightly more. Only cell lines with at least one adjacent compound are shown. doi:10.1371/journal.pone.0013055.g010 Combining Claims 9 and 10 we have DCDƒak d . As each vertex of C has degree at most d, there are at most adk d vertices in S, and the bound follows.
Theorem 11 (a,d)-HITTING SET is fixed-parameter tractable with parameter k and has a kernel of size at most adk d .
We note that although d must be a constant to obtain a polynomial time kernelization, a may be alternatively given as an additional parameter, without change to the kernelization.
This kernelization may be extended to an even more general version of the problem, where we not only specify lower bounds for the number of hits, but also upper bounds: (a,b,d)-HITTING SET Instance: A bipartite graph G~(S ] C,E) where for all c [ C we have d(c)ƒd, two hitting functions g 1 : C?½0,a and g 2 : C?½0,b and an integer k. Question: Is there a set S'(S with DS'Dƒk such that for every c [ C we have g 2 (c) §DN(c)\S'D §g 1 (c)?
Corollary 12 (a,b,d)-HITTING SET is fixed-parameter tractable with parameter k and has a kernel of size at most adk d .