Mechanistic Modeling Quantifies The Influence Of Tumor Growth Kinetics On The Response To Anti-Angiogenic Treatment
===================================================================================================================

* Thomas D. Gaddy
* Stacey D. Finley

## Abstract

Tumors exploit angiogenesis, the formation of new blood vessels from pre-existing vasculature, in order to obtain nutrients required for continued growth and proliferation. Targeting factors that regulate angiogenesis, including the potent promoter vascular endothelial growth factor (VEGF), is therefore an attractive strategy for inhibiting tumor growth. Systems biology modeling enables us to identify tumor-specific properties that influence the response to those anti-angiogenic strategies. Here, we build on our previous systems biology model of VEGF transport and kinetics in tumor-bearing mice to include a tumor compartment whose volume depends on the “angiogenic signal” produced when VEGF binds to its receptors on tumor endothelial cells. We trained and validated the model using *in vivo* measurements of xenograft tumor volume to produce a model that accurately predicts the tumor’s response to anti-angiogenic treatment. We applied the model to investigate how tumor growth kinetics influence the response to anti-angiogenic treatment targeting VEGF. Based on multivariate regression analysis, we found that certain intrinsic kinetic parameters that characterize the growth of tumors could successfully predict response to anti-VEGF treatment. This model is a useful tool for predicting which tumors will respond to anti-VEGF treatment, complementing pre-clinical *in vivo* studies.

## Introduction

Angiogenesis is the formation of new blood vessels from pre-existing vasculature and is important in both physiological and pathological conditions. Numerous promoters and inhibitors regulate angiogenesis. One key promoter of angiogenesis is the vascular endothelial growth factor-A (VEGF-A), which has been extensively studied and is a member of a family of pro-angiogenic factors that includes five ligands: VEGF-A, VEGF-B, VEGF-C, VEGF-D, and placental growth factor (PlGF). VEGF-A (or simply, VEGF) promotes angiogenesis by binding to its receptors VEGFR1 and VEGFR2 and recruiting co-receptors called neuropilins (NRP1 and NRP2). The VEGF receptors and co-receptors are expressed on many different cell types, including endothelial cells (ECs), cancer cells, neurons, and muscle fibers1. Together, VEGF and its receptors and co-receptors initiate the intracellular signaling necessary to promote vessel sprouting, and ultimately, the formation of fully matured and functional vessels. The new vasculature formed following VEGF signaling enables delivery of oxygen and nutrients and facilitates removal of waste products2.

Regulating angiogenesis presents an attractive treatment strategy for diseases characterized by either insufficient or excessive vascularization. In the context of excessive vascularization seen in many types of cancer, inhibiting angiogenesis can decrease tumor growth. Anti-angiogenic treatment targeting tumor vascularization is a particular focus area within cancer research3. One anti-angiogenic drug is bevacizumab, a recombinant monoclonal antibody that neutralizes VEGF (an “anti-VEGF” drug). Bevacizumab is approved as a monotherapy or in combination with chemotherapy for several cancers, including metastatic colorectal cancer, non-small cell lung cancer, and metastatic cervical cancer4. In 2008, the drug gained accelerated approval for treatment of metastatic breast cancer (mBC) through the US Food and Drug Administration (FDA). Though initial clinical trials initially showed that bevacizumab improved progression-free survival (PFS), subsequent results revealed that bevacizumab failed to improve overall survival (OS) in a wide range of patients and that the drug elicited significant adverse side effects5. Consequently, the FDA revoked its approval for the use of bevacizumab for mBC in late 20116.

The case of bevacizumab illustrates that although anti-angiogenic therapy can be effective, not all patients or cancer types respond to the treatment. This underscores the need for biomarkers that can help select patients who are likely to respond to anti-angiogenic treatment. Bevacizumab, and anti-angiogenic therapies in general, typically confer a modest overall survival benefit when combined with chemotherapy7. The prevailing hypothesis is that anti-angiogenic agents induce a restructuring and normalization of the tumor vasculature, allowing for faster and more efficient drug delivery7. However, some patients respond better than others, and it is therefore beneficial to establish markers that predict whether the treatment will benefit specific patient populations8.

Numerous studies have sought to identify biomarkers for anti-angiogenic treatment. Biomarkers can be used to determine which tumors will respond prior to any treatment being given (“predictive”), or to evaluate efficacy following treatment (“prognostic”)8. Biomarkers can also be used to determine optimal doses, to design combination therapies, and to identify resistance to therapies9. The concentration range of circulating angiogenic factors (CAFs), and VEGF in particular, is one possible predictor of the response to anti-angiogenic therapy8. Alternatively, expression of angiogenic receptors such as NRP1 and VEGFR1 on tumor cells, in the tumor interstitial space, or in plasma can serve as biomarker candidates5,10. Unfortunately, though some of these candidates are promising, a marker that predicts bevacizumab treatment outcome has not yet been validated5,8. In fact, relying on the concentrations of CAFs has produced inconclusive and inconsistent results9,8,11. Tumor growth kinetics have also been investigated as prognostic biomarkers of the response to anti-angiogenic treatment12,13,14,15,16. The most recent studies take advantage of improved imaging technology that can assess tumor volume, rather than only providing two-dimensional information12. The imaging analyses show that tumor growth kinetics may be a reliable indicator of treatment efficacy and are in good agreement with standardized approaches for assessing response treatment. However, utilizing tumor growth kinetics as a predictive biomarker has not been extensively studied.

In this work, we use a systems biology model to investigate the utility of tumor growth kinetics in predicting response to anti-VEGF treatment. We build upon our existing compartment model of VEGF distribution and kinetics in tumor-bearing mice17 by changing the dynamic tumor volume to be dependent on the pro-angiogenic complexes involving VEGF-bound receptors (the “angiogenic signal”). This new element of the model allows us to simulate anti-VEGF treatment and observe the effect of the treatment on tumor volume. We apply the new model to identify conditions and characteristics of tumor growth that may be predictive of a favorable response to anti-angiogenic treatment. Our work contributes to the identification of validated biomarkers that could be used to determine patients who are likely to respond to anti-angiogenic treatment.

## Results

### Model optimization and validation

Our model is able to recreate the growth dynamics of untreated breast tumor xenografts in mice and can predict the tumor volume in response to anti-angiogenic treatment. We first fit the model to control data from six experimental datasets quantifying tumor volume in mice bearing MDA-MB-231 xenograft tumors without any anti-VEGF treatment18,19,20,21,22. We used nonlinear least squares optimization to fit the model and estimate the optimal parameter values, minimizing the error between the model predictions and the experimental measurements. Specifically, we estimated the values of the tumor growth parameters, including k, k1, *ψ*, and *Ang* (see Methods section for more details). We performed the model fitting 20 times for each of the six datasets, obtaining 20 sets of optimized parameter values per dataset. Overall, the model does a good job of recreating the growth dynamics of untreated tumors (Figure 1A-F, blue lines). In fact, the final data point is within the 95% confidence interval for four out of the six datasets. One limitation is the fit to data from Volk *et al.*22, where the model fails to capture the sigmoidal shape of the experimental tumor growth curve (Figure 1F). We moved forward by validating the model with data not used in the fitting. Here, using the same fitted kinetic parameters as the control case, we simulated the treatment regimens described in the *in vivo* mouse studies (Supplemental Table S1). The model does an excellent job of matching the experimental data (Figure 1), capturing the effect of anti-VEGF treatment on tumor growth for the majority of datasets. Based on these results, the model is in agreement with experimental data of untreated tumor growth and was appropriately validated using treatment data.

![Figure 1.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F1.medium.gif)

[Figure 1.](http://biorxiv.org/content/early/2017/05/10/136531/F1)

Figure 1. Model fit and validation using full tumor growth time course for fitting.
Predicted tumor volume over time for the six datasets. **A,** Roland18. **B,** Zibara19. **C,** Tan21. **D,** Volk20. **E,** Volk 2011a22. **F,** Volk 2011b22. The model is able to reproduce experimental data for tumor growth without treatment and predict validation data not used in parameter fitting. Blue triangles and purple squares are control and treatment experimental data points, respectively. Solid lines are the mean of the model results obtained from the best fits, and the shading indicates the 95% confidence interval. Note different scales on both axes.

The optimized parameter values from the fits with the lowest errors are given in Figure 2. When fitting to the datasets from Volk et al., the model fitting and parameter estimation showed higher values and k/k1 ratios than the other three datasets. Some of the estimated parameter values varied widely, even up to 3-4 orders of magnitude, while others, such as, had a much narrower range.

![Figure 2.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F2.medium.gif)

[Figure 2.](http://biorxiv.org/content/early/2017/05/10/136531/F2)

Figure 2. Estimated model parameters obtained from fitting.
The estimated parameter values from the best fits are plotted for each dataset. **A,** k. **B,** k1. **C,. D,. E,** k/k1. Horizontal bar indicates the mean of the best fits obtained from fitting the model to each dataset. Statistical comparison of the estimated parameter sets is given in Figure 5.

The model can recreate both control and treatment dynamics even when parameter fitting is performed using a limited number of experimental measurements. We wanted to examine whether it is possible to accurately predict the response to anti-VEGF treatment when the model fitting only includes the initial tumor growth data. We selected the datasets that included at least three tumor volume measurements prior to administration of bevacizumab (two out of the six datasets fit this criterion). We fit those initial experimental data points for the control (no anti-VEGF treatment) and validated the fitted model using the anti-VEGF treatment data. We again performed the model fitting 20 times for each dataset. The optimized model fit using only the initial tumor growth data was able to predict the tumor volume following treatment (Figure 3). Although the 95% intervals were wider in this fitting as compared to the results obtained when all of the data points were used for model fitting (see Figure 1), the newly optimized model still predicted reasonable values for the tumor size. In general, the error was less than 50% between the simulated and experimental tumor volumes. Interestingly, the confidence intervals for the treatment cases were actually smaller than those for the control cases.

![Figure 3.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F3.medium.gif)

[Figure 3.](http://biorxiv.org/content/early/2017/05/10/136531/F3)

Figure 3. Model validation after fitting initial tumor growth data.
Predicted tumor volume over time for the two datasets with at least three pre-treatment measurements for tumor volume. **A,** Roland18. **B,** Volk 2011a22. Triangles (blue) and squares (purple) are control and treatment experimental data points, respectively. Only the blue data points are used for fitting. Solid line is the mean of the predicted results obtained from the best fits, and the shading shows the 95% confidence interval on the best fits. Note different scales on both axes. Confidence intervals for the predicted tumor volumes with treatment were reasonably small and contained the experimental data points.

We performed statistical analyses to compare the estimated parameter values obtained from fitting the full set of control data to the estimated parameters obtained from fitting only the initial measurements before treatment. The analyses revealed no significant difference between the estimated parameter values obtained from the two fittings (*p* > 0.9 for all parameters). Thus, the fitted model obtained using only the initial time points provides reliable predictions regarding the effect of anti-VEGF treatment on tumor volume and is consistent with the optimized model obtained from fitting all of the data. Overall, the model predicts tumor growth relatively well, even with a limited number of pre-treatment data points used for model fitting. This highlights the robustness of the model in that only three to five data points were needed in order to get a reasonable prediction for treatment response. However, for the remainder of the paper, we present simulations obtained using the optimized parameter sets estimated from fitting the full time course for the control data. Using the fitting from the full training data allows us to compare the response to treatment across all six independent experimental datasets.

### Predicting the response to anti-VEGF treatment

Having validated our model, we used the optimized parameter sets to predict the tumor volume in response to anti-VEGF treatment. We ran the model for each of the six datasets, using all 20 sets of optimized parameter values. For each set of parameters, the model was simulated for three cases: no treatment (control) and two treatment conditions (2 and 10 mg/kg bevacizumab). For the treatment cases, twice-weekly injections were simulated, starting when the tumor volume reached 0.1 cm3 (termed “*T**start*”). We selected this volume, since it is thought to be the critical time at which tumors typically start secreting higher levels of angiogenic factors in order to recruit the vasculature necessary to support further growth (∼1-2 mm in diameter). For all cases, the model was simulated for 6 weeks after *T**start*. We used the model to predict the relative tumor volume (RTV), the ratio of the final tumor volume for the control and treatment cases: ![Formula][1]</img>  where *V**treatment* and *V**control* are the tumor volumes at the end of the simulation with treatment and without treatment, respectively. Thus, the RTV represents the fold-change in tumor size due to treatment, where an RTV value less than one indicates that the treatment reduced the tumor volume, compared to the control. We use the RTV value to characterize the response to anti-VEGF treatment. The predicted response to bevacizumab treatment at a dose of 10 mg/kg using the best fit parameter values is shown in Figure 4. The range of predicted RTV values indicates that certain tumors are more responsive to anti-VEGF treatment than others. In particular, the predicted RTV values obtained using fitted parameter values from fitting to data from Volk are higher than the predicted response for the other datasets. We next performed a thorough statistical comparison of the RTV and the estimated parameter values obtained in the fitting.

![Figure 4.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F4.medium.gif)

[Figure 4.](http://biorxiv.org/content/early/2017/05/10/136531/F4)

Figure 4. Predicted response to anti-VEGF treatment.
The model was used to simulate bevacizumab treatment at a dose of 10 mg/kg. The relative tumor volume (RTV) predicted by the model is shown. Horizontal bar indicates the mean of the predicted RTV for the best fits from each dataset.

Our statistical analysis indicates a relationship between particular kinetic parameters that characterize tumor growth and the effectiveness of treatment. We used to statistical analyses to determine whether the sets of estimated parameters or the predicted RTV values were statistically significantly different (*p* < 0.05) across the six datasets (Figure 5). Based on this analysis, we found that all datasets with significantly different predicted RTV values either had significantly different k, , values or k/k1 ratios. The parameter represents the switch from exponential to linear growth9, and k/k1 is the ratio of the exponential to linear growth rates23. Interestingly, there was no statistically significant difference in the estimated values, the “basal angiogenic signal”, between any of the datasets. Overall, the statistical analysis reveals that certain kinetic parameters (k, , and k/k1) varied considerably between datasets and corresponded to significantly different treatment response (as indicated by the RTV value). The values of those parameters, which characterize the kinetics of tumor growth, may be used to predict the response to treatment.

![Figure 5.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F5.medium.gif)

[Figure 5.](http://biorxiv.org/content/early/2017/05/10/136531/F5)

Figure 5. Statistical analysis of the optimized parameter sets.
Standard ANOVA analysis, followed by pairwise comparisons, was used to determine whether the sets of optimized parameter values were statistically different. **A,** upper triangle: k; lower triangle: k1. **B,** upper triangle: *ψ* lower triangle: *Ang*. **C,** upper triangle: k/k1; lower triangle: RTV for bevacizumab dose of 10 mg/kg. The color indicates log10(*p*-value): **\***|**, (*p*-value ≤ 0.001); ******, (0.001 < *p*-value ≤ 0.01); *****, (0.01 < *p*-value < 0.05).

### Determination of relationship between tumor growth parameters and response to treatment

We applied PLSR, a multivariate regression analysis, to further quantify the importance of specific tumor growth characteristics in predicting the response to anti-VEGF treatment. We used the values of k, k1, *ψ*, *Ang*, and k/k1 as inputs (predictors) and the RTV at the two dosage levels for bevacizumab (2 and 10 mg/kg) as the responses. We determined the optimal PLSR model by varying the number of components from one to six and calculating the R2X, R2Y, and Q2Y values (see Methods section). We also varied the number of inputs, using different combinations of the estimated parameters. The final PLSR model (i.e., the model that best predicted the responses without over-fitting) had two components and included four inputs (k, k1, *ψ*, and *Ang*). This PLSR model is able to accurately predict the RTV at both dosage levels (Figure 6A) and also performs well with leave-one-out cross validation (Q2Y = 0.93). We analyzed the PLSR model to obtain insight regarding how the four inputs relate to the outputs. The variable importance of projection (VIP) scores for the four model inputs indicate that the value of *ψ* is the largest contributor to predicting the RTV (Figure 6B). This suggests that the value of *ψ* could be used to distinguish tumors that will respond to therapy or not.

![Figure 6.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F6.medium.gif)

[Figure 6.](http://biorxiv.org/content/early/2017/05/10/136531/F6)

Figure 6. Results from multivariate analysis.
PLSR analysis quantifies how the tumor growth parameters influence the response to treatment (RTV). **A,** PLSR model to predict RTV for two dosage levels of the anti-VEGF. **B,** VIP scores for the model inputs; a score greater than one indicate variables that are important for predicting the RTV. **C,** Scores of the model output; increasing in component one corresponds to decreasing efficacy. **D**, Loadings of the model inputs.

Although the PLSR components do not explicitly correspond to a physiological variable, plotting the loadings for the inputs and outputs provides some insight into the meaning of each component. A plot of the loadings for the outputs (Figure 6C) reveals that both components capture the treatment efficacy. Decreasing in component 1 and increasing in component 2 corresponds to increased efficacy of the anti-VEGF treatment. Here, one of the datasets from the 2011 paper by Volk *et al.*22, in which anti-VEGF treatment is the least effective in reducing tumor growth compared to the other datasets, has the highest loading in component 1 and lowest loading in component 2 (i.e., it appears in the lower right portion of the plot). In comparison, measurements from tumors in which anti-VEGF treatment leads to more growth inhibition appear in the upper left quadrant of the plot. Here, we consider both components, as together, they account for 80% of the variance in the output.

A plot of the loadings for the inputs reveals how the estimated tumor growth parameters are associated with treatment efficacy. Here, we focus only on the loadings for component 1, as this component accounts for 94% of the variance in the inputs. We find that *ψ* is positively correlated with low treatment efficacy, as it has a positive loading in component 1 (Figure 6D). This result suggests that a high value of *ψ* is associated with low treatment efficacy. In summary, the multivariate analysis provides a regression model that accurately predicts the relative tumor volume following anti-VEGF treatment, given the tumor growth parameters. Additionally, the analysis confirms the importance of *ψ* as a key predictor of the tumor’s response to anti-VEGF treatment.

### Effect of tumor receptor number on the response to treatment

The tumor microenvironment is predicted to influence the response to anti-VEGF treatment. After validating the model and investigating relationships between kinetic parameters describing tumor growth and response to treatment, we sought to investigate the effects of the tumor microenvironment. In particular, we examined the effect of neuropilin and VEGF receptor levels on relative tumor volume. VEGF receptor levels were varied from 0 to 10,000 receptors/cell, and NRP levels were varied from 0 to 100,000 receptors/cell. Using the best set of parameters from each dataset (i.e., the parameter set with the lowest error), we used the model to determine *T**start* for each combination of receptor levels. We then ran the model to simulate the tumor growth for six weeks past *T**start* to obtain the baseline control volumes. Treatment volumes were obtained by simulating twice-weekly bevacizumab injections at a dose of 10 mg/kg for six weeks after *T**start*. The RTV values were calculated for each combination of the tumor receptor densities. The model predicts that higher neuropilin levels led to decreased treatment efficacy, especially for high VEGFR levels (Figure 7). Regardless of neuropilin density, very high treatment efficacy (low RTV) was observed when tumor cells expressed few VEGFRs. The predicted RTV values obtained using the estimated parameters from certain datasets show that neuropilin expression has a noticeable impact on the response to treatment (Figure 7A-B). In comparison, neuropilin levels seem to have a diminished impact for the Volk dataset, indicated by surfaces that are almost identical, even with drastic changes in neuropilin receptor levels (Figure 7C). In summary, the model can be used to study tumor microenvironments that provide favorable conditions for anti-angiogenic treatment. Higher receptor expression is predicted to reduce anti-VEGF efficacy, although the relationship was not uniform across all datasets, indicating the importance of accounting for specific tumor microenvironmental conditions.

![Figure 7.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F7.medium.gif)

[Figure 7.](http://biorxiv.org/content/early/2017/05/10/136531/F7)

Figure 7. Effect of VEGF receptor expression on tumor cells.
Model predictions based on optimized parameter values obtained from fitting: **A,** Roland18. **B,** Zibara19. **C,** Volk 200820. Surface plots reveal the relationship between RTV and VEGFR1, VEGFR2, and neuropilin receptor density on tumor cells. Red color indicates higher RTV, representing tumor conditions that are less favorable for anti-VEGF treatment. Neuropilin density varies: 0 receptors/cell (*left*), 20,000 receptors/cell (*center*), and 100,000 receptors/cell (*right*). The colorbar indicates the RTV value, with the same range for all panels.

## Discussion

We have developed a compartmental model representing tumor-bearing mice in which the tumor volume is responsive to changes in VEGF concentration. The tumor volume explicitly depends on the “angiogenic signal”, which is the signal produced when VEGF binds to its receptors on tumor endothelial cells. In this way, the model can be applied to analyze the effect of anti-VEGF treatment on xenograft tumor growth, aiding in the analysis of pre-clinical data. Kinetic parameters are obtained by fitting the model to experimental data of breast xenograft tumor growth in mice and are validated with treatment data. By including a dynamic tumor volume that explicitly depends on the concentration of VEGF-bound receptors, we address a primary limitation of our previous work.

The change in tumor volume estimated by our computational model matches experimental data used to train the model. We determine the tumor growth kinetic parameters by fitting the computational model to tumor growth measurements obtained from *in vivo* tumor xenografts. Using the fitted model, we predict the tumor volume in response to anti-angiogenic treatment targeting VEGF. These model predictions compare well with measurements of treatment response observed in six different pre-clinical *in vivo* experiments. Our approach of training the model using control data and using the optimized model to predict treatment data is a significant advantage over previous modeling work. In model fitting performed by other groups, tumor growth parameters were estimated by simultaneously fitting both control and treatment groups23. In contrast, our computational model is able to accurately predict response to anti-VEGF treatment, data not used in the fitting.

The model provides unique insight into how certain kinetic parameters that characterize tumor growth correlate with response to anti-angiogenic treatment. Our results demonstrate how the parameters describing tumor growth could be used as a predictive biomarker for treatment response. In comparison, other studies have used volume-based growth tumor kinetics as a prognostic biomarker. Lee and coworkers found that the time to progression (defined as the time it takes the tumor to grow from its nadir in volume after treatment to a progressive disease state) was significantly correlated with overall survival12. In other work, researchers used tumor growth kinetics to determine the efficacy of anti-angiogenic treatment13–16. Excitingly, our approach is highly predictive, where volumetric measurements performed prior to treatment can give insight into how the tumor might respond to an anti-VEGF agent such as bevacizumab.

We performed various analyses to quantify how the tumor growth kinetic parameters influence the response to treatment. The PLSR and statistical analyses reveal that higher *ψ* values are related to decreased treatment efficacy. The parameter *ψ* represents the transition from exponential to linear growth9. Thus, according to our results, anti-angiogenic treatment targeting VEGF would be more effective in tumors that have a smoother transition from exponential to linear growth (low *ψ*). The fitted parameters related to the shape of the tumor growth curve (k and k1) were not shown to influence the response to treatment as significantly as *ψ*. However, for datasets with significantly different responses to anti-VEGF treatment, the k/k1 ratio is also significantly different. Our statistical analyses indicate an inverse relationship between the ratio k/k1 and effectiveness of treatment. Simeoni *et al.* posit that k and k1 may be indicative of the initial aggressiveness of the cell line and of the response of the animal to tumor progression (i.e., immunological or anti-angiogenic response), respectively23. According to this interpretation, treatment would be least effective for tumors with aggressive initial growth (high k) combined with a strong response from the animal (low k1). Additionally, we find that the basal angiogenic signal, *Ang*, is not predictive of anti-angiogenic treatment response. This agrees with experimental results indicating that the ability of basal levels of circulating angiogenic factors to predict treatment efficacy is limited8.

We used the model to investigate how the number of VEGF receptor and co-receptors on tumor cells influences the response to treatment. Currently, modified expression of VEGF receptors (VEGFR1, VEGFR2, or NRP1) appears to be among the most promising markers for bevacizumab treatment, though this has not consistently been replicated across different studies involving various cancer types5. In particular, low levels of soluble VEGFR1 expression in plasma and NRP1 expression on tumor cells are characteristics of a bevacizumab-responsive tumor5. Therefore, we wanted to use our model to predict the influence of tumor microenvironment on treatment efficacy. The model predicts that low levels of NRP1 or VEGFR lead to increased treatment efficacy for all datasets. The treatment is predicted to be most effective when both NRP and VEGFR have low expression; however, the effect of VEGFR levels appears to be more pronounced. That is, the treatment was still effective for high NRP levels, as long as VEGFR levels were low. These results are in agreement with other biomarker studies10,24. Although there was a consistent inverse relationship between receptor levels and treatment efficacy, this extent to which receptor numbers influenced the predicted relative tumor volume was not identical for all tumors. Datasets for tumors with higher k/k1 ratios and *ψ* values had higher RTV (i.e., the treatment was less effective), even for a wide range of receptor expression levels. This may indicate that intrinsic characteristics of the tumor related to its growth kinetics make anti-angiogenic treatment less effective, regardless of microenvironmental tumor conditions. As a result, solely using receptor expression as a predictive biomarker could lead to inconsistent results across tumor types.

Our model accurately predicts tumor growth profiles and matches experimental measurements. The focus of our model is on the molecular level interactions occurring between VEGF and its receptors. In our model, the number of VEGF-receptor (pro-angiogenic) signaling complexes formed directly influences tumor growth. We acknowledge that this representation of tumor growth omits the many intracellular signaling pathways and corresponding cellular-level responses (i.e., proliferation and migration) involved in new blood vessel formation. However, the model does indeed capture the dynamics of tumor growth, providing a mechanistic understanding of the growth kinetics that contribute to the response to anti-VEGF treatment.

We acknowledge some assumptions and limitations that may be addressed as more quantitative data become available. For example, we do not account for changes in tumor vascularity relative to the tumor volume. The tumor volume consists of interstitial space, vascular volume, and tumor cells. We account for tumor growth by assuming the tumor cell volume fraction increases, while the interstitial space volume fraction decreases, and the relative proportion of the vascular volume is constant (see Methods section for more detail). This means that the tumor vascularity does change as the overall tumor volume grows, but it remains in the same proportion relative to the whole tumor volume. Furthermore, we do not simulate remodeling of the blood compartment or changes in vascular permeability in response to anti-VEGF treatment. However, experimental data show a decrease in microvessel density following bevacizumab treatment25, and incorporating this observation would enhance the model. Additionally, anti-angiogenic treatment promotes normalization of the vasculature, which allows for more efficient delivery of chemotherapy to the tumor7. Accounting for changes in the microvascular density would allow us to simulate combination treatments that include chemotherapy and anti-angiogenic agents. Unfortunately, there is a lack of robust time-series data that can be used to predict changes in vascular density with treatment. This limitation may be addressed as additional quantitative measurements are published.

The model is highly successful in capturing the growth kinetics of exponential or linear growth curves. However, the model does not accurately predict sigmoidal tumor growth. The equation governing tumor growth used in our model is based on the foundational work of Simeoni *et al.*, who adapted a Gompertz model of tumor growth to investigate both the exponential and linear phases of growth23. Although this makes the tumor growth equation more flexible, it also limits the ability to simulate an eventual plateau in growth. The model’s inability to capture sigmoidal growth was particularly apparent when fitting the Volk datasets20,22. However, we have focused on exponential growth, as it has been implemented in many other mathematical models26,27 and shown to accurately fit tumor growth data28. Expansion of the tumor growth equation can be added in future studies.

## Conclusion

We constructed a computational model that simulates the kinetics of VEGF binding to its receptors and the influence of VEGF-bound receptor complexes on tumor volume in tumor-bearing mice. This model is a useful tool in the analysis of pre-clinical data. The model matches control tumor growth data used for fitting, and it was validated using data not used in fitting. The validated model accurately predicts the tumor growth upon administration of anti-angiogenic treatment that targets VEGF. Recently published imaging studies indicate that tumor growth kinetics are prognostic biomarkers for treatment efficacy. Interestingly, the fitted parameter values estimated in the present study point to the possibility of using tumor growth kinetics as a predictive biomarker for anti-angiogenic treatment. This model also helps to elucidate why biomarker candidates such as expression of VEGF receptors on tumor cells may not be reliable for all tumors. Although the model predicts that receptor levels influence response to treatment, the results are not uniform across all of the experimental datasets we analyzed. Thus, our modeling work lays the foundation for future studies to investigate the importance of tumor growth kinetics as a predictive and specific biomarker and can accelerate the discovery of biomarker candidates in pre-clinical studies.

## Materials and Methods

### Computational Model

#### Compartmental model

In this work, we expand on our previous three-compartment model17 by including VEGF-mediated tumor growth. We briefly describe the full model and detail the new additions that are the focus of this work. The model is comprised of three compartments representing the whole mouse: normal tissue (assumed to be skeletal muscle), blood, and tumor (Figure 8). The model includes human and mouse VEGF isoforms: human isoforms (VEGF121 and VEGF165) are secreted by tumor cells, and mouse isoforms (VEGF120 and VEGF164) are secreted by endothelial cells in the normal, blood, and tumor compartments and muscle fibers in the normal tissue. Geometric characteristics, receptor densities, kinetic parameters, and transport rates are all detailed in our previous paper17.

![Figure 8.](http://biorxiv.org/https://www.biorxiv.org/content/biorxiv/early/2017/05/10/136531/F8.medium.gif)

[Figure 8.](http://biorxiv.org/content/early/2017/05/10/136531/F8)

Figure 8. Model schematic.
The computational model includes three compartments: normal tissue, blood, and tumor volume. The compartments are connected via lymphatic flow from the interstitial space in the normal tissue to the blood and transendothelial macromolecular permeability. Molecular species include human and mouse VEGF isoforms, VEGF receptors and co-receptors (including the soluble receptor VEGFR1, sR1), and the protease inhibitor α-2-macroglobulin (a2m). Glycosaminoglycan (GAG) chains represent the extracellular matrix. The volume of the tumor depends on the concentration of receptor-bound VEGF complexes on tumor endothelial cells (denoted as [rec-VEGF]).

#### Tumor volume and growth

Previously, we assumed the tumor volume increased exponentially with time, based on measurements from tumor xenografts17. Under that assumption, cancer treatment, including anti-angiogenic therapy, has no effect on tumor growth.

In the present study, we address that limitation by introducing an equation for tumor growth wherein the volume of the tumor compartment is dependent on the “angiogenic signal” (*Ang*) produced when VEGF binds to its receptors on endothelial cells in the tumor.

The tumor compartment is assumed to consist of cancer cells, endothelial cells (vascular volume) and interstitial space, each of which has a defined volume fraction (i.e., volume relative to the total tumor volume). Our previous model assumed that as the total tumor volume increased, the relative proportions of cancer cells, vascular space, and interstitial space remain constant. Here, we still have the volume fraction for the vascular space remaining constant, based on a range of experimental data29,30,31. However, we used results from a recent imaging study to account for an increase in the relative volume of cancer cells as the tumor volume increases. Christensen and coworkers measure how tumor cell density increases as the tumor grows by tracking cancer cells in xenograft tumors in rats using near near-infrared (NIR) fluorescence dyes32. The authors quantify the fluorescence intensity in a tumor and use it to estimate the number of cancer cells as the tumor grows over time. The estimated cell count was normalized by the tumor volume to obtain the number of cancer cells per unit volume of tumor tissue as the tumor grows. We extracted the values obtained by Christensen and coworkers for MDA-MB-231 tumors and converted them to the cancer cell volume fraction using the volume of tumor cells, as we have done in our previous work17. Therefore, we have been able to incorporate into our model an increase in the cancer cell volume fraction over time. Assuming a tumor cell volume of 905 µm3, based on our previous examination of the literature17, we developed expressions describing the decay of interstitial space during tumor growth. We found that the relative decrease in interstitial space during tumor growth was adequately modeled by exponential decay.

Tumor growth was described by an adapted Gompertz model focusing on the exponential and linear phases of the tumor, as previously described9,23. The differential equation for the tumor volume is: ![Formula][2]</img> 

Here, *v(t)*is the tumor volume in cm3 at time *t*, k and k1 are parameters describing the rate of exponential and linear growth, respectively. The units of k and k1 are s-1 and cm3 tissue/s, respectively. The parameter represents the transition from exponential to linear tumor growth and is unitless. The *Ang* parameter represents the basal angiogenic signal (at time *t* = 0), and *Ang(t)*is the angiogenic signal at time *t*. The value of *Ang* at any time is calculated as the total concentration of pro-angiogenic VEGF-receptor complexes on tumor endothelial cells. This includes VEGFR1 and VEGFR2 bound to either mouse or human VEGF isoforms, with or without the NRP1 co-receptor. Thus, *Ang(t)*and *Ang* have units of concentration (mol/cm3 tissue). The values of the tumor growth parameters (k, k1, *ψ*, and *Ang*) were estimated by fitting the model to experimental data, as described below.

### Parameter estimation

#### Model training

Six datasets from independent *in vivo* experimental studies of MDA-MB-231 xenograft tumors in mice were used for parameter fitting and validation. Each of the datasets included growth profiles for untreated tumors (control), as well as tumors treated with the anti-VEGF agent bevacizumab. First, k, k1, *ψ*, and *Ang* (“free parameters”) were fit using the control tumor growth profiles for each dataset. Fitting was performed using the *lsqnonlin* function in MATLAB. We used this algorithm to minimize the sum of squared residuals (SSR): ![Formula][3]</img>  where *v**exp,i* is the *i*th experimentally measured tumor volume, *v**sim,i* is the *i*th simulated volume at the corresponding time point, and *n* is the total number of experimental measurements. The minimization is subject to Θ, the set of upper and lower bounds on each of the free parameters.

We found that weighting the residual by the experimental measurement biased the error towards early data points and reduced the model’s ability to fit the full course of tumor growth. Therefore, we minimized the residual, with no weighting, to fit the model to the experimental data.

We performed the parameter fitting 20 times for each dataset. To attempt to arrive at a global minimum for the error, we initialized each fitting run by randomly selecting a value for the free parameters within the specified upper and lower bounds. The bounds were set such that the range for each parameter was at least one order of magnitude: 10-8 to 10-2 for k and k1, 0.1 to 50 for *ψ*, and 10-16 to 10-14 for *Ang*. After performing the model fitting, we used the SSR to identify the optimal parameters. Parameter sets with the smallest errors were taken to be the “best” fits and were used for subsequent statistical analysis. The number of “best” parameter sets varied between datasets and ranged from 12 to 16 parameter sets. We first tested to see whether there were significant effects of the experimental data being fit on the estimated parameters values using one-way non-parametric ANOVA. This method makes no assumptions about the distributions of parameter values and tests whether samples originate from a common distribution. We then performed post-hoc analyses to make pairwise comparisons using the Kruskal-Wallis test. We corrected for multiple comparisons using statistical hypothesis testing (Dunn’s test). All statistical analyses were performed using GraphPad Prism.

Two of the experimental datasets contained at least three data points prior to administration of treatment18,20. These points were used in a separate model fitting to see whether limiting the data used for model training to only pre-treatment measurements could generate a fitted model that still accurately predicts the response to anti-angiogenic treatment.

#### Model validation

After fitting the control data, we validated the estimated parameters with data not used in the fitting. We applied the fitted model to simulate anti-angiogenic treatment and compared the predicted tumor growth profile to the experimental measurements for the treatment cases. Here, we simulated the dosing regimens used in each experiment with the same optimized parameters obtained by fitting the control data. For each dataset, we simulated intravenous injections lasting for one minute. All six experimental studies administered bevacizumab bi-weekly; however, the dosage varied between the studies. The dosing regimens are given in Supplemental Table S1. The binding affinity and clearance rate for bevacizumab were obtained from experimental studies in which VEGF was immobilized on a flow cell (FC) and bevacizumab was injected over the FC at varying concentrations33. Based on those measurements, the binding affinity was set to 4456 pM (*k**on* = 2.19 × 104 M-1s-1; *k**off* = 9.8 × 10-5 s-1), and 5.73 × 10-7 s-1 was used for the anti-VEGF clearance rate.

### Partial least squares regression analysis

Partial least squares regression (PLSR) modeling was used to determine the relationship between the fitted parameters characterizing tumor growth kinetics (inputs) and the response to treatment given by the RTV value (output). PLSR modeling seeks to maximize the correlation between the inputs and outputs. To accomplish this, the inputs and outputs are recast onto new dimensions called principal components (PCs), which are linear combinations of the inputs. The regression coefficients estimated by PLSR describe the relative importance of each input. Quantitative measures from the PLSR modeling, including the loadings and scores, provide some insight into the biological meaning of the PCs34. Additionally, we use the estimated regression coefficients to determine each input’s contribution across all responses. This metric is given by the “variable importance of projection” (VIP) for each predictor. The VIP value is the weighted sum of each input’s contribution to all of the responses. As such, the VIP values indicate the overall importance of the predictors. VIP values greater than one indicate variables that are important for predicting the output response.

For our analysis, the input matrix was 6 rows × 4 columns, where the 6 rows correspond to the best fit for each of the six datasets, and the 4 columns consisted of the estimated free parameters (k, k1, *ψ*, and *Ang*). The output matrix was 6 rows × 2 columns, where the rows corresponds to the predicted RTV using the best fit for each of the six datasets, and the columns are the two treatment doses (2 and 10 mg/kg). We used two metrics to evaluate the model fitness: R2Y and Q2Y, which each have a maximum value of 1. The R2Y value indicates how well the model fits the output data. The Q2Y metric specifies how much of the variation in the output data the model predicts35, and values greater than 0.5 indicate that the model can predict data not used in the fitting. We performed PLSR modeling using the nonlinear iterative partial least squares (NIPALS) algorithm36, implemented in MATLAB (Mathworks, Inc.).

### Numerical implementation

All model equations were implemented in MATLAB using the SimBiology toolbox. The final model is provided as SBML (Supplementary File S1). Parameter fitting was performed using the *lsqnonlin* function MATLAB. GraphPad Prism was used to run statistical analyses on parameter values.

## Author Contributions

SDF designed the research. TDG and SDF performed the research, analyzed the data, and wrote the manuscript. The authors thank members of the Finley research group for critical comments and suggestions.

## Additional Information

### Competing financial interests

The authors declare no competing financial interests.

## Acknowledgements

The research was supported by NSF CAREER award 1552065 (SDF) and a USC Provost’s Undergraduate Research Fellowship (TDG).

*   Received May 10, 2017.
*   Accepted May 10, 2017.


*   © 2017, Posted by Cold Spring Harbor Laboratory

This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/)

## References

1.  1.Olsson, A., Dimberg, A., Kreuger, J. & Claesson-Welsh, L. VEGF receptor signalling — in control of vascular function. Nat. Rev. Mol. Cell Biol. 7 359–371 (2006).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1038/nrm1911&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=16633338&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000237057500015&link_type=ISI) 

2.  2.Carmeliet, P. & Jain, R. K. Molecular Mechanisms and and clinical applications of angiogenesis. Nature 473 298–307 (2011).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1038/nature10144&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=21593862&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000290722400036&link_type=ISI) 

3.  3.Al-Husein, B., Abdalla, M., Trepte, M., DeRemer, D. L. & Somanath, P. R. Antiangiogenic therapy for cancer: An update. Pharmacotherapy 32 1095–1111 (2012).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1002/phar.1147&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=23208836&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

4.  4.FDA Approval for Bevacizumab - National Cancer Institute. at [http://www.cancer.gov/about-cancer/treatment/drugs/fda-bevacizumab](http://www.cancer.gov/about-cancer/treatment/drugs/fda-bevacizumab)
    
    

5.  5.Lambrechts, D., Lenz, H. J., De Haas, S., Carmeliet, P. & Scherer, S. J. Markers of response for the antiangiogenic agent bevacizumab. J. Clin. Oncol. 31 1219–1230 (2013).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamNvIjtzOjU6InJlc2lkIjtzOjk6IjMxLzkvMTIxOSI7czo0OiJhdG9tIjtzOjM3OiIvYmlvcnhpdi9lYXJseS8yMDE3LzA1LzEwLzEzNjUzMS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 

6.  6.Alberto, M. J., Escobar, M., Lopes, G., Glück, S. & Vogel, C. Bevacizumab in the Treatment of Metastatic Breast Cancer: Friend or Foe? Curr. Oncol. Rep. 14 1–11 (2012).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1007/s11912-011-0202-z&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=22012632&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

7.  7.Jain, R. K. Normalizing tumor microenvironment to treat cancer: bench to bedside to biomarkers. J. Clin. Oncol. 31 2205–2218 (2013).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamNvIjtzOjU6InJlc2lkIjtzOjEwOiIzMS8xNy8yMjA1IjtzOjQ6ImF0b20iO3M6Mzc6Ii9iaW9yeGl2L2Vhcmx5LzIwMTcvMDUvMTAvMTM2NTMxLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 

8.  8.Jain, R. K. et al. Biomarkers of response and resistance to antiangiogenic therapy. Nat. Rev. Clin. Oncol. 6 327–338 (2009).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1038/nrclinonc.2009.63&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=19483739&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000266528500008&link_type=ISI) 

9.  9.Sharan, S. & Woo, S. Quantitative Insight in Utilizing Circulating Angiogenic Factors as Biomarkers for Antiangiogenic Therapy: Systems Pharmacology Approach. CPT Pharmacometrics Syst. Pharmacol. 3 e139 (2014).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1038/psp.2014.36&link_type=DOI) 

10. 10.Van Cutsem, E. 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. 30 2119–2127 (2012).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamNvIjtzOjU6InJlc2lkIjtzOjEwOiIzMC8xNy8yMTE5IjtzOjQ6ImF0b20iO3M6Mzc6Ii9iaW9yeGl2L2Vhcmx5LzIwMTcvMDUvMTAvMTM2NTMxLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 

11. 11.Kopetz, S. et al. Phase II trial of infusional fluorouracil, irinotecan, and bevacizumab for metastatic colorectal cancer: Efficacy and circulating angiogenic biomarkers associated with therapeutic resistance. J. Clin. Oncol. 28 453–459 (2010).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamNvIjtzOjU6InJlc2lkIjtzOjg6IjI4LzMvNDUzIjtzOjQ6ImF0b20iO3M6Mzc6Ii9iaW9yeGl2L2Vhcmx5LzIwMTcvMDUvMTAvMTM2NTMxLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 

12. 12.Lee, J. H. et al. Volume-based growth tumor kinetics as a prognostic biomarker for patients with EGFR mutant lung adenocarcinoma undergoing EGFR tyrosine kinase inhibitor therapy: a case control study. Cancer Imaging 16 5 (2016).
    
    

13. 13.Seyal, A. R. et al. Performance of tumor growth kinetics as an imaging biomarker for response assessment in colorectal liver metastases: correlation with FDG PET. Abdom. Imaging 40 3043–3051 (2015).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1007/s00261-015-0546-1&link_type=DOI) 

14. 14.El Sharouni, S. Y., Kal, H. B. & Battermann, J. J. Accelerated regrowth of non-small-cell lung tumours after induction chemotherapy. Br. J. Cancer 89 2184–9 (2003).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1038/sj.bjc.6601418&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=14676792&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000188093600004&link_type=ISI) 

15. 15.Stein, W. D., Yang, J., Bates, S. E. & Fojo, T. Bevacizumab reduces the growth rate constants of renal carcinomas: a novel algorithm suggests early discontinuation of bevacizumab resulted in a lack of survival advantage. Oncologist 13 1055–1062 (2008).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTM6InRoZW9uY29sb2dpc3QiO3M6NToicmVzaWQiO3M6MTA6IjEzLzEwLzEwNTUiO3M6NDoiYXRvbSI7czozNzoiL2Jpb3J4aXYvZWFybHkvMjAxNy8wNS8xMC8xMzY1MzEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

16. 16.Rezai, P. et al. Change in the growth rate of localized pancreatic adenocarcinoma in response to gemcitabine, bevacizumab, and radiation therapy on MDCT. Int. J. Radiat. Oncol. Biol. Phys. 81 452–459 (2011).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1016/j.ijrobp.2010.05.060&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=21570199&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

17. 17.Finley, S. D., Dhar, M. & Popel, A. S. Compartment model predicts VEGF secretion and investigates the effects of VEGF trap in tumor-bearing mice. Front. Oncol. 3 196 (2013).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.3389/fonc.2013.00196&link_type=DOI) 

18. 18.Roland, C. L. et al. Inhibition of vascular endothelial growth factor reduces angiogenesis and modulates immune cell infiltration of orthotopic breast cancer xenografts. Mol. Cancer Ther. 8 1761–1771 (2009).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTA6Im1vbGNhbnRoZXIiO3M6NToicmVzaWQiO3M6ODoiOC83LzE3NjEiO3M6NDoiYXRvbSI7czozNzoiL2Jpb3J4aXYvZWFybHkvMjAxNy8wNS8xMC8xMzY1MzEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

19. 19.Zibara, K. et al. Anti-angiogenesis therapy and gap junction inhibition reduce MDA-MB-231 breast cancer cell invasion and metastasis in vitro and in vivo. Sci. Rep. 5 12598 (2015).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1038/srep12598&link_type=DOI) 

20. 20.Volk, L. D. et al. Nab-paclitaxel efficacy in the orthotopic model of human breast cancer is significantly enhanced by concurrent anti-vascular endothelial growth factor A therapy. Neoplasia 10 613–623 (2008).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1593/neo.08302&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=18516298&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000255885200010&link_type=ISI) 

21. 21.Tan, G. et al. Combination therapy of oncolytic herpes simplex virus HF10 and bevacizumab against experimental model of human breast carcinoma xenograft. Int. J. Cancer 136 1718–1730 (2015).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1002/ijc.29163&link_type=DOI) 

22. 22.Volk, L. D. et al. Synergy of nab-paclitaxel and bevacizumab in eradicating large orthotopic breast tumors and preexisting metastases. Neoplasia 13 327–338 (2011).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1593/neo.101490&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=21472137&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000289529400004&link_type=ISI) 

23. 23.Simeoni, M. et al. Predictive Pharmacokinetic-Pharmacodynamic Modeling of Tumor Growth Kinetics in Xenograft Models after Administration of Anticancer Agents Predictive Pharmacokinetic-Pharmacodynamic Modeling of Tumor Growth Kinetics in Xenograft Models after Administratio. Cancer Res. 64 1094–1101 (2004).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NjoiY2FucmVzIjtzOjU6InJlc2lkIjtzOjk6IjY0LzMvMTA5NCI7czo0OiJhdG9tIjtzOjM3OiIvYmlvcnhpdi9lYXJseS8yMDE3LzA1LzEwLzEzNjUzMS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 

24. 24.Jubb, A. M. et al. Impact of exploratory biomarkers on the treatment effect of bevacizumab in metastatic breast cancer. Clin. Cancer Res. 17 372–381 (2011).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTA6ImNsaW5jYW5yZXMiO3M6NToicmVzaWQiO3M6ODoiMTcvMi8zNzIiO3M6NDoiYXRvbSI7czozNzoiL2Jpb3J4aXYvZWFybHkvMjAxNy8wNS8xMC8xMzY1MzEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

25. 25.Willett, C. G. et al. Direct evidence that the VEGF-specific antibody bevacizumab has antivascular effects in human rectal cancer. Nat Med 10 145–147 (2004).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1038/nm988&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=14745444&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000188719600032&link_type=ISI) 

26. 26.Enderling, H. & Chaplain, M. A. J. Mathematical Modeling of Tumor Growth and Treatment. Current Pharmaceutical Design 20 4934–4940 (2014).
    
    

27. 27.Hahnfeldt, P., Panigrahy, D., Folkman, J. & Hlatky, L. Tumor Development under Angiogenic Signaling. Cancer Res. 59 4770 LP–4775 (1999).
    
    

28. 28.Benzekry, S. et al. Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth. PLOS Comput. Biol. 10 e1003800 (2014).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1371/journal.pcbi.1003800&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=25167199&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

29. 29.Lewin, M. et al. In vivo assessment of vascular endothelial growth factor-induced angiogenesis. Int J Cancer 83 798–802 (1999).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1002/(SICI)1097-0215(19991210)83:6<798::AID-IJC16>3.0.CO;2-W&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=10597197&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

30. 30.Bogin, L., Margalit, R., Mispelter, J. & Degani, H. Parametric imaging of tumor perfusion using flow- and permeability-limited tracers. J. Magn. Reson. Imaging 16 289–299 (2002).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1002/jmri.10159&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=12205585&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

31. 31.Cao, M., Liang, Y., Shen, C., Miller, K. D. & Stantz, K. M. Developing DCE-CT to Quantify Intra-Tumor Heterogeneity in Breast Tumors With Differing Angiogenic Phenotype. IEEE Trans. Med. Imaging 28 861–871 (2009).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1109/TMI.2008.2012035&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=19150783&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

32. 32.Christensen, J., Vonwil, D. & Prasad, V. S. Non-invasive in vivo imaging and quantification of tumor growth and metastasis in rats using cells expressing far-red fluorescence protein. PLoS One 10 1–14 (2015).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0131080&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=25919488&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 

33. 33.Yang, J. et al. Comparison of binding characteristics and in vitro activities of three inhibitors of vascular endothelial growth factor A. Mol. Pharm. 11 3421–3430 (2014).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1021/mp500160v&link_type=DOI) 

34. 34.Kreeger, P. K. Using partial least squares regression to analyze cellular response data. Sci. Signal. 6 tr7 (2013).
    
    [Abstract/FREE Full Text](http://biorxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6ODoic2lndHJhbnMiO3M6NToicmVzaWQiO3M6OToiNi8yNzEvdHI3IjtzOjQ6ImF0b20iO3M6Mzc6Ii9iaW9yeGl2L2Vhcmx5LzIwMTcvMDUvMTAvMTM2NTMxLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 

35. 35.Prasasya, R. D., Vang, K. Z. & Kreeger, P. K. A Multivariate Model of ErbB Network Composition Predicts Ovarian Cancer Cell Response to Canertinib. Biotechnol. Bioeng. 109 213–224 (2012).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1002/bit.23297&link_type=DOI) 
    
    [PubMed](http://biorxiv.org/lookup/external-ref?access_num=21830205&link_type=MED&atom=%2Fbiorxiv%2Fearly%2F2017%2F05%2F10%2F136531.atom) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=000297110300022&link_type=ISI) 

36. 36.Geladi, P. & Kowalski, B. R. Partial least-squares regression: a tutorial. Anal. Chim. Acta 185 1–17 (1986).
    
    [CrossRef](http://biorxiv.org/lookup/external-ref?access_num=10.1016/0003-2670(86)80028-9&link_type=DOI) 
    
    [Web of Science](http://biorxiv.org/lookup/external-ref?access_num=A1986E477500001&link_type=ISI)

 [1]: /embed/graphic-4.gif
 [2]: /embed/graphic-10.gif
 [3]: /embed/graphic-11.gif