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 k0, k1, ψ, and Ang0 (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.
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 k0/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.
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.
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 “Tstart”). 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 Tstart. We used the model to predict the relative tumor volume (RTV), the ratio of the final tumor volume for the control and treatment cases: where Vtreatment and Vcontrol 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.
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 k0, , values or k0/k1 ratios. The parameter represents the switch from exponential to linear growth9, and k0/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 (k0, , and k0/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.
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 k0, k1, ψ, Ang0, and k0/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 (k0, k1, ψ, and Ang0). 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.
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 Tstart for each combination of receptor levels. We then ran the model to simulate the tumor growth for six weeks past Tstart 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 Tstart. 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.
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 (k0 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 k0/k1 ratio is also significantly different. Our statistical analyses indicate an inverse relationship between the ratio k0/k1 and effectiveness of treatment. Simeoni et al. posit that k0 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 k0) combined with a strong response from the animal (low k1). Additionally, we find that the basal angiogenic signal, Ang0, 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 k0/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.
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:
Here, v(t)is the tumor volume in cm3 at time t, k0 and k1 are parameters describing the rate of exponential and linear growth, respectively. The units of k0 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 Ang0 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 Ang0 have units of concentration (mol/cm3 tissue). The values of the tumor growth parameters (k0, k1, ψ, and Ang0) 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, k0, k1, ψ, and Ang0 (“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): where vexp,i is the ith experimentally measured tumor volume, vsim,i is the ith 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 k0 and k1, 0.1 to 50 for ψ, and 10-16 to 10-14 for Ang0. 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 (kon = 2.19 × 104 M-1s-1; koff = 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 (k0, k1, ψ, and Ang0). 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).