Systems models of biological networks show promise for informing drug target selection/qualification, identifying lead compounds and factors regulating disease progression, rationalizing combinatorial regimens, and explaining sources of intersubject variability and adverse drug reactions. However, most models of biological systems are qualitative and are not easily coupled with dynamical models of drug exposure-response relationships. In this proof-of-concept study, logic-based modeling of signal transduction pathways in U266 multiple myeloma (MM) cells is used to guide the development of a simple dynamical model linking bortezomib exposure to cellular outcomes. Bortezomib is a commonly used first-line agent in MM treatment; however, knowledge of the signal transduction pathways regulating bortezomib-mediated cell cytotoxicity is incomplete. A Boolean network model of 66 nodes was constructed that includes major survival and apoptotic pathways and was updated using responses to several chemical probes. Simulated responses to bortezomib were in good agreement with experimental data, and a reduction algorithm was used to identify key signaling proteins. Bortezomib-mediated apoptosis was not associated with suppression of nuclear factor κB (NFκB) protein inhibition in this cell line, which contradicts a major hypothesis of bortezomib pharmacodynamics. A pharmacodynamic model was developed that included three critical proteins (phospho-NFκB, BclxL, and cleaved poly (ADP ribose) polymerase). Model-fitted protein dynamics and cell proliferation profiles agreed with experimental data, and the model-predicted IC50 (3.5 nM) is comparable to the experimental value (1.5 nM). The cell-based pharmacodynamic model successfully links bortezomib exposure to MM cellular proliferation via protein dynamics, and this model may show utility in exploring bortezomib-based combination regimens.
The fields of systems biology and pharmacokinetic (PK)/pharmacodynamic (PD) modeling have evolved largely in parallel, and there is an emerging consensus that an effective integration of these disciplines is needed to fully realize the promise of each in bringing new therapeutic molecules and combination regimens to the bedside (http://www.nigms.nih.gov/Training/Documents/SystemsPharmaWPSorger2011.pdf). Traditional PK/PD models of drug action use compartmental structures to integrate the time course of drug exposure, pharmacological properties (capacity, sensitivity, and transduction of drug-target interactions), and (patho)physiological turnover processes (Mager et al., 2003). Such semimechanistic models contain a minimal number of identifiable parameters to describe temporal profiles of macroscale therapeutic and adverse drug effects. When coupled with nonlinear mixed-effects modeling of relatively large clinical trials, a covariate analysis can be used to identify patient-specific characteristics (e.g., genetic polymorphisms) that explain the interindividual variability in model parameters (Pillai et al., 2005). Although a major component of model-informed drug development and therapeutics (Milligan et al., 2013), this approach can be limited by specific study designs and is rarely sufficient for recapitulating multiple, complex genotype-phenotype relationships.
Significant insights have been realized from the recognition that both drugs and disease processes give rise to complex and dynamic clinical phenotypes by altering natural interconnected biochemical networks and support the emergence of systems pharmacology models of drug action (Zhao and Iyengar, 2012; Huang et al., 2013; Jusko, 2013). Multiscale models that combine PK/PD principles and signaling networks can serve as a platform for integrating genomic/proteomic factors that regulate drug effects and clinical outcomes—so-called enhanced PD (ePD) models (Iyengar et al., 2012). Two major challenges to the development of ePD models include the lack of complete mathematical models of signal transduction networks (e.g., concentrations and reaction rate constants) and unknown quantitative relationships between individual genomic/proteomic differences and model parameters. Mechanistic network models are preferred over empirical structures (Birtwistle et al., 2013); however, such models may not be defined or calibrated to specific pharmacological and/or disease systems. Logic-based modeling techniques provide a global perspective of system properties in the absence of kinetic parameters through the integration of qualitative a priori knowledge of network connections (Albert and Wang, 2009). In this proof-of-concept study, a reduction algorithm is applied to a mathematical network to guide the development of a small signaling model for bortezomib, a potent proteasome inhibitor, which may ultimately serve as an ePD model for its use in multiple myeloma (MM).
Multiple myeloma is a B cell neoplasm associated with several comorbidities, including hypercalcemia, renal insufficiency, anemia, and bone lesions (Caers et al., 2008; Blade et al., 2010). The prognosis for advanced stages of MM is poor despite multiple treatment options, with a median survival of advanced-stage patients of less than 10 months (Richardson et al., 2003, 2005). Bortezomib is commonly prescribed alone or in combination with other antimyeloma agents (Oancea et al., 2004), and the addition of bortezomib has significantly improved overall survival of MM patients (Caers et al., 2008). However, almost all patients relapse and become refractory to all treatment options (Richardson et al., 2005; San Miguel et al., 2008), and a better understanding of MM disease progression and mechanisms of drug action is critical for improving the treatment of MM.
Bortezomib modulates both survival and apoptotic cellular pathways in MM cells (Fig. 1), and multiple mechanisms of action are proposed for inducing cell death, including the inhibition of proteasome and the phospho–nuclear factor κB (pNFκB) pathway (Hideshima et al., 2001, 2002, 2003a,b). However, recent reports suggest that bortezomib stimulates receptor-interacting protein (RIP; a signaling protein upstream of the NFκB pathway), leading to activation of pNFκB protein expression (Hideshima et al., 2009). Despite extensive qualitative information on bortezomib-induced intracellular signaling, mathematical models linking bortezomib exposure to intracellular protein dynamics have not been established. Systems-level modeling may facilitate the rational design of single and combinatorial dosing regimens, approaches to overcome drug resistance, and prevent suboptimal dosing (Berger and Iyengar, 2011; Zhao and Iyengar, 2012). Mechanistic models are emerging in which chemotherapy exposure is connected to ultimate responses through target occupancy and biomarker signal transduction (Yamazaki et al., 2011; Harrold et al., 2012; Kay et al., 2012; Kirouac et al., 2013). The purpose of this study is to integrate network systems analysis and PK/PD modeling principles to study bortezomib effects on signal transduction in MM cells. Although the current study is conducted in U266 cells, the strategic framework may be extended to other cell lines and applied to other therapeutic areas.
Materials and Methods
Tissue Culture Materials.
The human myeloma cell line U266 (TIB-196) was purchased from American Type Culture Collection (Manassas, VA). Inhibitor of κB kinase inhibitor (IKKi) (PS-1145) and Janus kinase (JAK) inhibitor I (JAKi) were purchased from Santa Cruz Biotechnology (Dallas, TX). IκBi (BAY11-7085) was purchased from Sigma-Aldrich (St. Louis, MO). The clinically available formulation of bortezomib was used for all experiments (Millennium Pharmaceuticals, Cambridge, MA). The primary monoclonal antibodies to phospho-p65, phosphor-IκBα (pIκBα), phospho–signal transducers and activator of transcription 3 (pStat3), Bcl-xL, poly (ADP ribose) polymerase (PARP), p65, Stat3, IκBα, and β-actin were purchased from Cell Signaling Technology (Danvers, MA). α-Tubulin primary and rabbit and mouse secondary antibodies were purchased from Santa Cruz Biotechnology. WST-1 reagent assay kit was purchased from Roche (Basel, Switzerland). Horseradish peroxidase conjugate was purchased from Bio-Rad (Hercules, CA), and enhanced chemiluminescence substrate was purchased from Thermo Scientific (Pittsburgh, PA). RPMI 1640 and fetal bovine serum (FBS) were purchased from American Type Culture Collection.
Cell Culture Experimental Design.
All experiments were conducted in U266 myeloma cells. Cells were cultured in RPMI 1640 medium supplemented with 15% FBS and 1% penicillin/streptomycin antibiotics. Treatment protocols included 1) IKKi (10 µM), 2) JAKi (10 µM), 3) IκBi (10 µM), 4) IKKi and JAKi (10 µM each), 5) IκBi and JAKi (10 µM each), 6) bortezomib (20 and 2 nM), and 7) dimethylsulfoxide vehicle control. All treatments were continuous for 48 hours, with the exception of IκBi experiments (10 hours). A transient exposure to bortezomib (20 nM) was also conducted, in which cells were incubated for either 1 or 2 hours followed by removal of drug from the culture media and subsequent incubation with vehicle control for up to 48 hours.
WST-1 Cell Proliferation Assay.
Approximately 10,000 cells per well were incubated with a range of concentrations (0.001–100 µM) of JAKi, IKKi, or a combination of IKKi and JAKi. Incubation experiments with bortezomib alone ranged from 0.001 to 100 nM. Cell viability was measured at 0, 24, 48, 72, and 96 hours according to manufacturer instructions. Absorbance was measured at 450 and 690 nm after 2-hour incubation with WST-1 reagent with a Microplate spectrophotometer (Molecular Devices, Sunnydale, CA).
Western Blot Analysis.
Relative protein expression levels of pStat3, pNFκB, pIκBα, cleaved PARP, and BclxL were measured following treatment regimens as specified in Supplemental Table 1. Cells (5 × 106/10 ml culture media) were incubated in 10-cm2 culture dishes. Cells were collected at the end of the treatment duration, washed with ice-cold phosphate-buffered saline, and lysed at 4°C in radioimmunoprecipitation assay RIPA buffer (Cell Signaling Technology) supplemented with protease and phosphatase inhibitor cocktail (Bio-Rad) and phenylmethylsulfonyl fluoride (Bio-Rad). Protein lysates were stored at −80°C until used. Equal amounts of proteins were separated on SDS-PAGE gels and transferred to nitrocellulose membranes. Immunoblotting was performed according to manufacturer instructions, and the relative intensity of bands was assessed by densitometric analysis of digitized images using ImageJ software (NIH, Bethesda, MD). All experiments were conducted in at least duplicate.
Bortezomib degradation in cell culture was assessed at different time points (0, 12, 24, 48, and 72 hours). Bortezomib concentrations were measured after incubation (20 nM) at 37°C in culture medium (RPMI 1640 supplemented with 15% FBS and 1% penicillin/streptomycin) using a liquid chromatography–tandem mass spectrometry assay previously validated in our laboratory (unpublished). Bortezomib degradation kinetics were fit with a monoexponential (λ) decay model, and the half-life was calculated as ln(2)/λ.
Boolean Network Assembly and Simulations.
The Boolean network model (Supplemental Fig. 1) was drawn with SmartDraw software (http://www.smartdraw.com) and was implemented in Odefy (http://www.helmholtz-muenchen.de/cmb/odefy) (Wittmann et al., 2009), a toolbox compatible with MATLAB (MathWorks, Natick, MA). Odefy converts Boolean relationships into a continuous framework using a Hill-type function. In the hypothetical case in which node B stimulates A, the relationship can be described byDefault parameter values for k (0.5), n (3), and τ (1) were used for all nodes. All model simulations were performed using the “HillCubeNorm” option within Odefy (MATLAB model code is provided in the Data Supplement). Each node in the model was initialized as ON (1) or OFF (0) based on baseline activity in the cell. For example, if the state is constitutively active in the cell, it was initialized as 1. Boolean transfer functions and initial state values are summarized in Supplemental Table 2. Examples of ordinary differential equations corresponding to each node, which were generated using Odefy, are listed in Supplemental Eqs. 1–5.
Three specific probes, JAKi, IκBi, and IKKi, were used to develop and update the network model. The effect of the JAKi was added such that it inhibits nodes JAK1 and JAK2. The IκBi inhibitor was added by direct inhibition of pIκB node activation. Similarly, the IKKi effect was added such that it inhibits node IKK. Model simulations were performed to predict bortezomib outcomes, and drug effect was incorporated as direct inhibition of the proteasome node and stimulation of the RIP node. For simulations of transient bortezomib exposure, the bortezomib node was activated for only one or two iterations.
The final full Boolean network model developed for U266 cells (Supplemental Fig. 1) was reduced using a model reduction algorithm (Veliz-Cuba, 2011) to identify critical proteins involved in bortezomib-mediated signal transduction. Figure 2 shows examples of the steps taken to reduce the network. First, nodes with only one input and one output are eliminated and the pathways are reconnected (Fig. 2A). For example, RAF was removed and RAS and MEK1 were directly connected (Fig. 2B). In another network reduction rule, nonfunctional edges are identified and removed based on the Boolean relationships (algorithm S from Veliz-Cuba, 2011). For example, Fig. 2C shows nodes A, B, and C, with the following Boolean functions: f(A) = B OR (B AND NOT C), f(B) = (A AND C), and f(C) = A. Hence, node C is nonfunctional in f(A) because of the “OR” relationship, and the B node is sufficient to elicit the same effect without invoking C. The resulting diagram is shown (block arrow, Fig. 2C) with the Boolean relationships: f(A) = B, f(B) = (A AND C), and f(C) = A. p21 is an example of nonfunctional edges (Fig. 2D), and Boolean relationships for p21 [f(p21) = p53 AND (NOT AKT OR NOT MDM OR NOT MYC OR NOT CKD4)] were refined using step 2 to f(p21) = p53 AND NOT AKT to remove nonfunctional edges. AKT is the sole connection to BAD; therefore, the p21 connection to AKT is retained (Fig. 2D).
Finally, nonfunctional nodes are identified and removed from the network (algorithm R from Veliz-Cuba, 2011). Before applying this rule, the network is rewired and the first rule is reapplied. Nodes that do not have any impact on other nodes are considered nonfunctional and are removed. In Fig. 2E, node C is regulating node B, but at the same time node A is regulating node C. Therefore, the connection from node C to B can be removed, and consequently node C can be removed. As an example, the only node modulated by MYC is CYCE, as the MYC connection to CYCD can be removed because of the AKT connection to CYCD (Fig. 2F). pStat3, mitogen-activated protein kinase, and extracellular signal-regulated kinase stimulate MYC independently, and mitogen-activated protein kinase was removed following the first network reduction rule (Fig. 2A). Similarly, extracellular signal-regulated kinase was removed as it does not modulate any other node besides MYC. pStat3 modulates several other nodes (e.g., BclxL and XIAP) and stimulates MYC. MYC can be removed following the one input–one output rule, and therefore, CYCE is directly modulated by pStat3 (Fig. 2F), resulting in four final nodes. This entire process was repeated until no further nodes could be removed from the Boolean network. Boolean simulations of the final reduced network were compared with simulations of the full network using Odefy (Wittmann et al., 2009).
Dynamical Model Development.
Several critical proteins identified using the model reduction algorithm were incorporated into a reduced cellular PD model based on mechanisms of bortezomib action and signal transduction pathways. The reduced dynamic model of bortezomib is depicted in Fig. 3. In brief, bortezomib stimulates protein expression of pNFκB and cleaved PARP (following delays) and inhibits expression of BclxL. A transit compartment (M1B) was added to account for the slight delay in stimulation of pNFκB protein expression (Mager and Jusko, 2001). Bortezomib stimulates the synthesis rate of M1B (ksyn_M1B) via a linear stimulation coefficient (Sm_M1B), and the rate of change in expression of M1B is defined as follows:(1)where CB is the concentration of bortezomib, and kdeg_M1B is a first-order degradation rate constant of M1B. pNFκB protein expression increases but returns toward the baseline after 24 hours of continuous bortezomib exposure. To capture this trend, a precursor-dependent indirect response model was proposed (Sharma et al., 1998), with stimulation of the first-order transfer rate constant from the precursor pool (pre_pNFκB) to the pNFκB compartment (ktr_NFκB). The rate of change of pre_pNFκB and pNFκB is described by the following equations:(2)(3)where ksyn_NFκB is a zero-order production rate constant of pre_pNFκB, and kdeg_NFκB is a first-order degradation rate constant of pNFκB. Changes in cleaved PARP expression follow a significant delay, and the protein expression profile was best characterized using a time-dependent transduction model (Mager and Jusko, 2001), with negative feedback from the last compartment to the synthesis rate of first cleaved PARP compartment to characterize the downward phase of protein expression. Bortezomib stimulates the synthesis rate (ksyn_PARP) of the first cleaved PARP compartment (cPARP1) with a linear stimulation coefficient (Sm_PARP):(4).The three subsequent transit compartments are described by the following general equation:(5)with n representing the compartment number (n = 2–4). pNFκB is a prosurvival transcription factor that upregulates antiapoptotic proteins (e.g., BclxL), whereas PARP is a proapoptotic factor that results in apoptosis when activated (cleaved PARP) and stimulates the degradation of BclxL. The rate of change of BclxL is defined as follows:(6)where ksyn_Bcl is a zero-order synthesis rate constant of BclxL, kdeg_Bcl is a first-order degradation rate constant of BclxL, and λ is a power coefficient for the cleaved PARP effect on BclxL. Cell proliferation is dependent on the balance between apoptotic and antiapoptotic signals. Therefore, cell proliferation (N) was modulated by BclxL and cleaved PARP expression profiles, with the first-order natural cell death rate constant (kd) inhibited by BclxL and stimulated by cleaved PARP. N represents cell proliferation:(7)with kg representing a first-order growth rate constant. Secondary equations defining zero-order production rate constants are listed in Supplemental Eqs. 6–9.
The dynamical model (Fig. 3) was first fitted to mean pNFκB and cleaved PARP data, followed by pooling of all relative changes in these protein expression patterns (naïve pooled data approach). Initially, pNFκB and cleaved PARP protein dynamics were fit separately. BclxL data were subsequently included, and all of the data including BclxL, pNFκB, and cleaved PARP protein dynamics were fitted simultaneously. Next, cell proliferation was introduced to the existing cellular dynamic model, and all of the data (three biomarkers and cell proliferation) were fitted simultaneously. The in vitro degradation half-life of bortezomib was included for all model runs. Parameters were estimated using MATLAB (fminsearch function, maximum likelihood algorithm, and ode23s) and a model development framework (Harrold and Abraham, 2014). All protein dynamics were described using a proportional error variance model (), and cell proliferation was fitted using an additive plus proportional error variance model (), where Yobs is the observation at time t, Ypred is the model-predicted value at time t, and σ and ε are estimated variance model parameters.
The final cell-based model and parameter estimates were used to simulate protein dynamics and cell proliferation after continuous exposure to a 10-fold-lower bortezomib concentration (2 nM). Simulations were compared with the experimentally measured cell proliferation and protein expression profiles. Only parameters associated with natural cell proliferation (kg and kd) were estimated. Cell proliferation was measured at 0, 24, 48, and 72 hours, and protein expression measurements were made at 0, 1, 4, 6, 8, 11, 24, 33, and 48 hours. The final model was also used to predict cell proliferation at 48 hours for a range of bortezomib concentrations (0.001–100 nM). The experimentally estimated IC50 value was compared with the IC50 obtained from the simulated data.
Network Development with IKK and JAK Inhibitors.
The final Boolean network model incorporates major survival and apoptotic pathways in U266 cells (Supplemental Fig. 1), and Supplemental Table 2 summarizes all node descriptions, Boolean internodal relationships, and initial values. The green connections in Supplemental Fig. 1 indicate how the initial network was updated based on cellular responses to the pathway probes. Nodes that are constitutively active under baseline cell conditions were set as being active (1) or otherwise inactive (0) in the network. An initial Boolean network was constructed, and the performance of the initial model was tested with two pathway-specific probes, IKKi and JAKi, to inhibit the NFκB and JAK/STAT pathways. In the Boolean network, the effect of each inhibitor was added such that IKKi inhibits node IKK, and JAKi inhibits nodes JAK1 and JAK2. Simulations of the initial Boolean network show that nodes representing pIκBα, pNFκB, pStat3, and BclxL expression decrease gradually upon inhibition of node IKK (Supplemental Fig. 2A, black dash-dotted lines). Inhibition of pStat3 expression appears to be due to interleukin-6 inhibition (major stimulus of JAK/STAT3 pathway) via pNFκB expression. However, in vitro Western blot experiments show that expressions of pNFκB, pStat3, and BclxL remained unchanged after 48-hour exposure to the IKKi, whereas pIκBα expression decreased (Supplemental Fig. 2B, black symbols and dotted lines). In contrast, simulations of the same proteins using the initial model following JAKi exposure (Supplemental Fig. 2A, green dashed lines) were comparable to observed experimental profiles over 48 hours (Supplemental Fig. 2B, green symbols and dashed lines). To account for the observed trend in pNFκB expression, the initial network was modified such that pStat3 also stimulates pNFκB (Supplemental Fig. 1). Although not yet confirmed in U266 or other MM cell lines, studies suggest cross-talk between NFκB and STAT3 signaling in tumors (Squarize et al., 2006; Lee et al., 2009; Saez-Rodriguez et al., 2011). To further test this hypothesis, model simulations were compared for the combination of IKKi and JAKi treatment with measured temporal expression profiles of pNFκB. Simulations of pNFκB following combination treatment with JAKi and IKKi using the initial Boolean network showed a steady decrease over time (data not shown). However, pNFκB expression measured by Western blot analysis only transiently decreased with a return to baseline values (Supplemental Fig. 2B, blue symbols and dash-dot lines). A detailed evaluation of pNFκB is described in the next section. Overall, cellular protein dynamic simulations using the final model (Supplemental Fig. 2C) reasonably agree with the experimental results (Supplemental Fig. 2B).
As the expression of pNFκB was not suppressed by the IKKi, either as a single agent or in combination with the JAKi, cellular responses to another NFκB pathway inhibitor (IκBi) were evaluated alone or in combination with the JAKi. Using the initial Boolean network, inhibition of node pIκBα by IκBi resulted in simulations showing inhibition of pNFκB expression (data not shown). However, continuous IκBi exposure to U266 cells for up to 10 hours (alone and in combination with JAKi) maintained pIκBα expression below baseline values (Supplemental Fig. 3A), whereas pNFκB expression transiently decreased followed by a steady increase above the baseline (Supplemental Fig. 3A). Based on the lack of pNFκB suppression with either the IKKi or IκBi, alone or in combination with JAKi, additional factors were assumed to govern pNFκB expression in U266 cells outside of IKK, pIκBα, and pStat3 signaling. This phenomenon was emulated by incorporating a dummy node “X” to constantly stimulate pNFκB in the final network (Supplemental Fig. 1). The addition of this factor stabilized the Boolean network and reconciled simulations (Supplemental Fig. 2, C, D, and F) with all experimental data (Supplemental Fig. 2, B, E, and G). These results highlight the need to inhibit NFκB activity directly rather than through upstream pathways. Both final model simulations and experimental data confirm inhibition of pIκBα expression following continuous 48-hour exposure to the IKKi (10 µM), as well as the lack of suppression of pNFκB and BclxL expression levels (Supplemental Fig. 2, B and C, black symbols and dotted lines). Analogously, continuous exposure to the JAKi (10 µM) suppressed expression of pStat3, whereas pNFκB and BclxL expression levels remained unchanged (Supplemental Fig. 2, B and C, green symbols and dashed lines). In addition, final model simulations of U266 cellular proliferation and apoptosis agreed well with experimental cellular responses (Supplemental Fig. 2, D–G). A range of IKKi and JAKi concentrations (0.001–100 µM) was used to evaluate cell viability, and neither induction of apoptosis nor decreased cell viability was observed after 72-hour continuous exposure alone or in combination (Supplemental Fig. 2, E and G).
Once the final Boolean network model was updated using protein dynamics and cellular responses after exposures to probe inhibitors, the model was used to evaluate bortezomib pharmacodynamics. Bortezomib effect was introduced into the network such that it directly inhibits the “proteasome” node and stimulates the “RIP” node (Supplemental Fig. 1). Stimulation of RIP activates a downstream cascade, leading to activation of pIκBα followed by pNFκB. In the Boolean network, a state initialized as 1 (active) cannot be further stimulated; therefore, pIκBα and pNFκB signal intensities remain unchanged (Supplemental Fig. 4A, blue dash-dot line). However, our immunoblot analysis revealed transient increases in pNFκB and pIκBα expression levels, followed by a gradual return to baseline values (Supplemental Fig. 4B, blue symbols and lines), which is in agreement with published profiles (Hideshima et al., 2009). Comparing steady-state values for pNFκB and pIκBα before and after treatment suggests no change in expression levels, which is in agreement with Boolean network simulations (Supplemental Fig. 4A). Furthermore, inhibition of proteasome results in accumulated cellular stress that activates the c-Jun N-terminal kinase (JNK) pathway as well as the apoptotic pathway. The mitochondrial apoptotic pathway activates caspase-3, leading to activation of p53 followed by p21, resulting in the inhibition of the Bcl-2 family proteins (e.g., BclxL and Bcl-2) (Supplemental Fig. 4B, blue symbols and lines). The activated JNK pathway inversely regulates pStat3, leading to downregulation of pStat3 expression (Supplemental Fig. 4B, blue symbols and lines). Simulations using the final Boolean network model correctly predicted the inhibition of pStat3 and BclxL in U266 cells after bortezomib exposure (Supplemental Fig. 4A, blue symbols and lines). In addition, final model simulations of U266 cellular outcomes to bortezomib exposure were in good agreement with in vitro measurements of cell viability and apoptosis measured via cleaved PARP expression (Fig. 4).
To further test the fidelity of the Boolean network model to predict bortezomib pharmacodynamics, we compared Boolean network simulations with experimental results of cellular responses upon transient bortezomib exposures. Simulations of transient drug exposure were achieved by maintaining the Bortezomib node as active (1 = ON) for a limited number of iterations, and the dynamics and steady states of protein and cellular response nodes were monitored. Transient experiments were conducted in which U266 cells were briefly exposed (i.e., for 1 and 2 hours) to bortezomib. Interestingly, some nodes (e.g., proteasome) returned to baseline values once bortezomib was removed from the system, whereas other nodes (e.g., cleaved PARP or Cl. PARP) achieved a new steady state even after drug was removed. Furthermore, steady-state outcomes for several nodes varied depending on the duration of simulated bortezomib exposure. For example, a relatively short duration of proteasome suppression produces a transient induction of apoptosis, whereas longer durations of bortezomib exposure resulted in steady states of apoptosis induction and inhibited cell growth that are similar to simulated outcomes following continuous drug exposure. A shorter duration of bortezomib exposure was associated with a delay in apoptosis induction as compared with continuous drug exposure (Supplemental Fig. 5A), and similar trends were observed experimentally (Supplemental Fig. 5B). Western blot analysis of BclxL and cleaved PARP expression shows a reduced magnitude in changes from baseline values between transient and continuous bortezomib exposure. The magnitude of cleaved PARP induction was reduced from 22- to 6-fold after 24 hours of bortezomib treatment. Model simulations of selected proteins were comparable with experimental results (Supplemental Fig. 5, A and B). However, stimulation of cleaved PARP is associated with apoptosis induction, downregulation of total NFκB expression, and subsequent suppression of pNFκB expression. In contrast, incubation of U266 cells with bortezomib (20 nM) for 2 hours failed to inhibit pNFκB expression (Supplemental Fig. 5C, bottom), despite activation of cleaved PARP for over 48 hours (Supplemental Fig. 5B, apoptosis, black triangles). A specific threshold of cleaved PARP activation might be required to inhibit NFκB expression, which is not incorporated into the current network model, and further studies are needed to test this hypothesis.
A Boolean network reduction approach (Veliz-Cuba, 2011) identified eight critical nodes (Supplemental Fig. 6A), five of which are survival pathway nodes (AKT, pNFκB, BclxL, proteasome, and pRB) and three are apoptotic pathway nodes (caspase-3, JNK, and p21). Boolean simulations of the reduced network are identical to the steady-state values of the full network, which further supports the reduction algorithm (Supplemental Fig. 6B). Although caspase-3 was identified as a critical protein through the Boolean network reduction, cleaved PARP expression was measured as a marker for apoptosis due to caspase-3 assay limitations.
Bortezomib Cellular Pharmacodynamic Model.
Three critical proteins (i.e., pNFκB, BclxL, and cleaved PARP) were selected from the reduced Boolean network model (Supplemental Fig. 6A) for measurement and inclusion in the cellular pharmacodynamic model. The final cellular model (Fig. 3) integrates the time courses of bortezomib exposure, protein dynamics, and cell proliferation. Bortezomib elicits its effects on pNFκB via stimulation of upstream proteins in the pNFκB pathway and on cleaved PARP via stress accumulation due to proteasome inhibition. Upon continuous exposure of U266 cells to bortezomib, pNFκB protein expression is stimulated after a slight delay with a return toward the baseline (Fig. 5A). The slight delay in the stimulation of pNFκB was well characterized by adding a simple transit compartment (M1B; Eqs. 1 and 2). pNFκB expression was well described by a precursor model, in which the delayed signal (M1B) stimulates the first-order transfer from the precursor to the pNFκB compartment (Eqs. 2 and 3). Parameters associated with M1B and pNFκB (kdeg_M1B and Sm_M1B) were not estimated with good precision; therefore, M1B-associated parameters were fixed to the estimates during an initial model run. This had no impact on model performance (i.e., model fits were not compromised and parameter values did not change substantially), but precision of the estimated parameters was significantly improved (Table 1). A relatively long delay was observed before cleaved PARP expression increased (Fig. 5B), and a transit-compartment model was selected to describe this delay. Four transit compartments were found to be adequate, and the total transit time was well estimated at 26.8 hours (Table 1, ktr_parp). BclxL expression starts decreasing after 12 hours (Fig. 5C), and this delay was also well characterized by the proposed model. The initial phase of BclxL was maintained at steady state by pNFκB, and the eventual decrease in BclxL expression was successfully described using cleaved PARP stimulation of the BclxL degradation rate constant. An exponential growth model well captured the natural cell proliferation under control conditions (Fig. 5D, triangles). The final signal transduction model reasonably captured the delay in cell death, and parameters were estimated with good precision (Table 1). Overall, the pharmacodynamic model well characterized the time courses of protein expression and cell proliferation profiles after exposure to bortezomib at 20 nM.
Sensitivity to Dose and Model Qualification.
The final pharmacodynamic model (Fig. 3) was developed using several cellular biomarkers after continuous exposure to a single, relatively high concentration of bortezomib (20 nM). Simulations of protein expression profiles and cell proliferation following exposure to different bortezomib concentrations (0.1–20 nM) were performed to evaluate the role of bortezomib dose in cellular outcomes (Supplemental Fig. 7). As bortezomib concentration is increased, the magnitude of pNFκB expression is also increased, and the time to peak response is decreased (Supplemental Fig. 7A). The extent of induced cleaved PARP is also increased (Supplemental Fig. 7B), along with greater suppression of BclxL expression (Supplemental Fig. 7C) as bortezomib concentration increases. For cell proliferation, cell growth is suppressed up to 2 nM (Supplemental Fig. 7D, red line), but greater concentrations result in cell death (Supplemental Fig. 7D). As an external predictive check, protein expression profiles and cell proliferation after low-dose bortezomib exposure (2 nM) were measured and compared with model simulations, which were obtained using the final model (Fig. 3) and parameter estimates (Table 1) from the high-dose bortezomib (20 nM) experiments. The low concentration of bortezomib (2 nM) was chosen, as it is close to its IC50 value (1.5 nM) from the in vitro exposure-response relationship at 48 hours. Thus, the selected concentration is relatively low but should still elicit a pharmacological response. Model-predicted pNFκB, cleaved PARP, and BclxL profiles and cell proliferation dynamics are in good agreement with the experimental data (Supplemental Fig. 8).
The final model was also used to predict the IC50 for bortezomib after 48 hours of continuous exposure. Although the model was developed based on temporal profiles after a single concentration of bortezomib, the model-predicted IC50 was comparable to that obtained from the in vitro concentration-effect experiments (3.5 versus 1.5 nM; Fig. 6).
Network-based approaches are increasingly used in drug discovery and development (Csermely et al., 2013; Harrold et al., 2013). Pharmacological networks are now commonly used for target identification (Sahin et al., 2009), evaluating signaling networks between normal and diseased conditions (Saez-Rodriguez et al., 2011), and understanding interactions among pathways (Thakar et al., 2007; Ge and Qian, 2009; Mai and Liu, 2009). For example, a Boolean network was used to identify novel targets to overcome trastuzumab resistance in breast cancer cell lines (Sahin et al., 2009), in which c-MYC was identified as a potential therapeutic target. In another example, Boolean network modeling was used to analyze immune responses due to the presence of virulence factors in lower respiratory tract infection (Thakar et al., 2007). Three distinct phases of Bordetellae infection were identified, which was not possible using traditional experimental methods (i.e., biochemical and molecular biology techniques). A similar approach is used in the present analysis, in which a Boolean network model is applied for evaluating bortezomib signal transduction pathways in U266 cells.
Our initial Boolean network model was developed based on the literature; however, available studies focus on either NFκB (Bharti et al., 2003b) or JAK/STAT3 (Bharti et al., 2003a; Park et al., 2011) pathways independently, and joint effects of these pathways have not been previously evaluated in U266 cells (e.g., apoptosis and cell proliferation). In addition, cellular protein dynamics and outcomes following exposure to specific probe compounds were also unavailable. Therefore, in vitro experiments and computational modeling were combined with pathway-specific inhibitors to evaluate the roles of both pathways individually and simultaneously. The Boolean network model confirmed that stress accumulation due to proteasome inhibition is a major pathway of myeloma cell death despite the activation of pNFκB protein expression, which contradicts a major proposed hypothesis of bortezomib effects.
Whereas logic-based models provide qualitative insight into network connectivity and drug-induced signal transduction, a quantitative dynamical model that includes important cellular biomarkers allows for the integration of critical factors responsible for cell death upon chemotherapy exposure (Yamazaki et al., 2011; Kay et al., 2012; Zhang et al., 2013). It is not yet practical to quantify all components in cellular signaling pathways; therefore, it is necessary to identify critical proteins for dynamical model development. A Boolean network reduction algorithm (Veliz-Cuba, 2011) identified critical factors regulating bortezomib cell death that guided development of a reduced PD model. Among the biomarkers (Supplemental Fig. 6), two antiapoptotic proteins and one proapoptotic protein (pNFκB, BclxL, and cleaved PARP) were integrated into a reduced pharmacodynamic model (Fig. 3). pNFκB was integrated in the model as bortezomib increases its expression, contradictory to proposed mechanisms. BclxL was selected as a major antiapoptotic protein that inhibits apoptosis, and cleaved PARP is a primary proapoptotic marker. The factor p21 was not incorporated in the current model as Boolean simulations suggest that bortezomib effects on apoptosis precede cell growth arrest.
The final cellular PD model (Fig. 3) was successfully qualified using the external data set following a low concentration of bortezomib (2 nM). Despite the fact that the cellular model was developed based on a single bortezomib concentration (20 nM) and linear coefficients, the final model reasonably predicted responses to the lower bortezomib exposure (Supplemental Fig. 8) and the dose-response curve for bortezomib at 48 hours (Fig. 6). Although there is a slight misfit for BclxL and cleaved PARP profiles (Supplemental Fig. 8, B and C), overall simulated profiles reasonably agree with experimental data. Most commonly, exposure-response relationships of protein expression data are sigmoidal in nature (Bharti et al., 2003a,b; Park et al., 2008, 2011). For example, the inhibition of pNFκB activity by curcumin is 10% at 1 μM and 40% at 10 μM (Bharti et al., 2003b). Similarly, cleaved PARP activation is about 1% at 10 μM curcumin but 100% at 50 μM curcumin (Park et al., 2008). It is therefore interesting that, although our model was developed using a single drug concentration and simple linear coefficients (e.g., Sm_M1B and Sm_PARP in Eqs. 1 and 4), it reasonably predicts profiles at a lower bortezomib concentration (Supplemental Fig. 8) and a sigmoidal concentration-effect relationship (Fig. 6). Furthermore, the predicted exposure-response relationship for BclxL can be overlaid with observed cell proliferation at 48 hours (Supplemental Fig. 9). This suggests that BclxL could serve as a potential biomarker to predict efficacy in MM and warrants further study. However, a potential limitation is that the proteins measured in this study reflect relative fold change and not absolute values of protein concentrations. Grounding on relative changes in protein expression maintains a degree of modularity, but more advanced proteomic methods are needed to quantitatively measure low-abundance signaling proteins. Cell growth and death parameters were also compared with available values from the literature. Natural cell growth and death were modeled using first-order growth rate constants (kg and kd). The estimated net growth rate constant (kg – kd) in the final model (0.00490 h−1) is comparable to net growth rate constants from U266 xenografts at 0.00487 h−1 (Siveen et al., 2014) and 0.00450 h−1 (Rhee et al., 2012) that were obtained by fitting an exponential model to digitized control curves (data not shown).
Preliminary experiments were conducted to measure pNFκB in the nucleus and cytoplasm of U266 cells after bortezomib exposure. Both cellular fractions showed a similar trend of increased relative pNFκB expression (data not shown); therefore, total cellular pNFκB was measured by Western blot analysis of cell lysates, revealing transient increased expression (Fig. 5A). Despite the stimulation of pNFκB protein expression, apoptosis is induced and the model is able to well capture this phenomenon. The NFκB pathway is one of the major cell survival pathways and is deregulated in many types of cancer (Karin, 2009). Bortezomib stimulates the upstream cascade of the NFκB pathway (Hideshima et al., 2009) leading to activation of pNFκB protein expression. Apoptosis is still induced through cellular stress resulting from the inhibition of proteasome (Fig. 4D) (Hideshima et al., 2001). The cellular model was able to well describe both mechanisms. Aberrant expression of pNFκB is responsible for the lack of efficacy and resistance in some cases of MM (Hideshima et al., 2002; Bharti et al., 2004). Greater expression of pNFκB results in greater expression of antiapoptotic proteins (e.g., BclxL) (Karin and Lin, 2002), with cross-talk between survival and apoptotic pathways determining cellular fate. Hence, greater drug concentrations are required to activate the apoptotic pathway. Initially, BclxL inhibition of cleaved PARP activation and cleaved PARP stimulation on removal of BclxL was incorporated. However, since BclxL exerts a negative feedback on cell death, as incorporated in the model, the BclxL effect on cleaved PARP was removed to simplify the model and reduce model redundancy.
In summary, systems pharmacology is an emerging field that seeks to couple systems biology and PD modeling, which could promote the discovery, development, and effective use of drugs based upon first principles. Computational tools based on graph theory (Liu et al., 2013), discrete dynamic relationships (Albert and Wang, 2009), and others can be used to identify critical system components within networks at multiple organizational levels, thus providing guidance for multiscale model construction and evaluation. We have used an integrated approach to investigate bortezomib signal transduction and exposure-response relationships in MM cells. A cell-based model of bortezomib in U266 human myeloma cells was developed that incorporates one of the major survival pathway proteins (pNFκB), antiapoptotic protein BclxL, and an apoptotic marker (cleaved PARP) to link bortezomib exposure to cell proliferation. The final model-predicted in vitro IC50 for bortezomib reasonably agrees with the experimental IC50 at 48 hours. Although the Boolean and dynamical models are relatively simple and specific to U266 cells and bortezomib, the overall strategic approach can be easily extended with slight modification to other myeloma cell lines (e.g., MM.1S and RPMI 8226) or bortezomib-based combination regimens (e.g., histone deacetylase inhibitors). This model may serve as a basis for studying bortezomib combinations with antimyeloma agents and to optimize xenograft combination studies, which is a focus of current research and will be reported separately. In addition, the model structure can be easily extended to describe responses in other cancer types that share similar pathways.
The authors thank Dr. John M. Harrold (University at Buffalo, SUNY) for developing MATLAB code for this project.
Participated in research design: Design and execution of experiments: Chudasama.
Conducted experiments: Chudasama.
Performed data analysis: Chudasama, Ovacik, Mager.
Wrote or contributed to the writing of the manuscript: Chudasama, Ovacik, Abernethy, Mager.
- Received March 28, 2015.
- Accepted July 9, 2015.
This work was supported by the National Institutes of Health National Institute of General Medical Sciences [Grant R01-GM57980] (to D.E.M.); the University at Buffalo-Pfizer strategic alliance (D.E.M.); an unrestricted training grant from Daiichi Sankyo Pharma Development (to D.E.M.); and an American Foundation for Pharmaceutical Education predoctoral fellowship (to V.L.C.).
- enhanced pharmacodynamics
- fetal bovine serum
- inhibitor of κB inhibitor
- inhibitor of κB kinase inhibitor
- janus kinase
- c-Jun N-terminal kinase
- multiple myeloma
- nuclear factor κB
- poly (ADP ribose) polymerase
- phospho–nuclear factor κB
- phospho–signal transducer and activator of transcription 3
- receptor-interacting protein
- signal transducer and activator of transcription
- U.S. Government work not protected by U.S. copyright