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

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

Background

Anti-vascular endothelial growth factor (anti-VEGF) is a drug commonly used as anticancer agent, although numerous studies show only modest results [1]. In combination with chemotherapy anti-VEGF improves 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 [16-20]; hence VEGF inhibition will enhance T cells function, and checkpoint blockade could therefore significantly advance antitumor therapy. The anticancer synergy between anti-VEGF and checkpoint inhibitors is currently being evaluated in clinical trials in renal cancer [21, 22].

In the case where anti-VEGF therapy decreases the perfusion of a second antitumor agent, the following question arises [23, 24]: Should the treatment with combination of anti-VEGF and a checkpoint inhibitor be given at the same time, or is it more beneficial to delay treatment of one of the two drugs, so that they are given non-overlappingly?

To address this question we develop a mathematical model using a system of partial differential equations (PDEs). The variables of the model include 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 [16-20]; VEGF also enhances the expression of PD-1 on CD 8+ T cells [27], and induces Treg proliferation [28].

Fig. 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.

Methods

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)

Notation	Description
D	density of dendritic cells
T 1	density of activated CD 4+ T cells
T 8	density of activated CD 8+ T cells
T r	density of activated Treg cells
E	density of endothelial cells
C	density of cancer cells
N C	density of necrotic cell
H	HMGB-1 concentration
I 12	IL-12 concentration
I 2	IL-2 concentration
T β	TGF- β concentration
W	oxygen concentration
G	VEGF concentration
P 1	concentration of PD-1 on CD 4+ T cells
P 8	concentration of PD-1 on CD 8+ T cells
L	PD-L1 concentration
Q	PD-1-PD-L1 concentration
A	anti-PD-L1 concentration
B	anti-VEGF concentration

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 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

Distribution of cells in space

Equation for DCs (D)

By necrotic cancer cells (N) we mean cancer cells undergoing the process of necrosis. Necrotic cancer cells release high mobility group box 1 protein (HMGB-1) [30]. We model the dynamics of N and HMGB-1 (H) by the following equations: where is the rate of cancer cells becoming necrotic and 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 and , 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 , or to . Here, the Michaelis-Menten law was used to account for the limited rate of receptor recycling which takes place in the process of DCs activation. The dynamics of DCs is given by the following equation:

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

Equations for CD 4+ T cells (T1) and CD 8+ T cells (T8)

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 [16-20]. 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 . Hence T1 satisfies the following equation:

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 [16-20], while IL-2 enhances the proliferation of activated CD 8+ T cells. Hence,

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

Equation for Tregs (T)

Naive CD 4+ T cells differentiate into T cells under the influence of TGF- β: (T) [37, 40] and VEGF [28]. The equation for T 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 is the carrying capacity of endothelial cells, λ(G)=λ(G−G0)+, 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 (C), 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 T1 and T8 cells, respectively. d is the natural death rate of cancer cells, and where W0 is the normal level of oxygen concentration in the blood.

Equation for IL-12 (I12)

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

Equation for IL-2 (I2)

IL-2 is produced by activated CD 4+ T cells (T1) [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 dW represents the take-up rate of oxygen by all the cells.

Equation for VEGF (G)

VEGF is produced by cancer cells [41, 42] and is depleted by anti-VEGF. Hence G satisfies the following equation:

where B is the effective anti-VEGF concentration in the tumor, and Here we assumed that the secretion rate of VEGF by cancer cells increases with the oxygen level, but falls off when oxygen level exceeds a certain level, W∗ [44].

Equations for PD-1 (P), PD-L1 (L) and PD-1-PD-L1 (Q)

PD-1 is expressed on the surface of activated 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 ε, for T cells. VEGF increases the PD-1 on T8 cells by a factor εG [27]. If we denote by ρ 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 T cells is given by

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

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 ρ. Then for some parameters ε,ε, where ε depends on the specific type of tumor.

To a change in T=T1+εT, given by , there corresponds a change of P1, given by . For the same reason, ∇·(uP1)=ρ∇·(uT) and ∇2P1=ρ∇2T when no anti-PD-1 drug is injected. Hence, P1 satisfies the following equation: Recalling Eqs. (3) and (5) for T1 and T, we get, where RHS=right-hand side. Similarly, When only anti-PD-1 drug (A) is injected, PD-1 is depleted at a rate μP1A, but when also anti-VEGF (B) is injected the depletion of PD-1 is reduced, due to restricted perfusion, by a factor ; 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:

Similarly,

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 α and d, respectively, we can write where P=P1+P8. The half-life of Q is less then 1 s (1.16×10−5day) [47], so that d is very large. Hence we may approximate the dynamical equation for Q by the steady state equation, αPL=dQ, so that

where σ=α/d.

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 γ throughout the tumor during the dosing period. We denote by μ the depletion rate of A caused by blocking PD-1. Hence,

where I(t)=1 during dosing, and I(t)=0 otherwise.

Equation for anti-VEGF (B)

We denote by γ the concentration of the injected drug B during the dosing period, and by μ the depletion rate of B blocking VEGF. The equation for B is then given by

where I(t)=1 during dosing, and I(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,T and C (in units of g/cm3) 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, where e 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 and at the tumor’s boundary, and, that upon crossing the tumor boundary, T1 and T8 are activated by IL-12 and T is activated by T. We then have the following flux conditions at the tumor’s boundary:

where we take and . We assume that the endothelial cells that are attracted by the VEGF into the tumor microenvironment have constant density at the tumor’s boundary, and take

The boundary condition for the oxygen is given by

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

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

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

Fig. 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

Table 3

Summary of parameter values

Notation	Description	Value used	References
δ D	Diffusion coefficient of DCs	8.64×10−7cm2day−1	[59]
δ T	Diffusion coefficient of T cells	8.64×10−7cm2day−1	[59]
δ E	Diffusion coefficient of endothelial cells	8.64×10−7cm2day−1	[59]
δ C	Diffusion coefficient of tumor cells	8.64×10−7cm2day−1	[59]
δ_I_12	Diffusion coefficient of IL-12	6.05×10−2cm2day−1	[15]
δ_I_2	Diffusion coefficient of IL-2	9.58×10−2cm2day−1	[15]
δ_T_β	Diffusion coefficient of TGF-β	8.52×10−2cm2day−1	[15]
δ W	Diffusion coefficient of oxygen	0.8 cm2day−1	Estimated
δ G	Diffusion coefficient of VEGF	8.64×10−2cm2 day−1	[57]
δ A	Diffusion coefficient of anti-PD-L1	4.73×10−2cm2day−1	[15]
δ B	Diffusion coefficient of anit-VGEF	4.70×10−2cm2day−1	Estimated
σ T	Flux rate of T1 and T8 cells at the boundary	1 cm−1	[59]
σ E	Flux rate of T1 and T8 cells at the boundary	1 cm−1	[59]
χ G	Chemoattraction coefficient of VEGF	10 cm5/g·day	[64, 65]
λ DC	Activation rate of DCs by tumor cells (mice)	17.5day	Estimated
λ DC	Activation rate of DCs by tumor cells (humans)	7.5day	Estimated
λ_T_1I_12	Activation rate of CD 4+ T cells by IL-12	11.65 day−1	Estimated
λ_T_1I_2	Activation rate of CD 4+ T cells by IL-2	0.25 day−1	[59]
λ_T_8I_12	Activation rate of CD 8+ T cells by IL-12	10.38 day−1	Estimated
λ_T_8I_2	Activation rate of CD 8+ T cells by IL-2	0.25 day−1	[59]
λ_T_rT_β	Activation rate of Tregs by TGF-β	0.415 day−1	Estimated
λ_T_rG	Activation rate of Tregs by VEGF	0.083 day−1	Estimated
λ E	Growth rate of endothelial cells (mice)	2.77×107cm3/g·day	Estimated
λ E	Growth rate of endothelial cells (humans)	2.08×107cm3/g·day	Estimated
λ CW	Growth rate of cancer cells (mice)	2.24 day−1	Estimated
λ CW	Growth rate of cancer cells (humans)	1.76 day−1	Estimated
λ_I_12D	Production rate of IL-12 by DCs	2.21×10−6day−1	[15]
λ_I_2T_1	Production rate of IL-2 by CD 4+ T cells	2.82×10−8day−1	[15]
λ_T_βC	Production rate of TGF-β by cancer cells	3.27×10−10day−1	Estimated
λ_T_βT_r	Production rate of TGF-β by Tregs	5.57×10−9day−1	[61]
λ WE	Production rate of oxygen by endothelial cells	7×10−2/day	[43]
λ GW	Production rate of VEGF by cancer cells	2.21×10−6day−1	Estimated
ε T	Relative expression of PD-1 in Tregs	0.8	[15]
ε C	Relative expression of PD-L1 in tumor cells	0.01	[15]
ε G	Relative rate of PD-1 promotion by VEGF	1.43×106cm3/g	Estimated

Table 4

Summary of parameter values

Notation	Description	Value used	References
d D	Death rate of DCs	0.1 day−1	[59]
d_T_1	Death rate of CD 4+ T cells	0.197day−1	[59]
d_T_8	Death rate of CD 8+ T cells	0.18day−1	[59]
d_T_r	Death rate of Tregs	0.2day−1	[61]
d E	Death rate of endothelial cells	0.69 day−1	[43]
d C	Death rate of tumor cells	0.17 day−1	[59]
d_I_12	Degradation rate of IL-12	1.38 day−1	[59]
d_I_2	Degradation rate of IL-2	2.376 day−1	[59]
d_T_β	Degradation rate of TGF-β	499.066 day−1	[15]
d W	Take-up rate of oxygen by cells	1.04 day−1	Estimated
d G	Degradation rate of VEGF	12.6 day−1	[43]
d A	Degradation rate of anti-PD-L1	0.34 day−1	Estimated
η 1	Killing rate of cancer cells by T1	30.19 cm3/g·day	Estimated
η 8	Killing rate of cancer cells by T8	60.375 cm3/g·day	Estimated
d B	Degradation rate of anti-VEGF	0.17 day−1	Estimated
μ PA	Blocking rate of PD-1 by anti-PD-1	4.33×107cm3/g·day	Estimated
μ AP	Degradation rate of anti-PD-1 in blocking PD-1	4.33×107cm3/g·day	Estimated
μ GB	Blocking rate of VEGF by anti-VEGF	2.19×107cm3/g·day	Estimated
μ BG	Degradation rate of anti-VEGF in blocking VEGF	1.31×108cm3/g·day	Estimated
K D	Half-saturation of CD 4+ T cells	4×10−4g/cm3	[15]
K_T_1	Half-saturation of CD 4+ T cells	2×10−3g/cm3	[15]
K_T_8	Half-saturation of CD 8+ T cells	1×10−3g/cm3	[15]
K_T_r	Half-saturation of Tregs	5×10−4g/cm3	[15]
K E	Half-saturation of endothelial cells	2.5×10−3g/cm3	[43]
K C	Half-saturation of tumor cells	0.4 g/cm3	[15]
K_I_12	Half-saturation of IL-12	8×10−10g/cm3	[15]
K_I_2	Half-saturation of IL-2	2.37×10−11g/cm3	[15]
K_T_β	Half-saturation of TGF-β	2.68×10−13g/cm3	[15]
K W	Half-saturation of oxygen	1.69×10−4g/cm3	[43]
K G	Half-saturation of VGEF	7×10−8g/cm3	[43]
K_TT_r	Inhibition of function of T cells by Tregs	5×10−4g/cm3	[15]
KTQ′	Inhibition of function of T cells by PD-1-PD-L1	1.68×10−18g2/cm6	Estimated
K PB	Block of anti-PD-1 perfusion by anti-VEGF (mice)	1×10−8g/cm3	Estimated
K PB	Block of anti-PD-1 perfusion by anti-VEGF (humans)	1×10−9g/cm3	Estimated
D 0	Density of immature DCs	2×10−5 g/ cm3	[59]
T 10	Density of naive CD 4+ T cells	4×10−4 g/ cm3	[15]
T 80	Density of naive CD 8+ T cells	2×10−4 g/ cm3	[15]
E M	Carrying capacity of endothelial cells	5×10−3 g/ cm3	[43]
C M	Carrying capacity of cancer cells	0.8 g/ cm3	[59]
G 0	Threshold VEGF concentration	3.65×10−10 g/ cm3	[43]
T^_1	Density of CD 4+ T cells from lymph node	4×10−3 g/ cm3	[15]
T^_8	Density of CD 8+ T cells from lymph node	2×10−3 g/ cm3	[15]
Ê	Density of endothelial cells from outside of tumor	5×10−3 g/ cm3	Estimated
W ∗	Hypoxia threshold oxygen level	1.69×10−4 g/ cm3	[44]
W 0	Normal threshold oxygen level	4.65×10−4 g/ cm3	[44]

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 γ and γ. The parameter γ 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 γ. We therefore conducted simulations with different choices of γ,γ 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 γ=3×10−8g/cm3·day and γ=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 γ=3.5×10−8g/cm3·day and γ=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 γ=0) anti-VEGF with γ=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

Growth of tumor volume under treatment with γ or γ, or combination (γ,γ). The anti-VEGF or/and the anti-PD-1 treatment started at day 0 and continued for 10 days. a γ=3×10−8g/cm3·day,γ=3×10−8g/cm3·day; b γ=3.5×10−8g/cm3·day,γ=0.5×10−8g/cm3·day. All other parameters are same as in Fig. 3

Fig. 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 γ=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 λ,λ,λ. 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 (γ,γ). 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

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

Schedule	(γB1,γA1)	(γB1,γA2)	(γB2,γA1)	(γB2,γA2)
S1	48.59	39.17	40.4	33.99
S2	16.05	13.26	15.25	12.78
S3	19.22	16.46	18.24	15.83

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 γ=1.2×10−10g/cm3·day,γ=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

Spatial profiles of the densities of dendritic cells, T cells, endothelial cells and cancer cells with γ=1.2×10−10g/cm3·day,γ=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 γ and γ. The horizontal axis scales the dose amount of anti-VEGF, γ, in units of g/cm3·day, and the vertical axis scales the dose amount of anti-PD-1, γ, 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 γ or γ 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

Efficacy maps: The time in weeks (T) at which the tumor volume decreases by 95% from its initial size under treatment with (γ,γ). 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/K) by . 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 γ and γ vary in the same range as in Fig. 8.

Fig. 9

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

Fig. 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 and η8, are negatively correlated with the tumor volume, while parameters that promote the inhibition of immune cells, such as K,K and K, and the proliferation of endothelial cells and Tregs, such as and λ, are positively correlated with the tumor volume. We also see that the inhibition of perfusion of the anti-PD-1 drug by anti-VEGF promotes cancer growth, namely, K 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 (the parameter which increases the inhibition of CD 8+ T cells by PD-1-PD-L1) has the largest impact on increasing tumor volume.

Discussion

Combination therapy for cancer has already shown improved benefits over a single agent therapy [50]. But significant challenges remain, as seen, for example, from the fact that of the successful phase II clinical trials with combination chemotherapy, less then 4% demonstrated, in Phase III, improvement of care within 5 years [11]. It was suggested in [13] that a selection of clinical trials needs to include scientific rationale underlying the data and hypothesis for the combination. Combination clinical trials should be preceded by analysis of the expected interactions between the diverse agents, addressing, for example the following questions:

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?

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?

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 [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 anti-VEGF, and they could play an important role in the design of clinical trials, where scheduling strategies may significantly affect the outcome.

Conclusions

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

Appendix

Parameter estimation

Half-saturation

In an expression of the form where Y is activated by X, the parameter K is called the half-saturation of X.

We assume that if the average density (or concentration) of X, converges to a steady state X0, 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 and δ are respectively the molecular weight and diffusion coefficient of VEGF, M is the molecular weight of p, and M=24kDa [56] and δ=8.64×10−2cm2 day−1 [57]. Since, M=149kDa (bevacizumab), we get δ=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 δ=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 , where by [15], d=0.1/day, C=K=0.4g/cm3,D=K=4×10−4g/cm3,D0=2×10−5g/cm3. We assume that K=4K where K=7×10−8g/cm3 [43]; hence λ=2.5dD/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 λ, taking λ=17.5/day in mice and λ=7.5/day in humans.

Eqs. (3) and (4).

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

From the steady state of Eq. (4), we have where by [15], /day, /day, . Hence .

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 and . The same considerations apply in the case of T.

Eq. (5).

We assume that TGF-β activates Tregs more than VEGF does, and take . From the steady state of Eq. (5), we get, , where [59], and /day [59]. Hence and /day.

Eq. (6).

By [43], d=0.69/day, E=5×10−3g/cm3,K=2.5×10−3g/cm3,G0=3.65×10−10. We take E=2K=5×10−3g/cm3. From the steady state of Eq. (6), we get λ=2d/(K−G0)=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 λ, taking λ=2.77×107cm3/g·day in mice and λ=2.08×107cm3/g·day in humans.

Eq. (7).

We take d=0.17day−1,C=0.8g/cm3 [59] and λ=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=K=1.69×10−4g/cm3 (see the estimates for Eq. (11)). From the steady state of Eq. (7) in the control case we have, where ; we take η8=2η1; hence and η1=30.19cm3/g·day. Since in the control case the tumor grows, we increase the growth rate of cancer cells taking λ=2.24/day in mice and λ=1.76/day in humans.

Eq. (10).

From the steady state of Eq. (10) we have, , where , and, by [15], and , and by [61], /day. Hence /day.

Eq. (11).

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

Eq. (12).

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

Eqs. (13)-(17).

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

and

In steady state with and G=K we have, by Eq. (17), K=σKK. We take . Hence, and where .

Eqs. (18)-(19).

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

Since the molecular mass of anti-PD-1 (32kDa [62]) approximates the molecular mass of PD-1 (20.5 - 40 kDa [62]), we take μ=μ=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 . We assume that 90% of B is depleted in blocking of VEGF, while the remaining 10% degrades naturally. Hence, The molecular mass of anti-VEGF is approximately 6 times larger than that of VEGF, so we take μ=6μ=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:

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

where , and . The mesh moves by , where 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.