## Abstract

Angiogenesis is a crucial step in tumour progression, as this process allows tumours to recruit new blood vessels and obtain oxygen and nutrients to sustain growth. Therefore, inhibiting angiogenesis remains a viable strategy for cancer therapy. However, anti-angiogenic therapy has not proved to be effective in reducing tumour growth across a wide range of tumours, and no reliable predictive biomarkers have been found to determine the efficacy of anti-angiogenic treatment. Using our previously established computational model of tumour-bearing mice, we sought to determine whether tumour growth kinetic parameters could be used to predict the outcome of anti-angiogenic treatment. A model trained with datasets from six *in vivo* mice studies was used to generate a randomized *in silico* tumour-bearing mouse population. We analysed tumour growth in untreated mice (control) and mice treated with an anti-angiogenic agent and determined the Kaplan–Meier survival estimates based on simulated tumour volume data. We found that the ratio between two kinetic parameters, *k*_{0} and *k*_{1}, which characterize the tumour's exponential and linear growth rates, as well as *k*_{1} alone, can be used as prognostic biomarkers of the population survival outcome. Our work demonstrates a robust, quantitative approach for identifying tumour growth kinetic parameters as prognostic biomarkers and serves as a template that can be used to identify other biomarkers for anti-angiogenic treatment.

## 1. Background

Tumour angiogenesis results in the vascularization of a tumour. This process facilitates tumour growth by allowing tumour cells to obtain oxygen and nutrients through the newly formed blood vessels. As excessive vascularization is often seen in many types of cancer, inhibiting angiogenesis is thought to decrease tumour growth. Therefore, anti-angiogenic treatment is pursued as an attractive therapeutic strategy in oncology [1,2].

Bevacizumab is a humanized monoclonal antibody against vascular endothelial growth factor A (VEGF), a key angiogenic promoter in tumours [1]. This drug has been approved as a monotherapy or in combination with chemotherapy for many cancers, including renal cell carcinoma, metastatic colorectal cancer, non-small cell lung cancer and metastatic cervical cancer [3]. It also gained accelerated approval for treatment of metastatic breast cancer through the US Food and Drug Administration (FDA) in 2008. However, subsequent results showed that bevacizumab failed to improve overall survival and that the drug elicited significant adverse side effects. Consequently, the FDA revoked its approval for use of bevacizumab for first-line metastatic breast cancer in late 2011 [4,5]. Several Phase II and III clinical stage studies have also revealed contradicting results regarding the benefit of add-on bevacizumab in the neoadjuvant treatment setting for breast cancer patients [6–11]. Altogether, these studies illustrate that angiogenic therapy may not be effective across a wide range of patients. Indeed, breast cancer is a genetically and clinically heterogeneous cancer type, which makes identifying optimal therapies a challenge [12].

More broadly, there is a need for biomarkers to predict the response to treatment and identify the tumours for which anti-angiogenic treatment will be effective. A number of mechanistic biomarkers have been investigated for their ability to predict response to anti-angiogenic treatment and to determine an optimal treatment strategy. Promising biomarker candidates include the concentration ranges of circulating angiogenic molecules (such as plasma levels of VEGF) [13,14], tissue markers (tumour microvessel density) [15–18] and imaging parameters (magnetic resonance imaging-measured *K*^{trans}) [15,19,20]. However, currently no validated and robust biomarkers are available that can guide selection of patients for whom anti-angiogenic therapy is most beneficial [5,15].

As an alternative, tumour growth kinetics may be used as biomarkers. There is a body of work that investigates how tumour growth kinetics can serve as prognostic biomarkers of the response to anti-angiogenic treatment [21–25]. Recently, a study showed that volume-based tumour growth kinetics may be a reliable indicator of treatment efficacy, and are in good agreement with standardized approaches for assessing response to treatment [21]. Moreover, we developed a computational systems biology model to further investigate the relationship between tumour growth kinetics and the response to anti-angiogenic therapy [26]. The model predicts VEGF distribution and kinetics in tumour-bearing mice, where the dynamic tumour volume is a function of the pro-angiogenic complexes involving VEGF-bound receptors (the ‘angiogenic signal’). By fitting the model to *in vivo* experimental data, we estimated the kinetic parameters that characterize tumour growth. We then used the trained model to predict the effect of anti-VEGF treatment on tumour volume, using only the estimated parameter values. The model predictions of tumour growth in response to anti-VEGF treatment closely matched experimental data. In this study, we concluded that there is a strong correlation between particular intrinsic kinetic parameters and the response to anti-VEGF treatment in terms of the end relative tumour volume (RTV).

Taking advantage of our established model framework and its strong predictive power, we now use this model to further investigate the utility of tumour growth kinetics to serve as a biomarker for anti-angiogenic treatment outcome. We performed an *in silico* randomized mouse study and estimated the survival of tumour-bearing mice in response to anti-VEGF treatment. Here, we introduced variability in the mouse population by allowing the tumour growth kinetic parameter values to vary within defined ranges. A total of 2400 mice with different tumour growth profiles were simulated in this study. By generating this large, heterogeneous *in silico* population of tumour-bearing mice, we can eliminate the likely bias caused by animals dropping out of experimental xenograft studies due to high tumour burden. In general, the average tumour size, particularly in the control group, can be underestimated in an experimental study, thereby underestimating the treatment effect, because large tumours are excluded from the analysis [27]. By contrast, computational modelling avoids these limitations and enables performance metrics (e.g. survival estimates) to be calculated [28]. Furthermore, computational systems biology is a powerful tool for studying how individual components contribute to the function and behaviour of a large system, and has been applied to study cancer at multiple scales [29–31]. Such computational models have been used to identify predictive biomarkers and to enhance the efficacy of anti-angiogenic therapies [13,32,33].

In our previous work, we did not explore how tumour growth parameters affect the response to treatment for individual tumours, nor did we examine time-based tumour growth inhibition. Therefore, in this study, we simulate the tumour volume over time in a heterogeneous population of mice and use more reliable and appropriate read-outs. Specifically, we performed time-to-event analysis [34] by determining the Kaplan–Meier survival curves based on the *in silico* population tumour growth data. We examined tumour growth kinetic parameters as prognostic biomarkers to distinguish the tumour response to anti-angiogenic treatment among the stratified groups.

## 2. Results

### 2.1. *In silico* mouse population tumour growth in the whole-body model

We performed an *in silico* randomized mouse study using our whole-body mouse model (figure 1). The model was previously fitted to each of six independent experimental datasets of control tumour volume in mice bearing MDA-MB-231 xenograft tumours and validated with a separate dataset [26]. The values of *k*_{0} and *k*_{1} (the rates of exponential and linear growth, respectively), and *Ang*_{0} (the basal angiogenic signal at time, *t* = 0) were estimated. A global sensitivity analysis indicated that *ψ* did not significantly influence tumour volume; thus, it was held constant. Here, we simulated the tumour growth of the six *in silico* populations of mice (henceforth referred to as ‘Roland’, ‘Zibara’, ‘Tan’, ‘Volk2008’, ‘Volk2011a’ and ‘Volk2011b’), with and without anti-VEGF treatment, in mice with different tumour growth kinetic parameters. For each population, the values of parameters *k*_{0} and *k*_{1} are randomly varied simultaneously with a uniform distribution within the ranges of their estimated values from our previous model fitting. Previously, a sensitivity analysis showed that the *Ang*_{0} parameter was an influential parameter to the model output when the model was fitted; however, further analysis using partial least-squares regression (PLSR) indicated that *Ang*_{0} was not a strong predictor of response to treatment [26]. Therefore, in each case, *Ang*_{0} is set as the median of the range of its estimated values. We generated 400 *in silico* mice for each of the six cases.

Our simulations show that among the six cases, the anti-VEGF treatment has differential effects in reducing tumour growth, when compared with the control group (figure 2). For all cases, we used a single treatment protocol different from protocols used in each of the six experimental studies, in order to compare the predicted results without bias (termed ‘protocol A’). For Roland, Tan, Volk2008 and Volk2011b (figure 2*a,c,d,f*), the treated tumour volumes are less than the untreated tumours. Meanwhile, for Zibara and Volk2011a (figure 2*b,e*), there is no apparent difference in the tumour volumes for the treated and control groups. Thus, the model simulations reveal distinct differences in the effect of anti-VEGF treatment.

We further studied the effect of anti-VEGF treatment on tumour growth using RTV, the ratio between the mean tumour volumes of the treated and control groups. We calculated the RTV at each time point for all simulated tumours (electronic supplementary material, figure S1). We also determined the RTV at the end of treatment (electronic supplementary material, figure S2). The RTV values in all cases are smaller than one, indicating that the anti-VEGF treatment limits tumour growth, similar to what has been observed experimentally [35–39]. For Zibara and Volk2011a, the endpoint RTV values are just slightly less than one (electronic supplementary material, figure S2B,E), which is an expected result based on the similar tumour growth curves between the control and treated groups (figure 2*b,e*). Comparing the endpoint RTV among all six cases, the effect of anti-VEGF treatment in limiting tumour growth is the strongest for Volk2011b (RTV = 0.459 ± 0.054), followed by Roland (0.454 ± 0.096), Volk2008 (0.615 ± 0.066) and Tan (0.638 ± 0.049). This treatment effect is the least significant in Zibara (0.979 ± 0.009) and Volk2011a (0.987 ± 0.013).

### 2.2. Kinetic parameters as potential predictor for stratified population response

We investigated the relationship between the parameters that characterize tumour growth kinetics and the effect of the anti-VEGF treatment. Previously, our PLSR analysis indicated that for nearly all pairwise comparisons, if the RTV values for two datasets were significantly different, their *k*_{0}/*k*_{1} ratios were also significantly different. This implies that *k*_{0}/*k*_{1} is a large contributor in predicting the endpoint RTV [26]. Additionally, plotting the RTV versus *k*_{0}, *k*_{1}, and *k*_{0}/*k*_{1} shows some relationship between the endpoint RTV and the tumour growth parameters (electronic supplementary material, figure S2). Therefore, we investigated whether these tumour growth parameters could stratify the simulated mouse populations, and distinguish their tumour growth and survival estimates. To address this question, we used our simulated tumour growth data for each case, noting the number of *in silico* mice at each time point. We record the time at which a mouse is ‘sacrificed’, which happens when the tumour volume reaches 2 cm^{3}, as typically done in experimental studies [40]. This approach for modelling population survival allows us to closely mimic the practice in preclinical animal studies, and provides easily interpretable insights for researchers and clinicians.

We used the simulated population survival data to determine if *k*_{0}, *k*_{1} or *k*_{0}/*k*_{1} can be used to discriminate between tumours for which anti-VEGF treatment is effective or not. We found that, in each case, a range of *k*_{0}/*k*_{1} ratios, as well as *k*_{1}, can be used to distinguish the population response to the anti-VEGF treatment (figure 3*b,c*). We term these ‘ratio_{thresh}’ and ‘*k*_{1,thresh}’ the values of the growth kinetic parameters that separate the simulated mouse population into groups with significantly different survival estimates. By contrast, we did not find any values of *k*_{0} alone that could be used to separate the simulated mouse population into groups whose survival estimates are statistically different for the Roland, Zibara and Volk2011b cases. For Tan and Volk2008, we only found one such *k*_{0} value in each case (figure 3*a*).

Interestingly, although the ranges of generated *k*_{0}/*k*_{1} ratios and *k*_{1} were different for each of the six sets of tumour growth data, we found that there is an overlap among the potential ratio_{thresh} or *k*_{1,thresh} values found in each of the six cases. The common range of ratio_{thresh} is 9.757 to 17.982, and that of *k*_{1,thresh} is 1.391 × 10^{−6} to 1.931 × 10^{−6}. This means that separating the treatment group by any *k*_{1,thresh} *or* ratio_{thresh} value within its respective range will produce two groups of treated mice that have statistically different survival estimates. Specifically, the treated group with *k*_{0}/*k*_{1} ratios larger than ratio_{thresh} has a better survival estimate than the treated group with smaller ratios. The treated group with *k*_{1} smaller than *k*_{1,thresh} has a better survival estimate than the treated group with larger *k*_{1}.

We used the median ratio_{thresh} value (13.689) to illustrate this distinction. We compare the survival estimates for a total of six groups: (i) all mice in the control group; (ii) all mice in the treatment group; (iii) control group with *k*_{0}/*k*_{1} < ratio_{thresh}; (iv) control group with *k*_{0}/*k*_{1} > ratio_{thresh}; (v) treatment group with *k*_{0}/*k*_{1} < ratio_{thresh}; and (vi) treatment group with *k*_{0}/*k*_{1} > ratio_{thresh}. We generated the Kaplan–Meier survival curves for these groups for each of the six cases investigated (figure 4). We also estimated the median survival of the six groups in each case (table 1), the Mantel–Haenszel hazard ratio (HR), with 95% confidence interval (CI), and the *p*-values from the Mantel–Cox log rank test for survival curve comparison (table 2). When comparing two groups, if the HR is less than one, the first group has a lower death rate (see Methods). Together these analyses emphasize that mice with larger *k*_{0}/*k*_{1} ratios survive for longer, with *p*-value < 0.05. Interestingly, for Zibara and Volk2011a, although the anti-VEGF treatment does not significantly reduce tumour growth and therefore does not yield a better survival estimate for the treated groups compared to their control groups (figures 2*b,e* and 4*b,e*), the stratified groups yield significantly different survival estimates. That is, the control and treated groups with *k*_{0}/*k*_{1} ratios larger than ratio_{thresh} have better survival estimates than those with smaller *k*_{0}/*k*_{1} ratios.

We performed a similar analysis using the median *k*_{1,thresh} value (1.661 × 10^{−6}) to show the distinction between the survival estimates (electronic supplementary material, figure S3). The control and treated groups with *k*_{1} smaller than *k*_{1,thresh} have better survival estimates than those with larger *k*_{1} values. We also estimated the median survival of the six groups separated using the median *k*_{1,thresh} (table 3), the Mantel–Haenszel HR and the *p*-values from the Mantel–Cox log rank test for survival curve comparison (table 4). From these analyses, mice with smaller *k*_{1} survive longer than those with larger *k*_{1}, and the HR is smaller than one (*p* < 0.05).

### 2.3. Alternative treatment strategies to improve survival estimates

We next sought to understand whether alternative treatment protocols can effectively reduce tumour volume for the Zibara and Volk2011a cases, because the baseline protocol did not significantly affect tumour volume. For the Zibara case, we simulated the original treatment protocol used in the experimental study (termed ‘protocol Z’). This protocol starts the 10 mg kg^{−1} biweekly treatment upon tumour engraftment (assuming the initial tumour volume to be 0.004 cm^{3}) [36]. The predicted tumour volumes are smaller in the treated group (electronic supplementary material, figure S4A), recapitulating the findings from the published experimental study. The predictions may suggest that, in this case, starting the treatment earlier is more effective in limiting the tumour growth. For mice with *k*_{0}/*k*_{1} ratios larger than the median ratio_{thresh}, or with *k*_{1} smaller than the median *k*_{1,thresh}, the HR between the treated and control groups is smaller than one, and the survival curves are significantly different (*p* < 0.0001) (tables 2 and 4).

For Volk2011a, we simulated treatment termed ‘protocol V11a’, which starts the 10 mg kg^{−1} biweekly treatment when the tumour volume reaches 0.5 cm^{3}, a start time extracted from the published preclinical study [39]. After 12 weeks, the simulated mean tumour volumes in the treated group are significantly smaller than the control tumours (electronic supplementary material, figure S4B). However, the survival estimates were not significantly different (*p* > 0.05). Again, the treated group with *k*_{0}/*k*_{1} ratios larger than the median ratio_{thresh}, or with *k*_{1} smaller than the median *k*_{1,thresh}*,* has a significantly better survival estimate than the opposite group (*p* < 0.0001) (tables 2 and 4). This phenomenon is similar to that observed in the Volk2011a case using protocol A, where the two groups separated according to the *k*_{0}/*k*_{1} ratio or *k*_{1} have distinct survival estimates, but there is no significant difference between the treated and control groups.

Finally, we explored whether another treatment protocol could significantly improve the survival estimates for the treated group compared to the control. We simulated protocol V11a-D, where biweekly treatment starts when the tumour volume reaches 0.5 cm^{3}, and the drug dosage is doubled to 20 mg kg^{−1}. This treatment protocol significantly limits the tumour growth (electronic supplementary material, figure S4C), and the survival curves are significantly better for the treated group compared to the control (*p* < 0.0001). Overall, the treated and control groups have an HR of 0.2016 (95% CI: 0.1343–0.3027) (table 2).

### 2.4. Validation of thresholds using an independent dataset

To validate the use of the range of ratio_{thresh} and *k*_{1,thresh} values that we found, we used a recently published independent set of data that measures tumour growth in mice with MDA-MB-231 xenografts, with or without bevacizumab treatment [41]. First, we fitted the model to the measured tumour volumes without treatment. We obtained 12 sets of estimated parameter values for *k*_{0}, *k*_{1} and *Ang*_{0} that allow the model to best fit to the control data. We then validated the fitted model by simulating anti-VEGF treatment and comparing to the experimental measurements. The predicted tumour growth with treatment matches closely to the experimental data (figure 5*a*).

Using the same approach as described above, we generated 400 sets of tumour volumes for an *in silico* mouse population with and without treatment (referred to as ‘Mollard’). To do so, we randomly varied *k*_{0} and *k*_{1} from the ranges of the 12 sets of estimated parameter values from model fitting to the Mollard dataset, with *Ang*_{0} held constant at the median of its estimated values. The simulated tumour volumes for the control and treated groups are shown in figure 5*b*.

We generated the population survival data based on the simulated tumour growth profiles. We tested whether the common range of ratio_{thresh} and *k*_{1,thresh} values identified using the six datasets described above are able to separate the population survival data for this validation case (Mollard). For all ratio_{thesh} values within the range, the survival estimate of the treated mice with *k*_{0}/*k*_{1} ratios larger than the threshold is better than those with smaller *k*_{0}/*k*_{1} ratios. Examples using the median ratio_{thresh} and the median *k*_{1,thresh} are shown in figure 5*c,d*. We calculated the HR values, as well as the *p*-value from the Mantel–Cox log rank test among the treated and control groups, separated using the median of the common ratio_{thresh} range (table 2) or the common *k*_{1,thresh} range (table 4). Thus, we were able to validate the threshold values.

### 2.5. Tumour growth dynamics among stratified populations

We explored the dynamics of the tumour growth for the groups separated by the threshold values to better understand why the anti-VEGF treatment has differential effects in the simulated mouse populations. As researchers have pointed out, log transformation of tumour growth data provides information on tumour growth rates (given by the slope of the curve) and is more suitable for detecting a transient biological or therapeutic effect [40,42,43]. Therefore, we compared the mean RTV time courses (electronic supplementary material, figure S1) and the mean tumour volume data plotted on the log scale (electronic supplementary material, figure S5) of the groups stratified by the median ratio_{thresh} (13.869) in each case.

For Roland, Tan and Volk2008, the mean RTV of the group with larger *k*_{0}/*k*_{1} ratios (electronic supplementary material, figure S1A,C,D) is initially larger, and then becomes smaller relative to the opposite group. This switch occurs because in the group with larger *k*_{0}/*k*_{1} ratios, the difference between the treated and control tumour volumes is smaller at early times, and then becomes larger (electronic supplementary material, figure S5). Meanwhile, the actual tumour volumes for this group are both relatively low. As a result, this group survives longer (figure 4). For the Mollard case, the differences between the treated and control tumour volumes in the group with larger *k*_{0}/*k*_{1} ratios are larger (figure 6*b*, dotted curves), giving rise to the larger mean RTV (figure 6*a*). However, the group with larger *k*_{0}/*k*_{1} ratios still survives longer because the actual tumour volumes are relatively low (figure 5*c*).

The tumour volume data plotted on the log scale also reveal that the tumour growth rates of the control and treated groups diverge at different time points. For Roland, the gap between the control and treated groups continually increases when plotted on a linear scale (figure 7*a*). However, the tumour volumes plotted on the log scale show that their growth rates are mostly different during days 14–40. The growth rates become similar during the later stage (after 40 days), as evidenced by the parallel curves on the log scale (figure 7*b*). Therefore, the increasingly large gap between the tumour volumes is a result of early differences in the tumour growth rates. A similar phenomenon is observed for Volk2011b, where the tumour growth rate of the treated group is suppressed transiently at early times but not in the later stage (electronic supplementary material, figure S5). In Zibara, Tan and Volk2008, the growth rates become different between day 30 and day 45, and only gradually become similar towards the end of the simulated time. Overall, analysis of the growth curves plotted on the log scale reveals that the anti-VEGF treatment has differential effects in limiting tumour growth, and the effects occur at different stages for the simulated cases. The treatment effect is predicted to be stronger for the group with *k*_{0}/*k*_{1} ratios larger than the median ratio_{thresh}.

## 3. Discussion

In this study, we focus on identifying tumour growth kinetic parameters as potential biomarkers for the outcome of anti-VEGF treatment. We developed a computational approach that incorporates model training, simulation of tumour growth within a heterogeneous population, and estimation and analysis of population response.

We simulated anti-VEGF treatment and compared the effect of treatment across tumour-bearing mice generated from our previous fitting to six independent preclinical studies. For most simulated tumours, the anti-VEGF agent significantly reduces tumour volume compared to the control. However, our simulations for Zibara and Volk2011a show that these populations do not respond to the treatment (figure 2*b,e*), which is different from the effect seen experimentally. This difference occurs for two reasons. First, our simulated treatment protocol A is universal across the six cases, and is different from what was used in each of the original six experimental studies. Second, in our simulations, *k*_{0} and *k*_{1} are varied simultaneously and independently of each other, possibly resulting in more variability than what occurs experimentally.

Our study demonstrates that the *k*_{0}/*k*_{1} ratio or *k*_{1} alone can be used to stratify the population response with or without anti-VEGF treatment. This finding agrees with our previous finding through PLSR analysis that the ratio is a key predictor of the tumour response to anti-VEGF treatment [26]. Building on that framework, we found that the survival estimate of mice with larger *k*_{0}/*k*_{1} ratios or smaller *k*_{1} is better compared to those with smaller ratios or higher *k*_{1}. Interestingly, the result for the ratio is the opposite of the conclusion we drew previously (that a larger ratio correlates with a poorer response to treatment). However, in that work, we focused only on whether the final RTV value was low. This highlights the fact that only evaluating the endpoint RTV of the treated and control group and neglecting the actual tumour volume data over time can lead to misinterpretation of the treatment effect. Indeed, researchers have recognized that while most preclinical studies focus on the endpoints of tumour growth, monitoring tumour growth kinetically may be more insightful [42,43].

We found that, in two cases (Volk2011a simulated with protocols A and V11a), no significant difference is observed in the survival estimates between the treated and control groups. However, even for these cases, two populations with significantly different survival estimates can be identified based on their *k*_{0}/*k*_{1} ratios (figure 4*b,e*) or *k*_{1} value (electronic supplementary material, figure S3B,E). This indicates that even when the treatment is not effective in reducing tumour volume, there is still a difference in tumour growth dynamics between the two populations stratified based on the tumour's growth kinetic parameters. Thus, we believe that the *k*_{0}/*k*_{1} ratio or *k*_{1} may be prognostic biomarkers to stratify populations for their survival estimate without anti-angiogenic treatment. Interestingly, the parameters provide mechanistic insight into tumour growth. In particular, they highlight that slower linear growth (larger ratio or smaller *k*_{1}) results in less aggressive overall tumour growth (electronic supplementary material, figure S5) and, therefore, a better survival outcome.

Another interesting aspect is the utility of *k*_{1} to serve as a prognostic biomarker. Although *k*_{1} was not revealed as a strong predictor of the final RTV previously in the PLSR analysis, it is inversely correlated with the *k*_{0}/*k*_{1} ratio, and therefore, in our study, it also can be used to stratify the population survival outcome. Here, performing the survival analysis addresses one limitation from our previous PLSR analysis, where we were able to identify which parameters were related to treatment efficacy, but could not identify the specific relationship between the kinetic parameter values and effectiveness of the treatment.

Compared to the mean RTV data, the tumour volume data provide more useful insight into the tumour growth characteristics of the stratified population. In particular, the tumour volume plotted on the log scale more clearly illustrates the source of the differences in the population survival estimates. Specifically, we found that larger *k*_{0}/*k*_{1} ratios often yield slower tumour growth in a population, and therefore, lead to a better survival estimate of the population. This conclusion could not be made if we were to only analyse the RTV data. In addition, the tumour volumes plotted on the log scale reveal that the effect of anti-VEGF treatment in tumour growth can be relatively transient or gradual.

Our study uses a predictive computational model of tumour growth. This is a pharmacokinetics–pharmacodynamics model with mechanistic detail that goes beyond what is found in other models. However, in the future, this model can be expanded to address limitations that are not currently accounted for. For example, we do not account for changes in tumour vascularity relative to tumour volume. We assume the vascular volume relative to total tumour volume remains constant, given the lack of robust quantitative data needed to develop a mathematical function describing how tumour vascularization changes over time. In addition, vascular normalization is an important process that has been shown to affect tumour growth and can be regulated by anti-VEGF agents [32]; however, this process is not included in our model. These aspects can be implemented into the model as more quantitative data become available and enable us to characterize the dynamics of vessel normalization. The model can then be further extended to account for other characteristics of tumour progression, including tumour perfusion and metastatic potential. The model can also be adapted to simulate the effect of cytotoxic drugs that target tumour cells, which in turn will affect the tumour volume. Furthermore, the range of threshold values for tumour stratification is constrained by the estimated parameter values from model training to each experimental dataset. It is possible that artefacts from experimental data quantification led to bias in the range of the fitted parameter values. This can be improved when more quantitative experimental data become available for additional model training. We note that the biomarker candidates identified in this study are best used to stratify populations for their survival outcome, whether the mice receive treatment or not, rather than to predict treatment efficacy. This is primarily because the datasets used for model training were tumour volumes measured over several weeks. Our results would be of broader applicability if only pretreatment data were adequate to train the model. We attempted such an approach in previous work [26]; however, the simulated volumes varied widely, preventing us from making conclusive predictions. Despite this perceived limitation, our modelling approach generates hypotheses about potential biomarkers, and spurs on experimental validation to ensure the utility of the biomarkers identified.

Our study demonstrates a time- and cost-effective way to generate large *in silico* mouse populations, predict anti-VEGF treatment outcome and stratify the populations. This approach provides useful information that could facilitate efficient experimental design, such as predicting the effect of different treatment protocols (varying the dosage and the timing of the injections). Additionally, our modelling approach can be adapted for analysis of the patient treatment outcome in clinical studies. With data from a small patient population, we can develop a patient-specific model and generate a larger *in silico* population. Analysis of the simulated tumour growth and survival data can be used to identify biomarkers that predict responders versus non-responders to anti-VEGF treatment, stratify the predicted population survival and test the response to various treatment schedules.

## 4. Conclusion

We examined tumour growth kinetic parameters as potential biomarkers of anti-angiogenic treatment outcome. Using a computational model that simulates VEGF-dependent tumour growth in tumour-bearing mice, we generated an *in silico* mouse population and related the kinetic parameters that characterize tumour growth to the response to anti-VEGF treatment. We found that the ratio between two tumour growth kinetic parameters, *k*_{0} and *k*_{1}, as well as *k*_{1} alone, can be prognostic biomarkers and that the simulated treatment protocol may have a better outcome for mice whose tumours have smaller linear growth rates. In fact, we found ranges of threshold values for the *k*_{0}/*k*_{1} ratio and *k*_{1} that distinguish tumours' response to the anti-VEGF treatment. This study demonstrates an approach for identifying tumour growth kinetic parameters as potential biomarkers, and this model framework can be adapted to predict the efficacy of other anti-angiogenic strategies.

## 5. Methods

### 5.1. Computational model

We use our previously calibrated and validated model of VEGF binding and distribution in a tumour-bearing mouse [26]. Electronic supplementary material, file S1, contains the full model description. Briefly, the model comprises three compartments representing the whole mouse (figure 1): normal tissue, blood and tumour tissue. We include human VEGF isoforms (VEGF_{121} and VEGF_{165}) secreted by tumour cells, as well as mouse isoforms (VEGF_{120} and VEGF_{164}) secreted by endothelial cells and muscle fibres. The model includes cell surface VEGF receptors, VEGFR1 and VEGFR2 and soluble VEGFR1 (sVEGFR1). We include neuropillin co-receptors (NRP1 and NRP2) that bind VEGF directly and also form tertiary complexes with the VEGFRs. The protease inhibitor α-2-macroglobulin binds VEGF in blood plasma. We consider the luminal and abluminal endothelial surfaces at the interface between the blood and each tissue compartment. The VEGF isoforms and sVEGFR1 are transported between compartments via transendothelial macromolecular permeability and lymphatic flow. Additionally, species are removed via clearance.

VEGF binding to its receptors on endothelial cells promotes intracellular signalling that mediates angiogenesis. Thus, we explicitly account for VEGF-mediated tumour growth by incorporating the concentration of ligated receptors localized on tumour endothelial cells into the tumour volume equation (figure 1). We simulate anti-VEGF treatment as intravenous injections lasting for 1 min by adding a net rate of secretion of the drug (bevacizumab) directly into the blood compartment.

### 5.2. Numerical implementation

Model equations were implemented in Matlab using the SimBiology toolbox. The model is provided as the SimBiology project file, SBML, Matlab m-file and full list of equations (electronic supplementary material, file S2). Parameter fitting was performed using the *lsqnonlin* function in Matlab. Kaplan–Meier survival estimation was performed using the *kmplot* function in Matlab, and GraphPad Prism was used for statistical survival analyses.

### 5.3. Simulation of *in silico* mouse population

We previously fitted the model to six independent control datasets to estimate the growth kinetic parameters (*k*_{0}, *k*_{1} and *Ang*_{0}). The parameter *ψ* was held constant, as it was not shown to significantly influence tumour growth, compared to the other parameters. The model-predicted tumour growth curves match closely to the experimental data (fitting error range: 0.0405–0.1833).

Here, we generated 400 sets of values for *k*_{0} and *k*_{1}, randomly selected from a uniform distribution within the range of the best fit parameter sets from our previous study (electronic supplementary material, table S1). The *Ang*_{0} value is set to be the median of the best fits in each case. These sets were used to calculate tumour growth with or without anti-VEGF treatment, simulating a population of mice for each dataset. To keep tumour growth profiles realistic, tumours that do not reach 0.1 cm^{3} within 10 days upon tumour engraftment (assuming an initial tumour volume of 0.004 cm^{3}) were excluded from the analyses.

We simulated anti-VEGF treatment for each dataset. Treatment protocol A is simulated universally across the six cases. In this protocol, weekly treatment starts when the tumour volume reached 0.1 cm^{3}, as the switch where angiogenesis is more strongly promoted occurs when the tumour reaches 1–2 mm in diameter. The treatment dosage is 10 mg kg^{−1}. The model was simulated for 12 weeks after treatment started. We also simulated alternative treatment protocols: Z denotes biweekly treatment at a dosage of 10 mg kg^{−1} starting when the tumour volume is 0.004 cm^{3}; V11a denotes biweekly treatment (twice a week) at a dosage of 10 mg kg^{−1}, starting when the tumour volume is 0.5 cm^{3}; and V11a-D denotes biweekly treatment at a dosage of 20 mg kg^{−1}, starting when the tumour volume is 0.5 cm^{3}. Information for all treatment protocols is given in electronic supplementary material, table S3.

### 5.4. Relative tumour volume

Based on the model-generated tumour growth data, the RTV is calculated at any simulated time point as follows: An RTV value less than one indicates that the treated tumour volume is smaller than the control.

### 5.5. Kaplan–Meier survival estimation

We applied time-to-event analysis to determine the survival of each mouse population [34]. An *in silico* mouse is recorded as ‘sacrificed’ when its tumour reaches 2 cm^{3} within the simulated time. Alternatively, a mouse is recorded as ‘censored’ at a particular time point, *t*, if its tumour volume simulation remains below 2 cm^{3} but ends before time *t*. All other mice are retained in the study and recorded as ‘alive’. Survival curves were estimated by the Kaplan–Meier method using the *kmplot* function in Matlab [44], and compared using the Mantel–Cox log rank test and the Mantel–Haenszel HR in GraphPad Prism.

The HR compares the rate of death in two groups, with the assumption that the population HR is consistent over time. It is calculated using the Mantel–Haneszel approach, which is more accurate than the log rank approach [45]. As an example, an HR of 0.5 between two groups means that the death rate of the first group is half of that of the second group.

### 5.6. Determination of threshold values

To determine threshold values for the *k*_{0}/*k*_{1} ratio, we ordered the simulated mouse tumour volume data for each of the six populations according to the *k*_{0}/*k*_{1} ratio. Then, we systematically tested each *k*_{0}/*k*_{1} ratio (called ‘ratio_{thresh}’) value to see if there is a significant difference between the survival estimates for the mice with the *k*_{0}/*k*_{1} ratio above and below ‘ratio_{thresh}’ in the log rank test (*p* < 0.05). We performed a similar analysis for *k*_{0} and *k*_{1} individually to determine any *k*_{0,thresh} and *k*_{1,thresh} values.

### 5.7. Validation of the predicted biomarker

Upon identifying a potential predictive biomarker for the efficacy of anti-VEGF treatment, we validated our findings using an independent set of data that was not used to determine the range of the threshold values. To do so, we fitted the control tumour growth for the independent dataset and generated an *in silico* mouse population based on the fitted parameters.

#### 5.7.1. Data extraction

For threshold validation, data from the published *in vivo* experimental study of MDA-MB-231 xenograft tumour growth in mice by Mollard *et al*. were used for parameter estimation and validation [41]. Experimental data were extracted using the WebPlotDigitizer program [46] and are shown in electronic supplementary material, table S1.

#### 5.7.2. Parameter estimation

We trained the model to fit the control tumour growth dataset from [41] using the same approach as described in our previous work [26]. The values of tumour growth parameters *k*_{0}, *k*_{1} and *Ang*_{0} were estimated. In their study, Mollard and co-workers only reported the tumour volumes relative to day 8. However, the absolute tumour volumes are needed to determine how the tumour interstitial volume varies as a function of the total tumour volume. Therefore, we compared the RTV at each time point in the work by Mollard and co-workers to that of all the available control datasets (electronic supplementary material, figure S6). We then chose to use the interstitial volume equation from the Zibara data, given that the RTV closely matches that of the data in Mollard. Finally, we fitted our tumour growth model to the Mollard control dataset.

Fitting was performed using the *lsqnonlin* function in Matlab to minimize the sum of squared residuals (SSR):
where *V*_{exp,I} is the *i*th experimental data point of tumour volume, *V*_{sim,I} is the *i*th simulated volume at the corresponding time point and *n* is the total number of experimental data points. The minimization is subject to *Θ*, the set of upper and lower bounds on each of the free parameters.

The bounds for each parameter spanned at least two orders of magnitude: 10^{−8} to 10^{−2} for *k*_{0} and *k*_{1} and 10^{−16} to 10^{−14} for *Ang*_{0}. After fitting to the control data, we validated the estimated parameters with data not used in the fitting for model validation. Specifically, we applied the fitted model to simulate anti-angiogenic treatment (bevacizumab) and compared to the experimental measurements for the treatment case. We simulated the dosing regimen used by Mollard *et al*.: three cycles of weekly intravenous injections lasting for 1 min starting from day 5. We used the combined SSR for the RTV between model prediction and the experimental data (both control and treatment) to identify the optimal parameters. Twelve parameter sets with the smallest errors were taken to be the ‘best’ sets (electronic supplementary material, table S1) and the ranges of the estimated parameter values were used for subsequent model simulations (electronic supplementary material, table S2).

We extracted the absolute tumour volume at day 8 from previously reported data from Mollard and co-workers [47] to determine the survival estimates for a mouse population simulated based on the fitted growth kinetics parameter values.

## Data accessibility

All data analysed during this study, including the model used for data generation, are included in this published article and its electronic supplementary material information files.

## Authors' contributions

S.D.F. designed the research. A.D.A. wrote Matlab scripts and performed initial model simulations. Q.W. performed model simulations and data analysis. Q.W. and S.D.F. wrote the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Funding

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

## Acknowledgements

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

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4188212.

- Received April 10, 2018.
- Accepted July 27, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.