BACKGROUND
in the last decade, a number of drugs targeting specific biologically relevant kinases have been developed that are becoming common in cancer research as a basis for personalized therapy.
the idea of treating cancer through inhibition of a specific tyrosine kinase was proven by the discovery that patients with chronic myeloid leukemia can be successfully treated by inhibiting the tyrosine kinase bcr-abl with the kinase inhibitor imatinib mesylate  <cit> .
however, the success rate of any one specific targeted drug for other forms of cancer, such as sarcoma, is limited as the tumors exhibit a wide variety of signaling pathways and are not uniformly dependent on the activity of a specific kinase  <cit> .
the numerous aberrations in molecular pathways that can produce cancer is one cause to necessitate the use of drug combinations for treatment of individual cancers.
combination therapy design requires a framework for inference of the individual tumor pathways, prediction of tumor sensitivity to targeted drug and algorithms for selection of the drug combinations under different constraints.
the current state of the art in predicting sensitivity to drugs is primarily based on assays measuring gene expression, protein abundance and genetic mutations of tumors; these methods often have low accuracy due to the breadth of available expression data coupled with the absence of information on the functional importance of many genetic mutations.
a commonly used method for predicting the success of targeted drugs for a tumor sample is based on the genetic aberrations in the tumor .
however, the accuracy of prediction of drug sensitivity based on mutation knowledge is limited in many forms of tumors as some of the mutations  may not be functionally important or tumors can develop without the known genetic mutations.
statistical tests have been used in  <cit>  to show that genetic mutations can be predictive of the drug sensitivity in non-small cell lung cancers but the classification rates of these predictors based on individual mutations for the aberrant samples are still low.
for specific diseases, some mutations have been able to predict the patients that will not respond to particular therapies: for instance  <cit>  reports a success rate of 87% in predicting non-responders to anti-egfr monoclonal antibodies using the mutational status of kras, braf, pik3ca and pten.
the prediction of tumor sensitivity to drugs has also been approached as a classification problem using gene expression profiles.
in  <cit> , gene expression profiles are used to predict the binarized efficacy of a drug over a cell line with the accuracy of the designed classifiers ranging from 64% to 92%.
in  <cit> , a co-expression extrapolation  approach is used to predict the binarized drug sensitivity in data points outside the training set with an accuracy of around 75%.
in  <cit> , a random forest based ensemble approach was used for prediction of drug sensitivity and achieved an r <dig> value of  <dig>  between the predicted ic50s and experimental ic50s.
supervised machine learning approaches using genomic signatures achieved a specificity and sensitivity of higher than 70% for prediction of drug response in  <cit> .
tumor sensitivity prediction has also been considered as  a drug-induced topology alteration  <cit>  using phospho-proteomic signals and prior biological knowledge of a generic pathway and  a molecular tumor profile based prediction  <cit> .
most interestingly, in the recent cancer cell line encyclopedia  study  <cit> , the authors characterize a large set of cell lines  with numerous associated data measurement sets: gene and protein expression profiles, mutation profiles, methylation data along with the response of around  <dig> of these cells lines across  <dig> anti-cancer drugs.
one of the goals of the study was to enable predictive modeling of cancer drug sensitivity.
for generating predictive models, the authors considered regression based analysis across input features of gene and protein expression profiles, mutation profiles and methylation data.
the performance  of the predictive models using  <dig> fold cross validation ranged between  <dig>  to  <dig> .
in particular, the correlation coefficient for prediction of sensitivity using genomic signatures for the drug erlotinib across > <dig> cell lines was < <dig> .
erlotinib is a commonly used tryosine kinase inhibitor selected primarily as an egfr inhibitor.
however, studies have shown  <cit>  that these targeted drugs often have numerous side targets that can play significant roles in the effectiveness of the inhibitor drugs.
the target inhibition profiles of drugs and sensitivity of trainings set of drugs can provide significant information for enhanced prediction of anti-cancer drug sensitivity as we have recently shown  <cit> .
by incorporating the drug-target interaction data and sensitivities of training drugs with genomic signatures, we were able to achieve a correlation coefficient of  <dig>   for prediction of erlotinib sensitivity using  <dig> fold cross validation.
the result illustrates the fundamental concept of the importance of drug-target interaction and functional data under which we develop the sensitivity prediction method presented in this paper.
by developing a framework around the functional and target information extracted from the primary tumor drug screen performed by our collaborators, we seek to develop a cohesive approach to sensitivity prediction and combination therapy design.
this necessitates the generation of the tumor pathway structure for individual patients to decide on the target inhibitors for therapy based on the personalized patient pathways.
we envision that the overall schematic of the design of personalized pathways and personalized therapy will be similar to the workflow shown in figure  <dig>
the explanations of the various steps in the design process are as follows:
a patient is diagnosed with cancer and a primary culture of the tumor is established.
cell viability after exposure to targeted drugs is measured through a drug screen.
use of this functional data rather than mutation or protein biomarkers provides a unique advantage.
a target inhibition map  is generated based on the ic50′s and the known targets of the drugs in the screen.
tim denotes a predictive model that provides the sensitivity for all possible target inhibitions.
specifically, a tim is composed of a set t={t <dig> t <dig> ⋯,tn} consisting of binary variables, each denoting inhibition of a target, and a function f relating the target inhibitions to the steady state sensitivity yt, i.e.
yt=f.
the inhibition vector  corresponding to a drug is known as the drug target inhibition profile .
a detailed example of tim is provided in additional file  <dig>  the coarse structure of the tim can also be represented by an abstract pathway which will be termed tim circuit.
the construction of the tim circuit is explained in the methods section.
further data is collected using sirna screens  <cit> , rna sequencing and protein phosphoarrays to reduce model parameter uncertainties.
based on the knowledge of the tim and tim-directed protein expression measurements, the dynamic model is created.
combination therapy is designed utilizing the personalized tim and the dynamic model .
various constraints such as avoiding resistance to drugs or minimizing toxicity can be applied to design the combination therapy.
a mouse xenograft model  <cit>  can be used to study development of resistance simultaneously.
the generated drug combinations are validated in vitro on the primary culture.
if needed, the circuit is revised or the drug combination with best response in vivo  and in vitro  is then provided to the patient.
the primary contributions of this paper are:  methods for extraction of numerically relevant drug targets from single-run drug screens,  design of the personalized tim circuit based on drug perturbation data,  algorithms for sensitivity prediction of a new drug or drug cocktail,  validation over canine osteosarcoma primary tumors and  pathway flow inference using sequential protein expression measurements.
the scope of the present article is concentrated around steps b, c and d of figure  <dig>
the perturbation data required for our proposed method originates from a drug screen consisting of  <dig> small molecule inhibitors with quantified kinase interaction behaviors.
this drug screen, denoted drug screen version  <dig> , consists of two sets of data:  the first set is the experimentally generated drug sensitivities provided as 50% inhibitory concentration  values.
the ic <dig> values denote the amount of a drug required to reduce the population of cancerous cells in vitro by half.
the sensitivity values are expected to change during each new cell line/ tumor culture experiment.
the generation of the sensitivities in step c can be done within  <dig> hours of initial biopsy using drug sensitivity assays which is a period of limited cell divisions for most primary cultures.
thus, the estimated personalized maps may be closer to real-time circuits in cancer cells - akin to the signaling found in an untreated patient within a day or two after biopsy, and not the evolving consensus pattern of signaling for growing and dividing tumor cells as subpopulations emerge with increased fitness in vitro.
in addition, the drug screen contains experimentally derived half-maximal concentration  values for the interaction of each drug and each kinase target.
the ec <dig> value is directly related to the notion of inhibition of a kinase target; in particular, the ec <dig> values correspond to the amount of a compound needed to deactivate via phosphorylation 50% of the population of the associated target.
hence, for a drug compound, a target with a lower ec <dig> is the one that will be heavily inhibited at low drug concentration levels.
thus, low ec <dig> targets are often considered to be the primary targets of a drug.
the remaining targets are considered to be the side targets of a drug, and are often ignored.
the utility of this ec <dig> data is its consistency throughout experiments; the ec <dig> values as curated from literature searches  <cit>  are fixed, regardless of change of tumor type or patient of origin.
this provides a great amount of prior information for analysis of the drug screen results, and its usage is supported from the experiments performed in  <cit> .
the overall goal of the methods presented in this paper is to create an input-output mathematical framework for the analysis of and inference on the functional data generated by the drug screens for the purpose of anti-cancer drug sensitivity prediction and inference of personalized tumor survival pathway.
the personalized tumor survival pathway refers to the visual circuit diagram generated from the inferred target inhibition map as explained in the methods section.
note that the circuit corresponding to a tim is only a coarse representation of the tim for visual understanding of the most probable target combinations whose inhibition can reduce the tumor survival.
since the experiments were conducted on in-vitro cell cultures with the output being cell viability measured in terms of ic <dig>  the survival here refers to tumor cell culture survival and not the overall survival of the patient.
RESULTS
tim generation for canine osteosarcoma tumor cultures and cross-validation estimates of prediction accuracy
the sensitivity prediction and circuit analysis performed on actual biological data are validations of the proposed methodology to be described in the methods section.
the experimental data on four tumor cultures and  <dig> targeted drug screen panel were generated in the keller laboratory at ohsu.
the cell lines applied to the drug screen were four canine osteosarcoma cell lines cultured from four distinct canines, denoted bailey, charley, sy, and cora.
the tumor cultures were collected by dr. bernard seguin of oregon state university from canines that are part of an ongoing clinical trial for osteosarcoma .
the tumor samples were collected from client-owned animals that have developed the disease naturally.
all procedures performed on these animals with regards to tumor collection were strictly for treatment purposes and nothing was done different because of the drug perturbation study.
all procedures were performed according to standard of care regardless of whether an animal had its tumor sampled.
for the generation of the experimental data, the canine osteosarcoma primary cell cultures were plated in  <dig> well plates at a seeding density of  <dig> cells per well over graded concentrations of  <dig> small-molecule kinase inhibitors.
each inhibitor was plated individually at four concentrations predicted to bracket the ic <dig> for that drug.
cells were cultured in rpmi  <dig> supplemented with 2mm glutamine, 2mm sodium pyruvate, 2mm hepes, 1% penicillin streptomycin, and 10% fetal bovine serum for  <dig> hours.
at the end of the  <dig> hour incubation, cell viability was assessed using the mts assay.
all values were normalized to the mean of seven wells on each plate containing no drug.
the ic <dig> for each drug was then determined by identification of the two concentrations bracketing 50% cell viability and application of the following formula: +da where cell viability value above 50% = a  and cell viability value below 50% = b .
the experimentally generated ic <dig> values are included as additional file  <dig>  the experimentally generated sensitivities  of the  <dig> drugs are then scaled to values between  <dig> and  <dig>
among the  <dig> drugs on the drug screen,  <dig> drugs have known target inhibition profiles; of these  <dig> drugs,  <dig> provide information only on the target mtor  and analysis of these drugs are trivial.
thus, the remaining  <dig> drugs are used to generate the tims.
these target profiles were extracted from several literature sources  based on experimental quantitative dissociation constants  which are treated as ec <dig> values  for each drug across kinase target assays with more than  <dig> targets.
the target profiles of the drugs are shown in additional file  <dig>  figures  <dig> and  <dig> represent the equivalent tim circuits generated from experimental data for bailey and sy respectively.
the tim circuits for charley and cora are included in additional file  <dig>
to emphasize the biological relevance provided by the tim framework employed in the analysis of the biological data, we present a more in-depth analysis of the tim circuit devised for the canine patient bailey .
the vast majority of human osteosarcomas contain genetic or post-translational abnormalities in one or both of the tumor suppressors p <dig>  <cit>  and prb  <cit> .
the first target identified in this circuit is pkc alpha .
pkc alpha modifies cdkn1a , which is the primary mediator of p <dig> tumor suppressor activity  <cit> .
psmb <dig> represents the proteasome .
previous studies  <cit>  and early preclinical data from the keller laboratory confirms in vitro sensitivity of many osteosarcomas to proteasome inhibitors and this sensitivity is hypothesized to be due to the integral role of the proteasome in p <dig> regulation  <cit> .
interestingly, cdk <dig> is also prominent in this circuit, which is a primary inhibitor of the tumor suppressor prb, which is also frequently abnormal in spontaneous human osteosarcoma  <cit> .
cdk <dig> is an important modifier of both p <dig> and prb and is also represented in this circuit  <cit> .
the importance of pi3k pathway in osteosarcoma has also been recently reported using high throughput genotyping  <cit> .
our tim circuit includes akt <dig> which is downstream of pi3k  <cit> .
also, ednra selected in the circuit has been known to interact with pkc and activate erk signaling  <cit> .
if the circuit models shown in figures  <dig> and  <dig> are used to predict sensitivities for comparison with experimentally generated data, we will get optimistic results as the models are trained using the entirety of the available data.
thus, we utilize leave one out  and 10-fold cross validation  approaches to test the validity of the tim framework that we present in this paper.
for the loo approach, a single drug among the  <dig> drugs with known inhibition profiles is removed from the dataset and a tim is built, using the sffs suboptimal search algorithm, from the remaining drugs.
the resulting tim is then used to predict the sensitivity of the withheld drug.
the predicted sensitivity value is then compared to it’s experimental value; the loo error for each drug is the absolute value of the experimental sensitivity y minus the predicted sensitivity , i.e.
|y−|.
the closer the predicted value is to the experimentally generated sensitivity, the lower the error for the withheld drug.
tables  <dig>   <dig>   <dig> and  <dig> provides the complete loo error tables and the average loo error  for each primary culture.
the average loo error over the  <dig> cell cultures is  <dig>  or  <dig> %.
for the 10-fold cross validation error estimate, we divided the available drugs into  <dig> random sets of similar size and the testing is done on each fold while being trained on the remaining  <dig> folds.
this is repeated  <dig> times and average error calculated on the testing samples.
we again repeated this experiment  <dig> times and the average of those mean absolute errors for the primary cell cultures are shown in table  <dig>  the detailed results of the 10-fold cross validation error analysis are included in additional file  <dig>  we note that both 10-fold cv and loo estimates for all the cultures have errors less than 9%, which is extremely low, especially considering the still experimental nature of the drug screening process performed in the keller laboratory and the available response of only  <dig> drugs with known target inhibition profile.
pred.
sens.
= predicted sensitivity, exp.
sens.
= experimental sensitivity, error = |pred.
sens.
−exp.
sens.|, avg err = average of the  <dig> leave one out errors.
pred.
sens.
= predicted sensitivity, exp.
sens.
= experimental sensitivity, error = |pred.
sens.
−exp.
sens.|, avg err = average of the  <dig> leave one out errors.
pred.
sens.
= predicted sensitivity, exp.
sens.
= experimental sensitivity, error = |pred.
sens.
−exp.
sens.|, avg err = average of the  <dig> leave one out errors.
pred.
sens.
= predicted sensitivity, exp.
sens.
= experimental sensitivity, error = |pred.
sens.
− exp.
sens.|, avg err = average of the  <dig> leave one out errors.
the  <dig> fold- cross valdiation error estimates were calculated for  <dig> runs and the average mean absolute error  across all the drugs is shown.
to provide a measure of the overlap between drugs, we considered a similarity measure Λ based on the ec <dig> of the drugs d <dig> and d <dig>  let the ec50′s of the drugs d <dig> and d <dig> be given by the n-length vectors e <dig> and e <dig> where n denotes the number of drug targets.
the entries for the targets that are not inhibited by the drugs  are set to  <dig>  let the vectors v <dig> and v <dig> represent the binarized targets of the drugs i.e.
it has a value of  <dig> if the target is inhibited by the drug and a value of zero if the target is not inhibited by the drug.
then, we define the similarity measure Λ as:
Λ=∑i=1nmin,e2)∗v1∗v2∑i=1nmax,e2)
note that Λ= <dig> and similarity between drugs with no overlapping targets is zero.
if two drugs have 50% targets overlapping with same ec50′s, then the similarity measure is  <dig> .
the similarities between the drugs are shown in additional file  <dig>  note that except two drugs rapamycin and temsirolimus that have a similarity measure of  <dig> , all other drugs have significantly lower similarities with each other.
the maximum similarity between two different drugs  is  <dig> .
this shows that any two drugs in the drug screen are not significantly overlapping and the prediction algorithm is still able to predict the response.
the low error rate illustrates the accuracy and effectiveness of this novel method of modeling and sensitivity prediction.
furthermore, these error rates are significantly lower than those of any other sensitivity prediction methodology we have found.
consistent with the analysis in  <cit> , the sensitivity prediction rates improve dramatically when incorporating more information about drug-protein interaction.
to more effectively compare the results generated via the tim framework with the results in  <cit> , we also present the correlation coefficients between the predicted and experimental drug sensitivity values in table  <dig>  the correlation coefficients for predicted and experimentally generated sensitivities for  <dig> drugs and more than  <dig> cell lines ranges from  <dig>  to  <dig>  when genomic characterizations are used to predict the drug sensitivities in the ccle study  <cit> .
in comparison, our approach based on sensitivity data on training set of drugs and drug-protein interaction information produced correlation coefficients > <dig>   for both leave one out  and 10-fold cross validation  approaches for error estimation.
ρ loo = average pearson correlation coefficient for test samples selected based on leave one out approach, ρ 10-fold cv = average pearson correlation coefficient for test samples selected based on  <dig> fold cross validation approach.
it should be noted that the sensitivity prediction is performed in a continuous manner, not discretely, and thus effective dosage levels can be inferred from the predictions made from the tim.
this shows that the tim framework is capable of predicting the sensitivity to anti-cancer targeted drugs outside the training set, and as such is viable as a basis for a solution to the complicated problem of sensitivity prediction.
in addition, we tested the tim framework using synthetic data generated from a subsection of a human cancer pathway taken from the kegg database  <cit> .
here, the objective is to show that the proposed tim method generates models that highly represent the underlying biological network which was sampled via synthetic drug perturbation data.
this experiment replicates in synthesis the actual biological experiments performed at the keller laboratory at ohsu.
to utilize the tim algorithm, a panel of  <dig> targeted drugs pulled from a library of  <dig> is used as a training panel to sample the randomly generated network.
additionally, a panel of  <dig> drugs is drawn from the library to serve as a test panel.
the training panel and the testing panel have no drugs in common.
each of the  <dig> training drugs is applied to the network, and the sensitivity for each drug is recorded.
the generated tim is then sampled using the test panel which determines the predicted sensitivities of the test panel.
the synthetic experiments were performed for  <dig> randomly generated cancer subnetworks for each of n= <dig> ⋯, <dig> active targets in the network.
the active targets are those which, when inhibited, may have some effect on the cancer downstream.
to more accurately mimic the boolean nature of the biological networks, a drug which does not satisfy any of the boolean network equations will have sensitivity  <dig>  a drug which satisfies at least one network equation will have sensitivity  <dig>  the inhibition profile of the test drugs is used to predict the sensitivity  of the new drug.
the average number of correctly predicted drugs for each n is reported in table  <dig>  this synthetic modeling approach generally produces respectable levels of accuracy, with accuracies ranging from 89% to 99%.
<dig> drugs for training mimics the drug screen setup used by our collaborators and testing  <dig> drugs for predicted sensitivity approximates a secondary drug screen to pinpoint optimal therapies.
the performance of the synthetic data shows fairly high reliability of the predictions made by the tim approach.
correct prediction = average number of correct prediction for  <dig> testing drugs for  <dig> different synthetic pathways.
we have also tested our algorithm on another set of randomly generated synthetic pathways.
the detailed results of the experiment are included in additional file  <dig>  a large number of testing samples were used for each pathway prediction and the results indicate an average error of less than 10% for multiple scenarios.
in comparison, the average error with random predictions was 44%.
the average correlation coefficient of the prediction to actual sensitivity for the  <dig> sets of experiments  was  <dig> .
the average correlation coefficient with random predictions was  <dig>  we also report the standard deviation of the errors and for a representative example, the  <dig> percentile of the error was - <dig>  and  <dig> percentile  <dig> , thus the 80% prediction interval for prediction μ was .
the results of the synthetic experiments on different randomly generated pathways shows that the approach presented in the paper is able to utilize a small set of training drugs from all possible drugs to generate a high accuracy predictive model.
methods
in this section, we provide an overview of the model design and inference from drug perturbation data for personalized therapy.
mathematical formulation
let us consider that we have drug ic <dig> data for a new primary tumor after application of m drugs in a controlled drug screen.
let the known multi-target inhibiting sets for these drugs be denoted by s <dig>  s <dig> ..., sm obtained from drug inhibition studies  <cit> .
the elements of set si are ei=  for i= <dig> ,⋯,m, where ei,j are real-valued elements describing the interaction of si with k= , the set of all kinase targets included in the drug screen.
the ei,j’s refer to the ec <dig> values discussed previously.
it should be noted that for all si, ei,j will most often be blank or an extremely high number denoting no interaction.
the initial problem we wish to solve is to identify the minimal subset of k, the set of all tyrosine kinase targets inhibited by the m drugs in the drug panel, which explains numerically the various responses of the m drugs.
denote this minimal subset of k as t. the rationale behind minimization of t is twofold.
first, as with any classification or prediction problem, a primary goal is avoidance of overfitting.
secondly, by minimizing the cardinality of the target set required to explain the drug sensitivities found in the exploratory drug screen, the targets included have supportable numerical relevance increasing the likelihood of biological relevance.
additional targets may increase the cohesiveness of the biological story of the tumor, but will not have numerical evidence as support.
this set t will be the basis of our predictive model approach to sensitivity prediction.
before formulation of the problem for elucidating t, let us consider the nature of our desired approach to sensitivity prediction.
from the functional data gained from the drug screen, we wish to generate a personalized tumor survival pathway model instead of a linear function approximator with minimal error.
we are working under the fundamental assumption that the tumor survival pathway is nonlinear in its behavior; this assumption is reasonable given the difficulty in treating multiple forms of cancer.
one frequent theory in personalized therapy is that effective treatment results from applying treatment across multiple important biological pathways.
these pathways generally consist of sequentially activated gene and protein nodes acting as a feedback network.
treatment of individual pathways may not be sufficient for majority of diseases, so multiple independent parallel pathways must be targeted to create an effective treatment.
we believe that one possible approach to the analysis of multiple pathway treatment is to begin with an underlying framework based on the boolean interactions of the multiple targets in the pathway architecture.
the approach is based on developing families of boolean equations that describe the multiple treatment combinations capable of acting as an effective intervention strategy.
for the initial step of developing the underlying boolean functions, an initial binarization of the data set must be performed.
however, the resulting model lends itself to numerous continuous approaches to sensitivity prediction which we will explore further in the paper.
binarization of drug targets and conversion of ic50′s to sensitivities
in this subsection, we present algorithms for generation of binarized drug targets  and continuous sensitivity score of each drug .
the inputs for the algorithms in this subsection are the ec50′s of the drug targets and the ic50′s of the drugs when applied to a tumor culture.
in order to perform the binarization, we must consider the nature of the data we are given.
in particular, we are provided with an ic <dig> for each drug , and an ec <dig> value for each kinase target inhibited by the drug.
under the assumption that the primary mechanism of tumor eradication is, in fact, the protein kinase inhibition enacted by these targeted drugs, a natural consequence would be the existence of a relationship between the ic <dig> and ec <dig> values.
this relationship is explained as such: suppose for a drug si the ic <dig> value of si and the ec <dig> of kinase target kj,  are of similar value, then it can be reasonably assumed that kinase target kj is possibly a primary mechanism in the effectiveness of the drug.
in other words, if 50% inhibition of a kinase target directly correlates with 50% of the tumor cells losing viability, then inhibition of the kinase target is most likely one of the causes of cell death.
hence, the target that matches the drug ic <dig> is binarized as a target hit for the drug.
the above assumption of direct correlation for all successful drugs is obviously an extremely restrictive assumption and will be unable to produce high accuracy predictions.
thus, the binarization scheme has to be modified to incorporate the following three factors:
first: noises in varying magnitude will be present in the drug screen data generated by our collaborators.
the noise is unavoidable, and as such, needs to be accounted for.
in addition, despite the high accuracy of the drug-protein interaction data procured from literature, we should still account for possible errors in the ec <dig> values for the numerous drugs.
second: the restrictive assumption considers that effective drugs operate on single points of failure within the patient’s signaling pathway.
in reality, high sensitivity to a drug is often attributed to a family of related kinases  or several independent kinases working synergistically over one or multiple pathways to induce tumor death.
this cooperative multivariate behavior needs to be taken into account while binarizing a drug to its multiple possible targets.
third: despite the high level of currently available knowledge on the biological effects of numerous targeted drugs, there remains the possibility of a drug having high sensitivity while having no known mechanisms explaining its sensitivity.
therefore, we must consider the situation where there are latent mechanisms not considered within the dataset that are proving to be effective in some combination of treatment.
this point does not necessarily eliminate the possibility of kinase mechanisms being an important factor.
we address all three concerns as follows:  by considering the log scaled ec <dig> values for each target and the log scaled ic <dig> value for each drug, we convert the multiplicative noise to additive noise.
in addition, we employ scalable bounds around the ic50′s to determine binarization values of the numerous kinase targets for each drug.
the bounds can be scaled to allow targets that may have ec50′s higher than the ic <dig> to be considered as a possible treatment mechanism.
we extend the bounds to low ec <dig> levels, and often down to  <dig>  to incorporate the possibility of target collaboration at various different ec <dig> levels.
while a high ic <dig> indicates the likelihood of drug side targets as therapeutic mechanisms, it does not preclude the possibility of a joint relationship between a high ec <dig> target and a low ec <dig> target.
hence, to incorporate the numerous possible effective combinations implied by the ic <dig> of an effective drug, the binarization range of targets for a drug is the range α · log≤ log≤β · log where 0≤α≪β.
for reliability and validity of the target set that we aim to construct, it is important to keep β in a reasonable range, i.e.
β should be a smaller constant such as  <dig> or  <dig>  for the situation where the above bounds do not result in at least one binarized target, the immediate option is to eliminate the drug from the data set before target selection.
this prevents incomplete information from affecting the desired target set.
as information concerning the drug screen agents gradually becomes complete with respect to other forms of data, such as gene interaction data, additional mechanisms for unexplained targets can be explored and incorporated more readily into the predictive model.
with binarization of the data set as explained, we now present the minimization problem that produces a numerically relevant set of targets, t.
consider the target set t= , where ti∈{ <dig> }.
here,  <dig> denotes inclusion in the target set t and  <dig> denotes exclusion.
for any target set t <dig>  one can find the representation under t <dig> of each drug si,i∈ <dig> ⋯,m as  = .
as the t <dig> will be the basis of the new representation for each drug, this will result in n <dig> columns which will be  <dig> for all si, where n <dig> is the number of ti= <dig>  i.e.
the number of targets not included in t <dig>  the resulting representation of each drug in t <dig> is then an n−n <dig> vector of ec <dig> values.
while the representation of each drug will change as the target set t changes, the ic <dig> values for each of the m drugs remains the same.
these experimental sensitivity values will be used to test the numerous different target sets to quantify the strength of the model for any target set.
to simplify scoring of the target set, we first convert the ic <dig> for each drug si to a continuous-valued sensitivity score yi∈  <cit>
yi= <dig> ific <dig> i<cmaxic∗/log,ifcmax≤x≤maxdosei
where maxdosei is the maximum dose of drug si given, cmaxi is the maximum achievable clinical dose of drug si, and c=1−log/log so that the scoring function is continuous.
maxdose is used to prevent inferences being made on data that is not available.
while it would be possible to attempt interpolation to infer an ic <dig> from the multiple available data points, such inference cannot be fully quantified.
hence, drugs which fail to achieve an ic <dig> within the allotted dosage are given the score of  <dig>  which means ineffective.
the cmax value is used to apply a variable score to the numerous drugs based on the inherent toxicity of the drug.
this will also prevent bias towards drugs with low ic50’s; some drugs may achieve efficacy at higher levels solely based on the drug ec <dig> values.
construction of the relevant target set
in this subsection, we present approaches for selection of a smaller relevant set of targets t from the set of all possible targets k. the inputs for the algorithms in this subsection are the binarized drug targets and continuous sensitivity score .
with the scaled sensitivities, we can develop a fitness function to evaluate the model strength for an arbitrary set of targets.
as has been established, for any set of targets t <dig>  drug si has a unique representation .
this representation can be used to separate the drugs into different bins based on the targets it inhibits under t <dig>  within each of these bins will be several drugs with identical target profiles but different scaled scores.
let the set of scores in each bin be denoted y for sj in an arbitrary bin, and we will assign to each bin the mean sensitivity score of the bin, e).
denote this value p. within each bin, we want to minimize the variation between the predicted sensitivity for the target combination, p, and the experimental sensitivities, y. this notion is equivalent to minimizing the inconsistencies of the experimental sensitivity values with respect to the predicted sensitivity values for all known target combinations for any set of targets, which in turn suggests the selected target set effectively explains the mechanisms by which the effective drugs are able to kill cancerous cells.
numerically, we can calculate the inter-bin sensitivity error using the following equation:
∑bins∑j∈bin|p−y|
this analysis has one notable flaw: if we attempt to only separate the various drugs into bins based on inter-bin sensitivity error, we can create an over-fitted solution by breaking each drug into an individual bin.
we take two steps to avoid this.
first, we attempt to minimize the number of targets during construction of t <dig>  second, we incorporate an inconsistency term to account for target behavior that we consider to be biologically inaccurate.
to expand on the above point, we consider there are two complementary rules by which kinase targets behave.
research has shown that the bulk of viable kinase targets behave as tumor promoters, proteins whose presence and lack of inhibition is related to the continued survival and growth of a cancerous tumor.
these targets essentially have a positive correlation with cancer progression.
this is in opposition to tumor suppressors, proteins that have been shown to have a negative correlation with the development of cancer.
to capture the behavior of oncogenes, we partially formulate our problem on two rules  <cit> :
rule 1: if  is the inhibiting set of targets for drug i and the drug is successful in inhibiting the circuit, then any set b such that si⊂b will also be successful in inhibiting the circuit.
rule 2: if  is the inhibiting set of targets for drug i and the drug is unsuccessful in inhibiting the circuit, then any set b such that b⊂si will also be unsuccessful in inhibiting the circuit.
rule  <dig> essentially says that if inhibiting a number of target proteins has blocked signaling pathways, then inhibiting more target proteins will not open any path that has already been blocked.
rule  <dig> captures the fact that if a set of target protein inhibitors is unsuccessful in blocking the paths of a circuit, then any reduced number of target protein inhibitors among the inhibiting proteins cannot block all the paths.
the above rules assume that the kinases in focus are oncogenes, genes that promote cancer growth and whose inhibition can prevent tumor development.
the majority of kinases in the drug screen panel behave as oncogenes, and as such, our approach utilizes the above rules.
target sets resulting in combination scores that do not follow the rule-based behavior incur an inconsistency penalty.
this penalty is calculated as follows:
∑drugs∑binsχ|p−y|
where χ is the indicator function which is  <dig> when the experimental drug score is inconsistent with the predicted subset/superset bin score.
we now present the complete target set score, and as such, the equation that we wish to solve:
mint∑bins∑j∈bin|p−y|+∑drugs∑binsχ|p−y|
which reduces to the minimization problem we wish to solve:
mint∑drugs∑bins−y|+χ|p−y|)
for brevity, we will denote the scoring function of a target set with respect to the binarized ec <dig> values s and the scaled sensitivity scores y, Γ. as the s and y sets will be fixed when target set generation begins, we reduce this notation further to Γ.
note that t⊆k where k denotes the set of all possible targets.
2|k| is the total number of possibilities for t which is extremely huge and thus prohibits exhaustive search.
thus the inherently nonlinear and computational intensive target set selection optimization will be approached through suboptimal search methodologies.
a number of methods can be applied in this scenario and we have employed sequential floating forward search   <cit>  to build the target sets.
we selected sffs as it generally has fast convergence rates while simultaneously allowing for a large search space within a short runtime.
additionally, it naturally incorporates the desired target set minimization aim as sffs will not add features that provide no benefit.
we present the sffs algorithm for construction of the minimizing target set in algorithm  <dig>
algorithm  <dig> sffs algorithm for generating relevant target sets
complexity of target set generation
the algorithm to generate the error score given a target set t is of order o, quadratic with respect to the number of drugs.
in general, the number of drugs remains relatively low.
the sffs algorithm has a single step runtime of |k|, making it linearly increasing with the number of kinase targets.
this number is often approximately  <dig>  the total computational cost of selecting a minimizing target set is o. it should be noted this algorithm is extremely parallelizable, and as such adding additional processors allows the effect of the addition of the numerous kinase targets to be computed significantly faster.
target combination sensitivity inference from a selected target set
in this subsection, we present algorithms for prediction of drug sensitivities when the binarized targets of the test drugs are provided.
the inputs for the algorithms in this subsection are the binarized drug targets, drug sensitivity score and the set of relevant targets for the training drugs.
construction of the target set that solves eq.
<dig> provides information concerning numerically relevant targets based on the drug screen data.
however, the resulting model is still limited in its amount of information.
given the binning behavior of the target selection algorithm, the predicted sensitivity values will include only those for which experimental data is provided, and again only a subset of those target combinations.
hence, in order to expand the current model from one of explanation to one that includes prediction, inferential steps have to be applied using the available information.
the first step in inference is prediction of sensitivity values for target combinations outside the known dataset.
consider that the set of drug representations, , consists of c unique elements.
in addition, the number of targets added to the minimizing target set is |t|=n.
the total possible target combinations is then 2n for binarized target inhibition, and there are thus 2n−c unknown target combination sensitivities.
we would like to be able to perform inference on any of the 2n−c unknown sensitivity combination, and we would like to utilize known sensitivities whenever possible.
to begin the inference step, let us first recall the  <dig> complementary rules for kinase target behavior upon which we base this model.
rule 1: if  is the inhibiting set of targets for drug i and the drug is successful in inhibiting the circuit, then any set b such that si⊂b will also be successful in inhibiting the circuit.
rule 2: if  is the inhibiting set of targets for drug i and the drug is unsuccessful in inhibiting the circuit, then any set b such that b⊂si will also be unsuccessful in inhibiting the circuit.
we consider a third rule that expresses target combination behavior as a function of its most similar target combinations.
rule 3: if  is the inhibiting set of a target combination with unknown sensitivity, then the sensitivity of  will be at least that of ’s closest subsets and will be at most ’s closest supersets.
rule  <dig> follows from the first two rules; rule  <dig> provides that any superset will have greater sensitivity, and rule  <dig> provides that any subset will have lower sensitivity.
to apply rule  <dig> in practical situations, we must guarantee that every combination  will have a subset and superset with an experimental value.
we will assume that the target combination that inhibits all targets in t will be very effective, and as such will have sensitivity  <dig>  in addition, the target combination that consists of no inhibition of any target, which is essentially equivalent to no treatment of the disease, will have no effectiveness, and as such will have a sensitivity of  <dig>  either of these can be substituted with experimental sensitivity values that have the corresponding target combination.
in numerous practical scenarios, the target combination of no inhibition has sensitivity  <dig>
with the lower and upper bound of the target combination sensitivity fixed, we now must perform the inference step by predicting, based on the distance between the subset and superset target combinations.
we perform this inference based on binarized inhibition, as the inference here is meant to predict the sensitivity of target combinations with non-specific ec <dig> values.
refining sensitivity predictions further based on actual drugs with specified ec <dig> values will be considered later.
let =  be the target combination of the subset of  with the highest sensitivity, and let = , the superset target combination with the lowest sensitivity.
let the sensitivity of  and  be yl and yu respectively.
let the hamming distance between cl and cu be h=⊕, and the hamming distance between  and  be d=⊕.
therefore, to transition from  to , it will require the inhibition of an additional d targets, denoted {t <dig> t <dig> ⋯,td}, and the remaining h−d, denoted {td+ <dig> ⋯,th} targets will remain uncontrolled.
for naive inference, we can consider that over the course of the addition of the h targets needed to transition from  to , the change in sensitivity due to the addition of each target is uniform.
with  as the lower bound of the drug sensitivity, the resulting naive sensitivity from the addition of d2≤h targets is
yi=yl+·d2hn
where n is a tunable inference discount parameter, where decreasing n increases yi and presents an optimistic estimate of sensitivity.
we can extend the sensitivity inference to a non-naive approach.
suppose for each target ti∈t, we have an associated target score αi.
the score can be derived from prior knowledge or pre-modeling analysis.
given this vector α= , we will define yi as follows:
yi=yl+·α+α+⋯+αα+α+⋯+α+α+⋯+αn
which can be written as
yi=yl+·∑i=1dα∑i=1hαn
as desired, if the majority of the mass of the weights of t <dig> t <dig> ⋯,th rest in t <dig> t <dig> ⋯,td, the sensitivity of yi will be close to yu.
with the inference function defined as above, we can create a prediction for the sensitivity of any binarized kinase target combination relative to the target set t; thus we can infer all of 2n−c unknown sensitivities from the experimental sensitivities, creating a complete map of the sensitivities of all possible kinase target-based therapies relevant for the patient.
as noted previously, this complete set of sensitivity combinations constitutes the tim.
the tim effectively captures the variations of target combination sensitivities across a large target set.
however, we also plan to incorporate inference of the underlying nonlinear signaling tumor survival pathway that acts as the underlying cause of tumor progression.
we address this using the tim sensitivity values and the binarized representation of the drugs with respect to target set, .
generation of tim circuits
in this subsection, we present algorithms for inference of blocks of targets whose inhibition can reduce tumor survival.
the resulting combination of blocks can be represented as an abstract tumor survival pathway which will be termed as the tim circuit.
the inputs for this subsection are the inferred tim from previous subsection and a binarization threshold for sensitivity.
the output is a tim circuit.
consider that we have generated a target set t for a sample cultured from a new patient.
with the ability to predict the sensitivity of any target combination, we would like to use the available information to discern the underlying tumor survival network.
due to the nature of the functional data, which is a steady-state snapshot and as such does not incorporate changes over time, we cannot infer models of a dynamic nature.
we consider static boolean relationships.
in particular, we expect two types of boolean relationships: logical and relationships where an effective treatment consists of inhibiting two or more targets simultaneously, and logical or relationships where inhibiting one of two or more sets of targets will result in an effective treatment.
here, effectiveness is determined by the desired level of sensitivity before which a treatment will not be considered satisfactory.
the two boolean relationships are reflected in the  <dig> rules presented previously.
by extension, a not relationship would capture the behavior of tumor suppressor targets; this behavior is not directly considered in this paper.
another possibility is xor  and we do not consider it in the current formulation due to the absence of sufficient evidence for existence of such behavior at the kinase target inhibition level.
thus, our underlying network consists of a boolean equation with numerous terms.
to construct the minimal boolean equation that describes the underlying network, we utilize the concept of tim presented in the previous section.
note that generation of the complete tim would require 2n−c≈2n inferences.
the inferences are of negligible computation cost, but for a reasonable n, the number of necessary inferences can become prohibitive as the tim is exponential in size.
we assume that generating the complete tim is computationally infeasible within the desired time frame to develop treatment strategies for new patients.
thus, we fix a maximum size for the number of targets in each target combination to limit the number of required inference steps.
let this maximum number of targets considered be m.
we then consider all non-experimental sensitivity combinations with fewer than m+ <dig> targets.
as we want to generate a boolean equation, we have to binarize the resulting inferred sensitivities to test whether or not a target combination is effective.
we denote the binarization threshold for inferred sensitivity values by θ∈  <cit> .
as θi→ <dig>  an effective combination becomes more restrictive, and the resulting boolean equations will have fewer effective terms.
there is an equivalent term for target combinations with experimental sensitivity, denoted θe.
we begin with the target combinations with experimental sensitivities.
for converting the target combinations with experimental sensitivity, we binarize those target combinations, regardless of the number of targets, where the sensitivity is greater than θe.
the terms that represent a successful treatment are added to the boolean equation.
furthermore, the terms that have sufficient sensitivity can be verified against the drug representation data to reduce the error.
to find the terms of the network boolean equation, we begin with all possible target combinations of size  <dig>  if the sensitivity of these single targets are sufficient relative to θi and θe, the target is binarized; any further addition of targets will only improve the sensitivity as per rule  <dig>  thus, we can consider this target completed with respect to the equation, as we have created the minimal term in the equation for the target.
if the target is not binarized at that level, we expand it by including all possible combinations of two targets including the target in focus.
we continue expanding this method, cutting search threads once the binarization threshold has been reached.
algorithm  <dig> algorithm for generation of minimal boolean equation
the method essentially resembles a breadth or depth-first search routine over n branches to a maximum depth of m. this routine has time complexity of o, and will select the minimal terms in the boolean equation.
the d term results from the cost of a single inference.
the time complexity of this method is significantly lower than generation of the complete tim and optimizing the resulting tim to a minimal boolean equation.
for the minimal boolean equation generation algorithm shown in algorithm  <dig>  let the function binary return the binary equivalent of x given the number of targets in t, and let sensitivity return the sensitivity of the inhibition combination x for the target set t.
with the minimal boolean equation created using algorithm  <dig>  the terms can be appropriately grouped to generate an equivalent and more appealing minimal equation.
to convey the minimal boolean equation to clinicians and researchers unfamiliar with boolean equations, we utilize a convenient circuit representation, as in figures  <dig> and  <dig>  these circuits were generated from two canine subjects with osteosarcoma, as discussed in the results section.
the circuit diagrams are organized by grouped terms, which we denote as blocks.
blocks in the tim circuit act as possible treatment combinations.
the blocks are organized in a linear or structure; treatment of any one block should result in high sensitivity.
as such, inhibition of each target results in its line being broken.
when there are no available paths between the beginning and end of the circuit, the treatment is considered effective.
as such, each block is essentially a modified and/or structure.
within the blocks, parallel lines denote an and relationship, and adjacent lines represent an or relationship.
the goal of an effective treatment then, from the perspective of the network circuit diagram, is to prevent the tumor from having a pathway by which it can continue to grow.
discussion
in this section, we discuss extensions of the tim framework presented earlier.
we provide foundational work for incorporating sensitivity prediction via continuous-valued analysis of ec <dig> values of new drugs as well as theoretical work concerning dynamical models generated from the steady state tims developed previously.
incorporating continuous target inhibition values
the analysis considered in the earlier sections was based on discretized target inhibition i.e.
each drug was denoted by a binary vector  representing the targets inhibited by the drug.
the framework can predict the sensitivities of new drugs with high accuracy as illustrated by the results on canine osteosarcoma tumor cultures.
however, the current framework can also be modified to incorporate the continuous nature of target inhibition and application of different concentrations of a new drug.
let us consider that a drug i with target set t <dig> and ec <dig> profile ei, <dig> ei, <dig> ⋯,ei,n is applied at concentration x nm.
for each ec <dig> value ei,j, we can fit a hill curve or a logistic function to estimate the inhibition of target j at concentration x nm.
for instance a logistic function will estimate the inhibition of target j as f=1/).
note that at concentration x=ei,j, f= <dig>  as desired.
this approach can be applied to arrive at a continuous target profile zi, <dig> zi, <dig> ⋯,zi,n of a drug that is dependent on the applied drug concentration.
the zi,j’s denote real numbers between  <dig> and  <dig> representing the inhibition ratio of target j. this approach can also be applied to generate drug target profiles for a combination of drugs at different concentrations.
to arrive at the sensitivity prediction for a new target inhibition profile, we can apply rules similar to rules  <dig>   <dig> and  <dig> along with searching for closest target inhibition profiles among the training data set.
the block analysis performed using discretized target inhibitions can provide smaller sub-networks to search for among the target inhibition profiles.
incorporating network dynamics in the tim formulation
the tim developed in the previous sections is able to predict the steady state behavior of target inhibitor combinations but cannot provide us with the dynamics of the model or the directionality  of the tumor pathways.
this limitation is a result of the experimental drug perturbation data being from the steady state.
our results show that the proposed approach is highly successful in locating the primary faults in a tumor circuit and predict the possible sensitivity of target combinations at the current time point.
however, extension of this model to incorporate the directional pathways will require protein or gene expression measurements.
the extension refers to steps f <dig> and f <dig> in figure  <dig>  these steps are not necessary to design the control policy but if performed can provide superior performance guarantees.
if we plan to infer a dynamic model from no prior knowledge, the number of required experiments will be huge and will primarily require time- series gene or protein expression measurements.
in this section, we will show that the circuit produced by our tim approach can be used to significantly reduce the search space of directional pathways.
to arrive at the potential dynamical models satisfying the inferred tim, we will consider the possible directional pathways that can generate the inferred tim and convert the directional pathways to discrete boolean network   <cit>  models.
the tim can be used to locate the feasible mutation patterns and constrain the search space of the dynamic models generating the tim.
for the duration of the network dynamics analysis, we will consider the two dynamic models shown in figure  <dig>
directional pathway to bn
to generate a discrete dynamical boolean network  model  <cit>  of a directional pathway, we will first consider the starting mutations or latent activations.
the number of states in the bn will be 2n+ <dig> for n targets.
each state will have n+ <dig> bits with first n bits referring to the discrete state of the n targets and the least significant bit  will correspond to the binarized phenotype ie.
tumor  or normal .
the rules of state transition are  a target state at time t+ <dig> becomes  <dig> if any immediate upstream neighbor has state  <dig> at time t for or relationships or all immediate upstream neighbors have state  <dig> at time t for and relationships.
note that the examples have or type of relations as they are the most commonly found relations in biological pathways .
for the bn without any drug, the targets that are mutated or have latent activations will transition to state  <dig> within one time step.
for a target with no inherent mutation or latent activation, the state will become  <dig> at time t+ <dig> if the immediate upstream activators of the target has state  <dig> at time t.
let us consider the simple example of a biological pathway shown in figure <dig>  the downstream target k <dig> can be activated by either of the upstream targets k <dig> or k <dig>  the tumor is in turn caused by the activation of k <dig>  for this directional pathway, we will assume that k <dig> and k <dig> are activated by their own mutations or have latent activations.
the corresponding bn transition diagram for this pathway is shown in figure  <dig>  for instance, if we consider the state  <dig> at time t, it denotes k <dig>  k <dig> being inactive and k <dig> being active and the phenotype being non-tumorous.
based on the directional pathway in figure  <dig>  activation of k <dig> causes tumor and thus the phenotype will change to tumor  at t+ <dig>  we are given that only k <dig> and k <dig> have mutations or latent activations, thus the activation k <dig> cannot be maintained without the activation of either k <dig> or k <dig> and thus we will have k3= <dig> at t+ <dig>  however, since k <dig> and k <dig> have mutations or latent activations, they will become  <dig> at time t+ <dig> which in turn will activate k <dig> at time t+ <dig>
dynamical model following target inhibition
the bn in figure  <dig> can also be represented by a 16× <dig> transition matrix q representing the state transitions.
to generate the dynamic model after inhibition of a specific target set s <dig> , we should consider that the transition i→j in the un-treated system will be converted to i→z in the treated system where z differs from j only in the target set s <dig> and all targets in s <dig> have value  <dig> for z. each target inhibition combination can be considered as multiplying a matrix tc to the initial transition matrix q. each row of tc contains only one non-zero element of  <dig> based on how the inhibition alters the state.
if we consider n targets, ntc’s in combination can produce a total of 2n possible transformation matrices t <dig> t <dig> ⋯,t2n.
the tim denotes the state of the lsb of the attractor for the 2n transition matrices t1q,t2q,⋯,t2nq starting from initial state 11⋯ <dig> .
for instance, if we consider that our drug inhibits the target k <dig> , the discrete dynamic model following application of the drug is shown in figure  <dig>  we should note that the equilibrium state of the network  <dig> has  <dig> for the tumor state.
this is because the tumor is activated by k <dig> and inhibition of k <dig> should eradicate the tumor.
on the other hand, since both k <dig> and k <dig> can cause tumor through activation of intermediate k <dig>  inhibition of only one of k <dig> and k <dig> will not block the tumor.
the bn following inhibition of k <dig> is shown in figure  <dig> where the attractor  <dig> denotes a tumorous phenotype.
experiment design to infer the dynamic pathway structure
the tim can be used to produce possible dynamic models based on assumptions of latent activations or mutations.
for instance, knowledge of the steady-state value of the target k <dig> following application of target inhibitor for k <dig>  will remove one of the possibilities.
following inhibition of k <dig>  the value of k <dig> will remain  <dig> for the case of figure  <dig> as k <dig> is upstream of k <dig>  conversely, the value of k <dig> will be  <dig> for the second case as k <dig> activates k <dig>
in the following paragraphs, we will consider a general pathway obtained from a tim having the structure shown in figure  <dig> but with unknown directionalities of the blocks and target positions.
for the current analysis, we will assume that there are no common targets in distinct blocks.
we will consider that the pathway has l blocks in series  and each block bi has ai parallel lines with each line j containing bji targets .
the total number of targets in the general map is ∑i=1l∑j=1aibji.
assuming that the n targets are distinct, the maximum number of distinct discrete dynamic models satisfying the structure is l!∏i=1l∏j=1ai!.
each parallel line j in block i can have bji!
possible directional orientations.
if the figure  <dig> represents a possible directional orientation, then only the targets k <dig> ,k <dig> ⋯ka <dig>  will have initial activations due to mutations or latent activations.
some other downstream target cannot have a mutation or latent activation otherwise the target inhibition combination k <dig> ,k <dig> ⋯ka <dig>  will not be effective.
for our analysis, we are assuming that we can inhibit specific targets of our choice and we can measure the steady state target expression following application of the target inhibitions.
we can locate the directionality of the blocks b <dig> to bl by using at most l− <dig> steady state measurements.
we can start by randomly picking any block bi and blocking the targets in that block, the blocks that will remain activated will be upstream of that block and the blocks that will be deactivated following the inhibition of block bi will be located down-stream of bi.
note that the number of experiments required is based on steady state measurements following particular perturbations.
time series measurements can reduce the number of experiments required but may not be always technically feasible.
the expected number of experiments required to detect the directionality of l serial blocks is
2l−13forl≥ <dig>
the next step will be locating the directionality of targets in each parallel line of the block.
we can start with an experiment where for each block bi, one target from each line up to a maximum of ai− <dig> lines will be inhibited.
we cannot inhibit all the lines in a block or else the downstream blocks will also be inhibited and no inference can be made on those blocks for that experiment.
while locating the directionality of the serial blocks bi, we have already validated the position of one target from each parallel line in a serial block.
if we consider a single block bi, each experiment can detect the location of ai− <dig> targets, thus the total number of experiments required to decipher the possible directionalities of the targets in the block bi is ≤max where si={ <dig> ⋯,ai}.
thus for the overall map, the worst case number of experiments new required to decipher the directionalities of all the targets is upper-bounded by
new≤maxi∈s{max}+l− <dig>
where s={ <dig> ⋯,l}.
utilizing equation  <dig>  the expected number of experiments nea required to decipher the directionalities of all the targets is upper-bounded by
nea≤maxi∈s{max⌉−˘1)}+2l− <dig>
CONCLUSIONS
in this article, we presented a novel framework for predicting the effectiveness of molecularly targeted drugs.
we used drug perturbation data to generate a map of the underlying genetic regulatory pathway.
using actual experimental data, we were able to show the effectiveness of our approach for drug sensitivity prediction.
the proposed tim approach produced a low average leave one out cross validation error of 5% when applied to perturbation data generated from four primary canine tumors using a set of  <dig> drugs.
we should note that the current  <dig> drug screen is a small one and technology has been developed for drug screens with a far greater number of drugs.
we are currently experimenting with pharmaceutical drug library consisting of more than  <dig> small molecule inhibitors.
we expect that the use of larger number of drugs will increase the accuracy further and generate maps with greater robustness.
the scope of the present article is concentrated around steps b, c and d of figure  <dig>
for future research, we will consider multiple data sources to increase the robustness of the designed maps.
as explained in figure  <dig>  we can use rapid sirna screens to validate single points of failures predicted by our tim approach.
furthermore, rnaseq and protein phosphoarray data can be used to further revise the circuit.
finally, time series data can be used to incorporate dynamics in the modeling framework.
for combination therapy design, we can use the tim framework to formulate control strategies with various constraints.
some possibilities are  minimal toxicity,  anticipating evolving drug resistance, and  success over a family of tims representing variations of a tumor.
for case , we can assume that the toxicity of a drug or drug combination is proportional to the number of targets being inhibited by the drug and search for the drug combination with high sensitivity but low set of target inhibitions.
for case , we would want to avoid resistance and thus would like to inhibit more than one independent blocking pathway such that for the scenario when resistance to one of the blocking pathways develops, the other independent pathway can still keep the tumor under check.
in other words, we would be interested in selecting a set of targets that can be divided into two or more non-intersecting sets such that the sensitivity of each set is higher than a threshold.
for case , the goal is to design control policies for the scenario when the exact pathway is not known but it belongs to a collection of pathways.
the uncertainty can arise when the experimental data is not sufficient enough to produce a unique pathway map or the current pathway may evolve into one of the different pathways obtained from tissues with same type of cancer.
this can approached from a worst case perspective  <cit>  or a bayesian perspective  <cit> .
in conclusion, the proposed framework provides a unique input-output based methodology to model a cancer pathway and predict the effectiveness of targeted drugs.
this framework can be developed as a viable approach for personalized cancer therapy.
to aide in the usage of our framework, we have developed a graphical user interface which implements in an easy to use way the algorithms and equations presented in this paper.
it is built in matlab, but is distributed as a compiled executable; as such, it is usable in a windows environment by downloading the matlab compile runtime  environment, which is free to download and requires no matlab installation.
it is available online at: http://cvial.ece.ttu.edu/ranadippal/research.html under the target inhibition map approach to inference of cancer pathways heading.
competing interests
author ck has received honoraria for scientific presentations at novartis, millennium/takeda pharmaceutical and glaxosmithkline, and has research joint ventures or sponsored research with novartis and johnson & johnson.
ck is also a paid consultant to the nci/ctep pediatric preclinical testing program .
ck is also a co-founder/stockholder of numira biosciences.
the remaining authors have no other conflicts to declare related to these studies except that a ec is a first degree relative of a bms scientist.
authors’ contributions
designed the algorithmic framework: nb, rp.
implemented the algorithms: nb.
performed the experiments: ld, ec, bs, ck.
analyzed the results: nb, rp.
all authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary analysis.
the file includes a detailed tim example, results of a set of synthetic experiments, theoretical analysis of errors and additional tim circuits.
click here for file
additional file 2
drug screenic <dig> dataset.
the experimentally measured ic50’s in nm for four cell cultures  across  <dig> drugs.
click here for file
additional file 3
drug screen panelec <dig>  the ec <dig> values in nm for  <dig> targets for  <dig> drugs.
click here for file
additional file 4
cross validation results.
the detailed results of the  <dig> fold cross validation error estimates for the tumor cultures bailey, charley, sy and cora.
click here for file
additional file 5
similarity between drugs.
the similarities between the  <dig> targeted drugs as measured by similarity measure Δ.
click here for file
acknowledgements
this work was supported in part by nsf grant ccf <dig> and nci grant 5r01ca133229-05/ <dig>

 in the last decade, a number of drugs targeting specific biologically relevant kinases have been developed that are becoming common in cancer research as a basis for personalized therapy. the idea of treating cancer through inhibition of a specific tyrosine kinase was proven by the discovery that patients with chronic myeloid leukemia can be successfully treated by inhibiting the tyrosine kinase bcr-abl with the kinase inhibitor imatinib mesylate. however, the success rate of any one specific targeted drug for other forms of cancer, such as sarcoma, is limited as the tumors exhibit a wide variety of signaling pathways and are not uniformly dependent on the activity of a specific kinase. the numerous aberrations in molecular pathways that can produce cancer is one cause to necessitate the use of drug combinations for treatment of individual cancers. combination therapy design requires a framework for inference of the individual tumor pathways, prediction of tumor sensitivity to targeted drug and algorithms for selection of the drug combinations under different constraints. the current state of the art in predicting sensitivity to drugs is primarily based on assays measuring gene expression, protein abundance and genetic mutations of tumors; these methods often have low accuracy due to the breadth of available expression data coupled with the absence of information on the functional importance of many genetic mutations. a commonly used method for predicting the success of targeted drugs for a tumor sample is based on the genetic aberrations in the tumor. however, the accuracy of prediction of drug sensitivity based on mutation knowledge is limited in many forms of tumors as some of the mutations may not be functionally important or tumors can develop without the known genetic mutations. statistical tests have been used in to show that genetic mutations can be predictive of the drug sensitivity in non-small cell lung cancers but the classification rates of these predictors based on individual mutations for the aberrant samples are still low. for specific diseases, some mutations have been able to predict the patients that will not respond to particular therapies: for instance reports a success rate of 87% in predicting non-responders to anti-egfr monoclonal antibodies using the mutational status of kras, braf, pik3ca and pten. the prediction of tumor sensitivity to drugs has also been approached as a classification problem using gene expression profiles. in, gene expression profiles are used to predict the binarized efficacy of a drug over a cell line with the accuracy of the designed classifiers ranging from 64% to 92%. in, a co-expression extrapolation approach is used to predict the binarized drug sensitivity in data points outside the training set with an accuracy of around 75%. in, a random forest based ensemble approach was used for prediction of drug sensitivity and achieved an r value of between the predicted ic50s and experimental ic50s. supervised machine learning approaches using genomic signatures achieved a specificity and sensitivity of higher than 70% for prediction of drug response in. tumor sensitivity prediction has also been considered as a drug-induced topology alteration using phospho-proteomic signals and prior biological knowledge of a generic pathway and a molecular tumor profile based prediction. most interestingly, in the recent cancer cell line encyclopedia study, the authors characterize a large set of cell lines with numerous associated data measurement sets: gene and protein expression profiles, mutation profiles, methylation data along with the response of around of these cells lines across anti-cancer drugs. one of the goals of the study was to enable predictive modeling of cancer drug sensitivity. for generating predictive models, the authors considered regression based analysis across input features of gene and protein expression profiles, mutation profiles and methylation data. the performance of the predictive models using fold cross validation ranged between to. in particular, the correlation coefficient for prediction of sensitivity using genomic signatures for the drug erlotinib across > cell lines was <. erlotinib is a commonly used tryosine kinase inhibitor selected primarily as an egfr inhibitor. however, studies have shown that these targeted drugs often have numerous side targets that can play significant roles in the effectiveness of the inhibitor drugs. the target inhibition profiles of drugs and sensitivity of trainings set of drugs can provide significant information for enhanced prediction of anti-cancer drug sensitivity as we have recently shown. by incorporating the drug-target interaction data and sensitivities of training drugs with genomic signatures, we were able to achieve a correlation coefficient of for prediction of erlotinib sensitivity using fold cross validation. the result illustrates the fundamental concept of the importance of drug-target interaction and functional data under which we develop the sensitivity prediction method presented in this paper. by developing a framework around the functional and target information extracted from the primary tumor drug screen performed by our collaborators, we seek to develop a cohesive approach to sensitivity prediction and combination therapy design. this necessitates the generation of the tumor pathway structure for individual patients to decide on the target inhibitors for therapy based on the personalized patient pathways. we envision that the overall schematic of the design of personalized pathways and personalized therapy will be similar to the workflow shown in figure the explanations of the various steps in the design process are as follows: a patient is diagnosed with cancer and a primary culture of the tumor is established. cell viability after exposure to targeted drugs is measured through a drug screen. use of this functional data rather than mutation or protein biomarkers provides a unique advantage. a target inhibition map is generated based on the ic50′s and the known targets of the drugs in the screen. tim denotes a predictive model that provides the sensitivity for all possible target inhibitions. specifically, a tim is composed of a set t={t t ⋯,tn} consisting of binary variables, each denoting inhibition of a target, and a function f relating the target inhibitions to the steady state sensitivity yt, i.e. yt=f. the inhibition vector corresponding to a drug is known as the drug target inhibition profile. a detailed example of tim is provided in additional file the coarse structure of the tim can also be represented by an abstract pathway which will be termed tim circuit. the construction of the tim circuit is explained in the methods section. further data is collected using sirna screens, rna sequencing and protein phosphoarrays to reduce model parameter uncertainties. based on the knowledge of the tim and tim-directed protein expression measurements, the dynamic model is created. combination therapy is designed utilizing the personalized tim and the dynamic model. various constraints such as avoiding resistance to drugs or minimizing toxicity can be applied to design the combination therapy. a mouse xenograft model can be used to study development of resistance simultaneously. the generated drug combinations are validated in vitro on the primary culture. if needed, the circuit is revised or the drug combination with best response in vivo and in vitro is then provided to the patient. the primary contributions of this paper are: methods for extraction of numerically relevant drug targets from single-run drug screens, design of the personalized tim circuit based on drug perturbation data, algorithms for sensitivity prediction of a new drug or drug cocktail, validation over canine osteosarcoma primary tumors and pathway flow inference using sequential protein expression measurements. the scope of the present article is concentrated around steps b, c and d of figure the perturbation data required for our proposed method originates from a drug screen consisting of small molecule inhibitors with quantified kinase interaction behaviors. this drug screen, denoted drug screen version, consists of two sets of data: the first set is the experimentally generated drug sensitivities provided as 50% inhibitory concentration values. the ic values denote the amount of a drug required to reduce the population of cancerous cells in vitro by half. the sensitivity values are expected to change during each new cell line/ tumor culture experiment. the generation of the sensitivities in step c can be done within hours of initial biopsy using drug sensitivity assays which is a period of limited cell divisions for most primary cultures. thus, the estimated personalized maps may be closer to real-time circuits in cancer cells - akin to the signaling found in an untreated patient within a day or two after biopsy, and not the evolving consensus pattern of signaling for growing and dividing tumor cells as subpopulations emerge with increased fitness in vitro. in addition, the drug screen contains experimentally derived half-maximal concentration values for the interaction of each drug and each kinase target. the ec value is directly related to the notion of inhibition of a kinase target; in particular, the ec values correspond to the amount of a compound needed to deactivate via phosphorylation 50% of the population of the associated target. hence, for a drug compound, a target with a lower ec is the one that will be heavily inhibited at low drug concentration levels. thus, low ec targets are often considered to be the primary targets of a drug. the remaining targets are considered to be the side targets of a drug, and are often ignored. the utility of this ec data is its consistency throughout experiments; the ec values as curated from literature searches are fixed, regardless of change of tumor type or patient of origin. this provides a great amount of prior information for analysis of the drug screen results, and its usage is supported from the experiments performed in. the overall goal of the methods presented in this paper is to create an input-output mathematical framework for the analysis of and inference on the functional data generated by the drug screens for the purpose of anti-cancer drug sensitivity prediction and inference of personalized tumor survival pathway. the personalized tumor survival pathway refers to the visual circuit diagram generated from the inferred target inhibition map as explained in the methods section. note that the circuit corresponding to a tim is only a coarse representation of the tim for visual understanding of the most probable target combinations whose inhibition can reduce the tumor survival. since the experiments were conducted on in-vitro cell cultures with the output being cell viability measured in terms of ic the survival here refers to tumor cell culture survival and not the overall survival of the patient.  tim generation for canine osteosarcoma tumor cultures and cross-validation estimates of prediction accuracy the sensitivity prediction and circuit analysis performed on actual biological data are validations of the proposed methodology to be described in the methods section. the experimental data on four tumor cultures and targeted drug screen panel were generated in the keller laboratory at ohsu. the cell lines applied to the drug screen were four canine osteosarcoma cell lines cultured from four distinct canines, denoted bailey, charley, sy, and cora. the tumor cultures were collected by dr. bernard seguin of oregon state university from canines that are part of an ongoing clinical trial for osteosarcoma. the tumor samples were collected from client-owned animals that have developed the disease naturally. all procedures performed on these animals with regards to tumor collection were strictly for treatment purposes and nothing was done different because of the drug perturbation study. all procedures were performed according to standard of care regardless of whether an animal had its tumor sampled. for the generation of the experimental data, the canine osteosarcoma primary cell cultures were plated in well plates at a seeding density of cells per well over graded concentrations of small-molecule kinase inhibitors. each inhibitor was plated individually at four concentrations predicted to bracket the ic for that drug. cells were cultured in rpmi supplemented with 2mm glutamine, 2mm sodium pyruvate, 2mm hepes, 1% penicillin streptomycin, and 10% fetal bovine serum for hours. at the end of the hour incubation, cell viability was assessed using the mts assay. all values were normalized to the mean of seven wells on each plate containing no drug. the ic for each drug was then determined by identification of the two concentrations bracketing 50% cell viability and application of the following formula: +da where cell viability value above 50% = a and cell viability value below 50% = b. the experimentally generated ic values are included as additional file the experimentally generated sensitivities of the drugs are then scaled to values between and among the drugs on the drug screen, drugs have known target inhibition profiles; of these drugs, provide information only on the target mtor and analysis of these drugs are trivial. thus, the remaining drugs are used to generate the tims. these target profiles were extracted from several literature sources based on experimental quantitative dissociation constants which are treated as ec values for each drug across kinase target assays with more than targets. the target profiles of the drugs are shown in additional file figures and represent the equivalent tim circuits generated from experimental data for bailey and sy respectively. the tim circuits for charley and cora are included in additional file to emphasize the biological relevance provided by the tim framework employed in the analysis of the biological data, we present a more in-depth analysis of the tim circuit devised for the canine patient bailey. the vast majority of human osteosarcomas contain genetic or post-translational abnormalities in one or both of the tumor suppressors p and prb. the first target identified in this circuit is pkc alpha. pkc alpha modifies cdkn1a, which is the primary mediator of p tumor suppressor activity. psmb represents the proteasome. previous studies and early preclinical data from the keller laboratory confirms in vitro sensitivity of many osteosarcomas to proteasome inhibitors and this sensitivity is hypothesized to be due to the integral role of the proteasome in p regulation. interestingly, cdk is also prominent in this circuit, which is a primary inhibitor of the tumor suppressor prb, which is also frequently abnormal in spontaneous human osteosarcoma. cdk is an important modifier of both p and prb and is also represented in this circuit. the importance of pi3k pathway in osteosarcoma has also been recently reported using high throughput genotyping. our tim circuit includes akt which is downstream of pi3k. also, ednra selected in the circuit has been known to interact with pkc and activate erk signaling. if the circuit models shown in figures and are used to predict sensitivities for comparison with experimentally generated data, we will get optimistic results as the models are trained using the entirety of the available data. thus, we utilize leave one out and 10-fold cross validation approaches to test the validity of the tim framework that we present in this paper. for the loo approach, a single drug among the drugs with known inhibition profiles is removed from the dataset and a tim is built, using the sffs suboptimal search algorithm, from the remaining drugs. the resulting tim is then used to predict the sensitivity of the withheld drug. the predicted sensitivity value is then compared to it’s experimental value; the loo error for each drug is the absolute value of the experimental sensitivity y minus the predicted sensitivity, i.e. |y−|. the closer the predicted value is to the experimentally generated sensitivity, the lower the error for the withheld drug. tables and provides the complete loo error tables and the average loo error for each primary culture. the average loo error over the cell cultures is or %. for the 10-fold cross validation error estimate, we divided the available drugs into random sets of similar size and the testing is done on each fold while being trained on the remaining folds. this is repeated times and average error calculated on the testing samples. we again repeated this experiment times and the average of those mean absolute errors for the primary cell cultures are shown in table the detailed results of the 10-fold cross validation error analysis are included in additional file we note that both 10-fold cv and loo estimates for all the cultures have errors less than 9%, which is extremely low, especially considering the still experimental nature of the drug screening process performed in the keller laboratory and the available response of only drugs with known target inhibition profile. pred. sens. = predicted sensitivity, exp. sens. = experimental sensitivity, error = |pred. sens. −exp. sens.|, avg err = average of the leave one out errors. pred. sens. = predicted sensitivity, exp. sens. = experimental sensitivity, error = |pred. sens. −exp. sens.|, avg err = average of the leave one out errors. pred. sens. = predicted sensitivity, exp. sens. = experimental sensitivity, error = |pred. sens. −exp. sens.|, avg err = average of the leave one out errors. pred. sens. = predicted sensitivity, exp. sens. = experimental sensitivity, error = |pred. sens. − exp. sens.|, avg err = average of the leave one out errors. the fold- cross valdiation error estimates were calculated for runs and the average mean absolute error across all the drugs is shown. to provide a measure of the overlap between drugs, we considered a similarity measure Λ based on the ec of the drugs d and d let the ec50′s of the drugs d and d be given by the n-length vectors e and e where n denotes the number of drug targets. the entries for the targets that are not inhibited by the drugs are set to let the vectors v and v represent the binarized targets of the drugs i.e. it has a value of if the target is inhibited by the drug and a value of zero if the target is not inhibited by the drug. then, we define the similarity measure Λ as: Λ=∑i=1nmin,e2)∗v1∗v2∑i=1nmax,e2) note that Λ= and similarity between drugs with no overlapping targets is zero. if two drugs have 50% targets overlapping with same ec50′s, then the similarity measure is. the similarities between the drugs are shown in additional file note that except two drugs rapamycin and temsirolimus that have a similarity measure of, all other drugs have significantly lower similarities with each other. the maximum similarity between two different drugs is. this shows that any two drugs in the drug screen are not significantly overlapping and the prediction algorithm is still able to predict the response. the low error rate illustrates the accuracy and effectiveness of this novel method of modeling and sensitivity prediction. furthermore, these error rates are significantly lower than those of any other sensitivity prediction methodology we have found. consistent with the analysis in, the sensitivity prediction rates improve dramatically when incorporating more information about drug-protein interaction. to more effectively compare the results generated via the tim framework with the results in, we also present the correlation coefficients between the predicted and experimental drug sensitivity values in table the correlation coefficients for predicted and experimentally generated sensitivities for drugs and more than cell lines ranges from to when genomic characterizations are used to predict the drug sensitivities in the ccle study. in comparison, our approach based on sensitivity data on training set of drugs and drug-protein interaction information produced correlation coefficients > for both leave one out and 10-fold cross validation approaches for error estimation. ρ loo = average pearson correlation coefficient for test samples selected based on leave one out approach, ρ 10-fold cv = average pearson correlation coefficient for test samples selected based on fold cross validation approach. it should be noted that the sensitivity prediction is performed in a continuous manner, not discretely, and thus effective dosage levels can be inferred from the predictions made from the tim. this shows that the tim framework is capable of predicting the sensitivity to anti-cancer targeted drugs outside the training set, and as such is viable as a basis for a solution to the complicated problem of sensitivity prediction. in addition, we tested the tim framework using synthetic data generated from a subsection of a human cancer pathway taken from the kegg database. here, the objective is to show that the proposed tim method generates models that highly represent the underlying biological network which was sampled via synthetic drug perturbation data. this experiment replicates in synthesis the actual biological experiments performed at the keller laboratory at ohsu. to utilize the tim algorithm, a panel of targeted drugs pulled from a library of is used as a training panel to sample the randomly generated network. additionally, a panel of drugs is drawn from the library to serve as a test panel. the training panel and the testing panel have no drugs in common. each of the training drugs is applied to the network, and the sensitivity for each drug is recorded. the generated tim is then sampled using the test panel which determines the predicted sensitivities of the test panel. the synthetic experiments were performed for randomly generated cancer subnetworks for each of n= ⋯, active targets in the network. the active targets are those which, when inhibited, may have some effect on the cancer downstream. to more accurately mimic the boolean nature of the biological networks, a drug which does not satisfy any of the boolean network equations will have sensitivity a drug which satisfies at least one network equation will have sensitivity the inhibition profile of the test drugs is used to predict the sensitivity of the new drug. the average number of correctly predicted drugs for each n is reported in table this synthetic modeling approach generally produces respectable levels of accuracy, with accuracies ranging from 89% to 99%. drugs for training mimics the drug screen setup used by our collaborators and testing drugs for predicted sensitivity approximates a secondary drug screen to pinpoint optimal therapies. the performance of the synthetic data shows fairly high reliability of the predictions made by the tim approach. correct prediction = average number of correct prediction for testing drugs for different synthetic pathways. we have also tested our algorithm on another set of randomly generated synthetic pathways. the detailed results of the experiment are included in additional file a large number of testing samples were used for each pathway prediction and the results indicate an average error of less than 10% for multiple scenarios. in comparison, the average error with random predictions was 44%. the average correlation coefficient of the prediction to actual sensitivity for the sets of experiments was. the average correlation coefficient with random predictions was we also report the standard deviation of the errors and for a representative example, the percentile of the error was - and percentile, thus the 80% prediction interval for prediction μ was. the results of the synthetic experiments on different randomly generated pathways shows that the approach presented in the paper is able to utilize a small set of training drugs from all possible drugs to generate a high accuracy predictive model. methods in this section, we provide an overview of the model design and inference from drug perturbation data for personalized therapy. mathematical formulation let us consider that we have drug ic data for a new primary tumor after application of m drugs in a controlled drug screen. let the known multi-target inhibiting sets for these drugs be denoted by s s..., sm obtained from drug inhibition studies. the elements of set si are ei= for i=,⋯,m, where ei,j are real-valued elements describing the interaction of si with k=, the set of all kinase targets included in the drug screen. the ei,j’s refer to the ec values discussed previously. it should be noted that for all si, ei,j will most often be blank or an extremely high number denoting no interaction. the initial problem we wish to solve is to identify the minimal subset of k, the set of all tyrosine kinase targets inhibited by the m drugs in the drug panel, which explains numerically the various responses of the m drugs. denote this minimal subset of k as t. the rationale behind minimization of t is twofold. first, as with any classification or prediction problem, a primary goal is avoidance of overfitting. secondly, by minimizing the cardinality of the target set required to explain the drug sensitivities found in the exploratory drug screen, the targets included have supportable numerical relevance increasing the likelihood of biological relevance. additional targets may increase the cohesiveness of the biological story of the tumor, but will not have numerical evidence as support. this set t will be the basis of our predictive model approach to sensitivity prediction. before formulation of the problem for elucidating t, let us consider the nature of our desired approach to sensitivity prediction. from the functional data gained from the drug screen, we wish to generate a personalized tumor survival pathway model instead of a linear function approximator with minimal error. we are working under the fundamental assumption that the tumor survival pathway is nonlinear in its behavior; this assumption is reasonable given the difficulty in treating multiple forms of cancer. one frequent theory in personalized therapy is that effective treatment results from applying treatment across multiple important biological pathways. these pathways generally consist of sequentially activated gene and protein nodes acting as a feedback network. treatment of individual pathways may not be sufficient for majority of diseases, so multiple independent parallel pathways must be targeted to create an effective treatment. we believe that one possible approach to the analysis of multiple pathway treatment is to begin with an underlying framework based on the boolean interactions of the multiple targets in the pathway architecture. the approach is based on developing families of boolean equations that describe the multiple treatment combinations capable of acting as an effective intervention strategy. for the initial step of developing the underlying boolean functions, an initial binarization of the data set must be performed. however, the resulting model lends itself to numerous continuous approaches to sensitivity prediction which we will explore further in the paper. binarization of drug targets and conversion of ic50′s to sensitivities in this subsection, we present algorithms for generation of binarized drug targets and continuous sensitivity score of each drug. the inputs for the algorithms in this subsection are the ec50′s of the drug targets and the ic50′s of the drugs when applied to a tumor culture. in order to perform the binarization, we must consider the nature of the data we are given. in particular, we are provided with an ic for each drug, and an ec value for each kinase target inhibited by the drug. under the assumption that the primary mechanism of tumor eradication is, in fact, the protein kinase inhibition enacted by these targeted drugs, a natural consequence would be the existence of a relationship between the ic and ec values. this relationship is explained as such: suppose for a drug si the ic value of si and the ec of kinase target kj, are of similar value, then it can be reasonably assumed that kinase target kj is possibly a primary mechanism in the effectiveness of the drug. in other words, if 50% inhibition of a kinase target directly correlates with 50% of the tumor cells losing viability, then inhibition of the kinase target is most likely one of the causes of cell death. hence, the target that matches the drug ic is binarized as a target hit for the drug. the above assumption of direct correlation for all successful drugs is obviously an extremely restrictive assumption and will be unable to produce high accuracy predictions. thus, the binarization scheme has to be modified to incorporate the following three factors: first: noises in varying magnitude will be present in the drug screen data generated by our collaborators. the noise is unavoidable, and as such, needs to be accounted for. in addition, despite the high accuracy of the drug-protein interaction data procured from literature, we should still account for possible errors in the ec values for the numerous drugs. second: the restrictive assumption considers that effective drugs operate on single points of failure within the patient’s signaling pathway. in reality, high sensitivity to a drug is often attributed to a family of related kinases or several independent kinases working synergistically over one or multiple pathways to induce tumor death. this cooperative multivariate behavior needs to be taken into account while binarizing a drug to its multiple possible targets. third: despite the high level of currently available knowledge on the biological effects of numerous targeted drugs, there remains the possibility of a drug having high sensitivity while having no known mechanisms explaining its sensitivity. therefore, we must consider the situation where there are latent mechanisms not considered within the dataset that are proving to be effective in some combination of treatment. this point does not necessarily eliminate the possibility of kinase mechanisms being an important factor. we address all three concerns as follows: by considering the log scaled ec values for each target and the log scaled ic value for each drug, we convert the multiplicative noise to additive noise. in addition, we employ scalable bounds around the ic50′s to determine binarization values of the numerous kinase targets for each drug. the bounds can be scaled to allow targets that may have ec50′s higher than the ic to be considered as a possible treatment mechanism. we extend the bounds to low ec levels, and often down to to incorporate the possibility of target collaboration at various different ec levels. while a high ic indicates the likelihood of drug side targets as therapeutic mechanisms, it does not preclude the possibility of a joint relationship between a high ec target and a low ec target. hence, to incorporate the numerous possible effective combinations implied by the ic of an effective drug, the binarization range of targets for a drug is the range α · log≤ log≤β · log where 0≤α≪β. for reliability and validity of the target set that we aim to construct, it is important to keep β in a reasonable range, i.e. β should be a smaller constant such as or for the situation where the above bounds do not result in at least one binarized target, the immediate option is to eliminate the drug from the data set before target selection. this prevents incomplete information from affecting the desired target set. as information concerning the drug screen agents gradually becomes complete with respect to other forms of data, such as gene interaction data, additional mechanisms for unexplained targets can be explored and incorporated more readily into the predictive model. with binarization of the data set as explained, we now present the minimization problem that produces a numerically relevant set of targets, t. consider the target set t=, where ti∈{ }. here, denotes inclusion in the target set t and denotes exclusion. for any target set t one can find the representation under t of each drug si,i∈ ⋯,m as =. as the t will be the basis of the new representation for each drug, this will result in n columns which will be for all si, where n is the number of ti= i.e. the number of targets not included in t the resulting representation of each drug in t is then an n−n vector of ec values. while the representation of each drug will change as the target set t changes, the ic values for each of the m drugs remains the same. these experimental sensitivity values will be used to test the numerous different target sets to quantify the strength of the model for any target set. to simplify scoring of the target set, we first convert the ic for each drug si to a continuous-valued sensitivity score yi∈ yi= ific i<cmaxic∗/log,ifcmax≤x≤maxdosei where maxdosei is the maximum dose of drug si given, cmaxi is the maximum achievable clinical dose of drug si, and c=1−log/log so that the scoring function is continuous. maxdose is used to prevent inferences being made on data that is not available. while it would be possible to attempt interpolation to infer an ic from the multiple available data points, such inference cannot be fully quantified. hence, drugs which fail to achieve an ic within the allotted dosage are given the score of which means ineffective. the cmax value is used to apply a variable score to the numerous drugs based on the inherent toxicity of the drug. this will also prevent bias towards drugs with low ic50’s; some drugs may achieve efficacy at higher levels solely based on the drug ec values. construction of the relevant target set in this subsection, we present approaches for selection of a smaller relevant set of targets t from the set of all possible targets k. the inputs for the algorithms in this subsection are the binarized drug targets and continuous sensitivity score. with the scaled sensitivities, we can develop a fitness function to evaluate the model strength for an arbitrary set of targets. as has been established, for any set of targets t drug si has a unique representation. this representation can be used to separate the drugs into different bins based on the targets it inhibits under t within each of these bins will be several drugs with identical target profiles but different scaled scores. let the set of scores in each bin be denoted y for sj in an arbitrary bin, and we will assign to each bin the mean sensitivity score of the bin, e). denote this value p. within each bin, we want to minimize the variation between the predicted sensitivity for the target combination, p, and the experimental sensitivities, y. this notion is equivalent to minimizing the inconsistencies of the experimental sensitivity values with respect to the predicted sensitivity values for all known target combinations for any set of targets, which in turn suggests the selected target set effectively explains the mechanisms by which the effective drugs are able to kill cancerous cells. numerically, we can calculate the inter-bin sensitivity error using the following equation: ∑bins∑j∈bin|p−y| this analysis has one notable flaw: if we attempt to only separate the various drugs into bins based on inter-bin sensitivity error, we can create an over-fitted solution by breaking each drug into an individual bin. we take two steps to avoid this. first, we attempt to minimize the number of targets during construction of t second, we incorporate an inconsistency term to account for target behavior that we consider to be biologically inaccurate. to expand on the above point, we consider there are two complementary rules by which kinase targets behave. research has shown that the bulk of viable kinase targets behave as tumor promoters, proteins whose presence and lack of inhibition is related to the continued survival and growth of a cancerous tumor. these targets essentially have a positive correlation with cancer progression. this is in opposition to tumor suppressors, proteins that have been shown to have a negative correlation with the development of cancer. to capture the behavior of oncogenes, we partially formulate our problem on two rules: rule 1: if is the inhibiting set of targets for drug i and the drug is successful in inhibiting the circuit, then any set b such that si⊂b will also be successful in inhibiting the circuit. rule 2: if is the inhibiting set of targets for drug i and the drug is unsuccessful in inhibiting the circuit, then any set b such that b⊂si will also be unsuccessful in inhibiting the circuit. rule essentially says that if inhibiting a number of target proteins has blocked signaling pathways, then inhibiting more target proteins will not open any path that has already been blocked. rule captures the fact that if a set of target protein inhibitors is unsuccessful in blocking the paths of a circuit, then any reduced number of target protein inhibitors among the inhibiting proteins cannot block all the paths. the above rules assume that the kinases in focus are oncogenes, genes that promote cancer growth and whose inhibition can prevent tumor development. the majority of kinases in the drug screen panel behave as oncogenes, and as such, our approach utilizes the above rules. target sets resulting in combination scores that do not follow the rule-based behavior incur an inconsistency penalty. this penalty is calculated as follows: ∑drugs∑binsχ|p−y| where χ is the indicator function which is when the experimental drug score is inconsistent with the predicted subset/superset bin score. we now present the complete target set score, and as such, the equation that we wish to solve: mint∑bins∑j∈bin|p−y|+∑drugs∑binsχ|p−y| which reduces to the minimization problem we wish to solve: mint∑drugs∑bins−y|+χ|p−y|) for brevity, we will denote the scoring function of a target set with respect to the binarized ec values s and the scaled sensitivity scores y, Γ. as the s and y sets will be fixed when target set generation begins, we reduce this notation further to Γ. note that t⊆k where k denotes the set of all possible targets. 2|k| is the total number of possibilities for t which is extremely huge and thus prohibits exhaustive search. thus the inherently nonlinear and computational intensive target set selection optimization will be approached through suboptimal search methodologies. a number of methods can be applied in this scenario and we have employed sequential floating forward search to build the target sets. we selected sffs as it generally has fast convergence rates while simultaneously allowing for a large search space within a short runtime. additionally, it naturally incorporates the desired target set minimization aim as sffs will not add features that provide no benefit. we present the sffs algorithm for construction of the minimizing target set in algorithm algorithm sffs algorithm for generating relevant target sets complexity of target set generation the algorithm to generate the error score given a target set t is of order o, quadratic with respect to the number of drugs. in general, the number of drugs remains relatively low. the sffs algorithm has a single step runtime of |k|, making it linearly increasing with the number of kinase targets. this number is often approximately the total computational cost of selecting a minimizing target set is o. it should be noted this algorithm is extremely parallelizable, and as such adding additional processors allows the effect of the addition of the numerous kinase targets to be computed significantly faster. target combination sensitivity inference from a selected target set in this subsection, we present algorithms for prediction of drug sensitivities when the binarized targets of the test drugs are provided. the inputs for the algorithms in this subsection are the binarized drug targets, drug sensitivity score and the set of relevant targets for the training drugs. construction of the target set that solves eq. provides information concerning numerically relevant targets based on the drug screen data. however, the resulting model is still limited in its amount of information. given the binning behavior of the target selection algorithm, the predicted sensitivity values will include only those for which experimental data is provided, and again only a subset of those target combinations. hence, in order to expand the current model from one of explanation to one that includes prediction, inferential steps have to be applied using the available information. the first step in inference is prediction of sensitivity values for target combinations outside the known dataset. consider that the set of drug representations,, consists of c unique elements. in addition, the number of targets added to the minimizing target set is |t|=n. the total possible target combinations is then 2n for binarized target inhibition, and there are thus 2n−c unknown target combination sensitivities. we would like to be able to perform inference on any of the 2n−c unknown sensitivity combination, and we would like to utilize known sensitivities whenever possible. to begin the inference step, let us first recall the complementary rules for kinase target behavior upon which we base this model. rule 1: if is the inhibiting set of targets for drug i and the drug is successful in inhibiting the circuit, then any set b such that si⊂b will also be successful in inhibiting the circuit. rule 2: if is the inhibiting set of targets for drug i and the drug is unsuccessful in inhibiting the circuit, then any set b such that b⊂si will also be unsuccessful in inhibiting the circuit. we consider a third rule that expresses target combination behavior as a function of its most similar target combinations. rule 3: if is the inhibiting set of a target combination with unknown sensitivity, then the sensitivity of will be at least that of ’s closest subsets and will be at most ’s closest supersets. rule follows from the first two rules; rule provides that any superset will have greater sensitivity, and rule provides that any subset will have lower sensitivity. to apply rule in practical situations, we must guarantee that every combination will have a subset and superset with an experimental value. we will assume that the target combination that inhibits all targets in t will be very effective, and as such will have sensitivity in addition, the target combination that consists of no inhibition of any target, which is essentially equivalent to no treatment of the disease, will have no effectiveness, and as such will have a sensitivity of either of these can be substituted with experimental sensitivity values that have the corresponding target combination. in numerous practical scenarios, the target combination of no inhibition has sensitivity with the lower and upper bound of the target combination sensitivity fixed, we now must perform the inference step by predicting, based on the distance between the subset and superset target combinations. we perform this inference based on binarized inhibition, as the inference here is meant to predict the sensitivity of target combinations with non-specific ec values. refining sensitivity predictions further based on actual drugs with specified ec values will be considered later. let = be the target combination of the subset of with the highest sensitivity, and let =, the superset target combination with the lowest sensitivity. let the sensitivity of and be yl and yu respectively. let the hamming distance between cl and cu be h=⊕, and the hamming distance between and be d=⊕. therefore, to transition from to, it will require the inhibition of an additional d targets, denoted {t t ⋯,td}, and the remaining h−d, denoted {td+ ⋯,th} targets will remain uncontrolled. for naive inference, we can consider that over the course of the addition of the h targets needed to transition from to, the change in sensitivity due to the addition of each target is uniform. with as the lower bound of the drug sensitivity, the resulting naive sensitivity from the addition of d2≤h targets is yi=yl+·d2hn where n is a tunable inference discount parameter, where decreasing n increases yi and presents an optimistic estimate of sensitivity. we can extend the sensitivity inference to a non-naive approach. suppose for each target ti∈t, we have an associated target score αi. the score can be derived from prior knowledge or pre-modeling analysis. given this vector α=, we will define yi as follows: yi=yl+·α+α+⋯+αα+α+⋯+α+α+⋯+αn which can be written as yi=yl+·∑i=1dα∑i=1hαn as desired, if the majority of the mass of the weights of t t ⋯,th rest in t t ⋯,td, the sensitivity of yi will be close to yu. with the inference function defined as above, we can create a prediction for the sensitivity of any binarized kinase target combination relative to the target set t; thus we can infer all of 2n−c unknown sensitivities from the experimental sensitivities, creating a complete map of the sensitivities of all possible kinase target-based therapies relevant for the patient. as noted previously, this complete set of sensitivity combinations constitutes the tim. the tim effectively captures the variations of target combination sensitivities across a large target set. however, we also plan to incorporate inference of the underlying nonlinear signaling tumor survival pathway that acts as the underlying cause of tumor progression. we address this using the tim sensitivity values and the binarized representation of the drugs with respect to target set,. generation of tim circuits in this subsection, we present algorithms for inference of blocks of targets whose inhibition can reduce tumor survival. the resulting combination of blocks can be represented as an abstract tumor survival pathway which will be termed as the tim circuit. the inputs for this subsection are the inferred tim from previous subsection and a binarization threshold for sensitivity. the output is a tim circuit. consider that we have generated a target set t for a sample cultured from a new patient. with the ability to predict the sensitivity of any target combination, we would like to use the available information to discern the underlying tumor survival network. due to the nature of the functional data, which is a steady-state snapshot and as such does not incorporate changes over time, we cannot infer models of a dynamic nature. we consider static boolean relationships. in particular, we expect two types of boolean relationships: logical and relationships where an effective treatment consists of inhibiting two or more targets simultaneously, and logical or relationships where inhibiting one of two or more sets of targets will result in an effective treatment. here, effectiveness is determined by the desired level of sensitivity before which a treatment will not be considered satisfactory. the two boolean relationships are reflected in the rules presented previously. by extension, a not relationship would capture the behavior of tumor suppressor targets; this behavior is not directly considered in this paper. another possibility is xor and we do not consider it in the current formulation due to the absence of sufficient evidence for existence of such behavior at the kinase target inhibition level. thus, our underlying network consists of a boolean equation with numerous terms. to construct the minimal boolean equation that describes the underlying network, we utilize the concept of tim presented in the previous section. note that generation of the complete tim would require 2n−c≈2n inferences. the inferences are of negligible computation cost, but for a reasonable n, the number of necessary inferences can become prohibitive as the tim is exponential in size. we assume that generating the complete tim is computationally infeasible within the desired time frame to develop treatment strategies for new patients. thus, we fix a maximum size for the number of targets in each target combination to limit the number of required inference steps. let this maximum number of targets considered be m. we then consider all non-experimental sensitivity combinations with fewer than m+ targets. as we want to generate a boolean equation, we have to binarize the resulting inferred sensitivities to test whether or not a target combination is effective. we denote the binarization threshold for inferred sensitivity values by θ∈. as θi→ an effective combination becomes more restrictive, and the resulting boolean equations will have fewer effective terms. there is an equivalent term for target combinations with experimental sensitivity, denoted θe. we begin with the target combinations with experimental sensitivities. for converting the target combinations with experimental sensitivity, we binarize those target combinations, regardless of the number of targets, where the sensitivity is greater than θe. the terms that represent a successful treatment are added to the boolean equation. furthermore, the terms that have sufficient sensitivity can be verified against the drug representation data to reduce the error. to find the terms of the network boolean equation, we begin with all possible target combinations of size if the sensitivity of these single targets are sufficient relative to θi and θe, the target is binarized; any further addition of targets will only improve the sensitivity as per rule thus, we can consider this target completed with respect to the equation, as we have created the minimal term in the equation for the target. if the target is not binarized at that level, we expand it by including all possible combinations of two targets including the target in focus. we continue expanding this method, cutting search threads once the binarization threshold has been reached. algorithm algorithm for generation of minimal boolean equation the method essentially resembles a breadth or depth-first search routine over n branches to a maximum depth of m. this routine has time complexity of o, and will select the minimal terms in the boolean equation. the d term results from the cost of a single inference. the time complexity of this method is significantly lower than generation of the complete tim and optimizing the resulting tim to a minimal boolean equation. for the minimal boolean equation generation algorithm shown in algorithm let the function binary return the binary equivalent of x given the number of targets in t, and let sensitivity return the sensitivity of the inhibition combination x for the target set t. with the minimal boolean equation created using algorithm the terms can be appropriately grouped to generate an equivalent and more appealing minimal equation. to convey the minimal boolean equation to clinicians and researchers unfamiliar with boolean equations, we utilize a convenient circuit representation, as in figures and these circuits were generated from two canine subjects with osteosarcoma, as discussed in the results section. the circuit diagrams are organized by grouped terms, which we denote as blocks. blocks in the tim circuit act as possible treatment combinations. the blocks are organized in a linear or structure; treatment of any one block should result in high sensitivity. as such, inhibition of each target results in its line being broken. when there are no available paths between the beginning and end of the circuit, the treatment is considered effective. as such, each block is essentially a modified and/or structure. within the blocks, parallel lines denote an and relationship, and adjacent lines represent an or relationship. the goal of an effective treatment then, from the perspective of the network circuit diagram, is to prevent the tumor from having a pathway by which it can continue to grow. discussion in this section, we discuss extensions of the tim framework presented earlier. we provide foundational work for incorporating sensitivity prediction via continuous-valued analysis of ec values of new drugs as well as theoretical work concerning dynamical models generated from the steady state tims developed previously. incorporating continuous target inhibition values the analysis considered in the earlier sections was based on discretized target inhibition i.e. each drug was denoted by a binary vector representing the targets inhibited by the drug. the framework can predict the sensitivities of new drugs with high accuracy as illustrated by the results on canine osteosarcoma tumor cultures. however, the current framework can also be modified to incorporate the continuous nature of target inhibition and application of different concentrations of a new drug. let us consider that a drug i with target set t and ec profile ei, ei, ⋯,ei,n is applied at concentration x nm. for each ec value ei,j, we can fit a hill curve or a logistic function to estimate the inhibition of target j at concentration x nm. for instance a logistic function will estimate the inhibition of target j as f=1/). note that at concentration x=ei,j, f= as desired. this approach can be applied to arrive at a continuous target profile zi, zi, ⋯,zi,n of a drug that is dependent on the applied drug concentration. the zi,j’s denote real numbers between and representing the inhibition ratio of target j. this approach can also be applied to generate drug target profiles for a combination of drugs at different concentrations. to arrive at the sensitivity prediction for a new target inhibition profile, we can apply rules similar to rules and along with searching for closest target inhibition profiles among the training data set. the block analysis performed using discretized target inhibitions can provide smaller sub-networks to search for among the target inhibition profiles. incorporating network dynamics in the tim formulation the tim developed in the previous sections is able to predict the steady state behavior of target inhibitor combinations but cannot provide us with the dynamics of the model or the directionality of the tumor pathways. this limitation is a result of the experimental drug perturbation data being from the steady state. our results show that the proposed approach is highly successful in locating the primary faults in a tumor circuit and predict the possible sensitivity of target combinations at the current time point. however, extension of this model to incorporate the directional pathways will require protein or gene expression measurements. the extension refers to steps f and f in figure these steps are not necessary to design the control policy but if performed can provide superior performance guarantees. if we plan to infer a dynamic model from no prior knowledge, the number of required experiments will be huge and will primarily require time- series gene or protein expression measurements. in this section, we will show that the circuit produced by our tim approach can be used to significantly reduce the search space of directional pathways. to arrive at the potential dynamical models satisfying the inferred tim, we will consider the possible directional pathways that can generate the inferred tim and convert the directional pathways to discrete boolean network models. the tim can be used to locate the feasible mutation patterns and constrain the search space of the dynamic models generating the tim. for the duration of the network dynamics analysis, we will consider the two dynamic models shown in figure directional pathway to bn to generate a discrete dynamical boolean network model of a directional pathway, we will first consider the starting mutations or latent activations. the number of states in the bn will be 2n+ for n targets. each state will have n+ bits with first n bits referring to the discrete state of the n targets and the least significant bit will correspond to the binarized phenotype ie. tumor or normal. the rules of state transition are a target state at time t+ becomes if any immediate upstream neighbor has state at time t for or relationships or all immediate upstream neighbors have state at time t for and relationships. note that the examples have or type of relations as they are the most commonly found relations in biological pathways. for the bn without any drug, the targets that are mutated or have latent activations will transition to state within one time step. for a target with no inherent mutation or latent activation, the state will become at time t+ if the immediate upstream activators of the target has state at time t. let us consider the simple example of a biological pathway shown in figure the downstream target k can be activated by either of the upstream targets k or k the tumor is in turn caused by the activation of k for this directional pathway, we will assume that k and k are activated by their own mutations or have latent activations. the corresponding bn transition diagram for this pathway is shown in figure for instance, if we consider the state at time t, it denotes k k being inactive and k being active and the phenotype being non-tumorous. based on the directional pathway in figure activation of k causes tumor and thus the phenotype will change to tumor at t+ we are given that only k and k have mutations or latent activations, thus the activation k cannot be maintained without the activation of either k or k and thus we will have k3= at t+ however, since k and k have mutations or latent activations, they will become at time t+ which in turn will activate k at time t+ dynamical model following target inhibition the bn in figure can also be represented by a 16× transition matrix q representing the state transitions. to generate the dynamic model after inhibition of a specific target set s, we should consider that the transition i→j in the un-treated system will be converted to i→z in the treated system where z differs from j only in the target set s and all targets in s have value for z. each target inhibition combination can be considered as multiplying a matrix tc to the initial transition matrix q. each row of tc contains only one non-zero element of based on how the inhibition alters the state. if we consider n targets, ntc’s in combination can produce a total of 2n possible transformation matrices t t ⋯,t2n. the tim denotes the state of the lsb of the attractor for the 2n transition matrices t1q,t2q,⋯,t2nq starting from initial state 11⋯. for instance, if we consider that our drug inhibits the target k, the discrete dynamic model following application of the drug is shown in figure we should note that the equilibrium state of the network has for the tumor state. this is because the tumor is activated by k and inhibition of k should eradicate the tumor. on the other hand, since both k and k can cause tumor through activation of intermediate k inhibition of only one of k and k will not block the tumor. the bn following inhibition of k is shown in figure where the attractor denotes a tumorous phenotype. experiment design to infer the dynamic pathway structure the tim can be used to produce possible dynamic models based on assumptions of latent activations or mutations. for instance, knowledge of the steady-state value of the target k following application of target inhibitor for k will remove one of the possibilities. following inhibition of k the value of k will remain for the case of figure as k is upstream of k conversely, the value of k will be for the second case as k activates k in the following paragraphs, we will consider a general pathway obtained from a tim having the structure shown in figure but with unknown directionalities of the blocks and target positions. for the current analysis, we will assume that there are no common targets in distinct blocks. we will consider that the pathway has l blocks in series and each block bi has ai parallel lines with each line j containing bji targets. the total number of targets in the general map is ∑i=1l∑j=1aibji. assuming that the n targets are distinct, the maximum number of distinct discrete dynamic models satisfying the structure is l!∏i=1l∏j=1ai!. each parallel line j in block i can have bji! possible directional orientations. if the figure represents a possible directional orientation, then only the targets k,k ⋯ka will have initial activations due to mutations or latent activations. some other downstream target cannot have a mutation or latent activation otherwise the target inhibition combination k,k ⋯ka will not be effective. for our analysis, we are assuming that we can inhibit specific targets of our choice and we can measure the steady state target expression following application of the target inhibitions. we can locate the directionality of the blocks b to bl by using at most l− steady state measurements. we can start by randomly picking any block bi and blocking the targets in that block, the blocks that will remain activated will be upstream of that block and the blocks that will be deactivated following the inhibition of block bi will be located down-stream of bi. note that the number of experiments required is based on steady state measurements following particular perturbations. time series measurements can reduce the number of experiments required but may not be always technically feasible. the expected number of experiments required to detect the directionality of l serial blocks is 2l−13forl≥ the next step will be locating the directionality of targets in each parallel line of the block. we can start with an experiment where for each block bi, one target from each line up to a maximum of ai− lines will be inhibited. we cannot inhibit all the lines in a block or else the downstream blocks will also be inhibited and no inference can be made on those blocks for that experiment. while locating the directionality of the serial blocks bi, we have already validated the position of one target from each parallel line in a serial block. if we consider a single block bi, each experiment can detect the location of ai− targets, thus the total number of experiments required to decipher the possible directionalities of the targets in the block bi is ≤max where si={ ⋯,ai}. thus for the overall map, the worst case number of experiments new required to decipher the directionalities of all the targets is upper-bounded by new≤maxi∈s{max}+l− where s={ ⋯,l}. utilizing equation the expected number of experiments nea required to decipher the directionalities of all the targets is upper-bounded by nea≤maxi∈s{max⌉−˘1)}+2l−  in this article, we presented a novel framework for predicting the effectiveness of molecularly targeted drugs. we used drug perturbation data to generate a map of the underlying genetic regulatory pathway. using actual experimental data, we were able to show the effectiveness of our approach for drug sensitivity prediction. the proposed tim approach produced a low average leave one out cross validation error of 5% when applied to perturbation data generated from four primary canine tumors using a set of drugs. we should note that the current drug screen is a small one and technology has been developed for drug screens with a far greater number of drugs. we are currently experimenting with pharmaceutical drug library consisting of more than small molecule inhibitors. we expect that the use of larger number of drugs will increase the accuracy further and generate maps with greater robustness. the scope of the present article is concentrated around steps b, c and d of figure for future research, we will consider multiple data sources to increase the robustness of the designed maps. as explained in figure we can use rapid sirna screens to validate single points of failures predicted by our tim approach. furthermore, rnaseq and protein phosphoarray data can be used to further revise the circuit. finally, time series data can be used to incorporate dynamics in the modeling framework. for combination therapy design, we can use the tim framework to formulate control strategies with various constraints. some possibilities are minimal toxicity, anticipating evolving drug resistance, and success over a family of tims representing variations of a tumor. for case, we can assume that the toxicity of a drug or drug combination is proportional to the number of targets being inhibited by the drug and search for the drug combination with high sensitivity but low set of target inhibitions. for case, we would want to avoid resistance and thus would like to inhibit more than one independent blocking pathway such that for the scenario when resistance to one of the blocking pathways develops, the other independent pathway can still keep the tumor under check. in other words, we would be interested in selecting a set of targets that can be divided into two or more non-intersecting sets such that the sensitivity of each set is higher than a threshold. for case, the goal is to design control policies for the scenario when the exact pathway is not known but it belongs to a collection of pathways. the uncertainty can arise when the experimental data is not sufficient enough to produce a unique pathway map or the current pathway may evolve into one of the different pathways obtained from tissues with same type of cancer. this can approached from a worst case perspective or a bayesian perspective. in conclusion, the proposed framework provides a unique input-output based methodology to model a cancer pathway and predict the effectiveness of targeted drugs. this framework can be developed as a viable approach for personalized cancer therapy. to aide in the usage of our framework, we have developed a graphical user interface which implements in an easy to use way the algorithms and equations presented in this paper. it is built in matlab, but is distributed as a compiled executable; as such, it is usable in a windows environment by downloading the matlab compile runtime environment, which is free to download and requires no matlab installation. it is available online at: http://cvial.ece.ttu.edu/ranadippal/research.html under the target inhibition map approach to inference of cancer pathways heading. competing interests author ck has received honoraria for scientific presentations at novartis, millennium/takeda pharmaceutical and glaxosmithkline, and has research joint ventures or sponsored research with novartis and johnson & johnson. ck is also a paid consultant to the nci/ctep pediatric preclinical testing program. ck is also a co-founder/stockholder of numira biosciences. the remaining authors have no other conflicts to declare related to these studies except that a ec is a first degree relative of a bms scientist. authors’ contributions designed the algorithmic framework: nb, rp. implemented the algorithms: nb. performed the experiments: ld, ec, bs, ck. analyzed the results: nb, rp. all authors read and approved the final manuscript. supplementary material additional file 1 supplementary analysis. the file includes a detailed tim example, results of a set of synthetic experiments, theoretical analysis of errors and additional tim circuits. click here for file additional file 2 drug screenic dataset. the experimentally measured ic50’s in nm for four cell cultures across drugs. click here for file additional file 3 drug screen panelec the ec values in nm for targets for drugs. click here for file additional file 4 cross validation results. the detailed results of the fold cross validation error estimates for the tumor cultures bailey, charley, sy and cora. click here for file additional file 5 similarity between drugs. the similarities between the targeted drugs as measured by similarity measure Δ. click here for file acknowledgements this work was supported in part by nsf grant ccf and nci grant 5r01ca133229-05/ 
the proposed framework provides a unique input-output based methodology to model a cancer pathway and predict the effectiveness of targeted anti-cancer drugs. we illustrate the high prediction accuracy of our framework on synthetic data generated from the kyoto encyclopedia of genes and genomes and an experimental dataset of four canine osteosarcoma tumor cultures following application of targeted small-molecule drugs. we achieve a low leave one out cross validation error of <10% for the canine osteosarcoma tumor cultures using a drug screen consisting of targeted drugs. this framework can be developed as a viable approach for personalized cancer therapy. the prediction of tumor sensitivity to targeted drugs remains a major challenge in the design of optimal therapeutic strategies. the current sensitivity prediction approaches are primarily based on genetic characterizations of the tumor sample.
 the success of targeted anti-cancer drugs are frequently hindered by the lack of knowledge of the individual pathway of the patient and the extreme data requirements on the estimation of the personalized genetic network of the patient’s tumor. the prediction of tumor sensitivity to targeted drugs remains a major challenge in the design of optimal therapeutic strategies. the current sensitivity prediction approaches are primarily based on genetic characterizations of the tumor sample. we propose a novel sensitivity prediction approach based on functional perturbation data that incorporates the drug protein interaction information and sensitivities to a training set of drugs with known targets.  we illustrate the high prediction accuracy of our framework on synthetic data generated from the kyoto encyclopedia of genes and genomes and an experimental dataset of four canine osteosarcoma tumor cultures following application of targeted small-molecule drugs. we achieve a low leave one out cross validation error of <10% for the canine osteosarcoma tumor cultures using a drug screen consisting of targeted drugs.  the proposed framework provides a unique input-output based methodology to model a cancer pathway and predict the effectiveness of targeted anti-cancer drugs. this framework can be developed as a viable approach for personalized cancer therapy. 