Abstract
The objective of this study is to evaluate the heterogeneity in pharmacodynamic response in four in vitro multiple myeloma cell lines to treatment with bortezomib, and to assess whether such differences are associated with drug-induced intracellular signaling protein dynamics identified via a logic-based network modeling approach. The in vitro pharmacodynamic-efficacy of bortezomib was evaluated through concentration-effect and cell proliferation dynamical studies in U266, RPMI8226, MM.1S, and NCI-H929 myeloma cell lines. A Boolean logic–based network model incorporating intracellular protein signaling pathways relevant to myeloma cell growth, proliferation, and apoptosis was developed based on information available in the literature and used to identify key proteins regulating bortezomib pharmacodynamics. The time-course of network-identified proteins was measured using the MAGPIX protein assay system. Traditional pharmacodynamic modeling endpoints revealed variable responses of the cell lines to bortezomib treatment, classifying cell lines as more sensitive (MM.1S and NCI-H929) and less sensitive (U266 and RPMI8226). Network centrality and model reduction identified key proteins (e.g., phosphorylated nuclear factor-κB, phosphorylated protein kinase B, phosphorylated mechanistic target of rapamycin, Bcl-2, phosphorylated c-Jun N-terminal kinase, phosphorylated p53, p21, phosphorylated Bcl-2–associated death promoter, caspase 8, and caspase 9) that govern bortezomib pharmacodynamics. The corresponding relative expression (normalized to 0-hour untreated-control cells) of proteins demonstrated a greater magnitude and earlier onset of stimulation/inhibition in cells more sensitive (MM.1S and NCI-H929) to bortezomib-induced cell death at 20 nM, relative to the less sensitive cells (U266 and RPMI8226). Overall, differences in intracellular signaling appear to be associated with bortezomib pharmacodynamic heterogeneity, and key proteins may be potential biomarkers to evaluate bortezomib responses.
Introduction
Multiple myeloma (MM) is a hematologic cancer of rapidly proliferating terminally differentiated plasma cells originating in the bone marrow. It is the second most common hematologic cancer, accounting for 10%–15% of such cases (Röllig et al., 2015; Rajkumar and Kumar, 2016). The current 5-year survival rate for people with the disease is approximately 47% (https://www.cancer.net/cancer-types/multiple-myeloma/statistics) (Naymagon and Abdul-Hay, 2016). Symptomatic myeloma manifests as tissue or organ dysfunction, primarily renal insufficiency, hypercalcemia, bone disease, and anemia. The presence of at least 10% of clonal bone marrow plasma cells and abnormal monoclonal protein in serum or urine confers diagnosis (Durie et al., 2003, 2006). Drug treatment often includes chemotherapeutic agents, such as proteasome inhibitors (e.g., bortezomib and carfilzomib), immunomodulators (e.g., lenalidomide and thalidomide), monoclonal antibodies (e.g., daratumumab), histone deacetylase inhibitors (e.g., panabinostat), steroids (e.g., dexamethasone), and alkylating agents (e.g., melphalan) (Naymagon and Abdul-Hay, 2016).
The evolution of MM goes through asymptomatic stages starting from monoclonal gammopathy of undetermined significance, progressing to smoldering myeloma (SM), and finally to symptomatic myeloma. The risk of progression from monoclonal gammopathy of undetermined significance to SM, and SM to myeloma is approximately 1% and 10% per year (Kyle et al., 2007; Landgren et al., 2009). The development of MM is a complex multistep transformation process involving the progressive acquisition of genetic abnormalities (Morgan et al., 2012). The primary genetic events include chromosomal translocations on chromosome 14, specifically at the Ig switch region, and hyperdiploidy. Disease progression continues with the occurrence of further secondary genetic abnormalities (e.g., chromosome deletions/gains), epigenetic events (e.g., hypomethylation/hypermethylation), and acquired mutations (Bergsagel and Kuehl, 2001; Fonseca et al., 2002). This results in the emergence of several clonal and subclonal populations of abnormal plasma cells that differ in their genetic composition (Walker et al., 2012). Clonal progression appears to resemble the branching evolution model proposed by Darwin in describing the “origin of species” (Nowell, 1976). The clones/subclones exist in equilibrium and evolve in complexity under the influence of the tumor microenvironment and treatment regimens (Keats et al., 2012). The profiles of clonal dynamics in small subsets of patients using high-throughput sequencing techniques have revealed the existence of molecularly different clones at diagnosis, during treatment, and at relapse (Egan et al., 2012; Keats et al., 2012; Magrangeas et al., 2013). The composition of the clonal populations has also been found to change under therapeutic treatment, microenvironment influence, and the natural progression of the disease (Fuhler et al., 2011; Egan et al., 2012; Magrangeas et al., 2013). Overall, the clonal evolution of myeloma is the basis for its heterogeneous nature.
Systems biology/pharmacology approaches provide platforms for studying the complexity associated with intracellular signaling and identifying mechanisms of drug action, targetable pathways, biomarkers, feedback mechanisms, and combinatorial drug regimens (Van der Graaf and Benson, 2011; Iyengar et al., 2012; Harrold et al., 2013). Most chemotherapeutic agents work via modulating intracellular signaling, and the development of signal transduction networks is a mechanistic approach to linking pharmacokinetics and pharmacodynamic efficacy of drugs and regimens (Kirouac and Onsum, 2013). The identification of appropriate biomarkers using these approaches is needed given the heterogeneity and diversity in responses to current treatment in patients and the desire to individualize pharmacotherapy.
The objective of this study is to examine heterogeneity in pharmacodynamic responses of four MM cell lines to bortezomib treatment and evaluate the association between heterogeneous responses and drug-induced protein signaling using an approach that integrates traditional pharmacokinetic/pharmacodynamic modeling, systems pharmacology, and a high-throughput multiplexing assay. Four genetically diverse mammalian MM cell lines: U266, RPMI8226, MM.1S, and NCI-H929 were used to evaluate the influence of cellular heterogeneity on bortezomib efficacy/sensitivity. Differences between cell lines were tested using model fitting and the estimation of pharmacodynamic response parameters, such as the IC50 for bortezomib obtained from concentration-effect data and Kmax and KC50 (nonlinear second-order capacity and sensitivity cytotoxicity terms) obtained from modeling temporal cell proliferation dynamics.
To achieve a systems-level understanding of the signaling mechanisms and key proteins stimulated by bortezomib in myeloma cells, a published network-based approach investigating bortezomib protein signaling dynamics in U266 myeloma cells (Chudasama et al., 2015) was adopted and extended to add additional signaling mechanisms and establish a cell line–unbiased network for myeloma cell growth and proliferation. The network was qualified with experimental protein expression measurements for the newly added pathways and prior qualification criteria of the base model. Network simulations, centrality measures, and a reduction algorithm identified critical proteins to further test the hypothesis that sensitivity differences among cell lines upon bortezomib treatment are associated with differences in expression patterns of intracellular signaling proteins. The temporal profiles of network-identified proteins were measured using a high-throughput magnetic bead–based multiplexing protein assay. This study highlights the need to understand the mechanistic basis of heterogeneous responses in myeloma to support drug efficacy assessment, and the inclusion of cell-based biomarker measurements might be used ultimately to tailor individual pharmacotherapy.
Materials and Methods
Cell Culture and Reagents
The human myeloma cell lines U266 (catalog number TIB–196), RPMI8226 (catalog number CCL–155), MM.1S (catalog number CRL–2974), and NCI–H929 (catalog number CRL–9068) were purchased from American Type Culture Collection (ATCC; Manassas, VA). Cells were maintained in culture in RPMI 1640 medium (catalog number 30–2001; ATCC) supplemented with fetal bovine serum (catalog number 30–2020; ATCC) and penicillin-streptomycin (BioWhittaker; Lonza, Walkersville, MD). NCI-H929 cells were also supplemented with 2-mercaptoethanol (Sigma-Aldrich; St. Louis, MO). Bortezomib was purchased from LC Laboratories (Woburn, MA) with chemical purities greater than 99%. Bortezomib was dissolved in dimethylsulfoxide (DMSO; Sigma-Aldrich) to obtain a 0.1 µM stock solution, which was stored at −80°C. The cell proliferation reagent WST-1 assay kit was purchased from Roche Diagnostics (Indianapolis, IN). Cellular protein isolation reagents included radioimmunoprecipitation assay buffer (10×) (Cell Signaling Technology, Danvers, MA) supplemented with Halt Protease and Phosphatase Inhibitor Cocktail (100×) (Thermo Fisher Scientific, Rockford, IL) and phenylmethylsulfonyl fluoride (Thermo Fisher Scientific). Total protein concentrations were measured using the DC Protein Assay (BIO-RAD, Hercules, CA). Reagents for protein expression analysis using the MAGPIX System included magnetic bead microspheres, detection antibody, streptavidin-phycoerythrin, assay buffer, and amplification buffer, all of which were purchased from EMD Millipore Corporation (Saint Charles, MO).
Concentration-Effect Studies
The initial cell seeding density and the incubation time after the addition of cell proliferation reagent WST-1 were optimized for each cell line in preliminary experiments based on the manufacturer assay kit protocol instructions. Approximately 10,000 (U266), 2500 (RPMI8226), 30,000 (MM.1S), and 25,000 (NCI-H929) cells per well (96-well plate) were incubated with bortezomib at concentrations ranging from 0.01 to 100 nM after a preincubation period of 24 hours for acclimatization of cells to the 96-well plate. The different concentrations were obtained by diluting the stock solution of bortezomib in DMSO with cell line–specific cell culture medium. At 0, 24, 48, and 72 hours after bortezomib exposure, cells were incubated with the cell proliferation reagent WST-1 for 2 hours. At the end of the incubation period, cell viability was measured at 450 and 690 nm (reference wavelength) with a Microplate Spectrophotometer (Molecular Devices, Sunnydale, CA). Cell viability was reported as a percentage of untreated cells (control) and reported as six replicates.
Cell Proliferation Studies
Approximately 10,000 (U266), 2500 (RPMI8226), 30,000 (MM.1S), and 25,000 (NCI-H929) cells per well (96-well plate) were incubated with bortezomib at concentrations of 2, 4, 10, and 20 nM after a preincubation period of 24 hours for the acclimatization of cells to the 96-well plate. The different concentrations were obtained by diluting the stock solution of bortezomib in DMSO with cell line–specific cell culture medium. At 0, 2, 4, 8, 12, 24, 30, 48, 72, and 96 hours after bortezomib exposure, cells were incubated with the cell proliferation reagent WST-1 for 2 hours. At the end of the incubation period, cell viability was measured at 450 and 690 nm with a Microplate Spectrophotometer (Molecular Devices). Cell viability was normalized to the untreated cells (control) at 0 hour and reported as six replicates.
MAGPIX System–Based Protein Assay
Experimental Design.
All four myeloma cell lines (U266, RPMI8226, MM.1S, and NCI-H929) in culture media were incubated in several 10-cm2 culture dishes at a density of 5 × 106 cells/10 ml of culture media. Three treatment groups were included: vehicle control (culture medium); 2 nM bortezomib (only RPMI8226 and MM.1S cell lines); and 20 nM bortezomib (all cell lines). Proteins were measured at 0, 2, 6, 12, 24, 32, 48, and 72 hours after bortezomib exposure. At each time point, cells in culture media were collected, temporarily stored on ice (to stop cellular reactions), and centrifuged to remove media and isolate them. Cells were then lysed at 4°C in radioimmunoprecipitation assay buffer supplemented with Halt Protease and Phosphatase Inhibitor Cocktail (100×) and phenylmethylsulfonyl fluoride for 30 minutes. After lysis, the samples were centrifuged at 11,000 rpm for 20 minutes at 4°C, and the supernatant was stored at −80°C until further analysis. Total protein concentrations in each sample were analyzed with the DC Protein Assay (BIO-RAD), which is similar to the Lowry assay (Lowry et al., 1951). The entire experimental design was repeated three times in total for replicated measurements.
The MAGPIX Protein Assay is based on the Luminex xMAP Technology (Luminex Corporation, Austin, TX), which is a magnetic bead–based assay platform. Multiple fluorescent dye–based proprietary color-coded microspheres (magnetic beads) coated with specific capture antibodies are used to capture several target proteins in a test sample. Upon capturing the target protein, a biotinylated detection antibody is introduced, followed by incubation with streptavidin-phycoerythrin conjugate (reporter molecule). A light-emitting diode–based analysis system measures the fluorescence intensity of the reporter to quantify proteins. The magnetic bead–based system enables multiple protein assay measurements from the same sample.
Time Course of Protein Expression.
The four myeloma cell lines treated with bortezomib were analyzed for the longitudinal expression of the following proteins: phosphorylated NFκB (pNFκB), total NFκB, phosphorylated AKT (pAKT), total AKT, phosphorylated mTOR (pmTOR), total mTOR, Bcl-2, phosphorylated JNK (pJNK), total JNK, phosphorylated p53 (pp53), total p53, p21, pBAD, caspase 8, and caspase 9. These proteins were chosen based on the key proteins identified using the Boolean network model, which includes proteins governing cell survival, proliferation, and apoptosis in myeloma cells. Proteins were purchased as both a premixed panel (human early apoptosis magnetic bead 7-plex kit containing pAKT, Bcl-2, pJNK, pp53, pBAD, caspase 8, and caspase 9) and individual assays, which were mixed per manufacturer instructions after considering compatibility with the magnetic bead detection regions and assay buffer. The phosphorylated and total proteins were measured separately (different 96-well plates) since the analytes share the same magnetic bead regions and cannot be mixed.
Protein Assay Procedure.
The assay procedure was performed according to the manufacturer protocol. A total protein concentration of 20 μg was used for analyzing intracellular proteins in cell lysates. Briefly, cell lysates were first diluted with the assay buffer (1:1) to achieve a final total protein amount of 20 μg. The diluted cell lysates (in 96-well plate) were then incubated overnight (16–20 hours) with magnetic beads in a shaker at 4°C in dark. The samples were later washed twice with assay buffer before incubation with the detection antibody for 1 hour with shaking in the dark at room temperature. Thereafter, the detection antibody was removed followed by the addition of streptavidin-phycoerythrin and incubation for 15 minutes at room temperature (with shaking under dark conditions). Amplification buffer was then added to the wells (without removing streptavidin-phycoerythrin) with another 15 minutes of incubation at room temperature (with shaking under dark conditions). After incubation, streptavidin-phycoerythrin and amplification buffer were removed, and the beads were resuspended in assay buffer prior to being read on the MAGPIX instrument. Three replicates were analyzed for all samples. Manufacturer-provided positive and negative control cell lysates were also analyzed along with experimental samples in addition to the assay buffer control (blank).
Data Analysis.
The reported median fluorescence intensity values were used for data analysis. The averaged median fluorescence intensity value of the assay buffer control (blank) in duplicate for each protein was subtracted from the median fluorescence intensity of the samples. The data were normalized by the housekeeping protein (glyceraldehyde-3-phosphate dehydrogenase). All of the data were also normalized to the protein expression in the 0 hour–untreated control cells to obtain the relative protein expression in bortezomib-treated cells. In addition, the phosphorylated proteins were normalized by their total protein equivalents.
Mathematical Modeling
a. In Vitro Concentration-Effect Model.
The concentration-effect curves for the four cell lines at three time points (24, 48, and 72 hours) were modeled individually, and the data were modeled using a standard inhibitory effect function (Hill, 1910):(1)with R representing cell viability (response) expressed as the percentage of control, R0 is the baseline response in the absence of bortezomib, Imax is the maximum inhibition of cell viability, C is the concentration of bortezomib, IC50 is the IC50 of bortezomib, and γ is the Hill coefficient.
b. Cell Proliferation Model.
Cell viability profiles upon treatment with bortezomib (2, 4, 10, and 20 nM) were modeled simultaneously within each cell line. Briefly, myeloma cells proliferate at an exponential rate represented by a first-order growth rate constant (kg), and bortezomib causes cell death represented by a nonlinear cytotoxicity term comprising the parameters Kmax and KC50. This term is similar to the Michaelis-Menton–type equation, where Kmax is the maximum killing rate constant, and KC50 is the Michaelis-Menton constant equivalent to the concentration of bortezomib that results in half-maximal Kmax (Yano et al., 1998):(2)Because the data were normalized to untreated-control cells at 0 hour, the first-order growth rate constants for each cell line were fixed to that estimated from modeling the control data alone. The model also incorporates a first-order degradation rate constant (kdeg) for bortezomib in the cell culture medium (Chudasama et al., 2015):(3)with C0 as the concentration of bortezomib at time 0 hour.
All model fitting and parameter estimation was performed using the maximum likelihood estimation algorithm in ADAPT5 (D’Argenio et al., 2009). The variance model was defined as follows:(4)with σ1 and σ2 as estimated variance model parameters, and Y (θ, ti) is the ith predicted value from the pharmacodynamic model. Model performance was assessed based on the following goodness-of-fit criteria: model convergence, precision of parameter estimates, and visual inspection of predicted versus observed values and residual plots.
c. Network Model Construction
A Boolean logic–based network model describing intracellular signaling in the U266 myeloma cell line served as the primary base model (Chudasama et al., 2015). Upon extensive literature review, the network was extended with the inclusion of additional signaling pathways [e.g., mTOR pathway and insulin-like growth factor type 1 (IGF-1)/insulin receptor substrate (IRS) signaling axis] and expansion of existing signaling pathways [e.g., cell cycle and renin-angiotensin system (RAS) pathways]. The overall process of building/expanding the network primarily entailed the following two steps: 1) building an interaction network; and 2) translating the interactions to obtain a logic-based network by adding Boolean relationships to pathway regulation (Wynn et al., 2012). The base network was modified to include an additional 24 nodes and several interactions. Two nodes [i.e., reactive oxygen species (ROS) and MnSOD2 (Superoxide dismutase 2, mitochondrial)] were also adopted from a bortezomib-vorinostat combination interaction network (Nanavati, 2016). Existing nodes and their relationships were revised based on the interactions with newly added pathways and literature-reported mRNA/protein expression studies (Wynn et al., 2012; Saadatpour and Albert, 2013). Once the interaction network was in place, information flow through the interaction network was specified with logic functions (also known as transfer functions or logic gates) that define subsequent state space as the network is updated at discrete time steps. These logic functions use one or more of the Boolean operators “AND” and “OR,” which define the processing of multiple signals acting on a single node. The “AND” operator simply reflects that the activation of a particular node depends on the simultaneous activation of all the nodes regulating it, whereas the “OR” operator requires activation of only one of the regulating nodes for activation to occur (Wynn et al., 2012). Network building also required assignment of initial states or initial conditions (0 or 1) for each node of the network (Chudasama et al., 2015). The initial conditions were carefully assigned based on their functional relevance (inactive or active) in an untreated proliferating myeloma cell. Many of the growth and proliferative pathway nodes were assigned an initial condition of 1 (active), and nodes pertaining to cell death/apoptosis were assigned an initial condition of 0 (inactive). Boolean relationships/rules governing any node are implemented in time steps, such that the state of a given node at time “t + 1” is dependent on the state of all the nodes controlling it at time “t.”
d. Network Simulations
The discrete Boolean network was converted into a system of ordinary differential equations to conduct continuous dynamic simulations of the network. Odefy is a MATLAB-compatible toolbox that implements the HillCube technique based on multivariate polynomial interpolation and Hill kinetics to transform Boolean models into a system of first-order ordinary differential equations (Wittmann et al., 2009; Krumsiek et al., 2010). Time takes on arbitrary units as Boolean networks operate in discretized time scales. Simulations were implemented in MATLAB as described by Krumsiek et al. (2010) using the default parameter values available in the “HillCubeNorm” simulation option (normalizes Hill functions between 0 and 1). Simulations of the network were performed under untreated conditions (control, bortezomib “off”/inactive/0) and drug treatment conditions (bortezomib: “on”/active/1) via direct inhibition of the proteasome node and stimulation of the RIP (Receptor Interacting Protein) node (Chudasama et al., 2015).
Network Qualification: In Silico Perturbations
Network qualification of the base model was achieved by comparing model simulations with experimental measurements of proteins via immunoblotting (Chudasama et al., 2015). The model was simulated under the presence of a dual mTORC inhibitor (acting via inhibiting the mTORC1 and mTORC2 nodes simultaneously) alone and in combination with bortezomib. Simulations were compared with pAKT measurements and cell viability measurements published in the literature (Hoang et al., 2010). The model was also qualified with all the specific probes and proteins that were used to qualify the base model (Chudasama et al., 2015).
e. Network Topology and Centrality Measures
The final network diagram was built in yEd (yWorks, GmbH, Tübingen, Germany), which was also used to evaluate network topology and for performing centrality analysis. Centrality analysis is known to be useful for identifying key players regulating biologic processes (Koschützki and Schreiber, 2008; Pavlopoulos et al., 2011). The total number of nodes and edges, degree centrality (based on number of outgoing connected edges), and betweenness centrality of the network were determined using yEd.
Network Reduction.
To identify critical proteins involved in bortezomib pharmacodynamics, the Boolean network reduction method proposed by Veliz-Cuba (2011) was manually applied to the final network. The method implements three basic reduction steps applied recursively: 1) nodes with only one input and one output are removed to connect the upstream and downstream nodes directly; 2) algorithm S: removal of nonfunctional edges; and 3) Algorithm R: removal of a nonfunctional nodes. This approach was also applied to the base model (Chudasama et al., 2015). Dynamical simulations of the reduced network under untreated-control and bortezomib treatment conditions were performed using Odefy, and simulations of the full and reduced networks were compared to ensure consistency in the steady-state simulations of the reduced nodes.
f. Statistical Analysis
Spearman’s rank correlation coefficient matrix signatures were calculated between cell viability values under bortezomib treatment and relative protein expressions for all four cell lines in SigmaPlot (version 12.0; Systat Software Inc., San Jose, CA). Hierarchical clustering of the correlation matrix was also performed with the default Euclidean distance as the similarity measure in MATLAB (R2016b). The correlation matrices and the clustered dendrogram of the correlation matrix for each cell line were plotted in MATLAB (R2016b).
Results
Concentration-Effect Relationships
The observed and model-fitted (eq. 1) concentration-effect profiles for four myeloma cell lines are shown in Fig. 1. Bortezomib exhibits concentration- and time-dependent cell-killing effects, with increasing cell death with increasing drug concentrations and longer durations of drug exposure across all cell lines. The estimated model parameters are listed in Table 1. The Imax and γ terms were fixed to 1 and 5 whenever the estimated values were close to 1 (Imax) or greater than 5 (γ). All parameters were estimated with good precision, and the model described the data well (Fig. 1). A comparison of the parameters across cell lines showed a consistently lower IC50 value for bortezomib in MM.1S cells at all time points relative to U266, RPMI8226, and NCI-H929 cells. The estimated Imax values for all cells were around 1, indicating maximum inhibition at high bortezomib concentrations. At the 24-hour time point, only MM.1S cells showed an Imax value of 1.
The differences in bortezomib sensitivity across cell lines were tested statistically using standard assessments of model predictive performance (Sheiner and Beal, 1981). Briefly, data were co-modeled so that all cell lines were permitted to share one common IC50 parameter versus cell line specific estimates. The mean squared prediction errors and the differences in mean prediction errors were calculated and tested for statistical significance at different α levels (Supplemental Table 1). The calculated mean squared prediction errors for the individual models were lower than that of the co-model with a shared potency, indicating that the individual IC50 values were significantly different. The difference in mean squared prediction error between the individual and the co-model analysis was significant at an α level of 0.05 for MM.1S cells at 24 and 72 hours and for NCI-H929 at 24 hours. This suggests greater bortezomib sensitivity for MM.1S cells relative to the other three cell lines. Overall, the concentration-effect analysis identified relative bortezomib potency as MM.1S > NCI-H929 ≥ RPMI8226 > U266 at all three time points.
Cell Proliferation Studies
The time-course of cell proliferation also exhibited concentration- and time-dependent cell kill effects with increasing bortezomib concentrations and longer exposures in all four cell lines (Fig. 2). At the lowest bortezomib concentration (2 nM), no difference in cell viability across cell lines was observed until 48 hours (Supplemental Fig. 1). This difference at later time points was magnified with 4 nM bortezomib, with separation in curves beginning as early as 24 hours. MM.1S and NCI-H929 cells showed consistently greater reduction in viability over time in comparison with U266 and RPMI8226 cells. At 4 nM, complete cell death was observed only in MM.1S cells at 72 and 96 hours. At higher bortezomib concentrations (10 and 20 nM), a clear separation was observed, with MM.1S and NCI-H929 cells demonstrating near complete cell death at 12 hours post-treatment, and U266 and RPMI8266 cells showing complete cell death at 30 hours post-treatment. The differences in the time-course of cell proliferation across cell lines were tested statistically using a three-way analysis of variance test. The effects of cell lines (U266, RPMI8226, MM.1S, and NCI-H929), bortezomib concentrations [4, 10, and 20 nM (2 nM was not included as no separation between the cell lines was observed at this concentration)], and time (0–96 hours) were examined on the mean normalized cell viability (dependent variable). A statistically significant effect of all three variables was observed with P values <0.001. A pairwise multiple-comparison test (Holm-Sidak method) identified significant differences among the cell lines U266 versus RPMI8226/MM.1S/NCI-H929 and RPMI8226 versus MM.1S/NCI-H929; between bortezomib concentrations of 4 versus 10/20 nM, and among many time points (except among a few earlier and later time points), confirming the differences in sensitivity among the four cell lines.
Cell proliferation data for each cell line were modeled separately, with the model capturing the temporal profiles for all four concentrations (Fig. 2), and model parameters were estimated with reasonable precision (Table 2). The exponential in vitro kg values for all four cell lines were similar, indicating that the cell lines do not differ in their unperturbed growth kinetics. The estimated Kmax values for NCI-H929 and MM.1S cells were greater than those for U266 and RPMI8226 cells. The estimates for KC50 values are rank ordered: NCI-H929 > MM.1S > U266 > RPMI8226; which is similar to those for Kmax values. The ratio of Kmax and KC50 was similar in all four MM cell lines. Overall, cell proliferation studies suggest increased sensitivity of MM.1S and NCI-H929 cells to bortezomib treatment.
Network Development
The final comprehensive Boolean network model incorporating intracellular signaling processes that regulate MM cell survival, growth, and apoptosis is depicted in Supplemental Fig. 2. The Boolean relationships for all nodes, their initial states, and the associated literature references are reported in Supplemental Table 2. Although the base network included several pathways relevant to growth and survival [e.g., NFκB, janus kinase (JAK)/signal transducers and activators of transcription (STAT), and phosphoinositide 3-kinase (PI3K)/AKT], new mechanisms relevant to myeloma signaling dynamics were added. Overall, 24 new nodes were included belonging to the mTOR pathway (e.g., mTORC1, mTORC2, Rheb, TSC, DEPTOR, and p70S6K) (Acosta-Jaquez et al., 2009; Tchevkina and Komelkov, 2012); IGF-1- and IRS-mediated signaling (Qiang et al., 2002; Abroun et al., 2004; Ishikawa et al., 2006; Vivanco et al., 2007; Chiron et al., 2013); interferon regulatory factor 4 (IRF4) (Shaffer et al., 2009; Rui et al., 2011); and cell cycle regulation–associated nodes [e.g., Cdc25A/B, checkpoint kinase (Chk) 1/2, and CYCB] (Otto and Sicinski, 2017). Some of the nodes already in the base that required revision of their logic relationships included, but were not limited to: AKT, pNFκB, RAS, pRB (retinoblastoma protein), MYC, p21, p53, and BAD (Tan et al., 1999; Harada et al., 2001; Hermeking, 2003; Giacinti and Giordano, 2006; Guertin et al., 2006; Ishikawa et al., 2006; Dhanasekaran and Reddy, 2008; Shaffer et al., 2009; Rui et al., 2011; Busino et al., 2012; Maurer et al., 2014; Otto and Sicinski, 2017). Notably, the pNFκB node in the base model incorporated stimulation by a dummy node “X” to account for bortezomib stimulation of pNFκB (contrary to reports in the literature bortezomib-mediated inhibition of pNFκB) (Hideshima et al., 2009; Chudasama et al., 2015). This dummy node was replaced by a more mechanistic explanation of pNFκB being stimulated under stress conditions in a cell, which was elaborated in the bortezomib-vorinostat combination interaction network (an iteration of the base model developed for the specific combination treatment of myeloma) (Huang et al., 2003; Dai et al., 2005; Nanavati, 2016). This mechanism required the adoption of two nodes (ROS and MnSOD2) into the interaction network.
Network Qualification
The approach used for qualification of the network was the same as that used for the base model. Briefly, considering that the mTOR pathway was added to the base network, it was necessary to qualify the model with probes specific to this pathway. mTOR inhibitors in combination with other antimyeloma agents have exhibited efficacy in myeloma. Hoang et al. (2010) tested a dual mTOR inhibitor, pp242 (inhibits both mTORC1 and mTORC2), in combination with bortezomib in different MM cell lines. They observed that pp242 is synergistic with bortezomib primarily because it inhibits mTORC2-mediated stimulation of AKT. The authors report the percentage of apoptosis and immunoblots for pAKT (Ser473) and total AKT in RPMI8226 and MM.1S cells, respectively. Our network was simulated in the presence of a dual mTOR pathway inhibitory probe, and the simulations for AKT and apoptosis nodes were in agreement with the experimental data. As expected, dynamical simulations of logic-based networks cannot differentiate synergism in an observed effect. The model simulation for apoptosis under combination perturbation overlapped the profile of single-agent bortezomib treatment. The model simulations are shown in Supplemental Fig. 3. Simulations were performed under the influence of pathway-specific probes [IKKi (inhibitor of kB kinase inhibitor) and JAKi (JAK inhibitor) alone and in combination] to determine the dynamics of downstream nodes [pIκBα (phospho-inhibitor of NFκB), pNFκB, pStat3, and Bcl-xL], and these were compared with experimental protein measurements via immunoblotting. Simulations of the final network were consistent with the base model under similar conditions.
Bortezomib Pharmacodynamics
The qualified network was used to investigate the intracellular signaling mechanisms altered by bortezomib that eventually trigger myeloma cell death. The effect of bortezomib was applied on the “Proteasome” and “RIP” nodes with bortezomib directly inhibiting “Proteasome” and stimulating the “RIP” node (Supplemental Fig. 2). The model was simulated to observe the effect of bortezomib on the dynamics of nodes (proteins) relative to the simulations obtained in the absence of drug treatment (mimics the untreated-control proliferative myeloma cell). Simulations were in agreement with bortezomib mechanisms of action, including the following: stimulation of JNK (stress pathway), stimulation of p53 and p21 (cell-cycle inhibitors; arrests cell division), stimulation of pNFκB (prosurvival protein), stimulation of BAD (proapoptotic protein), inhibition of Bcl-xL and Bcl-2 (antiapoptotic proteins), and activation of caspase-9 (intrinsic mitochondrial apoptosis pathway). In addition, several new mechanisms were also highlighted, such as: stimulation of AKT and mTOR (prosurvival), stimulation of DNAdam and Chk1/Chk2 (DNA damage response pathway), and activation of caspase-8 (extrinsic apoptosis pathway). These mechanisms suggest that bortezomib induces cell stress by inhibiting proteolysis causing DNA breaks, generation of ROS, and activation of the extrinsic apoptosis pathway (via stimulation of JNK), along with mitochondrial apoptosis, coercing the cells to antagonize the effect of drug by stimulating their prosurvival pathways. These outcomes were substantiated with evidence from the literature (Mitsiades et al., 2002; Qin et al., 2005; Que et al., 2012; Kubiczkova et al., 2014; Mimura et al., 2014), adding confidence in the network performance and the ability to use networks to study system dynamics.
Network Topology and Centrality Measures
The final network consists of a total of 97 nodes and 202 edges (Supplemental Fig. 2). The nodes are color coded based on their role in the cell (e.g., nodes highlighted in green are prosurvival proteins, whereas nodes highlighted in orange are proapoptotic proteins). The nodes with a yellow fill color represent the newly added nodes in comparison with the base model. Overall, the nodes of the model are schematically positioned such that the nodes belonging to each pathway are grouped together. For example, the top panel of Supplemental Fig. 2 includes nodes associated with prosurvival pathways (e.g., mTOR, PI3K/AKT, NF-κB, JAK/STAT, RAS, IGF-1/IRS, and IRF4) in that order starting from left to right with the growth factors on the top of the figure. The bottom panel of Supplemental Fig. 2 includes the following: intrinsic/extrinsic apoptosis nodes, DNA damage and stress pathway nodes, and cell cycle/proliferation–associated nodes in that order from left to right.
Centrality analysis of the network was performed based on the following two methods, 1) the number of connected edges (measure of how well the nodes interact or are associated/connected with each other); and 2) node betweenness centrality (a measure of how often a node lies on the path between other nodes) (Koschützki and Schreiber, 2008; Pavlopoulos et al., 2011). The top 10 ranking nodes identified from the analysis are reported in Table 3. The nodes are scored relatively with the highest-ranking node having a score of 1. The two assessments of node centrality revealed overlapping nodes, such as p53, AKT, pNFκB, JNK, MYC, and caspase 3, representing proteins from important prosurvival, proapoptotic, and stress-related pathways. These results emphasize the significance of some of these well-known intracellular signaling proteins and pathways, and these nodes may represent potential biomarkers for bortezomib pharmacodynamics.
Network Reduction
Network reduction attempts to identify nodes that are central and critical to the network while preserving its dynamical properties. The reduced network obtained after applying the reduction algorithms is shown in Fig. 3 and consists of 19 nodes and 37 edges. Steady-state simulations of the reduced network agreed with that of the larger final network (Fig. 3). Many of the nodes from the reduced network overlapped with the nodes identified from the centrality measures, re-emphasizing the importance of critical prosurvival (AKT, pNFκB, and mTOR), stress-associated (JNK), proapoptotic (caspase 9), and cytostatic (p53, p21) proteins in regulating the pharmacodynamic effects of bortezomib. This provides a concise list of proteins that might serve as potential biomarkers of bortezomib effects across myeloma cell lines to study heterogeneous drug responses.
Protein Dynamics
The overall intracellular signaling protein dynamics across the four myeloma cell lines are shown in Fig. 4, which provides a snapshot of all 10 proteins in untreated-control and bortezomib-treated cells. The data represent log-transformed values of the mean relative expression obtained from three separate replicate experiments. The four columns represent the individual myeloma cell lines: MM.1S, RPMI8226, U266, and NCI-H929. The three row panels on the y-axis are relative expressions (from 0-hour untreated-control cells) in untreated-control, and bortezomib 2 nM and 20 nM treatments. The x-axis denotes time in hours. Given the logarithmic transformation of the data normalized to the 0-hour untreated-control cells in each cell line, the time-course data for the untreated-control arm varies around a mean value of zero (Fig. 4, green). Any stimulation/increase in expression in Fig. 4 is represented as red, whereas inhibition/decrease in expression is denoted as blue as specified in the key.
The log-relative expression of proteins in the untreated-control arm of all cell lines is around zero, with some minor variability in a few proteins. This variability is insignificant considering that these are replicated experiments relative to treatment arms. The 2 nM bortezomib treatment did not result in any stimulation or inhibition of protein expression in RPMI8226 cells. However, stimulation and inhibition of proteins was observed in MM.1S cells; especially p21 and caspase 8 stimulation at later time points (32–72 hours). Prosurvival proteins, such as pNFκB and pAKT, exhibit an initial delay, followed by minor stimulation toward later time points. pJNK exhibits a gradual activation throughout the duration of treatment. pmTOR, Bcl-2, and pp53 remain unchanged, whereas pBAD follows the untreated-control curve up to 32 hours, followed by minor inhibition at 48 and 72 hours. Caspase 9 shows stimulation that dampens after 24 hours. The line plots of the time-course of relative expression of proteins upon administration of 2 nM bortezomib in MM.1S and RPMI8226 cells are shown in Supplemental Figs. 4 and 5. All four cell lines exhibit substantial stimulation or inhibition in protein profiles at the higher bortezomib concentration (20 nM). MM.1S and NCI-H929 cells exhibit relatively greater and earlier onset of activation/inhibition of proteins in comparison to U266 and RPMI8226 cells.
Protein Expression Dynamics (20 nM Bortezomib)
pNFκB.
The temporal dynamics of the prosurvival protein pNFκB relative expression across cell lines is depicted in Fig. 5. An initial delay (12 hours for U266 and RPMI8226; 6 hours for MM.1S and NCI-H929) in stimulation of the protein is observed in all the four cell lines. MM.1S and U266 cells show a gradual rise over time, whereas NCI-H929 and RPMI8226 cells display maximal activation at 24 hours and between 24 and 48 hours, followed by dampening of the signal at 72 hours. MM.1S cells exhibit a 150-fold increase in stimulation in relative expression compared with control cells, which is much greater than those of the other cell lines (20-fold in NCI-H929, 40-fold in U266, and about 2-fold in RPMI8226).
pAKT.
AKT is a serine/threonine-specific protein kinase that plays a major role in cell growth and survival. A stimulation of pAKT is observed in all four myeloma cell lines after an initial delay of 2–6 hours in activation. U266 and RPMI8226 cells displayed a slower rate of stimulation in contrast to the other two cell lines. The overall peak fold change in expression was similar across cell lines and ranged from 8-fold to 15-fold. The trend followed attainment of a peak activation followed by a sustained signal with minor diminishment (Supplemental Fig. 6).
pmTOR.
The dynamics of another proliferative pathway associated protein, pmTOR, are shown in Supplemental Fig. 7. The overall trend indicated delayed stimulation of the protein under bortezomib treatment in U266 and MM.1S cells. No change from untreated-control was observed in RPMI8226 and NCI-H929 cells. All profiles were also associated with very high variability among replicates, making it difficult to assess the magnitude of activation.
Bcl-2.
Bcl-2 is an antiapoptotic protein belonging to the Bcl-2 family of apoptosis-regulating proteins. Bortezomib treatment induces a suppression of Bcl-2 in all four myeloma cell lines, as expected. The suppression was very gradual in U266 and RPMI8226 cells compared with a faster rate of inhibition in MM.1S and NCI-H929 cells. Unlike U266 and RPMI8226 cells, a near complete inhibition was observed by 12 hours in MM.1S and NCI-H929 cells. U266 and RPMI8226 cells showed only 50% inhibition even until 72 hours (Fig. 6).
pJNK.
JNKs belong to the mitogen-activated protein kinase family, and these proteins are activated under stress stimuli. Bortezomib is known to induce the expression of pJNK, as confirmed in the four cell lines investigated in this study (Supplemental Fig. 8). The induction was associated with a slight delay (2–6 hours) in all cell lines; a maximum fold change in induction was observed for U266 cells (close to 30-fold to 40-fold). Other cell lines exhibited an induction ranging from 6-fold to 10-fold. Peak activation was observed near 24–48 hours in U266 and RPMI8226 cells, whereas MM.1S and NCI-H929 cells achieved peak activation at 6 hours after bortezomib treatment. U266 cells showed sustained induction after the peak signal throughout the time-course, whereas the rest of the cells displayed a dampening of signal after the peak.
pp53.
The dynamics of p53 in the four myeloma cell lines are shown in Fig. 7. U266 and RPMI8226 cells presented a clear 3-fold to 4-fold stimulation of the protein, with sustained stimulation in RPMI8226 cells and reduced stimulation at later time points in U266 cells. No change in protein expression was observed for MM.1S and NCI-H929 cells throughout the time-course.
p21.
p21 is a cyclin-dependent kinase (CDK) inhibitor that inhibits CDK1 and CDK2. The major role of p21 is arresting or inhibiting cell cycle progression. The tumor suppressor p53 activates and tightly regulates the expression of p21. The profiles for p21 relative expression in the four cell lines suggested stimulation, but profiles were associated with very high S.D. and variability. The overall directionality was toward activation, with a clear stimulation in U266 cells peaking at 12 hours, followed by loss of the protein at later time points. The peak fold change in relative expression was 6-fold. The other cell lines presented unclear profiles, with an indication of activation. A clear stimulation of p21 was observed in MM.1S cells at the lower-concentration bortezomib (2 nM) treatment. The p21 relative expressions in the four cell lines are plotted in Supplemental Fig. 9.
pBAD.
BAD is a proapoptotic protein of the Bcl-2 family. Apoptosis induction in cells is associated with a stimulation of BAD mediated by dephosphorylation of pBAD. BAD and pBAD are in equilibrium, and cellular apoptosis is related with an activation of BAD and therefore an inhibition of pBAD. pBAD (Ser112) was measured, and an inhibition was observed in all four cell lines upon bortezomib treatment. The inhibitory curves are very similar to that of Bcl-2, with gradual inhibition in U266 and RPMI8226 cells. The rate of inhibition is much faster in MM.1S and NCI-H929 cells. The onset of activation of inhibition is also earlier in these cells compared to U266 and RPMI8226 cells (Supplemental Fig. 10).
Caspase 8.
Caspase 8 is a marker for activation of the extrinsic apoptotic pathway. Bortezomib has been associated with the stimulation of caspase 8 via pJNK (Mitsiades et al., 2002; Qin et al., 2005). This stimulation was confirmed via the protein dynamics in all four myeloma cell lines (Fig. 8). A strong stimulation of caspase 8 was observed, especially in MM.1S cells, with a peak stimulation of 150-fold at 12 hours, followed by a reduction and stabilization at 50-fold expression. A similar profile was observed in NCI-H929 cells with peak stimulation (20-fold to 25-fold) at 12 hours and stabilization at 5-fold at later time points. The magnitude of peak stimulation at 12 hours was much lower in the case of U266 (8-fold to 10-fold) and RPMI8226 (2-fold to 2.5-fold) cells. No change in expression from control cells was noted at later points in U266 and RPMI8226 cells.
Caspase 9.
Caspase 9 is a marker for activation of the mitochondrial intrinsic apoptotic pathway. All four cell lines exhibit activation of caspase 9, as expected, based on the mechanism of action of bortezomib. Similar to caspase 8, the expression of caspase 9 activation was higher in MM.1S and NCI-H929 cells relative to U266 and RPMI8226 cells. MM.1S and NCI-H929 cells exhibited a similar magnitude of relative expression (about 10-fold to 15-fold from untreated-control cells). The stimulation peaked at 6 hours, and the relative expression was lost after 12–24 hours. The fold-change peak activation was close to 8-fold to 10-fold in U266 cells (at 12 hours) and 2-fold to 3-fold in RPMI8226 cells (at 6–12 hours). The change in relative expression from untreated-control cells was lost in these cells lines after 24 hours (Supplemental Fig. 11).
Correlation between Cell Viability and Protein Expression
Spearman’s rank correlation was performed to test the relationship between cell viability under treatment with different protein expression signatures. Spearman’s rank correlation coefficients were calculated between temporal values of experimental cell viability and relative protein expression for each cell line individually (Supplemental Fig. 12). An inverse correlation was observed between cell viability under treatment and pNFκB, pAKT, pJNK, and caspase 8 in all four cell lines. This was expected since treatment induced reduction in cell viability results from the activation of pJNK (stress) and caspase 8. The reduced viability might also activate survival mechanisms in cells to combat cell death and attempt cell growth in the presence of a chemotherapeutic agent by way of the activation of survival pathways (e.g., pNFκB and pAKT). A positive correlation was found between cell viability and Bcl-2 and pBAD across cell lines, as anticipated, since cell death is associated with the inhibition of these proteins. The overall correlation matrix signature (correlation between cell viability-protein expression, and protein-protein expressions) was similar across cell lines. Visual inspection suggests that the matrix was very similar for MM.1S and NCI-H929 cells, suggesting similarities in protein activation upon bortezomib treatment. The matrix signatures were not as comparable between U266 and RPMI8226 cells, especially because of differences in the strength of correlations and not the directionality per se. Finer hierarchical clustering analysis of the four correlation matrices to detect similarities in the tendency of the proteins to group together based on the Euclidean distance measure was also tested (Supplemental Fig. 13). Overall, fewer clustering similarities were noted between cell lines. Broadly, proteins clustered as upstream (pNFκB, pAKT, pJNK, pmTOR, and pp53) and downstream (pBAD, Bcl-2, p21, and caspase 9) proteins, with the exception of caspase 8 in all cell lines.
Discussion
MM is a hematological cancer known to present interclonal and intraclonal heterogeneity among cancer cells, leading to interindividual and intraindividual diversity in response to treatment regimens (Keats et al., 2012; Bolli et al., 2014). The objective of this study is to explore the differences in sensitivity to bortezomib treatment across four myeloma cell lines using mathematical model–based approaches.
A standard inhibitory effect function (eq. 1) was used as the first approach to investigate differences in the sensitivity of cell lines toward bortezomib. MM.1S cells exhibited the lowest estimated IC50 values (Fig. 1; Table 1) at all time points. These results compared well with a similar study by Shabaneh et al. (2013), who reported lower IC50 values (at 48 hours of bortezomib treatment) in MM.1S (2.3 nM) and NCI-H929 (1.9 nM) cells when compared with RPMI8226 cells (5.9 nM). Co-modeling was used to statistically test differences in potency (Table 2), supporting the conclusion that MM.1S cells appear to represent a more sensitive cell line compared with U266, RPMI8226, and NCI-H929 cells.
The temporal profiles of cell proliferation identified two groups of cell lines based on sensitivity at higher bortezomib concentrations: more sensitive cell lines (MM.1S and NCI-H929); and less sensitive cell lines (U266 and RPMI8226) (Supplemental Fig. 1). This difference in bortezomib sensitivity among cell lines was shown to be statistically significant. Shabaneh et al. (2013) also showed increased bortezomib sensitivity of MM.1S and NCI-H929 cell lines in comparison with RPMI8226 based on the measured percentage of apoptotic cells. Mathematical modeling suggested exponential growth in all four untreated cells with a second-order saturable cytotoxicity function to describe the bortezomib treatment data. The model also included the in vitro degradation kinetics of bortezomib in media since the duration of the study was up to 96 hours (Chudasama et al., 2015). The model fits (Fig. 2) estimated similar values for the exponential kg values for the four cell lines (Table 2), in accordance with prior estimates of their doubling times (U266, 55–70 hours; RPMI8226, 60–70 hours; MM.1S, 72 hours; NCI-H929, 70 hours) (Greenstein et al., 2003; https://www.dsmz.de/catalogues/details/culture/). At high drug concentrations, Kmax governs cell death, and greater Kmax values were estimated for more sensitive MM.1S and NCI-H929 cell lines. At concentrations lower than the KC50, the Kmax/KC50 ratio governs cell death, and similar values for the ratio across cell lines suggest little to no difference across the cell lines at 2 and 4 nM concentrations.
The literature suggests that the ratio of proteasome load/capacity may determine differences in apoptotic sensitivity (Bianchi et al., 2009). The rank ordered measured values of proteasome activity (which reflects proteasome capacity), and proteasomal degradation (which reflects proteasome load) was reported by Bianchi et al. (2009) in four MM cell lines (U266, RPMI8226, KMS.18, and MM.1S). These data were digitized to calculate proteasome load/capacity values, and a comparison between load/capacity and IC50 values at 24, 48, and 72 hours is shown in Supplemental Fig. 14, a–c. Although the data are sparse, the comparison shows that RPMI8226 and U266 cells behave differently from MM.1S cells, in agreement with the findings of Bianchi et al. (2009). Shabaneh et al. (2013) also experimentally derived proteasome load/capacity through functional assays and obtained a reasonable correlation coefficient (r2) of 0.74 with seven MM cell lines. After digitizing the proteasome load/capacity values and comparing with IC50 values (Supplemental Fig. 14, d–f), NCI-H929 and MM.1S cells separate from RPMI8226 cells, albeit with poor correlation coefficients. The analysis suggests the need for more research to determine whether this ratio represents a suitable biomarker of variability in responses to bortezomib treatment.
To test the hypothesis that protein-signaling differences are associated with the heterogeneous bortezomib pharmacodynamic responses, a logic-based network–modeling approach was employed to identify critical intracellular signaling proteins stimulated by bortezomib in myeloma cells. Dynamic simulations of our comprehensive network model identified several interesting mechanisms of bortezomib action. The most noteworthy of these was the activation of the extrinsic apoptotic pathway mediated by the stimulation of caspase 8. This stimulation of the extrinsic pathway is primarily directed by the stimulation of JNK by bortezomib, which in turn activates caspase 8 via the death receptors (Mitsiades et al., 2002; Qin et al., 2005; Kubiczkova et al., 2014). Bortezomib-mediated cell death results from both intrinsic (mitochondrial regulation of apoptosis; biomarker, caspase 9) and extrinsic (death receptor–mediated signaling; biomarker, caspase 8) pathways of apoptosis (Mitsiades et al., 2002; Kubiczkova et al., 2014). The stimulation of the prosurvival AKT pathway by bortezomib is also consistent with reports (Que et al., 2012; Mimura et al., 2014). The counterintuitive stimulation of prosurvival pathways suggests cancer cells push survival pathways into overdrive to combat stress conditions generated in the cells.
The arbitrary time units in the network simulations make assessment of temporal patterns of cellular processes difficult, making a case for longitudinal measurements of cellular processes/proteins along with continuous ordinary differential equation–based pharmacodynamic models. Hence, we used the network approach to qualitatively identify critical proteins mediating bortezomib pharmacodynamics via model simulations, topology-based centrality measures (Table 2), and a reduction algorithm (Fig. 3; Table 3). The investigation of graded protein responses will require the development of a quantitative model, which is not within the scope of qualitative Boolean network models.
Network-guided experimental snapshot of protein dynamics showed that cell lines MM.1S and NCI-H929 exhibit a greater magnitude of stimulation/inhibition of proteins compared with U266 and RPMI8226 cells. This observation was consistent at the 2 nM treatment, with MM.1S cells exhibiting stimulation/inhibition of proteins, whereas no change was observed in RPMI8226 cells. MM.1S and NCI-H929 cells also showed an earlier onset of activation of the key proteins. Protein dynamics alongside of an in vitro pharmacodynamic response highlighted a potential association between increased sensitivity to bortezomib treatment with earlier onset of activation and a greater magnitude of change in protein expression (Fig. 4).
The magnitude of stimulation of the prosurvival protein pNFκB (Fig. 5) was greater (150-fold) in the more sensitive MM.1S cells compared with all the other cell lines. Between NCI-H929 and U266 cells, the less sensitive U266 cells showed greater activation, although both cell lines displayed a considerably greater stimulation from control cells when compared against RPMI8226 cells. In the case of pAKT and pmTOR (Supplemental Figs. 6 and 7)—also cell proliferation associated proteins—little difference could be discerned in the protein profiles, making a corresponding association with drug sensitivity difficult. For the stress-associated protein pJNK (Supplemental Fig. 8), U266 cells showed the maximum fold change in induction, whereas all the other cell lines displayed similar magnitudes of activation.
Sensitivity to bortezomib appeared to be better associated with apoptosis-related proteins. The proteins belonging to the Bcl-2 family, Bcl-2 and pBAD (Fig. 6; Supplemental Fig. 10), showed a faster rate and stronger (near complete inhibition) suppression in the more sensitive MM.1S and NCI-H929 cells, whereas U266 and RPMI8226 cells presented with a slower rate with only 50% inhibition. Similar dynamics were reflected in the protein expression profiles of caspase 8 and caspase 9 (Fig. 8; Supplemental Fig. 11). MM.1S cells displayed greater stimulation (150-fold) of caspase 8 (similar to that of pNFκB), and NCI-H929 cells exhibited moderately greater stimulation (20-fold to 25-fold), both of which were greater than the stimulation in U266 and RPMI8226 cells. Caspase 9 demonstrated similar trends among cell lines. The gradual loss of signal for caspases 8 and 9 in all cell lines over time could depend on the stability of these proteins after apoptosis is triggered in cells and in the increasingly acidic media containing dead cells (Stennicke and Salvesen, 1997; MacKenzie and Clark, 2012).
Heterogeneity in the temporal dynamics of pp53 (tumor suppressor) was observed across all cell lines (Fig. 7). No change in the expression of pp53 was observed in the more sensitive cell lines MM.1S and NCI-H929 (carrying a wild type, TP53), whereas stimulation was observed in U266 and RPMI8826 cells (carrying a mutated TP53) (Pichiorri et al., 2010). The possible reasons for this observation could be: 1) the presence of MDM2 associated inhibition of p53 expression (Pichiorri et al., 2010; Bi and Chng, 2014; Teoh et al., 2014); and/or 2) the regulation of p53 by microRNA (miR; miR-25 and miR-30d) (Kumar et al., 2011). The dynamics of p21 (Supplemental Fig. 9), a cell-cycle inhibitor, were unclear and were not amenable to speculating on the differences between cell lines. The lack of clarity in the dynamics could be due to the short half-life of p21 (minutes to a few hours) in cancer cells (Gareau et al., 2011; Scoumanne et al., 2011). Additionally, statistical correlation analysis also revealed similarities between the two groups of cells (more sensitive/less sensitive) in terms of temporal protein–protein associations, along with similarities in cell viability–protein expression (Supplemental Fig. 12).
This research suggests an association between sensitivity differences with intracellular signaling protein dynamics triggered upon bortezomib exposure to cells. However, there are some limitations to the current study. The number of cell lines tested in this study is relatively small for establishing correlations between sensitivity and protein biomarker dynamics. Some of the protein profiles were associated with considerable variability among replicates, clouding the certainty in terms of magnitude and directionality of protein stimulation, especially p21 dynamics. A panel of more cell lines and replicates would add to this work to better assess apparent associations. In addition, genetically modified MM cell lines could be used to test and validate the network model. Nevertheless, this proof-of-concept study suggests that intracellular signaling heterogeneity represents a possible source of pharmacodynamic heterogeneity and highlights the need to monitor heterogeneity by identifying critical protein biomarkers with the potential to track therapeutic responses.
Authorship Contributions
Participated in research design: Ramakrishnan and Mager.
Conducted experiments: Ramakrishnan.
Performed data analysis: Ramakrishnan and Mager.
Wrote or contributed to the writing of the manuscript: Ramakrishnan and Mager.
Footnotes
- Received January 18, 2018.
- Accepted April 5, 2018.
This work was supported by the National Institutes of Health [Grant GM57980].
↵This article has supplemental material available at jpet.aspetjournals.org.
Abbreviations
- AKT
- protein kinase B
- ATCC
- American Type Culture Collection
- BAD
- B-cell lymphoma 2–associated death promoter
- Bcl-2
- B-cell lymphoma 2
- CDK
- cyclin-dependent kinase
- Chk
- checkpoint kinase
- DMSO
- dimethylsulfoxide
- γ
- Hill coefficient
- Imax
- maximum inhibition of cell viability
- IGF-1
- insulin-like growth factor type 1
- IRF4
- interferon regulatory factor 4
- IRS
- insulin receptor substrate
- JAK
- janus kinase
- JNK
- c-Jun N-terminal kinase
- kg
- growth rate constant
- KC50
- Michaelis-Menton constant equivalent to the concentration of bortezomib that results in half-maximal maximum killing rate constant
- Kmax
- maximum killing rate constant
- miR
- microRNA
- MM
- multiple myeloma
- mTOR
- mechanistic target of rapamycin
- NFκB
- nuclear factor-κB
- PI3K
- phosphoinositide 3-kinase
- pp53
- phosphorylated p53
- R0
- baseline response in the absence of bortezomib
- RAS
- renin-angiotensin system
- ROS
- reactive oxygen species
- SM
- smoldering myeloma
- STAT
- signal transducers and activators of transcription
- Copyright © 2018 by The American Society for Pharmacology and Experimental Therapeutics