The complexity of cancer signaling networks limits the efficacy of most single-agent treatments and brings about challenges in identifying effective combinatorial therapies. In this study, we used chronic active B-cell receptor (BCR) signaling in diffuse large B-cell lymphoma as a model system to establish a computational framework to optimize combinatorial therapy in silico. We constructed a detailed kinetic model of the BCR signaling network, which captured the known complex cross-talk between the NFκB, ERK, and AKT pathways and multiple feedback loops. Combining this signaling model with a data-derived tumor growth model, we predicted viability responses of many single drug and drug combinations in agreement with experimental data. Under this framework, we exhaustively predicted and ranked the efficacy and synergism of all possible combinatorial inhibitions of eleven currently targetable kinases in the BCR signaling network. Ultimately, our work establishes a detailed kinetic model of the core BCR signaling network and provides the means to explore the large space of possible drug combinations. Cancer Res; 77(8); 1818–30. ©2017 AACR.
Using chronic active B-cell receptor signaling in diffuse large B-cell lymphoma as a model system, we developed a kinetic-modeling based computational framework for predicting effective combination therapy in silico. By integrative modeling of signal transduction, drug kinetics, and tumor growth, we were able to directly predict drug-induced cell viability responses at various dosages, which were in agreement with published cell line experimental data. We implemented computational screening methods that identified potent and synergistic combinations in silico and validated our independent predictions in vitro.
Quick Guide to Equations and Assumptions
A kinetic model was constructed for the BCR signaling network according to the following rules. For protein–protein binding interactions, we assumed that this type of reaction were under equilibrium and solved the steady-state level of protein complex analytically in the model,
where [TA] and [TB] stand for the total concentration of protein species A and B, respectively, [AB] represents the concentration of the protein complex, is the inverse of the dissociation constant Kd. For kinase catalyzed reactions, we adopted the classic Michaelis–Menten kinetics,
where is the rate of catalytic product formation, kcat is the enzyme turnover rate, Km is the Michaelis–Menten constant, [E] and [S] are concentration of enzyme and substrate, respectively. Each interaction in the BCR signaling network was depicted by corresponding terms according to above rules. The full kinetic model consists of 28 state variables, each representing concentration of a specific form of a protein species, depicted by 10 steady-state equations and 18 ordinary differential equations (see Supplementary Text S1).
A data-driven tumor growth model was formulated to link signaling response to cell viability output. We assume that tumor cell population is at exponential growth phase,
where N0 and N are cell number at time 0 and time t, is basal death rate, while growth rate is dependent on three downstream survival and proliferation signals NFκB, pAKT, and pERK (normalized by untreated control) through a Hill function,
Here, r* is basal growth rate, n is Hill coefficient, w1, w2, w3 are the estimated weights of three signaling outputs NFκB, pAKT, and pERK. Therefore, viability response defined as the ratio of cell number monitored under treated condition (for a time span of ) and untreated control can be predicted as following,
The activation of intracellular signaling pathways in response to environmental stimulus leads to important cell decisions such as proliferation. The amplitude and duration of pathway activation are precisely and robustly controlled by complex regulatory loops to maintain cellular homeostasis. In cancer, activating mutations or deletion of signaling repressors frequently result in sustained and exaggerated pathway activation that drives uncontrolled tumor survival and proliferation. Targeted therapies that use small-molecule inhibitors to repress specific signaling pathway members, for example, kinases, can directly block oncogenic pathway activation and lead to tumor cell death. These targeted therapies are expected to provide improved efficacy and reduced toxicity compared with chemotherapy.
However, clinical application of targeted therapy is facing several challenges such as low response rate and frequently acquired drug resistance. The limited efficacy of single-agent–targeted therapy is at least partially due to pathway cross-talks and compensatory circuits within signaling networks targeted by these agents (1). Cross-talks and compensatory circuits allow signals to bypass drug inhibition and reactivate downstream effectors. By simultaneously repressing multiple nodes in a signaling network, combination therapy has the potential to completely extinguish oncogenic signaling and induce more potent and durable treatment response. Thus, novel drug combinations where two or more drugs work cooperatively to suppress corrupted signaling networks need to be identified to achieve maximum therapeutic efficacy. The complexity of signaling networks makes it difficult to simply guess which combinations will be effective and synergistic and which ones will not. Moreover, given the large number of possible drug combinations against complex signaling networks, comprehensive experimental screening, including exploration of multiple dosages, is not practically feasible. Besides, results acquired from such screening may be specific to the cell line tested, thus lacking general applicability to highly variable primary tumors found in patients.
Computational models of signaling networks that can accurately reconstruct signaling dynamics in silico may represent a useful alternative to experimental screening and trial-and-error experimental investigation. Once proven reliable, these models can be used to exhaustively test the efficacy of a large number of single drug and drug combinations by quantifying signaling output under corresponding network perturbations. Even though computational modeling has been widely used to study the dynamics of signaling network in the past decades, the development of cancer signaling models and its application to predicting effective combinatorial therapies is still lacking. Here we demonstrate the feasibility of this approach using chronic activation of B-cell receptor (BCR) signaling in diffuse large B cell lymphoma (DLBCL) as a model system. We adopted a systems biology approach and established a computational framework to optimize anti-DLBCL combinatorial therapy in silico. The proposed approach is broadly applicable and can be used for other malignancies driven by aberrantly active signaling pathways.
The deregulation of BCR signaling is central to the pathogenesis of many B-cell malignancies. It is especially central in the activated B-cell–like subtype of DLBCL (ABC DLBCL). ABC DLBCLs exhibit chronic active BCR signaling, and are addicted to constitutive activation of downstream survival and proliferation signals such as NFκB (2). It has recently been found that a subset of the germinal center B-cell like subtype of DLBCLs are also dependent on BCR signaling through activation of the PI3K/AKT pathway (3). Multiple small-molecule inhibitors against BCR signaling were developed and proved effective in killing BCR-dependent DLBCLs in vitro and in vivo (4–6). However when tested in clinical trials, single-agent treatments again demonstrated limited responsiveness and efficacy (7), suggesting an urgent need for the design of effective combination therapies.
In this work, we present a kinetic modeling–based computational framework for predicting and optimizing combinatorial therapy against chronic active BCR signaling (Fig. 1). We constructed a detailed kinetic model of the BCR signaling network parameterized by published signaling responses and protein concentrations. Mathematical models of proximal BCR signaling and downstream transcriptional network have been reported (8–10). But, to our knowledge, this is the first kinetic model to reconstruct the entire core BCR signaling network in silico. Using published drug response data in a BCR signaling–dependent cell line, we trained a tumor growth model, which in combination with the kinetic model, allowed us to simulate viability response upon various targeted treatments. Under this framework, we exhaustively tested the efficacy and synergism of all possible combinations of inhibition of eleven currently targetable kinases in the BCR signaling network. We discuss how these results pave the way for the discovery of effective drug combinations.
Materials and Methods
Protein–protein interactions in the BCR signaling network
Antigen-induced BCR crosslinking allows SRC family kinases, mainly LYN, to phosphorylate the immuno-receptor tyrosine–based activation motifs (ITAM) of the intracellular BCR subunits Igα (CD79A) and Igβ (CD79B; ref. 11). Dually phosphorylated ITAM motifs then recruit SYK and activate it via tyrosine phosphorylation (12). Activated SYK phosphorylates adapter BLNK, which recruits BTK to the plasma membrane to facilitate its phosphorylation and subsequent activation by SYK and LYN (13). Activated BTK further phosphorylates PLCγ2, which catalyzes the hydrolysis of phosphatidylinositol-4,5-bisphosphate (PI(4,5)P2) into diacylglycerol (DAG) and inositol trisphosphate (IP3; ref. 14). DAG together with elevated intracellular calcium induced by IP3 triggered endoplasmic reticulum (ER) calcium-release activates PKCβ (15), which then stimulates two divergent pathways that activate NFκB and ERK, respectively. Phosphorylation of CARD11 by PKCβ leads to the assembly of the CBM complex composed of CARD11, BCL10, and MALT1 (16). CBM acts as a scaffolding complex that facilitates IKK phosphorylation by TAK1, which in turn phosphorylates IKB and induces its degradation, releasing NFκB into the nucleus to elicit transcriptional activity (17). In addition, protease activity of MALT1 positively regulates NFκB signaling by cleaving and inactivating inhibitors against NFκB activation such as A20 and RELB (18, 19). In the meantime, PKCβ and DAG activate RASGRP, which triggers the canonical MAPK signaling cascade, leading to eventual phosphorylation and activation of ERK (20). On the other hand, SYK and LYN phosphorylate BCAP and CD19, respectively, which activate PI3K by membrane recruitment (21, 22). PI(3,4,5)P3 synthesized by PI3K further facilitates PDK1-catalyzed AKT phosphorylation by binding to both proteins via their plextrin homology (PH) domains (23). Importantly, LYN negatively regulates PI3K signaling by activating SHIP1, which hydrolyzes PI(3,4,5)P3 into PI(4,5)P2 (24).
Besides major signal transduction pathways as described above, our model includes key regulatory interactions in the BCR signaling network such as pathway cross-talks and feedback loops. The PI3K pathway positively regulates NFκB and ERK signaling by enhancing BTK membrane recruitment via PI(3,4,5)P3 binding. At the same time, it conversely attenuates ERK signaling via AKT catalyzed RAF phosphorylation (25). It has recently been found that MEK negatively regulates PI3K/AKT signaling by recruiting PTEN to the plasma membrane (26), which dephosphorylates PI(3,4,5)P3 into PI(4,5)P2. BTK amplifies BCR signaling by two coupled positive feedback loops. It recruits PIP5K to the plasma membrane, which produces PI(4,5)P2 to sustain both PI(3,4,5)P3 synthesis and PI(4,5)P2 hydrolysis (27). In addition, BTK phosphorylates BCAP, further facilitating the membrane recruitment of PI3K (22). The activity of BTK is attenuated by active PKCβ via disruption of its membrane localization, constituting a negative feedback loop (28). Besides, another indirect feedback from PKC to SYK was added into the model as knockdown of PKCδ was shown to mediate hyperphosphorylation of SYK (29). Furthermore, multiple negative feedback loops exist within the MAPK signaling cascade to fine-tune its activation amplitude and duration (30).
Cell viability assay
In the published drug response experiment (31), cells were treated at time 0 and incubated for 48 hours. Viability response was normalized to the plate positive control (bortezomib) and negative control(DMSO) as described previously (31). These normalized data were used for comparison with simulation results.
For our independent validation of synergistic and antagonistic drug combinations, the DLBCL cell lines TMD8, HBL1, and OCI-LY10 were generously donated from the Melnick laboratory at Weill Cornell Medicine (New York, NY) in 2014. These cell lines were verified using the cell line authentication service at Biosynthesis (BSI) using their STR Profiling and Comparison Analysis Service (32) in 2015. TMD8, HBL1 cell lines were cultured using RPMI1640 media supplanted with 10% FBS, 5% HEPES, 5% l-glutamine, and 5% penicillin/streptomycin. OCI-LY10 cell line was cultured using Iscove's Modified Dulbecco's Media supplanted with 10% FBS, 5% l-glutamine, and 5% PS. All cell culture media and reagents were from Life Technologies. R406, dasatinib, and MI-2 were purchased from Selleck Chemicals. Cells were grown at concentrations sufficient to keep untreated cells in exponential growth during the time of drug exposure. Cells were treated with 6 doses of each drug or combination in triplicate. Drug combinations were administered in constant ratio. Cell viability was determined by an ATP luminescent method (CellTiter-Glo, Promega). Luminescence was measured with the Syngery4 microplate reader (BioTek). Cell viability in drug-treated cells was normalized to vehicle-treated controls.
Kinetic model of BCR signaling network
Upon antigen-induced BCR crosslinking, most kinases involved in the BCR signaling pathway are recruited to the plasma membrane through activation of adaptor proteins such as BLNK, BCAP, and CD19 to form an integrated signalosome (33). Therefore, we expect most interactions depicted in our kinetic model to happen in a localized fashion. Thus, we do not consider any compartmentalization and model the cell as a well-mixed system.
As protein–protein binding is in general a very fast process at the second time scale, we assumed that this type of reaction were under equilibrium and solved the steady-state level analytically in the model. For a reversible protein–protein binding reaction,
where [A] and [B] stand for the concentration of freed form of A and B; [TA] and [TB] stand for the total concentration of A and B; [AB] represents the concentration of the complex. By solving the above three equations, we have
where , is the inverse of the dissociation constant Kd.
Under a few circumstances where a protein may bind to more than one partner, the interactions were considered independently for simplification.
For kinase catalyzed reaction, we adopted the classic Michaelis–Menten kinetics,
where is the rate of catalytic product formation, kcat is the enzyme turnover rate, Km is the Michaelis–Menten constant, [E] and [S] are concentration of enzyme and substrate, respectively.
Reactions in the BCR signaling network were written into corresponding equations according to rules discussed above. The full model consists of 28 state variables each representing concentration of a specific form of a protein species, depicted by 10 steady-state equations and 18 ordinary differential equations (ODE; see Supplementary Text S1). Kinetic parameters in the model are summarized in Supplementary Table S1. Bounds for kinetic parameter estimation are listed in Supplementary Table S2. Parameters of total protein concentrations are summarized in Supplementary Table S3. We performed parameter sensitivity analysis where each parameter was perturbed independently across four orders of magnitudes and viability response was recorded. We found overall robustness and identified the most sensitive parameters as parameters regulating main axis of the NFκB pathway (see Supplementary Fig. S1).
We imposed temporal pLYN and pSYK stimulus in normal BCR signaling modeled by two double exponential functions,
where and are total concentration of LYN and SYK, respectively. and are the phosphorylated fraction of LYN and SYK, respectively. and were estimated by fitting to the phosphorylation time course data using a genetic algorithm. A,B,C,D were estimated by fitting to the normalized pLYN and pSYK time course.
In contrast, we imposed constitutive pLYN and pSYK stimulus in simulations of diseased ABC DLBCLs, with negative feedback from SYK and PKC, respectively.
The 0.2 coefficient is to account for LYN attenuation effect due to CD79B mutation in TMD8 (2). The parameter for negative feedback is estimated by fitting to single-drug viability response.
Tumor growth model
Assume a tumor cell population is at exponential growth phase (34, 35),
where N0 and N are cell number at time 0 and time t, is basal death rate, while growth rate is dependent on three downstream survival and proliferation signals NFκB, pAKT, and pERK (normalized by untreated control) through a Hill function,
Here r* is basal growth rate, n is Hill coefficient, w1, w2, w3 are the estimated weights of three signaling outputs NFκB, pAKT, and pERK. Therefore, viability response defined as the ratio of cell number monitored under treated condition (for a time span of ) and untreated control can be predicted as following,
Parameters required in this function were trained with viability data of three single-drug viability responses, namely NFκB, AKT, and MEK inhibitor, respectively. First the level of the three downstream survival and proliferation signals were predicted via simulation of the signaling model, and then input into the tumor growth functions to compute the viability output. Parameters were chosen by minimizing the sum of residuals between the viability predictions and experimental data.
To simulate an inhibitor's effect at a given dosage, percent activity of the targeted kinase was calculated via the medium effect equation,
The drug's IC50 was taken from literature (Supplementary Table S4), while m was assumed to be 1 under a first-order approximation. Then, the activity of the targeted kinase (i.e., parameters representing catalytic or activation rate of targeted kinase) was reduced to the corresponding percentage in the kinetic model. We list perturbed parameters in each simulated inhibitor treatment in Supplementary Table S4.
Under the Bliss Independence model, the additive effect of two inhibitors is computed as the multiplication of the effect of individual inhibitors,
where indicates the fraction unaffected. To evaluate mode of interaction between two inhibitors, we computed viability response at 10×10 virtual dosages by varying the percent inhibition of each targeted kinase independently from 0% to 90% at 10% interval. These viability values were used to estimate parameter that minimizes the following metric,
where x,y are virtual dosages for inhibitor 1 and 2, respectively. indicates synergism, additive, and antagonism, respectively.
The numerical simulations of ODEs were performed in Matlab _R2012a while data analysis was performed in R. Specifically, we used ode15s function in Matlab for solving ODEs and ga function in Matlab for parameter estimation using genetic algorithms. We have uploaded scripts for ODE simulation and parameter estimation in github (https://github.com/weiduweillcornell/BCR-signaling-model).
Kinetic modeling of BCR signaling network reproduces normal BCR signaling in silico
We first curated the central BCR signaling network by gathering experimentally validated protein–protein interactions from literature. The reconstructed network is shown in Fig. 2, and includes three major signaling pathways downstream of BCR, namely NFκB, PI3K/AKT, and RAF/RAS/ERK. We chose to include these three pathways because they have been reported to closely regulate cell survival and proliferation in B cells and B-cell malignancies (36). Detailed information on the interactions included in the constructed BCR signaling network is provided in Materials and Methods.
Instead of directly applying mass action kinetics to characterize elementary reactions in the network, we chose to adopt more streamlined mathematical representations derived from mass action law under reasonable assumptions (see Materials and Methods). This strategy greatly reduced the number of variables, equations, and most importantly parameters required in the mathematical model. As elementary protein–protein binding reactions generally reach equilibrium within seconds, we modeled them by deriving the steady-state relationships from mass action law (see Materials and Methods). For enzymatic reactions such as phosphorylation or dephosphorylation, we adopted a classic Michaelis–Menten kinetic framework. To parameterize the model, we first retrieved protein concentrations in B lymphocytes quantified by mass spectrometry from the MOPED protein expression database (37). We modeled LYN and SYK as two independent input signals that triggered a downstream response. Their activation kinetics were approximated by double exponential functions (see Materials and Methods) where parameters in the functions were estimated by fitting to the phosphorylation time course data (38). We then used genetic algorithms to optimize the remaining 72 kinetic parameters within bounded biologically reasonable ranges (Supplementary Table S2; refs. 39–41) by minimizing residual sum of squares between simulated phosphorylation time courses and published Western blot data (38). Experimental data and simulated results were each normalized to their respective maximum value for comparison. Ten sets of kinetic parameters were identified from 5,000 independent runs that fit almost equally well (Supplementary Table S5). Simulated trajectory under these 10 parameter sets together with phosphorylation time course data are shown in Fig. 3A. To examine whether our fitting method is prone to overfitting problem, we manually injected different levels of white noise into the training time course data and refitted the model using our genetic algorithm. The simulation was done using optimized parameter values as initial population of parameter choice and 10 independent optimization processes were performed for each noise level. We showed that the residual sum of square greatly increases as the noise level increases (Supplementary Fig. S2). This indicates that the model is not prone to fitting to noise. Parameter sensitivity analysis was performed as described in Materials and Methods (Supplementary Fig. S1). We note that predicted activation time course of BLNK, PLC, and PKC deviate from experimental data. We believe that the most likely reason behind this relatively poor fit may be currently unknown regulatory interactions involving additional proteins in the BCR signaling pathway. However the relatively poor fit in these components in the upstream of NFκB pathway would not significantly influence the model's predictions as our parameter sensitivity analysis indicates viability response is most sensitive to parameters in the downstream part of the NFκB pathway (Supplementary Fig. S1). Of note, we were independently able to find 39 kinetic parameter values from the literature, and we compared these values with the range of estimated 10 parameter sets (Fig. 3B). By shuffling the literature-retrieved parameters 10,000 times, we found that the literature-retrieved parameter values fall within the estimated parameter ranges significantly more often than random (P = 0.05; Fig. 3C). We note, however, that many discrepancies were found between estimated and published parameters (Fig. 3B). We speculate that many of these discrepancies are likely due to in vitro nature of the experiments used to quantify kinetic parameters.
Combining BCR signaling model with a tumor growth model predicted cell viability response upon single and combinatorial drug treatments in a BCR signaling–dependent ABC DLBCL cell line
We next sought to simulate the effect of various small-molecule inhibitors on ABC DLBCL cell viability and to compare simulation results with published drug response data in a BCR signaling–dependent ABC DLBCL cell line TMD8 (31). We selected TMD8 because of the extensive drug combinatorial data available on this cell line (31). We first made several modifications to the model to accommodate the differences between normal BCR signaling and aberrant BCR signaling in ABC DLBCL. Instead of applying a temporal upstream stimulus, we assumed constitutive LYN and SYK phosphorylation as observed both in ABC DLBCL cell lines and in primary DLBCL patient samples (see Materials and Methods; refs. 2, 42). In addition, we accounted for genetic alterations in members of the BCR signaling network in TMD8 compared with normal B cells. Specifically, TMD8 was shown to carry CD79B mutation that attenuates LYN activity by approximately 80% (2). Correspondingly, we decreased the enzymatic activity of LYN in the model to the same extent (see Materials and Methods).
To predict cell viability response from signaling output, we formulated a tumor growth model in which the growth rate of tumor cells is dependent on the weighted sum of the three downstream survival and proliferation signals NFκB, ERK, and AKT through a Hill function (see Materials and Methods). A similar formalism has been used to characterize tumor growth of ERBB-amplified breast cancer driven by ERK and AKT activation (43). We used published viability response data in TMD8 to parameterize the tumor growth model, where cells were treated with IKK, AKT, and MEK inhibitors at multiple dosages (31). Specifically, using the median effect equation (44), we calibrated the percent activity left on the targeted kinase for each inhibitor at a given dosage based on the inhibitor's IC50 value (see Materials and Methods; Supplementary Table S4). We then reduced the activity of the targeted kinase to the same level in the model and simulated steady-state signaling output. Parameters in the tumor growth model were estimated by minimizing residual sum of squares between predicted viability response and experimental data.
We first simulated single-drug viability response of inhibitors covering the NFκB, PI3K/AKT, and MAPK pathway and compared with experimental data. We observed that in silico simulation with the BCR signaling model and the tumor growth model recapitulated the viability response of the three training single-drug response, namely IKK, AKT, MEK inhibitors (Fig. 4A; ref. 31). This is not surprising as the growth model was fitted on the basis of training data. As independent predictions, we also simulated drug response of inhibitors targeting other kinases in the network, for example, CAL-101 against PI3K, ibrutinib against BTK, and found that predicted viability response matched favorably with TMD8 drug response data (Fig. 4B; ref. 31). At the same time, we found simulated viability response of SYK inhibition to deviate from experimental data (gray line), yet this discrepancy can be partially rescued by adding a negative feedback from SYK to LYN (blue line). It has been reported that SYK functions as a negative regulator of BCR signaling by phosphorylating Ig-α (29, 45, 46). As Ig-α primarily interacts with LYN, we assumed in the model that SYK indirectly negatively regulates LYN.
Beyond single-drug viability response, we also simulated combinatorial drug response of ibrutinib in combination with various other kinases targeting the BCR pathway, and observed the predicted response contour to match favorably with experimental results (Fig. 5). These results demonstrate that our model can correctly capture the interaction between inhibitors as well.
Overall, these results suggest that viability response of small-molecule inhibitors targeting the BCR signaling network can be predicted via in silico simulation of the BCR signaling model in combination with the tumor growth model.
Cross-talk between PI3K and NFκB pathway mediates efficacy of PI3K inhibition in TMD8
In both the drug response data and model's simulation, we observed that PI3K inhibition is significantly more effective at inhibiting tumor growth than blockage of its downstream effector AKT. A similar phenomenon was reported in other studies, where PI3K inhibition was shown to attenuate NFκB transcriptional activity (3, 47). We hypothesized that the efficacy of PI3K inhibition is primarily attributed to suppression of NFκB signaling, which is mediated by upstream cross-talk between the PI3K and NFκB pathways. To test this hypothesis, we abolished the cross-talk between PI3K and NFκB by knocking out in silico PI(3,4,5)P3-mediated membrane recruitment of BTK in the signaling model. Under this condition, we resimulated the viability response of PI3K inhibition, which showed significant reduction compared with both experimental data and simulation with the full signaling model (Fig. 4, brown line). This result supports the notion that the upstream crosstalk between PI3K and NFκB pathway is critical in mediating tumor growth inhibition by PI3K inhibitor. It also provides further rational support for the clinical use of PI3K inhibitors in DLBCL that are dependent in NFκB signaling (3, 47).
Computational optimization of targeted therapy against chronic active BCR signaling
Using the above modeling framework, we sought to identify targeted therapies against the BCR signaling network that most effectively inhibit tumor growth. We exhaustively tested all drug pairs based on 11 small-molecule inhibitors currently available that target various kinases in the network, yielding 55 treatment strategies in total. In each scenario, viability response was simulated at 10×10 virtual dosages where each targeted kinase was inhibited at 0% to 99% evenly spaced in log10 space. We calculated area under the combinatorial viability response surface as an overall indicator of drug combination potency. The smaller the value is, the more potent the drug target combination is (Fig. 6A). We found that under the same inhibition potency, efficacy of different treatment strategies was highly variable, ranging from almost no growth inhibition to up to 80% reduction (Fig. 6B). Specifically, inhibiting downstream of the NFκB signaling pathway, especially through MALT1 and IKK inhibitor, exhibited the most prominent efficacy, and combined MALT1 and IKK blockage yielded highest tumor growth inhibition. In comparison, tumor cell growth was relatively insensitive to blockage of MAPK pathway in our simulations. In summary, this computational screening result suggests that various treatment strategies against a signaling network can yield highly variable therapeutic responses and that in silico simulation can help identify targets that confer intrinsic vulnerability.
We then sought to identify drug combinations that are synergistic via computational simulations. For a given two-drug combination, the combinatorial drug response at 10× 10 virtual dosage as discussed above were used to estimate mode of drug interaction under the Bliss independence model (see Materials and Methods; Fig. 7A). Computational screening predicted dual blockage of LYN and SYK as the most synergistic combination. To test this prediction, we treated ABC DLBCL cell lines TMD8, HBL1, and OCI-LY10 with LYN inhibitor dasatinib and SYK inhibitor R406, at multiple doses. Comparing combinatorial drug response data to theoretical additive response predicted by the Bliss independence model (see Materials and Methods), we confirmed synergism between dasatinib and R406 (Fig. 7B). We also tested and confirmed the predicted antagonism between dual SYK and MALT1 inhibition using SYK inhibitor R406 and MALT1 inhibitor MI-2 across three ABC DLBCL cell lines (Fig. 7B).
It is increasingly acknowledged that aberrant BCR signaling plays a central role in the development and maintenance of many B-cell malignancies (48). Although a large panel of small-molecule inhibitors against BCR signaling have been developed, rational methodologies that can predict effective combinatorial therapy and guide the design of specific treatment strategy in individual patients have been lacking. To bridge this gap, we aimed to construct the first kinetic model of the core BCR signaling network and use this model to investigate targeted therapy against BCR signaling. We showed that simulations with the signaling model reconstructed dynamics of normal B-cell signaling in silico. Combining the signaling model with a data-trained tumor growth model successfully predicted viability response of multiple drug combinations, and identified novel synergistic drug combinations such as LYN and SYK inhibitor.
In this work, we chose to develop a highly detailed kinetic model of the BCR signaling network, modeling direct protein–protein interactions using highly quantitative ODEs wherever possible. We note that there are other simpler modeling techniques for analysis of signaling pathways, for example, simplifying signaling cascades as coupled Hill functions (43) or using Boolean network models. To investigate the predictive power of a simpler model, we constructed a Boolean model of a simplified BCR signaling network consisting of only the 11 targetable nodes together with signaling pathway end points and a node representing cell viability (Supplementary Fig. S3). In this model, each node has only two states (0, inactive; 1, active), and the state value in a particular step is determined by the values of all its regulators in the previous step. A node will be active if the majority of its regulators are activating, except that the Viability node is active when any of the signaling outputs is active. Assuming the two input nodes LYN and SYK are constitutively active, all possible initial states (2^12) are exhaustively simulated until reaching attractors. These initial states ended up in two attractors, both attractors consisting of states in which Viability is on (Supplementary Fig. S4). Then all two drug combinations were tested. To simulate drug-mediated inhibition, each targeted node was made constitutively inactive. There are four drug combinations that result in a single global attractor in which Viability is off. These four drug combinations are BTK-MEK, BTK-RAF, SYK-MEK, and SYK-RAF. These four drug combinations ranked low in the drug efficacy predictions made by the full ODE model (Fig. 6B). Moreover, experimental data in TMD8 cell line (31) indicates that the BTK-MEK and BTK-RAF combinations are not very effective in decreasing cell viability—the MEK and RAF inhibitors are not responsive (along the y-axis) and synergism is lacking (Supplementary Fig. S5). Overall, these results suggest that a simple Boolean network model is not able to capture the same results as the full ODE model and that its prediction accuracy is lacking when compared with experimental data.
The clinical application of targeted therapy is frequently challenged by highly variable drug response among cancer patients. Heterogeneous response to BCR signaling–targeted therapy was observed in ABC DLBCL cell lines (47) and in clinical trials of ABC DLBCL patients (7). This is likely due to differential expression of proteins in the BCR signaling pathway that impact pathway activities in individual patient. In this study, signal transduction is explicitly modeled on the molecular level, which provides a straightforward framework to incorporate protein level variation to develop patient-specific predictive models. Specifically, expression levels of proteins within the BCR signaling network can be measured experimentally using protein expression profiling techniques such as reverse-phase protein array (RPPA) or other approaches. Then protein expression changes in patients relative to cell line can be incorporated into the model to predict optimal treatment strategy for individual patients. We believe that the use of patient-specific predictive models can greatly improve the performance of targeted therapy in cancer.
In this study, we focused on the exploration of combinatorial efficacy of kinase inhibitors against BCR signaling pathway components. In the future, it may be compelling to combine kinase inhibitors with traditional chemotherapeutic drugs in the clinical setting. Regrettably, capturing the combinatorial effect of kinase inhibitors and chemotherapeutic drugs likely requires a complex theoretical framework that may involve a more detailed mechanistic characterization of cell cycle and DNA damage response than the simple growth model we use here.
Finally, we note that besides DLBCL, aberrant BCR signaling was shown to play a role in other B-cell malignancies such as chronic lymphocytic leukemia (CLL; ref. 49) and mantle cell lymphoma (MCL; ref. 50). In phase II studies of BTK inhibitor ibrutinib, 71% and 68% overall response rate (ORR) was reported in CLL and MCL patients, respectively (51), suggesting targeting BCR signaling as promising treatment strategy. Correspondingly, the overall framework and predictions reported in this work may also be useful to identify drug combinations for CLL and MCL-targeted treatment.
Disclosure of Potential Conflicts of Interest
L. Cerchietti reports receiving a commercial research grant from Celgene. A.M. Melnick is a consultant at Epizyme and reports receiving a commercial research grant from Janssen, GSK, and Eli Lilly. No potential conflicts of interest were disclosed by the other authors.
Conception and design: W. Du, A.M. Melnick, O. Elemento
Development of methodology: W. Du, O. Aly
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): R. Goldstein, Y. Jiang, O. Aly
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): W. Du, R. Goldstein, L. Cerchietti, A.M. Melnick
Writing, review, and/or revision of the manuscript: W. Du, O. Aly, A.M. Melnick, O. Elemento
Study supervision: L. Cerchietti, O. Elemento
This work was supported by the CAREER grant from National Science Foundation (DB1054964), NIH grant R01CA194547, the Starr Cancer Consortium and Hirschl Trust. A.M. Melnick is supported by NCI R01 CA143032, R01 CA104348, the Burroughs Welcome Foundation, and the Chemotherapy Foundation.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
We thank all Elemento and Melnick laboratory members for productive discussions.
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).
- Received February 17, 2016.
- Revision received December 10, 2016.
- Accepted January 23, 2017.
- ©2017 American Association for Cancer Research.