 Research article
 Open Access
 Published:
How to schedule VEGF and PD1 inhibitors in combination cancer therapy?
BMC Systems Biology volume 13, Article number: 30 (2019)
Abstract
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 antiVEGF, which is known to affect the perfusion of other drugs. As a second drug we take antiPD1. Both drugs are known to increase the activation of anticancer T cells. Our simulations show that in the case where antiVEGF reduces the perfusion, a nonoverlapping schedule is significantly more effective than a simultaneous injection of the two drugs, and it is somewhat more beneficial to inject antiPD1 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
Antivascular endothelial growth factor (antiVEGF) is a drug commonly used as anticancer agent, although numerous studies show only modest results [1]. In combination with chemotherapy antiVEGF 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 antiVEGF 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 nonsmall cell lung cancer (NSCLC) after antiVEGF therapy, highlighting the importance of drug scheduling in combination therapy when antiVEGF 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 antiVEGF 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 anticancer T cells [16–20]; hence VEGF inhibition will enhance T cells function, and checkpoint blockade could therefore significantly advance antitumor therapy. The anticancer synergy between antiVEGF and checkpoint inhibitors is currently being evaluated in clinical trials in renal cancer [21, 22].
In the case where antiVEGF therapy decreases the perfusion of a second antitumor agent, the following question arises [23, 24]: Should the treatment with combination of antiVEGF 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 nonoverlappingly?
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 IL12 and IL2. The network of interactions among these species is shown in Fig. 1. This figure includes also oxygen concentration, and programmed cell death protein 1 (PD1) and its ligand PDL1. 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–20]; VEGF also enhances the expression of PD1 on CD 8^{+} 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 antiVEGF decreases the perfusion of antiPD1 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 antiVEGF is delayed or advanced relative to antiPD1, 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 nonoverlappingly.
We finally consider briefly the case where antiVEGF increases the perfusion of antiPD1, and show that in this case there are more benefits when the two drugs are given simultaneously.
Methods
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 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 antitumor 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 (HMGB1) [30]. We model the dynamics of N_{C} and HMGB1 (H) by the following equations:
where \(\lambda _{N_{C}C}\) is the rate of cancer cells becoming necrotic and \(\lambda _{HN_{C}}\) is the production rate of HMGB1 by necrotic cells. We note that since molecules like HMGB1, 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 HMGB1 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}Cd_{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 HMGB1 [32, 33]. Hence, the activation rate of immature dendritic cells, with density D_{0}, is proportional to \(D_{0}\frac {H}{K_{H}+H}\), or to \(D_{0}\frac {C}{K_{C}+C}\). Here, the MichaelisMenten 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].
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 (T_{1}) under IL12 (I_{12}) environment [34, 35], and this process is inhibited by Tregs [36, 37] and by VEGF [16–20]. The proliferation of activated T_{1} cells is enhanced by IL2. Both processes of activation and proliferation of T_{1} are assumed to be inhibited by PD1PDL1 (Q) [38, 39]; we represent this inhibition by the factor \(\frac {1}{1+Q/K_{TQ}}\). Hence T_{1} satisfies the following equation:
where T_{10} is the density of the naive CD 4^{+} T cells.
Similarly, inactive CD 8^{+} T cells are activated by IL12 [34, 35], a process resisted by Tregs [36, 37] and VEGF [16–20], while IL2 enhances the proliferation of activated CD 8^{+} T cells. Hence,
where T_{80} is the density of the inactive CD 8^{+} T cells.
Equation for Tregs (T _{r})
Naive CD 4^{+} 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
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: 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 (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 IL12 (I _{12})
The proinflammatory cytokine IL12 is secreted by activated DCs [34, 35]; hence it satisfies the equation:
Equation for IL2 (I _{2})
IL2 is produced by activated CD 4^{+} T cells (T_{1}) [35]. Hence,
Equation for TGFβ: (T _{β})
TGFβ is produced by tumor cells [36] and Tregs [37], so that
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:
where d_{W}W represents the takeup rate of oxygen by all the cells.
Equation for VEGF (G)
VEGF is produced by cancer cells [41, 42] and is depleted by antiVEGF. Hence G satisfies the following equation:
where B is the effective antiVEGF 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 PD1 (P), PDL1 (L) and PD1PDL1 (Q)
PD1 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 PD1 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 PD1 on T_{8} cells by a factor ε_{G}G [27]. If we denote by ρ_{P} the ratio between the mass of one PD1 protein to the mass of one T cell, then the total concentration of PD1 on T_{1} and T_{r} cells is given by
and the concentration of PD1 on T_{8} cells is given by
PDL1 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 PDL1 per cell is the same for T_{1} and T_{8} cells, and denote the ratio between the mass of one PDL1 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 \(\frac {\partial T}{\partial t}\), there corresponds a change of P_{1}, given by \(\rho _{P}\frac {\partial T}{\partial t}\). For the same reason, ∇·(uP_{1})=ρ_{P}∇·(uT) and ∇^{2}P_{1}=ρ_{P}∇^{2}T when no antiPD1 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=righthand side. Similarly,
When only antiPD1 drug (A) is injected, PD1 is depleted at a rate μ_{PA}P_{1}A, but when also antiVEGF (B) is injected the depletion of PD1 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 P_{1} satisfies the following equation:
Similarly,
PDL1 combines with PD1 on the plasma membrane of T cells, to form a complex PD1PDL1 (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 halflife 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 antiPD1 (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 PD1. Hence,
where I_{A}(t)=1 during dosing, and I_{A}(t)=0 otherwise.
Equation for antiVEGF (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 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, T_{1} and T_{8} are activated by IL12 and T_{r} is activated by T_{β}. We then have the following flux conditions at the tumor’s boundary:
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
The boundary condition for the oxygen is given by
where W_{0} is the normal level of oxygen concentration in the blood. We impose noflux boundary condition for all the remaining variables:
It is tacitly assumed here that the receptors PD1 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.
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 IL12 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 antiVEGF 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 antiPD1 as a single agent reduces the volume of tumor more than antiVEGF 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 antiVEGF as a single agent reduces the tumor volume more than the antiPD1 as single agent, in agreement with Fig. 4j in [27]. In Fig. 5 we see that (with γ_{A}=0) antiVEGF with γ_{B}=2×10^{−8}g/cm^{3}·day reduces significantly the PD1 expression on CD 8^{+} T cells, in agreement with experimental results in [27].
Clinical trials in silico
In clinical trials the treatment and followup 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 antiPD1 and antiVEGF are both administered continuously in the first 3 weeks of the 9week cycle. In schedule S2 antiPD1 is given continuously in the first 3 weeks, followed by antiVEGF in the next 3 weeks, with no drugs for the remaining 3 weeks of the cycle. Schedule S3 is similar to schedule S2 with antiVEGF in the first 3 weeks and antiPD1 in the next 3 weeks. We refer to schedules S2 and S3 as nonoverlapping 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 nonoverlapping schedules S2 and S3 reduce the tumor volume significantly faster than schedule S1, and schedule S2 is somewhat more effective than schedule S3.
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^{−10}g/cm^{3}·day,γ_{B}=9.5×10^{−9}g/cm^{3}·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.
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 antiVEGF, γ_{B}, in units of g/cm^{3}·day, and the vertical axis scales the dose amount of antiPD1, γ_{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.
We conclude that nonoverlapping treatment is far more beneficial than simultaneous treatment, and there is some advantage in injecting antiPD1 before antiVEGF.
So far we considered the case where antiVEGF blocks perfusion. We next consider briefly the case where antiVEGF 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 \(\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 nonoverlapping 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 \(\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 pvalues 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 \(\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 K_{DG},K_{TQ} and K_{TG}, 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 antiPD1 drug by antiVEGF 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 CD 8^{+} T cells) has the largest impact in reducing tumor volume, and K_{TQ} (the parameter which increases the inhibition of CD 8^{+} T cells by PD1PDL1) 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 interactions between the diverse agents, addressing, for example the following questions:

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.
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.
If the two drugs are injected intermittently, should they be injected simultaneously or nonoverlappingly, 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 (antiPD1), and the second drug is a cancer vaccine [14] or a BRAFinhibitor [15]. In the present paper we addressed the question (iii) when the two drugs are antiPD1 and antiVEGF. AntiVEGF 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 antiVEGF decreases the perfusion of the antiPD1, the time it takes to reduce tumor volume by 95% is much shorter when the injections of the two drugs are nonoverlapping than when the injections are given at the same time. Furthermore, in the nonoverlapping treatment, if we inject the antiPD1 first we get the 95% reduction somewhat faster than if we inject antiVEGF first.
On the other hand, in the case when antiVEGF increases the perfusion of antiPD1, a treatment with simultaneous infections is more beneficial than a treatment with nonoverlapping injections.
The optimal scheduling in combination therapy for cancer was considered in several publications, both in terms of toxicity [51] and efficacy [52–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 antiVEGF, 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 antiVEGF and antiPD1. AntiVEGF increases the activation of anticancer cells, but it also modifies the perfusion of other drugs. Using a mathematical model we showed that a nonoverlappling adminstration of the two drugs is more effective in reducing tumor volume than simultaneous administration of the drugs in the case where antiVEGF degrades perfusion, as it occurs in some cancers, while the opposite is true when perfusion is enhanced by antiVEGF. The mathematical methodology developed in this paper could be extended to treatment with other combinations of drugs.
Appendix
Parameter estimation
Halfsaturation
In an expression of the form \(Y\frac {X}{K_{X}+X}\) where Y is activated by X, the parameter K_{X} is called the halfsaturation of X.
We assume that if the average density (or concentration) of X,
converges to a steady state X_{0}, then
is not “too small” and not “too close” to 1, and for definiteness we take
so that
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]:
where M_{G} and δ_{G} are respectively the molecular weight and diffusion coefficient of VEGF, M_{p} is the molecular weight of p, and M_{G}=24kDa [56] and δ_{G}=8.64×10^{−2}cm^{2} day^{−1} [57]. Since, M_{B}=149kDa (bevacizumab), we get δ_{B}=4.70×10^{−2}cm^{2} day^{−1}. The diffusion coefficient of oxygen in the extracellular matrix (ECM) in the range of 7×10^{−6}−2×10^{−5}cm^{2}/s [58]; we take it to be δ_{W}=0.8cm^{2}/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], d_{D}=0.1/day, C=K_{C}=0.4g/cm^{3},D=K_{D}=4×10^{−4}g/cm^{3},D_{0}=2×10^{−5}g/cm^{3}. We assume that K_{DG}=4K_{G} where K_{G}=7×10^{−8}g/cm^{3} [43]; hence λ_{DC}=2.5d_{D}D/D_{0}=5/day. For simplicity we had assumed that the source of inactive dendritic cells, D_{0}, 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/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], \(\lambda _{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\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
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 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 \(\lambda _{T_{1}I_{12}}\) and \(\lambda _{T_{8}I_{12}}\). The same considerations apply in the case of T_{r}.
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], d_{E}=0.69/day, E_{M}=5×10^{−3}g/cm^{3},K_{E}=2.5×10^{−3}g/cm^{3},G_{0}=3.65×10^{−10}. We take E_{M}=2K_{E}=5×10^{−3}g/cm^{3}. From the steady state of Eq. (6), we get λ_{E}=2d_{E}/(K_{G}−G_{0})=1.98×10^{7}cm^{3}/g·day. Assuming that the threshold G_{0} (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×10^{7}cm^{3}/g·day in mice and λ_{E}=2.08×10^{7}cm^{3}/g·day in humans.
Eq. (7).
We take d_{C}=0.17day^{−1},C_{M}=0.8g/cm^{3} [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/cm^{3}, and W=W_{0}=K_{W}=1.69×10^{−4}g/cm^{3} (see the estimates for Eq. (11)). From the steady state of Eq. (7) in the control case we have,
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.19cm^{3}/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, λ_{WE}E−d_{W}W=0, where λ_{WE}=7×10^{−2}/day [43], W=K_{W}=1.69×10^{−4}g/cm^{3} [43], E=K_{E}=2.5×10^{−3}g/cm^{3}. Hence, d_{W}=λ_{WE}E/W=1.04/day.
Eq. (12).
From steady state of Eq. (12) we get, λ_{GW}C−d_{G}G=0, where d_{G}=12.6/day [43], G=K_{G}=7×10^{−8}g/cm^{3} [43], C=K_{C}=0.4g/cm^{3}. 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/K_{G}=1.43×10^{6}cm^{3}/g. From Eqs. (13)(16) we get,
and
In steady state with \(\phantom {\dot {i}\!}P=K_{P}=K_{P_{1}}+K_{P_{8}},L=K_{L},Q=K_{Q}\) and G=K_{G} we have, by Eq. (17), K_{Q}=σK_{P}K_{L}. 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
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 d_{A}=0.34day^{−1} and assume that 10% of A is used in blocking PD1, while the remaining 90% degrades naturally. Hence,
Since the molecular mass of antiPD1 (32kDa [62]) approximates the molecular mass of PD1 (20.5  40 kDa [62]), we take μ_{PA}=μ_{AP}=4.33×10^{7} cm^{3}/g·day.
By [63], the halflife of antiVEGF is 2.824.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,
The molecular mass of antiVEGF 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_{i}^{k}\) and \(C_{i}^{k}\) denote numerical approximations of ith grid point and \(C(r_{i}^{k},n\tau)\), respectively, where τ is the size of the timestep. The discretization of Eq. (29) is derived by the fully implicit finite difference scheme:
where \(C_{r}=\frac {h_{1}^{2}C_{i+1}^{k+1}h_{1}^{2}C_{i1}^{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_{i1}^{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_{i1}^{k+1}(h_{1}^{2}h_{1}^{2})u_{i}^{k+1}}{h_{1}(h_{1}^{2}h_{1}h_{1})},h_{1}=r_{i1}^{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 nondimensionalized the model and then applied the above method.
References
 1
Meadows KL, Hurwitz HI. Antivegf therapies in the clinic. Cold Spring Harb Perspect Med. 2012; 2(10):006577.
 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.
 3
Turley RS, Fontanella AN, Padussis JC, Toshimitsu H, Tokuhisa Y, Cho EH, Hanna G, Beasley GM, Augustine CK, Dewhirst MW, et al.Bevacizumabinduced alterations in vascular permeability and drug delivery: a novel approach to augment regional chemotherapy for intransit melanoma. Clin Cancer Res. 2012; 18(12):3328–39.
 4
Pastuskovas CV, Mundo EE, Williams SP, Nayak TK, Ho J, Ulufatu S, Clark S, Ross S, Cheng E, ParsonsReponte K, et al.Effects of antivegf on pharmacokinetics, biodistribution, and tumor penetration of trastuzumab in a preclinical breast cancer model. Mol Cancer Ther. 2012; 11(3):752–62.
 5
DaldrupLink HE, Okuhata Y, Wolfe A, Srivastav S, Øie S, Ferrara N, Cohen RL, Shames DM, Brasch RC. Decrease in tumor apparent permeabilitysurface area product to a mri macromolecular contrast medium following angiogenesis inhibition with correlations to cytotoxic drug accumulation. Microcirculation. 2004; 11(5):387–96.
 6
Cesca M, Morosi L, Berndt A, Nerini IF, Frapolli R, Richter P, Decio A, Dirsch O, Micotti E, Giordano S, et al.Bevacizumabinduced inhibition of angiogenesis promotes a more homogeneous intratumoral distribution of paclitaxel, improving the antitumor response. Mol Cancer Ther. 2016; 15(1):125–35.
 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.
 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 cpt11. Br J Cancer. 2003; 88(12):1979.
 9
Liu Y, Suzuki M, Masunaga Si, Chen YW, Kashino G, Tanaka H, Sakurai Y, Kirihata M, Ono K. Effect of bevacizumab treatment on pboronophenylalanine distribution in murine tumor. J Radiat Res. 2012; 54(2):260–7.
 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 antivegf therapy: implications for scheduling of antiangiogenic drugs. Cancer Cell. 2012; 21(1):82–91.
 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.
 12
Sharma MR, Stadler WM, Ratain MJ. Randomized phase ii trials: a longterm investment with promising returns. J Natl Cancer Inst. 2011; 103(14):1093–100.
 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.
 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.
 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. https://doi.org/10.1186/s1291801704469.
 16
Ohm JE, Gabrilovich DI, Sempowski GD, Kisseleva E, Parman KS, Nadaf S, Carbone DP. Vegf inhibits tcell development and may contribute to tumorinduced immune suppression. Blood. 2003; 101(12):4878–86.
 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.
 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.
 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.
 20
Li YL, Zhao H, Ren XB. Relationship of vegf/vegfr with immune and cancer cells: staggering or forward?. Cancer Biol Med. 2016; 13(2):206–14.
 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 antigenspecific tcell migration in metastatic renal cell carcinoma. Nat Commun. 2016; 7:12624.
 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.
 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.
 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.
 25
Alfaro C, Suarez N, Gonzalez A, Solano S, Erro L, Dubrot J, Palazon A, HervasStubbs S, Gurpide A, LopezPicazo 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.
 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.
 27
Voron T, Colussi O, Marcheteau E, Pernot S, Nizard M, Pointet AL, Latreche S, Bergaya S, Benhamouda N, Tanchot C, et al.Vegfa modulates expression of inhibitory checkpoints on cd8+ t cells in tumors. J Exp Med. 2015; 212(2):139–48.
 28
Terme M, Pernot S, Marcheteau E, Sandoval F, Benhamouda N, Colussi O, Dubreuil O, Carpentier AF, Tartour E, Taieb J. Vegfavegfr pathway blockade inhibits tumorinduced regulatory tcell proliferation in colorectal cancer. Cancer Res. 2013; 73(2):539–49.
 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 antitumour effect in vivo. Clin Exp Immunol. 2013; 172(3):500–6.
 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.
 31
Zandarashvili L, Sahu D, Lee K, Lee YS, Singh P, Rajarathnam K, et al.Realtime kinetics of highmobility group box 1 (hmgb1) oxidation in extracellular fluids studied by in situ protein nmr spectroscopy. J Biol Chem. 2013; 288(17):11621–7.
 32
Palucka J, Banchereau J. Cancer immunotherapy via dendritic cells. Nat Rev Cancer. 2012; 12(4):265–77.
 33
Saenz R, Futalan D, Leutenez L, Eekhout F, Fecteau JF, Sundelius S, et al.Tlr4dependent activation of dendritic cells by an hmgb1derived peptide adjuvant. J Transl Med. 2014; 12(211):1–11.
 34
Janco JMT, Lamichhane P, Karyampudi L, Knutson KL. Tumorinfiltrating dendritic cells in cancer pathogenesis. J Immunol. 2015; 194(7):2985–91.
 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 factorbeta signaling pathway in cutaneous melanoma. Ann Dermatol. 2013; 25(2):135–44.
 37
Whiteside TL. The role of regulatory t cells in cancer immunology. Immunotargets Ther. 2015; 4:159–71.
 38
Shi L, Chen S, Yang L, Li Y. The role of pd1 and pdl1 in tcell immune suppression in patients with hematological malignancies. J Hematol Oncol. 2013; 6(74):10–118617568722674.
 39
Muppidi MR, George S. Immune checkpoint inhibitors in renal cell carcinoma. J Target Ther Cancer 2015. 2015; 4:47–52.
 40
Umansky V, Blattner C, Gebhardt C, Utikal J. The role of myeloidderived suppressor cells (mdsc) in cancer progression. Vaccines. 2016; 4(36):1–16.
 41
Szomolay B, Eubank TD, Roberts RD, Marsh CB, Friedman A. Modeling the inhibition of breast cancer growth by gmcsf. J Theor Biol. 2012; 303:141–151.
 42
Chanmee T, Ontong P, Konno K, Itano N. Tumorassociated macrophages as major players in the tumor microenvironment. Cancers. 2014; 6(3):1670–90.
 43
Hao W, Friedman A. Serum upar as biomarker in breast cancer recurrence: A mathematical model. PLoS ONE. 2016; 11(4):0153508.
 44
Chen D, Roda JM, Marsh CB, Eubank TD, Friedman A. Hypoxia inducible factorsmediated inhibition of cancer by gmcsf: a mathematical model. Bull Math Biol. 2012; 74(11):2752–77.
 45
Cheng X, Veverka V, Radhakrishnan A, Waters LC, Muskett FW, Morgan SH, et al.Human pdl1/b7h1/cd274 protein. Sino Biological Inc. http://www.sinobiological.com/PDL1B7H1CD274Proteing533.html.
 46
Francisco LM, Sage PT, Sharpe AH. The pd1 pathway in tolerance and autoimmunity. Immunol Rev. 2010; 236:219–42.
 47
Mautea RL, Gordona SR, Mayere AT, McCrackena MN, Natarajane A, Ring NG, et al.Engineering highaffinity pd1 variants for optimized immunotherapy and immunopet imaging. Proc Natl Acad Sci USA. 2015; 112(47):6506–14.
 48
D’Acunto B. Computational Methods for PDE in Mechanics. Series on Advances in Mathematics for Applied SciencesVol.67. Singapore: Word Scientific; 2004.
 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.
 50
W. Humphrey R, M BrockwayLunardi 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.
 51
Miles D, von Minckwitz G, Seidman AD. Combination versus sequential singleagent therapy in metastatic breast cancer. Oncologist. 2002; 7(Supplement 6):13–9.
 52
FujimotoOuchi K, Tanaka Y, Tominaga T. Schedule dependency of antitumor activity in combination therapy with capecitabine/5?deoxy5fluorouridine and docetaxel in breast cancer models. Clin Cancer Res. 2001; 7(4):1079–86.
 53
Frei III E, Eder JP. Principles of dose, schedule, and combination therapy. Cancer Med. 2003; 7:590–9.
 54
Cabibbo G, Tremosini S, Galati G, Mazza G, GadaletaCaldarola G, Lombardi G, Antonucci M, Sacco R. Transarterial chemoembolization and sorafenib in hepatocellular carcinoma. Expert Rev Anticancer Ther. 2014; 14(7):831–45.
 55
Young ME. Estimation of diffusion coefficients of proteins. Biotech Bioeng. 1980; XXII:947–55.
 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.
 57
Liao KL, Bai XF, Friedman A. Mathematical modeling of interleukin27 induction of antitumor t cells response. PLoS ONE. 2014; 9(3):91844.
 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.
 59
Friedman A, Hao W. The role of exosomes in pancreatic cancer microenvironment. Bull Math Biol. 2017; 80(5):1111–113.
 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.
 61
Lo W. C, Arsenescu V, Arsenescu RI, Friedman A. Inflammatory bowel disease: How effective is tnfalpha suppression?. PLoS One. 2016; 11(11):0165782.
 62
Abcam/. Antipd1 antibody (ab89828). http://www.abcam.com/pd1antibodyab89828.html.
 63
Agrawal S, Joshi M, Christoforidis JB. Vitreous inflammation associated with intravitreal antivegf pharmacotherapy. Mediat Inflamm. 2013; 2013(943409):1–6.
 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.
 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.
Acknowledgements
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.
Funding
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
Affiliations
Contributions
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 (http://creativecommons.org/licenses/by/4.0/), 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(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 AntiPD1
 AntiVEGF
 Combination therapy
 Scheduling
 PDE model