Biphasic Mathematical Model of Cell–Drug Interaction That Separates Target-Specific and Off-Target Inhibition and Suggests Potent Targeted Drug Combinations for Multi-Driver Colorectal Cancer Cells

Quantifying the response of cancer cells to a drug, and understanding the mechanistic basis of the response, are the cornerstones for anti-cancer drug discovery. Classical single target-based IC50 measurements are inadequate at describing cancer cell responses to targeted drugs. In this study, based on an analysis of targeted inhibition of colorectal cancer cell lines, we develop a new biphasic mathematical model that accurately describes the cell–drug response. The model describes the drug response using three kinetic parameters: ratio of target-specific inhibition, F1, potency of target-specific inhibition, Kd1, and potency of off-target toxicity, Kd2. Determination of these kinetic parameters also provides a mechanistic basis for predicting effective combination targeted therapy for multi-driver cancer cells. The experiments confirmed that a combination of inhibitors, each blocking a driver pathway and having a distinct target-specific effect, resulted in a potent and synergistic blockade of cell viability, improving potency over mono-agent treatment by one to two orders of magnitude. We further demonstrate that mono-driver cancer cells represent a special scenario in which F1 becomes nearly 100%, and the drug response becomes monophasic. Application of this model to the responses of >400 cell lines to kinase inhibitor dasatinib revealed that the ratio of biphasic versus monophasic responses is about 4:1. This study develops a new mathematical model of quantifying cancer cell response to targeted therapy, and suggests a new framework for developing rational combination targeted therapy for colorectal and other multi-driver cancers.


Introduction
Some cancers rely on a single proliferative driver and its associated signaling pathway. Abl in chronic myeloid leukemia (CML) [1], ErbB2 in some breast cancers [2], and EGFR in some non-small cell lung cancer [3] are a few examples of such mono-driver cancers. Mono-driver cancers can be effectively treated by targeted therapy blocking the function of the proliferative drivers. Small molecule inhibitors or monoclonal antibodies blocking Abl, ErbB2, or EGFR have become the standard of care for these cancers [1][2][3]. Table 1. Inhibition of HT-29 viability by protein kinase inhibitors. Means ± SEM are reported for six sets of data from two independent experiments in triplicate. The table only lists the most common kinase targets for each inhibitor. Most of the inhibitors have been tested against most of the kinome, and complete inhibitory information is available at the Library of Integrated Network-based Cellular Signatures (http://lincs.hms.harvard.edu/kinomescan/).

Drug
Target IC 50  HG6-64-1 is a specific BRAF inhibitor [33] and AZD-6244 (selumetinib) is a Mek kinase inhibitor [34]. The HT-29 genome contains two BRAF mutations, V600E and T119S [26]. The V600E mutation activates BRAF, which phosphorylates and activates Mek. The inhibition of HT-29 viability by these two inhibitors suggests that the BRAF-activated MAP kinase pathway is a driver for HT-29 proliferation and/or viability. Dasatinib is a broad inhibitor for kinases in the Abl, Src, and PDGFR families [35]. Due to the fact that the cells did not respond to the Abl inhibitor imatinib, or the PDGFR inhibitor sunitinib, and Src is over-expressed in HT-29 cells with a Z score of 2.99 [26] and activated [27], the most likely dasatinib target is Src. BMS-754807 is an inhibitor for the insulin receptor (IR) and the insulin-like growth factor 1 receptor (IGF-1R) [36]. These results suggest that BRAF/MAP kinases, Src kinases and IR/IGF-1R kinases all contribute to the viability of HT-29 cells.

IC 50 Values Do Not Accurately Reflect How HT-29 Cells Respond to Targeted Therapy
Traditionally, cell response to a drug is measured by IC 50 , the concentration that inhibits cell viability by 50%, which is commonly calculated by fitting the data into a single target inhibition  Figure S1 were the best fitting single target kinetic curves for HG6-64-1, BMS-754807, dasatinib and AZD-6244. The poor fitting is evident judged visually or by the root mean square error (RMSE), and is clearly not the result of experimental inaccuracy. The kinetically best fitting IC 50 values were often significantly different from the visually determined IC 50 values reported in Table 1.
In an attempt to better understand the dose responses, the data were also fitted to a modified Hill equation: I = I max × D n /(IC 50 * n + D n ). In this equation, the relative inhibition (I) is expressed as a function of drug concentration (D) determined by three parameters: I max , IC 50 * and n. I max is the maximal relative inhibition a drug can achieve at saturating concentration, IC 50 * is the drug concentration that achieves half of maximal inhibition, and n is the Hill coefficient, or slope. We did not use the fourth parameter of the full Hill equation, I 0 (inhibition in the absence of an inhibitor), because there is no inhibition at the drug concentration of zero. Fitting the data into the above equation produced the graphs shown in Figure S1.
The equation fit the data reasonably well, producing IC 50 * values that reflected how potently the drugs inhibited cell viability. However, each drug generated a highly variable n value significantly less than 1 (ranging from 0.34 for BMS-754807 to 0.84 for HG6-64-1), suggesting varying levels of negative cooperativity. As there is no obvious explanation for such apparent negative cooperativity, it is presumed to reflect some undefined biological property of the system. Such a mathematically accommodating, but biologically undefined, parameter is considered one of the limitations of the Hill equation when it is applied to drug dose response curves.
Furthermore, three of the four drugs also produced an I max that is significantly less than 1, suggesting that even when a drug saturates its target, it still does not fully inhibit cell viability. This indicated that the target kinase of the drug is only partially responsible for cell viability. We further noticed that some drugs, such as dasatinib and BMS-754807, displayed a biphasic inhibition pattern, with a nM inhibitory phase and µM inhibitory phase. The nM inhibition was assumed to correlate to the target-specific inhibition, while the µM inhibition likely reflected an off-target effect. Based on these considerations, we decided to further modify the Hill equation to set n at 1, and add a second phase inhibition. Thus the equation was modified into a biphasic equation: . In this equation, the inhibition (I) has two phases: F 1 and F 2 as fractions of total cell viability with their respective K d s (K d1 and K d2 ). The biphasic assumption dictates that F 1 + F 2 = 1 (or 100%), thus F 2 is not an independent variable, but a derivative of F 1 . The biphasic curve fitting graphs are shown in Figure 1A-D. In every case, the biphasic equation fits the data well, judged visually or by the RMSE. The RMSE values for the biphasic equation fitting were often 5-10 times lower than the monophasic equation fitting, and similar to those of the full Hill equation fitting. The biphasic kinetic curve fitting revealed the best fitting F1, F2, Kd1 and Kd2 for each drug. The inhibition by HG6-64-1 had an F1 of 57% with a Kd1 of 19.8 nM, and an F2 of 43% with Kd2 > 1000 µM (largely resistant). This analysis indicated that HG6-64-1 inhibited 57% of cell viability with a Kd of 19.8 nM, but left the remaining portion of viability unaffected. Similarly, the responses to the other drugs also consisted of an F1 with a low nM Kd1 and an F2 with a >10 µM Kd2. These analyses suggest that the cell responses can be accurately described by three constants: F1/F2 ratio, Kd1 and Kd2 in the biphasic model.

Biological Meaning of Biphasic Inhibition
Two plausible scenarios are consistent with the biphasic inhibition patterns. One is that HT-29 cells contain two subpopulations, one being more sensitive (F1 and Kd1) than the other (F2 and Kd2). If this scenario is true, treatment with a drug should eliminate/reduce the subpopulation sensitive to the drug. To test this scenario, the HT-29 cells were treated with 2 µM dasatinib, BMS-754807, or no drug for three days, and the residual cells were recovered and retested for their sensitivity toward the two drugs. The recovered cells from the three treatments displayed virtually identical sensitivity to either dasatinib or BMS-754807 (Figure 2A,B). This result argues against the two-population possibility.
The second scenario is that each drug causes the F1 inhibition by blocking a high affinity driver and causes the F2 inhibition through a non-specific off-target mechanism. This scenario is consistent with two basic features of this experimental system: (1) the kinase inhibitors used inhibit their targets with low nM Kds, but also inhibit other protein kinases non-specifically at µM concentrations; and (2) HT-29 cells are sensitive to multiple kinase inhibitors, suggestive of multiple independent drivers. The combination of these two features makes it possible that each drug partially blocks viability when it inhibits a high affinity target, and inhibits the remaining portion of viability when the drug concentration reaches a cytotoxic level. It is well established that most kinase inhibitors cause cytotoxicity at µM concentrations. Elimination of the first scenario leaves the second as the most plausible explanation.
. Each set of inhibition data represents two independent experiments in triplicate, with standard errors presented as error bars. The inhibitory parameters obtained from curve fitting are presented in the graphs.
The biphasic kinetic curve fitting revealed the best fitting F 1 , F 2 , K d1 and K d2 for each drug. The inhibition by HG6-64-1 had an F 1 of 57% with a K d1 of 19.8 nM, and an F 2 of 43% with K d2 > 1000 µM (largely resistant). This analysis indicated that HG6-64-1 inhibited 57% of cell viability with a K d of 19.8 nM, but left the remaining portion of viability unaffected. Similarly, the responses to the other drugs also consisted of an F 1 with a low nM K d1 and an F 2 with a >10 µM K d2 . These analyses suggest that the cell responses can be accurately described by three constants: F 1 /F 2 ratio, K d1 and K d2 in the biphasic model.

Biological Meaning of Biphasic Inhibition
Two plausible scenarios are consistent with the biphasic inhibition patterns. One is that HT-29 cells contain two subpopulations, one being more sensitive (F 1 and K d1 ) than the other (F 2 and K d2 ). If this scenario is true, treatment with a drug should eliminate/reduce the subpopulation sensitive to the drug. To test this scenario, the HT-29 cells were treated with 2 µM dasatinib, BMS-754807, or no drug for three days, and the residual cells were recovered and retested for their sensitivity toward the two drugs. The recovered cells from the three treatments displayed virtually identical sensitivity to either dasatinib or BMS-754807 (Figure 2A,B). This result argues against the two-population possibility.
The second scenario is that each drug causes the F 1 inhibition by blocking a high affinity driver and causes the F 2 inhibition through a non-specific off-target mechanism. This scenario is consistent with two basic features of this experimental system: (1) the kinase inhibitors used inhibit their targets with low nM K d s, but also inhibit other protein kinases non-specifically at µM concentrations; and (2) HT-29 cells are sensitive to multiple kinase inhibitors, suggestive of multiple independent drivers. The combination of these two features makes it possible that each drug partially blocks viability when it inhibits a high affinity target, and inhibits the remaining portion of viability when the drug concentration reaches a cytotoxic level. It is well established that most kinase inhibitors cause cytotoxicity at µM concentrations. Elimination of the first scenario leaves the second as the most plausible explanation.

Src, IR/IGF-1R/AKT Signaling and the MAP Kinase Pathway Are All Partially Responsible for Driving HT-29 Cell Proliferation and Viability
To verify that AZD-6244, BMS-754807, dasatinib, and HG6-64-1 are indeed blocking distinct drivers, we determined the effects of each drug on the activation status of key signaling proteins, including insulin receptor (IR), IGF-1R, PDK1, AKT, AKT substrates PRAS40 and GSK3β, BRAF, Mek, Erk and Src ( Figure 2C). For this analysis, we also included an AKT inhibitor, GSK690693 [37], even though GSK690693 surprisingly did not inhibit HT-29 viability ( Table 1).
The immunoblots in Figure 2C revealed the following. (1) AZD-6244 (Mek inhibitor) strongly inhibited the phosphorylation of Erk, a Mek substrate. It may have also inhibited, to a lesser degree, the activation of AKT. (2) Western blots of the phosphorylated IR and IGF-1R were not informative, due to the detection of a large number of protein bands smaller than the expected sizes of these two receptors ( Figure S4). Despite repeated attempts under different conditions in consultation with the antibody manufacturer, these lower bands, either degradation products or nonspecific recognition, persisted. However, BMS-754807 (inhibitor of IR and IGF-1R) clearly inhibited Akt phosphorylation on both T308 and S473, an indication of Akt activation [38]. The blockade of Akt activation was confirmed by the undetectable phosphorylation of PRAS40 and reduced phosphorylation level of GSK3β, two direct substrates of Akt [38]. In addition, the phosphorylation level of BRAF, MEK and ERK also decreased in response to BMS-754807. These results indicated that blocking IR/IGF-1R by BMS-754807 resulted primarily in the down-regulation of AKT, but also the MAP kinase pathway.
The immunoblots in Figure 2C revealed the following. (1) AZD-6244 (Mek inhibitor) strongly inhibited the phosphorylation of Erk, a Mek substrate. It may have also inhibited, to a lesser degree, the activation of AKT. (2) Western blots of the phosphorylated IR and IGF-1R were not informative, due to the detection of a large number of protein bands smaller than the expected sizes of these two receptors ( Figure S4). Despite repeated attempts under different conditions in consultation with the antibody manufacturer, these lower bands, either degradation products or nonspecific recognition, persisted. However, BMS-754807 (inhibitor of IR and IGF-1R) clearly inhibited Akt phosphorylation on both T308 and S473, an indication of Akt activation [38]. The blockade of Akt activation was confirmed by the undetectable phosphorylation of PRAS40 and reduced phosphorylation level of GSK3β, two direct substrates of Akt [38]. In addition, the phosphorylation level of BRAF, MEK and ERK also decreased in response to BMS-754807. These results indicated that blocking IR/IGF-1R by BMS-754807 resulted primarily in the down-regulation of AKT, but also the MAP kinase pathway. (3) Dasatinib completely blocked Src phosphorylation on Y527 by another PTK, Csk [39], and significantly inhibited Src autophosphorylation on Y416. This is consistent with earlier reports that both Src and Csk are sensitive to dasatinib [40]. Dasatinib also appeared to partially inhibit Mek and Erk activation, suggestive of crosstalk between Src and the MAP kinase pathway. (4) HG6-64-1 [33] (a BRAF inhibitor) moderately inhibited BRAF autophosphorylation, but completely blocked Mek phosphorylation by BRAF, and reduced Erk phosphorylation by Mek. It also reduced the phosphorylation of AKT and PRAS40. These results suggested considerable cross talk between the MAP kinase, IR/IGF-1R/Akt and Src pathways.
GSK690693 is a specific Akt inhibitor [37]. Although it did not inhibit HT-29 cell viability (Table 1), it completely blocked the phosphorylation of PRAS40 and strongly inhibited phosphorylation of GSK3β ( Figure 2C), confirming its ability to block Akt kinase function. However, it also strongly activated Akt phosphorylation on S473, as previously reported [37,41,42]. Phosphorylation of Akt on S473 is catalyzed by PDK2, a kinase activity attributed to several different kinases [43], especially rapamycin complex 2 (mTORC2) [44]. Activation of Akt S473 phosphorylation by GSK690693 was thought to be a feedback mechanism by mTORC2. mTORC2 also phosphorylates other AGC kinases at their corresponding hydrophobic motif residues [43]. The activation of PDK2/mTORC2 by GSK690693 may explain its lack of inhibition of HT-29 viability.
The above results overall demonstrated that AZD-6244, BMS-754807, dasatinib, and HG6-64-1 indeed blocked their respective target kinases, and caused distinct patterns of signaling inhibition. These results further support the multi-driver proliferation hypothesis in HT-29 cells.

HT-29 Cell Viability Can Be Effectively and Synergistically Blocked by Combinations of Drugs, Each Blocking an Independent Driver
If the multi-driver proliferation mechanism is true, then the combination of inhibitors with non-overlapping F 1 targets should synergistically block the proliferation. The effects of pair-wise drug combinations (1:1 ratio) on HT-29 viability were determined. In Figure 3A, the IC 50 s for dasatinib and BMS-754807 were 500 nM and 790 nM, respectively, and the IC 60 s for either drug was above 2 µM. The IC 50 and IC 60 for the drug combination were 28 nM and 50 nM, respectively. The dose reduction index (DRI), a measure of how many folds the drug doses may be reduced by the combination to achieve a given inhibition level [45], was 11-fold at 50% inhibition and 31-fold at 60% inhibition ( Figure 3A), meaning that the combination is 31 times more effective than the two drugs alone in achieving 60% inhibition. The combinations of dasatinib with AZD-6244 or HG6-64-1, and the combination of BMS-754807 with AZD-6244, also showed dramatic synergy ( Figure 3B-D). These results suggest that the drugs in these pairs (dasatinib/BMS-754807, dasatinib/HG6-64-1, dasatinib/AZD-6244, and BMS-754807/AZD-6244) exert their effects on HT-29 cells largely independently, which is consistent with the multi-driver proliferation hypothesis.
We also tested the combination of HG6-64-1 and AZD-6244, inhibitors of two consecutive steps in the MAP kinase pathway. As their F 1 effects on HT-29 cell viability would largely overlap, their combination would be redundant rather than complementary/synergistic. The combination indeed produced very little synergistic benefit ( Figure 3E), still retaining approximately 40% viability when both drugs were at 20 µM. Thus, the portion of viability that is resistant to HG6-64-1 is also largely resistant to AZD-6244, and independent of the MAP kinase pathway.
Two drugs, BGJ398 (FGFR inhibitor) [46] and crizotinib (Met inhibitor) [47] did not inhibit HT-29 viability in the nM range, but inhibited HT-29 viability at around 10 µM (Table 1, Figure 3F), likely through off-target toxicity. The combination of these two drugs did not result in synergistic inhibition ( Figure 3F), indicating that only target-specific effects by drugs on different targets would have synergistic effects.  The combination of three drugs, dasatinib, BMS-754807 and AZD-6244, had an IC50 of 30 ± 5 nM on HT-29 cells, marginally lower than the IC50s of the pair-wise combinations, suggesting that blocking any two of the three pathways is largely sufficient in blocking HT-29 viability. This overlap between the effects of these drugs is consistent with the cross-talk observed in the Western blots ( Figure 2C).

Other CRC Cancer Cell Lines Follow a Similar Pattern in Their Responses to Individual Kinase Inhibitors and Inhibitor Combinations
We previously reported [31] that the CRC cell lines SK-CO-1 and NCI-H747 were mildly inhibited by both BMS-754807 and AZD-6244 and much more potently by the two inhibitors combined. We tested if the responses of these two cells to BMS-754807 and AZD-6244 were also biphasic ( Figure 4). The responses of SK-CO-1 to BMS-754807 and AZD-6244 clearly did not fit the monophasic equation ( Figure 4A,B) (dashed lines), indicated by a poor fitting curve with high RMSEs. Biphasic analyses of the same data ( Figure 4A,B, solid lines) yielded much better fitting curves. The two drugs had F1 values of 44% and 55% respectively, and Kd1s 32 nM and 168 nM, respectively. This suggested that the targets of BMS-754807 and AZD-6244 are respectively responsible for 44% and 55% of SK-CO-1 viability. The results suggested that the MAP kinase The combination of three drugs, dasatinib, BMS-754807 and AZD-6244, had an IC 50 of 30 ± 5 nM on HT-29 cells, marginally lower than the IC 50 s of the pair-wise combinations, suggesting that blocking any two of the three pathways is largely sufficient in blocking HT-29 viability. This overlap between the effects of these drugs is consistent with the cross-talk observed in the Western blots ( Figure 2C).

Other CRC Cancer Cell Lines Follow a Similar Pattern in Their Responses to Individual Kinase Inhibitors and Inhibitor Combinations
We previously reported [31] that the CRC cell lines SK-CO-1 and NCI-H747 were mildly inhibited by both BMS-754807 and AZD-6244 and much more potently by the two inhibitors combined. We tested if the responses of these two cells to BMS-754807 and AZD-6244 were also biphasic ( Figure 4). The responses of SK-CO-1 to BMS-754807 and AZD-6244 clearly did not fit the monophasic equation ( Figure 4A,B) (dashed lines), indicated by a poor fitting curve with high RMSEs. Biphasic analyses of the same data ( Figure 4A,B, solid lines) yielded much better fitting curves. The two drugs had F 1 values of 44% and 55% respectively, and K d1 s 32 nM and 168 nM, respectively. This suggested that the targets of BMS-754807 and AZD-6244 are respectively responsible for 44% and 55% of SK-CO-1 viability. The results suggested that the MAP kinase pathway and the BMS-754807-sensitive pathway, the IR/IGF-1R initiated Akt signaling pathway, each counted for approximately half of SK-CO-1 viability. pathway and the BMS-754807-sensitive pathway, the IR/IGF-1R initiated Akt signaling pathway, each counted for approximately half of SK-CO-1 viability. As shown in 4C, SK-CO-1 cells were much more potently inhibited by the combination of AZD-6244 and BMS-754807 than either drug alone. The IC50 and IC70 were more than an order of magnitude lower than the corresponding values for either drug alone. The synergy was especially pronounced at higher inhibition levels. The IC70 values were >7 µM for either drug alone, but only 261 nM for the drug combination. Very similar biphasic inhibitory patterns and synergistic inhibition were observed for NCI-H747 cells ( Figure 4D-F). These results suggest that multi-driver proliferation is a common mechanism in these cells and combination targeted therapy may be broadly effective for CRC. The synergy and potency of double and triple drug combinations for three CRC cell lines are summarized in Figure 5A,B. As shown in 4C, SK-CO-1 cells were much more potently inhibited by the combination of AZD-6244 and BMS-754807 than either drug alone. The IC 50 and IC 70 were more than an order of magnitude lower than the corresponding values for either drug alone. The synergy was especially pronounced at higher inhibition levels. The IC 70 values were >7 µM for either drug alone, but only 261 nM for the drug combination. Very similar biphasic inhibitory patterns and synergistic inhibition were observed for NCI-H747 cells ( Figure 4D-F). These results suggest that multi-driver proliferation is a common mechanism in these cells and combination targeted therapy may be broadly effective for CRC. The synergy and potency of double and triple drug combinations for three CRC cell lines are summarized in Figure 5A,B. Cancers 2020, 12, x FOR PEER REVIEW 10 of 22

Mono-Driver Cancer Cell Responses to Targeted Therapy Is Monophasic
To validate that the biphasic responses and combination synergy are unique to multi-driver cancer cells, we determined how a mono-driver cancer system would behave in these analyses. Targeted therapy using EGFR inhibitors has become the standard of care for non-small cell lung cancers that contain activating EGFR mutations [3], even though not all patients with EGFR mutations respond to EGFR inhibitors, such as erlotinib [48]. HCC-827 is a model cell line for such cancers [49]. It contains four mutations that activate EGFR, which drives the proliferation [26]. HCC-827 cells are indeed very sensitive to several EGFR inhibitors, gefitinib, afatinib and erlotinib, with IC50s between 7.8 nM and 15.4 nM ( Table 2). They are also inhibited by EGFR inhibitors neratinib and lapatinib, although with higher IC50s. The cells are also sensitive to dasatinib and bosutinib with IC50s of 140 and 273 nM, respectively. Mek >20,000 BX-912 PDK1 >20,000

Mono-Driver Cancer Cell Responses to Targeted Therapy Is Monophasic
To validate that the biphasic responses and combination synergy are unique to multi-driver cancer cells, we determined how a mono-driver cancer system would behave in these analyses. Targeted therapy using EGFR inhibitors has become the standard of care for non-small cell lung cancers that contain activating EGFR mutations [3], even though not all patients with EGFR mutations respond to EGFR inhibitors, such as erlotinib [48]. HCC-827 is a model cell line for such cancers [49]. It contains four mutations that activate EGFR, which drives the proliferation [26]. HCC-827 cells are indeed very sensitive to several EGFR inhibitors, gefitinib, afatinib and erlotinib, with IC 50 s between 7.8 nM and 15.4 nM ( Table 2). They are also inhibited by EGFR inhibitors neratinib and lapatinib, although with higher IC 50 s. The cells are also sensitive to dasatinib and bosutinib with IC 50 s of 140 and 273 nM, respectively. In contrast to the CRC cells, the HCC-827 responses to EGFR and Src inhibitors followed single-target equation. For example, the gefitinib inhibition data fit well to the monophasic equation, with an IC 50 of 14.1 nM ( Figure 6A, dashed line). Fitting the same data to the biphasic equation resulted in a very similar graph with an F 1 of 92%, and a K d1 of 11.2 nM ( Figure 6A, solid line), indicating that the target-specific F 1 response is predominantly responsible for gefitinib inhibition of HCC-827. Similar results were obtained for erlotinib inhibition (IC 50 = 19.6 nM, F 1 = 93%, K d1 = 16.1 nM) ( Figure 6B) and dasatinib (IC 50 = 165 nM, F 1 = 100%, K d1 = 165 nM) ( Figure 6C). These results indicated that any of these drugs inhibiting its specific target (EGFR or Src) alone was largely sufficient to block HCC-827 viability. The fact that inhibiting either EGFR or Src would be sufficient to block the viability indicates that EGFR and Src are dependent on each other in supporting the viability of this cell line. Therefore, EGFR and Src are not independent drivers for HCC-827 proliferation and viability.
In contrast to the CRC cells, the HCC-827 responses to EGFR and Src inhibitors followed singletarget equation. For example, the gefitinib inhibition data fit well to the monophasic equation, with an IC50 of 14.1 nM ( Figure 6A, dashed line). Fitting the same data to the biphasic equation resulted in a very similar graph with an F1 of 92%, and a Kd1 of 11.2 nM ( Figure 6A, solid line), indicating that the target-specific F1 response is predominantly responsible for gefitinib inhibition of HCC-827. Similar results were obtained for erlotinib inhibition (IC50 = 19.6 nM, F1 = 93%, Kd1 = 16.1 nM) ( Figure  6B) and dasatinib (IC50 = 165 nM, F1 = 100%, Kd1 = 165 nM) ( Figure 6C). These results indicated that any of these drugs inhibiting its specific target (EGFR or Src) alone was largely sufficient to block HCC-827 viability. The fact that inhibiting either EGFR or Src would be sufficient to block the viability indicates that EGFR and Src are dependent on each other in supporting the viability of this cell line. Therefore, EGFR and Src are not independent drivers for HCC-827 proliferation and viability.  At µM concentrations, dasatinib caused a slight but consistent uptick in HCC-827 cell viability ( Figure 6C). This is likely a reflection of the complex effects of dasatinib on cell signaling and biochemistry. As a broad spectrum kinase inhibitor, it is conceivable that it could directly or indirectly stimulate cell viability by inhibiting some kinases. It has been previously reported that dasatinib induced elevated levels of phosphorylated AKT protein as well as ERK1/2 proteins in FaDu and A431 cells [50], and the activation of these growth stimulating kinases could result in the stimulation of cell viability. Whether or not the apparent stimulation of HCC-827 cells is due to this mechanism is yet to be determined.
The synergy between Src and EGFR inhibitors on HCC-827 was also tested ( Figure 6D,E). Neither the neratinib/dasatinib combination nor the gefitinib/bosutinib combination displayed significant synergy, indicating that blocking EGFR, Src or both achieved largely the same effect on HCC-827 viability, further demonstrating that Src and EGFR are not independent drivers in HCC-827.
A Western blot analysis demonstrated that ( Figure 6F) gefitinib blocked EGFR phosphorylation and activation, and the EGFR downstream pathways, such as the phosphorylation of Mek, Erk, and Akt substrate PRAS40. The gefitnib treatment did not affect the Src phosphorylation on Tyr416 or Tyr527. Dasatinib not only blocked Src phosphorylation on Tyr416 and its activation, but also blocked EGFR phosphorylation, Mer, Erk and PRAS40 phosphorylation. Therefore, blocking Src function also prevented EGFR signaling, indicating that the function of EGFR in activating the MAP kinase pathway and the Akt pathway is dependent on Src function. This finding is consistent with reports that EGFR activation is dependent on recruiting Src to phosphorylate several phosphorylation sites on EGFR [51,52]. These results suggest that mutation-activated EGFR is the true driver of HCC-827 proliferation, while Src is an essential component in the EGFR-initiated signaling pathways.

Inhibition of Cancer Cell Viability by Dasatinib Is Mostly Biphasic
Even though the biphasic inhibition pattern is strongly supported in our hands, such an inhibitory pattern has not been reported in the literature to our knowledge. The only report [53] that comes close to it is that 28% of cancer cell-drug response curves were better described by a multi-phasic model. To determine if the biphasic inhibition model is broadly applicable, we applied the analysis to dasatinib dose response data that are available in the Genomics of Drug Sensitivity in Cancer database [54].
In this dataset, the responses of 426 cancer cell lines to dasatinib at nine concentrations from 5 nM to 1.28 µM (2× dilution series) were determined. The data were collected through a high throughput screen, with each cell line/drug concentration tested in one well in the 384-well format [54]. We removed dose responses that contained obvious scattering, and the remaining 362 datasets were analyzed. The responses to dasatinib were divided into four categories ( Figure 7A): Category LV (low viability) contains 13 cell lines characterized by low viability readings regardless of the drug concentration ( Figure 7B). All 13 cell lines are suspension blood cancer cell lines. The second category (category NI, no inhibition) contains 175 cells, and they are not inhibited by dasatinib up to 1.28 µM. Due to the lack of concentration-dependent inhibition, both of these categories were not analyzed further for this study.
The remaining 174 cell lines responded to dasatinib in dose-dependent manners with varying levels of sensitivity. Both monophasic and biphasic analyses were performed on these cell lines (Data S1). In the biphasic analysis, the F 1 value varied from 100% to about 10%. Thirty-four (34/174) cell lines had an F 1 above 85%, and these cells were inhibited more than~80% at 1.28 µM dasatinib ( Figure 7C). Their responses fit the monophasic equation and the biphasic equation similarly well, with comparable RMSE values (0.086 + 0.004 for MP analysis versus 0.074 + 0.004 for BP analysis). Kinetic analyses of the ETK-1 cell line are shown as an example of this group in Figure 7D. In this group, the K d1 values closely aligned to the IC 50 values from the monophasic analysis ( Figure 7E). This group is called the MP category, where inhibiting the dasatinib target is largely sufficient to block cell viability. Cells in this category are likely fully or nearly fully dependent on a proliferative driver that is sensitive to dasatinib, such as abl, Src of PDGFR family kinases.  Figure 7G,H. As a group, the IC50 values for these cells range from ~100 nM to ~4 µ M, and the increase in IC50 values inversely correlates to the decrease in F1 ( Figure 7I). Because the highest dasatinib concentration used in the screen was 1.28 µ M, too low to cause significant off-target toxicity, the data do not yield a reliable Kd2. It is also interesting to note that the Kd1 for this group mostly remained below 100 nM, even though the IC50 are up to ~4 µ M ( Figure 7J), indicating that the affinity The remaining 140 cell lines (BP category) display clear biphasic characteristics with F1 < 85%. They all have shallow response curves, displaying clear inhibition at low dasatinib concentrations, but still retain considerable cell viability even at 1.28 µM ( Figure 7F). These response datasets fit the biphasic equation much better than the monophasic equation (RMSE = 0.142 ± 0.004 for MP analysis, and 0.038 ± 0.002 for BP analysis). Two examples of the monophasic and biphasic fitting comparisons are shown in Figure 7G,H. As a group, the IC 50 values for these cells range from~100 nM to~4 µM, and the increase in IC 50 values inversely correlates to the decrease in F 1 ( Figure 7I). Because the highest dasatinib concentration used in the screen was 1.28 µM, too low to cause significant off-target toxicity, the data do not yield a reliable K d2 . It is also interesting to note that the K d1 for this group mostly remained below 100 nM, even though the IC 50 are up to~4 µM ( Figure 7J), indicating that the affinity between dasatinib and its target is largely independent of the cell host. The ratio of monophasic versus biphasic responses is 34:140, suggesting that most cancer cells are likely to be multi-driver dependent. This analysis indicates that even though the biphasic response pattern to targeted therapy has not been widely appreciated, it is actually very common.

Discussion
Most cancers rely on multiple drivers for full proliferation, and require drug combinations for effective targeted therapy. Currently, most drug combinations are discovered by empirical testing [18,19], a rational and a mechanism-based platform would greatly aid in the formulation of effective drug combinations. In the current study, we investigated how multi-driver cancer cells respond to targeted therapy, and developed a mathematical model that quantitatively describes the response of multi-driver cancer cells to targeted therapy. This analysis suggested that multi-driver cancer cells have a biphasic response to targeted agents: one caused by the drug blocking a high affinity target, and the other due to off-target toxicity. The target-specific response correlates to the role a drug target plays in the proliferation of a cancer cell. Quantifying the target-specific effects of a given drug on cancer cells provides a mechanistic basis for a rational approach for identifying effective drug combinations for multi-driver cancers. We demonstrated that combinations of drugs, each blocking a specific driver, potently and synergistically block the viability of cancer cells that are refractory to inhibition by individual drugs.

Multi-Driver Cancers Are a Major Challenge to Targeted Therapy
Despite the dramatic successes targeted cancer therapy has enjoyed, it faces two major challenges: acquired resistance [55,56] and lack of sensitivity. A mono-driver cancer can be effectively treated initially with targeted therapy, but it acquires resistance to the same treatment after the initial sensitive period. The acquired resistance is due to the selection for cells that contain resistance-conferring mutations in the original driver, or new alternative signaling pathways being activated [57,58]. The acquired resistance can be overcome by developing drugs that can block the mutated driver or newly activated driver. Combination therapy that targets more than one step in a proliferative signaling pathway also helps prevent the development of acquired resistance. The second challenge is that most cancers are naturally refractory to targeted therapy. The lack of sensitive response suggests that a cancer is not addicted to a predominant driver or pathway, and multiple proliferative drivers are sustaining the proliferation. This is consistent with genetic evidence that most cancer cells contain multiple driver genes and driver mutations [4,6]. It is logical that only a combination of drugs, each blocking an independent driver, would be effective in stopping multi-driver cancers.
Despite considerable preclinical and clinical effort at developing combination therapy, effective combination cancer therapy is still elusive [7,23]. A recent analysis of many FDA-approved drug combinations revealed that most of the benefits are derived from patient-to-patient variability without drug additivity or synergy [23]. This was also found to be true for combinations tested in patient-derived tumor xenografts [23]. Thus, how to develop effective combination therapy to treat multi-driver cancers that are naturally insensitive to individual drugs is still a challenge. Underlying that challenge is the fact that how multi-driver proliferation responds to targeted therapy is still poorly understood [6].
Multi-driver proliferation can account for many of the existing challenges in targeted therapy. For example, Src kinases are shown to play important roles in the proliferation of various cancers [14,15], but clinical trials of Src inhibitors for such cancers have produced disappointing results [16,17]. The discrepancy between cell-based studies and clinical trials has been a major enigma in targeted therapy. Another example is the BRAF V600E mutant which is found in many cancer types. However, targeting BRAF is an effective treatment in melanoma [13], but not for colorectal cancers containing the same mutation [13,59]. A plausible explanation is that Src or BRAF are likely one of the multi-drivers in some cancers. Even though Src kinases are associated with many cancer types, they have not been shown to be an exclusive driver in any cancer. As we demonstrated in this study, BRAF V600E and Src are part of the multi-driver mechanism in the HT-29 cancer cell line. Furthermore, in most of the dasatinib-sensitive cancer cells, the dasatinib target, likely Src in many cases, is one of the drivers in multi-driver cells. For such multi-driver cancers, finding effective drug combinations to stop most, if not all proliferative drivers, is needed for effective combination targeted therapy.

Multi-Driver Proliferation Challenges Traditional Pharmacological Metrics for Analyzing Cancer Cell-Drug Response
Traditionally, the effect of a cancer drug on cancer cells is measured by the IC 50 , which works well for traditional drugs, and measures how potently a drug kills cancer cells. As previously noted [24] and further demonstrated in the current study, IC 50 does not capture how a targeted therapy agent affects a multi-driver cancer cell. Current mathematic models for calculating IC 50 are derivatives of the Hill equation which was developed to describe O 2 binding to hemoglobin. They are all based on the assumption of cooperative binding of a ligand to a single target. Two assumptions of the Hill equation are not consistent with the response of multi-driver cancer cells to targeted therapy. First, there is no evidence that kinase inhibitors bind to kinase with cooperativity. Second, kinase inhibitors tend to inhibit many kinase targets that are playing various roles in regulating cell viability and proliferation, thus inhibition of multi-driver cancer cells by targeted therapy is clearly a multi-target system. These considerations are clearly reflected in the dose response data of colorectal cancer cells. An analysis of the response of multi-driver cancer cells to kinase-based targeted therapy resulted in the evolution of the Hill equation into a biphasic mathematical model. Using HT-29 and other CRC cell lines, we demonstrated that multi-driver cancer cells respond to targeted therapy blocking one driver in two phases: a target-specific phase and an off-target toxicity phase. The biphasic response can be represented by three parameters, F 1 /F 2 ratio, K d1 and K d2 . Each parameter has a distinct mechanistic basis: F 1 /F 2 ratio represents how much a role a driver may play in the cell proliferation and viability; K d1 represents the affinity a drug has for the driver; and K d2 represents the potency of a drug's general toxicity to a host cell. The three-parameter mathematical model accurately describes the drug response of all the cancer cells we have examined so far. If a cancer cell is mono-driver dependent, the drug response becomes an extreme scenario, in which F 1 becomes 100% or nearly 100% and the second phase is eliminated. Analysis of cancer cell responses to dasatinib reveals that the ratio of biphasic versus monophasic responses is about 4:1.
The biphasic response is consistent with multi-driver proliferation and the nature of targeted therapeutics. When a cancer cell has multiple independent proliferative drivers, if we had a perfect targeted drug that exclusively interacts with the driver alone, the drug would only partially inhibit cell proliferation, resulting in a shallow or partial inhibition. However, most targeted therapy agents are protein kinase inhibitors, and there are over 500 similarly structured protein kinases that are potential targets for a kinase inhibitor. Each kinase inhibitor displays a varying level of promiscuity toward the protein kinome. Most protein kinase inhibitors in clinical use have toxicity at some dose. Thus, biphasic inhibition is a natural manifestation of the unique characteristics of multi-driver proliferation, the precision focus of targeted therapy, and the promiscuity of targeted drugs.

Biphasic Analysis Provides a Basis for Formulating Combination Targeted Therapy for Multi-Driver Cancers
As noted above, the biphasic analysis breaks cancer cell-drug response into three quantitative parameters, F 1 /F 2 , K d1 , and K d2 , each representing a mechanistic aspect of the cell-drug interaction. Quantifying these mechanistic aspects, combined with knowledge of drug-target specificity and signaling pathways of targeted kinases, offers a rational basis for predicting effective drug combinations. When more than one drug with distinct F 1 effects is identified against a cancer cell, the combination of these drugs has been found to be dramatically effective in blocking cell viability. Dose reduction of an order of magnitude or more is often achieved ( Figure 5A). This has led to the identification of two-drug (BMS-754807 and AZD-6244) and three-drug (BMS-754807, AZD-6244 and dasatinib) combinations that are broadly effective against colorectal cancer cells.
The CRC cell lines, HT-29, SK-CO-1, and NCI-H747 are potently inhibited by the combination of AZD-6244 and BMS-754807 with IC 50 s below 70 nM. The triple combination of AZD-6244/BMS-754807/dasatinib inhibited all three cell lines with IC 50 s at 40 nM or less ( Figure 5B). These results, together with a previous report from this lab [31] indicate that most CRC cell lines are sensitive to the AZD-6244/BMS-754807 combination (7/11), and even more sensitive to the AZD-6244/BMS-754807/dasatinib combination (8/11). This level of potency is comparable to those of the most successful targeted therapy drugs for mono-driver cancers. Further in vivo and clinical investigations are needed to test the effectiveness of these drug combinations on CRC.
It is yet to be tested how broadly applicable the biphasic analysis model is in predicting effective combination targeted therapy in vitro and in vivo beyond CRC. Multi-driver proliferation is believed to be broadly applicable to most cancers, because most cancers contain multiple driver mutations [4][5][6], and only a relatively small number of cancer subtypes can be treated effectively by mono-agent targeted therapy. However, this hypothesis needs to be directly and experimentally tested in broad cancer types. We have investigated two triple negative breast cancer cell lines and found that both relied on multiple drivers, and the biphasic mathematical model accurately describes their responses to targeted therapy, and predicts highly potent and strongly synergistic drug combinations [60]. Further investigations are needed to evaluate the clinical applicability of this approach.
Even though this current study focused on kinase-based target therapy agents, it is noted that targeted therapy can often be combined with chemotherapeutic agents to produce remarkable synergy [21]. It is possible that the synergistic targeted therapy combinations may be combined with chemotherapy to achieve the most favorable outcomes in clinical applications.

Biphasic Analysis Provides the Tool to Quantitatively Assess the Target-Specific and Off-Target Effects of Targeted Therapeutics
Quantification of the toxicity effect using biphasic analysis is also potentially crucial in developing targeted therapy, especially for multi-driver cancers. Even though it is well recognized that protein kinase inhibitors tend to have off-target effects, which are likely the basis of clinical dose-limiting toxicity, it is difficult to distinguish the target-specific effects from the off-target effects. The IC 50 -based evaluation measures the combined sum of both effects. Consequently, off-target effect contributes to a lower IC 50 , which makes a drug appear more potent in vitro, but diminishes the effectiveness of a drug in the clinical applications. This Trojan horse characteristic of the off-target effect may have contributed to the discrepancy between in vitro studies and the clinical results of targeted therapy. The biphasic analysis offers a quantitative tool to evaluate the target-specific and off-target effects. The quantitative understanding of the two effects of a drug offers a better guide to laboratory and clinical applications of targeted therapy agents.

Cell Culture and Viability Assays
Cell culture and viability assays were performed as described previously [18]. Briefly, cells were plated in 96-well plates at 25,000 cells per well, and were incubated with drugs at indicated concentrations in 100 µL medium containing 1% DMSO for 72 h. At the end of the incubation, 10 µL staining solution containing 5 mg/mL 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT dye) (Thermo Fisher Scientific) was added into each well, and the plates were incubated at 37 • C for another 4 h. The medium in the wells was removed by aspiration, and DMSO (100 µL) was added into each well. After the purple formazon was dissolved, absorbance at 490 and 750 nm for each well was determined using a VersaMax Tunable Microplate Reader. The A490-A750 values were taken as indicators of cell viability. All cell growth and drug inhibition experiments were performed in triplicate at least three times.

Drug Synergy Analysis and Combination Index Calculation
Drug synergy analysis was performed as described [31]. To determine if two drugs were synergistic in their inhibition of cell proliferation, cell viability was determined after incubating the cells in the presence of each drug alone or both drugs (1:1 ratio) at 16 concentrations ranging from 0.6 nM to 20 µM. The concentration of a drug or drug combination that inhibited cell viability at a certain percentage (IC x ) was determined manually from the graphs of drug dose response curves. The dose reduction index (DRI) was calculated according to Chou [45], using the Chou-Talalay method [45]. Briefly, the IC 40 , IC 50 , IC 60 , etc., were determined from the dose response curves, and the DRI at a given level of inhibition was calculated by the following equation: where IC x−1 , IC x−2 , and IC x−1+2 are the concentrations of drug 1, drug 2, and drugs 1 + 2, causing X% inhibition of cell viability, respectively.

Curve Fitting by Single Target Equation and Biphasic Equation
Curve fitting was performed in Microsoft Excel. The equation for the single target binding is I = [D]/([D] + IC 50 ), where I is the relative inhibition as the dependent variable, [D] is the drug concentration as the independent variable, IC 50 is the drug concentration that inhibits viability by 50%, a constant to be determined through curve fitting. In curve fitting, the root mean square error (RMSE) between the experimentally determined cell viability and the calculated cell viability from the equation at each drug concentration was minimized using the IC 50 as a variable. The IC 50 that resulted in the best fit with the smallest RMSE was taken as the best fitting IC 50 . The RMSE minimization was carried out using the Solver add-in program in Microsoft Excel.
The biphasic kinetic curve fitting was performed similarly, except that the equation used is a biphasic kinetic equation: I = F 1 × [D]/([D] + K d1 ) + F 2 × [D]/([D] + K d2 ). In this equation, the inhibition of cell viability (I) as a function of variable drug concentration is determined by four constants: F 1 and F 2 , the phase 1 and phase 2 as fractions of total cell viability, K d1 and K d2 , and the affinity of the drug for phase 1 and phase 2 targets. The equation assumed that the responses were biphasic, thus F 1 and F 2 added up to 1 (or 100%), so F 2 = 1 − F 1 . The RMSE between the actual cell viability (1 − I) and the calculated cell viability from the biphasic equation as a function of the drug concentration was minimized, allowing the F 1 , K d1 and K d2 as variables. The F 1 , K d1 and K d2 that resulted in the smallest RMSE were taken as the best fitting parameters. F 2 was calculated as 1 − F 1 .
The dose response data of HT-29 cells were also analyzed by a modified Hill equation: I = I max × D n /(IC 50 * n + D n ). The meaning of the parameters and assumptions of this equation are discussed in the Results Section 2.2. Curve fitting to this equation, which minimized RMSE by varying I max , n and IC 50 *, was performed in Microsoft Excel using the Solver function.

Western Blot Analysis of Drug Effects on Cell Signaling
To determine the effects of the kinase inhibitors on the signaling network in a given cell line, the cells were cultured to 70% confluency, and treated by the drugs at the indicated concentrations for 2 h. After the treatment, the culture medium was aspirated, and the cells were then resuspended in the SDS-PAGE loading buffer containing a protease inhibitor cocktail (Sigma-Aldrich, St Louis, MO, USA) and protein phosphatase inhibitors (PhosSTOP, Sigmal-Aldrich), and then lysed in a 95 • C heat bath immediately. SDS-PAGE and Western blotting were performed as described [31]. The loading was normalized based on the β-actin expression levels in the lysates determined by an antibody specific for this protein. All antibodies were purchased from Cell Signaling Technology (Danvers, MA, USA). The density of each protein band in the Western blots was measured using the ImageJ software (https://imagej.nih.gov/ij/).

Monophasic and Biphasic Analyses of Genomics of Drug Sensitivity in Cancer Cell Responses to Dasatinib
The dasatinib inhibition data (GDSC1-raw-data) were downloaded from the Genomics of Drug Sensitivity in Cancer website (https://www.cancerrxgene.org/downloads/bulk_download). The raw viability data were normalized to the controls as relative viability. In cases where the viability reading in the presence of low concentration drugs were higher than the controls, the highest reading was used as a control for normalizing the data. Monophasic and biphasic analyses were carried out as described in Section 4.4.

Conclusions
Targeted therapy is an effective treatment for cancers that rely on a single proliferative driver. However, most cancers are sustained by multiple independent proliferative drivers. Such multi-driver cancer cells are only mildly inhibited by targeted therapy agents. This study revealed that the inhibition of several colorectal cancer cell lines by targeted therapy does not follow a traditional single target binding equation, and the potency cannot be represented by IC 50 . A biphasic mathematical model was developed that quantifies the drug response using three inhibitory parameters, ratio of target-specific inhibition, F 1 , potency of target-specific inhibition, K d1 , and potency of off-target toxicity, K d2 . This biphasic model and the mechanistic insights it reveals provide a new mechanistic perspective and mathematical tool for predicting effective combination targeted blockades against multi-driver cancer cells. Further in vivo and clinical investigations are required to validate this rational approach.
The study identified an intriguing drug combination toward CRC. The CRC cell lines, HT-29, SK-CO-1, and NCI-H747, were potently inhibited by the combination of AZD-6244/ BMS-754807 (IC 50 s < 70 nM), and triple combination of AZD-6244/BMS-754807/dasatinib (IC 50 s < 40 nM). This level of potency is comparable to those of the most successful targeted therapy drugs for mono-driver cancers.