Skip to main content
Download PDF

Predictive model identifies strategies to enhance TSP1-mediated apoptosis signaling

Cell Communication and Signaling volume 15, Article number: 53 (2017) Cite this article

Abstract

Background

Thrombospondin-1 (TSP1) is a matricellular protein that functions to inhibit angiogenesis. An important pathway that contributes to this inhibitory effect is triggered by TSP1 binding to the CD36 receptor, inducing endothelial cell apoptosis. However, therapies that mimic this function have not demonstrated clear clinical efficacy. This study explores strategies to enhance TSP1-induced apoptosis in endothelial cells. In particular, we focus on establishing a computational model to describe the signaling pathway, and using this model to investigate the effects of several approaches to perturb the TSP1-CD36 signaling network.

Methods

We constructed a molecularly-detailed mathematical model of TSP1-mediated intracellular signaling via the CD36 receptor based on literature evidence. We employed systems biology tools to train and validate the model and further expanded the model by accounting for the heterogeneity within the cell population. The initial concentrations of signaling species or kinetic rates were altered to simulate the effects of perturbations to the signaling network.

Results

Model simulations predict the population-based response to strategies to enhance TSP1-mediated apoptosis, such as downregulating the apoptosis inhibitor XIAP and inhibiting phosphatase activity. The model also postulates a new mechanism of low dosage doxorubicin treatment in combination with TSP1 stimulation. Using computational analysis, we predict which cells will undergo apoptosis, based on the initial intracellular concentrations of particular signaling species.

Conclusions

This new mathematical model recapitulates the intracellular dynamics of the TSP1-induced apoptosis signaling pathway. Overall, the modeling framework predicts molecular strategies that increase TSP1-mediated apoptosis, which is useful in many disease settings.

Background

Angiogenesis, the formation of new capillaries from pre-existing blood vessels, plays a critical role in tumor progression. Angiogenesis enables the tumor to generate its own blood supply and obtain oxygen and nutrients from the microenvironment. This process is regulated by a dynamic interplay between the angiogenic promoters, such as vascular endothelial growth factor (VEGF) and fibroblast growth factor (FGF), as well as angiogenic inhibitors, such as thrombospondin-1 (TSP1) [1,2,3,4,5].

Due to its importance in tumor development, invasion, and metastasis, angiogenesis has become a prominent target for cancer therapies. In addition to strategies targeting pro-angiogenic species, such as inhibiting VEGF signaling using antibodies and tyrosine kinase inhibitors, anti-angiogenic species hold promise in reducing tumor angiogenesis. TSP1 is a well-known, potent endogenous angiogenesis inhibitor. TSP1 expression is lost in multiple cancer types; however, its re-expression can delay cancer progression, promote tumor cell apoptosis, and decrease microvascular density. For these reasons, it has been of interest to mimic TSP1’s functions in regulating angiogenesis [3, 6,7,8,9].

TSP1 is a multifunctional matricellular protein that acts to inhibit angiogenesis in multiple ways [2, 10, 11], which include altering the availability of pro-angiogenic factors and promoting anti-angiogenic signaling through its receptors CD36 and CD47. Several studies have shown that TSP1 mediates its anti-proliferative and pro-apoptotic effects in a highly specific manner on endothelial cells. TSP1 primarily promotes these effects by binding to the CD36 receptor [3, 12, 13], which is associated with capillary blood vessel regression [10, 12, 14, 15]. TSP1 interaction with CD36 leads to recruitment of the Src-related kinase Fyn, activation of p38MAPK, and processing of caspase-3, a vital protease that mediates apoptosis [12, 15, 16]. TSP1-CD36 signaling also causes transcriptional activation of Fas ligand (FasL), a death ligand that also promotes pro-apoptotic signaling, ultimately inhibiting angiogenesis. This apoptosis pathway is further enhanced as pro-angiogenic factors induce increased levels of Fas receptors, sensitizing the cells to FasL [17].

Unfortunately, therapies that mimic TSP1 activity have not demonstrated definitive clinical efficacy. For example, ABT-510, a TSP1 peptide mimetic that binds to CD36, was previously tested in a Phase II study in 2007 for treatment of metastatic melanoma. However, the drug failed to reach its primary endpoint (18-week treatment failure rate), resulting in termination of the study [18]. ABT-510 also showed little clinical effect in a Phase II trial for advanced renal cell carcinoma [19]. These disappointing results indicate that there is a need to better understand the effects of anti-angiogenic agents and develop effective treatment strategies. This requires a detailed and quantitative understanding of the dynamic concentrations of the factors involved in angiogenesis signaling.

Computational systems biology offers powerful tools for studying complex biological processes that involve a large number of molecular species and signaling reactions that occur on multiple time- and spatial-scales. Systems biology aims to study how individual components of biological systems give rise to the function and behavior of the system [20]. Additionally, computational modeling aids in the development of therapeutic strategies that specifically target tumor angiogenesis to optimally inhibit tumor progression, complementing pre-clinical and clinical angiogenesis research [21].

Substantial research has focused on the pro-angiogenic factors and their extracellular interactions [21,22,23]. However, a consideration of the intracellular mechanisms of anti-angiogenic factors is also needed in order to fully understand the dynamics of the signaling networks regulated by angiogenesis promoters and inhibitors. In this study, we focus on TSP1-mediated apoptosis signaling through the CD36 receptor. Although some aspects of the TSP1-CD36 pathway have been studied experimentally, the signaling network has not been quantitatively and systematically analyzed. We constructed the first computational model that describes the intracellular signaling network induced by TSP1-CD36 binding in endothelial cells, a complex network comprised of biochemical reactions that lead to cell apoptosis. We applied the model to predict the effects of modulating protein expression and enzyme activity on apoptosis signaling. The model quantifies the effects of these perturbations and predicts promising targets, both in terms of the averaged response of a population of endothelial cells and individual cells within the population. Thus, the model is a quantitative framework to predict strategies to enhance TSP1-mediated apoptosis. Ultimately, the model can be used to identify novel pharmacologic targets and optimize therapeutic strategies that promote apoptosis and, subsequently, inhibit angiogenesis.

Methods

Mathematical model

We constructed a computational model of TSP1-mediated apoptosis signaling via the CD36 receptor in endothelial cells. The molecular interactions depicted in Fig. 1 were translated into biochemical reaction equations, with the assumption that the reactions follow well-established kinetic laws, including mass-action or Michaelis-Menten kinetics (Additional file 1: Table S1). A system of nonlinear ordinary differential equations (ODEs) was formulated to describe the rate of change of the species’ concentration. The model is comprised of 53 ODEs to predict the concentrations of the 53 species in the signaling network over time. The SimBiology toolbox (MATLAB) was used to implement the biochemical reaction equations, and the MATLAB stiff solver ODE15s was used the numerically solve the system of ODEs. The model file is provided in the Additional file 2.

Fig. 1
figure 1

Model schematic of TSP1-mediated apoptosis signaling. TSP1 binding to the CD36 receptor recruits p59fyn, which induces activation of the caspase-3 cascade. The kinase p38MAPK is subsequently phosphorylated and translocated to the nucleus. NF-κB translocates into the nucleus and is activated in the presence of phosphorylated p38MAPK. This leads to transcriptional activation of FasL. FasL protein binds to its receptor Fas, forming the DISC complex, which binds to c-FLIP (FL) and procaspase-8 (pro8) to form the p43-FLIP complex. This complex activates IKK, which releases NF-κB from its inhibitor IκB. Blue arrows indicate transport reactions

Full size image

Solving the set of ODEs with the baseline initial conditions provides the averaged response of a population of cells. Additionally, we accounted for heterogeneity in a population of cells by solving the ODE model 2000 times, each with a different set of initial conditions. We refer to this as the “population-based model”.

Cytosolic and nuclear compartments

The model is comprised of two compartments, cytosolic and nuclear, both assumed to be well mixed. Specific molecules, such as NF-κB, may move from one compartment to another at a defined translocation rate. The volume of the nuclear compartment is estimated to be 14.32% of the cytosolic compartment [24], and the concentrations of species transported between the two compartments are converted using this ratio.

Initial protein concentrations

The initial conditions used in the model are given in Additional file 1: Table S2. Very few references for initial concentrations of proteins are available. Therefore, we adapted values from a previously published model [25] and adjusted the initial concentrations of several species in order for the model to match experimental measurements. For the CD36 and Fas receptors, we used flow cytometry to quantify the average numbers of receptors on cultured human microvascular endothelial cells (data described below; similar to previous work [26], and converted the receptor numbers to concentrations using the total cell volume of 1 picoliter [27].

When simulating the population-based model, we randomly selected the initial conditions from a gamma distribution. The gamma distribution is characterized by two parameters: the shape factor, a, and the scale parameter, b. These parameters are related to the mean, m, and standard deviation, sd, of the distribution: 1/a = sd 2/m 2; b = sd 2/m. Thus, a × b = m. We set m to be the baseline value of the initial condition for each species (given in Additional file 1: Table S2) and assumed a shape factor of 5.5 (based on previous work [28]).

Rate constants

All baseline model parameter values are listed in Additional file 1: Table S3.

Production of soluble species

The basal rate at which each species is synthesized (K syn_all ) is set to be 10−4 μM/min, with the exception of FasL, procaspase-8, and procaspase-3, whose production rates are described below.

The model accounts for FasL production mediated by TSP1, and we described the production of FasL mRNA production (DNA transcription) using Michaelis-Menten kinetics:

$$ \mathrm{V}={\mathrm{V}}_{\max \_\mathrm{FasL}}\ast \mathrm{NF}\hbox{-} \kappa \mathrm{B}\_\mathrm{p}/\left({\mathrm{K}}_{\mathrm{m}\_\mathrm{FasL}}+\mathrm{NF}\kappa \mathrm{B}\_\mathrm{p}\right) $$

where V max_FasL and K m_FasL are the Michaelis-Menten kinetic rate constants for FasL mRNA production, and NF-κB_p is the activated transcriptional factor that catalyzes this process. The molecular details involved in FasL protein production encompass the mRNA translocation and translation, and protein secretion. The rates involved in these reactions are not readily available in published literature. Therefore, we estimated the values in model fitting in order to match experimental data.

The synthesis rate of procaspase-8 and procaspase-3 were assumed to be dependent on the concentration of DISC present in the system, as a partial effect of Fas ligation. The synthesis rate is described as:

$$ \mathrm{V}=\mathrm{F}\ast \mathrm{DISC}+{\mathrm{K}}_{\mathrm{syn}\_\mathrm{all}} $$

where F is a hand-tuned coefficient, DISC is the complex formed by FasL binding to Fas, and K syn_all is the basal level synthesis rate assigned to all the other species except for FasL.

Protein degradation

Protein species are assumed to be degraded at the same rate, 10−3 min−1, unless there was a degradation rate available in the literature or from a previous model. This allows the system to balance and reach steady-state in the absence of TSP1 stimulation. The degradation rates of caspase-8, caspase-3, the p43:FLIP:IKK_a complex, and cytosolic NF-κB have unique values adapted from previous modeling work by Neumann et al. [25].

Receptor-ligand interactions

The affinity of TSP1 and its receptor CD36 has been measured experimentally: the K d value is 230 nM [29]. We assumed that FasL binds to Fas with an affinity of 0.4 nM. In all cases, the dissociation rate for the receptors is 1.2 × 10−2 min−1. Receptors are internalized and inserted at the cell membrane such that the total number of receptors (ligated plus unbound) is constant.

FasL cascade

The model includes DISC formation upon FasL binding with Fas, and the downstream caspase-8 and NF-κB activation reactions. The molecular details were adapted from the model established by Neumann et al. [25]. We altered this portion of their model by adding reversible binding reactions to ensure the reaction network is consistent with the other parts of our model. We tuned the universal dissociation rate K off to be 1.2 × 10−2 min−1 to match the data presented in their paper. The simulations of the implemented minimal model are shown in Additional file 3: Figure S1.

Sensitivity analysis

There is limited quantitative experimental data available to specify the values of the kinetic parameters. However, the parameters must be set to appropriate values in order for the model to generate reliable predictions. We first used sensitivity analysis to reduce the number of parameters to be estimated. Specifically, to identify the influential kinetic parameters before each step of model fitting, we conducted global sensitivity analysis using the extended Fourier Amplitude Sensitivity Test (eFAST) method [30], as we have done in previous work [21, 22]. All inputs were allowed to vary simultaneously one order of magnitude above and below the baseline value, and the effects of multiple inputs on the model outputs were quantified. An additional global sensitivity analysis was performed after model training, in order to quantify the robustness of the model with respect to varying the kinetic parameters (Additional file 3: Figure S2). Sensitivity analysis was also used to determine the effects of initial protein concentrations to inform perturbation simulations (Additional file 3: Figure S3).

Quantification of experimental data

The experimental data used in model fitting and validation are extracted from previously published studies [12, 15]. Jimenez et al. stimulated human microvascular endothelial cells (HMVECs) with 10 nM TSP1 and used immunoblotting of cell lysates to measure TSP1-induced association of activated p59fyn with CD36 over 30 min (TSP1:CD36:p59fyn complex and activated p59Fyn, pp59), and activated p38 (pp38) over 60 min. We analyzed the immunoblots using ImageJ (https://imagej.nih.gov) and extracted quantitative data needed for model fitting and validation. The local background from bands was subtracted and their intensity was quantified. The intensities of subject species were normalized to the corresponding control band intensities. One data point (TSP1:CD36:p59fyn concentration at 30 min) was excluded due to low image quality.

Two sets of data for caspase-3 activity upon TSP1 stimulation measured with fluorescence assays were extracted: from Jimenez et al. and Nor et al.. HMVECs stimulated with 5 nM TSP1 over 5 h (300 min), and human dermal vascular endothelial cells (HDMECs) stimulated with 0.388 nM TSP1 over 12 h (720 min) were measured in these studies, respectively. We quantified the caspase-3 activity at each timepoint directly from the published figures using ImageJ.

Receptor quantification

We measured CD36 and Fas receptor numbers on HMVECs, following methods previously established [31]. Briefly, HMVECs from Lonza were cultured in flasks and maintained in Endothelial Cell Growth Media-2 (EGM-2) supplemented by the EGM-2 Single Quot Kit (Lonza). Cells were maintained at 37 °C in 95% air and 5% CO2, and we only use cells at passage numbers 2–4. To dissociate cells from the culture plate, cells were incubated with a non-enzymatic cell dissociation solution, CellStripper (Corning), for 5 min at 37 °C. Cells were centrifuged at 500×g for 5 min to obtain a final concentration of 4 × 106 cells/mL in stain buffer (PBS, bovine serum albumin, and sodium azide).

Aliquots of cells (25 μL, ~105 cells) were labeled with phycoerythrin (PE)-conjugated monoclonal antibodies (Biolegend) and allowed to incubate on ice for 45 min. The volumes of antibody solution used (15 μL for CD36 receptor; 10 μL for Fas receptor) were the optimal volumes as determined by saturation experiments. Cells were then washed with ice-cold stain buffer, centrifuged twice at 500×g and re-suspended in 200 μL stain buffer prior to single-cell analysis via flow cytometry to quantify the number of receptors per cell.

Flow cytometry was performed on a MACSQuant flow cytometer (Miltenyi), and FlowJo (BD Biosciences) software was used to analyze the data. To identify dead cells, 5 μg/mL Sytox Blue (Thermofisher) was added to all samples, and tubes were vortexed immediately prior to placement in the flow cytometer. Cells exhibiting very low Sytox Blue fluorescence were identified as live cells, and gating was performed to collect 10,000 live cells for each sample. Finally, the gated cells were examined in a plot of forward scatter area (FSC-A) versus side scatter area (SSC-A) to identify the single-cell population.

To determine the number of receptors per cell, the fluorescence of Quantibrite PE beads (BD Biosciences) was measured. We measured the fluorescence for beads with different numbers of binding sites (as specified by the manufacturer). We applied linear regression to the fluorescence measurements, constructing a calibration curve to convert the geometric mean of PE fluorescence to the number of bound molecules. The average number of receptors on a cell in the population was then calculated using the linear regression and the cell fluorescence data. For each experiment, two biological replicates were used, and the experiments were repeated 3–4 times. We report the mean and standard error of the mean of the measurements of all samples: CD36 = 24,372 ± 2365 receptors/cell and Fas = 7860 ± 395 receptors/cell. The receptor distributions of representative samples of each receptor are shown in Additional file 3: Figure S4.

Parameter estimation

The estimation of the kinetic parameters was achieved using the “lsqnonlin” function in MATLAB, as done in our previous work [22, 32, 33]. This algorithm solves the nonlinear least squares problem using the trust-region-reflective optimization algorithm, minimizing the weighted sum of the squared residuals (WSSR). The minimization is subject to the upper and lower bounds of the free parameters. One hundred runs were performed in each fitting step, and a global sensitivity analysis was performed with the best fit parameter values (the parameters that produce the lowest WSSR). The step-wise iteration was repeated four times to ensure fine-tuning of the parameter values. Parameter values used in the implemented model are from the best fit (lowest WSSR) from the last step. We also report the mean and standard deviation of the estimated parameter values in Additional file 1: Table S3 and Additional file 3: Figure S5.

Definition of apoptotic cells

Cleaved poly(ADP-ribose) polymerase (cPARP) is the output of the model used as an indicator of apoptosis, since loss of intact PARP results in failure to repair DNA damage. Our model simulations show that the dynamics of cPARP follows a switch-like action; however, the range of cPARP varies widely depending on the initial concentration of PARP. Previous study [34] has shown that low doxorubicin (DXR) dosage with 10 nM TSP1 stimulation resulted in approximately 50% of the cells becoming apoptotic in 24 h. Therefore, we simulated this treatment condition using the population-based model, and determined the cPARP concentration that results in 50% cell apoptosis. We then use this concentration, 1.05 μM, as the defined threshold that needs to be reached for cell apoptosis to occur. Thus, the definition of which cells are apoptotic is based on literature evidence.

Simulated perturbations to TSP1-mediated apoptosis

We applied the model to simulate seven specific perturbations to the intracellular signaling network, to find strategies that enhance apoptosis signaling. Below, we list the motivation and literature evidence for each of the seven perturbations. We also describe how the perturbation was simulated in our mathematical model. The abbreviations listed in parentheses are also used in the results figures. Generally, perturbations are simulated in the ODE model by adjusting the baseline initial conditions (Additional file 1: Table S2) and parameter values (Additional file 1: Table S3).

  1. 1.

    XIAP downregulation (“XIAP”): Experimental studies show that downregulation of X-linked inhibitor of apoptosis protein (XIAP) can promote apoptotic signaling [35,36,37,38]. We simulated this effect by reducing XIAP concentration to 0.5-fold of the baseline value.

  2. 2.

    Low dosage doxorubicin treatment (“DXR”): Experimental studies [34, 39] have shown that a low dose of doxorubicin upregulates the expression of Fas receptor and other protein species. We simulated this effect by increasing the initial Fas receptor level by 3-fold and K syn_all by 10-fold.

  3. 3.

    Phosphatase inhibition (“Ptase”): Studies have shown that inhibiting MAPK phosphatase (MKP) activity can promote apoptosis signaling [40, 41]. We simulated this effect by decreasing the association rate (K on_dephos ) of the phosphatase with phosphorylated p38MAPK (pp38) and the dephosphorylation rate (K dephos ) by 10-fold.

  4. 4.

    Kinase promotor (“Kp”): Literature evidence suggests that the tumor microenvironment likely upregulates many kinases’ activities in the tumor-related endothelial cells [42,43,44]. We simulated the kinase promoter by increasing the phosphorylation rates of p59fyn, p38MAPK, and IκB by 10-fold.

  5. 5.

    Procaspase-3 upregulation (“pro3”): Global sensitivity analysis (Additional file 3: Figure S3) and baseline model simulations (Fig. 3) indicate that upregulation of procaspase-3 increases the cPARP level upon TSP1 stimulation. Therefore, we investigated the effect of procaspase-3 upregulation by increasing procaspase-3 concentration by 3-fold.

  6. 6.

    Fas upregulation (“Fas”): Experimental investigation by Quesada et al. suggests that upregulation of the receptor Fas promotes TSP1-induced apoptosis [34]. We simulated the effect of Fas upregulation by increasing Fas concentration by 3-fold.

  7. 7.

    Translocation rate increase (“Ktrsp”): Based on the structure of the signaling network, we hypothesized that increasing the cytoplasm-to-nucleus transport would enhance apoptosis. We simulated the effect of faster cytosol-to-nucleus translocation by increasing the translocation rate (K trsp ) by 10-fold.

Analysis of sensitivity and specificity of a cPARP-based classifier

Here, we consider binary classification of the cells’ response to TSP1 stimulation based on cPARP levels at 24 h: apoptotic or non-apoptotic. We aim to classify the cells as apoptotic or non-apoptotic using certain model variables as predictors (i.e., the initial species’ concentrations). The goal is to determine whether the initial amounts of one or more species are accurate predictors of what the response to TSP1 stimulation would be. That is, whether the cell will become apoptotic or not. We constructed the ROC curve to determine which model variables are accurate predictors. Here, the “actual response” is the classification of a cell as apoptotic or non-apoptotic based on its cPARP level predicted by the mechanistic model of TSP1-mediated apoptosis signaling presented above.

For a binary classification system such as this, there are four possible predicted outcomes for a given cPARP cutoff value: true positive, a cell predicted to be apoptotic is actually apoptotic; false positive, a cell predicted to be apoptotic is actually non-apoptotic; true negative, the predicted and actual response are both non-apoptotic; false negative, a cell predicted to be non-apoptotic is actually apoptotic. This analysis determines the fraction of positives predicted correctly (sensitivity or the true positive rate) and the fraction of true negatives predicted (specificity or the true negative rate) for different cutoff values of cPARP.

To evaluate tradeoffs between sensitivity and specificity, we constructed a receiver operator characteristic (ROC) curve for cPARP. The ROC curve plots the true positive rate versus the false positive rate (1-specificity). An ideal input maximizes true positives, with minimal false positives (i.e., the (0,1) point on the ROC graph). An ROC curve that lies on the 45-degree angle line indicates that the input does not classify the output any better than a random guess, where the area under the ROC curve (AUC) is 0.5. Thus, having an AUC value significantly greater than 0.5 indicates that the input can be used to classify the data. We performed the ROC analysis using the custom “roc” function in MATLAB.

Results

Model training and validation

We constructed a model of the signaling network of TSP1-mediated apoptosis in endothelial cells based on literature evidence. TSP1 binds to CD36, activating caspase-3, the core executioner protease. Caspase-3 promotes apoptosis by cleaving PARP in endothelial cells. Activation of caspase-3 also mediates intracellular signaling leading to the production of FasL, a death ligand that binds to its receptor Fas on endothelial cells and further promotes apoptosis through activation of caspase-3 [12, 15, 17]. The signaling network, illustrated in Fig. 1, includes several important feedback loops involved in TSP1-mediated apoptosis, including the caspase cascade (caspase-3 activates its activator, caspase-8) and Fas signaling (TSP-1 promotes the production of Fas, which also activates caspase-3). We implemented the signaling network mathematically to generate an ODE model, assuming that the reactions follow mass-action or Michaelis-Menten kinetics rate laws.

The model was trained using quantitative experimental data and validated with an independent set of measurements. We extracted experimental data from the literature in order to calibrate the model and estimate the kinetic parameters. Specifically, the fold-changes in the caspase-3 activity and the levels of three intracellular species (TSP1:CD36:p59fyn, pp59fyn, and p38MAPK) upon TSP1 stimulation were quantified from Western blot data and used to train the ODE model.

We used a step-wise strategy comprised of global sensitivity analysis and parameter estimation to ensure that the model could match the training data (see Methods). As a result of this approach, we obtained 12 sets of parameters that enable the model to closely reproduce the training data (Fig. 2a-e, Additional file 1: Table S3).

Fig. 2
figure 2

Model training and validation. The ODE model was trained to match experimental measurements of activated species in the TSP1-mediated apoptosis signaling pathway. (a) TSP1:CD36:p59fyn [12]; (b) pp59fyn [12]; (c) pp38MAPK [12]; (d) caspase-3 activity [12]; and (e) caspase-3 activity [15]. (f) An independent set of data for caspase-3 activity under the condition of p38MAPK inhibition [12] was used to validate the model prediction. Solid line: mean of 12 best fits. Shaded area: 95% confidence interval. Squares: experimental data

Full size image

After fitting the model to the experimental data, we used a separate set of measurements to validate the model predictions. Here, we applied the trained model to predict the dynamics of caspase-3 activity when p38MAPK is inhibited, mimicking an experimental study from Jimenez et al. [12]. This inhibitory effect on p38MAPK is simulated by setting the phosphorylation rate of NF-κB by active p38MAPK (pp38MAPK) to be zero. The 95% confidence interval of the model predictions produced with all 12 sets of parameters shown in Fig. 2f is not visible. This indicates that the parameter sets produce similar dynamics of caspase-3 activity with p38MAPK inhibition. The model qualitatively matches this independent set of data (Fig. 2f, Additional file 1: Table S3), where caspase-3 activity is reduced at 300 min, compared to the case without p38 MAPK inhibition (Fig. 2d). Overall, the model fitting and validation produces a trained model that generates reliable predictions related to the dynamics of TSP1 simulation. Results from a representative set of parameter values are shown in Additional file 3: Figure S6, where the baseline model is simulated to produce the dynamics of all 53 species upon 24-h simulation with 10 nM TSP1. Notably, TSP1 decays rapidly, and cPARP has a sigmoidal shape.

A global sensitivity analysis was performed to reveal the robustness of the trained model. The sensitivity of all 53 species in the model with respect to changes in the parameter values (Additional file 3: Figure S2) and species with non-zero initial concentrations (Additional file 3: Figure S3) was computed. These results show that the model output, cPARP, is largely influenced by the concentrations of its immediate effectors (procaspase-3, XIAP, and PARP), as well as critical parameters identified and estimated during model training. Upstream or intermediate species, such as those involved in p38 signaling and FasL signaling (Additional file 3: Figure S3, bottom panel), are sensitive to changes in a variety of initial concentrations and parameters values.

Altering the concentrations of intracellular signaling species influences the apoptotic response

We first applied the trained and validated model to investigate the effects of varying the concentrations of cell surface receptors and intracellular signaling species, in combination with different TSP1 stimulation levels. In this study, we specifically focus on predicting the concentration of cleaved PARP (cPARP) as an indicator of cell apoptosis. Caspase-3 promotes apoptosis by cleaving PARP, and cleavage of PARP by caspases is considered a hallmark of apoptosis [45]. Sensitivity analysis revealed that the concentrations of procaspase-3, XIAP, and PARP most significantly influence the cPARP level throughout the simulated time course (Additional file 3: Figure S3). This analysis suggests that varying the concentrations of those intracellular species can impact TSP1-mediated apoptosis signaling. Since the receptor concentration influences the initial dynamics of TSP1 stimulation, we also hypothesized that increasing the receptor’s availability (i.e., increasing the receptor:ligand ratio) can amplify the signaling induced by ligand-receptor binding. Therefore, we ran the model and individually altered the expression level (initial conditions) of the CD36 or Fas receptors, or intracellular species procaspase-3, XIAP, and PARP, by 10-fold above and below the baseline values. We applied this relatively large alteration in the protein expression levels to explore the extent of changes in the model output. The initial conditions were varied for each of the 12 fitted parameter sets, and we compared the cPARP level at various time points for each case.

Across the simulated time points, there is a dose-dependent response to TSP1, where increasing the concentration of TSP1 increases the predicted cPARP concentration. Interestingly, altering the expression levels of the CD36 or Fas receptors does not affect the cPARP level, compared to the baseline model (Fig. 3a and b). This result holds true for all TSP1 concentrations investigated, and is in accordance with the findings from the global sensitivity analysis, which identified CD36 and Fas as non-influential to the cPARP level (Additional file 3: Figure S3).

Fig. 3
figure 3

Dose-dependent response of apoptosis signaling with varied initial concentrations. Initial concentrations of (a) CD36, (b) Fas, (c) procaspase-3, (d) XIAP, and (e) PARP were varied 10-fold above (right column) and below (left column) the baseline values (center column). The model was used to simulate cPARP level in response to four different TSP1 concentrations: 0.1, 1, 10, and 100 nM. The predicted cPARP level at 24 h was generated using the 12 best sets of parameter values for each condition. The mean cPARP concentration is plotted; error bars show the standard deviation

Full size image

In contrast, varying the initial concentrations of the influential species led to significant changes in cPARP levels. Increasing the amount of procaspase-3, the unprocessed form of caspase-3, by 10-fold leads to increased cPARP level at every simulated time point, as compared to the baseline model (Fig. 3c, right panel). Downregulation of procaspase-3 to 0.1-fold of the baseline value slightly decreased cPARP level at intermediate time points (6 to 12 h), but did not affect the cPARP level at 24 h, compared to the baseline model.

Regulation of the caspase-3 inhibitor XIAP reduces apoptosis signaling. That is, increasing the level of XIAP by 10-fold dramatically decreased cPARP concentration to less than 20% of the baseline level, as shown in the right panel of Fig. 3d. However, decreasing XIAP by 0.1-fold results in a larger and faster increase in cPARP level compared to the baseline model (Fig. 3d, left panel). For example, after 24 h, the decreased XIAP resulted in 41% and 34% more cPARP than the baseline level, with 0.1 nM and 100 nM TSP1, respectively.

Lastly, the model predicts that increasing PARP levels significantly influences cPARP levels (Fig. 3e). When PARP is increased by 10-fold, the cPARP level at all time points is approximately nine times higher than the amount produced in the baseline model. In summary, the apoptotic response stimulated by TSP1 is sensitive to varying the concentrations of certain intracellular species.

Perturbing the signaling pathway influences the population response to TSP1 stimulation

Next, we implemented perturbations in the model and predicted the response of individual cells in a population. We accounted for heterogeneity in the cell population by varying the initial concentrations of protein species. Cellular heterogeneity is observed for multiple dimensions of single cell measurements, and detailed molecular differences can be used to distinguish cell-to-cell variation [46]. In this population-based model, the initial concentrations of all starting species are drawn from a gamma distribution [28, 47] (see Methods for details). Here, we focus on extrinsic noise (i.e., variability in the protein levels), as opposed to intrinsic variations (fluctuations in the rates of the biochemical reactions), since several studies have demonstrated that the experimentally-observed cell-to-cell heterogeneity is largely due to extrinsic rather than intrinsic noise [28, 48,49,50,51]. We performed the simulations with the baseline model and the parameter set that best fit the data out of the 12 parameter sets obtained from model training. We then ran the model 2000 times, representing 2000 independent cells, and analyzed the population-level response to TSP1 stimulation. We simulate the response to seven conditions (as described in the Methods) at two TSP1 concentrations (0.1 and 10 nM). Particularly, we investigated whether the predicted results from solving the deterministic model with the fixed initial conditions (Fig. 3) hold true when accounting for heterogeneity at the population level.

We characterized the population-level response based on the cells’ cPARP concentration. We extracted predicted intracellular cPARP concentrations for the 2000 cells at distinct time points up to 24 h of TSP1 stimulation, and generated histograms. This provides a direct visualization of the distribution of cPARP levels in the cell population. Based on literature data, we defined the threshold of intracellular cPARP required for apoptosis to occur within each cell to be 1.05 μM (see Methods). We used the model to predict the percentage of apoptotic cells, based on the predicted cPARP concentrations. Cells that have high cPARP level (greater than 1.05 μM) are classified as apoptotic, since their cPARP level exceeds the threshold value. Cells whose intracellular cPARP concentration is below the threshold value are classified as non-apoptotic. We also analyzed when cells that have high cPARP level appear at the simulated time points.

In the baseline model, the apoptotic response initiates within 6 h after 10 nM TSP1 stimulation (Fig. 4a). The size of the cPARP-positive population increases throughout the 24-h stimulation. By 24 h, the apoptotic cells make up 41% of the total population (Fig. 4e). Below, we compare the population-level response for the baseline model to the response when particular species in the intracellular signaling network are perturbed.

Fig. 4
figure 4

Distribution of cPARP concentration in population-level model. a-d: Histogram showing the percentage of the 2000 cells with a given cPARP concentration, in response to 10 nM TSP1 stimulation. (a) Baseline model; (b) XIAP downregulation; (c) DXR treatment; and (d) Increased nuclear translocation rate. A different color is assigned to each time point. The cPARP threshold is marked by a solid line and the region in the x-y plane beyond the threshold is shaded as light purple. (e-h): The predicted percentage of non-apoptotic (black) and apoptotic (purple) cells in response to 10 nM TSP1 stimulation. (e) Baseline model; (f) XIAP downregulation; (g) DXR treatment; and (h) Increased nuclear translocation rate

Full size image

According to the model with fixed initial conditions, downregulation of XIAP strongly promotes apoptotic signaling (Fig. 3d, left panel). To explore whether this conclusion still holds with a heterogeneous cell population, we decreased XIAP expression level by 0.5-fold, a physiologically reasonable change to the protein expression, and simulate the population-based response to 10 nM TSP1 stimulation under this condition. The results show that with XIAP downregulation, cells with high cPARP level appear at 4 h (Fig. 4b, f). By 24 h, 57% of the cell population is apoptotic. Additionally, the apoptotic population exceeds the non-apoptotic population by 0.4-fold (Fig. 4f).

We also considered the effect of doxorubicin (DXR) on the apoptosis signaling network. A published experimental study suggest that a low dosage of DXR sensitizes cells to pro-apoptotic signaling [34]. Specifically, Quesada and coworkers observed that Fas receptor expression increases approximately 3-fold following DXR treatment. Thus, we simulated the effect of such a DXR treatment by increasing Fas expression by 3-fold. Additionally, we increased the synthesis rates of certain intracellular species, as DXR has been shown to increase protein expression [39]. Here, we increased K syn_all by 10-fold. Under this simulated DXR treatment condition, cells with high cPARP level appear within 2 h, a faster apoptotic response than in the baseline model or with XIAP downregulation (Fig. 4c). However, the population of cells with high cPARP is 52%, which is not as large as what is predicted with XIAP downregulation (57%). Additionally, with DXR treatment, the progressive increase in the percentage of apoptotic cells through 24 h is more gradual than with XIAP downregulation (Fig. 4f, g). By 24 h, 52% of the cells are apoptotic, and there are 0.2-fold more apoptotic cells than non-apoptotic cells.

Upon examining the structure of the signaling network, we hypothesized that increasing the translocation rate of phosphorylated p38MAPK (pp38) and NF-κB into the nucleus can promote apoptotic signaling. Therefore, we simulated the model with K trsp increased by 10-fold. The apoptotic response is slower and smaller in scale than in the baseline model (Fig. 4d). Additionally, the positive population does not appear until 10 h after starting TSP1 stimulation, and by 24 h, less than 30% of the cells are apoptotic (Fig. 4h).

We also simulated the population response with procaspase-3 upregulation, increased Fas expression, phosphatase inhibition, and kinase-activity upregulation. The results are shown in Additional file 3: Figure S7. These perturbations to the signaling network do not dramatically affect the population response. That is, the speed and magnitude of the response in each case are similar to the baseline model (Fig. 4a and e). Apoptotic cells begin to appear by 6 h, and after 24 h of TSP1 stimulation, at least 40% of the cells are apoptotic (Additional file 3: Figure S7). Additionally, we predicted the population-based response for the baseline model and the seven network perturbations when the cells are stimulated with 0.1 nM TSP1. These results are shown in Additional file 3: Figure S8. Next, we present our detailed analysis of the predicted results for the two TSP1 stimulation levels and compare the apoptotic response.

Strategies to enhance apoptotic response have differential effects on the magnitude and time scale of TSP1-mediated signaling

We applied the model to distinguish the effects of possible strategies to promote TSP1-induced apoptosis under different levels of TSP1 stimulation. Here, we compared three quantities: the percentage of cells that have reached the cPARP threshold by 24 h, the maximum cPARP level reached in each cell, and the time to reach the threshold cPARP concentration (T t ) within the 24-h simulation time. Combined with the results shown in Fig. 5, these quantities provide insight into TSP1-mediated apoptosis signaling. Below, we compare the results to the baseline case.

Fig. 5
figure 5

Predicted population-based response to TSP1 stimulation. Comparison of three quantities that characterize the population-level response: (a) the percentage of cells that have reached the cPARP threshold by 24 h, (b) the maximum cPARP level reached in each cell, and (c) the time to reach threshold (T t ). Asterisks in panels (b) and (c) indicates the p-value from ANOVA: ****, p < 0.0001; ***, p < 0.001; **, p < 0.01; *, p < 0.1; ns, not significant

Full size image

In the baseline model, 0.1 nM TSP1 stimulation leads to 37% of the cells being apoptotic at 24 h. The model predicts that 49% more cells are apoptotic with XIAP downregulation than the baseline (Fig. 5a, left panel). Additionally, phosphatase inhibition and DXR treatment increased the apoptotic population by 14% and 9%, respectively, compared to the baseline model. Upregulation of procaspase-3 or kinase activity both increase the apoptotic population by approximately 7%. Increasing Fas expression did not have an effect on the percentage of apoptotic cells, while increasing the translocation rate decreased the apoptotic cells by 27%.

With 0.1 nM TSP1 stimulation, the median values for the maximum cPARP attained under each simulated condition follow the same order of effectiveness as observed with the percentage of apoptotic cells (Fig. 5b, left panel). Additionally, we performed statistical analyses to compare the maximum cPARP between the baseline model and each perturbation (Fig. 5b, left panel, asterisks above each column). The maximum cPARP reached is highly significantly different from the baseline model when XIAP is downregulated (the maximum cPARP is higher; p < 0.0001) or when the translocation rate is increased (cPARP decreases; p < 0.0001). Notably, the effects of altering XIAP or the nuclear translocation rate are more significantly different from the baseline than all of the other strategies.

Interestingly, the effectiveness of the approaches on time to reach threshold does not follow the same order as that of the percentage of apoptotic cells or of maximum cPARP (Fig. 5c, left panel). Based on our statistical analysis, compared to the baseline level, the time to reach threshold is significantly shorter (with high significance level) with XIAP downregulation, phosphatase inhibition, DXR treatment or procaspase-3 over-expression (Fig. 5c, left panel, asterisks above each column). In contrast, the time to reach threshold is significantly longer when the translocation rate is increased by 10-fold.

With 10 nM TSP1 stimulation, XIAP downregulation leads to 40% more apoptotic cells than the baseline model (Fig. 5a, right panel). Phosphatase inhibition and DXR treatment also lead to a strong response, with 17% and 27% more apoptotic cells, respectively. Increasing procaspase-3 or kinase activity both increased the relative size of the apoptotic population by approximately 10%, compared to the baseline.

XIAP downregulation, phosphatase inhibition, DXR treatment, and increase of procaspase-3 expression significantly increased the maximum cPARP level and shortened the time to reach threshold (Fig. 5b and c, right panels). Increasing the translocation rate significantly decreased the maximum cPARP reached, and prolonged the time to reach threshold.

In summary, these results show that XIAP downregulation is more effective than the other approaches in both increasing the maximum cPARP level attained in the cell population and shortening the time to reach the cPARP threshold. This holds true when cells are stimulated with either 0.1 nM or 10 nM TSP1. Increasing the level of TSP1 stimulation to 10 nM makes all of the pro-apoptotic strategies more effective. However, it is interesting to note that the ordering of the strategies from most effective to least effective changes with the level of TSP1 stimulation.

Initial protein expression levels influence apoptosis response

In order to explore the cause for the different responses to the apoptosis signaling within the population of 2000 cells, we compared the initial conditions of the apoptotic cells and non-apoptotic cells at 24 h. Statistical analysis of the distributions of the initial protein concentrations (normalized to their baseline values) indicate that XIAP concentration in the apoptotic population is significantly lower than in the non-apoptotic population (p < 0.0001), while PARP concentration is significantly higher in the apoptotic population (p < 0.0001) (Additional file 1: Table S5). In fact, the relationship between the apoptotic response and the XIAP and PARP initial conditions is easily visualized (Fig. 6). This illustrates that apoptotic cells (Fig. 6a, purple markers) have higher PARP and lower XIAP than non-apoptotic cells (Fig. 6a, black markers). The NF-κB concentration in the cytosolic compartment is also significantly lower in the apoptotic population compared to apoptotic cells (p = 0.04). Thus, statistical analysis shows that the initial concentrations of certain species distinguish apoptotic from non-apoptotic cells. Moreover, the distribution of the initial concentrations of XIAP and PARP are significantly different for apoptotic versus non-apoptotic cells (Fig. 6b and c; Additional file 1: Table S5).

Fig. 6
figure 6

Relationship between initial conditions and predicted apoptotic response. (a) Initial concentrations of PARP and XIAP for apoptotic cells (purple) and non-apoptotic cells (grey). The difference between the initial amounts of PARP and XIAP is larger for apoptotic cells, which reach higher cPARP levels. (b) Histogram showing distribution of initial XIAP concentration for apoptotic cells (purple) and non-apoptotic cells (grey). The apoptotic cells have relatively lower initial XIAP concentrations. (c) Histogram showing distribution of initial PARP concentrations for apoptotic cells (purple) and non-apoptotic cells (grey). The apoptotic cells have higher initial PARP concentrations

Full size image

Building on the results from the statistical analysis of the initial concentrations, we investigated whether the initial concentrations can be used to classify cells as apoptotic or non-apoptotic upon 24 h of TSP1 stimulation. Here, we generated a receiver operator characteristic (ROC) curve. The ROC curve illustrates the ability of an input descriptor to classify a response with high sensitivity and specificity (see Methods). For inputs, we used the initial concentration of each species that does not start at zero (16 species), the ratios of XIAP and PARP concentrations ([XIAP]/[PARP] and [PARP]/[XIAP]), and the absolute value of the difference between the concentrations of PARP and XIAP (i.e., |[PARP]-[XIAP]|). Thus, we considered 19 total inputs that may possibly predict the response to TSP1-mediated apoptosis. We used the area under the ROC curve (AUC) to compare the ability of the inputs to predict the response to TSP1.

Constructing the ROC curve for the 19 inputs shows that the absolute difference between the initial concentrations of XIAP and PARP predicts which cells will be apoptotic. Specifically, the difference between XIAP and PARP can classify the cells with high sensitivity and specificity (Fig. 7). In this case, the AUC is 0.97, with a 95% confidence interval of [0.96–0.98], providing a quantitative measure of the predictive ability of the difference between the XIAP and PARP concentrations. This AUC is significantly different than 0.5 (p < 0.0001), which indicates that using the difference between the concentrations of XIAP and PARP is more predictive than classifying cells as apoptotic or not based purely on chance. The AUC values for classifying the apoptotic response using the initial concentration of XIAP or PARP alone are 0.65 and 0.95, respectively, and in both cases, the AUC values are significantly different from 0.5 (p < 0.0001). Thus, although using the concentration of XIAP or PARP alone is better than randomly guessing which cells will respond, these concentrations are less reliable predictors when considered individually. The AUC values for the remaining 16 inputs range from 0.48 to 0.54, and are not significantly different than randomly selecting which cells will be apoptotic. Quantitative results from constructing the ROC curve for all of the inputs are provided in Additional file 1: Table S6. Overall, the results of this analysis demonstrate that the initial concentrations of PARP and XIAP, and especially the difference between their concentrations, can be used to predict which cells will respond to TSP1 signaling.

Fig. 7
figure 7

ROC curves for classifying the apoptotic response. Blue: ROC curve for the difference between the initial concentration of PARP and XIAP input; AUC is 0.97. Purple: ROC curve for the initial PARP concentration input; AUC is 0.95. Cyan: ROC curve for the initial XIAP concentration input; AUC is 0.65. All of three AUC values are significantly different than 0.5 (given by the gray dashed line, p < 0.0001), indicating that these inputs are reliable predictors of the cells’ response to TSP1 stimulation

Full size image

Discussion

We have developed a molecular-detailed model of the TSP1-induced apoptotic signaling. The model represents the reaction network of interactions involved in the intracellular signaling pathway and includes multiple feedback loops. Our model captures the feedback from transcriptional regulation by NF-κB onto Fas signaling, which allows us to expand the dynamics of the network to a longer time scale. Predictions from the trained model match experimental data. We further validated the model using a separate set of experimental measurements. We implemented the model using ODEs, which can be solved once to simulate the average response of a population of cells. We also accounted for randomness in the protein concentrations by solving the ODEs 2000 times with varied initial conditions to predict the individual responses of 2000 cells. In this case, the cells each have different initial concentrations of the molecular species, representing heterogeneity of the cell population. This heterogeneity influences the responses to TSP1 treatment and the effectiveness of strategies aiming to increase apoptosis signaling.

The model provides mechanistic and quantitative explanations of the effects of several approaches to promote TSP1-mediated apoptosis. Using the model, we proposed and simulated the effects of perturbing the signaling network, including altering receptor and protein levels, rates of protein synthesis and transport, the activity of specific phosphatases, and the overall kinase activity within the cells. These simulations exploit the power of mathematical modeling to generate quantitative predictions that would otherwise be time- and cost- consuming to explore experimentally. Overall, our model provides quantitative insight into the effects of targeting particular aspects of the TSP1-mediated apoptosis signaling pathways.

Global sensitivity analyses reveal insight into the robustness of the structure of the signaling network. We find that the model output cPARP is most significantly influenced by the concentrations of its immediate effectors, procaspase-3, XIAP, and PARP. Cleaved PARP is insensitive to changes in the initial conditions of upstream signaling species in the network; however, its effectors are influenced by the upstream signaling molecules. One such example is procaspase-8. The species cFLIP-L (FL) forms a complex with procaspase-8 involving DISC, and promotes activation of caspase-8 or NFκB. Caspase-8 then cleaves caspase-3, which in turn cleaves PARP. Therefore, cPARP is indirectly influenced by upstream signaling. However, the multiple feedback loops regulating this network appear to attenuate or buffer the effects of the upstream signal on the end output cPARP.

The model predicts that downregulation of XIAP is a promising strategy to enhance TSP1’s apoptotic signaling effect, sensitizing cells to low-dosage TSP1 treatment. XIAP is a potent inhibitor of caspase activity [35]. Experimental studies have shown that specific over-expression of XIAP can rescue cells from apoptosis, and antisense downregulation of XIAP led to a dramatic decrease in resistance to radiation-induced cell death [36,37,38, 52, 53]. Our model simulation shows that in TSP1-mediated apoptotic signaling, modulating XIAP level also has a similar effect. Importantly, the model provides quantitative and mechanistic insight into the effects of targeting XIAP. Downregulation of XIAP directly leads to an increased level of active caspase-3, a crucial mediator in the apoptosis signaling network. Our simulation results demonstrated that this approach is the most effective in promoting endothelial cell apoptosis. The further analysis based on the ROC curve supports the model simulations showing the importance of XIAP in influencing the dynamics of TSP1-mediated apoptosis. Interestingly, the ROC curve reveals a relationship between XIAP and PARP that we had not identified using model simulations alone. The difference between the initial concentrations of those proteins can accurately predict which cells will undergo apoptosis in response to TSP1 stimulation, even before the cells are exposed to TSP1. This highlights one potential clinical application of our work, to predict the response to pro-apoptosis signaling. With our model, it may be possible to determine whether a TSP1-based anti-angiogenic treatment that targets tumor endothelial cells would be effective. For example, endothelial cells isolated from a tissue sample obtained from a cancer patient’s tumor biopsy can be analyzed to determine the intracellular XIAP and PARP levels. Those measurements can be input into the model and used to inform whether the treatment would be effective. Although validating the model for this purpose is beyond the scope of this study, our work provides a foundation to pursue such investigation.

In another example, a strategy to increase apoptosis signaling is supported by experimental studies. We used the model to quantify the effect of inhibiting MAPK phosphatase (MKP) activity. We simulated the effect of this approach by decreasing the binding affinity between the phosphatase and its substrate, pp38, and the dephosphorylation rate. Therefore, the active pp38 level remains high as the phosphatase activity is inhibited, enabling downstream signaling. The results indicate that this approach can effectively promote the apoptotic response. Our predictions agree with experimental results that show that modulating MKPs is a viable option to promote apoptosis mediated by p38MAPK [40, 41].

The model also generates non-intuitive results. The simulations show that increasing the rate of translocation from the cytosol to the nucleus delayed and attenuated the apoptotic response. This observation is counterintuitive, as faster translocation of species immediately upstream of FasL production is expected to accelerate signal transduction. However, the model simulation suggests that with faster translocation, the pool of caspase-3 and p38MAPK is depleted in upstream signal transduction, before DISC formation occurs (data not shown). The effect of a pan-kinase activity promoter is another example of non-intuitive predictions. We simulated the kinase promoter by increasing the phosphorylation rates of p59fyn, p38MAPK, and IκB by 10-fold. Intriguingly, this approach did not affect either the response time (time to reach apoptosis threshold) or the percentage of apoptotic cells. The explanation is that increased kinase activity depleted certain species before their accumulation is achieved, in a similar manner to the faster cytosol-to-nucleus translocation case.

In addition to proposing strategies to increase apoptosis signaling, the model provides mechanistic insight into experimental observations. At low dosages, doxorubicin treatment has been shown to promote TSP1-mediated apoptosis. Using our model, we propose the potential mechanism of action of this effect. A study by Quesada et al. showed that endothelial cell apoptosis was due to a synergistic effect of the upregulation of FasL induced by TSP1 and upregulation of Fas by doxorubicin [34]. Our model simulations show that altering Fas receptor expression level alone does not affect the apoptotic response; rather, the availability of intracellular signaling species significantly influences apoptosis signaling as well. Specifically, the model simulations show that increasing the synthesis of particular molecular species is required to qualitatively match the effects of DXR observed experimentally. Based on these modeling results, we hypothesize that the low dosage doxorubicin treatment not only regulates the protein expression of Fas, but also other species. Our predictions agree with other experimental studies that show translation of multiple proteins increases when the cells are exposed to stress [25, 30]. Our model simulation demonstrates the effect of a low dosage of doxorubicin through the proposed mechanism of action, which can be further tested with experiments.

This work provides a biophysically realistic network that generates reliable predictions of the population-level responses. We incorporated variability into the ODE model to represent a heterogeneous population of cells. This implementation provides a framework to test how molecular-targeted strategies influence individual cells in a population. Solving the ODE model only once for a single set of initial conditions indicates that upregulation of procaspase-3 greatly increases the magnitude of the apoptotic response, while procaspase-3 downregulation does not affect cPARP level. This implies that a threshold value of procaspase-3 expression may be needed in order to enhance apoptosis signaling. On the other hand, the population-based model simulations showed that overexpressing procaspase-3 by 3-fold has only a mild effect of increasing the apoptotic response. These conflicting results demonstrate the importance of accounting for heterogeneity within a cell population. Therefore, we emphasize the utility of the population-based model to make predictions, as the deterministic model that represents the dynamics of an average single cell can be misleading in certain cases.

To our knowledge, this is the first mechanistic model to investigate TSP1-mediated apoptosis. The apoptotic signaling in this study is essential to the inhibitory effect of TSP1 [15]. Notably, TSP1 not only induces apoptosis through the CD36 receptor, but also has many other anti-angiogenic functions, possibly making it more potent than a single apoptosis inducer for regulation of angiogenesis. The framework established in our study can be readily adapted and combined with models describing other interactions of TSP1 [54], or pro-angiogenic signaling networks [22, 23] for future modeling studies. The model can also be expanded to include cell-cell interactions with both exogenous and endogenous FasL signaling. Future work can improve the model framework in various ways, for example, by adding the downstream function of PARP to address the balance between survival and apoptotic effects as PARP loses its repair function upon cleavage; specifying initial concentrations from different cell types; including more detailed reaction mechanisms (such as NF-κB activation by p38MAPK); or accounting for the mitochondrial feed-forward loop, which provides the link to the intrinsic apoptotic pathway. Finally, we note that strategies to increase apoptosis of diseased cells may also impact normal endothelial cells. We can expand the model to study both normal and diseased cells, which may have differential receptor expression levels, as observed for VEGF receptors in normal versus tumor cells [26, 55]. In this way, we can predict strategies that more specifically target the diseased population of cells.

Our model complements other studies that focus on apoptosis signaling promoted by death ligands and their receptors [24, 25, 56]. We adapted the model for Fas-mediated apoptotic and NF-κB signaling established by Neumann et al., which served as a foundation for the FasL signaling cascade in our work. In other work, Albeck et al. established a model to describe TNF-related apoptosis-inducing ligand (TRAIL) induced apoptosis in HeLa cells (an ovarian cancer cell line), with a focus on the “variable-delay, snap-action” switching mechanism of extrinsic apoptosis. Both the Albeck model and our model exhibit the cells’ response of switching to a high apoptotic response. However, one key difference between our model and theirs is how the switching arises. In the Albeck model, all-or-none switching of the cell death response is achieved by a network that does not include feedback. They concluded that this snap-action arises from interplay between the biochemistry of protein-protein interaction and translocation between the cytosol and mitochondria. We have simplified this extrinsic apoptosis pathway in our model; however, the snap-action behavior of apoptosis is still evident, shown by the fast accumulation of cPARP. This switching is a direct result of the reaction between caspase-3 and PARP, amplified by the feedback loops. Another difference between the two models is how stochasticity is implemented. The time to reach threshold (T t ) in our model is analogous to the delay time (T d ) in Albeck model. Albeck and co-workers impose a distribution in the delay time by randomly selecting T d from a defined range. In contrast, in our model, cell-to-cell differences in the switching delay emerge solely based on the variation in intracellular protein concentrations. This is a more realistic framework that represents an actual population of cells. Thus, our model is a tool to analyze the population-level responses to TSP1 stimulation and perturbations to the signaling network.

Conclusions

In summary, our model quantitatively describes the TSP1-mediated intracellular signaling via the CD36 receptor, which leads to endothelial cell apoptosis. This model predicts that downregulation of XIAP is the most promising way to effectively promote TSP1-mediated apoptosis. In addition, we propose an alternative mechanism of action for the effect of low dosage doxorubicin treatment in sensitizing cells to TSP1 stimulation. This model framework can be ultimately used to generate and optimize TSP1-based therapeutic strategies for promoting apoptosis.

Abbreviations

NF-κB:

Nuclear factor κB

cPARP:

Cleaved poly(ADP-ribose) polymerase

DXR:

Doxorubicin

FGF:

Fibroblast growth factor

MAPK:

Mitogen-activated protein kinase

ODE:

Ordinary differential equations

PARP:

Poly(ADP-ribose) polymerase

TSP1:

Thrombospondin-1

VEGF:

Vascular endothelial growth factor

WSSR:

Weighted sum of the squared residuals

XIAP:

X-linked inhibitor of apoptosis protein

References

  1. Bergers G, Benjamin LE. Angiogenesis: tumorigenesis and the angiogenic switch. Nat Rev Cancer. 2003;3:401–10.

    Article  CAS  PubMed  Google Scholar 

  2. Ren B, Yee KO, Lawler J, Khosravi-Far R. Regulation of tumor angiogenesis by thrombospondin-1. Biochim Biophys Acta. 2006;1765:178–88.

    CAS  PubMed  Google Scholar 

  3. Rusnati M, Urbinati C, Bonifacio S, Presta M, Taraboletti G. Thrombospondin-1 as a paradigm for the development of Antiangiogenic agents endowed with multiple mechanisms of action. Pharmaceuticals. 2010;3:1241–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Shibuya M. Differential roles of vascular endothelial growth factor receptor-1 and receptor-2 in angiogenesis. J Biochem Mol Biol. 2006;39:469–78.

    CAS  PubMed  Google Scholar 

  5. Zhang X, Kazerounian S, Duquette M, Perruzzi C, Nagy JA, Dvorak HF, et al. Thrombospondin-1 modulates vascular endothelial growth factor activity at the receptor level. FASEB J. 2009;23:3368–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Armstrong LC, Bornstein P. Thrombospondins 1 and 2 function as inhibitors of angiogenesis. Matrix Biol J Int Soc Matrix Biol. 2003;22:63–71.

    Article  CAS  Google Scholar 

  7. Kazerounian S, Yee KO, Lawler J. Thrombospondins: from structure to therapeutics. Cell Mol Life Sci. 2008;65:700–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Lawler J. Thrombospondin-1 as an endogenous inhibitor of angiogenesis and tumor growth. J Cell Mol Med. 2002;6:1–12.

    Article  CAS  PubMed  Google Scholar 

  9. Mirochnik Y, Kwiatek A, Volpert O. Thrombospondin and apoptosis: molecular mechanisms and use for design of complementation treatments. Curr Drug Targets. 2008;9:851–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Volpert OV. Modulation of endothelial cell survival by an inhibitor of angiogenesis thrombospondin-1: a dynamic balance. Cancer Metastasis Rev. 2000;19:87–92.

    Article  CAS  PubMed  Google Scholar 

  11. Isenberg JS, Frazier WA, Roberts DD. Thrombospondins: from structure to therapeutics. Cell Mol Life Sci. 2008;65:728.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Jiménez B, Volpert OV, Crawford SE, Febbraio M, Silverstein RL, Bouck N. Signals leading to apoptosis-dependent inhibition of neovascularization by thrombospondin-1. Nat Med. 2000;6:41–8.

    Article  PubMed  Google Scholar 

  13. Dawson DW, Pearce SFA, Zhong R, Silverstein RL, Frazier WA, Bouck NP. CD36 mediates the in vitro inhibitory effects of Thrombospondin-1 on endothelial cells. J Cell Biol. 1997;138:707–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Folkman J. Angiogenesis and apoptosis. Semin Cancer Biol. 2003;13:159–67.

    Article  CAS  PubMed  Google Scholar 

  15. Nor JE, Mitra RS, Sutorik MM, Mooney DJ, Castle VP, Polverini PJ. Thrombospondin-1 induces endothelial cell apoptosis and inhibits angiogenesis by activating the Caspase death pathway. J Vasc Res. 2000;37:209–18.

    Article  CAS  PubMed  Google Scholar 

  16. Githaka JM, Vega AR, Baird MA, Davidson MW, Jaqaman K, Touret N. Ligand-induced growth and compaction of CD36 nanoclusters enriched in Fyn induces Fyn signaling. J Cell Sci. 2016;129:4175–89.

    Article  CAS  PubMed  Google Scholar 

  17. Volpert OV, Zaichuk T, Zhou W, Reiher F, Ferguson TA, Stuart PM, et al. Inducer-stimulated Fas targets activated endothelium for destruction by anti-angiogenic thrombospondin-1 and pigment epithelium–derived factor. Nat Med. 2002;8:349–57.

    Article  CAS  PubMed  Google Scholar 

  18. Markovic SN, Suman VJ, Rao RA, Ingle JN, Kaur JS, Erickson LA, et al. A phase II study of ABT-510 (thrombospondin-1 analog) for the treatment of metastatic melanoma. Am J Clin Oncol. 2007;30:303–9.

    Article  CAS  PubMed  Google Scholar 

  19. Ebbinghaus S, Hussain M, Tannir N, Gordon M, Desai AA, Knight RA, et al. Phase 2 study of ABT-510 in patients with previously untreated advanced renal cell carcinoma. Clin Cancer Res. 2007;13:6689–95.

    Article  CAS  PubMed  Google Scholar 

  20. Kholodenko B, Yaffe MB, Kolch W. Computational approaches for analyzing information flow in biological networks. Sci Signal. 2012;5:re1.

    Article  PubMed  Google Scholar 

  21. Finley SD, Chu L-H, Popel AS. Computational systems biology approaches to anti-angiogenic cancer therapeutics. Drug Discov Today. 2015;20:187–97.

    Article  CAS  PubMed  Google Scholar 

  22. Rohrs JA, Sulistio CD, Finley SD. Predictive model of thrombospondin-1 and vascular endothelial growth factor in breast tumor tissue. Npj Syst Biol Appl. 2016;2:16030.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Jain H, Jackson T. Mathematical modeling of cellular cross-talk between endothelial and tumor cells highlights counterintuitive effects of VEGF-targeted therapies. Bull Math Biol. 2017. doi: 10.1007/s11538-017-0273-6.

  24. Albeck JG, Burke JM, Spencer SL, Lauffenburger DA, Sorger PK. Modeling a snap-action, variable-delay switch controlling extrinsic cell death. PLOS Biol. 2008;6:e299.

    Article  PubMed Central  Google Scholar 

  25. Neumann L, Pforr C, Beaudouin J, Pappa A, Fricker N, Krammer PH, et al. Dynamics within the CD95 death-inducing signaling complex decide life and death of cells. Mol Syst Biol. 2010;6:352.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Imoukhuede PI, Popel AS. Quantitative fluorescent profiling of VEGFRs reveals tumor cell and endothelial cell heterogeneity in breast cancer xenografts. Cancer Med. 2014;3:225–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Rubin DB, Drab EA, Bauer KD. Endothelial cell subpopulations in vitro: cell volume, cell cycle, and radiosensitivity. J Appl Physiol. 1989;67:1585–90.

    CAS  PubMed  Google Scholar 

  28. Birtwistle MR, Rauch J, Kiyatkin A, Aksamitiene E, Dobrzyński M, Hoek JB, et al. Emergence of bimodal cell population responses from the interplay between analog single-cell signaling and protein expression noise. BMC Syst Biol. 2012;6:109.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Pearce SFA, Wu J, Silverstein RL. Recombinant GST/CD36 fusion proteins define a Thrombospondin binding domain EVIDENCE FOR a SINGLE CALCIUM-DEPENDENT BINDING SITE ON CD36. J Biol Chem. 1995;270:2981–6.

    Article  CAS  PubMed  Google Scholar 

  30. Saltelli A, Bolado R. An alternative way to compute Fourier amplitude sensitivity test (FAST). Comput Stat Data Anal. 1998;26:445–60.

    Article  Google Scholar 

  31. Chen S, Guo X, Imarenezor O, Imoukhuede PI. Quantification of VEGFRs, NRP1, and PDGFRs on endothelial cells and fibroblasts reveals serum, intra-family Ligand, and cross-family Ligand regulation. Cell Mol Bioeng. 2015;8:383–403.

    Article  CAS  Google Scholar 

  32. Finley SD, Dhar M, Popel AS. Compartment model predicts VEGF secretion and investigates the effects of VEGF trap in tumor-bearing mice. Front Oncol. 2013;3 doi: 10.3389/fonc.2013.00196.

  33. Gaddy TD, Wu Q, Arnheim AD, Finley SD. Mechanistic modeling quantifies the influence of tumor growth kinetics on the response to anti-angiogenic treatment. PLoS Computat Biol. 2017. In press.

  34. Quesada AJ, Nelius T, Yap R, Zaichuk TA, Alfranca A, Filleur S, et al. In vivo upregulation of CD95 and CD95L causes synergistic inhibition of angiogenesis by TSP1 peptide and metronomic doxorubicin treatment. Cell Death Differ. 2005;12:649–58.

    Article  CAS  PubMed  Google Scholar 

  35. Eckelman BP, Salvesen GS, Scott FL. Human inhibitor of apoptosis proteins: why XIAP is the black sheep of the family. EMBO Rep. 2006;7:988–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Chawla-Sarkar M, Bae SI, Reu FJ, Jacobs BS, Lindner DJ, Borden EC. Downregulation of Bcl-2, FLIP or IAPs (XIAP and survivin) by siRNAs sensitizes resistant melanoma cells to Apo2L/TRAIL-induced apoptosis. Cell Death Differ. 2004;11:915–23.

    Article  CAS  PubMed  Google Scholar 

  37. Liu S, Zhang P, Chen Z, Liu M, Li X, Tang H. MicroRNA-7 downregulates XIAP expression to suppress cell growth and promote apoptosis in cervical cancer cells. FEBS Lett. 2013;587:2247–53.

    Article  CAS  PubMed  Google Scholar 

  38. Tong Q-S, Zheng L-D, Wang L, Zeng F-Q, Chen F-M, Dong J-H, et al. Downregulation of XIAP expression induces apoptosis and enhances chemotherapeutic sensitivity in human gastric cancer cells. Cancer Gene Ther. 2005;12:509–14.

    Article  CAS  PubMed  Google Scholar 

  39. Hammer E, Bien S, Salazar MG, Steil L, Scharf C, Hildebrandt P, et al. Proteomic analysis of doxorubicin-induced changes in the proteome of HepG2cells combining 2-D DIGE and LC-MS/MS approaches. Proteomics. 2010;10:99–114.

    Article  CAS  PubMed  Google Scholar 

  40. Haagenson KK, Wu GS. The role of MAP kinases and MAP kinase phosphatase-1 in resistance to breast cancer treatment. Cancer Metastasis Rev. 2010;29:143–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Grossi V, Peserico A, Tezil T, Simone C. p38α MAPK pathway: a key factor in colorectal cancer therapy and chemoresistance. World J Gastroenterol. 2014;20:9744–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Samples J, Willis M, Klauber-DeMore N. Targeting angiogenesis and the tumor microenvironment. Surg Oncol Clin N Am. 2013;22:629–39.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Pedersen AK, Mendes Lopes de Melo J, Mørup N, Tritsaris K, Pedersen SF. Tumor microenvironment conditions alter Akt and Na+/H+ exchanger NHE1 expression in endothelial cells more than hypoxia alone: implications for endothelial cell function in cancer. BMC Cancer. 2017;17:542.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Cheng H-W, Chen Y-F, Wong J-M, Weng C-W, Chen H-Y, Yu S-L, et al. Cancer cells increase endothelial cell tube formation and survival by activating the PI3K/Akt signalling pathway. J Exp Clin Cancer Res. 2017;36:27.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Chaitanya GV, Alexander JS, Babu PP. PARP-1 cleavage fragments: signatures of cell-death proteases in neurodegeneration. Cell Commun Signal CCS. 2010;8:31.

    Article  CAS  PubMed  Google Scholar 

  46. Altschuler SJ, Wu LF. Cellular heterogeneity: when do differences make a difference? Cell. 2010;141:559–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Friedman N, Cai L, Xie XS. Linking stochastic dynamics to population distribution: an analytical framework of gene expression. Phys Rev Lett. 2006;97:168302.

    Article  PubMed  Google Scholar 

  48. Meyer R, D’Alessandro LA, Kar S, Kramer B, She B, Kaschek D, et al. Heterogeneous kinetics of AKT signaling in individual cells are accounted for by variable protein concentration. Front Physiol. 2012;3 doi: 10.3389/fphys.2012.00451.

  49. Iwamoto K, Shindo Y, Takahashi K. Modeling cellular noise underlying heterogeneous cell responses in the epidermal growth factor signaling pathway. PLoS Comput Biol. 2016;12:e1005222.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Kellogg RA, Tay S. Noise facilitates transcriptional control under dynamic inputs. Cell. 2015;160:381–92.

    Article  CAS  PubMed  Google Scholar 

  51. Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature. 2009;459:nature08012.

    Google Scholar 

  52. Holcik M, Yeh C, Korneluk RG, Chow T. Translational upregulation of X-linked inhibitor of apoptosis (XIAP) increases resistance to radiation induced cell death. Nat Oncogene. 2000;19 doi: 10.1038/sj.onc.1203765.

  53. Virág L, Szabó C. The therapeutic potential of poly(ADP-ribose) polymerase inhibitors. Pharmacol Rev. 2002;54:375–429.

    Article  PubMed  Google Scholar 

  54. Bazzazi H, Isenberg JS, Popel AS. Inhibition of VEGFR2 activation and its downstream signaling to ERK1/2 and calcium by Thrombospondin-1 (TSP1): in silico investigation. Front Physiol. 2017;8 doi: 10.3389/fphys.2017.00048.

  55. Imoukhuede PI, Popel AS. Quantification and cell-to-cell variation of vascular endothelial growth factor receptors. Exp Cell Res. 2011;317:955–65.

    Article  CAS  PubMed  Google Scholar 

  56. Bentele M, Lavrik I, Ulrich M, Stößer S, Heermann DW, Kalthoff H, et al. Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis. J Cell Biol. 2004;166:839–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

The authors thank members of the Finley research group for critical comments and suggestions and Dr. James Finley for insight regarding the ROC analysis. The authors thank Dr. P.I. Imoukhuede and Si Chen from the University of Illinois, Urbana Champaign for training and insight regarding receptor quantification.

Funding

The authors acknowledge the support of the US National Science Foundation (CAREER Award 1552065). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files.

Author information

Authors and Affiliations

  1. Department of Biomedical Engineering, University of Southern California, Los Angeles, California, USA

    Qianhui Wu & Stacey D. Finley

  2. Department of Chemical Engineering and Materials Science, University of Southern California, Los Angeles, California, USA

    Stacey D. Finley

Authors
  1. Qianhui Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Stacey D. Finley
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

SDF designed the research. QW trained model to experimental data and generated model simulations and analyses. SDF performed the ROC analysis. Both authors contributed to writing the manuscript and have read and approved the final manuscript.

Corresponding author

Correspondence to Stacey D. Finley.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Table S1.

Biochemical reactions and parameter values. Table S2. Initial protein concentrations. Table S3. Estimated parameter values and model fitness. Table S4. Results from ANOVA analysis for T t and max cPARP. Table S5. Results from non-parametric t test for XIAP and PARP. Table S6. Results from ROC analysis. (XLSX 32 kb)

Additional file 2:

Ordinary differential equations and initial conditions (PDF 47 kb)

Additional file 3: Figure S1.

Comparison of minimal model of FasL signaling to experimental data. Figure S2. Global sensitivity analysis of parameters. Figure S3. Global sensitivity analysis of non-zero initial concentrations. Figure S4. Receptor number distributions of CD36 and Fas. Figure S5. Distribution of estimated parameter values. Figure S6. Dynamics of all species in baseline model. Figure S7. Population-level response to 10 nM TSP1 stimulation. Figure S8. Population-level response to 0.1 nM TSP1 stimulation. (PDF 13089 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wu, Q., Finley, S.D. Predictive model identifies strategies to enhance TSP1-mediated apoptosis signaling. Cell Commun Signal 15, 53 (2017). https://doi.org/10.1186/s12964-017-0207-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12964-017-0207-9

Keywords

Cell Communication and Signaling

ISSN: 1478-811X