The impact of tumor receptor heterogeneity on the response to anti-angiogenic cancer treatment

Abstract

Multiple promoters and inhibitors mediate angiogenesis, the formation of new blood vessels, and these factors represent potential targets for impeding vessel growth in tumors. Vascular endothelial growth factor (VEGF) is a potent angiogenic factor targeted in anti-angiogenic cancer therapies. In addition, thrombospondin-1 (TSP1) is a major endogenous inhibitor of angiogenesis, and TSP1 mimetics are being developed as an alternative type of anti-angiogenic agent. The combination of bevacizumab, an anti-VEGF agent, and ABT-510, a TSP1 mimetic, has been tested in clinical trials to treat advanced solid tumors. However, the patients’ responses are highly variable and show disappointing outcomes. To obtain mechanistic insight into the effects of this combination anti-angiogenic therapy, we have constructed a novel whole-body systems biology model including the VEGF and TSP1 reaction networks. Using this molecular-detailed model, we investigated how the combination anti-angiogenic therapy changes the amounts of pro-angiogenic and anti-angiogenic complexes in cancer patients. We particularly focus on answering the question of how the effect of the combination therapy is influenced by tumor receptor expression, one aspect of patient-to-patient variability. Overall, this model complements the clinical administration of combination anti-angiogenic therapy, highlights the role of tumor receptor variability in the heterogeneous responses to anti-angiogenic therapy, and identifies the tumor receptor profiles that correlate with a high likelihood of a positive response to the combination therapy. Our model provides novel understanding of the VEGF–TSP1 balance in cancer patients at the systems-level and could be further used to optimize combination anti-angiogenic therapy.

---

Insight, innovation, integration

This work reports a novel multi-compartment model integrating the reaction networks of VEGF and TSP1, two potent angiogenic factors, in human cancer patients. With the model, we gain new insight into the impact of tumor receptor heterogeneity on the response to a clinically tested combination anti-angiogenic therapy targeting both VEGF and TSP1 signaling and we identify potential tissue biomarkers. The model predictions provide mechanistic explanations of several important results reported in clinical trials of anti-angiogenic therapy, including the heterogeneous response to treatment and the importance of VEGFR1 as a predictive biomarker. Overall, this computational framework can be applied to improve combination anti-angiogenic therapy by increasing our understanding of the heterogeneous response, identifying potential biomarkers, and optimizing the pre-clinical and clinical studies.

Introduction

Angiogenesis is a hallmark of cancer that facilitates tumor progression in many aspects.¹ Tumor growth relies on the formation of new blood vessels to enable waste exchange and to deliver oxygen and nutrients to the tumor. In addition, angiogenesis increases the likelihood of metastasis by enabling tumor cells to enter the bloodstream and disperse to other sites in the body.² Considering the outstanding importance of angiogenesis for tumor development, anti-angiogenic therapy was designed to starve the tumor of its nutrient supply and limit its growth.³

Tumor angiogenesis is controlled by both pro- and anti-angiogenic signaling.^4–6 Common anti-angiogenic therapy uses single agents to reduce the pro-angiogenic signals, and the primary anti-angiogenic agent being used in the clinic inhibits signaling mediated by vascular endothelial growth factor-A (VEGF), a potent promoter of angiogenesis. However, this approach is not effective in all cancers. For example, bevacizumab, a monoclonal antibody that binds VEGF, is no longer approved for the treatment of metastatic breast cancer because it was not shown to be effective and safe for patients.⁷ Sunitinib, a tyrosine kinase inhibitor that targets VEGF receptors and other growth factor receptors, has also shown limited success.⁸ These limitations of anti-VEGF treatment prompt the need to optimize anti-angiogenic therapy. One alternate approach is to enhance the signal of anti-angiogenic factors. Thrombospondin-1 (TSP1) is one of the most studied endogenous inhibitors of angiogenesis and has been shown to inhibit vascular growth and tumorigenesis in preclinical trials.^9–11 Inspired by the effect of TSP1, TSP1 mimetics were developed for tumor treatment.¹² One such drug, named ABT-510, reached Phase II clinical trials. However, ABT-510 failed to show clear evidence of efficacy and is no longer tested as a single-agent drug in clinical development.^13,14

The disappointing outcomes of clinical studies of anti-angiogenic drugs as single agents prompt the development of combination anti-angiogenic therapy. Administering a combination of anti-angiogenic agents that simultaneously target multiple angiogenic signals is expected to achieve efficient and durable suppression of angiogenesis by strongly shifting the relative balance of inducers and inhibitors of angiogenesis to oppose the “angiogenic switch”.^15,16 The combination of agents targeting different pathways might also prevent tumors from leveraging complementary pathways to escape anti-angiogenic treatment. With this in mind, the ABT-510 TSP1 mimetic was clinically tested in combination with bevacizumab in patients with advanced solid tumors. However, patients displayed a heterogeneous response to this combination therapy:¹⁷ one patient had a partial response and only 32% of the patients had prolonged stable disease (≥6 months). Unfortunately, the mechanisms driving these disappointing results were not elucidated in the trial. Three fundamental questions remain: do the levels of TSP1 and VEGF balance one another in tumor tissue, how does combination therapy influence this VEGF–TSP1 balance, and does inter-patient heterogeneity lead to significantly different responses. Answering these questions contributes to our understanding of the action of combination therapy. In addition, considering the angiogenic balance at the levels of tissue, organs and the whole body can help us optimize anti-angiogenic therapy. In this study, we address these questions using a computational systems biology model.

Mathematical modeling serves as a useful tool to study the response to anti-angiogenic treatment, complementing experimental and clinical studies. Various modeling approaches have been applied to study anti-angiogenesis therapies, including differential equation based models, Boolean network models, multi-compartment models, hybrid cellular automaton models, multiscale agent-based models, image-based models and bioinformatics-based modeling.¹⁸ In particular, the effects of anti-VEGF agents on VEGF and its receptors have been intensively investigated with computational modeling. Targeting VEGF binding, VEGF secretion and VEGFR signaling as anti-angiogenic strategies have been illustrated in different studies.¹⁸ Recently, mathematical modeling was used to understand the impact of the cross-talk between tumor cells and endothelial cells¹⁹ and the vascular phenotypes²⁰ on the effects of anti-angiogenic treatment, which generates new insights of anti-angiogenic therapy using mathematical modeling. Each model provides quantitative insight into the response to various anti-angiogenic strategies, at different levels of detail and scales, ranging from intracellular signaling to whole-body dynamics. However, these existing models do not provide information regarding the balance of pro- and anti-angiogenic signals and cannot be used to study the mechanistic effects of combination therapy targeting both sides of the angiogenic balance. In this study, one of our goals is to understand the anti-angiogenic therapy targeting VEGF and TSP1 simultaneously, which has not been the focus of previous models. To achieve this goal, we build a novel ordinary differential equation (ODE) based multi-compartment model.

The model presented here significantly builds upon our previous modeling efforts, and is mainly based on two published models. One model is the whole-body model of the VEGF–receptor system, which was previously used to illustrate the counterintuitive increasing of VEGF after anti-VEGF treatment.²¹ Second, we build on our recent model of TSP1 and VEGF interactions in tumor tissue, which predicts the effects of various strategies mimicking TSP1's anti-angiogenic properties.²² This model of the VEGF–TSP1 balance in tumor tissue, however, omits the trafficking of soluble species and only considers drugs delivered directly into the tumor. Here, we significantly expand these two previous works^22,23 to generate a novel whole-body model of the VEGF–TSP1 interaction network. The expanded model has three compartments to incorporate the drug pharmacokinetics (PK) and species distribution in the human body. The model also incorporates pharmacodynamics (PD). Since the binding of angiogenic factors to their receptors triggers a cascade of intracellular reactions, including phosphorylation of the receptors, we use the number of ligand–receptor complexes as an approximation of the receptor activation level to capture the status of pro- and anti-angiogenic signaling in tissue. Altogether, our new model enables a complete PK/PD study of the clinically tested combination anti-angiogenic therapy targeting both VEGF and TSP1.

We apply the model to understand how the angiogenic balance of VEGF and TSP1 is modulated by bevacizumab and ABT-510 combination therapy. Then we use the model to investigate the impact of inter-patient heterogeneity, specifically the tumor receptor heterogeneity, on the response to combination anti-angiogenic therapy. Compared to other inter-patient variability, tumor receptor heterogeneity is one of the most well supported in published literature. It has been observed experimentally that tissue samples from patients with different types of cancer and different samples from patients with the same type of cancer have different receptor expressions.^24–27 Additionally, the VEGFR2 heterogeneity was shown to affect the response to an anti-angiogenic cyclophosphamide treatment in an in vitro experimental setting.²⁸ The patient-to-patient VEGFR1 and neuropilin variability has been associated with the response to bevacizumab treatment as intra-tumoral biomarkers.^29,30 Thus, understanding the effects of tumor receptor variability is clinically relevant.^28–31 Our model predicts that the metric of angiogenic receptor expression can serve as a predictive tissue biomarker to distinguish the patients in which the combination anti-angiogenic therapy will elicit a strong therapeutic response. Overall, we establish a new computational framework to predict the effects of anti-angiogenic therapies and understand clinical observations.

Methods

The compartmental whole-body model

This mechanistic model characterizes the extracellular distribution of angiogenic species in the human body. We follow the compartmental model structure used in previous works.^23,32 In this approach, a tissue is assumed to be a collection of capillaries, surrounded by parenchymal cells. The interstitial space lies between the parenchymal cells and the capillaries, and is comprised of the extracellular matrix (ECM), parenchymal basement membranes (PBM) and endothelial basement membranes (EBM). The soluble species are assumed to diffuse within the available interstitial space very fast compared to the timescale of the biochemical reactions,³³ thus all of the structures are modeled in a spatially-averaged manner as a simplification. A human cancer patient is represented by a three-compartment model: normal tissue (“normal,” represented by skeletal muscle), the vasculature (“blood”), and diseased tissue (“tumor”) (Fig. 1A). The soluble species are introduced to the system by being secreted by cells and are removed from the system through degradation and clearance from the blood. Receptors are uniformly distributed on the cell surfaces and can be internalized by the cell and recycled back to the surface. It is largely unknown how the expression level of angiogenic receptors changes over time during anti-angiogenic treatment. Additionally, the timescale over which receptor levels typically vary (on the order of minutes) is much shorter than the timescale over which we are simulating in the model (days to weeks). Therefore, we assumed the total number for each type of receptor is conserved at every simulated time point (the receptor recycling rate is the same as the internalization rate). Transport of soluble species between the compartments is mediated by transcapillary permeability and lymphatic flow.

Fig. 1 Compartmental model of VEGF–TSP1 system in human cancer patients and TSP1 and VEGF interactions. (A) The model includes three compartments: normal tissue, blood, and tumor tissue. The endothelial cells and parenchymal cells secrete soluble species (VEGF, TSP1, MMP3, and proMMP9) into the compartments. Receptors are localized on the luminal and abluminal surfaces of endothelial and parenchymal cells. Free and ligand-bound receptors can be internalized. Transport between compartments occurs by trans-endothelial permeability and lymphatic flow. Soluble species are degraded in the tissue or cleared from the blood. Specifically, the compartment model includes: (B) the molecular interactions of two active VEGF isoforms (VEGF₁₂₁ and VEGF₁₆₅), receptors (R1 and R2), and co-receptors (N1 and N2); (C) the interactions between TSP1, its receptors (CD36, CD47, β1, and LRP1), VEGF and VEGFR2; (D) the activation of proMMP9 via cleavage by MMP3; cleavage of TSP1; and proteolysis of VEGF₁₆₅ (free or bound to glycosaminoglycan, GAG, chains) to form VEGF₁₁₄ by active MMPs; (E) the sequestration of VEGF₁₆₅ and TSP1 by GAG chains in the extracellular matrix and the cellular basement membranes; (F) the sequestration of VEGF by α2M in blood; (G) the binding between ABT-510 and receptors and sequestration of VEGF by bevacizumab.

Model construction is based on certain assumptions. Firstly, we formulate the compartmental model assuming that the tumor volume is constant. Admittedly, there is a change in the number of healthy or diseased cells in human patients undergoing anti-angiogenic therapy, as a primary goal of treatment is to reduce tumor volume. However, as we have shown in our previous work,³⁴ that since the tumor is nearly 2000 times smaller than the normal compartment, the tumor volume must change by at least two orders of magnitude (to ∼3300 cm³) for it to significantly influence the distributions of the soluble factors, which is a central focus of our PK/PD compartment model. Since this size of tumor is not physiologically realistic, we assume constant tumor volume and instead focus on the distribution of the soluble factors and the formation of pro- and anti-angiogenic complexes in each of the compartments. Accordingly, we assume that the total surface area of the microvessels is constant, as the tissue vascularity is characterized as the ratio of the microvascular surface area to the tumor volume. Finally, we assume that the vascular permeability between compartments is fixed. We implement this simplification because there is a scarcity of quantitative data available to formulate a mathematical equation to capture the relationship between the anti-angiogenic therapy and vascular permeability. Rather than impose further uncertainty in the model, we maintain a constant value for the permeability.

Rule-based model of VEGF–TSP1 reaction network

The characterization of the species’ dynamics is based on the principles of mass action kinetics and biological transport. BioNetGen, a rule-based modeling approach, is used to construct the model.³⁵ The biological reaction rules are defined in BioNetGen, which automatically generates the set of ordinary differential equations (ODEs) that describe how the species’ concentrations evolve over time. A detailed description of the derivation of the ODE model through BioNetGen is documented in Supplementary file S1 (ESI†), which provides an explanation of the propagation of the rules and reactions in the model generation. Here, we briefly summarize the defined rules that govern the molecular interactions and corresponding reactions (Fig. 1B–G). Following our previous works, the model includes two active VEGF isoforms (VEGF₁₆₅ and VEGF₁₂₁). The inactive form, VEGF₁₁₄, is the product of proteolytic cleavage of VEGF₁₆₅. Two predominant VEGF receptors, VEGFR1 and VEGFR2 (R1 and R2), and neuropilin co-receptors, NRP1 and NRP2 (N1 and N2), are considered (Fig. 1B). TSP1 binds to its receptors, CD36, CD47, low density lipoprotein receptor-related protein 1 (LRP1) and α_xβ₁ intergrins (β₁, a generic form representing several species) (Fig. 1C). We also include matrix metalloproteinase species (MMP3, MMP9 and proMMP9), which promote VEGF cleavage. TSP1 impedes the activation of MMP9 (Fig. 1D) as a means of inhibiting pro-angiogenic signaling. Glycosaminoglycan (GAG) chains reside in the interstitial space, representing the extracellular matrix, as well as in the cellular basement membranes. GAG chains are able to bind and sequester TSP1 and VEGF₁₆₅ (Fig. 1E). The α-2-macroglobulin (α2M) species, a protease inhibitor, is confined to the blood compartment, where it binds to VEGF (Fig. 1F).

Model parameterization

There are 157 parameters presented in our model, including geometric parameters, kinetic parameters, receptor numbers, secretion and degradation rates, transport rates, and parameters for the drug properties. The model parameter values are reported in Supplementary file S1 (ESI†), with literature references, and are described below.

Geometric parameters (27 parameters). The geometric parameters characterize the fundamental structure of the model, defining the volume of compartments, the interstitial space volume, and tissue surface areas of endothelial and parenchymal cells. These parameters are based on experimental measurements taken directly from in vivo mouse tumor models.³² We assume that the geometric characteristics of xenograft tumors in mice recapitulate human tumors, rather than relying on data from in vitro cell culture. The set of geometric parameters has been used in multiple previous studies,^21,23,34 and we adopt the parameters without changing their values.

Kinetic parameters (47 parameters). The kinetic parameters specify the association and dissociation rates for the binding of molecular species. For the VEGF axis, the kinetic parameters are based on experimental measurements for the biochemical interactions of VEGF and its receptors, which have been implemented in previous models.^23,36 Likewise, the kinetic parameters for TSP1 axis are based on experimental measurements that estimate the rates of interactions between TSP1 and its receptors and other binding partners, which were systematically reported in our published model of VEGF and TSP1 in tumor tissue.²²

Receptor numbers (32 parameters). The receptor numbers for the VEGF axis (VEGFR1, VEGFR2 and neuropilin-1 and -2) are taken from quantitative measurements of receptor expression on cells from mouse xenograft studies. The cell surface expression of these receptors was measured via flow cytometry.^24,25 For the TSP1 receptor numbers, we referred to the qualitative measurements reported in the Human Protein Atlas,²⁷ assuming “high”, “medium”, and “low” expression levels correspond to 10 000, 5000, and 2500 receptors per cell, respectively.

Secretion and clearance rates (34 parameters). The clearance and degradation rates are based on the reported protein half-life values. The secretion rates of VEGF were fit based on modeling in vivo population PK data in our previous study.²³ The secretion rates of TSP1, MMP3 and proMMP9 (10 parameters) are fitted in this study to match the experimental measurements shown in Table 1.

Table 1 Comparison of predicted and experimental steady state concentrations of VEGF, TSP1 and MMPs

Species	Range of experimental measurements^a	Predicted concentration	Source and references
a Median value is shown in parentheses, if provided in literature. If the experimental data reflects the total concentration in tissue, we assume 50% of the total protein amount is localized in extracellular space. b TSP1 concentration includes both the active TSP1 and the cleaved TSP1.
Tumor tissue
VEGF	8.0–389 pM	148.6 pM	Multiple cancer types^34
TSP1^b	1.0–6.2 nM (2.0)	2.0 nM	Breast cancer patients^71
MMP3	1.8–65.1 nM (5.1)	4.9 nM	Oral squamous cell carcinoma^72
MMP9 total	1.0–287.8 nM (9.0)	9.3 nM	Oral squamous cell carcinoma^72
MMP9 active	0–22.4 nM (0.8)	0.2 nM	Oral squamous cell carcinoma^72
Plasma
VEGF	0.4–3.0 pM	1.6 pM	Multiple cancer types^73
TSP1^b	0.8–2.1 nM (1.2)	1.2 nM	Breast cancer patients^71
MMP3	1.9–2.0 nM	1.9 nM	Breast cancer patients^74
MMP9 total	0.6–0.7 nM	0.7 nM	Breast cancer patients^74
Normal tissue
VEGF	0.3–3.0 pM	1.0 pM	Healthy subjects^34
TSP1^b	0.2 nM	0.9 nM	Healthy subjects^75
MMP3	0.9–66.4 nM (4.1)	4.2 nM	Oral squamous cell carcinoma^72
MMP9 total	0.8–27.7 nM (5.0)	3.2 nM	Oral squamous cell carcinoma^72
MMP9 active	0–4.1 nM (0.02)	0.05 nM	Oral squamous cell carcinoma^72

---

Transport rates (6 parameters). We assume passive transport for all soluble species in our model. That is, the transport of angiogenic species between compartments is only mediated by lymphatic flow and vascular permeability. The rate of lymphatic flow is based on the reported experimental measurements.³⁷ The vascular permeability of VEGF is determined in our previous work,³² by considering the Stokes–Einstein radius for the VEGF protein. Since other angiogenic factors are similar in size to VEGF,^38,39 we assume that all of the newly introduced angiogenic species have the same vascular permeability as VEGF. This assumption can be relaxed in future work.

Properties of the anti-angiogenic drugs (11 parameters). The degradation and clearance rates of bevacizumab and ABT-510 are converted from their reported half-life values measured in clinical trial. The vascular permeability of bevacizumab and its binding rates are the same as those used in our previous modeling work,³⁴ which capture the clinically-measured pharmacokinetic data⁴⁰ (Fig. S1A, ESI†). The vascular permeability of ABT-510 was set to be the same as VEGF. Since ABT-510 is a TSP1-derived peptidomimetic that specifically binds to the CD36 receptor, we assume it has the same affinity to CD36 as TSP1. The bioavailability of ABT-510 is tuned to be 30% in subcutaneous injection,⁴¹ in order to match pharmacokinetic data of TSP1⁴² (Fig. S1B, ESI†).

In summary, parameters are either taken from our previous modeling studies or estimated based on experimental measurements. Only 11 parameters in total are fitted in this study (the secretion rates of TSP1, MMP3 and proMMP9 in the compartments and the ABT-510 bioavailability). Given the presence of the uncertainty of the parameters, we performed a global sensitivity analysis to understand the robustness of the baseline model predictions (see below).

Sensitivity analysis. We perform the extended Fourier amplitude sensitivity test (eFAST), a global variance-based sensitivity analysis, to understand how the uncertainty of parameters (“model inputs”) affects the baseline model predictions (“model outputs”).⁴³ We analyzed the effects of three groups of parameters (receptor numbers, kinetic parameters, and vascular permeability) on nine different model outputs (the concentrations of TSP1, VEGF, proMMP9, MMP9, and MMP3; the TSP1–VEGF, proMMP9–MMP3, and MMP3–TSP1 complexes; and the angiogenic ratio) in each of the three compartments. In each case, the parameter values were allowed to vary 10-fold above and below the baseline values (a total range of two orders of magnitude) to account for uncertainty in the model parameters. In the eFAST method, the inputs are varied together, at different frequencies. The Fourier transform of the outputs is calculated to identify the influence of each parameter, based on the amplitude of each input's frequency. Two different sensitivity indices are generated in the eFAST analysis: the first-order FAST indices, S_i, and the total FAST indices, S_Ti. The first-order indices (S_i) measure the local sensitivity of individual inputs, while the total indices (S_Ti) represent the global sensitivity by accounting for second- and higher-order interactions between multiple inputs. The eFAST method is implemented using MATLAB code developed by Kirschner and colleagues.⁴³ We have performed this analysis to characterize the robustness of our previous models.^22,23

Simulating receptor variability. To investigate the impact of tumor receptor heterogeneity, we perform a Monte Carlo analysis. In this analysis, parameter values are randomly varied, by drawing from a defined distribution for each parameter. Here, we vary the receptor expression parameter values in the tumor compartment. In total, the densities of 16 receptors are varied in the simulations: four VEGF receptors on tumor cells (R1_Tum, R2_Tum, N1_Tum, N2_Tum), four TSP1 receptors on tumor cells (CD36_Tum, CD47_Tum, LRP1_Tum, β1_Tum), four VEGF receptors on tumor endothelial cells (R1_disEC, R2_disEC, N1_disEC, N2_disEC), and four TSP1 receptors on tumor endothelial cells (CD36_disEC, CD47_disEC, LRP1_disEC, β1_disEC). The densities of these receptors were randomly chosen from a uniform distribution within a range of 10-fold above and below the baseline value. We generated 1000 different combinations of receptor density profiles, representing 1000 unique cancer patients. We ran the model for each of the receptor profiles with anti-angiogenic treatment to examine how the response to treatment varies across the 1000 parameter sets.

Simulation of therapy

Administration of combination therapy. The combination therapy of bevacizumab and ABT-510 is simulated by mimicking the administration strategy used in clinical trials.¹⁷ We first allowed the model to reach steady state (this occurs within 24 hours) before the start of treatment. We then simulated one cycle of the combination therapy: bevacizumab was administered once at the beginning of the cycle; ABT-510 was administered every 12 hours for 14 days. Bevacizumab was given at a dose of 10 mg kg⁻¹ through intravenous infusion lasting 90 minutes, while ABT-510 was administered at 100 mg twice daily through subcutaneous injection. The bolus of ABT-510 was given directly to a subcutaneous compartment⁴⁴ (assumed to be a reservoir with a volume of 30 cm³), and it is subsequently transported into blood. The transportation between the subcutaneous and blood compartments is unidirectional and is assumed to occur at the same rate as the transport between the normal and blood compartments. ABT-510 binds to TSP1 receptor CD36 to induce an anti-angiogenic signal⁴⁵ (Fig. 1G). Bevacizumab is a VEGF antibody that sequesters VEGF to keep it from binding to its receptors, thereby inhibiting the pro-angiogenic signal (Fig. 1G).

Characterization of the response to treatment. In our study, the response to anti-angiogenic treatment is characterized based on the angiogenic balance in tumor tissue. We define the angiogenic balance as the ratio of the concentrations of the pro-angiogenic complexes to the anti-angiogenic complexes. The ratio indicates the activation level of the pro-angiogenic receptors relative to the activation level of anti-angiogenic receptors. Specifically, the pro-angiogenic complexes include the ligand-bound VEGF receptors that are not interacting with active TSP1 receptors. Here we assume a ligand-bound VEGF receptor coupled with active TSP1 receptor is not a pro-angiogenic complex, since the downstream signaling of the ligand-bound VEGF receptor could be inhibited by TSP1.^46,47 The anti-angiogenic complexes are the active TSP1 receptors, those bound to TSP1 or the TSP1 mimetic. The fold-change of the angiogenic ratio in the tumor compartment (F) characterizes the response to anti-angiogenic treatment:
where [VEGF:Rec] is the total concentration of pro-angiogenic complexes, [TSP1:Rec] is the total concentration of anti-angiogenic complexes, T denotes tumor compartment, and t indicates the simulated time: before treatment (t = 0 day) and after one cycle of treatment (t = 14 days).

Quantification of combination effect

We determine the combined effect of the two anti-angiogenic agents, bevacizumab and ABT-510. Commonly used dose–effect-based approaches to quantify drug combination effects rely on the mathematical framework known as Loewe Additivity.⁴⁸ It calculates the combination of dose a for drug X and dose b for drug Y that can produce the same effect as dose A of drug X alone and dose B of drug Y alone. The combination effect of drug X and drug Y can be expressed as:
where CI is the combination index. Additionally, reference drug values are usually set based on a desired end-point (i.e., the IC₅₀ value). However, here, there is no specified characteristic of the efficacy of the anti-angiogenic drugs (individually or in combination). Therefore, we selected a specific end point (a desired “fold-change of angiogenic ratio”, 0.75). This end point is firstly used to determine the drug doses a and b in combination. The ratio between doses a and b is a constant, equal to the ratio used in clinical protocol. Then, we change the doses of bevacizumab and ABT-510 to get the dose-effect curve of single-agent treatments and find doses A and B corresponding to the desired end point 0.75. A CI value less than one indicates synergy between the two drugs, a CI equal to one represents simple additivity, and CI greater than one indicates antagonism.⁴⁹

Partial least squares regression

Partial least squares regression (PLSR) relates input variables to output variables by maximizing the correlation between the variables. This is accomplished by projecting the output variables onto new dimensions (principal components, PCs), which are linear combinations of the inputs. We use PLSR to investigate the relationship between tumor receptor numbers and the response to anti-angiogenic therapy (the fold-change in the angiogenic ratio in the tumor tissue, termed F). The PLSR model is trained and validated with the 1000 different combinations of tumor receptor densities (inputs) and the corresponding responses to treatment generated by the three-compartment systems biology model (outputs). The receptor values are normalized by dividing by the lower bound of the sampling range. Since the receptor numbers are sampled from a uniform distribution within 10-fold above and below the baseline value, the normalized receptor values range from 1 to 100. The nonlinear iterative partial least squares (NIPALS) algorithm was used to implement the PLSR analysis.⁵⁰ Ten-fold cross validation is applied to quantify the predictive capability of the PLSR model. Furthermore, we calculated the “variable importance of projection” (VIP) for each input variable,⁵¹ which is the weighted sum of each input's contribution to the output. As such, the VIP evaluates the overall importance of an input variable in predicting the output. VIP values greater than one indicate variables that are important for predicting the output response.

Results

Baseline model is calibrated to experimental data

We report the species’ concentrations (VEGF, TSP1, MMP3 and MMP9) predicted by the baseline model and compare them with experimental measurements in Table 1. Our predictions quantitatively match the experimental data, where all of the predicted concentrations are within the measured ranges. We also compare the predicted pharmacokinetics of bevacizumab with the data collected in clinical trials (Fig. S1, ESI†) and confirm that our model can capture the dynamics of these two agents at different dosages.

We performed a global sensitivity analysis to quantify the robustness of the baseline model predictions. The eFAST method is used to calculate the sensitivity indices S_i and S_Ti (see Methods section). The S_i estimates how an input influences an output individually, while the S_Ti indicates how influential an input is in combination with other parameters. We find that the predicted species’ concentrations in the tumor compartment are mainly affected by the densities of receptors on tumor cells and diseased endothelial cells (Fig. S2A, ESI†). Additionally, the predicted concentrations for species in the blood and normal compartments are heavily affected by the receptor densities on normal endothelial cells and normal parenchymal cells. The vascular permeability significantly affects the model predictions (Fig. S2B, ESI†). Intuitively, the predicted concentrations in the tumor compartment are mainly affected by the permeability between tumor and vascular system, while the model predictions for blood and normal compartment are predominantly affected by permeability between normal tissue and the vascular system. The kinetic parameters are shared by reactions in all three compartments, thus they influence the predictions in all three compartments to varying degrees (Fig. S2C, ESI†).

Overall, the eFAST results indicate the parameters that should be fitted or for which additional experimental data is needed. Vascular permeability is known to be affected by VEGF and anti-angiogenic drugs. However, to our knowledge, there are no robust, quantitative measurements available that can be used to specify the mathematical relationship between VEGF or anti-angiogenic agents and permeability. Therefore, we keep the permeability constant, based on the size of the molecular species (see Methods). All of the kinetic parameter values, with the exception of the coupling rates between CD36 and VEGFR2, CD36 and β1, and CD47 and VEGFR2, are based on the reported experimental measurements. The three coupling rates for which we do not have experimental values are estimated to have little influence on the model predictions: kc_CD36:R2, kc_CD36:β1 and kc_CD47:R2 are either not significant or have very low S_i and S_Ti values (Fig. S2C, last three columns, ESI†). Therefore, the uncertainty of those parameters does not hamper the robustness of the predictions.

Overall, the sensitivity analysis demonstrates that the model is mostly affected by parameters whose values can be specified from experimental data. Below, we present results from applying the calibrated model to predict the effects of anti-angiogenic therapy. We explicitly study the effects of variability in the tumor receptor numbers, as those are shown to influence the predicted concentrations of the angiogenic factors, both in our sensitivity analysis and in experimental and clinical studies.^28–31

The angiogenic balance is predicted to shift during combination therapy

We apply the compartmental model to predict how the angiogenic signal changes with anti-angiogenic treatment. We simulate one cycle (14 days) of combination therapy, following the protocol used in clinical trial.¹⁷ Bevacizumab was given at a dose of 10 mg kg⁻¹ through intravenous infusion lasting 90 minutes. ABT-510 was administered at 100 mg every 12 hours for 14 days. The concentrations of unbound (“free”) VEGF and TSP1 over time, along with dynamics of the “angiogenic ratio”, are reported to show the change of angiogenic signals during combination therapy (Fig. 2). The angiogenic ratio is the number of pro-angiogenic complexes over the number of anti-angiogenic complexes (see Methods). If the angiogenic ratio is greater than one, the compartment is considered to be in a pro-angiogenic state.

Fig. 2 The effects of one cycle of combination anti-angiogenic therapy. The system is allowed to reach steady state before the start of treatment (Day 0). One cycle of treatment lasts 14 days: bevacizumab is given once every two weeks, and ABT-510 is given twice daily. (A) VEGF is shown to decrease after the treatment start in all three compartments. (B) TSP1 is predicted to stay at a stable level. (C) The angiogenic ratio (the ratio of pro-angiogenic complexes to anti-angiogenic complexes) is predicted to oscillate throughout the course of treatment.

After one cycle of the combination therapy, the concentration of free VEGF is predicted to decrease in all three compartments (Fig. 2A). The tumor VEGF level changes from 149 pM to 14 pM after one treatment cycle. The levels of free VEGF in plasma and normal tissue decreased from 1.6 pM to 1.2 pM and from 1.0 pM to 0.9 pM, respectively. In contrast, the free active TSP1 level in all three compartments is highly stable throughout the whole treatment cycle (Fig. 2B), remaining at 117 pM, 96 pM, and 66 pM in the tumor, blood, and normal compartments, respectively. This indicates that the combination therapy does not affect the free TSP1 level in human body. The experimental measurements from clinical study also show that the level of TSP1 does not significantly change,¹⁷ which validates this prediction. However, the combination therapy is predicted to significantly reduce the angiogenic ratio in all compartments. Before treatment, the angiogenic ratio (Fig. 2C) is predicted to highly favor angiogenesis in tumor and blood, where the ratios are 14.4 and 6.2, respectively. In comparison, the angiogenic ratio is almost balanced in the rest of the body (1.3 in normal compartment). The values of the angiogenic ratios at the end of the treatment cycle are 6.6, 2.0, and 0.6 in the tumor, blood, and normal compartments. Thus, following treatment, the angiogenic balance is significantly reduced, opposing angiogenesis in normal tissue, but still favoring angiogenesis in tumor tissue and blood. In all compartments, the angiogenic ratio is predicted to oscillate during the treatment. Since the ABT-510 drug has a short half-life (1.2 hours in circulation), it requires a high frequency of administration (twice-daily injections), which causes the angiogenic ratio to oscillate following each injection. Specifically, the angiogenic ratio decreases to a low level for a short time after the administration of ABT-510 and goes back to a medium level before the next dose of ABT-510.

To get detailed insight into how the angiogenic ratio varies in the tumor, we also report the change of all tumor angiogenic complexes in combination therapy and single agent treatment (Fig. S3, ESI†). VEGF can bind to four different receptors, including VEGF receptor 1 (R1), VEGF receptor 2 (R2), neuropilin-1 (N1) and neuropilin-2 (N2). Thus, four different types of pro-angiogenic complexes are formed. TSP1 binds to four different receptors, including CD47, CD36, LRP1 and β1, to produce four types of anti-angiogenic complexes. Our predictions show that bevacizumab promotes a decrease in the pro-angiogenic complex involving VEGFR1 (Fig. S3A, ESI†), and ABT-510 drives an increase in the anti-angiogenic complex involving CD36 (Fig. S3B, ESI†). Note that the ABT-510 drug has a short half-life (1.2 hours in circulation), requiring a high frequency of administration (twice-daily injections). Since ABT-510 binds to CD36 directly, the amount of ligand-bound CD36 goes up after the administration of ABT-510 and decreases to a low level before the next dose of ABT-510. This explains the pronounced fluctuations of CD36 during the treatment cycle. The absolute levels of other angiogenic complexes only show small changes. In Fig. S3C (ESI†), the curve shows aspects of the curves in Fig. S3A and B (ESI†), which indicates that the result of these two agents in combination involves each of the drug's individual effects.

The therapeutic response varies due to tumor receptor heterogeneity

To examine the impact of tumor receptor heterogeneity, we varied the number of VEGF and TSP1 receptors in the tumor tissue, both on tumor cells and diseased endothelial cells. Our sensitivity analysis shows that these receptor numbers significantly influence the tumor angiogenic ratio. Additionally, pre-clinical and clinical investigations indicate that tumor receptors affect anti-angiogenic treatment.^28–31 We set the density of each receptor by sampling within two orders of magnitude of the baseline value. This variability in the receptor number influences the dynamics of free VEGF, free TSP1, and the angiogenic ratio, as shown in Fig. S4 (ESI†). We use the “fold-change” of free VEGF, free TSP1, and the angiogenic ratio to characterize the response to combination treatment. The fold-change is the post-treatment level compared to the pre-treatment level. When the fold-change is less than one, the quantity has decreased after treatment. We further define a therapeutic response to be when the fold-change of free VEGF or the angiogenic ratio is less than one, indicating a shifting of the angiogenic balance to oppose angiogenesis. The model predicts that combination treatment elicits a therapeutic effect for all of the receptor levels simulated. That is, the fold-change of free VEGF and the angiogenic ratio are predicted to be lower than one in all three compartments (Fig. 3A), indicating that free VEGF and the angiogenic ratio decrease after treatment.

Fig. 3 The impact of tumor receptor heterogeneity on the response to anti-angiogenic therapy. We simulated 1000 tumor receptor profiles and predict: (A) the fold-changes of the angiogenic ratio, VEGF, and TSP1 in tumor (red), blood (green), and normal (blue) compartments. (B) The distributions of the fold-change of the tumor angiogenic ratio (F) for single-agent treatments and combination therapy (green: ABT-510 alone; purple: bevacizumab alone; orange: combination treatment).

Model predictions show that the response level of the fold-change in tumor varies depending on the tumor receptor profile. Different combinations of sampled parameter values represent patients with different tumor receptor profiles. The model predicts a wide range of fold-changes for the angiogenic ratio and free VEGF in the tumor compartment (Fig. 3A), which indicates that different tumor receptor profiles lead to different responses to the combination therapy. Furthermore, the skewed distribution of the fold-changes of the tumor angiogenic ratio (F) (Fig. 3B, orange) shows that the majority (83%) of the sampled tumor profiles lead to intermediate responses with fold-change ranging from 0.5 to 1. A smaller percentage of profiles (17%) have a strong response to treatment (F < 0.5). Interestingly, there are some outliers for which combination treatment significantly decreases the pro-angiogenic signal, with fold-changes of free VEGF or angiogenic ratio much lower than one. These results imply that particular tumor receptor profiles could have a strong therapeutic response to combination therapy. In single-agent anti-angiogenic tumor therapy (ABT-510 or bevacizumab treatment alone), heterogeneous responses are also observed (Fig. 3B, green and purple). Thus, the response to single-agent treatment also varies depending on the tumor receptor profile.

Tumor receptor heterogeneity influences treatment outcomes in other aspects. Firstly, variability in the tumor receptor numbers could affect the targeting of combination anti-angiogenic therapy. Fig. 3A shows the fold-changes in the angiogenic ratio in the blood and normal compartments are ∼0.3 and ∼0.45, respectively. In comparison, in the tumor, the fold-change varies from ∼0.3–0.9, with the median being 0.6. That is, across all 1000 simulations, the median effect of the combination therapy is a reduction of the angiogenic ratio by 40% in tumor, while the angiogenic ratio is reduced by 70% in blood and by 55% in normal tissue. Thus, for most cases, the fold-change of the angiogenic ratio in healthy tissue and in blood is lower than in the tumor (Fig. 3A). This indicates the combination therapy is not always efficiently targeting tumor angiogenesis. That is, the combination therapy is targeting physiological angiogenesis instead of tumor angiogenesis. This insight is not captured by solely examining the fold-change of VEGF.

Secondly, tumor receptor heterogeneity influences the effect of combination therapy on individual angiogenic pathways. The model predicts that the fold-change of the VEGFR2- and neuropilin-1-containing pro-angiogenic complexes in tumor could be higher than one for some sampled cases (Fig. S5, ESI†). This indicates that for some tumor receptor profiles, the pro-angiogenic signaling pathway of VEGFR2 or NRP-1 will be enhanced after treatment instead of inhibited. Similar results are seen when simulating bevacizumab treatment (data not shown). Interestingly, these unintended outcomes following bevacizumab treatment predicted by the model have been observed clinically: a substantial increase in phosphorylated VEGFR2 following bevacizumab treatment⁵² and variable changes in phosphorylated VEGFR2 in the clinical trial of the combination anti-angiogenic therapy.¹⁷

Tumor receptor heterogeneity affects the drug combination effect

We used the model to compare the response to single- and dual-agent therapy and find that the response to combination therapy is related to the response to single-agent treatment (Fig. 4A). The color change from red to yellow along the diagonal direction (from upper right to lower left) shows that the response to combination treatment (Fc) improves when the response to either single-agent treatment is higher. This indicates that a better response to single-agent treatment may correspond to a better response to combination therapy. We then examined the superiority of the drug combination compared to the single agents. The color of the points in Fig. 4B indicates the difference between the response to combination therapy and the best response amongst the two single-agent treatments. Although this difference changes depending on the tumor receptor profile, the difference is higher than zero for all sampled cases. This indicates that the response to combination therapy is better than single agent treatment. Thus, when considering the fold-change of tumor angiogenic ratio (F), the combination therapy of bevacizumab and ABT-510 is superior to single-agent therapy for all of the sampled tumor receptor profiles.

Fig. 4 The impact of tumor receptor heterogeneity on the combination effect. (A) The comparison of the response to combination therapy (Fc) and responses to single agent treatments (Fa: response to ABT-510 alone; Fb: response to bevacizumab alone). (B) The comparison of the response to combination therapy (Fc) and the best response between the individual treatments (the smaller of Fa and Fb). (C) The calculated combination index (CI) for four representative tumor receptor profiles.

We explicitly calculated the combination effect to determine if the action of the two drugs is synergistic or additive. We selected four sets of tumor receptor profiles (Fig. 4A, triangles) and calculated the combination index (CI) of the Loewe additivity (see Methods). These are representative cases corresponding to patients that respond well to ABT-510 alone (“1”), bevacizumab alone (“2”), or combination therapy (“3”), and one patient that does not respond well to either drug or the combination (“4”). The calculated CI values are only slightly different for these four tumor profiles (Fig. 4C), where the combination effect is nearly additive for all four patients (CI ranges from 0.88 to 0.96). Thus, the bevacizumab and ABT-510 have an additive effect in shifting the tumor angiogenic balance for the four selected patients.

The therapeutic response depends on tumor cell VEGR1, CD36 and CD47 expression

Given the observed impact of tumor receptor heterogeneity in the predictions from the mechanistic model, we explored if we can use tissue biomarkers (tumor receptor levels) to determine the response to combination anti-angiogenic therapy a priori. We used partial least squares regression (PLSR) to characterize the association between tumor receptor levels and the response to combination therapy (the fold-change in the angiogenic ratio in the tumor). Specifically, we established a predictive PLSR model that estimates the response to combination therapy using the tumor angiogenic receptor profile (Fig. S6a, ESI†). We used 18 variables as the inputs in PLSR modeling: the 16 VEGF and TSP1 receptor numbers on diseased endothelial cells and tumor cells that were varied and two nonlinear terms comprised of ratios of receptor numbers. The fold-changes estimated by the PLSR model closely match the actual values given by the mechanistic model (the points lie along the 45°-angle line). The high Q²Y value of 0.94, indicates that in ten-fold cross validation, the PLSR model predicts over 94% of the total variation in the fold-change of tumor angiogenic ratio. The R²Y value of 0.93 shows the most of the variation in the response is also captured by the PLSR model, while the R²X value indicates that only a small portion of the receptor information (15%) is useful for predicting the therapeutic response.

The variable importance of projection (VIP) score was calculated to find which variables, across the entire PLSR model, most significantly contribute to determining the response to combination therapy.⁵¹ VEGFR1, CD36 and CD47 levels on tumor cells were found to be important for explaining the response, as they have VIP scores higher than one (Fig. 5A). The two nonlinear terms of tumor cell receptor numbers are also identified as important variables. We found that all of the receptors on tumor endothelial cells have VIP score lower than 0.5, indicating that tumor endothelial cell receptor expression was not as influential as receptor expression on tumor cells in predicting the response to combination therapy. Plotting the values of the receptors with VIP scores greater than one versus the response to combination therapy shows that lower VEGFR1 expression, lower CD47 expression and higher CD36 expression on tumor cells correlate with higher response (Fig. 5B–D). This indicates the patients with low tumor cell VEGFR1 or CD47 expression or high tumor cell CD36 expression are likely to show a stronger response to combination therapy than others. The associations between response and other receptors are not pronounced (Fig. S7, ESI†). The therapeutic response is also related to the two nonlinear terms in a monotonic fashion (Fig. 5E and F). Interestingly, although the response to treatment exhibits a trend that correlates with certain individual receptors, the data highly deviates from the LOWESS smoothing of the data (Fig. 5, orange lines). This indicates that none of the 16 receptor variables can accurately predict the response alone, and the variability in the expression of a single receptor cannot account for all of the variability in the response to combination anti-angiogenic therapy. Hence, the ratio of receptor numbers is a valuable quantity in predicting the response to combination therapy.

Fig. 5 The association between tumor angiogenic receptors and the response to combination anti-angiogenic therapy. (A) The calculated variable importance of projection (VIP) score estimated by the PLSR model. VEGFR1, CD36 and CD47 on tumor cells have VIP scores higher than one (orange bars), indicating that they each significantly contribute to the response to treatment, the fold-change in the angiogenic ratio in the tumor. (B–F) The associations between normalized receptor level and the response to combination therapy. The orange line is the locally weighted scatterplot smoothing (LOWESS).

Discussion

We present a systems biology model of the angiogenic balance in breast cancer patients. We significantly expand upon our previous work,²² extending the single-compartment model to three compartments. This whole-body model allows us to account for the PK characteristics and distribution of anti-angiogenic drugs, which was not possible with our previous single-compartment model. This new model enables us to investigate the combination anti-angiogenic therapy targeting the VEGF and TSP1 axes. Additionally, the model parameters represent physiological and biological characteristics of cancer patients, allowing us to answer clinically relevant questions. Specifically, we investigated the impact of tumor receptor variability, one key aspect of inter-patient heterogeneity.

Overall, in the present study, we build our model based on biological knowledge and calibrate the model with both in vivo and in vitro data. We use the model to simulate anti-angiogenic therapy, generating new predictions and quantifying the impact of tumor receptor heterogeneity on the response to anti-angiogenic tumor therapy. Several model predictions are supported by published clinical observations, as discussed below.

This work uses a molecular-based PK/PD model to bridge the heterogeneous response of anti-angiogenic therapy with tumor receptor heterogeneity. Tumor angiogenic receptor heterogeneity has been observed in various settings. Qualitative information on the differential expression of angiogenic receptors from different tumor samples is reported in the Human Protein Atlas Database.^26,27 Quantitative measurements also show the cell-to-cell variations in the density of VEGF receptors across different cell lines.^24,25 This heterogeneity is thought to influence prognosis and response to treatment, which has been observed in various experimental settings. The VEGFR2 heterogeneity was proved to influence the response to anti-angiogenic cyclophosphamide treatment in a pre-clinical model.²⁸ The inter-patient VEGFR1 and neuropilin heterogeneity has been associated with the response to bevacizumab treatment.^29,30 Our modeling results provide a possible quantitative characterization of the mechanisms leading to the clinically observed heterogeneous responses to combination anti-angiogenic therapy,^3,8,31 complementing the pre-clinical and clinical studies. The model predicts that after one cycle of treatment, the fold-change of the angiogenic ratio in tumor tissue varies widely depending on the tumor receptor levels. This indicates that tumor receptor heterogeneity could contribute to the heterogeneous response and should be considered when developing more effective and personalized anti-angiogenic therapy.

The characterization of the response in our study (the fold-change of angiogenic ratio) can be studied and measured in the experimental setting. The binding of angiogenic factors to their receptors initiates a cascade of intracellular reactions, including activation or phosphorylation of the receptors. The phosphorylated receptors can be measured (via flow cytometry or immunohistochemistry analyses of a tumor sample) and used to calculate the ratio of the pro- and anti-angiogenic signals. Thus, the fold-change of angiogenic ratio reflects how the anti-angiogenic therapy changes the amount of phosphorylated pro-angiogenic receptors relative to the activated anti-angiogenic receptors, which together mediate tumor angiogenesis. Multiple experimental studies qualitatively explore the activated (phosphorylated) levels of the angiogenic receptors in the context of the TSP1–VEGF axis;^47,53–56 however, quantitative data needed to validate our model are still missing. Therefore, a direct comparison with the clinical measurements is not possible. However, our modeling results are meaningful, as they identify important potential factors to measure (the tumor angiogenic receptors level before treatment), and can help design pre-clinical and clinical studies.

Our work helps identify biomarkers that predict response to anti-angiogenic treatment. The heterogeneous response to anti-angiogenic therapy is a major drawback for optimizing the combination therapy of bevacizumab and ABT-510, and generally, when considering all anti-angiogenic agents.^3,31 However, despite substantial research,⁵⁷ there is currently no reliable biomarker for identifying patients for which anti-angiogenic therapy will be effective. The PLSR analysis in our study showed that the levels of CD36, CD47, and VEGFR1 on tumor cells are potential biomarkers to select the patients who will likely respond to combination treatment. These receptors are shown to have a pronounced impact on the response: low expression of VEGFR1 and CD47 and high expression of CD36 all associate with high response to combination treatment. Interestingly, the relationship between VEGFR1 and the effect of bevacizumab has also been observed in clinical trial. Fountzilas et al. report that intra-tumoral VEGFR1 overexpression was an indicator for poor survival in metastatic breast cancer patients receiving bevacizumab and chemotherapy,⁵⁸ which supports our finding and provides a qualitative foundation for our study. In addition to identifying these biomarkers, our study emphasizes the importance of the integration of information from multiple receptors to accurately predict the response. We showed that the expression of a single receptor is not able to fully account for the heterogeneous response to therapy. Instead, a comprehensive analysis that incorporates multiple biomarkers is more appropriate. In the experimental and clinical practice of biomarker discovery, it is commonly accepted that a single biomarker is not able to predict anti-angiogenic treatment efficacy.⁵⁷ Our findings support this stance and demonstrate that a computational approach could help identify biomarkers by providing systematic understanding of the response to treatment. Specifically, the PLSR analysis indicated that the CD36/(CD47 + LRP1) ratio has the highest VIP score and is most important in predicting the fold-change in the angiogenic ratio in the tumor, indicating the importance of multiple receptors as a potential tissue biomarker.

Given the molecular detail of the model, we gain mechanistic insight into the model predictions. For example, we can understand how the individual tumor receptor populations change in response to anti-angiogenic therapy, insight that explains why VEGFR1, CD36, and CD47 are potential biomarkers. VEGFR1 has been supported as a possible in situ predictive biomarker for anti-VEGF agents in both preclinical and clinical settings.^57–59 Our model predicts that the pro-angiogenic complexes involving VEGFR1 play an important role in the response to bevacizumab treatment: complexes involving VEGFR1 significantly decrease after the administration of bevacizumab while the levels of other pro-angiogenic complexes are relatively stable before and after bevacizumab treatment. The high expression of VEGFR1 can efficiently compete for VEGF, thereby limiting the benefits of VEGF neutralization by bevacizumab. On the other side, high expression of CD36 provides more available receptors for ABT-510 to bind to, which enables the formation of more anti-angiogenic complexes to shift the angiogenic balance. Since CD47 has a high affinity for TSP1, increasing the expression of CD47 will significantly elevate the baseline level of anti-angiogenic complexes. This means that the tumor has a higher baseline anti-angiogenic signal and consequently shows an attenuated response to the formation of more anti-angiogenic complexes by ABT-510. The insights generated by our model can guide future experimental or clinical studies of predictive biomarkers for anti-angiogenic treatment.

We do not include clinical response in this study due to lack of data needed to connect PK/PD data with clinical outcome. Thus, we currently cannot distinguish “complete response” and “partial response” in our model prediction at this point. However, we are able to get insights into outcomes from clinical trials involving anti-angiogenic agents from the predictions of our model at a qualitative level. ABT-510 failed in Phase II clinical trials for advanced renal cell carcinoma and metastatic melanoma.^13,14 According to our model simulations, low CD36 expression on tumor cell could lead to the low observed response to ABT-510. Interestingly, the Human Protein Atlas shows CD36 is not detected in human renal cancer samples (in 10 out of 10 samples) and melanoma samples (in 11 out of 11 samples).²⁷ The heterogeneous expression of CD36 might be a possible reason for the failure of ABT-510. In the combination of bevacizumab and ABT-510 for advanced solid tumors, only 32% percent of the patients have prolonged stable disease for more than 6 months,¹⁷ which shows no evidence of being better than bevacizumab alone. Our model predicts that the difference between combination therapy and the best response of individual bevacizumab or ABT-510 treatment will become zero when the response to ABT-510 is very low. This means that the combination will not be superior to bevacizumab if ABT-510 alone is highly ineffective, which provides a possible explanation for the clinical observations.

Although we specifically simulate the combination of bevacizumab and ABT-510 in this work, our model can be adapted to study other anti-angiogenic agents that target VEGF and TSP1 signaling and can help explore different anti-angiogenic strategies. Given the many different anti-angiogenic agents^3,12 and possible drug regimens, our realistic model could be a valuable help in the search for an optimal combination anti-angiogenic strategy. Additionally, this model has practical interest, which can help circumvent difficulty in experimental studies. Traditional experiment-based study of comparing the effects of combination therapy to individual agents is limited by the difficulty of accumulating enough data to make meaningful comparisons. In contrast, our model can easily simulate multiple anti-angiogenic treatment strategies and quantitatively evaluate the added benefit of combination therapy as compared to single-agent therapy. In this study, we show that the combination of bevacizumab and ABT-510 leads to a stronger effect than administering either agent alone, in all of the sampled cases. However, the combination index revealed that the enhanced effect does not mean a strong synergistic effect. Our model would be a help to understand the synergism between two drugs and to design the combination therapy that achieves high efficacy with low drug dosage.

We acknowledge some limitations of the model that can be addressed in future studies. In this study, we use the angiogenic ratio to represent the angiogenic state, assuming that each angiogenic complex equally contributes to angiogenesis. However, it is known that the angiogenic receptors have unique functions and influence angiogenesis in different ways. For example, VEGFR1 primarily modulates blood vessel angiogenesis by ligand-trapping and receptor dimerization, while VEGFR2 is the predominant receptor that promotes pro-angiogenic VEGF signaling pathways.^60,61 On the anti-angiogenic side, CD36 and CD47 are reported to inhibit angiogenesis by antagonizing survival pathways and activating apoptotic pathways.⁶² This knowledge of the relevant intracellular signaling pathways can be combined with our model to study specific functional responses to anti-angiogenic therapy. To date, several multiscale modeling studies have linked PK/PD models with the downstream signaling and cell-surface reactions.^63,64 We can implement a similar approach, combining our model with downstream signaling models^65–67 to enable broader application of our model. In addition, anti-angiogenic therapy is usually used in combination with chemotherapy in the clinic. The variable response of cell lines to the chemotherapeutic agent doxorubicin has been recently modeled computationally.⁶⁸ It is possible to combine these models and consider both anti-angiogenic and chemotherapeutic agents. This might help improve our current work and generate new insights. Secondly, it is worth noting that there are alternate ways of characterizing the magnitude of the response to anti-angiogenic therapy. We focus on the angiogenic ratio, where the balance of pro- and anti-angiogenic factors has been experimentally shown to be a more accurate approach to study angiogenesis than analyzing the level of angiogenic factors individually.⁶⁹ We further calculate the fold-change of the angiogenic ratio to characterize the response to anti-angiogenic treatment, while the area under the angiogenic ratio curve for the tumor compartment could be an alternative way to quantify the response to treatment. In that case, the response to ABT-510 is more substantial due to the strong shifting of the angiogenic ratio in the short time after each bolus (data not shown). However, the conclusions of our study remain unchanged. Another metric to consider is whether anti-angiogenic treatment shifts the angiogenic balance in tumor to the level observed in healthy tissue, related to vessel normalization.⁷⁰ Since a universal definition of the response to anti-angiogenic therapy is still missing, the fold-change of the angiogenic balance of pro- and anti-angiogenic receptor complexes remains an important indicator, which this research is well suited to address. The primary focus for this study is the tumor receptor heterogeneity, and we do not investigate other possible inter-patient variations in this work, such as variation in the drugs’ PK properties, different tumor sizes, or variable secretion of the angiogenic factors variation. Currently, these variations are poorly quantified and studied in pre-clinical and clinical experiments. Using our model as a tool to understand the impact of these other aspects of heterogeneity might be a great interest for future research.

Concluding thoughts

Our model illustrates the effect of combination therapy of bevacizumab and ABT-510 on changing the balance between two opposing angiogenic signals. The model provides a quantitative description of the impact of tumor receptor heterogeneity on the response to combination anti-angiogenic therapy and aids in the discovery of predictive biomarkers. We expect that the insights generated by our model predictions will shed light on previously obscure clinical observations and that the model will be used to facilitate the optimization of new clinical trials.

Funding

The authors acknowledge the support of the US National Science Foundation (CAREER Award 1552065 to S. D. F.), the American Cancer Society (130432-RSG-17-133-01-CSM to S. D. F.) and the USC Provost's PhD Fellowship (D. L.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The authors thank members of the Finley research group for critical comments and suggestions.

References

1. D. Hanahan and R. A. Weinberg, Hallmarks of cancer: the next generation, Cell, 2011, 144(5), 646–674,  DOI:10.1016/j.cell.2011.02.013.
2. N. Nishida, H. Yano, T. Nishida, T. Kamura and M. Kojiro, Angiogenesis in cancer, Vasc. Health Risk Manage., 2006, 2(3), 213–219.
3. N. S. Vasudev and A. R. Reynolds, Anti-angiogenic therapy for cancer: current progress, unresolved questions and future directions, Angiogenesis, 2014, 17(3), 471–494.
4. D. Hanahan and J. Folkman, Patterns and emerging mechanisms of the angiogenic switch during tumorigenesis, Cell, 1996, 86(3), 353–364.
5. T. Sakurai and M. Kudo, Signaling pathways governing tumor angiogenesis, Oncology, 2011, 81, 24–39.
6. Z. Huang and S. D. Bao, Roles of main pro- and anti-angiogenic factors in tumor angiogenesis, World J. Gastroenterol., 2004, 10, 463–470.
7. K. L. Meadows and H. I. Hurwitz, Anti-VEGF therapies in the clinic, Cold Spring Harbor Perspect. Med., 2012, 2(10), a006577.
8. M. Wehland, J. Bauer, M. Infanger and D. Grimm, Target-based anti-angiogenic therapy in breast cancer, Curr. Pharm. Des., 2012, 18(27), 4244–4257.
9. B. Ren, K. O. Yee, J. Lawler and R. Khosravi-Far, Regulation of tumor angiogenesis by thrombospondin-1, Biochim. Biophys. Acta, 2006, 1765(2), 178–188.
10. J. Lawler, Thrombospondin-1 as an endogenous inhibitor of angiogenesis and tumor growth, J. Cell. Mol. Med., 2002, 6(1), 1–12 . Available from: http://www.ncbi.nlm.nih.gov/pubmed/12003665.
11. L. C. Armstrong and P. Bornstein, Thrombospondins 1 and 2 function as inhibitors of angiogenesis, Matrix Biol., 2003, 63–71.
12. A. Jeanne, C. Schneider, L. Martiny and S. Dedieu, Original insights on thrombospondin-1-related antireceptor strategies in cancer, Front. Pharmacol., 2015, 6, 252.
13. S. N. Markovic, V. J. Suman, R. A. Rao, J. N. Ingle, J. S. Kaur and L. A. Erickson, et al., A phase II study of ABT-510 (thrombospondin-1 analog) for the treatment of metastatic melanoma, Am. J. Clin. Oncol., 2007, 30(3), 303–309 . Available from: http://ezproxy.library.usyd.edu.au/login?url=http://ovidsp.ovid.com/ovidweb.cgi?T=JS&CSC=Y&NEWS=N&PAGE=fulltext&D=med5&AN=17551310%5Cnhttp://dd8gh5yx7k.search.serialssolutions.com/?sid=OVID:medline&id=pmid:17551310&id=doi:&issn=0277-3732&isbn=&volume=30&i.
14. S. Ebbinghaus, M. Hussain, N. Tannir, M. Gordon, A. A. Desai and R. A. Knight, et al., Phase 2 study of ABT-510 in patients with previously untreated advanced renal cell carcinoma, Clin. Cancer Res., 2007, 13(22 Pt 1), 6689–6695 . Available from: http://www.ncbi.nlm.nih.gov/pubmed/18006769.
15. L. Yadav, N. Puri, V. Rastogi, P. Satpute and V. Sharma, Tumour angiogenesis and angiogenic inhibitors: a review, J. Clin. Diagn. Res., 2015, 9, XE01–XE05.
16. G. Bergers and L. E. Benjamin, Angiogenesis: Tumorigenesis and the angiogenic switch, Nat. Rev. Cancer, 2003, 3(6), 401–410 . Available from: http://www.nature.com/doifinder/10.1038/nrc1093.
17. H. E. Uronis, S. M. Cushman, J. C. Bendell, G. C. Blobe, M. A. Morse and A. B. Nixon, et al., A phase I study of ABT-510 plus bevacizumab in advanced solid tumors, Cancer Med., 2013, 2(3), 316–324 . Available from: http://www.ncbi.nlm.nih.gov/pubmed/23930208%5Cnhttp://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC3699843.
18. S. D. Finley, L.-H. Chu and A. S. Popel, Computational systems biology approaches to anti-angiogenic cancer therapeutics, Drug Discovery Today, 2015, 20(2), 187–197 . Available from: http://linkinghub.elsevier.com/retrieve/pii/S135964461400395X.
19. H. Jain and T. Jackson, Mathematical Modeling of Cellular Cross-Talk Between Endothelial and Tumor Cells Highlights Counterintuitive Effects of VEGF-Targeted Therapies, Bull. Math. Biol., 2017, 1–46.
20. L. G. Hutchinson, E. A. Gaffney, P. K. Maini, J. Wagg, A. Phipps and H. M. Byrne, Vascular phenotype identification and anti-angiogenic treatment recommendation: A pseudo-multiscale mathematical model of angiogenesis, J. Theor. Biol., 2016, 398, 162–180,  DOI:10.1016/j.jtbi.2016.03.002.
21. M. O. Stefanini, F. T. H. Wu, F. Mac Gabhann and A. S. Popel, Increase of plasma VEGF after intravenous administration of bevacizumab is predicted by a pharmacokinetic model, Cancer Res., 2010, 70(23), 9886–9894.
22. J. A. Rohrs, C. D. Sulistio and S. D. Finley, Predictive model of thrombospondin-1 and vascular endothelial growth factor in breast tumor tissue, NPJ Syst. Biol. Appl., 2016, 2, 16030 . Available from: http://dx.doi.org/10.1038/npjsba.2016.30%5Cnhttp://10.1038/npjsba.2016.30%5Cnhttp://www.nature.com/articles/npjsba201630#supplementary-information.
23. S. D. Finley, P. Angelikopoulos, P. Koumoutsakos and A. S. Popel, Pharmacokinetics of Anti-VEGF Agent aflibercept in cancer predicted by data-driven, molecular-detailed model, CPT: Pharmacometrics Syst. Pharmacol., 2015, 4(11), 641–649.
24. P. I. Imoukhuede and A. S. Popel, Quantitative fluorescent profiling of VEGFRs reveals tumor cell and endothelial cell heterogeneity in breast cancer xenografts, Cancer Med., 2014, 3(2), 225–244.
25. P. I. I. Imoukhuede and A. S. Popel, Quantification and cell-to-cell variation of vascular endothelial growth factor receptors, Exp. Cell Res., 2011, 317(7), 955–965 . Available from: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3073416/%5Cnhttp://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3073416&tool=pmcentrez&rendertype=abstract.
26. M. Uhlen, L. Fagerberg, B. M. Hallstrom, C. Lindskog, P. Oksvold and A. Mardinoglu, et al., Tissue-based map of the human proteome, Science, 2015, 347(6220), 1260419 . Available from: http://www.sciencemag.org/cgi/doi/10.1126/science.1260419.
27. M. Uhlen, A Human Protein Atlas for Normal and Cancer Tissues Based on Antibody Proteomics, Mol. Cell. Proteomics, 2005, 4(12), 1920–1932 . Available from: http://www.mcponline.org/cgi/doi/10.1074/mcp.M500279-MCP200.
28. S. G. Patten, U. Adamcic, K. Lacombe, K. Minhas, K. Skowronski and B. L. Coomber, VEGFR2 heterogeneity and response to anti-angiogenic low dose metronomic cyclophosphamide treatment, BMC Cancer, 2010, 10, 683 . Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3009683&tool=pmcentrez&rendertype=abstract.
29. A. M. Jubb, K. D. Miller, H. S. Rugo, A. L. Harris, D. Chen and J. D. Reimann, et al., Impact of exploratory biomarkers on the treatment effect of bevacizumab in metastatic breast cancer, Clin. Cancer Res., 2011, 17(2), 372–381.
30. E. Van Cutsem, S. de Haas, Y.-K. Kang, A. Ohtsu, N. C. Tebbutt and J. Ming Xu, et al., Bevacizumab in combination with chemotherapy as first-line therapy in advanced gastric cancer: a biomarker evaluation from the AVAGAST randomized phase III trial, J. Clin. Oncol., 2012, 30(17), 2119–2127 . Available from: http://www.ncbi.nlm.nih.gov/pubmed/22565005.
31. L. Moserle, G. Jiménez-Valerio and O. Casanovas, Antiangiogenic therapies: going beyond their limits, Cancer Discovery, 2014, 4, 31–41.
32. M. O. Stefanini, F. T. H. Wu, F. Mac Gabhann and A. S. Popel, A compartment model of VEGF distribution in blood, healthy and diseased tissues, BMC Syst. Biol., 2008, 2, 77 . Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2562372&tool=pmcentrez&rendertype=abstract.
33. F. Mac Gabhann, J. W. Ji and A. S. Popel, Multi-scale computational models of pro-angiogenic treatments in peripheral arterial disease, Ann. Biomed. Eng., 2007, 982–994.
34. S. D. Finley and A. S. Popel, Effect of Tumor Microenvironment on Tumor VEGF During Anti-VEGF Treatment: Systems Biology Predictions, J. Natl. Cancer Inst., 2013, 105(11), 804–811.
35. J. R. Faeder, M. L. Blinov and W. S. Hlavacek, Rule-based modeling of biochemical systems with BioNetGen, Methods Mol. Biol., 2009, 500(1), 113–167.
36. F. Mac Gabhann and A. S. Popel, Targeting neuropilin-1 to inhibit VEGF signaling in cancer: comparison of therapeutic approaches, PLoS Comput. Biol., 2006, 2(12), 1649–1662.
37. M. O. Stefanini, F. T. H. Wu, F. Gabhann Mac and A. S. Popel, The presence of VEGF receptors on the luminal surface of endothelial cells affects VEGF distribution and VEGF signaling, PLoS Comput. Biol., 2009, 5(12), e1000622.
38. K. Tan, M. Duquette, J. H. Liu, Y. Dong, R. Zhang and A. Joachimiak, et al., Crystal structure of the TSP-1 type 1 repeats: a novel layered fold and its biological implication, J. Cell Biol., 2002, 159(2), 373–382.
39. T. Klein and R. Bischoff, Physiology and pathophysiology of matrix metalloproteases, Amino Acids, 2011, 41, 271–290.
40. M. S. Gordon, K. Margolin, M. Talpaz, G. W. Sledge, E. Holmgren and R. Benjamin, et al., Phase I safety and pharmacokinetic study of recombinant human anti-vascular endothelial growth factor in patients with advanced cancer, J. Clin. Oncol., 2001, 19(3), 843–850.
41. J. D. Jensen, L. W. Jensen and J. K. Madsen, The pharmacokinetics of recombinant human erythropoietin after subcutaneous injection at different sites, Eur. J. Clin. Pharmacol., 1994, 46(4), 333–337.
42. R. Hoekstra, F. Y. de Vos, F. A. Eskens, J. A. Gietema, A. van der Gaast and H. J. M. Groen, et al., Phase I safety, pharmacokinetic, and pharmacodynamic study of the thrombospondin-1-mimetic angiogenesis inhibitor ABT-510 in patients with advanced cancer, J. Clin. Oncol., 2005, 23(22), 5188–5197.
43. S. Marino, I. B. Hogue, C. J. Ray and D. E. Kirschner, A methodology for performing global uncertainty and sensitivity analysis in systems biology, J. Theor. Biol., 2008, 254, 178–196.
44. L. Kagan, Pharmacokinetic modeling of the subcutaneous absorption of therapeutic proteins, Drug Metab. Dispos., 2014, 42(11), 1890–1905.
45. J. C. Anderson, J. R. Grammer, W. Wang, L. B. Nabors, J. Henkin and J. E. Stewart, et al., ABT-510, a modified type 1 repeat peptide of thrombospondin, inhibits malignant glioma growth in vivo by inhibiting angiogenesis, Cancer Biol. Ther., 2007, 6(3), 454–462.
46. S. Kaur, G. Martin-Manso, M. L. Pendrak, S. H. Garfield, J. S. Isenberg and D. D. Roberts, Thrombospondin-1 inhibits VEGF receptor-2 signaling by disrupting its association with CD47, J. Biol. Chem., 2010, 285(50), 38923–38932.
47. L. Y. Chu, D. P. Ramakrishnan and R. L. Silverstein, Thrombospondin-1 modulates VEGF signaling via CD36 by recruiting SHP-1 to VEGFR2 complex in microvascular endothelial cells, Blood, 2013, 122(10), 1822–1832.
48. J. Foucquier and M. Guedj, Analysis of drug combinations: current methodological landscape, Pharmacol. Res. Perspect., 2015, 3(3), e00149.
49. T. C. Chou, Drug combination studies and their synergy quantification using the chou-talalay method, Cancer Res., 2010, 70(2), 440–446.
50. P. Geladi and B. R. Kowalski, Partial least-squares regression: a tutorial, Anal. Chim. Acta, 1986, 185, 1–17.
51. P. K. Kreeger, Using partial least squares regression to analyze cellular response data, Sci. Signaling, 2013, 6(271), tr7 . Available from: http://www.ncbi.nlm.nih.gov/pubmed/23592846.
52. S. B. Wedam, J. A. Low, S. X. Yang, C. K. Chow, P. Choyke and D. Danforth, et al., Antiangiogenic and antitumor effects of bevacizumab in patients with inflammatory and locally advanced breast cancer, J. Clin. Oncol., 2006, 24(5), 769–777.
53. K. Osz, M. Ross and J. Petrik, The thrombospondin-1 receptor CD36 is an important mediator of ovarian angiogenesis and folliculogenesis, Reprod. Biol. Endocrinol., 2014, 12(1), 21.
54. J. Sun, B. D. Hopkins, K. Tsujikawa, C. Perruzzi, I. Adini and R. Swerlick, et al., Thrombospondin-1 modulates VEGF-A-mediated Akt signaling and capillary survival in the developing retina, Am. J. Physiol. Heart Circ. Physiol., 2009, 296(5), H1344–H1351.
55. S. Russell, M. Duquette, J. Liu, R. Drapkin, J. Lawler and J. Petrik, Combined therapy with thrombospondin-1 type I repeats (3TSR) and chemotherapy induces regression and significantly improves survival in a preclinical model of advanced stage epithelial ovarian cancer, FASEB J., 2015, 29(2), 576–588.
56. S. Kazerounian, M. Duquette, M. A. Reyes, J. T. Lawler, K. Song and C. Perruzzi, et al., Priming of the vascular endothelial growth factor signaling pathway by thrombospondin-1, CD36, and spleen tyrosine kinase, Blood, 2011, 117(17), 4658–4666.
57. M. Wehland, J. Bauer, N. E. Magnusson, M. Infanger and D. Grimm, Biomarkers for Anti-Angiogenic Therapy in Cancer, , 2013, 9338–9364.
58. G. Fountzilas, H. P. Kourea, M. Bobos, D. Televantou, V. Kotoula and C. Papadimitriou, et al., Paclitaxel and Bevacizumab as First Line Combined Treatment in Patients with Metastatic Breast Cancer: The Hellenic Cooperative Oncology Group Experience with Biological Marker Evaluation, Anticancer Res., 2011, 3018, 3007–3018.
59. D. Lambrechts, H. J. Lenz, S. De Haas, P. Carmeliet and S. J. Scherer, Markers of response for the antiangiogenic agent bevacizumab, J. Clin. Oncol., 2013, 31, 1219–1230.
60. M. Shibuya, Vascular Endothelial Growth Factor (VEGF) and Its Receptor (VEGFR) Signaling in Angiogenesis: A Crucial Target for Anti- and Pro-Angiogenic Therapies, Genes Cancer, 2011, 2(12), 1097–1105 . Available from: http://gan.sagepub.com/lookup/doi/10.1177/1947601911423031.
61. N. Rahimi, VEGFR-1 and VEGFR-2: two non-identical twins with a unique physiognomy, Front. Biosci., 2006, 11, 818–829 . Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1360224&tool=pmcentrez&rendertype=abstract.
62. P. R. Lawler and J. Lawler, Molecular basis for the regulation of angiogenesis by thrombospondin-1 and -2, Cold Spring Harbor Perspect. Med., 2012, 2(5), a006627.
63. C. Nanavati, D. Ruszaj and D. E. Mager, Cell signaling model connects vorinostat pharmacokinetics and tumor growth response in multiple myeloma xenografts, CPT: Pharmacometrics Syst. Pharmacol., 2017, 6(11), 756–764.
64. G. Dwivedi, L. Fitz, M. Hegen, S. W. Martin, J. Harrold and A. Heatherington, et al., A Multiscale Model of Interleukin-6–Mediated Immune Regulation in Crohn's Disease and Its Application in Drug Discovery and Development, CPT: Pharmacometrics Syst. Pharmacol., 2014, 3(1), e89 . Available from: http://doi.wiley.com/10.1038/psp.2013.64.
65. Q. Wu and S. D. Finley, Predictive model identifies strategies to enhance TSP1-mediated apoptosis signaling, Cell Commun. Signaling, 2017, 15(1), 53.
66. W. H. Tan, A. S. Popel and F. Mac Gabhann, Computational model of VEGFR2 pathway to ERK activation and modulation through receptor trafficking, Cell. Signalling, 2013, 25(12), 2496–2510.
67. J. Kanodia, D. Chai, J. Vollmer, J. Kim, A. Raue and G. Finn, et al., Deciphering the mechanism behind Fibroblast Growth Factor (FGF) induced biphasic signal-response profiles, Cell Commun. Signaling, 2014, 12, 34.
68. M. T. McKenna, J. A. Weis, S. L. Barnes, D. R. Tyson, M. I. Miga and V. Quaranta, et al., A Predictive Mathematical Modeling Approach for the Study of Doxorubicin Treatment in Triple Negative Breast Cancer, Sci. Rep., 2017, 7(1), 5725.
69. E. Roudier, C. Gineste, A. Wazna, K. Dehghan, D. Desplanches and O. Birot, Angio-adaptation in unloaded skeletal muscle: new insights into an early and muscle type-specific dynamic process, J. Physiol., 2010, 588(22), 4579–4591,  DOI:10.1113/jphysiol.2010.193243.
70. R. K. Jain, Normalizing tumor microenvironment to treat cancer: bench to bedside to biomarkers, J. Clin. Oncol., 2013, 2205–2218.
71. G. J. Byrne, K. E. Hayden, G. McDowell, H. Lang, C. C. Kirwan and L. Tetlow, et al., Angiogenic characteristics of circulating and tumoural thrombospondin-1 in breast cancer, Int. J. Oncol., 2007, 31(5), 1127–1132.
72. E. A. Baker, D. J. Leaper, J. P. Hayter and A. J. Dickenson, The matrix metalloproteinase system in oral squamous cell carcinoma, Br. J. Oral Maxillofac. Surg., 2006, 44(6), 482–486 . Available from: http://www.sciencedirect.com/science/article/pii/S0266435605003220.
73. C. Kut, F. Mac Gabhann and A. S. Popel, Where is VEGF in the body? A meta-analysis of VEGF distribution in cancer, Br. J. Cancer, 2007, 97(7), 978–985 . Available from: http://www.nature.com/doifinder/10.1038/sj.bjc.6603923.
74. F. Vasaturo, F. Solai, C. Malacrino, T. Nardo, B. Vincenzi and M. Modesti, et al., Plasma levels of matrix metalloproteinases 2 and 9 correlate with histological grade in breast cancer patients, Oncol. Lett., 2012, 5(1), 316–320.
75. B. Hoier, M. Passos, J. Bangsbo and Y. Hellsten, Intense intermittent exercise provides weak stimulus for vascular endothelial growth factor secretion and capillary growth in skeletal muscle, Exp. Physiol., 2013, 98(2), 585–597,  DOI:10.1113/expphysiol.2012.067967.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8ib00019k

---