Skip to main content
  • Research article
  • Open access
  • Published:

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



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.


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.


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.


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 anti-cancer 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 [1620]; 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 CD 4+ (Th1) and CD 8+ 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 (antigen-presenting) dendritic cells [25, 26], and it suppresses the functions of activated T cells [1620]; VEGF also enhances the expression of PD-1 on CD 8+ T cells [27], and induces Treg proliferation [28].

Fig. 1
figure 1

Interaction of immune cells with cancer cells. Sharp arrows indicate proliferation/activation, blocked arrows indicate killing/blocking, and the inverted arrow indicates recruitment/chemoattraction. C: cancer cells, D: dentritic cells, T1: CD 4+ Th1 cells, T8: CD 8+ T cells, Treg: T regulatory cells, Endo: endothelial cells, Ox: Oxygen from the blood. T1 and T8 cells and Tregs express PD-1 and PD-L1; tumor expresses PD-L1

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.

Table 1 List of variables (in units of g/ cm3)

We assume that the total density of cells within the tumor remains constant throughout the tumor tissue, for all time:

$$ D+T_{1}+T_{8}+T_{r}+E+C=\text{constant}, $$

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 CD 4+ and CD 8+ T cells, remain constant throughout the tumor tissue. As cancer cells proliferate, they “push away,” or displace, other cells. There is 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.

Fig. 2
figure 2

Distribution of cells in space

Equation for DCs (D)

By necrotic cancer cells (NC) 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 NC and HMGB-1 (H) by the following equations:

$$ {\begin{aligned} &\frac{\partial N_{C}}{\partial t}+\underbrace{\nabla \cdot (\mathbf{u} N_{C})}_{\text{velocity}}-\underbrace{\delta_{N_{C}}\nabla^{2} N_{C}}_{\text{difusion}}=\underbrace{\lambda_{N_{C}C}C}_{\text{derived\ from\ life\ cancer\ cells }}-\underbrace{d_{N_{C}}N_{C}}_{\text{removal}},\\ &\frac{\partial H}{\partial t}-\underbrace{\delta_{H}\nabla^{2} H}_{\text{difusion}}=\underbrace{\lambda_{HN_{C}}N_{C}}_{\text{released from nerotic cancer cells}}-\underbrace{d_{H}H,}_{\text{degradation}}\\ \end{aligned}} $$

where \(\lambda _{N_{C}C}\) is the rate of cancer cells becoming necrotic and \(\lambda _{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 \(\lambda _{N_{C}C}C-d_{N_{C}}N_{C}=0\) and \(\lambda _{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 D0, is proportional to \(D_{0}\frac {H}{K_{H}+H}\), or to \(D_{0}\frac {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:

$$ {\begin{aligned} \frac{\partial D}{\partial t}+\underbrace{\nabla \cdot (\mathbf{u} D)}_{\text{velocity}}-\underbrace{\delta_{D}\nabla^{2} D}_{\text{difusion}}=\underbrace{\lambda_{DC}D_{0}\frac{C}{K_{C}+C}}_{\text{activation by HMGB-1}}\cdot \underbrace{\frac{1}{1+G/K_{DG}}}_{\text{inhibition by VEGF}}-\underbrace{d_{D}D,}_{\text{death}} \end{aligned}} $$

where δD is the diffusion coefficient, dD is the death rate of DCs, and 1/(1+G/KDG) represents the impairment of maturation of dendritic cells by VEGF [25, 26].

Equations for CD 4+ T cells (T 1) and CD 8+ T cells (T 8)

Naive CD 4+ T cells differentiate into active Th1 cells (T1) under IL-12 (I12) environment [34, 35], and this process is inhibited by Tregs [36, 37] and by VEGF [1620]. The proliferation of activated T1 cells is enhanced by IL-2. Both processes of activation and proliferation of T1 are assumed to be inhibited by PD-1-PD-L1 (Q) [38, 39]; we represent this inhibition by the factor \(\frac {1}{1+Q/K_{TQ}}\). Hence T1 satisfies the following equation:

$$ {\begin{aligned} \frac{\partial T_{1}}{\partial t} &+\nabla \cdot (\mathbf{u} T_{1})-\delta_{T}\nabla^{2} T_{1}\\&\,\,\,\,= \underbrace{\left(\lambda_{T_{1}I_{12}}T_{10}\cdot \frac{I_{12}}{K_{I_{12}}+I_{12}}\right.}_{\text{activation by IL-12}}\cdot \underbrace{\frac{1}{1+T_{r}/K_{TT_{r}}}}_{\text{inhibition by Tregs}}\cdot \underbrace{\frac{1}{1+G/K_{TG}}}_{\text{inhibition by VEGF}} \\ &+\underbrace{ \left. \lambda_{T_{1}I_{2}}T_{1}\frac{I_{2}}{K_{I_{2}}+I_{2}}\right)}_{\text{IL-2-induced proliferation}}\times \underbrace{\frac{1}{1+Q/K_{TQ}}}_{\text{inhibition by PD-1-PD-L1}}-\underbrace{d_{T_{1}}T_{1},}_{\text{death}} \end{aligned}} $$

where T10 is the density of the naive CD 4+ T cells.

Similarly, inactive CD 8+ T cells are activated by IL-12 [34, 35], a process resisted by Tregs [36, 37] and VEGF [1620], while IL-2 enhances the proliferation of activated CD 8+ T cells. Hence,

$$ {\begin{aligned} \frac{\partial T_{8}}{\partial t}&+\nabla \cdot (\mathbf{u} T_{8})-\delta_{T}\nabla^{2} T_{8}\\&\,\,\,\,=\underbrace{\left(\lambda_{T_{8}I_{12}}T_{80}\cdot \frac{I_{12}}{K_{I_{12}}+I_{12}}\right.}_{\text{activation by IL-12}}\cdot \underbrace{\frac{1}{1+T_{r}/K_{TT_{r}}}}_{\text{inhibition by Tregs}}\cdot \underbrace{\frac{1}{1+G/K_{TG}}}_{\text{inhibition by VEGF}}\\ &+\underbrace{\left.\lambda_{T_{8}I_{2}}T_{8}\frac{I_{2}}{K_{I_{2}}+I_{2}}\right)}_{\text{IL-2-induced proliferation}} \times \underbrace{\frac{1}{1+Q/K_{TQ}}}_{\text{inhibition by PD-1-PD-L1}}-\underbrace{d_{T_{8}}T_{8},}_{\text{death}} \end{aligned}} $$

where T80 is the density of the inactive CD 8+ T cells.

Equation for Tregs (T r)

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

$$ \begin{aligned} \frac{\partial T_{r}}{\partial t}+\nabla \cdot (\mathbf{u} T_{r})-\delta_{T}\nabla^{2} T_{r}=&T_{10}\underbrace{\left(\lambda_{T_{r}T_{\beta}}\frac{T_{\beta}}{K_{T_{\beta}}+T_{\beta}}\right.}_{\text{TGF}-\beta-\text{induced proliferation}}\\&+\underbrace{\lambda_{T_{r}G}\left.\frac{G}{K_{G}+G}\right)}_{\text{promotion by VEGF}}-\underbrace{d_{T_{r}}T_{r}.}_{\text{death}} \end{aligned} $$

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

$$ \begin{aligned} \frac{\partial E}{\partial t}+\nabla \cdot (\mathbf{u} E)-\delta_{E}\nabla^{2} E=&\underbrace{\lambda_{E}(G)E\left(1-\frac{E}{E_{M}}\right)}_{\text{proliferation}}\\&-\underbrace{\nabla \cdot (\chi_{G} E\nabla G)}_{\text{recruited by VEGF}}-\underbrace{d_{E}E,}_{\text{death}} \end{aligned} $$

where EM is the carrying capacity of endothelial cells, λE(G)=λEG(GG0)+, and G0 is a threshold below which endothelial cells do not proliferate [43]. Here we used the notion: X+=X if X≥0 and X+=0 if X<0.

Equation for cancer cells (C)

We assume a logistic growth for cancer cells, with carrying capacity (CM), 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:

$$ \begin{aligned} \frac{\partial C}{\partial t}+\nabla \cdot (\mathbf{u} C)-\delta_{C}\nabla^{2} C=&\underbrace{\lambda_{C}(W)C\left(1-\frac{C}{C_{M}}\right)}_{\text{proliferation}}\\*&-\underbrace{\eta_{1}T_{1}C-\eta_{8} T_{8}C}_{\text{killing by T cells}}-\underbrace{d_{C}C,}_{\text{death}} \end{aligned} $$

where η1 and η8 are the killing rates of cancer cells by T1 and T8 cells, respectively. dC is the natural death rate of cancer cells, and

$$\lambda_{C}(W)=\left\{ \begin{array}{ll} \lambda_{CW}\frac{W}{W_{0}} & \,\,\text{if}\:\: W\le W_{0}\\ \lambda_{CW} & \,\,\text{if}\:\: W> W_{0}, \end{array} \right. $$

where W0 is the normal level of oxygen concentration in the blood.

Equation for IL-12 (I 12)

The proinflammatory cytokine IL-12 is secreted by activated DCs [34, 35]; hence it satisfies the equation:

$$\begin{array}{@{}rcl@{}} \frac{\partial I_{12}}{\partial t}-\delta_{I_{12}}\nabla^{2} I_{12}&=&\underbrace{\lambda_{I_{12}D}D}_{\text{production by DCs}}-\underbrace{d_{I_{12}}I_{12}.}_{\text{degradation}} \end{array} $$

Equation for IL-2 (I 2)

IL-2 is produced by activated CD 4+ T cells (T1) [35]. Hence,

$$\begin{array}{@{}rcl@{}} \frac{\partial I_{2}}{\partial t}-\delta_{I_{2}}\nabla^{2} I_{2}&=&\underbrace{\lambda_{I_{2}T_{1}}T_{1}}_{\text{production by \(T_{1}\)}}-\underbrace{d_{I_{2}}I_{2}.}_{\text{degradation}} \end{array} $$

Equation for TGF-β: (T β)

TGF-β is produced by tumor cells [36] and Tregs [37], so that

$$ \begin{aligned} \frac{\partial T_{\beta}}{\partial t}-\delta_{T_{\beta}}\nabla^{2} T_{\beta}=&\underbrace{\lambda_{T_{\beta}C}C}_{\text{production by cancer cells}}\\&+\underbrace{\lambda_{T_{\beta}T_{r}}T_{r}}_{\text{production by Tregs}}-\underbrace{d_{T_{\beta}}T_{\beta}.}_{\text{degradation}} \end{aligned} $$

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:

$$ \begin{aligned} \frac{\partial W}{\partial t}-\delta_{W}\nabla^{2} W=&\underbrace{\lambda_{WE}E}_{\text{source from blood}}-\underbrace{d_{W}W,}_{\text{consumption by cells}} \end{aligned} $$

where dWW 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:

$$ \begin{aligned} \frac{\partial G}{\partial t}-\delta_{G}\nabla^{2} G=&\underbrace{\lambda_{G}(W)C}_{\text{production by cancer cells}}- \underbrace{\mu_{GB}GB}_{\text{inhibition by anti-VEGF}}\\ &-\underbrace{d_{G}G}_{\text{degradation}}, \end{aligned} $$

where B is the effective anti-VEGF concentration in the tumor, and

$$\lambda_{G}(W)=\lambda_{GW}\times \left\{ \begin{array}{ll} \frac{W}{W^{*}} & \,\,\text{if}\:\: 0\le W\le W^{*}\\ 1-0.7\frac{W-W^{*}}{W_{0}-W^{*}} & \,\,\text{if}\:\:W^{*}< W\le W_{0}\\ 0.3 & \,\,\text{if}\:\: W> W_{0}.\\ \end{array}\right. $$

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 CD 4+ T cells, activated CD 8+ T cells, and Tregs [38, 45]. We assume that the number of PD-1 receptors per cell is the same for T1 and T8 cells, but is smaller, by a factor εT, for Tr cells. VEGF increases the PD-1 on T8 cells by a factor εGG [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 T1 and Tr cells is given by

$$ P_{1}=\rho_{P}\left(T_{1}+\varepsilon_{T} T_{r}\right), $$

and the concentration of PD-1 on T8 cells is given by

$$ P_{8}=\rho_{P}T_{8}(1+\varepsilon_{G}G). $$

PD-L1 is expressed on activated CD 4+ T cells, activated CD 8+ T cells [38], Tregs [46], and cancer cells [38, 39]. We assume that the number of PD-L1 per cell is the same for T1 and T8 cells, and denote the ratio between the mass of one PD-L1 protein to the mass of one cell by ρL. Then

$$L=\rho_{L}(T_{1}+T_{8}+\varepsilon_{T}T_{r}+\varepsilon_{C} C), $$

for some parameters εT,εC, where εC depends on the specific type of tumor.

To a change in T=T1+εTTr, given by \(\frac {\partial T}{\partial t}\), there corresponds a change of P1, given by \(\rho _{P}\frac {\partial T}{\partial t}\). For the same reason, ·(uP1)=ρP·(uT) and 2P1=ρP2T when no anti-PD-1 drug is injected. Hence, P1 satisfies the following equation:

$$ \begin{aligned} \frac{\partial P_{1}}{\partial t}+\nabla \cdot (\mathbf{u} P_{1})-\delta_{T} \nabla^{2} P_{1}=&\rho_{P}\left[\frac{\partial(T_{1}+\varepsilon_{T}T_{r})}{\partial t} \right.\\ &+\nabla \cdot (\mathbf{u} (T_{1}+\varepsilon_{T}T_{r}))\\ &\left.-\delta_{T}\nabla^{2} (T_{1}+\varepsilon_{T}T_{r}){\vphantom{\frac{\partial(T_{1}+\varepsilon_{T}T_{r})}{\partial t}}}\right]. \end{aligned} $$

Recalling Eqs. (3) and (5) for T1 and Tr, we get,

$$ \begin{aligned} \frac{\partial P_{1}}{\partial t}+&\nabla \cdot (\mathbf{u} P_{1})-\delta_{T} \nabla^{2} P_{1}\\&=\rho_{P} \left[\text{RHS of Eq.\:(3)}+\varepsilon_{T},\times \text{RHS of Eq.\:(5)} \right] \end{aligned} $$

where RHS=right-hand side. Similarly,

$$ {\begin{aligned} \frac{\partial P_{8}}{\partial t}+&\nabla \cdot (\mathbf{u} P_{8})-\delta_{T} \nabla^{2} P_{8}\\ &=\rho_{P} \left[(1+\varepsilon_{G}G)\times \text{RHS of Eq.\:(4)} + \varepsilon_{G} T_{8} \times\text{RHS of Eq.\:(12)}\right.\\ & \left.+ \varepsilon_{G}T_{8}\nabla\cdot\mathbf{u}G+\varepsilon_{G}(\delta_{G}-\delta_{T})T_{8}\nabla^{2}G-2\varepsilon_{G}\delta_{T}\nabla T_{8}\cdot \nabla G \right]. \end{aligned}} $$

When only anti-PD-1 drug (A) is injected, PD-1 is depleted at a rate μPAP1A, but when also anti-VEGF (B) is injected the depletion of PD-1 is reduced, due to restricted perfusion, by a factor \(\frac {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 P1 satisfies the following equation:

$$ \begin{aligned} \frac{\partial P_{1}}{\partial t}+&\nabla \cdot (\mathbf{u} P_{1})-\delta_{T} \nabla^{2} P_{1}\\ &=\rho_{P}\left[\text{RHS of Eq.\:(3)} +\varepsilon_{T}\times \text{RHS of Eq.\:(5)}\right]\\ &-\underbrace{\mu_{PA} P_{1}A}_{\text{ depletion\ by\ anti-PD-1 }}\times \underbrace{\frac{1}{1+B/K_{PB}}.}_{\text{ blockade\ of\ perfusion\ by\ anti-VEGF }} \end{aligned} $$


$$ {\begin{aligned} \frac{\partial P_{8}}{\partial t}+&\nabla \cdot (\mathbf{u} P_{8})-\delta_{T} \nabla^{2} P_{8}\\&=\rho_{P}\left[(1+\varepsilon_{G}G)\times \text{RHS of Eq.\:(4)}\right.\\ & +\varepsilon_{G}T_{8} \times\text{RHS of Eq.\:(12)}\\ &\left. + \varepsilon_{G} T_{8}\nabla\cdot\mathbf{u}G+\varepsilon_{G}(\delta_{G}-\delta_{T})T_{8}\nabla^{2}G -2\varepsilon_{G}\delta_{T}\nabla T_{8}\cdot \nabla G \right] \\ & -\underbrace{\mu_{PA} P_{8}A}_{\text{ depletion\ by\ anti-PD-1 }}\times \underbrace{\frac{1}{1+B/K_{PB}}.}_{\text{ blockade\ of\ perfusion\ by\ anti-VEGF }} \end{aligned}} $$

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 dQ, respectively, we can write

$$ P+L \overset{\alpha_{PL}}{\underset{d_{Q}}{\rightleftharpoons}} Q, $$

where P=P1+P8. The half-life of Q is less then 1 s (1.16×10−5day) [47], so that dQ is very large. Hence we may approximate the dynamical equation for Q by the steady state equation, αPLPL=dQQ, so that

$$ Q=\sigma PL, $$

where σ=αPL/dQ.

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,

$$ {\begin{aligned} \frac{\partial A}{\partial t}-\delta_{A} \nabla^{2} A&=&\gamma_{A}I_{A}(t) -\underbrace{\mu_{AP} PA}_{\text{depletion through blocking PD-1}}-\underbrace{d_{A}A,}_{\text{degradation}} \end{aligned}} $$

where IA(t)=1 during dosing, and IA(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

$$ \begin{aligned} \frac{\partial B}{\partial t}-\delta_{B} \nabla^{2} B=&\gamma_{B}I(t) -\underbrace{\mu_{BG}GB}_{\text{depletion by VEGF}}-\underbrace{d_{B}B,}_{\text{degradation}} \end{aligned} $$

where IB(t)=1 during dosing, and IB(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/cm3.

We take the steady state density of endothelial cells to be E=2.5×10−3 g/ cm3 [43]. The steady state densities of the immune cells D, T1,T8,Tr and C (in units of g/cm3) are taken to be (see Appendix: Parameter estimation)

$$ {\begin{aligned} D=4\times 10^{-4}, \ \ T_{1}=2\times 10^{-3},\ \ T_{8}=1\times 10^{-3}, \ \ T_{r}= 5\times 10^{-4}, \ \ C=0.4. \end{aligned}} $$

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

$$ 0.4064\times \nabla\cdot \mathbf{u}=\sum_{j=2}^{7} \left[\text{RHS of Eq.\:(j)} \right]. $$

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≤rR(t). In particular, u=u(r,t)er, where er 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

$$ \frac{dR(t)}{dt}=u(R(t),t). $$

Boundary conditions

We assume that the naive CD 4+ T cells and CD 8+ T cells which migrated from the lymph nodes into the tumor microenvironment have constant densities \(\hat T_{1}\) and \(\hat T_{8}\) at the tumor’s boundary, and, that upon crossing the tumor boundary, T1 and T8 are activated by IL-12 and Tr is activated by Tβ. We then have the following flux conditions at the tumor’s boundary:

$$ {\begin{aligned} &\frac{\partial T_{1}}{\partial r}+\sigma_{T}(I_{12})(T_{1}-\hat T_{1})=0,\:\: \frac{\partial T_{8}}{\partial r}+\sigma_{T}(I_{12})(T_{8}-\hat T_{8})=0,\\ &\frac{\partial T_{r}}{\partial r}+\sigma_{T_{r}}(T_{\beta})(T_{r}-\hat T_{1})=0\quad \text{at}\:\: r=R(t), \end{aligned}} $$

where we take \(\sigma _{T}(I_{12})=\sigma _{0} \frac {I_{12}}{K_{I_{12}}+I_{12}}\) and \(\sigma _{T_{r}}(T_{\beta })=\sigma _{0} \frac {T_{\beta }}{K_{T_{\beta }}+T_{\beta }}\). We assume that the endothelial cells that are attracted by the VEGF into the tumor microenvironment have constant density \(\hat E\) at the tumor’s boundary, and take

$$ \begin{aligned} &\frac{\partial E}{\partial r}+\sigma_{E}\frac{G}{K_{G}+G}(E-\hat E)=0 \quad \text{at}\:\: r=R(t). \end{aligned} $$

The boundary condition for the oxygen is given by

$$ \begin{aligned} &\frac{\partial W}{\partial r}+\sigma_{W}(W-W_{0})=0 \quad \text{at}\:\: r=R(t), \end{aligned} $$

where W0 is the normal level of oxygen concentration in the blood. We impose no-flux boundary condition for all the remaining variables:

$$ \begin{aligned} \text{No-flux for} &\:D,C,I_{12},I_{2},T_{\beta},G,P_{1},P_{8}, A,\:\text{and} \:B\: \text{at}\: r=R(t). \end{aligned} $$

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/cm3):

$$ {\begin{aligned} &D=2\times 10^{-3}, \: T_{1}\!=4\times 10^{-3},\:T_{8}\!=2\times 10^{-3}, \: T_{r}\!= 3\times 10^{-4}, \:E=2.45\times 10^{-3},\\ & C= 0.39565, \: I_{12}=2.88\times 10^{-9},\: I_{2}=4.74\times 10^{-11},\: T_{\beta}=2.57\times 10^{-13},\\ & W=1.52\times 10^{-4},\: G=6.3\times 10^{-8},\:P_{1}=1.06\times 10^{-9},\:P_{8}=5.43\times 10^{-11}, \\ & R(0)=0.01\:\text{cm}. \end{aligned}} $$

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.


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 half-saturation value, KX, as assumed in the estimation of some of the model parameters (in Appendix). This shows the consistency in the parameters estimation.

Fig. 3
figure 3

Average densities/concentrations, in g/cm3, of all the variables in the model with control case (no drugs). All parameter values are the same as in 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 T1 and T8: 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−8g/cm3·day and γA=3×10−8g/cm3·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−8g/cm3·day and γA=0.5×10−8g/cm3·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−8g/cm3·day reduces significantly the PD-1 expression on CD 8+ T cells, in agreement with experimental results in [27].

Fig. 4
figure 4

Growth of tumor volume under treatment with γB or γA, or combination (γ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−8g/cm3·day,γA=3×10−8g/cm3·day; b γB=3.5×10−8g/cm3·day,γA=0.5×10−8g/cm3·day. All other parameters are same as in Fig. 3

Fig. 5
figure 5

Anti-VEGF decreases PD-1 expression on CD 8+ T cells. The treatment with anti-PD-1 drug started at day 0 and continued for 30 days with γB=2×10−8g/cm3·day (a) Growth of tumor volume. (b) Expression level of PD-1 on CD 8+ T cells

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.

Fig. 6
figure 6

Tumor volume under schedules S1, S2 and S3 for 4 pairs (γB,γA). Here γB1=9.5×10−9g/cm3·day,γA1=1.2×10−10g/cm3·day,γB2=10×10−9g/cm3·day,γA2=1.5×10−10g/cm3·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

Figure 7 displays the density profiles of the immune cells, endothelial cells and cancer cells at different times along the radius of the tumor, with γA=1.2×10−10g/cm3·day,γB=9.5×10−9g/cm3·day. Figure 7a shows the profiles under schedule S2 at times t=0,8,15 and 16 weeks, and Fig. 7b shows the profiles under schedule S3 at times t=0,8,18 and 19 weeks. We see that the initially assumed flat density profiles are evolving to develop a distinct monotonic form. The density of cancer cells (C) is increasing from the tumor core to the tumor boundary where the level of oxygen (W) is more favorable to their abnormal proliferation. Since the total density of all the cells is constant (by Eq. (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.

Fig. 7
figure 7

Spatial profiles of the densities of dendritic cells, T cells, endothelial cells and cancer cells with γA=1.2×10−10g/cm3·day,γB=9.5×10−9g/cm3·day. a Profiles under schedule S2 at day t=0,8,15 and 16 weeks. b Profiles under schedule S3 at days t=0,8,18 and 19 weeks

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/cm3·day, and the vertical axis scales the dose amount of anti-PD-1, γA, in unit of g/cm3·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.

Fig. 8
figure 8

Efficacy maps: The time in weeks (Tcrit) at which the tumor volume decreases by 95% from its initial size under treatment with (γB,γA). a Efficacy map under schedule S1; b Efficacy map under schedule S2; c Efficacy map under schedule S3. The color columns show the time at which the tumor volume was reduced by 95%

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/KPB) by \(\left (1+\frac {B}{K_{B}+B}\right)\). 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.

Fig. 9
figure 9

Tumor volume under the schedules S1, S2 and S3 for 4 pairs (γB,γA). Here γB1=9.5×10−9g/cm3·day,γA1=1.2×10−10g/cm3·day,γB2=9.5×10−9g/cm3·day,γA2=1.5×10−10g/cm3·day. a Tumor volume under schedule S1; b Tumor volume under schedule S2; c Tumor volume under schedule S3

Sensitivity analysis

We performed sensitivity analysis with respect to the tumor volume at day 30 for parameters \(\phantom {\dot {i}\!}\lambda _{DC}, \lambda _{T_{8}I_{12}}, \lambda _{T_{r}G}, \lambda _{E}, \lambda _{GW}, \eta _{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.

Table 3 Summary of parameter values
Table 4 Summary of parameter values
Fig. 10
figure 10

Statistically significant PRCC values (p-value <0.01) for tumor volume at day 30

We see that parameters that promote the killing of cancer cells, such as \(\phantom {\dot {i}\!}\lambda _{DC}, \lambda _{T_{8}I_{12}}\) and η8, are negatively correlated with the tumor volume, while parameters that promote the inhibition of immune cells, such as KDG,KTQ and KTG, and the proliferation of endothelial cells and Tregs, such as \(\phantom {\dot {i}\!}\lambda _{E}, \lambda _{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, KPB is positively correlated to the tumor volume. Among the various parameters, η8 (the killing rate of cancer cells by CD 8+ T cells) has the largest impact in reducing tumor volume, and KTQ (the parameter which increases the inhibition of CD 8+ T cells by PD-1-PD-L1) has the largest impact on increasing tumor volume.


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 interactions between the diverse agents, addressing, for example the following questions:

  1. 1.

    Are the two drugs positively correlated at any amount of dose, and, if so, what should be the most beneficial ratio between the two agents?

  2. 2.

    If the drugs are not always positively correlated, what are the zones of antagonism, i.e., what are the amounts of, and the ratios between, the agents in the combination that may decrease the treatment benefits, and should not be used in clinical trials?

  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?

In [14] and [15] we addressed the questions (i) and (ii) in the case where one of the drugs is a checkpoint inhibitor (anti-PD-1), and the second drug is a cancer vaccine [14] or a BRAF-inhibitor [15]. In the present paper we addressed the question (iii) when the two drugs are anti-PD-1 and anti-VEGF. Anti-VEGF is known to block the perfusion of chemotherapy in melanoma [3], breast cancer [4, 5, 9] and ovarian cancer [6], but to increase perfusion in colon cancer [7].

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 [5254]. 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.


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.


Parameter estimation


In an expression of the form \(Y\frac {X}{K_{X}+X}\) where Y is activated by X, the parameter KX is called the half-saturation of X.

We assume that if the average density (or concentration) of X,

$$\frac{\int X dx}{\int dX}, $$

converges to a steady state X0, then

$$\frac{X_{0}}{K_{X}+X_{0}} $$

is not “too small” and not “too close” to 1, and for definiteness we take

$$\frac{X_{0}}{K_{X}+X_{0}}=\frac{1}{2} $$

so that

$$ K_{X}=X_{0}. $$

To estimate parameters, we assume that all the average densities and concentrations converge to their respective steady states, and thus use Eq. (28) in each of the steady state equations of the model.

Diffusion coefficients

We use the following relation for estimating the diffusion coefficient of a protein p [55]:

$$\delta_{p}=\frac{M^{1/3}_{G}}{M^{1/3}_{p}}\delta_{G}, $$

where MG and δG are respectively the molecular weight and diffusion coefficient of VEGF, Mp is the molecular weight of p, and MG=24kDa [56] and δG=8.64×10−2cm2 day−1 [57]. Since, MB=149kDa (bevacizumab), we get δB=4.70×10−2cm2 day−1. The diffusion coefficient of oxygen in the extracellular matrix (ECM) in the range of 7×10−6−2×10−5cm2/s [58]; we take it to be δW=0.8cm2/day.

Eq. (2).

From the steady state of Eq. (2) (more precisely, by setting to zero the RHS of Eq. (2)), we get \(\lambda _{DC}D_{0}\frac {C}{K_{C}+C}\cdot \frac {1}{1+G/K_{DG}}=d_{D}D\), where by [15], dD=0.1/day, C=KC=0.4g/cm3,D=KD=4×10−4g/cm3,D0=2×10−5g/cm3. We assume that KDG=4KG where KG=7×10−8g/cm3 [43]; hence λDC=2.5dDD/D0=5/day. For simplicity we had assumed that the source of inactive dendritic cells, D0, is constant. However, this source actually decreases as more dendritic cells become activated. We take this into count by increasing λDC, taking λDC=17.5/day in mice and λDC=7.5/day in humans.

Eqs. (3) and (4).

We assume that in steady state, Q/KTQ=2 (the value of KTQ is derived in the estimates for Eqs. (13)-(17)). We also assume that KTG=4KG where KG=7×10−8g/cm3 [43]. From the steady state of Eq. (3), we get

$$\left(\lambda_{T_{1}I_{12}}T_{10}\cdot\frac{1}{2}\cdot\frac{1}{2}\cdot\frac{4}{5}+\lambda_{T_{1}I_{2}}T_{1}\cdot\frac{1}{2}\right)\cdot\frac{1}{3}-d_{T_{1}}T_{1}=0, $$

where, by [15], \(\lambda _{T_{1}I_{2}}=0.25\)/day, \(d_{T_{1}}=0.197\)/day, T10=4×10−4g/cm3 and \(T_{1}=K_{T_{1}}=2\times 10^{-3} \mathrm {g}/\text {cm}^{3}\). Hence \(\lambda _{T_{1}I_{12}}=11.65/\text {day}\).

From the steady state of Eq. (4), we have

$$\left(\lambda_{T_{8}I_{12}}T_{80}\cdot\frac{1}{2}\cdot\frac{1}{2}\cdot\frac{4}{5}+\lambda_{T_{1}I_{2}}T_{8}\cdot\frac{1}{2}\right)\cdot\frac{1}{3}-d_{T_{8}}T_{8}=0 $$

where by [15], \(\lambda _{T_{8}I_{2}}=0.25\)/day, \(d_{T_{8}}=0.18\)/day, \(T_{80}=2\times 10^{-4} \mathrm {g}/\text {cm}^{3}, T_{8}=K_{T_{8}}=1\times 10^{-3} \mathrm {g}/\text {cm}^{3}\). Hence \(\lambda _{T_{8}I_{12}}=10.38/\text {day}\).

As in the case of Eq. (2), we had ignored the decrease in the sources T10 and T80 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 \(\lambda _{T_{1}I_{12}}\) and \(\lambda _{T_{8}I_{12}}\). The same considerations apply in the case of Tr.

Eq. (5).

We assume that TGF-β activates Tregs more than VEGF does, and take \(\lambda _{T_{r}T_{\beta }}=5\lambda _{T_{r}G}\). From the steady state of Eq. (5), we get, \((\lambda _{T_{r}T_{\beta }}\cdot \frac {1}{2}+\lambda _{T_{r}G}\cdot \frac {1}{2})T_{10}-d_{T_{r}}T_{r}=0\), where \(T_{10}=1\times 10^{-3} \mathrm {g}/\text {cm}^{3}, T_{r}=K_{T_{r}}=5\times 10^{-4} \mathrm {g}/\text {cm}^{3}\) [59], and \(d_{T_{r}}=0.2\)/day [59]. Hence \(\lambda _{T_{r}G}=0.083/\text {day}\) and \(\lambda _{T_{r}T_{\beta }}=0.415\)/day.

Eq. (6).

By [43], dE=0.69/day, EM=5×10−3g/cm3,KE=2.5×10−3g/cm3,G0=3.65×10−10. We take EM=2KE=5×10−3g/cm3. From the steady state of Eq. (6), we get λE=2dE/(KGG0)=1.98×107cm3/g·day. Assuming that the threshold G0 (which was taken to be constant) is actually increasing with the progression of the cancer in the control case (no drugs), we increase λE, taking λE=2.77×107cm3/g·day in mice and λE=2.08×107cm3/g·day in humans.

Eq. (7).

We take dC=0.17day−1,CM=0.8g/cm3 [59] and λCW=1.6/day [60]. In the steady state of the control case (no drugs), we assume that C is approximately 0.4 g/cm3, and W=W0=KW=1.69×10−4g/cm3 (see the estimates for Eq. (11)). From the steady state of Eq. (7) in the control case we have,

$$\frac{1}{2}\lambda_{CW}K_{W}/W_{0}-\eta_{1}K_{T_{1}}-\eta_{8}K_{T_{8}}-d_{C}=0, $$

where \(K_{T_{1}}=2\times 10^{-3} \mathrm {g}/\text {cm}^{3},K_{T_{8}}=1\times 10^{-3} \mathrm {g}/\text {cm}^{3}\); we take η8=2η1; hence \(\eta _{8}=(\lambda _{C}K_{W}/(2W_{0})-d_{C})/(2K_{T_{8}})=60.375 \text {cm}^{3}/\mathrm {g}\cdot \text {day}\) and η1=30.19cm3/g·day. Since in the control case the tumor grows, we increase the growth rate of cancer cells taking λCW=2.24/day in mice and λCW=1.76/day in humans.

Eq. (10).

From the steady state of Eq. (10) we have, \(\lambda _{T_{\beta } C}C+\lambda _{T_{\beta } T_{r}}T_{r}=d_{T_{\beta }}T_{\beta }\), where \(T_{r}=K_{T_{r}}=5\times 10^{-4} \mathrm {g}/\text {cm}^{3},C=K_{C}=0.4 \mathrm {g}/\text {cm}^{3}\), and, by [15], \(T_{\beta }=K_{T_{\beta }}=2.68\times 10^{-13} \mathrm {g}/\text {cm}^{3}\) and \(d_{T_{\beta }}=499.07 \text {day}^{-1}\), and by [61], \(\lambda _{T_{\beta } T_{r}}=5.57\times 10^{-9}\)/day. Hence \(\lambda _{T_{\beta } C}=3.27\times 10^{-10}\)/day.

Eq. (11).

From steady state of Eq. (11) we get, λWEEdWW=0, where λWE=7×10−2/day [43], W=KW=1.69×10−4g/cm3 [43], E=KE=2.5×10−3g/cm3. Hence, dW=λWEE/W=1.04/day.

Eq. (12).

From steady state of Eq. (12) we get, λGWCdGG=0, where dG=12.6/day [43], G=KG=7×10−8g/cm3 [43], C=KC=0.4g/cm3. Hence, λGW=2.21×10−6/day.

Eqs. (13)-(17).

By [15], ρP=2.49×10−7,ρL=5.22×10−7,εT=0.8,εC=0.01. We assume that εG=0.1/KG=1.43×106cm3/g. From Eqs. (13)-(16) we get,

$$\begin{array}{@{}rcl@{}} K_{P_{1}}=P_{1}&=&\rho_{P}(T_{1}+\varepsilon_{T} T_{r})\\ &=&(2.49\!\times\! 10^{-7})\times\left[2\!\times\! 10^{-3} \,+\,0.8\!\times\! (5 \times 10^{-4})\right]\\ &=&5.98\times 10^{-10}\mathrm{g}/\text{cm}^{3}, \end{array} $$
$$\begin{array}{@{}rcl@{}} K_{P_{8}}=P_{8}&=&\rho_{P}T_{8}(1+\varepsilon_{G} G)\\ &=&\left(2.49\times 10^{-7}\right)\times\left(1\times 10^{-3}\right)\times (1+0.1)\\ &=&2.74\times 10^{-10}\mathrm{g}/\text{cm}^{3}, \end{array} $$


$$\begin{array}{@{}rcl@{}} K_{L}=L&=&\rho_{L}(T_{1}+T_{8}+\varepsilon_{T} T_{r}+\varepsilon_{C} C)\\ &=&\left(5.22\times 10^{-7}\right)\times\left[3.4\times 10^{-3} +0.01\times 0.4\right]\\ &=&3.86\times 10^{-9}\mathrm{g}/\text{cm}^{3}. \end{array} $$

In steady state with \(\phantom {\dot {i}\!}P=K_{P}=K_{P_{1}}+K_{P_{8}},L=K_{L},Q=K_{Q}\) and G=KG we have, by Eq. (17), KQ=σKPKL. We take \(K_{TQ}=K_{Q}=\frac {1}{2}\sigma K_{P}K_{L}\). Hence, \(Q/K_{TQ}=PL/\left (\frac {1}{2}K_{P}K_{L}\right)\) and

$$\frac{1}{1+Q/K_{TQ}}=\frac{1}{1+PL/\left(\frac{1}{2}K_{P}K_{L}\right)}=\frac{1}{1+PL/K'_{TQ}}, $$

where \(K^{\prime }_{TQ}=\frac {1}{2}K_{P}K_{L}=\frac {1}{2}\times \left (5.98\times 10^{-10}+2.74\times 10^{-10}\right)\times \left (3.86\times 10^{-9}\right)=1.68\times 10^{-18} \mathrm {g}^{2}/\text {cm}^{6}\).

Eqs. (18)-(19).

We take dA=0.34day−1 and assume that 10% of A is used in blocking PD-1, while the remaining 90% degrades naturally. Hence,

$${{} \begin{aligned} \mu_{AP}\,=\,\frac{d_{A}}{9P}\,=\,\frac{0.34}{9\times (8.715\!\times\! 10^{-10})}\,=\,4.33\!\times\! 10^{7}\: \text{cm}^{3}/\mathrm{g}\cdot \text{day}. \end{aligned}} $$

Since the molecular mass of anti-PD-1 (32kDa [62]) approximates the molecular mass of PD-1 (20.5 - 40 kDa [62]), we take μPA=μAP=4.33×107 cm3/g·day.

By [63], the half-life of anti-VEGF is 2.82-4.58 days; we take it to be 4 days, so that \(d_{B}=\frac {\text {ln} 2}{4}=0.17 \text {day}^{-1}\). We assume that 90% of B is depleted in blocking of VEGF, while the remaining 10% degrades naturally. Hence,

$$\mu_{BG}=\frac{9d_{B}}{G}=\frac{9\time 0.17}{7\times 10^{-8}}=2.19\times 10^{7}\: \text{cm}^{3}/\mathrm{g}\cdot \text{day}. $$

The molecular mass of anti-VEGF is approximately 6 times larger than that of VEGF, so we take μGB=6μBG=1.31×108 cm3/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:

$$ \begin{aligned} \frac{\partial C(r,t)}{\partial t}=\delta_{C}\Delta C(r,t)-div(\mathbf{u}C)+ F, \quad \end{aligned} $$

where F represents the term in the right hand side of Eq. (7). Let \(r_{i}^{k}\) and \(C_{i}^{k}\) denote numerical approximations of i-th grid point and \(C(r_{i}^{k},n\tau)\), respectively, where τ is the size of the time-step. The discretization of Eq. (29) is derived by the fully implicit finite difference scheme:

$$ \begin{aligned} \frac{C_{i}^{k+1}-C_{i}^{k}}{\tau}=&\delta_{C}\left(C_{rr}+\frac{2}{r_{i}^{k}}C_{r}\right)-\left(\frac{2}{r_{i}^{k+1}}u_{i}^{k+1}+u_{r}\right)C_{i}^{k+1}\\&-u_{i}^{k+1}C_{r}+F_{i}^{k+1}, \end{aligned} $$

where \(C_{r}=\frac {h_{-1}^{2}C_{i+1}^{k+1}-h_{1}^{2}C_{i-1}^{k+1}-(h_{1}^{2}-h_{-1}^{2})C_{i}^{k+1}}{h_{1}(h_{-1}^{2}-h_{1}h_{-1})},C_{rr}=2\frac {h_{-1}C_{i+1}^{k+1}-h_{1}C_{i-1}^{k+1}+(h_{1}-h_{-1})C_{i}^{k+1}}{h_{1}(h_{1}h_{-1}-h_{-1}^{2})}\),\(u_{r}=\frac {h_{-1}^{2}u_{i+1}^{k+1}-h_{1}^{2}u_{i-1}^{k+1}-(h_{1}^{2}-h_{-1}^{2})u_{i}^{k+1}}{h_{1}(h_{-1}^{2}-h_{1}h_{-1})},h_{-1}=r_{i-1}^{k+1}-r_{i}^{k+1}\) and \(h_{1}=r_{i+1}^{k+1}-r_{i}^{k+1}\). The mesh moves by \(r_{i}^{k+1}=r_{i}^{k}+u_{i}^{k+1}\tau \), where \(u_{i}^{k+1}\) 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.


  1. Meadows KL, Hurwitz HI. Anti-vegf therapies in the clinic. Cold Spring Harb Perspect Med. 2012; 2(10):006577.

    Article  CAS  Google Scholar 

  2. Nerini IF, Cesca M, Bizzaro F, Giavazzi R. Combination therapy in cancer: effects of angiogenesis inhibitors on drug pharmacokinetics and pharmacodynamics. Chin J Cancer. 2016; 35(1):61.

    Article  Google Scholar 

  3. Turley RS, Fontanella AN, Padussis JC, Toshimitsu H, Tokuhisa Y, Cho EH, Hanna G, Beasley GM, Augustine CK, Dewhirst MW, et al.Bevacizumab-induced alterations in vascular permeability and drug delivery: a novel approach to augment regional chemotherapy for in-transit melanoma. Clin Cancer Res. 2012; 18(12):3328–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Pastuskovas CV, Mundo EE, Williams SP, Nayak TK, Ho J, Ulufatu S, Clark S, Ross S, Cheng E, Parsons-Reponte K, et al.Effects of anti-vegf on pharmacokinetics, biodistribution, and tumor penetration of trastuzumab in a preclinical breast cancer model. Mol Cancer Ther. 2012; 11(3):752–62.

    Article  CAS  PubMed  Google Scholar 

  5. Daldrup-Link HE, Okuhata Y, Wolfe A, Srivastav S, Øie S, Ferrara N, Cohen RL, Shames DM, Brasch RC. Decrease in tumor apparent permeability-surface area product to a mri macromolecular contrast medium following angiogenesis inhibition with correlations to cytotoxic drug accumulation. Microcirculation. 2004; 11(5):387–96.

    Article  CAS  PubMed  Google Scholar 

  6. Cesca M, Morosi L, Berndt A, Nerini IF, Frapolli R, Richter P, Decio A, Dirsch O, Micotti E, Giordano S, et al.Bevacizumab-induced inhibition of angiogenesis promotes a more homogeneous intratumoral distribution of paclitaxel, improving the antitumor response. Mol Cancer Ther. 2016; 15(1):125–35.

    Article  CAS  PubMed  Google Scholar 

  7. Kabbinavar F, Hurwitz HI, Fehrenbacher L, Meropol NJ, Novotny WF, Lieberman G, Griffing S, Bergsland E. Phase ii, randomized trial comparing bevacizumab plus fluorouracil (fu)/leucovorin (lv) with fu/lv alone in patients with metastatic colorectal cancer. J Clin Oncol. 2003; 21(1):60–5.

    Article  CAS  PubMed  Google Scholar 

  8. Wildiers H, Guetens G, De Boeck G, Verbeken E, Landuyt B, Landuyt W, De Bruijn E, Van Oosterom A. Effect of antivascular endothelial growth factor treatment on the intratumoral uptake of cpt-11. Br J Cancer. 2003; 88(12):1979.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Liu Y, Suzuki M, Masunaga S-i, Chen Y-W, Kashino G, Tanaka H, Sakurai Y, Kirihata M, Ono K. Effect of bevacizumab treatment on p-boronophenylalanine distribution in murine tumor. J Radiat Res. 2012; 54(2):260–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  10. Van der Veldt AA, Lubberink M, Bahce I, Walraven M, de Boer MP, Greuter HN, Hendrikse NH, Eriksson J, Windhorst AD, Postmus PE, et al.Rapid decrease in delivery of chemotherapy to tumors after anti-vegf therapy: implications for scheduling of anti-angiogenic drugs. Cancer Cell. 2012; 21(1):82–91.

    Article  CAS  PubMed  Google Scholar 

  11. Maitland ML, Hudoba C, Snider KL, Ratain MJ. Analysis of the yield of phase ii combination therapy trials in medical oncology. Clin Cancer Res. 2010; 16(21):5296–302.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Sharma MR, Stadler WM, Ratain MJ. Randomized phase ii trials: a long-term investment with promising returns. J Natl Cancer Inst. 2011; 103(14):1093–100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Paller CJ, Bradbury PA, Ivy SP, Seymour L, LoRusso PM, Baker L, Rubinstein L, Huang E, Collyar D, Groshen S, et al.Design of phase i combination trials: recommendations of the clinical trial design task force of the nci investigational drug steering committee. Clin Cancer Res. 2014; 20(16):4210–7.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Lai X, Friedman A. Combination therapy of cancer with cancer vaccine and immune checkpoint inhibitor: A mathematical model. PLoS ONE. 2017; 12(2):0178479.

    Google Scholar 

  15. Lai X, Friedman A. Combination therapy of cancer with braf inhibitor and immune checkpoint inhibitor: A mathematical model. BMC Syst Biol. 2017; 11(70):1–18.

    Google Scholar 

  16. Ohm JE, Gabrilovich DI, Sempowski GD, Kisseleva E, Parman KS, Nadaf S, Carbone DP. Vegf inhibits t-cell development and may contribute to tumor-induced immune suppression. Blood. 2003; 101(12):4878–86.

    Article  CAS  PubMed  Google Scholar 

  17. Mulligan JK, Rosenzweig SA, Young MRI. Tumor secretion of vegf induces endothelial cells to suppress t cell functions through the production of pge2. J Immunother (Hagerstown, Md.: 1997). 2010; 33(2):126–35.

    Article  CAS  Google Scholar 

  18. Gavalas NG, Tsiatas M, Tsitsilonis O, Politi E, Ioannou K, Ziogas AC, Rodolakis A, Vlahos G, Thomakos N, Haidopoulos D, et al.Vegf directly suppresses activation of t cells from ascites secondary to ovarian cancer via vegf receptor type 2. Br J Cancer. 2012; 107(11):1869–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Ziogas AC, Gavalas NG, Tsiatas M, Tsitsilonis O, Politi E, Terpos E, Rodolakis A, Vlahos G, Thomakos N, Haidopoulos D, et al.Vegf directly suppresses activation of t cells from ovarian cancer patients and healthy individuals via vegf receptor type 2. Int J Cancer. 2012; 130(4):857–64.

    Article  CAS  PubMed  Google Scholar 

  20. Li Y-L, Zhao H, Ren X-B. Relationship of vegf/vegfr with immune and cancer cells: staggering or forward?. Cancer Biol Med. 2016; 13(2):206–14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Wallin JJ, Bendell JC, Funke R, Sznol M, Korski K, Jones S, Hernandez G, Mier J, He X, Hodi FS, et al.Atezolizumab in combination with bevacizumab enhances antigen-specific t-cell migration in metastatic renal cell carcinoma. Nat Commun. 2016; 7:12624.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Kuusk T, Albiges L, Escudier B, Grivas N, Haanen J, Powles T, Bex A. Antiangiogenic therapy combined with immune checkpoint blockade in renal cancer. Angiogenesis. 2017; 20(2):205–215.

    Article  CAS  PubMed  Google Scholar 

  23. Ott PA, Hodi FS, Buchbinder EI. Inhibition of immune checkpoints and vascular endothelial growth factor as combination therapy for metastatic melanoma: an overview of rationale, preclinical evidence, and initial clinical data. Front Oncol. 2015; 5(202):1–7.

    Google Scholar 

  24. Einstein DJ, McDermott DF. Combined blockade of vascular endothelial growth factor and programmed death 1 pathways in advanced kidney cancer. Clin Adv Hematol Oncol. 2017; 15(6):478–88.

    PubMed  Google Scholar 

  25. Alfaro C, Suarez N, Gonzalez A, Solano S, Erro L, Dubrot J, Palazon A, Hervas-Stubbs S, Gurpide A, Lopez-Picazo JM, et al.Influence of bevacizumab, sunitinib and sorafenib as single agents or in combination on the inhibitory effects of vegf on human dendritic cell differentiation from monocytes. Br J Cancer. 2009; 100(7):1111–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Gabrilovich DI, Chen HL, Girgis KR, Cunningham HT, Meny GM, Nadaf S, Kavanaugh D, Carbone DP. Production of vascular endothelial growth factor by human tumors inhibits the functional maturation of dendritic cells. Nat Med. 1996; 2(10):1096–103.

    Article  CAS  PubMed  Google Scholar 

  27. Voron T, Colussi O, Marcheteau E, Pernot S, Nizard M, Pointet A-L, Latreche S, Bergaya S, Benhamouda N, Tanchot C, et al.Vegf-a modulates expression of inhibitory checkpoints on cd8+ t cells in tumors. J Exp Med. 2015; 212(2):139–48.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Terme M, Pernot S, Marcheteau E, Sandoval F, Benhamouda N, Colussi O, Dubreuil O, Carpentier AF, Tartour E, Taieb J. Vegfa-vegfr pathway blockade inhibits tumor-induced regulatory t-cell proliferation in colorectal cancer. Cancer Res. 2013; 73(2):539–49.

    Article  CAS  PubMed  Google Scholar 

  29. Yasuda S, Sho M, Yamato I, Yoshiji H, Wakatsuki K, Nishiwada S, Yagita H, Nakajima Y. Simultaneous blockade of programmed death 1 and vascular endothelial growth factor receptor 2 (vegfr2) induces synergistic anti-tumour effect in vivo. Clin Exp Immunol. 2013; 172(3):500–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Sims GP, Rowe DC, Rietdijk ST, Herbst R, Coyle AJ. Hmgb1 and rage in inflammation and cancer. Annu Rev Immunol. 2010; 28:367–88.

    Article  CAS  PubMed  Google Scholar 

  31. Zandarashvili L, Sahu D, Lee K, Lee YS, Singh P, Rajarathnam K, et al.Real-time kinetics of high-mobility group box 1 (hmgb1) oxidation in extracellular fluids studied by in situ protein nmr spectroscopy. J Biol Chem. 2013; 288(17):11621–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Palucka J, Banchereau J. Cancer immunotherapy via dendritic cells. Nat Rev Cancer. 2012; 12(4):265–77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Saenz R, Futalan D, Leutenez L, Eekhout F, Fecteau JF, Sundelius S, et al.Tlr4-dependent activation of dendritic cells by an hmgb1-derived peptide adjuvant. J Transl Med. 2014; 12(211):1–11.

    Google Scholar 

  34. Janco JMT, Lamichhane P, Karyampudi L, Knutson KL. Tumor-infiltrating dendritic cells in cancer pathogenesis. J Immunol. 2015; 194(7):2985–91.

    Article  CAS  PubMed Central  Google Scholar 

  35. Ma Y, Shurin1 GV, Peiyuan Z, Shurin MR. Dendritic cells in the cancer microenvironment. J Cancer. 2013; 4(1):36–44.

  36. Perrot CY, Javelaud D, Mauviel A. Insights into the transforming growth factor-beta signaling pathway in cutaneous melanoma. Ann Dermatol. 2013; 25(2):135–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Whiteside TL. The role of regulatory t cells in cancer immunology. Immunotargets Ther. 2015; 4:159–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Shi L, Chen S, Yang L, Li Y. The role of pd-1 and pd-l1 in t-cell immune suppression in patients with hematological malignancies. J Hematol Oncol. 2013; 6(74):10–118617568722674.

    Google Scholar 

  39. Muppidi MR, George S. Immune checkpoint inhibitors in renal cell carcinoma. J Target Ther Cancer 2015. 2015; 4:47–52.

    Google Scholar 

  40. Umansky V, Blattner C, Gebhardt C, Utikal J. The role of myeloid-derived suppressor cells (mdsc) in cancer progression. Vaccines. 2016; 4(36):1–16.

    Google Scholar 

  41. Szomolay B, Eubank TD, Roberts RD, Marsh CB, Friedman A. Modeling the inhibition of breast cancer growth by gm-csf. J Theor Biol. 2012; 303:141–151.

    Article  CAS  PubMed  Google Scholar 

  42. Chanmee T, Ontong P, Konno K, Itano N. Tumor-associated macrophages as major players in the tumor microenvironment. Cancers. 2014; 6(3):1670–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Hao W, Friedman A. Serum upar as biomarker in breast cancer recurrence: A mathematical model. PLoS ONE. 2016; 11(4):0153508.

    Google Scholar 

  44. Chen D, Roda JM, Marsh CB, Eubank TD, Friedman A. Hypoxia inducible factors-mediated inhibition of cancer by gm-csf: a mathematical model. Bull Math Biol. 2012; 74(11):2752–77.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. Cheng X, Veverka V, Radhakrishnan A, Waters LC, Muskett FW, Morgan SH, et al.Human pd-l1/b7-h1/cd274 protein. Sino Biological Inc.

  46. Francisco LM, Sage PT, Sharpe AH. The pd-1 pathway in tolerance and autoimmunity. Immunol Rev. 2010; 236:219–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Mautea RL, Gordona SR, Mayere AT, McCrackena MN, Natarajane A, Ring NG, et al.Engineering high-affinity pd-1 variants for optimized immunotherapy and immuno-pet imaging. Proc Natl Acad Sci USA. 2015; 112(47):6506–14.

    Article  CAS  Google Scholar 

  48. D’Acunto B. Computational Methods for PDE in Mechanics. Series on Advances in Mathematics for Applied Sciences-Vol.67. Singapore: Word Scientific; 2004.

    Google Scholar 

  49. Marino S, Hogue I, Ray C, Kirschner D. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol. 2008; 254(1):178–96.

    Article  PubMed  PubMed Central  Google Scholar 

  50. W. Humphrey R, M Brockway-Lunardi L, T Bonk D, Dohoney KM, Doroshow JH, Meech SJ, Ratain MJ, Topalian SL, M Pardoll D. Opportunities and challenges in the development of experimental drug combinations for cancer. J Natl Cancer Inst. 2011; 103(16):1222–6.

    Article  PubMed Central  CAS  Google Scholar 

  51. Miles D, von Minckwitz G, Seidman AD. Combination versus sequential single-agent therapy in metastatic breast cancer. Oncologist. 2002; 7(Supplement 6):13–9.

    CAS  PubMed  Google Scholar 

  52. Fujimoto-Ouchi K, Tanaka Y, Tominaga T. Schedule dependency of antitumor activity in combination therapy with capecitabine/5?-deoxy-5-fluorouridine and docetaxel in breast cancer models. Clin Cancer Res. 2001; 7(4):1079–86.

    CAS  PubMed  Google Scholar 

  53. Frei III E, Eder JP. Principles of dose, schedule, and combination therapy. Cancer Med. 2003; 7:590–9.

    Google Scholar 

  54. Cabibbo G, Tremosini S, Galati G, Mazza G, Gadaleta-Caldarola G, Lombardi G, Antonucci M, Sacco R. Transarterial chemoembolization and sorafenib in hepatocellular carcinoma. Expert Rev Anticancer Ther. 2014; 14(7):831–45.

    Article  CAS  PubMed  Google Scholar 

  55. Young ME. Estimation of diffusion coefficients of proteins. Biotech Bioeng. 1980; XXII:947–55.

    Article  Google Scholar 

  56. Shui YB, Wang X, Hu JS, Wang SP, Garcia CM, et al.Vascular endothelial growth factor expression and signaling in the lens. Invest Ophthalmol Vis Sci. 2003; 44(9):3911–9.

    Article  PubMed  Google Scholar 

  57. Liao KL, Bai XF, Friedman A. Mathematical modeling of interleukin-27 induction of anti-tumor t cells response. PLoS ONE. 2014; 9(3):91844.

    Article  CAS  Google Scholar 

  58. Androjna C, Gatica JE, Belovich JM, Derwin KA. Oxygen diffusion through natural extracellular matrices: implications for estimating critical thickness values in tendon tissue engineering. Tissue Eng A. 2008; 14(4):559–69.

    Article  CAS  Google Scholar 

  59. Friedman A, Hao W. The role of exosomes in pancreatic cancer microenvironment. Bull Math Biol. 2017; 80(5):1111–113.

    Article  PubMed  CAS  Google Scholar 

  60. Chen D, Bobko AA, Gross AC, Evans R, Marsh CB, Khramtsov VV, et al.Involvement of tumor macrophage hifs in chemotherapy effectiveness: mathematical modeling of oxygen, ph, and glutathione. PLoS ONE. 2014; 9(10):107511.

    Article  CAS  Google Scholar 

  61. Lo W. -C, Arsenescu V, Arsenescu RI, Friedman A. Inflammatory bowel disease: How effective is tnf-alpha suppression?. PLoS One. 2016; 11(11):0165782.

    Google Scholar 

  62. Abcam/. Anti-pd1 antibody (ab89828).

  63. Agrawal S, Joshi M, Christoforidis JB. Vitreous inflammation associated with intravitreal anti-vegf pharmacotherapy. Mediat Inflamm. 2013; 2013(943409):1–6.

    Article  CAS  Google Scholar 

  64. Kim Y, Lawler S, Nowicki MO, Chiocca EA, Friedman A. A mathematical model for pattern formation of glioma cells outside the tumor spheroid core. J Theor Biol. 2009; 260(3):359–71.

    Article  PubMed  Google Scholar 

  65. Kim Y, Wallace J, Li F, Ostrowski M, Friedman A. Transformed epithelial cells and fibroblasts/myofibroblasts interaction in breast tumor: a mathematical model and experiments. J Theor Biol. 2010; 61(3):401–21.

    Google Scholar 

Download references


This work is supported by the Mathematical Biosciences Institute and the National Science Foundation (DMS 0931642), and by the Renmin University of China and the National Natural Science Foundation of China (Grant No. 11501568), and the International Postdoctoral Exchange Fellowship Program 2016 by the Office of China Postdoctoral Council.


XL was supported by National Natural Science Foundation of China (Grant No. 11501568).

Availability of data and materials

The dataset supporting the conclusions of this article is included within the article.

Author information

Authors and Affiliations



AF and XL developed and simulated the model, and wrote the final manuscript. Both authors read and approved the final manuscript.

Corresponding author

Correspondence to Avner Friedman.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lai, X., Friedman, A. How to schedule VEGF and PD-1 inhibitors in combination cancer therapy?. BMC Syst Biol 13, 30 (2019).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: