How to schedule VEGF and PD-1 inhibitors in combination cancer therapy?

Background One of the questions in the design of cancer clinical trials with combination of two drugs is in which order to administer the drugs. This is an important question, especially in the case where one agent may interfere with the effectiveness of the other agent. Results In the present paper we develop a mathematical model to address this scheduling question in a specific case where one of the drugs is anti-VEGF, which is known to affect the perfusion of other drugs. As a second drug we take anti-PD-1. Both drugs are known to increase the activation of anticancer T cells. Our simulations show that in the case where anti-VEGF reduces the perfusion, a non-overlapping schedule is significantly more effective than a simultaneous injection of the two drugs, and it is somewhat more beneficial to inject anti-PD-1 first. Conclusion The method and results of the paper can be extended to other combinations, and they could play an important role in the design of clinical trials with combination therapy, where scheduling strategies may significantly affect the outcome.


Background
Anti-vascular endothelial growth factor (anti-VEGF) is a drug commonly used as anticancer agent, although numerous studies show only modest results [1]. In combination with chemotherapy anti-VEGF improves anticancer therapy, although the outcome depends on the cancer type and on the scheduling of the treatment [2]. The scheduling issue arises from the fact that anti-VEGF decreases perfusion of chemotoxic agents in some cancers, including melanoma [3], breast cancer [4,5] and ovarian cancer [6], while it increases perfusion of chemotoxic agents in other cancers, such as colon cancer [7,8] and head and neck cancer [9]. More recently Asrid et al [10] reported a rapid decrease in the delivery of chemotherapy to the tumor in patients of non-small cell lung cancer (NSCLC) after anti-VEGF therapy, highlighting the importance of drug scheduling in combination therapy when anti-VEGF is one of the drugs. Less than 4% of positive phase II cancer clinical trials with combination chemotherapy demonstrate improvement of care in phase III [11]. Hence, the decision to go from phase II to phase III needs to identify more effectively which combinations will have a higher probability of success in phase III [12]. It was suggested in [13] that the design of clinical trials with combination therapy should be based, among other factors, on the scientific rationale underlying data and hypothesis for the combination.
In a previous work [14] we considered a combination therapy with a checkpoint inhibitor and cancer vaccine, and explored the synergy between the two drugs, taking into account potential negative side effects. In another paper [15] we considered the combination of BRAF/MEK inhibitor and checkpoint inhibitor, and showed that although the two drugs are positively correlated for most combinations of the doses, there is an exceptional range of doses where the two drugs are mutually antagonistic.
In the present paper we consider a combination therapy of anti-VEGF and a checkpoint inhibitor, and focus on the scheduling issue of these drugs. The rationale for using such a combination originates from the fact that VEGF impairs the function of anti-cancer T cells [16][17][18][19][20]; hence VEGF inhibition will enhance T cells function, and checkpoint blockade could therefore significantly advance antitumor therapy. The anticancer synergy between anti-VEGF and checkpoint inhibitors is currently being evaluated in clinical trials in renal cancer [21,22].
In the case where anti-VEGF therapy decreases the perfusion of a second antitumor agent, the following question arises [23,24]: Should the treatment with combination of anti-VEGF and a checkpoint inhibitor be given at the same time, or is it more beneficial to delay treatment of one of the two drugs, so that they are given non-overlappingly?
To address this question we develop a mathematical model using a system of partial differential equations (PDEs). The variables of the model include CD4 + (Th1) and CD8 + T cells, regulatory T cells (Tregs), dendritic cells (DCs), endothelial cells, and cancer cells. The model also includes VEGF and TGF-β produced by cancer cells, and cytokines IL-12 and IL-2. The network of interactions among these species is shown in Fig. 1. This figure includes also oxygen concentration, and programmed cell death protein 1 (PD-1) and its ligand PD-L1. As indicated in Fig. 1. VEGF impairs the maturation of (antigenpresenting) dendritic cells [25,26], and it suppresses the functions of activated T cells [16][17][18][19][20]; VEGF also enhances the expression of PD-1 on CD8 + T cells [27], and induces Treg proliferation [28].
The mathematical model is based on the network shown in Fig. 1. The basic assumption in the model is that anti-VEGF decreases the perfusion of anti-PD-1 in the tumor microenvironment as in the case in melanoma, breast cancer, and ovarian cancer. The range of the injected amounts of drugs in the model was chosen so that the simulations will be in agreement with experimental results in mice models [27,29].
We can use the model to assess the efficacy of the combination therapy when the treatment with anti-VEGF is delayed or advanced relative to anti-PD-1, for different schedules if treatment. In particular we show that the treatment is significantly more effective if instead of administering the two drugs at the same time, we administer them non-overlappingly.
We finally consider briefly the case where anti-VEGF increases the perfusion of anti-PD-1, and show that in this case there are more benefits when the two drugs are given simultaneously.

Mathematical model
The mathematical model is based on the network shown in Fig. 1. The list of variables is given in Table 1.
We assume that the total density of cells within the tumor remains constant throughout the tumor tissue, for all time: and that the density of debris of dead cells is also constant. We further assume that the densities of immature dendritic cells, and of naive CD4 + and CD8 + T cells, remain constant throughout the tumor tissue. As cancer cells proliferate, they "push away, " or displace, other cells. There is anti-VEGF concentration also migration of endothelial cells and immune cells into the tumor. Since the total density of cells was assumed to be constant at each point and time, by Eq.(1), these increases in the population of cells create a pressure (p) among the cells with an associated velocity field u. Under some additional assumptions on the material structure of the tumor, one can actually connect u to p (for example, by Darcy's law in porous media), but we shall not need to do this in our model. The vector u is a function of space and time, taken in units of cm/day. We also assume that all the cytokines and anti-tumor drugs are diffusing within the tumor tissue, and that also the cells are undergoing diffusion (i.e. dispersion), although with much smaller coefficients.
Although in our model we use densities of cells, it is interestingly to visualize how individual cells interact within the tumor. Figure 2 displays a distribution of cells in space, based on Fig. 1. We note, in particular, that cancer cells move toward the tumor boundary where the oxygen level can supports their abnormal proliferation; hence, by Eq. (1), the other types of cells are "pushed" toward the tumor core.

Equation for DCs (D)
By necrotic cancer cells (N C ) we mean cancer cells undergoing the process of necrosis. Necrotic cancer cells release high mobility group box 1 protein (HMGB-1) [30]. We model the dynamics of N C and HMGB-1 (H) by the following equations: where λ N C C is the rate of cancer cells becoming necrotic and λ HN C is the production rate of HMGB-1 by necrotic cells. We note that since molecules like HMGB-1, or other proteins, are several orders of magnitude smaller than cells, their diffusion coefficients are several orders of magnitude larger than the diffusion coefficients of cells, and they are only marginally influenced by the cells velocity u, so we do not include a velocity term in their equations. The degradation of HMGB-1 is fast (∼0.01/day) [31], and we assume that the process of necrosis is also fast. We may then approximate the two dynamical equations by the steady state equations λ N C C C − d N C N C = 0 and λ HN C N C − d H H = 0, so that H is proportional to C. Dendritic cells are activated by HMGB-1 [32,33]. Hence, the activation rate of immature dendritic cells, with density D 0 , is proportional to D 0 H K H +H , or to D 0 C K C +C . Here, the Michaelis-Menten law was used to account for the limited rate of receptor recycling which takes place in the process of DCs activation. The dynamics of DCs is given by the following equation: where δ D is the diffusion coefficient, d D is the death rate of DCs, and 1/(1 + G/K DG ) represents the impairment of maturation of dendritic cells by VEGF [25,26].

Equation for Tregs (T r )
Naive CD4 + T cells differentiate into T r cells under the influence of TGF-β: (T β ) [37,40] and VEGF [28]. The equation for T r takes the following form:

Equation for endothelial cells (E)
Endothelial cells are chemoattracted by VEGF, and their proliferation is increased by VEGF [41,42]. The equation for the density of endothelial cells is given by ∂E ∂t where E M is the carrying capacity of endothelial cells, λ E (G) = λ EG (G−G 0 ) + , and G 0 is a threshold below which endothelial cells do not proliferate [43]. Here we used the notion:

Equation for cancer cells (C)
We assume a logistic growth for cancer cells, with carrying capacity (C M ), in order to account for competition for space among these cells. The proliferation rate depends on the density of oxygen (W ) [42]. The equation for C takes the following form: where η 1 and η 8 are the killing rates of cancer cells by T 1 and T 8 cells, respectively. d C is the natural death rate of cancer cells, and where W 0 is the normal level of oxygen concentration in the blood.

Equation for oxygen (W)
Oxygen is infused through the blood [41,42]. We identify the blood concentration with the density of endothelial cells and, accordingly, write the equation for W in the following form: consumption by cells (11) where d W W represents the take-up rate of oxygen by all the cells.

Equation for VEGF (G)
VEGF is produced by cancer cells [41,42] and is depleted by anti-VEGF. Hence G satisfies the following equation: where B is the effective anti-VEGF concentration in the tumor, and Here we assumed that the secretion rate of VEGF by cancer cells increases with the oxygen level, but falls off when oxygen level exceeds a certain level, W * [44].

Equations for PD-1 (P), PD-L1 (L) and PD-1-PD-L1 (Q)
PD-1 is expressed on the surface of activated CD4 + T cells, activated CD8 + T cells, and Tregs [38,45]. We assume that the number of PD-1 receptors per cell is the same for T 1 and T 8 cells, but is smaller, by a factor ε T , for T r cells. VEGF increases the PD-1 on T 8 cells by a factor ε G G [27]. If we denote by ρ P the ratio between the mass of one PD-1 protein to the mass of one T cell, then the total concentration of PD-1 on T 1 and T r cells is given by and the concentration of PD-1 on T 8 cells is given by PD-L1 is expressed on activated CD4 + T cells, activated CD8 + T cells [38], Tregs [46], and cancer cells [38,39]. We assume that the number of PD-L1 per cell is the same for T 1 and T 8 cells, and denote the ratio between the mass of one PD-L1 protein to the mass of one cell by ρ L . Then for some parameters ε T , ε C , where ε C depends on the specific type of tumor.
To a change in T = T 1 + ε T T r , given by ∂T ∂t , there corresponds a change of P 1 , given by ρ P ∂T ∂t . For the same reason, ∇ · (uP 1 ) = ρ P ∇ · (uT) and ∇ 2 P 1 = ρ P ∇ 2 T when no anti-PD-1 drug is injected. Hence, P 1 satisfies the following equation: Recalling Eqs. (3) and (5) for T 1 and T r , we get, where RHS=right-hand side. Similarly, When only anti-PD-1 drug (A) is injected, PD-1 is depleted at a rate μ PA P 1 A, but when also anti-VEGF (B) is injected the depletion of PD-1 is reduced, due to restricted perfusion, by a factor 1 1+B/K PB ; since all other species of the model evolve within the tumor, we ignored the effect of restricted perfusion in setting up their dynamics. We conclude that P 1 satisfies the following equation: blockade of perfusion by anti-VEGF (15) Similarly, blockade of perfusion by anti-VEGF (16) PD-L1 combines with PD-1 on the plasma membrane of T cells, to form a complex PD-1-PD-L1 (Q) [38,39]. Denoting the association and disassociation rates of Q by α PL and d Q , respectively, we can write where P = P 1 + P 8 . The half-life of Q is less then 1 s (1.16 × 10 −5 day) [47], so that d Q is very large. Hence we may approximate the dynamical equation for Q by the steady state equation, α PL PL = d Q Q, so that where σ = α PL /d Q .

Equation for anti-PD-1 (A)
We assume that the drug enters the tumor from the boundary and quickly diffuses, so that its concentration is a constant γ A throughout the tumor during the dosing period. We denote by μ AP the depletion rate of A caused by blocking PD-1. Hence, where I A (t) = 1 during dosing, and I A (t) = 0 otherwise.

Equation for anti-VEGF (B)
We denote by γ B the concentration of the injected drug B during the dosing period, and by μ BG the depletion rate of B blocking VEGF. The equation for B is then given by where I B (t) = 1 during dosing, and I B (t) = 0 otherwise.

Equation for cells velocity (u)
We assume that the density of each cell type in the growing tumor tends to a steady state, and take the density of the extracellular matrix (ECM) in steady state to be 0.6g/cm 3 . We take the steady state density of endothelial cells to be E = 2.5 × 10 −3 g/cm 3 [43]. The steady state densities of the immune cells D, T 1 , T 8 , T r and C (in units of g/cm 3 ) are taken to be (see Appendix: Parameter estimation) With these choices, the constant in Eq. (1) equal to 0.4064. We further assume that all cells have approximately the same diffusion coefficient. Adding Eqs. (2)-(7), we get To simplify the computations, we assume that the tumor is spherical and denote its radius by r = R(t). We also assume that all the densities and concentrations are radially symmetric, that is, they are functions of (r, t), where 0 ≤ r ≤ R(t). In particular, u = u(r, t)e r , where e r is the unit radial vector.

Equation for free boundary (R)
We assume that the free boundary r = R(t) moves with the velocity of the cells, so that

Boundary conditions
We assume that the naive CD4 + T cells and CD8 + T cells which migrated from the lymph nodes into the tumor microenvironment have constant densitiesT 1 andT 8 at the tumor's boundary, and, that upon crossing the tumor boundary, T 1 and T 8 are activated by IL-12 and T r is activated by T β . We then have the following flux conditions at the tumor's boundary: where we take σ T (I 12 ) = σ 0 We assume that the endothelial cells that are attracted by the VEGF into the tumor microenvironment have constant densityÊ at the tumor's boundary, and take The boundary condition for the oxygen is given by where W 0 is the normal level of oxygen concentration in the blood. We impose no-flux boundary condition for all the remaining variables: No-flux for D, C, I 12 , I 2 , T β , G, P 1 , P 8 , A, and B at r = R(t).
It is tacitly assumed here that the receptors PD-1 become active only after the T cells are already inside the tumor.

Initial conditions
We take the following initial values (in units of g/cm 3 ): Note that the total density of cells at time t = 0 satisfies Eq. (1) with the already chosen constant 0.4064. We also mention that the choice of initial conditions does not affect the simulation results after a few days.

Results
The simulations of the model were performed by Matlab, based on the moving mesh method for solving partial differential equations with free boundary [48] (see the section on computational method).
Simulation results for mouse model Figure 3 shows the profiles of the average densities/concentrations of all the variables of the model in the control case (no drugs) and the tumor volume in the first 30 days. We note that each species X reaches a nearly steady state that is approximately the halfsaturation value, K X , as assumed in the estimation of some of the model parameters (in Appendix). This shows the consistency in the parameters estimation.  Tables 3 and 4, for a mouse model From Fig. 3 we see that in the first few days E, and corresponding W, are increasing before they begin to gradually decrease and reach a steady state. The production rate of G is given by λ(W )C where λ(W ) is bimodal and the oscillation in the profile of G is affected by the bimodal profile of W and the initial growth in the profile of C. The initial increase in G results in initial decrease in D, and hence also a decrease in IL-12 and the T cells. The profile of C is affected by the growth in W and by the decrease in T 1 and T 8 : C begins to grow after a short time, but the growth rate decreases to near 0 as time increases. Before applying our model to clinical situations, we need to determine the order of magnitude of γ B and γ A . The parameter γ B is proportional to the amount of anti-VEGF administered to a patient, but it is difficult to determine this proportionality coefficient; the same is true for γ A . We therefore conducted simulations with different choices of γ B , γ A in order to find values for which the simulations agree qualitatively with mice experiments. One set of simulations is displayed in Fig. 4. Figure 4a shows that with γ B = 3 × 10 −8 g/cm 3 · day and γ A = 3 × 10 −8 g/cm 3 · day the anti-PD-1 as a single agent reduces the volume of tumor more than anti-VEGF as a single agent, in agreement with Fig. 1 in [29]. Figure 4b shows that with γ B = 3.5 × 10 −8 g/cm 3 · day and γ A = 0.5 × 10 −8 g/cm 3 · day the anti-VEGF as a single agent reduces the tumor volume more than the anti-PD-1 as single agent, in agreement with Fig. 4j in [27]. In Fig. 5 we see that (with γ A = 0) anti-VEGF with γ B = 2×10 −8 g/cm 3 ·day reduces significantly the PD-1 expression on CD8 + T cells, in agreement with experimental results in [27].

Clinical trials in silico
In clinical trials the treatment and follow-up periods are much longer than in mouse models, and the tumor growth is significantly slower. Accordingly, we shall modify some of the parameters in order to slow the growth of the tumor; these parameters are λ DC , λ E , λ CW . We shall also decrease the range of the drug, from order of magnitude 10 −8 to 10 −9 − 10 −10 .
We consider clinical treatment in cycles of 9 weeks, whereby each drug is given continuously for 3 weeks during a cycle with drug holiday for 6 weeks. We introduce three different schedules. In schedule S1 anti-PD-1 and anti-VEGF are both administered continuously in the first 3 weeks of the 9-week cycle. In schedule S2 anti-PD-1 is given continuously in the first 3 weeks, followed by anti-VEGF in the next 3 weeks, with no drugs for the remaining 3 weeks of the cycle. Schedule S3 is similar to schedule S2 with anti-VEGF in the first 3 weeks and anti-PD-1 in the next 3 weeks. We refer to schedules S2 and S3 as non-overlapping schedules. Figure 6 shows the profile of the tumor volume under schedules S1, S2 and S3 for four different dose pairs (γ B , γ A ). Table 2 summarizes the time, in weeks, at which the tumor volume decreased by 95% from its initial size. We see that the non-overlapping schedules S2 and S3 reduce the tumor volume significantly faster than schedule S1, and schedule S2 is somewhat more effective than schedule S3.  (γ B , γ A ). The anti-VEGF or/and the anti-PD-1 treatment started at day 0 and continued for 10 days. a γ B = 3 × 10 −8 g/cm 3 · day, γ A = 3 × 10 −8 g/cm 3 · day; b γ B = 3.5 × 10 −8 g/cm 3 · day, γ A = 0.5 × 10 −8 g/cm 3 · day. All other parameters are same as in Fig. 3 a b  (1)), the immune cells and the endothelial cells are "pushed" back toward the core of the tumor, so that their profile is monotone decreasing from the tumor core toward the tumor boundary. Figure 8 shows "efficacy maps, " namely, the time at which tumor volume was reduced by 95%, under schedules S1, S2 and S3, for a range of the parameters γ B and γ A . The horizontal axis scales the dose amount of anti-VEGF, γ B , in units of g/cm 3 · day, and the vertical axis scales the dose amount of anti-PD-1, γ A , in unit of g/cm 3 · day. The color columns show the time at which the tumor volume decreased by 95% from its initial size. We see that as γ A or γ B increases, the time when the tumor volume was reduced by 95% is decreased. Also, as in the special case of Fig. 6, schedules S2 and S3 reduce significantly the time for the 95% reduction, and the tumor volume reduction by schedule S2 is a little faster than by schedule S3. a b c Fig. 6 Tumor volume under schedules S1, S2 and S3 for 4 pairs (γ B , γ A ). Here γ B1 = 9.5 × 10 −9 g/cm 3 · day, γ A1 = 1.2 × 10 −10 g/cm 3 · day, γ B2 = 10 × 10 −9 g/cm 3 · day, γ A2 = 1.5 × 10 −10 g/cm 3 · day. a Tumor volume under schedule S1; b Tumor volume under schedule S2; c Tumor volume under schedule S3 Table 2 The time, in weeks, at which the tumor volume decreased to 5% of its initial size We conclude that non-overlapping treatment is far more beneficial than simultaneous treatment, and there is some advantage in injecting anti-PD-1 before anti-VEGF.
So far we considered the case where anti-VEGF blocks perfusion. We next consider briefly the case where anti-VEGF increases perfusion (as in colon cancer [7,8]). In this case we have to modify Eqs. (13) and (14) by replacing the term 1/(1 + B/K PB ) by 1 + B K B +B . In contrast to Figs. 6, 9 shows that simultaneous injection reduces the tumor volume (by 95%) faster than non-overlapping injections: S1 is preferable to S2 and S3, while S2 is preferable to S3. These conclusions hold also when γ B and γ A vary in the same range as in Fig. 8.

Sensitivity analysis
We performed sensitivity analysis with respect to the tumor volume at day 30 for parameters λ DC , λ T 8 I 12 , λ T r G , λ E , λ GW , η 8 , K DG , K TQ , K TG , K PB . Following the method of [49], we performed Latin hypercube sampling and generated 5000 samples to calculate the partial rank correlation coefficients (PRCC) and the p-values with respect to the tumor volume at day 30. In sampling all the parameters, we took the range of each parameter from 1/2 to twice its value in Tables 3 or 4. The results are shown in Fig. 10.
We see that parameters that promote the killing of cancer cells, such as λ DC , λ T 8 I 12 and η 8 , are negatively correlated with the tumor volume, while parameters that promote the inhibition of immune cells, such as K DG , K TQ and K TG , and the proliferation of endothelial cells and Tregs, such as λ E , λ T r G and λ GW , are positively correlated with the tumor volume. We also see that the inhibition of perfusion of the anti-PD-1 drug by anti-VEGF promotes cancer growth, namely, K PB is positively correlated to the tumor volume. Among the various parameters, η 8 (the killing rate of cancer cells by CD8 + T cells) has the largest impact in reducing tumor volume, and K TQ (the parameter which increases the inhibition of CD8 + T cells by PD-1-PD-L1) has the largest impact on increasing tumor volume.

Discussion
Combination therapy for cancer has already shown improved benefits over a single agent therapy [50]. But significant challenges remain, as seen, for example, from the fact that of the successful phase II clinical trials with combination chemotherapy, less then 4% demonstrated, in Phase III, improvement of care within 5 years [11]. It was suggested in [13] that a selection of clinical trials needs to include scientific rationale underlying the data and hypothesis for the combination. Combination clinical trials should be preceded by analysis of the expected a b Fig. 7 Spatial profiles of the densities of dendritic cells, T cells, endothelial cells and cancer cells with γ A = 1.2 × 10 −10 g/cm 3 · day, γ B = 9.5 × 10 −9 g/cm 3   3. If the two drugs are injected intermittently, should they be injected simultaneously or non-overlappingly, and, in the latter case, which drug should be injected first?
We developed a mathematical model and simulations for in silico clinical trials that addressed the complexity of interactions among the two drugs. The model is represented by a system of PDEs that includes the most relevant cells and cytokines associated with the treatment. The simulations show that in the case where anti-VEGF decreases the perfusion of the anti-PD-1, the time it takes to reduce tumor volume by 95% is much shorter when the injections of the two drugs are non-overlapping than when the injections are given at the same time. Furthermore, in the non-overlapping treatment, if we inject the anti-PD-1 first we get the 95% -reduction somewhat faster than if we inject anti-VEGF first.
On the other hand, in the case when anti-VEGF increases the perfusion of anti-PD-1, a treatment with simultaneous infections is more beneficial than a treatment with non-overlapping injections. The optimal scheduling in combination therapy for cancer was considered in several publications, both in terms of toxicity [51] and efficacy [52][53][54]. The present paper considered the efficacy question by a mathematical model. The method and results of the paper can be extended to other combinations, for example to a chemotherapeutic agent combined with anti-VEGF, and they could play an important role in the design of clinical trials, where scheduling strategies may significantly affect the outcome.

Conclusions
The design of cancer clinical trials with combination of two drugs should take into account the potential interactions between the two agents, which may suggest an optimal scheduling strategy in the administration of the drugs. Mathematical models can play a useful role in this process. This was illustrated in the case where the drugs are anti-VEGF and anti-PD-1. Anti-VEGF increases the activation of anti-cancer cells, but it also modifies the perfusion of other drugs. Using a mathematical model we showed that a non-overlappling adminstration of the two drugs is more effective in reducing tumor volume than simultaneous administration of the drugs in the case where anti-VEGF degrades perfusion, as it occurs in some cancers, while the opposite is true when perfusion is enhanced by anti-VEGF. The mathematical methodology developed in this paper could be extended to treatment with other combinations of drugs.

Eqs. (3) and (4).
We assume that in steady state, Q/K TQ = 2 (the value of K TQ is derived in the estimates for Eqs. (13)- (17)). We also assume that K TG = 4K G where K G = 7×10 −8 g/cm 3 [43]. From the steady state of Eq. (3), we get where, by [15], λ T 1 I 2 = 0.25/day, d T 1 = 0.197/day, T 10 = 4 × 10 −4 g/cm 3 and T 1 = K T 1 = 2 × 10 −3 g/cm 3 . Hence λ T 1 I 12 = 11.65/day. From the steady state of Eq. (4), we have where by [15], λ T 8 I 2 = 0.25/day, d T 8 = 0.18/day, T 80 = 2 × 10 −4 g/cm 3 , T 8 = K T 8 = 1 × 10 −3 g/cm 3 . Hence λ T 8 I 12 = 10.38/day. As in the case of Eq. (2), we had ignored the decrease in the sources T 10 and T 80 as more of these cells become activated, but we also ignored the contribution of the flux of T cells at the tumor's boundary. We assume that the flux of T cells compensates for decrease in the source of the inactive T cells, and retain the above values of λ T 1 I 12 and λ T 8 I 12 . The same considerations apply in the case of T r .
The molecular mass of anti-VEGF is approximately 6 times larger than that of VEGF, so we take μ GB = 6μ BG = 1.31 × 10 8 cm 3 /g · day.

Computational method
We employ a moving mesh method [48] to numerically solve the free boundary problem for the tumor proliferation model. To illustrate this method, we take Eq. (7) as an example and rewrite it in the following form: where F represents the term in the right hand side of Eq. (7). Let r k i and C k i denote numerical approximations of i-th grid point and C(r k i , nτ ), respectively, where τ is the size of the time-step. The discretization of Eq. (29) is derived by the fully implicit finite difference scheme: where C r = , h −1 = r k+1 i−1 − r k+1 i and h 1 = r k+1 i+1 − r k+1 i . The mesh moves by r k+1 where u k+1 i is solved by the velocity equation. To deal with the different scales of the variables in the simulations, we first non-dimensionalized the model and then applied the above method.