Pharmacodynamic Models of Differential Bortezomib Signaling Across Several Cell Lines of Multiple Myeloma.

The heterogeneous polyclonal nature of multiple myeloma complicates the identification of protein biomarkers predictive of drug response. In this study, a pharmacodynamic systems modeling approach was used to link in vitro bortezomib exposure and myeloma cell death. The exposure-response was integrated through a network of important protein biomarker dynamics activated by bortezomib in four myeloma cell lines. The pharmacodynamic models reasonably characterized the protein and myeloma cell dynamics simultaneously following bortezomib (20 nM) treatment. The models were used to identify differences in pathway dynamics across cell lines from model-estimated protein biomarker turnover parameters and global sensitivity analyses. Additionally, a statistical correlation analysis between drug sensitivity and model-fitted protein activation profiles (i.e., cumulative area under the protein expression-time curves) supported the identification of shared biomarkers associated with sensitivity differences among the cell lines. Both types of analysis identified similar important proteins associated with bortezomib pharmacodynamics, such as phosphorylated Nuclear Factor kappa-light-chain-enhancer of activated B cells (pNFkappaB), phosphorylated protein kinase B (pAKT), and caspase-8 (Cas 8).

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Multiple myeloma is a constantly changing disease due to clonal heterogeneity. The clonal populations of genetically diverse myeloma cells result in complex disease dynamics leading to heterogeneity in response to therapeutic interventions.

WHAT QUESTION DID THIS STUDY ADDRESS?

Heterogeneous responses to chemotherapy necessitate the need for adequate risk assessment and therapy management. This study includes a quantitative pharmacology approach to investigate the heterogeneous intracellular signaling mechanisms governing drug action and enables the identification of potential biomarkers associated with drug sensitivity.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

Sensitivity analyses of the final cell line specific dynamic models and an exploratory statistical analysis suggest potential protein biomarkers of bortezomib pharmacodynamic response.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

Upon adequate translation across a larger panel of cell lines, in vivo xenograft systems, and patient‐derived tumor samples, this approach could enable risk stratification based on biomarker profiles, aiding the timely identification of aggressive disease progression and enabling appropriate therapeutic interventions (i.e., precision medicine).

Multiple myeloma is a hematological cancer characterized by heterogeneity in its presentation and progression.1, 2 The changing genetic progression of the disease, similar to branching Darwinian evolution, introduces heterogeneity through the development of clonal populations of cells.3, 4 These subpopulations of cells proliferate with differing clonogenic potential and sensitivity to chemotherapy, which results in substantial interpatient variability in responses to the same treatment regimen as well as intrapatient variability across multiple cycles of treatment.3, 4 Multiple myeloma is incurable at present and is characterized by repeating cycles of remission and relapse, with an overall 5‐year survival rate of just 47%.5 This presents a significant challenge to treatment, which requires continual assessment of the clonal composition to develop appropriate treatment strategies.6 Thus, approaches are needed to translate tumor heterogeneity for appropriate risk assessment and the recommendation of safe and effective treatment plans. Integration of disease‐related and drug‐related (i.e., mechanism‐based) biomarkers in clinical trials is one such approach with considerable potential. Achieving this goal requires the identification of specific sets of markers that may influence cellular components and could be targeted by therapeutic agents, appropriate regimens, and/or combination therapies.

The objective of this study is to investigate the role of intracellular signaling protein expression within different myeloma cell lines that differ in their pharmacodynamic responses to bortezomib exposure.7 The differences in cell‐line responses were established via experimental and traditional pharmacokinetic/pharmacodynamic modeling, which showed that MM.1S and NCI‐H929 cells exhibited greater intracellular signaling, faster cell killing, and lower half‐maximal inhibitory concentration (IC50) values as compared to U266 and RPMI8226 myeloma cells. Bortezomib is a targeted proteasome inhibitor that influences many intracellular signaling proteins that regulate proliferation, cellular stress, and apoptosis pathways in myeloma cells.8 The drug triggers the nuclear factor‐kappa B (NFκB) and phosphoinositide 3 (PI3)/protein kinase B (AKT) pathways, stress‐associated Jun NH2‐terminal kinase (JNK) and p53 pathways, cell cycle inhibitory p21 and p27 pathways, and extrinsic (caspase‐8) and intrinsic (caspase‐9) apoptotic pathways.9 In this study, systems models integrating cell line specific signaling mechanisms were developed for the less sensitive U266 and RPMI8226 cell lines, which were subsequently tailored to describe bortezomib pharmacodynamics in the more sensitive MM.1S and NCI‐H929 cell lines. The four models are drug specific, individual cell‐type models that may facilitate the identification of common biomarkers. Global sensitivity analyses of the models enabled a comparison of proteins/pathways governing responses in each cell line, and results were in agreement with a comparative statistical correlation analysis. Such biomarker‐based pharmacodynamic models can provide insights into the determinants of system properties that are not possible with empirical models.

Materials and Methods

Pharmacodynamic systems models

The time‐courses of intracellular proteins following bortezomib exposure in four multiple myeloma cell lines were obtained from a published study.7 Briefly, all cell lines were treated with bortezomib (20 nM) and 10 intracellular proteins regulating bortezomib pharmacodynamics were measured at several time points. Bortezomib exposure was linked to myeloma cell death using protein biomarker dynamics, and the protein interactions in the models were guided by a qualitative logic‐based Boolean network describing signal transduction in myeloma cells.7 The models consist of a series of ordinary differential equations for multiple indirect‐response turnover profiles describing relative protein expression, with drug/protein(s) modulating synthesis or degradation rate constants under treatment conditions.10 Time‐dependent transduction models were used for temporal delays in network interactions as needed.11 The four model structures are shown in Figure 1 a–d and the detailed mathematical equations are provided in the Supplementary Materials.

Figure 1

Pharmacodynamic systems models for (a) U266, (b) RPMI8226, (c) MM.1S, and (d) NCI‐H929 cell lines. Bortezomib‐mediated intracellular protein dynamics link drug exposure with drug efficacy. Nodes highlighted in green represent proteins in prosurvival pathways, and nodes highlighted in red represent pro‐apoptotic proteins. Nodes with a dashed outline are transit compartments. Stimulatory and inhibitory connections between nodes are represented by ““ and “.”

The proteins were integrated into the pharmacodynamic end point of myeloma cell growth, and bortezomib‐induced cell death is defined as:

with cell growth and death described by first‐order growth (k ) and natural cell death (k ) rate constants. Cell growth was governed by proliferation pathway proteins pAKT and mammalian target of rapamycin (pmTOR). A cell cycle inhibitory protein (p21) was used to inhibit cell growth, whereas caspases stimulate cell death. Myeloma cell viability was reported as values normalized to untreated control cell viability at time 0 hour, thus the initial condition of Eq. (1) was set to 1. The temporal dynamics of proteins up to 24 hours were included in model development because cells showed complete cell death by that time. The degradation kinetics of bortezomib in media was described as mono‐exponential decay with a half‐life of 144 hours.12

Naïve pooled data analysis was used for model fitting, and parameter estimation was achieved using the maximum likelihood estimation algorithm in ADAPT5.13 Final models were selected based on goodness‐of‐fit criteria, which included model convergence, parameter estimation, and visual inspection of predicted vs. observed values and residual plots.

Global sensitivity analysis

Global Sobol sensitivity analysis was implemented in the open source Sensitivity Analysis Library in Python (SALib).14, 15 SALib was used to determine the sensitivity of the area under the effect curve (AUEC) of bortezomib‐induced cell death in 24 hours to changes in the individual protein turnover rate constant parameters, keeping all the other model parameters (stimulatory coefficients and power coefficients) fixed to their estimated values. The AUEC was calculated as the difference between the area under the cell proliferation and cell death curves. The rate constants were varied in the range of twofold for RPMI8226, MM.1S, and NCI‐H929, and ± 20% for U266 from their estimated values, depending on estimation accuracy of the Sobol indices, to generate a computationally feasible number of simulations.16 Samples of parameter sets were generated using the Saltelli sampling method that generates N × (2D + 2) samples, in which a value of 10,000 was used for N and D as the number of parameters being varied. Model output was simulated for all the generated parameter sets, and the large number of parameter sets ensured that adequate confidence intervals were obtained for the parameters with significant Sobol index values (above 0.05).17 The total‐order, first‐order, and second‐order indices along with their 95% confidence interval values were computed with SALib.

Statistical analysis

A Spearman's rank correlation coefficient was calculated between cell viability values under bortezomib treatment at 24 hours and the calculated cumulative relative protein expression AUEC (0–24 hours) of the model‐fitted profiles (calculated using Phoenix 64) for all proteins in SigmaPlot (version 12.0; Systat Software, Inc., San Jose, CA).

Results

Pharmacodynamic systems models of “less sensitive” cells

Bortezomib‐mediated responses were similar in U266 and RPMI8226 cells. In addition, the comparable sensitivity to bortezomib in terms of cell proliferation dynamics suggests similar mechanistic interactions among signaling proteins.7 The models for U266 and RPMI8226 (Figure 1 a,b) were used to describe protein profiles and cell proliferation data simultaneously, and the experimentally observed data and model fitted profiles are shown in Figure 2 (line plots depicting the model fits of the proteins are provided in Figures S1 and S2 in Data S1). The models successfully characterized the dynamics of the system, and all model parameters were estimated with good precision (Table 1), with the exception of the p21 degradation rate constant (k dp21) in U266 cells. The two models differed only in the transit compartments and power coefficient parameterization as described in the Methods section.

Figure 2

Model‐fitted time‐course of protein dynamics and cell proliferation in U266 (a,b) and RPMI8226 (c,d) cells. (a,c) The left panel depicts the experimentally observed and model fitted data. Log relative expression of the 10 proteins are on the Y‐axis, and time (in hours) is on the x‐axis. The dark red color denotes maximum activation and the dark blue denotes maximum inhibition. (b,d) Solid symbols represent experimentally observed cell proliferation for untreated‐control (black) and 20 nM bortezomib treatment (blue). Solid lines are model‐fitted curves for untreated‐control (black) and 20 nM bortezomib treatment (blue).

Table 1

Final parameter estimates for the four pharmacodynamic models

Parameter	Unit	Description	Value (% CV)
U266
k tpNFκB	Hour−1	pNFκB transit rate constant	0.03528 (3.17)
S mpNFκB	nM−1	pNFκB stimulatory coefficient	3.984 (7.28)
k tpAKT	Hour−1	pAKT transit rate constant	0.08517 (1.54)
S mpAKT	nM−1	pAKT stimulatory coefficient	0.5079 (0.989)
k dpmTOR	Hour−1	pmTOR degradation rate constant	0.01434 (3.16)
k dBcl‐2	Hour−1	Bcl‐2 degradation rate constant	0.01094 (6.19)
γBcl‐2	–	Bcl‐2 power coefficient	5.000*
k tpJNK	Hour−1	pJNK transit rate constant	0.07654 (0.127)
S mpJNK	nM−1	pJNK stimulatory coefficient	2.757 (0.256)
k dpp53	Hour−1	pp53 degradation rate constant	0.3222 (9.70)
k dp21	Hour−1	p21 degradation rate constant	4.709 (289)
γp21	–	p21 power coefficient	2.500*
k tpBAD	Hour−1	pBAD transit rate constant	0.08116 (2.99)
k dCas8	Hour−1	Caspase 8 degradation rate constant	0.3287 (3.17)
k dCas9	Hour−1	Caspase 9 degradation rate constant	0.9148 (11.0)
γCas9	–	Caspase 9 power coefficient	1.278 (0.689)
k g	Hour−1	Myeloma cell growth rate constant	0.01882 (0.0856)
k d	Hour−1	Myeloma death rate constant	0.003776 (3.19)
σ2Protein	–	Error coefficient for proteins	0.4724 (0.00210)
σ2MMcell	–	Error coefficient for myeloma cell	0.3645 (0.155)
RPMI8226
k tpNFκβ	Hour−1	pNFκβ transit rate constant	0.07685 (10.9)
S mpNFκβ	nM−1	pNFκβ stimulatory coefficient	0.04046 (13.3)
k tpAKT	Hour−1	pAKT transit rate constant	0.04528 (0.311)
S mpAKT	nM−1	pAKT stimulatory coefficient	1.011 (0.064)
k dpmTOR	Hour−1	pmTOR degradation rate constant	0.0007853 (0.676)
k dBcl‐2	Hour−1	Bcl‐2 degradation rate constant	0.008890 (10.9)
k dpJNK	Hour−1	pJNK degradation rate constant	0.07269 (2.70)
S mpJNK	nM−1	pJNK stimulatory coefficient	0.2965 (0.602)
k dpp53	Hour−1	pp53 degradation rate constant	0.04395 (3.14)
γpp53	–	pp53 power coefficient	1.746 (0.348)
k dp21	Hour−1	p21 degradation rate constant	0.01834 (0.911)
γp21	–	p21 power coefficient	5.000*
k tpBAD	Hour−1	pBAD transit rate constant	0.2257 (72.6)
γpBAD	–	pBAD power coefficient	0.3705 (2.41)
k dCas8	Hour−1	Caspase 8 degradation rate constant	0.1656 (34.8)
γCas8	–	Caspase 8 power coefficient	5.640 (16.2)
k dCas9	Hour−1	Caspase 9 degradation rate constant	0.008027 (9.52)
γCas9	–	Caspase 9 power coefficient	5.000*
k g	Hour−1	Myeloma cell growth rate constant	0.01680 (0.116)
k d	Hour−1	Myeloma death rate constant	0.05166 (11.0)
σ2Protein	–	Error coefficient for proteins	0.4324 (0.042)
σ2MMcell	–	Error coefficient for myeloma cell	0.4118 (0.187)
MM.1S
k tpNFκβ	Hour−1	pNFκβ transit rate constant	0.09507 (0.154)
S mpNFκβ	nM−1	pNFκβ stimulatory coefficient	5.054 (0.236)
k tpAKT	Hour−1	pAKT transit rate constant	0.1503 (1.32)
S mpAKT	nM−1	pAKT stimulatory coefficient	0.5823 (0.498)
k dpmTOR	Hour−1	pmTOR degradation rate constant	0.01159 (0.561)
k dBcl‐2	Hour−1	Bcl‐2 degradation rate constant	0.09135 (0.104)
γBcl‐2	–	Bcl‐2 power coefficient	5.000*
k tpJNK	Hour−1	pJNK transit rate constant	0.5659 (0.161)
S mpJNK	nM−1	pJNK stimulatory coefficient	0.3985 (0.0272)
k dp21	Hour−1	p21 degradation rate constant	0.1072 (2.16)
S mp21	nM−1	p21 stimulatory coefficient	0.09941 (0.412)
γp21	–	p21 power coefficient	1.200*
k dpBAD	Hour−1	pBAD degradation rate constant	0.02842 (0.741)
k dCas8	Hour−1	Caspase 8 degradation rate constant	0.4254 (0.676)
γCas8	–	Caspase 8 power coefficient	3.228 (0.0345)
k dCas9	Hour−1	Caspase 9 degradation rate constant	0.1439 (0.319)
k g	Hour−1	Myeloma cell growth rate constant	0.01901 (0.0395)
k d	Hour−1	Myeloma death rate constant	0.0003691 (2.95)
σ2Protein	–	Error coefficient for proteins	0.6237 (0.0162)
σ2MMcell	–	Error coefficient for myeloma cell	0.8254 (0.0721)
NCI‐H929
k tpNFκβ	Hour−1	pNFκβ transit rate constant	0.1123 (0.000363)
S mpNFκβ	nM−1	pNFκβ stimulatory coefficient	3.808 (0.00102)
k tpAKT	Hour−1	pAKT transit rate constant	0.1633 (0.275)
S mpAKT	nM−1	pAKT stimulatory coefficient	0.4531 (0.0167)
k dpmTOR	Hour−1	pmTOR degradation rate constant	0.005520 (0.0543)
k dBcl‐2	Hour−1	Bcl‐2 degradation rate constant	0.1105 (0.0888)
γBcl‐2	–	Bcl‐2 power coefficient	6.167 (0.0200)
k tpJNK	Hour−1	pJNK transit rate constant	0.1394 (0.0268)
S mpJNK	nM−1	pJNK stimulatory coefficient	0.2635 (0.0216)
k dp21	Hour−1	p21 degradation rate constant	0.07804 (0.274)
S mp21	nM−1	p21 stimulatory coefficient	3.627 (0.00830)
k dpBAD	Hour−1	pBAD transit rate constant	0.3298 (0.829)
γpBAD	–	pBAD power coefficient	1.103 (0.0217)
k dCas8	Hour−1	Caspase 8 degradation rate constant	0.2219 (0.0710)
γCas8	–	Caspase 8 power coefficient	3.163 (0.0116)
k dCas9	Hour−1	Caspase 9 degradation rate constant	1.544 (5.89)
k g	Hour−1	Myeloma cell growth rate constant	0.02278 (0.000380)
k d	Hour−1	Myeloma death rate constant	0.004006 (19.2)
σ2Protein	–	Error coefficient for proteins	0.7509 (0.000995)
σ2MMcell	–	Error coefficient for myeloma cell	0.2987 (0.000831)

% CV, percentage of coefficient of variation; JNK, Jun NH2‐terminal kinase; pAKT, proliferative protein kinase B; pmTOR, phosphorylated mammalian target of rapamycin; pNFκB, phosphorylated nuclear factor‐kappa B.

U266

Bortezomib‐mediated proteasome inhibition causes an accumulation of unwanted ubiquitinated proteins, resulting in severe stress in cells and leading to the activation of pJNK,18 pNFκB via the stimulation of upstream receptor‐interacting protein (RIP),12 and pAKT via unknown mechanisms.19, 20 Hence, the models included the direct stimulation of these proteins by bortezomib, characterized by linear stimulatory coefficients (Table 1). The relatively low value of the coefficient for pAKT suggests that stimulation of the AKT pathway is not a predominant bortezomib effect. Transit compartments (one for pJNK and pAKT and two for pNFκB) adequately described the delay in the stimulation of these proteins. A minor stimulation of the mTOR pathway was observed and was characterized by stimulation via the upstream regulator, pAKT.21

The pJNK stimulates the synthesis of the tumor suppressor p53,22 which in turn stimulates the synthesis of the cell cycle inhibitor p21. These proteins are primarily activated under stress/DNA damage conditions in order to halt cellular proliferation. Thus, proliferative pathways (pAKT) also play a role in the inhibition of these proteins.18, 22, 23, 24, 25, 26 The upregulation of pp53 and p21 was adequately described by the model; however, no expression of p21 was observed at 24 hours. One potential explanation is the relatively short activation half‐life of p21 (minutes to a few hours) and protein stability in acidic media resulting from dead cells at 24 hours.27, 28 The model was unable to capture the absence of p21 because the equation defining myeloma cell death is governed by the presence of p21 and caspases, which are needed to be active/present for maintaining continual cell death observed in in vitro experiments. This is the likely reason for the poor estimation of the p21 turnover parameter, k dp21.

Cas 8, the primary marker for the extrinsic apoptosis pathway, is activated by stress‐induced upregulation of pJNK and inhibited by cellular FLICE inhibitory protein (c‐FLIP), an apoptosis regulatory protein stimulated by pNFκB.8, 18, 29, 30 The profiles suggest a similar time to activation as pJNK. Hence, the synthesis of Cas 8 was driven by pJNK and pNFκB. Bortezomib also activates the intrinsic mitochondrial pathway of apoptosis, which involves the Bcl‐2 family of proteins and several caspases, such as Cas 9. Bcl‐2 is an anti‐apoptotic protein that resides on the mitochondrial membrane and inhibits release of cytochrome c, an activator of Cas 9 in the cytosol. BAD is a pro‐apoptotic protein present in the cytosol, which upon Ser128 phosphorylation by JNK, translocates to the mitochondria to form a pro‐apoptotic complex with Bcl‐xL (another anti‐apoptotic protein similar in function to Bcl‐2), leading to the release of cytochrome c and activation of apoptosis.31 In this current analysis, Ser112 phosphorylated BAD was measured, which is an inactive form of BAD that remains sequestered by the 14‐3‐3 family of proteins. JNK is also known to phosphorylate 14‐3‐3ς, releasing the sequestered Ser112 BAD to promote apoptosis.31, 32 Hence, pJNK was used to stimulate the degradation of Ser112 BAD in the model, and the delay in its inhibition was achieved by incorporating two transit compartments. Bcl‐2 dynamics is primarily governed by pBAD (BAD being the direct regulator), pNFκB (direct transcription), and pp53 (via PUMA),33, 34, 35, 36 and these profiles were characterized reasonably well by the model (Figure 2). Finally, the marker for intrinsic apoptosis (i.e., Cas 9) is regulated by Bcl‐2 (inhibiting activation of Cas 9), Cas 8, and pNFκβ (via XIAP) in the model. Cas 8 cleaves cytosolic BID, and its fragment (tBID) translocates to the mitochondria and oligomerizes with BAX (Bcl‐2‐associated X protein) or BAK (Bcl‐2 homologous antagonist/killer), causing release of cytochrome c from mitochondria and subsequent activation of Cas 9.35 The dynamics of Cas 9 was also adequately captured by the model (Figure 2).

RPMI8226

The relative stimulation of direct drug targets (pJNK and pNFκB) was lower than U266 cells, with smaller values of the estimated linear stimulatory coefficients (Table 1). Proteins that were similarly stimulated or inhibited in RPMI8226 and U266 cells included pAKT, pmTOR, pp53, p21, Bcl‐2, and pBAD. Unlike U266 cells, the fold‐change in expression of Cas 8 and 9 (twofold) was relatively less. Similar to p21 in U266 cells, the model was unable to characterize the loss of Cas 9 at 24 hours (due to short half‐life of active caspases37, 38) because continued presence of Cas 9 was essential to sustain in vitro cell death until 96 hours and prevent the prediction of regrowth at later time points.

Pharmacodynamic systems models of “more sensitive” cells

The U266 and RPMI8226 models served as base model structures for MM.1S and NCI‐H929 models. One of the key differences between the more sensitive and less sensitive myeloma cells was the lack of activation of pp53 upon bortezomib treatment as compared to untreated‐controls in MM.1S and NCI‐H929 cells. Although this could result from counteracting stimulatory and inhibitory signals, this required removal of the pp53 compartment from the models for these cells. However, the model then required inclusion of a pp53 independent stimulation of p21. Evidence suggests a direct activation of p21 by bortezomib via regulation of microRNA‐29b and HDAC in myeloma cells.39, 40 Hence, a direct stimulation of p21 by bortezomib was included in the models for MM.1S and NCI‐H929 cells with respective linear stimulatory coefficients. In addition, given the high relative expression change in select proteins in MM.1S cells (e.g., pNFκB), two protein‐protein interactions had to be removed to reconcile model fitting to the data: (i) pNFκB mediated stimulation of Bcl‐2, and (ii) pAKT mediated inhibition of p21. The absence of these two interactions may not affect the response outcome (i.e., cell proliferation).

The temporal profiles of protein expression between MM.1S and NCI‐H929 cells were similar in onset and extent of stimulation/inhibition, with the exception of a few proteins (e.g., Cas 8 and pNFκB), which exhibited significantly greater stimulation in MM.1S cells. The cell lines were also similar in cell proliferation responses, suggesting similar signaling interactions for regulating drug response. The models for MM.1S and NCI‐H929 cells (Figure 1 c,d) described the quantitative profiles of protein expression and myeloma cell proliferation simultaneously. A comparison of the experimentally observed data and model fitted profiles are shown in Figure 3 (line plots depicting the model fits of the proteins are provided in Figures S3 and S4 in Data S1). All system parameters associated with the protein and myeloma cell dynamics were estimated with good precision (Table 1).

Figure 3

Model‐fitted time‐course of protein dynamics and cell proliferation in MM.1S (a,b) and NCI‐H929 (c,d) cells. (a,c) The left panel depicts the experimentally observed and model fitted data. Log relative expression of the 10 proteins are on the y‐axis, and time (in hours) is on the x‐axis. The dark red color denotes maximum activation and the dark blue denotes maximum inhibition. (b,d) Solid symbols represent experimentally observed cell proliferation for untreated‐control (black) and 20 nM bortezomib treatment (blue). Solid lines are model‐fitted curves for untreated‐control (black) and 20 nM bortezomib treatment (blue).

MM.1S

The model for MM.1S cells included stimulation of pJNK,18 pNFκB,12 pAKT,19, 20 and p2139, 40 directly by bortezomib. A low value for the stimulatory coefficient for p21 (Table 1) suggests a relatively lower stimulation of the protein, which was adequately characterized, although no expression of the protein was observed after 6 hours. This could be due to the proposed explanation of a short half‐life and stability issues with the protein.27, 28 Delay in stimulation of pJNK, pNFκB, and pAKT required the incorporation of transit compartments (two each for pJNK and pNFκB, and one for pAKT). The mild stimulation of pmTOR mediated by pAKT was characterized well.8, 21

MM.1S cells exhibited a high magnitude of stimulation of Cas 8 (about 150‐fold), and the model was able to describe the profile well via regulation by pJNK and pNFκB.8, 18, 29, 30 The earlier onset and steeper inhibition of pBAD and Bcl‐2 in comparison to the less sensitive cells were characterized well. A ~10–20‐fold stimulation in expression of Cas 9 at 6 hours was also adequately described by the model. However, the model was unable to characterize the observed absence of Cas 9 at later time points, similar to RPMI8226 cells. The strong stimulation of Cas 8 in comparison to Cas 9 suggests that bortezomib relies more on the extrinsic apoptotic pathway activation in this cell line.

NCI‐H929

Bortezomib stimulation of three direct targets (i.e., pJNK, pNFκB, and pAKT) in NCI‐H929 cells was quantitatively similar to that of MM.1S cells, which explains the similarly estimated values of the linear stimulatory coefficients (Table 1). The coefficient for stimulation of p21 was greater in NCI‐H929 cells (S mp21 = 3.63 nM−1), as reflected in the greater relative expression of p21 in this cell line at 6 hours; however, the overall trend in p21 expression over time was unclear, possibly due to stability issues. Several proteins (pJNK, pAKT, pNFκB, pmTOR, Bcl‐2, and pBAD) were similarly stimulated or inhibited in the more sensitive cells, and the model described the dynamics of these proteins satisfactorily. Distinct from MM.1S cells, NCI‐H929 cells exhibited less stimulation of Cas 8 and 9. The model underestimated the profile of Cas 9 considerably, likely due to the lack of more time intensive data around the 6‐hour time point when the stimulation signal is most evident. Similar to MM.1S cells, the model for NCI‐H929 cells also fitted a sustained presence of Cas 9 and p21 at later time points.

Global Sobol sensitivity analysis

The results of the sensitivity analysis are shown in Figure 4. The sensitivity analysis for the U266 model was performed on turnover parameters for nine proteins, with the turnover rate of p21 fixed to its estimated value (due to uncertainty in its estimation). Interestingly, both U266 and RPMI8226 cells show the turnover rate constants of pJNK, pNFκB, and pAKT as the most significant parameters governing cell death. In addition, the rank‐order of the first five proteins was also similar in these cell lines. The total‐order sensitivity indices are essentially the sum of first‐order, second‐order, and further interaction orders of sensitivity indices.17 The total‐order and first‐order sensitivity index values were close for U266, suggesting absence of protein‐protein interactions regulating the output. However, the values differed for RPMI8226 indicating a significant pJNK and pAKT second‐order interaction term, although the confidence on this interaction index was poor (46%).

Figure 4

Sobol sensitivity indices of protein turnover parameters for all pharmacodynamic systems models and the correlation between drug sensitivity and biomarker activation. Total‐order (dark blue bar) and first‐order (light blue bar) sensitivity indices are shown for (a) U266, (b) RPMI8226, (c) MM.1S, and (d) NCI‐H929 models. The model is more sensitive to perturbations in parameters with greater sensitivity indices, and proteins are ordered on the basis of decreasing total‐order sensitivity indices.

Turnover rate constants of pJNK, pNFκB, and pAKT were identified as significant sensitivity indices for the MM.1S model, whereas the NCI‐H929 model identified pNFκB and pJNK turnover rate constants as significant. Although only the first few parameters have significant Sobol sensitivity indices in both models, the rank‐order of all parameters differed between the two cell lines. Additionally, for the NCI‐H929 model, the total‐order and first‐order index values were similar, except for one significant second‐order interaction between pJNK and pNFκB (SI = 0.0834; 31% confidence). Similarly, for the MM.1S model, the values differed indicating a significant pJNK and pNFκB second‐order interaction term, although the confidence on this interaction index was poor (31%). Another second‐order interaction that came close to significance was between pAKT and pNFκB (SI = 0.0457).

Correlation analysis between drug sensitivity and biomarker activation

This analysis was not powered for statistical significance owing to the small sample size of only four cell lines. The Spearman's rank correlation coefficients and a scatter plot showing the association between drug sensitivity and model‐fitted protein expression AUEC values for pNFκB, pAKT, and Cas 8 are shown in Figure 5 a,b. The results suggest inverse correlations (i.e., protein expression associated with decreasing cell viability) with proteins pAKT, pNFκB, and Cas 8. Cas 9, Bcl‐2, pBAD, and pmTOR exhibited moderate inverse correlations with cell viability. pJNK showed a slight positive correlation with viability (Figure 5), and cell viability was not correlated with model‐fitted p21 AUEC values (data not shown). A correlation with pp53 could not be assessed due to the sample size of only two cell lines. Although the correlations were not statistically significant (due to a small sample size of only four myeloma cell lines), they agree with the sensitivity analysis results and support the model‐based approach to identification of potential protein biomarkers.

Figure 5

Correlation analyses of protein expression and multiple myeloma cell viability. (a) Spearman rank correlation coefficients between cell viability under bortezomib treatment at 24 hours and model‐fitted area under the protein expression‐time curves (AUEC) until 24 hours are shown across the four myeloma cell lines. Negative and positive correlations are shown in red and blue bars. Correlation coefficients for pp53 and p21 are not shown. (b) Scatter plot showing the association between myeloma cell viability under bortezomib treatment at 24 hour and model‐fitted AUEC values until 24 hours for pNFκB, pAKT, and Cas 8. Cell lines are differentiated by color: U266 (red), RPMI8226 (blue), MM.1S (pink), and NCI‐H929 (dark green).

Discussion

The final pharmacodynamic models were constructed using intensive time‐series biomarker measurements of 10 important signaling proteins to explain myeloma cell growth and death under bortezomib exposure. Initially, temporal dynamics of proteins until 72 hours were included in model development. The observed data showed no expression of proteins p21, Cas 8, and Cas 9 at later time points (48–72 hours) along with sustained cell death until 96 hours in the treatment group. Cell growth inhibition and cell death were governed by p21 and Cas 8 and 9 in the model, and this led to misspecifications in model fits for p21, Cas 8, and Cas 9 (sustained activation of these proteins at later time points was noted) and the prediction of regrowth at later time points. Complete in vitro cell death was observed by 24–30 hours, and the observation of sustained cell death until 96 hours was actually an artifact of the experimental system because these studies permit evaluation of a limited number of cells. Viable intracellular protein signaling only happens until all the cells are dead (24–30 hours), and any protein signaling observed at later time points has uncertainty associated with its origin. Factors such as inherent stability of proteins in acidic media containing dead cells and activation half‐lives of proteins may play a role. Hence, protein signaling data were truncated at 24 hours for model development purposes. This change allowed for better characterization of protein profiles. The overall rich temporal profiles (data until 24 hours) allowed for the inclusion of a number of mechanisms pertaining to cell growth, stress, cytostasis, and apoptosis (intrinsic and extrinsic) pathways. In addition, the data provided confidence in the estimation of parameters simultaneously and precluded the need for fixing model parameters.

The final models (Figure 1) provide a detailed mechanistic framework and highlight critical biomarkers that may regulate bortezomib‐mediated drug response. The cell growth function (Eq. (1)) includes pAKT, pmTOR, and p21 on cell proliferation and Cas 8 and Cas 9 on cell death in a multiplicative manner. However, one limitation of the model is that the individual contributions of these factors cannot be identified using the current data. It was assumed that these signaling proteins act independently (and therefore multiplicative) in order to avoid additional unidentifiable model parameters. The quantitative determination of individual factor contributions will require further testing.

The expression (or lack thereof) of pp53 was a noteworthy distinction between the more and less sensitive cell lines. In contrast to the mutated TP53 isoform in U266 and RPMI8226 cells, MM.1S and NCI‐H929 cells contain the wild‐type TP53 gene.41, 42 Exposure to increasing concentrations of nutlin (inhibitor of the MDM2‐p53 interaction) in myeloma cells expressing Wild Type‐TP53 (MM.1S and NCI‐H929) caused an induction in the expression of p53, which was not observed in myeloma cells expressing mutated‐TP53 (U266 and RPMI8226 cells).42, 43 In addition, MM.1S and NCI‐H929 cells also overexpress MDM2 mRNA.41 These observations suggest an MDM2‐associated inhibition of p53 expression in MM.1S and NCI‐H929 cells.41, 42, 43, 44 Another mechanism for the lack of expression of p53 in these cells could be microRNA (miR‐25 and miR‐30d)‐based inhibition of p53 expression. NCI‐H929 cells exposed to an anti‐miR‐25 and anti‐miR‐30d show increased expression of p53.45 These mechanisms might explain the absence of pp53 stimulation in the more sensitive cells. Thus, the node for pp53 and its associated interactions were removed from the system models in these cell lines. Ideally, a unified model structure with cell line‐specific parameters would be preferable; however, discordance in p53 signaling precluded the development of such a model. It is unclear whether a unified model structure is possible, and more research is needed to address this issue.

The dynamics of p21 was unclear in more sensitive cells (Figure 3). The cell lines presented an increase in expression at 6 hours, albeit with large variability, followed by no expression at later time points, possibly due to the stability concerns and short half‐life of the protein.27, 28 The increased expression of p21 was evident in MM.1S cell lines treated with 2 nM bortezomib,7 which strongly suggested an activation of p21 in these cells following bortezomib exposure. Accordingly, the observed p53 independent stimulation of p21 was specified in the model as direct bortezomib stimulation, a mechanism supported by a bortezomib‐mediated increase in expression of miR‐29b downregulating HDAC and upregulating p21.39, 40 The difference in the system models for MM.1S and NCI‐H929 cells were specific to only minor secondary protein interactions, and the overall model structures were similar. Collectively, the models established that bortezomib mechanisms are fundamentally conserved, from a systems perspective, and are amenable to incorporation of functional differences due to the heterogeneity in signaling protein expression.

Global sensitivity analysis highlighted the importance of upstream proteins (e.g., pJNK, pNFκB, and pAKT; Figure 4) directly regulated by bortezomib, emphasizing the role of signal relay from upstream to downstream proteins eventually resulting in cellular apoptosis. The relay process introduces several opportunities for the presentation of signaling heterogeneity. One of the possible reasons for observing the same proteins as sensitive across models could be the nature of the data. Relative expressions from untreated controls do not render estimation of real synthesis and degradation rate constants.46 Absolute quantities of intracellular signaling proteins are difficult to measure due to their low abundance and challenges in isolating these proteins.

The exploratory Spearman rank correlation analysis between model‐fitted protein activation AUCs and drug sensitivity across all cell lines (Figure 5) supplemented the sensitivity analysis and highlighted similar proteins that are associated with differences in bortezomib sensitivity. Some of the key proteins include: pAKT, pNFκB, Cas 8, Bcl‐2, and pBAD. An increase in the expression of pAKT and pNFκB upon drug treatment is probably due to inherent cellular mechanisms that enhance prosurvival pathways. Decreasing cell viability is also a result of increasing expression of Cas 8 and increasing inhibition of Bcl‐2 and pBAD. Both approaches suggest pAKT, pNFκB, and Cas 8 as key proteins associated with bortezomib pharmacodynamics despite heterogeneous expression among the cell lines. The proteins identified by these approaches could serve as potential biomarkers, and their expression may be predictive of responses and provide an opportunity for directing therapy toward achieving a desired outcome.

The pharmacodynamic models were developed for a single bortezomib treatment concentration and cannot account for potential dose‐dependent nonlinearities and delays in signal transduction. The models also assume linear stimulation of proteins by bortezomib and other proteins. The models differed in the number of transit compartments for pJNK, pNFκB, pAKT, and pBAD, which could suggest unknown and/or dissimilar upstream processes. These processes could affect responses at the in vivo translational stage, in the presence of combination therapy, and may affect the generalizability of the model across cell lines. However, these limitations do not preclude the hypothesis that biomarker‐based modeling approaches can be used to explain pharmacological heterogeneity. The developed system models are mechanistic and provide a framework that enables comparison across cell lines. Bortezomib has been shown to exhibit a synergistic interaction when given sequentially with vorinostat (an HDAC inhibitor).47 The final bortezomib models (Figure 1) could be coupled with a mechanistic model of vorinostat pharmacodynamics48 to further explore the biomarkers associated with this combination and may provide a platform for assessing other bortezomib‐based chemotherapy combinations (e.g., Bcl‐2 inhibitors49). In addition, the models can be potentially extended to other myeloma cell lines and to include a range of doses for improving their translational potential.

Funding

This work was supported, in part, by an unrestricted educational grant from Daiichi Sankyo.

Conflict of Interest

As an Associate Editor for CPT: Pharmacometrics & Systems Pharmacology, Donald E. Mager was not involved in the review or decision process for this article.

Author Contributions

V.R. and D.E.M. wrote the manuscript. V.R. and D.E.M. designed the research. V.R. performed the research. V.R. and D.E.M. analyzed the data.

Data S1. Pharmacodynamic Systems Models and Figures S1‐S4.

Click here for additional data file.

Data S2. Model code.

Click here for additional data file.