Mechanistic insights into the heterogeneous response to anti-VEGF treatment in tumors

Vascular endothelial growth factor (VEGF) is a strong promoter of angiogenesis in tumors, and anti-VEGF treatment, such as a humanized antibody to VEGF, is clinically used as a monotherapy or in combination with chemotherapy to treat cancer patients. However, this approach is not effective in all patients or cancer types. To better understand the heterogeneous responses to anti-VEGF and the synergy between anti-VEGF and other anticancer therapies, we constructed a computational model characterizing angiogenesis-mediated growth of in vivo mouse tumor xenografts. The model captures VEGF-mediated cross-talk between tumor cells and endothelial cells and is able to predict the details of molecular- and cellular-level dynamics. The model predictions of tumor growth in response to anti-VEGF closely match the quantitative measurements from multiple preclinical mouse studies. We applied the model to investigate the effects of VEGF-targeted treatment on tumor cells and endothelial cells. We identified that tumors with lower tumor cell growth rate and higher carrying capacity have a stronger response to anti-VEGF treatment. The predictions indicate that the variation of tumor cell growth rate can be a main reason for the experimentally observed heterogeneous response to anti-VEGF. In addition, our simulation results suggest a new synergy mechanism where anticancer therapy can enhance anti-VEGF simply through reducing the tumor cell growth rate. Overall, this work generates novel insights into the heterogeneous response to anti-VEGF treatment and the synergy of anti-VEGF with other therapies, providing a tool that be further used to test and optimize anticancer therapy.


INTRODUCTION
growth factor (VEGF), which includes five isoforms generated from alternative splicing of mRNA: VEGF-A, VEGF-B, VEGF-C, VEGF-D, and placental growth factor (PlGF). 3 The angiogenic effects of VEGF are primarily initiated through the VEGF-A ligand, 4,5 which binds to two tyrosine kinase receptors, VEGFR1 and VEGFR2, and recruits co-receptors called neuropilins (NRP1 and NRP2). The VEGF receptors and co-receptors were thought to only be expressed on endothelial cells (ECs). However, in recent studies, they are reported to be expressed on many other cell types, including macrophages, hematopoietic stem cells, neurons, muscle fibers, and various tumor cell types. [5][6][7] VEGF binding to its receptors and co-receptors initiates cellular processes, including proliferation, migration, survival, and vessel maturation, which then leads to new blood vessel formation in tumor sites. VEGF is reported to be associated with malignancy, metastasis and poor survival rate for a variety of solid tumors, including breast, prostate, lung, brain, ovarian and pancreatic. [8][9][10][11][12] Therefore, VEGF and its receptors have been intensively studied as the target of anti-angiogenic therapy to hinder tumor angiogenesis and growth. 10,13 The most notable anti-angiogenic drug is bevacizumab (Avastin R ; a humanized anti-VEGF-A antibody). 14 It was proven to be effective in prolonging survival and delaying tumor progression in patients with various cancers, 15 leading to it being approved for use as a monotherapy or in combination with chemotherapy for several cancers, 16 including metastatic colorectal cancer (mCRC), non-small cell lung cancer metastatic renal cell carcinoma, epithelial ovarian cancer and metastatic breast cancer (mBC). However, the variations between tumors (known as intertumor heterogeneity) can lead to a range of responses to anti-VEGF treatment. Later studies revealed that not all patients or cancer types respond to bevacizumab treatment. In late 2011, the US Food and Drug Administration (FDA) revoked the approval of bevacizumab for the treatment of metastatic breast cancer because it was not shown to be effective and safe for mBC patients. 14,17 To date, bevacizumab is still used for mCRC patients, although differences between colorectal tumors also raise concerns as to heterogeneous responses that are observed clinically. 18,19 The heterogeneous response to bevacizumab has inspired substantial computational modeling work to systematically study tumor angiogenesis and further optimize anti-VEGF treatment. 20,21 Models of intracellular signaling of VEGF characterize the biochemical events inside the cell that mediate the signal transduction necessary for cell proliferation and migration. 22 These models help in the identification of intracellular targets for mitigating or blocking VEGF signaling. In addition, models can be used to capture extracellular interactions, 23 including ligand-receptor binding, receptor-coreceptor coupling and drug-target binding events. These models were mainly used to illustrate the distribution of angiogenic factors and therapeutic agents in tumor tissue and in the whole body. 24,25 In addition, the ligand-receptor details are further linked to higher-level processes to study tumor angiogenesis and response to anti-VEGF treatment. For example, the VEGF ligand-receptor dynamics are linked with different cellular processes, including proliferation, apoptosis, and maturation, and influence cellular crosstalk and the development of vascular phenotypes. 12,26 The VEGF ligand-receptor details are also connected with tumor volume dynamics to study the in vivo response to anti-VEGF treatment in several modeling studies. [27][28][29] Overall, the computational models provide insights into tumor angiogenesis and the effects of anti-angiogenic therapy across different scales. In this study, we aim to understand how the variations of properties related to tumor growth lead to the heterogeneous response to anti-VEGF treatment, which has not been the focus of previous works.
The work presented here significantly builds upon two previous modeling efforts. 12,28 Specifically, we followed the modeling approach used in a previous study 12 to describe tumor angiogenesis, tumor growth, and response to treatment across molecular and cellular scales. Our model captures the cross-talk between endothelial cells and tumor cells and is able to predict the details of molecular signaling and cellular dynamics, in a more simplified manner than the previous paper. A novel aspect of our study is to fit and validate the simplified model to multiple tumor growth experimental datasets from various sources, while in previous work by Jain and Jackson, 12 the model was fit to a single dataset from Fujita et al. 30 without model validation. Mouse xenograft models of human tumors, in which human cancer cell lines are used to establish a tumor in the mouse, are broadly utilized in preclinical studies to evaluate the efficacy of anti-VEGF. Importantly, the observed tumor growth and response vary across mice, even when injected with the same cancer cell line and under the same experimental conditions. Thus, inspired by our previous work, 28,31 we trained the constructed model with quantitative measurements from multiple preclinical mouse xenograft studies. 30,[32][33][34][35][36][37] In model validations, the predictions of tumor growth in response to anti-VEGF treatment closely match experimental data, which shows the great utility and predictive power of the model. With the constructed model, we investigated several unexplored aspects of anti-VEGF treatment. Specifically, we study the effects of anti-VEGF treatment on endothelial and tumor cells individually, how tumor-specific parameters (tumor cell growth rate and carrying capacity) influence the response to anti-VEGF treatment, and how anti-VEGF treatment synergizes with other cancer therapies. F I G U R E 1 Computational model of in vivo angiogenesis and tumor growth. Tumor and endothelial cells express VEGF receptors (R1 and R2). VEGF binding to R2 promotes cell proliferation. Tumor cells secrete VEGF and promote hypoxia, which further enhances tumor cell VEGF secretion. Endothelial cells inhibit hypoxia, and hypoxia inhibits tumor cell proliferation. The anti-VEGF and chemotherapeutic agents reach the tumor through a constant transport rate (from the blood).
Finally, we applied the model to study different experimentally tested combination therapy strategies of anti-VEGF and chemotherapy. Altogether, this work generates mechanistic insights into the heterogeneous response and the synergistic effect of anti-VEGF and provides useful guidance for future preclinical and clinical practices.

Modeling VEGF-mediated angiogenesis and tumor growth
Our model describes in vivo angiogenesis and tumor growth in mouse xenografts. Figure 1 presents a schematic of the model. The model characterizes a simplified homogeneous tumor site consisting of two cellular species: tumor cells (T) and endothelial cells (E). The free VEGF receptors (R1 and R2) are expressed on both the tumor cell surface (R T1 and R T2 ) and on the endothelial cell surface (R E1 and R E2 ). Following the approach used in the previous study, 12 we used a scaled ratio of endothelial cell number (E) over tumor cell number (T) as the indicator of the degree of hypoxia in the tumor site. It is reported that a growing tumor, is comprised of 9.2% vessel cells and has an oxygen level of 10 mm Hg oxygen. 38,39 Assuming a linear relationship between oxygen level and the percent-age of endothelial cells, 12  ) × 100%. In the model, tumor cells grow through cell proliferation and will drive the formation of hypoxia (lower the H level). An increase in the hypoxia degree will promote the secretion of VEGF (V) from tumor cells and inhibit the growth of tumor cells. For the endothelial cells, they are assumed to undergo programmed cell death, apoptosis, and will be removed from the system at a constant rate. The proliferation of endothelial cells relieves hypoxia and consequently promotes the growth of tumor cells. The secreted VEGF can bind to both R1 and R2 and form ligand-receptor complexes on both the tumor cell surface (C T1 and C T2 ) and the endothelial cell surface (C E1 and C E2 ). VEGF receptors have unique functions and influence angiogenesis in different ways. Previous work assumes that R2 is the signaling receptor of endothelial cells and R1 is the signaling receptors of tumor cells. 12 However, functionally active VEGFR2 has been identified on cancer cells, 6,7 and VEGFR2 is commonly considered as the predominant receptor that promotes pro-angiogenic VEGF signaling pathway. 3,40 The VEGFR1 is reported to primarily regulate angiogenesis by ligand-trapping and receptor dimerization. 3,40 To better reflect this, we assumed that, for both cell types, the R2 is the primary signaling receptor and R1 only serves as nonsignaling binding sites in the model. Therefore, only the formed VEGFR2 complexes (C T2 and C E2 ) induce cell proliferation. In addition, we assumed the tumor cell growth and endothelial cells growth are subject to different physical constraints. The tumor cells growth is limited by the carrying capacity, which is a tumor-specific property. The endothelial cells' growth is limited by the number of tumor cells (T) to reflect that the endothelial cells are growing together with tumor cells within the tumor site.

Mathematical equations of the model
The model is represented by a set of coupled non-linear ordinary differential equations (ODEs). The details of the model equations are included in Supplementary File 1. We implemented the model with MATLAB (The MathWorks, Natick, MA, USA). The model files are provided in supplementary files as well. The ODEs were solved using the ode15s solver in MATLAB to predict the species' quantities over time.

Model calibration and validation
There are 38 parameters presented in our model, including molecular dynamic parameters, cellular response parameters, and drug pharmacokinetic/pharmacodynamic (PK/PD) parameters. The full list of parameters is reported in the Supplementary File 1. Most of the parameters are derived from previous studies. In our work, we focus on estimating four parameters, including two global parameters for VEGF-mediated angiogenesis (e grow and μ E ) and two dataset-specific parameters for tumor growth kinetics (t grow and t cap ). The four parameters are estimated by fitting the model to tumor volume data. Experimental data from various mouse xenograft studies were extracted using the WebPlotDigitizer program 41 and are shown in Supplementary File 2. The data extracted from published studies indicate that tumor growth is variable in each dataset, even though they were generated using the same breast cancer cell line and similar experimental protocols. Several reasons can contribute to the discrepancies, including the mouse strain, intrinsic properties of the cancer cells, number of the tumor cells injected, the location of the injection site, and the method used to quantify tumor volume. To account for these dataset-specific variations, the tumor cell growth rate (t grow ) and carrying capacity (t cap ) are selected as dataset-specific parameters, meaning they have different values for different datasets to capture the variation in tumor growth dynamics across studies. We must estimate the dataset-specific parameters individually for each dataset first, and then use the estimated parame-ters to predict the response to anti-VEGF. Parameters were estimated with the Particle Swarm Optimization (PSO) algorithm, which minimizes the sum of squared residuals (SSR). 42 We experimented with different numbers of dataset-specific parameters and performed sensitivity analysis (Fig. S1). We settled upon fitting t grow and t cap , as other combinations led to overfitting or underfitting the data.
We collected seven published datasets of breast cancer that used the MDA-MB-231 cell line to form xenografts in nude mice (data and descriptions are included in Supplementary File 2). [32][33][34][35][36][37] These data are used to establish the baseline model. The process of parameter estimation and model validation is illustrated in Fig. S2. To estimate the parameters and validate the model, we separate the seven datasets into training data (six datasets) and test data (one dataset, e.g. ELHajjar et al.). Firstly, in step 1, we use the six training datasets to estimate the values of global parameters and the values of dataset-specific parameters for the training datasets. In step 2, we use the control data from the test dataset to estimate the values of datasetspecific parameters for that dataset. To do so, we use the global parameters estimated from the first step. In step 3, we validate the estimated model parameters by predicting the anti-VEGF response with the parameter values from the first and second steps and comparing to the treatment group data in the test dataset. In this way, we are able to confirm the global parameters and dataset-specific parameters for the test dataset.
To examine the bias introduced from the selection of training data, we applied leave-one-out cross validation to investigate the predictive capability of the model. Specifically, we trained and validated the model seven times. Each time, a different dataset from these seven datasets was selected as the test dataset with the other six datasets used for training. After the cross-validation, we further validate the model on one dataset of breast cancer for the MDA-MB-435 cell line 36 and one dataset of head and neck squamous cell carcinomas with the UM-SCC-17B cell line 30 to investigate the generalizability of the model. For this validation, the global parameters are set to the values estimated with MDA-MB-231 datasets. The dataset-specific parameters are estimated with the control group data. The treatment group data is used to validate estimated datasetspecific parameters.
In addition to the anti-VEGF therapy, our model includes chemotherapy to simulate the administration of Nab-paclitaxel (ABX). For the chemotherapy module, we use a dataset of the ABX plasma concentration-time profile 43 in a mouse xenograft model to estimate the PK parameters. The parameters for the cell-killing effect were estimated by fitting the model to tumor volume datasets in which ABX was administered in mouse xenografts formed using the MDA-MB-231 cell line. 36 Then we compare the model prediction for various combination therapies with experimental data from breast cancer xenografts formed using the MDA-MD-435 cell line 36 to validate the prediction of the combination therapy.

Quantification of treatment response
The details of model construction are included in Supplementary File 1. To predict the tumor volume change, the predicted numbers of tumor cells (TCs) and endothelial cells (ECs) are converted to tumor volume following the approach used in the previous study 12 : tumor volume in μl is V = (volume of 1 million ECs)×E + (volume of 1 million TCs)×T, where E and T are the predicted number of ECs and TCs, respectively. The volume of 1 EC is reported to be 2.2 × 10 −6 µl 44 and that of 1 TC is 1 × 10 −6 µl. 45 With the developed model, we can generate the tumor volume dynamics for the control and treatment groups.
We used the simulated tumor volume dynamics to quantify two different indicators of response to anti-VEGF therapy. Firstly, we computed the area under the curve (AUC) of simulated tumor volume and termed this the tumor volume response (R AUC ). Specifically, for different combinations of growth rate and carrying capacity, we simulated the control and treatment condition for 90 days (assuming a three-month observation window, similar to experimental studies). For the treatment condition, we simulated administering 4 mg/kg bevacizumab twice a week for two weeks, once the tumor volume reached 0.5 ml. To quantify the effect of treatment, we compare the AUC of the simulated tumor volume for the control case to the treatment case. We calculated the response as = − .
In addition to R AUC , we defined the half-time response, R Half-time . To calculate R Half-time , we simulate the control and treatment conditions until steady state is reached and record the time required for the tumor volume to reach 50% of the final tumor volume. The difference in the halftime is defined as half-time response: R Half-time = Halftime treatment -Half-time ctrl , which indicates the effect of treatment in delaying the time to reach half of the final tumor volume.
For the survival estimates, we simulate a heterogeneous in silico mouse population, by generating 100 sets of values randomly selected from a uniform distribution within the range of the parameter estimations for the two tumor-specific parameters (t grow and t cap ). We predicted the tumor volume dynamics for each virtual mouse (a unique combination of t grow and t cap ) and applied timeto-event analysis to determine the survival of the whole mouse population. Following the approach used in previous study, 31 an in silico mouse is recorded as 'sacrificed' when its tumor reaches 2 cm 3 within the simulated time. Alternatively, a mouse is recorded as 'censored' at a particular time point, t, if its predicted tumor volume remains below 2 cm 3 but ends at a time t that is before the end of the predefined simulation time. All other mice are retained in the study and recorded as 'alive'. Survival curves were estimated by the Kaplan-Meier method in GraphPad Prism (GraphPad Software, San Diego, CA, USA).

Model validation
The cross-validation is used to examine the prediction capacity of the model calibrated with different splits of training and validation data. We applied leave-one-out cross validation with the seven collected datasets. Each time, we held out one dataset as a test dataset and use the other six datasets as training datasets. The training datasets are used to estimate the two global parameters (e grow and μ E ) and dataset-specific parameters (tumor cell growth rate t grow and carrying capacity t cap ) for training datasets.
With the estimated values for global parameters, we use the control group data from the test dataset to estimate the dataset-specific parameters for the test dataset. We also use the predicted tumor volume dynamic under anti-VEGF treatment compared to experimental data from the test dataset for model calibration. Figure 2 shows the results of cross-validation.  36 and one dataset of head and neck squamous cell carcinomas from the UM-SCC-17B cell line. 30 As described above, the control data is used to estimate dataset-specific parameters, while the treatment group data is used to validate the model. Figure 3 shows that the model prediction also quantitatively matches the experimental data for these tumor xenografts. With a carefully constructed model validated with experimental data, we can apply the model to investigate the effects of anti-VEGF strategies, alone, and in combination with chemotherapeutic agents.

Anti-VEGF inhibits tumor growth primarily through blocking VEGF pro-angiogenic effect
With the constructed model, we simulate in vivo tumor growth under different hypothetical therapeutic con-ditions to help us better understand the mechanism of action of anti-VEGF treatment on tumor and endothelial cells. 4 VEGF binds to its receptors on both tumor cells and endothelial cells to promote cell proliferation. To understand the cell type specific effect, we simulated five different conditions, given in Figure 4 treatment condition affects the cell such that ligand-bound R2 does not contribute to cell proliferation. Overall, MuteE is fully blocking the VEGF pro-angiogenic effect and MuteT is fully blocking the VEGF pro-tumorigenic effect. The MuteET is fully blocking both the pro-angiogenic and pro-tumorigenic effects, which sets the maximum efficacy for the anti-VEGF. We reported the predicted tumor volume changes ( Figure 4A) and the percentage change of endothelial cells to reflect the dynamics of the degree of vasculature ( Figure 4B).
The results show that anti-VEGF treatments reduce the degree of vasculature to inhibit the tumor growth. We present three cases with representative anti-VEGF regimens in Figure 4, and Figure S5 includes the simulation results for other dataset conditions. The results show that the curves of anti-VEGF are below the control group, which indicates the anti-VEGF inhibits tumor growth (Row A), corresponding to the reduced degree of vasculature (Row B). Continuous bevacizumab treatment can delay the tumor growth dynamics and reduce the final tumor volume at a steady state (Figure 4, column I).
In comparison, a short period of bevacizumab treatment cannot reduce the final tumor volume (Figure 4, columns II and III). This is because the drug is degraded and cannot promote a significant effect if given over a short period of time. The overlapping of the MuteT (magenta) and Control (blue) curves implies that MuteT has a similar effect as no treatment at all, which indicates that blocking the pro-tumorigenic effect of VEGF does not significantly contribute to inhibiting tumor growth. The MuteE curve (brown) is significantly lower than Control (blue) in all figures and very close to the MuteET (black), which proves that blocking VEGF-driven endothelial cell proliferation can strongly reduce tumor angiogenesis and inhibit the tumor growth. MuteET is the lowest tumor volume curve in all figures, which means completely blocking both tumor cells and endothelial cells has the strongest therapeutic response. The curves of the anti-VEGF group (red) are above the MuteET group (black), which indicates the anti-VEGF only partially blocks VEGF signaling in these conditions. Together, the results suggest that anti-VEGF inhibits the tumor growth primarily through targeting ECs, blocking the pro-angiogenic effect of the VEGF.

F I G U R E 4 Predicted outcome of anti-VEGF treatment and hypothetical treatments target tumor and endothelial cells. A, Predicted tumor volume. B, Number of endothelial cells relative to the total number of cells. I, Simulation with the estimated dataset-specific parameters and
anti-VEGF treatment protocol for Volk2011a. 36 II., Simulation with the estimated dataset specific parameters and anti-VEGF treatment protocol for Tan. 34 III., Simulation with the estimated dataset-specific parameters and anti-VEGF treatment protocol for Zibara. 33 Grey shaded area indicates the time period for anti-VEGF treatment.

The tumor cells' growth rate affects the response to anti-VEGF treatment
We quantified the impact of two tumor-specific parameters, tumor cells growth rate t grow and carrying capacity t cap , on the response of anti-VEGF treatment. This is related to our previous works, where we identified tumor growth kinetics as biomarkers for the outcomes of anti-VEGF treatments. The estimated tumor-specific parameters could take on a range of values, depending on the dataset being fit: the tumor cell growth rates ranged from 0.1 to 10 s -1 , and carrying capacities ranged from 1 × 10 9 to 2 × 10 10 cells. Thus, we predicted the tumor volume and vasculature degree changes for various combinations of the tumor-specific parameters with and without anti-VEGF treatment (Fig. S6). For the anti-VEGF treatment, we simulated two weeks of treatment in which bevacizumab was administered twice per week (two bolus injections at Day 1 and Day 4) at a dose of 4 mg/kg (a dosage commonly used in experimental studies). 36 The treatment starts when the tumor volume reaches 0.5 cm 3 . With the simulated dynamics, we quantified two responses of the anti-VEGF (R AUC and R Half-time ), as mentioned in Section 2.4. We note that in the model, t grow and t cap are assumed to be constant tumor-specific properties and are not affected by treatment. Specifically, the tumor cell growth rate (t grow ) discussed here is a model parameter representing the basal tumor cell growth rate.
The model predicts that the tumor cell growth rate (t grow ) and carrying capacity (t cap ) significantly affect the outcome of anti-VEGF treatment ( Figure 5). The color change from red to blue along the diagonal direction (from upper left to lower right) in Figure 5A shows that the tumor volume AUC response is stronger in tumors with lower tumor cell growth rate and higher carrying capacity. Additionally, the color gradient is more pronounced along the vertical direction (color changing from red to blue in Figure 5A), which implies the variation in tumor cell growth rate can significantly influence the response in conditions with different carrying capacities. This is more evident in Figure 5B (the half time response, R Half-time ), which shows the color mainly changes along the vertical axis. Therefore, the variations of tumor cells' growth rate might be an important contributor to the heterogeneous response of anti-VEGF observed in preclinical and clinical studies. F I G U R E 5 Impact of dataset-specific parameters on the response to anti-VEGF treatment. We vary the dataset-specific parameters and calculate the response to anti-VEGF treatment (R AUC and R Half-time ) as mentioned in Section 2.4. The circles in different colors represent the estimated parameter combinations for different datasets used in previous sections.

F I G U R E 6
Predicted synergy of anti-VEGF and chemotherapy. A, Results from estimating the ABX pharmacokinetic parameters (k X12 , k X21 , and k X23 ) by fitting the model to experimental data. B, Comparing the simulated tumor volume dynamics with the experimental data. The cell killing parameters of ABX are estimated by fitting the model to experimental data (Group ABX and Group bevacizumab + ABX).

Chemotherapy can synergize with anti-VEGF through reducing tumor cell growth rate
Since we observed that the lower tumor cell growth rate generates a stronger response to anti-VEGF, we hypothesized that anticancer therapy that reduces tumor cell growth rate can synergize with anti-VEGF treatment. We included a simplified chemotherapy submodule into the model to simulate the treatment of nab-paclitaxel (ABX). ABX is a cell-cycle-phase specific drug, which targets cells in the G 2 /M phase of the cell cycle, while having a limited effect on cells in G 1 /G 0 phase. 46 Therefore, it has higher killing effect on highly proliferating cells. Here, we assume that ABX only affects the tumor cells. This assumption can be relaxed, and we consider the implications of this assumption in Discussion section. In our model, ABX will trigger a transition of regular tumor cells to damaged cells. The damaged tumor cells do not proliferate and will either recover to regular tumor cells or go through an irreversible cell death process. This cell damage and death model is adapted from previous work 47 and is detailed in Supplementary File 1. Overall, the administered ABX modulates the net growth rate of tumor cells.
The pharmacokinetic parameters of ABX (three parameters) are estimated by fitting the chemotherapy submodule to experimentally measured paclitaxel plasma concentrations. 43 The fitting result is shown in Figure 6A. Then the parameters that characterize ABX's cell killing effects (six parameters) are estimated by fitting the tumor growth model combined with the chemotherapy module to tumor volume data for treatments involving ABX (Figure 6B, blue and brown curves). 36 The model simulations precisely capture the observed synergy of bevacizumab and ABX in the experimental data ( Figure 6B). The combination treatment group ( Figure 6B, brown) shows a significant regression of tumor volume, while the ABX group alone (blue) leads to much less tumor volume regression, and bevacizumab alone (red) only delays the increase of tumor volume.
Given this result for the effects of ABX and anti-VEGF combination therapy, we applied the fitted model to predict the response to different ABX/bevacizumab combination treatment regimens. The combination treatment results shown in Figure 6 are for two cycles of ABX and beva-cizumab given together throughout the treatment period. Here, we sought to investigate the effects of alternative combination treatment protocols. The three selected combination treatment protocols (Figure 7) were previously tested in an in vivo mouse study with mice bearing breast cancer MDA-MB-435 xenografts. 36 The model predicts the combination of intermittent ABX with continuous bevacizumab that extends beyond the final dose of ABX (Group 3) has the best overall therapeutic outcome. We also simulated intermittent ABX followed by a period of continuous bevacizumab treatment (Group 1) and a protocol in F I G U R E 8 Kaplan-Meier survival estimates for the heterogeneous virtual mouse population. We generated a heterogeneous population of 100 mice by varying the dataset-specific parameters. The mice were subjected to different treatments: bevacizumab alone, ABX alone, and three combinations of bevacizumab and ABX. The survival for each group is compared to no treatment.
which bevacizumab is given continuously along with intermittent ABX where both treatments end at the same time (Group 2). Compared to Group 1, Group 2 leads to a lower tumor burden for approximately 80 days, but the tumor growth rebounds and has an accelerated growth after bevacizumab treatment is discontinued at day 100. These predictions ( Figure 7A) qualitatively match the trend shown in the experimental data of breast cancer MDA-MB-435 mouse xenografts ( Figure S7). The model also predicts the increase of the degree of vasculature ( Figure 7B) in Group 2 after day 100, which potentially explains the rapid rebound of the tumor volume observed in the experimental data. Specifically, the treatment protocols that reduce both tumor growth and control the size of the endothelial cell population (Groups 1 and 3) have slower tumor regrowth after treatment. Group 3 in particular exhibits strong control of overall tumor growth, as the numbers of endothelial and tumor cells are significantly reduced with treatment.
In summary, bevacizumab alone strongly reduces the number of endothelial cells, while ABX reduces the number of tumor cells. Since they target two distinct cell populations, the combination of ABX and bevacizumab has the strongest effect, and the two drugs are even able to synergize. Altogether, these simulation results show that an anticancer treatment that robustly inhibits growth of both the tumor and endothelial cell populations can lead to tumor regression.

Combination therapy can extend survival time in a virtual mouse population
We next investigated the effects of different combinations of bevacizumab and ABX for a range of tumor-specific properties. We generated a population of 100 virtual mice by uniformly sampling the tumor-specific parameters within the range of the estimated values ( Figure S3). With the virtual mouse population, we simulate various treatment strategies and apply time-to-event analysis to produce the survival data. As shown with the Kaplan-Meier survival estimates (Figure 8), the model predicts that bevacizumab treatment alone only slightly improves the median survival, while chemotherapy alone and other combination therapies can significantly improve the median survival, even increasing survival time by more than two months compared to control. Adding bevacizumab after ABX treatment has ended (Group 1) is less effective compared to Groups 2 and 3 where ABX and bevacizumab are given simultaneously. This suggests that the concurrent combination therapy is important for extending the overall survival time of the population. Among the three combination therapies, the combination of intermittent ABX with simultaneous and continuous bevacizumab extending beyond the end of ABX treatment (Group 3) shows the best therapeutic outcome. Interestingly, the model predicts that ABX and bevacizumab have a synergistic effect on the median survival days. Treatment Group 3 extends the median survival days to 109 days, which improves the survival days by 40 days compared to ABX alone (Group 1). This improvement is significantly larger than the 6-day increase in survival time between the bevacizumab and the control groups.

DISCUSSION
In this work, we constructed a computational model to understand the experimentally observed heterogeneous response of anti-VEGF treatment and its synergy with other anticancer strategies. There are several novel aspects of our modeling approach that help achieve this goal.
Our previous works identify that tumor growth kinetics has a great potential to serve as biomarkers of the outcome of anti-VEGF treatment. The model used in previous work includes 258 species, 129 global parameters, and three tumor-specific parameters estimated with control group data. 28,31 Here, we greatly simplified the model structure and reduced the number of parameters, increasing its utility and ease of understanding. In this work, we build on the structure from Jain et al., 12 removing the intracellular signaling details while retaining the VEGF ligandreceptor interaction network and cellular cross-talk. Our model consists of 20 species and 38 parameters. Importantly, the model is able to explain the same sets of experimental data explored in previous works 28,31 and link the heterogeneous response to easily interpretable growth kinetic parameters. Specifically, we reduced the number of tumor-specific parameters to two, both of which have clear biological meaning. This enables us to investigate the impact of tumor-specific growth properties on the response to anti-VEGF treatment and directly illustrate the connection between tumor growth kinetics and the response to anti-VEGF. Additionally, this model is sufficiently detailed at the cellular level to investigate the effects of interventions that target a specific cell type. Lastly, we included a simple module modeling the administration of nab-paclitaxel to investigate the synergy of a chemotherapeutic agent with anti-VEGF treatment, which is also a significant extension of previous models. Achieving these features in one model is a significant advantage over previous modeling works.
The model provides unique insights into the anti-VEGF therapeutic outcome. With the developed model, we generated a novel understanding of the two-pronged effect of the anti-VEGF on both endothelial cells and tumor cells, which is hard to fully evaluate in the experimental setting. It is reported that VEGF not only promotes the endothelial cells' growth but also increases tumor cell proliferation. 6,7 Therefore, the anti-VEGF can block both the pro-angiogenic and pro-tumorigenic effects of VEGF by influencing endothelial and tumor cells, respectively. Interestingly, our model predictions show that the therapeutic outcome of anti-VEGF is mainly from blocking the pro-angiogenic effect of VEGF, while blocking VEGF-mediated cancer cell proliferation only slightly contributes to inhibiting the tumor volume increase. The parameter values in our model show that the basal growth rate of tumor cells (0.1-14.7 s -1 ) is much higher than the basal growth rate of endothelial cells (0.07 s -1 considering the basal cell proliferation and apoptosis rates together). Activation of the VEGF receptors significantly changes the endothelial cells' growth rate but only slightly contributes to the tumor cells' net growth rate. This explains the prediction of the model that VEGF has a strong impact on inhibiting endothelial cells. It is important to note that this prediction is not influenced by the assumption that tumor cells are more sensitive to VEGF than endothelial cells. In this work, we assumed that the VEGF receptors on tumor cells have the same affinity to VEGF as the receptors expressed on endothelial cells. Therefore, the VEGFR2 fractional occupancy for both cell types is the same value. However, for these two cell types, the same fractional occupancy can have distinct effects on proliferation depending on the values of μ T and μ E . In the baseline model, we set the parameter μ T to be ten times larger than μ E . We explored other assumptions, setting μ T to be 10 times smaller than μ E or even equal to μ E . The model prediction that VEGF more strongly affects endothelial cells rather than tumor cells holds true independently of the relative values of μ T and μ E . In addition, the value of μ T is shown to not strongly influence the model predictions, based on our sensitivity analysis results ( Figure S7; Supplementary File 1). This prediction is useful for understanding the role of VEGF receptors on the tumor cell surface. A direct implication is that therapeutic agents competing for VEGF binding sites on the tumor cell surface might have unintended effects. As we found in previous modeling works, 23,24 anti-angiogenic agents binding to non-signaling binding sites of VEGF might release the VEGF ligand and force it to bind to signaling sites, leading to counter-therapeutic outcomes. Therefore, this model prediction can provide guidance to the preclinical and clinical development of anticancer agents blocking VEGF receptors. 48 The clinically observed heterogeneous tumor response of anti-VEGF treatments is a major drawback to their use. The wide range of responses has prompted the investigation of biomarkers to identify patients for which antiangiogenic therapy will be effective. Our work nicely supplements studies exploring the use of tumor growth kinetics as biomarkers of anti-VEGF response. 31,[49][50][51] Various experimental studies investigated using tumor growth kinetics as a prognostic biomarker, in which the tumor volume dynamics are correlated with survival rate and used to determine the efficacy of anti-angiogenic treatment. 28,31,49,50,52,53 In this work, we quantitatively profiled the impact of two tumor-specific properties of growth kinetics, the tumor cell growth rate and carrying capacity, on the response to anti-VEGF treatment. The model predicted that individual mice with lower tumor cell growth rate and higher carrying capacity have a stronger response to anti-VEGF treatment, which indicates that VEGF treatment might be more beneficial for patients with tumors that have a slow progression but the ability to reach a larger tumor burden. Indeed, the variations of tumor growth rate and carrying capacity are observed in various preclinical and clinical trials, 54 which can be a contributor to the observed heterogeneous response of anti-VEGF. Furthermore, these two growth kinetic parameters can be estimated with early scans of the patients' tumor in clinical practice. 55 Therefore, the model could be used to predict the prognosis and select the patients that will benefit from the treatment.
Additionally, this work proposes a new potential mechanism for the synergy between chemotherapy and anti-VEGF treatment. The combination therapy of anti-VEGF with other drugs has been intensively studied in recent years. 56,57 Various anticancer therapies, including chemotherapy, immunotherapy, and radiation therapy, show a synergistic effect with anti-VEGF treatment in different settings. In this work, we explored the hypothesis that chemotherapy is able to synergize with anti-VEGF simply through reducing the tumor cells' growth rate. To investigate this hypothesis, we simulated administering nab-paclitaxel to modulate the net growth rate of tumor cells combined with simulating anti-VEGF treatment, and we successfully reproduced the synergy observed in experimental data. Commonly, the synergy between chemotherapy therapy and anti-VEGF is explained by the chemotherapeutic agent normalizing the tumor vasculature and modulating the vessel permeability to improve the delivery of the anti-VEGF to the tumor. 58,59 Our work provides a new angle to examine the observed synergy with anti-VEGF. Although we are specifically simulating the administration of nab-paclitaxel, this mechanism can be applied to other anticancer therapies that modulate the tumor cell growth rate. Therefore, the insight generated by the model can guide future experimental or clinical studies of an anti-VEGF in combination with various anticancer therapies.
The model developed in this work can help explore different regimens for combination therapy including anti-VEGF treatment and other agents. For a combination of the same drug agents, different dosing schedules can lead to highly different outcomes. 36 Therefore, arranging dosing sequences in a proper way to maximize the synergy effect has been an active research area of combination therapy. Such studies can lead to important discoveries for the administration protocol and ultimately, greater biological understanding of the effects of the treatment protocols. 59,60 In this work, we compared three combination therapies with different dosing schedules. The model predicts that concurrent therapy of chemotherapy and continuous anti-VEGF treatment (Group 3 in Results Section 3.4) leads to the best treatment outcomes. Importantly, our model not only generates the tumor volume dynamics that match the data, but also provides details on how the degree of vasculature varies, along with the Kaplan-Meier survival estimates. Measuring the dynamics of tumor vasculature and survival estimates is a time-and resourceintensive endeavor. However, these details are necessary for understanding the mechanism behind the tumor volume change. Therefore, our model has great practical interest for preclinical studies, providing a time-and costeffective way to explore different combination strategies to establish proof of concept before moving to in vivo studies.
We acknowledge some limitations of the model that can be addressed in future studies. In our model, we assumed the oxygen level has a linear relationship with the percentage of endothelial cells. However, this assumption could be oversimplified and not physiologically realistic. For example, oxygen diffusion and the physical properties of tumor tissues strongly influence the development of oxygen gradients. Given this complexity, in our studies, we do not directly compare the model prediction to the oxygen measurements reported in experimental studies. There are existing models that can characterize the in vivo oxygen level change in tumor xenografts. 59,61,62 It is possible to combine these models with our model to address this limitation. In addition, endothelial cell maturation and vascular normalization are important processes that have been shown to influence tumor progression and were reported to be affected by the anti-angiogenic treatment; however, they are not included in our model. These aspects have been investigated in previous computational works, 26,63 but the models are not rigorously calibrated due to the scarcity of experimental data. This limitation can be addressed as more quantitative data become available, enabling us to account for vessel maturation. Lastly, to constrain the scope of our work, the influence of tumor angiogenesis on the drug transport between the blood and the tumor is not considered in our model. We assumed a constant transportation rate of the therapeutic agents during tumor progression. Published PK/PD models with more details 64,65 can be incorporated into our model to further improve the accuracy of the drug dose-response prediction. Finally, in the chemotherapy module, we assumed that normal cells (endothelial cells) are not affected by chemotherapy drugs, which might not be true in real condition. We can relax this assumption and easily adapt the model for broader research purposes. Overall, our model provides quantitative mechanistic insight into the effects of VEGF-targeted treatment on tumor cells and endothelial cells. Importantly, the model explains the heterogeneous response of tumor xenografts to anti-VEGF treatment, generates a novel understanding of the synergy of combination therapy and provides a framework for testing various treatment regimens.

A C K N O W L E D G M E N T S
The authors thank members of the Finley research group for critical comments and suggestions. This work was supported by the US National Science Foundation (CAREER Award 1552065 to SDF), the American Cancer Society (130432-RSG-17-133-01-CSM to SDF), and the USC Provost's PhD Fellowship (DL).