Combination therapy for cancer with oncolytic virus and checkpoint inhibitor: A mathematical model.

Oncolytic virus (OV) is a replication competent virus that selectively invades cancer cells; as these cells die under the viral burden, the released virus particles proceed to infect other cancer cells. Oncolytic viruses are designed to also be able to stimulate the anticancer immune response. Thus, one may represent an OV by two parameters: its replication potential and its immunogenicity. In this paper we consider a combination therapy with OV and a checkpoint inhibitor, anti-PD-1. We evaluate the efficacy of the combination therapy in terms of the tumor volume at some later time, for example, 6 months from initial treatment. Since T cells kill not only virus-free cancer cells but also virus-infected cancer cells, the following question arises: Does increasing the amount of the checkpoint inhibitor always improve the efficacy? We address this question, by a mathematical model consisting of a system of partial differential equations. We use the model to construct, by simulations, an efficacy map in terms of the doses of the checkpoint inhibitor and the OV injection. We show that there are regions in the map where an increase in the checkpoint inhibitor actually decreases the efficacy of the treatment. We also construct efficacy maps with checkpoint inhibitor vs. the replication potential of the virus that show the same antagonism, namely, an increase in the checkpoint inhibitor may actually decrease the efficacy. These results have implications for clinical trials.

Introduction

PD-1 is an immunoinhibitory receptor predominantly expressed on activated T cells [1, 2]. Its ligand PD-L1 is upregulated on the same activated T cells, and in some human cancer cells [2, 3]. The complex PD-1-PD-L1 is known to inhibit T cell function [1]. Immune checkpoints are regulatory pathways in the immune system that inhibit its active response against specific targets. In the case of cancer, the complex PD-1-PD-L1 functions as an immune checkpoint for anti-tumor T cells. CTLA-4 is another immunoinhibitory receptor expressed on activated T cells; when it combines with its ligand B7 on dendritic cells, the complex CTLA-4-B7 acts as a checkpoint inhibitor for anti-tumor T cells [4, 5]. There has been much progress in recent years in developing checkpoint inhibitors, primarily anti-PD-1 and anti-PD-L1 [6], and anti-CTLA-4 [7, 8].

Oncolytic virus (OV) is a genetically engineered virus that can selectively invade into and replicate within cancer cells while not harming normal healthy cells. OV therapy has been explored as an approach to combat cancer, and clinical trials were carried out on different types of cancer [9-12]. However, therapeutic efficacy remains a challenge [13, 14]. One of the factors that limits OV therapy is the antigenicity of the infected cells; the macrophages of the innate immune system recognize these cells and destroy them together with the virus particles inside them. For this reason, experimental studies considered combination of OV therapy with immune suppressive drugs [15-18].

In another direction, some studies consider OV with viruses designed to both replicate within cancer cells and stimulate cytotoxic T cells; such viruses include vesicular stomatitis virus [19, 20], Newcastle Disease Virus [21], vaccinia [22, 23], measle virus [24], and others [25, 26]. Advances in the design of various oncolytic viruses are reported in [27, 28]. The underlying assumption in these studies is that the virus will survive long enough, under the pressure of the innate immune attack, to activate a sufficiently large number of cytotoxic T cells that will eradicate or significantly reduce the cancer. To make this approach more effective, it was suggested to combine the OV drug with checkpoint inhibitors. Several mouse experiments, with different types of cancer cells, reported that both CTLA-4 and PD-L1 checkpoints blockade enhanced the OV therapy [29-33]. There are also several clinical trials with OV and checkpoint inhibitors [34-37].

In previous work the authors considered combination therapies with checkpoint inhibitor and, as a second agent, tumor vaccine [38] or BRAF inhibitor [39]. In the present paper the second agent is oncolytic virus. This poses a dilemma, since T cells kill not only virus-free cancer cells but also virus-infected cancer cells (thus reducing the anti-cancer effect of the virus), while checkpoint inhibitors enhance the T cells activities. Thus, it is natural to ask whether increasing the amount of the checkpoint inhibitor does always result in a decrease in tumor volume. We develop a mathematical model to address this question. We denote by γ the dose amount of the injected OV and by γ the dose amount of the checkpoint inhibitor, and define the efficacy of the treatment by (γ, γ) in terms of the tumor volume at some arbitrary time, for example, 24 weeks from the beginning of the treatment. We use the mathematical model to develop an efficacy map, and we find that there are regions in (γ, γ) plane where an increase in γ results in actual decrease in the efficacy. We denote by the replication rate of viruses within infected cancer cells. We then construct efficacy maps for (, γ) and find regions where an increase in γ again results in decreased efficacy. In such regions, the indiscriminent killing of infected and uninfected cancer cells has pro-cancer effect. These have implications for clinical trials.

The mathematical model includes CD4+ Th1 cells and CD8+ T cells, macrophages, and dendritic cells. Dendritic cells are activated indirectly by the virus and by necrotic cancer cells, while macrophages are activated by virus-infected cancer cells. Macrophages engulf and destroy infected cancer cells, but they also kill, at a lesser rate, uninfected cancer cells. When a cancer cell is infected by an extracellular virus, the extracellular virus becomes an intracellular virus within the infected cell. Intracellular viruses multiply within the cancer cells and cause them to lyse, thereby releasing all their viruses to the extracellular environment. T cells are activated by IL-12 produced by dendritic cells, and they also proliferate by IL-2 produced by Th1 cells. Fig 1 shows the network of interactions among the cells, with PD-1 and PD-L1 on T cells and PD-L1 also on tumor cells.

Fig 1

Interaction of tumor cells with virus and immune cells.

Sharp arrows indicate proliferation/activation, blocked arrows indicate killing/blocking, and dashed lines indicate proteins on T cells. C: uninfected cancer cells, C: infected cancer cells, V: extracellular virus, V: intracellular virus, D: dendritic cells, T1: CD4+ Th1 cells, T8: CD8+ T cells, I2: IL-2, I12: IL-12, P: PD-1, L: PD-L1, Q: PD-1-PD-L1 complex.

We assume that the treatment with combination therapy extends over a period of 16 weeks, and we evaluate the results of the treatment at the end of 24 weeks. We can use the model to compute the tumor volume at the end of 24 weeks for each pair of parameters of (, λ) and doses (γ, γ).

The mathematical model is represented by a system of partial differential equations based on Fig 1.

Mathematical model

The mathematical model is based on the diagram in Fig 1. The list of variables is given in Table 1, where the density of cells and concentration of cytokines are all in unit of g/cm3. The time unit is 1 day.

Table 1

List of variables (in units of g/cm3).

Notation	Description	Notation	Description
C	density of cancer cells	T8	density of activated CD8+ T cells
Ci	density of infected cancer cells	I12	IL-12 concentration
V	density of extracellular virus	I2	IL-2 concentration
Vi	density of intracellular virus	P	PD-1 concentration
M	density of macrophages	L	PD-L1 concentration
D	density of dentritic cells	Q	PD-1-PD-L1 concentration
T1	density of activated CD4+ T cells	A	anti-PD-L1 concentration

We assume that the total density of cells within the tumor remains constant in space and time, so that

Since cancer cells proliferate while T cells and macrophages enter the tumor, the assumption (1) implies that there is an internal pressure among cells, and this gives rise to a velocity u of cells.

Equation for uninfected cancer cells (. We assume a logistic growth for cancer cells, and that cancer cells are killed primarily by CD8+ T cells at a rate η8 T8 C where η8 is a constant. Cancer cells become infected by V at a rate proportional to CV. Therefore C satisfies the following equation: where δ is the dispersion coefficient and d is the death rate by apoptosis.

Equation for infected cancer cells (. We assume that CD8+ T cells kill infected cancer cells at a rate , where is a constant larger than η8. We also assume that macrophages kill infected cancer cells by phagocytosis [40] at a rate proportional to C M. The death rate of infected cancer cells is larger than the death rate of uninfected cancer cells by a factor V which represents the effect of viral burden. We take the rate by which cancer cells become infected to be β CV. We finally assume that the density of the V cells is proportional to the density of the C cells within which the V reside, so that they have the same dispersion coefficient. Hence the equation for infected cancer cells is given by

Equation for extracellular virus (. We assume that virus at amount γ is injected into the tumor at successive days t1, t2, …, t. Thus at each day t we have to increase V by an amount γ, that is, V(t + 0) − V(t − 0) = γ. This increase can be written in the form where δ(s) is the Dirac measure. We assume that when an infected cell dies the intracellular viral particles are released into the tumor microenvironment; however, when an infected cell is killed by macrophages or T8 cells, the virus particles inside it are cleared out. Extracellular virus are endocytosed by macrophages, and the rate of their depletion is proportional to MV. Hence, the equation for V takes the following form: where N is the average number of viral particles released at death of an infected cancer cell. Note that the coefficient β is related to the coefficient β in Eq (2) by the equation β = β m, where m is the ratio of the mass of one virus to one cancer cell.

Equation for intercellular virus (. Viruses multiply in a cancer cell by exploiting the DNA of the cell as a ‘resource’. We represent the proliferation of the viruses in the cell by C. The equation for V is the following:

The last two terms represent a loss of V due to death of their host C by the macrophages and CD8+ T cells. Note that V moves with the same velocity as C.

Equation for macrophages (. The growth rate of the proinflammatory macrophages is promoted by infected cancer cells and is represented by a term . Hence M satisfies the following equation: where λ is a source of macrophages prior to the treatment with OV.

Equation for dentritic cells (. Oncolytic virus is often armed to elicit adaptive immune response [19, 24]. In particular, we assume that inactive dendritic cells with density D0 are activated by intracellular armed viruses at a rate proportional to D0V. Dendritic cells are also activated by HMGB-1 [41, 42], which is produced by necrotic cancer cells (NCs) [43]. We assume that the concentration of HMGB-1 is proportional to the density of NCs and that the density of NCs is proportional to the density of cancer cells. Hence, the activation rate of inactive dendritic cells is proportional to , where the Michaelis-Menten law is used to account for the limited receptor recycling time which occurs in the process of DC activation. The dynamics of DCs is given by

Equation for CD4. Naive CD4+ T cells are activated by IL-12 while in direct contact with dendritic cells. IL-2 induces proliferation of activated T1 cells [44, 45]. Both processes are inhibited by the complex PD-1-PD-L1 (Q) [46], by a factor . Hence T1 satisfies the following equation:

Equation for CD8. IL-12 activates CD8+ T cells and IL-2 induces proliferation of CD8+ T cells [44, 45]. Hence, similarly to the equation for T1, T8 satisfies the following equation:

Equation for IL-12 (. IL-12 is produced by activated DCs, so that

The diffusion coefficient of I12 is several orders of magnitude larger than the diffusion coefficient of cells. Hence the transport term ∇ ⋅ (uI12) is negligible compared to the diffusion term , and it was therefore omitted.

Equation for IL-2 (. IL-2 is produced by activated CD4+ T cells. Hence,

Here again the transport term was omitted.

Equation for PD-1 (. PD-1 is expressed on the surface of activated CD4+ T cells and activated CD8+ T cells. Hence, P is given by P = ρ(T1 + T8), where ρ is the ratio of the mass of all the PD-1 proteins in one T cell to the mass of one T cell. Thus, P satisfies the equation or, by Eqs (8) and (9), where . Note that P undergoes the same advection velocity as the T cells. We assume that PD-1 is depleted (or blocked) by A at rate μPA, so that

In the sequal we take the dimension of μ to be cm3/g ⋅ day so that A is given in unit of g/cm3.

PD-L1 is expressed on the surface of activated CD4+ T cells, activated CD8+ T cells, and on tumor cells. Hence, the concentration of PD-L1 (L) is proportional to (T1 + T8) and C: where ρ is the ratio of the mass of all the PD-L1 proteins in one T cell to the mass of one T cell, and ε depends on the specific type of tumor.

PD-1 and PD-L1 form a complex PD-1-PD-L1 (Q), with association and disassociation rates α and d, respectively:

The half-life of Q is less then 1 second (i.e. 1.16 × 10−5 day) [47]. Hence, we may assume that the dynamics in (14) is in quasi-steady state, so that α PL = d Q, or where σ = α/d.

Equation for anti-PD-1 (. We assume that anti-PD-1 is injected intraperitoneally in the amount γ at the same days t1, t2, …, t as in the injection of virus. The PK/PD effect of the drug is assumed to be , where H(s) = 0 if s ≤ 0, H(t) = 1 if s > 1. The drug A is depleted in the process of blocking PD-1. Hence,

Equation for cells velocity (: Cells disperse within the tissue, and its random motility may vary from one cell type to another. If the differences in the dispersion coefficients are ignored, then by adding the equations for all the cells and using Eq (1), we get

To simplify the model we assume that the differences between the dispersion coefficients of the different cell types are small (but see comments in “Parameter estimation” (in “Diffusion coefficients”) and “Sensitivity analysis”), and proceed to use the above equation for ∇ ⋅ u.

We assume that the average density of each cell type eventually stabilizes with the following values: For cancer cells, 0.4 g/cm3; for dendritic cells 0.4 × 10−4 g/cm3; for macrophage, 0.2 g/cm3; for T1 cells, 2 × 10−3 g/cm3; and for T8 cells, 1 × 10−3 g/cm3. Recalling Eq (1) we find that θ = 0.6034, so that

To simplify the computations, we assume that the tumor is spherical with moving boundary r = R(t), and that all the densities and concentrations are radially symmetric, that is, 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 (: We assume that the free boundary r = R(t) moves with the velocity of cells, so that

Boundary conditions We assume that naive CD4+ T cells and CD8+ T cells which migrated from the lymph nodes into the tumor microenvironment have constant densities at the tumor boundary, and that they are activated by dentritic cells and IL-12 upon entering the tumor. We represent this process by the flux conditions at the boundary: where .

We impose zero-flux boundary condition on all the remaining variables:

It is implicity assumed that receptors P become expressed only after T1 and T8 cells were already inside the tumor.

Initial conditions We prescribe the following values (in unit g/cm3) at day t = 0:

Note that the initial values satisfy Eq (1). We took the initial values for cells to be different from their above assumed asymptotic values. The choice of the initial conditions have little effect on 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).

Fig 2 shows the profiles of the average densities/concentrations of all the variables of the model in the first 30 days in the control case, that is, without treatment. The simulation results show that the steady states of all the cytokines and cells are approximately equal to the half-saturation values that we assumed in estimating the parameters of the model.

Fig 2

Average densities/concentrations, in g/cm3, of all the variables in the model in the control case.

All parameter values are the same as in Tables 2 and 3. Initial values are as in (21).

Table 2

Summary of parameter values.

Notation	Description	Value used	References
δD	diffusion coefficient of DCs	8.64 × 10−7 cm2 day−1	[38]*
δ_T_1	diffusion coefficient of CD4+ T cells	8.64 × 10−7 cm2 day−1	[38]*
δ_T_8	diffusion coefficient of CD8+ T cells	8.64 × 10−7 cm2 day−1	[38]*
δC	diffusion coefficient of tumor cells	8.64 × 10−7 cm2 day−1	[38]*
δM	diffusion coefficient of macrophages	8.64 × 10−7 cm2 day−1	[38]*
δ_I_12	diffusion coefficient of IL-12	6.0472 × 10−2 cm2 day−1	[39]
δ_I_2	diffusion coefficient of IL-2	9.9956 × 10−2 cm2 day−1	[39]
δA	diffusion coefficient of IL-2	4.73 × 10−2 cm2 day−1	[39]
αT	flux rate of T cells on the boundary	1 cm−1	estimated
λC	growth rate of cancer cells	0.65 day−1	estimated
λ_V_i	growth rate of intracellular virus	6 × 10−4 day−1	estimated
λM	growth rate of macrophages	0.009 day−1	[51]*
λ_MC_i	activation rate of macrophages by Ci	0.04 cm3/g	estimated
λDV	activation rate of DCs by virus infection	5.2 × 1010 cm3/g ⋅ day	estimated
λDC	activation rate of DCs by tumor cells	5.2 day−1	estimated
λ_T_1I_12	activation rate of CD4+ T cells by IL-12	9.32 day−1	[39]
λ_T_1I_2	activation rate of CD4+ T cells by IL-2	0.25 day−1	[39]
λ_T_8I_12	activation rate of CD8+ T cells by IL-12	8.30 day−1	[39]
λ_T_8I_2	activation rate of CD8+ T cells by IL-2	0.25 day−1	[39]
λI12D	production rate of IL-12 by DCs	2.76 × 10−6 day−1	[39]
λ_I_2T_1	production rate of IL-2 by CD4+ T cells	2.82 × 10−8 day−1	[39]
βC	infection rate of cancer cells by virus	9 × 104 cm3/g ⋅ day	estimated
βV	rate of transition from Ve to Vi by infection	0.09 cm3/g ⋅ day	estimated
μCiM	killing rate of Ci by M	4.8 × 10−2 cm3/g ⋅ day	estimated
μVeM	clearance rate of Ve by M	2 cm3/g ⋅ day	estimated
μ_V_i	death rate of infected cell due to viral burden	5 × 107 day−1	estimated
N	burst size of Vi from natural death of Ci	100	estimated
η8	killing rate of tumor cells by CD8+ T cells	1.38 × 102 day−1 ⋅ cm3/g	estimated
η_8C_i	killing rate of infected cancer cells by CD8+ T cells	7.59 × 103 day−1 ⋅ cm3/g	estimated
μPA	blocking rate of PD-1 by anti-PD-1	6.87 × 104 cm3/g ⋅ day	[39]
ρP	expression of PD-1 in T cells	2.49 × 10−7	[39]
ρL	expression of PD-L1 in T cells	5.22 × 10−7	[39]
ε	relative expression of PD-L1 in tumor cells	0.01	[39]
dC	death rate of uninfected tumor cells	0.17 day−1	[50][*]
dM	death rate of macrophages	0.015 day−1	[51]
dD	death rate of DCs	0.1 day−1	[50]
d_T_1	death rate of CD4+ T cells	0.197 day−1	[50]
d_T_8	death rate of CD8+ T cells	0.18 day−1	[50]
d_I_12	degradation rate of IL-12	1.38 day−1	[50]
d_I_2	degradation rate of IL-2	2.376 day−1	[50]

* In this reference the value was estimated but not obtained directly from experimental results.

Table 3

Summary of parameter values.

KC	half-saturation of tumor cells	0.4 g/cm3	[50]
KD	half-saturation of DCs	0.4 × 10−4 g/cm3	[39]
K_I_12	half-saturation of IL-12	1.5 × 10−10 g/cm3	[50]
K_I_2	half-saturation of IL-2	2.37 × 10−11 g/cm3	[50]
K_T_1	half-saturation of CD4+ T cells	2 × 10−3 g/cm3	[39]
K_T_8	half-saturation of CD8+ T cells	1 × 10−3 g/cm3	[39]
KTQ′	inhibition of function of T cells by PD-1-PD-L1	1.365 × 10−18 g/cm3	[39]*
θ	total cell density	0.6034 g/cm3	**
D0	density of immature DCs	2 × 10−5 g/cm3	[50]
T10	density of naive CD4+ T cells	4 × 10−4 g/cm3	[39]*
T80	density of naive CD8+ T cells	2 × 10−4 g/cm3	[39]*
CM	carrying capacity of cancer cells	0.8 g/cm3	[50]
T^_1	density of CD4+ T cells from lymph node	4 × 10−3 g/cm3	[39]*
T^_8	density of CD8+ T cells from lymph node	2 × 10−3 g/cm3	[39]*

* In this reference the value was estimated but not obtained directly from experimental results.

** The value is determined by Eq (1) with steady state densities of the cells.

We proceed to simulate the treatment of cancer by OV and anti-PD-1 as single agents, and by a combination of the two drugs. Following mice experiments reported in [49], we apply the OV injections in days 0,2,4, and anti-PD-1 injection in days 4,7,11. From Figs 1(b) and 2(b) in [49] we see that although all the mice were identical and were treated with the same amounts of dose, their responses were varied: For some mice the tumor volume grew faster with OV than with anti-PD-1 as single agents, while for others this was the reverse, which means that the “effective” dose amounts varied with each subject. We account for this, in our model, by taking for each mouse somewhat different values of γ and γ which represent the effective doses for this subject. We also note that γ and γ should be approximately proportional to the amount of dose injected in the experiments. We determined the proportionality coefficients, or rather the orders of the magnitude of γ and γ, so that the doubling time of the tumor volume (under treatment with a single agent) is approximately 20 days, which is the case for a large number of the mice in Figs 1(b), 2(b) of [49].

Fig 3(a)–3(c) show some simulations of the model with different values of γ and γ. The profiles are similar to many of those given in [49]. In Fig 3(a) treatment with anti-PD-1 as single agent reduces tumor growth more than treatment with OV as single agent, and in Fig 3(b) and 3(c), it is the reverse, in agreement with profiles in [49]. In all cases, the combination reduces the tumor growth more than a single agent.

Fig 3

The growth of tumor volume.

OV is given at days t = 0, 2, 4 with the amount γ and anti-PDE-1 is given at days t = 4, 7, 11 with the amount γ. (a) γ = 0.1 × 10−10 g/cm3, γ = 8 × 10−7 g/cm3. (b) γ = 0.2 × 10−10 g/cm3, γ = 3 × 10−7 g/cm3. (c) γ = 0.5 × 10−10 g/cm3, γ = 7 × 10−7 g/cm3. Parameter values are the same as in Fig 2.

We can characterize the anticancer effectiveness of a virus by (i) its ability to replicate within cancer cells, as represented by the parameters in Eq (5), and (ii) by its ability to stimulate the anticancer immune response, as represented by the activation rate λ in Eq (7). For any pair (, λ) we may associate a “virtual virus” having these two parameters.

We proceed to use the mathematical model to conduct in silico clinical trials. As an example, a treatment will be given for a period of 16 weeks, and the patient’s tumor volume will be measured at the end of 24 weeks from the initial treatment. The virus is injected into the tumor at the beginning of weeks 1,3,5,7,9,11,13 and 15, at an amount γ, and the anti-PD-1 is given at the beginning of weeks 1,4,7,10, 13, 16 at an amount γ. We denote by V24(γ, γ) the volume of the tumor at the end of 24 weeks, and define the efficacy of the treatment by the formula: thus the efficacy is increased if the tumor volume V24(γ, γ) is decreased.

Fig 4 shows an efficacy map for the parameters , λ = 5.2 × 1010 cm3/g ⋅ day. For clarity we marked tumor volumes V24(γ, γ) on the equi-efficacy curves. We see that as γ increases so does the efficacy. However, the same is not true of γ: there are regions where the efficacy decreases as γ increases. To understand what happens in such regions we take two points with the same γ: (3.7 × 10−7, 5.4 × 10−8) and (3.7 × 10−7, 6.6 × 10−8). Fig 5 shows that the tumor volume for the larger γ is somewhat larger than the tumor volume for the smaller γ. Fig 6 explains what has actually occurred. With the higher dose, more infected cancer cells were killed, and the virus population decreased. Hence the number of activated dendritic cells decreased and then also the number of T cells decreased, which resulted in an increase in the number of uninfected cancer cells.

Fig 4

The tumor volume at week 24 for different pair of (γ, γ).

Here , λ = 5.2 × 1010 cm3/g ⋅ day, γ = 1 × 10−7−4 × 10−7 g/cm3 and γ = 0.6 × 10−8−9 × 10−8 g/cm3. All other parameter values are the same as in Tables 2 and 3.

Fig 5

Growth of tumor volume.

Here, γ = 3.7 × 10−7 g/cm3. Other parameter values are the same as in Fig 4.

Fig 6

The average densities and tumor volume.

Blue: γ = 3.7 × 10−7 g/cm3, γ = 5.4 × 10−8 g/cm3. Red: γ = 3.7 × 10−7 g/cm3, γ = 6.6 × 10−8 g/cm3. Other parameter values are the same as in Fig 4.

We find the same phenomenon in Fig 7, which is an efficacy map for (γ, ), for specific values of γ = 2.5 × 10−7 g/cm3/day and λ = 5.2 × 1010 cm3/g ⋅ day. The tumor volume decreases as increases, but there are values of for which the tumor volume increases when γ is increased.

Fig 7

The tumor volume at week 24.

Here, γ = 2.5 × 10−7 g/cm3 and λ = 5.2 × 1010 cm3/g⋅day. All other parameter values are the same as in Tables 2 and 3.

We note that λ and γ are positively correlated, since both are increasing the activity of effector T cells. We therefore expect that, unlike situation in Fig 4, the tumor volume will decrease as γ increases.

One is tempted to replace PDE system by a simpler system of ODEs where the diffusion and advection terms are dropped. However, since the diffusion of cells is several orders of magnitude smaller than diffusion of cytokines and extracellular virus, the ODE system cannot adequately represent the PDE model. For example, Fig 6 show antagonism between the OV and anti-PD-1, whereas this antagonism disappears in the ODE model.

Conclusion

Oncolytic virus (OV) is a genetically modified virus that can selectively invade cancer cells and replicate inside them. When an infected cell dies, its virus particles are released and proceed to infect other cancer cells. OV therapy, as a single agent, had not been successful because macrophages recognize infected cells and kill them together with their viruses. Recent studies use new designs of OV that can stimulate cytotoxic T cells to kill cancer cells before the viral population is significantly depleted by the macrophages. Some of these studies introduce enhancement of the T cells by blocking their checkpoints. Mice experiments demonstrated that both CTLA-4 and PD-L1 checkpoints blockade enhance the OV treatment [29-33]. There are recent clinical trials with OV and checkpoint inhibitors [34-37]. In particular, in clinical trials for melanoma, reported in [36], patients were treated with OV (T-VEC) and anti-CTLA-4 (ipilimumab) for a period of 13 weeks and were observed for an average period of 20 months.

Since T cells kill not only virus-free cancer cells but also virus-infected cancer cells, they may disrupt the anti-cancer effect of the OV. Hence an increase in the dose of the checkpoint inhibitor may actually have a pro-cancer effect. In order to clarify this situation we developed a mathematical model that includes the immune cells (macrophages, dendritic cells, and effective T cells), and characterized the OV by two parameters: the replication potential () and its immunogenicity potential (λ). We first simulated a treatment corresponding to mice experiments, where the dose γ of anti-PD-1 and the dose γ of OV were administered for 11 days, and the tumor volume was observed for 30 days. We found quantitative agreement with experimental results [49].

We then proceeded to use the model to run in silico clinical trials, where the treatment with a combination (γ, γ) was given for 14 weeks, and we observed the results of the treatment for 24 weeks. We simulated the tumor volume V24(γ, γ) at the end of 24 weeks for a range of (γ, γ). We found that there are regions in the (γ, γ)-plane where an increase in γ results in an increase in the tumor volume. Thus, there are regions of antagonism between the two drugs, where an increase in the anti-PD-1 decreases the efficacy of the treatment. We also simulated V24 (, γ) for a fixed parameter γ and variable (, γ). We again found regions of antagonism where an increase in γ results in an increase in the tumor volume.

These results have implications for clinical trials. Indeed for a clinical trial to be successful, the regions of antagonism between the doses of the checkpoint inhibitor and the OV doses should be determined early on, and avoided; these regions may depend on the specific oncolytic virus which is used in the clinical trial.

Materials and methods

Parameter estimation

Many of parameters in Tables 2 and 3 were taken directly from [38, 39, 50, 51]. When a value taken from these references was not obtained directly from experimental papers or was not carefully estimated from experimental results, we added an asterisk “*” next to the reference. These parameters were used in sensitivity analysis (see Figs 8 and 9) in order to see how the tumor value is affected by a random increase or decrease of these parameters by a factor of 2; but the dispersion coefficient of macrophages was increased by up to a factor of 4 because they are highly mobile.

Fig 8

Statistically significant PRCC values (p-value< 0.01) for R(t) at day 60.

Fig 9

Statistically significant PRCC values (p-value< 0.01) for R(t) at day 60.

Diffusion coefficients

The diffusion coefficients of cytokines were computed in [39] based on the formula where A is a constant and p is any protein with diffusion coefficient δ and molecular weight M. The diffusion coefficients of cells may vary depending on the cell type. For simplicity we take them equal, and choose the common value as in [39], but we show in the section on sensitivity analysis that the tumor growth is affected very little by taking different diffusion coefficients for different cell types.

Half-saturation

In an expression of the form where Y is activated by X, the parameter K is called the half-saturation of X. If X reaches a steady state X0, we expect that X0/(K + X0) will not be “too close” to 0 and not “too close” to 1. For definiteness we take X0/(K + X0) to be 1/2, so that the steady state X0 is derived from experimental or clinical data.

Eq (2)

We take λ = 0.65/day, which is slightly smaller than in [38], and β = 9 × 104 cm3/g ⋅ day, which is slightly larger than in [52], and we take η8 = 1.38 × 102 cm3/g ⋅ day, which is slightly larger than in [39].

Eq (3)

We assume that T8 cells kill C much more efficiently than they kill C, and take . We assume that OV is designed to stimulate the immune system while it is in infected cells. Therefore the viral burden does not increase the death rate of infected cells very much. We therefore take so that V is very small compared to 1, e.g. it increases the death of the infected cancer cell by 2% when the viral load is 10−9 g/cm3. Macrophages engulf infected cancer cells [40], but we assume that the rate is extremely small compared to the rate by which T8 kill infected cancer cells. We accordingly take μ = 4.8 × 10−2/day.

Eqs (4) and (5)

We take N = 100, which is in the range considered in [51]. We assume that the ratio of mass of one virus to one cell is m = 10−6. Hence β = βm = 0.09 cm3/g ⋅ day. We assume that the clearance rate of V by macrophages is much larger than the rate by which V invades uninfected cancer cells, and take μ = 2 cm3/g ⋅ day. We assume that the ratio of C/V in the first 12 days averages 3 × 106, and that the replication of an intracellular virus occurs approximately every 22-23 hours, so that the growth rate per day is 1.5 × 29V. Hence and the growth rate of V is then determined by the equation , so that .

Eq (6)

Without OV, we assume λ = d M in steady state, where d = 0.015/day and M = K = 0.2 g/cm3. Hence λ = 0.003/day. We note however that in estimating λ, we ignored the contribution of ∇ ⋅ (uM), whose integral over the tumor {r < R(t)} is , which is a positive quantity. Hence, is actually decreased when we equate to zero the right-hand side (RHS) of Eq (6); we therefore need to increase λ; we take λ = 0.09/day. Since initially tumor is with radius R(0) = 0.01 cm, macrophages had already arrived into the tumor tissue so that the additional increase in macrophages, , is assumed to be ‘relatively’ small; we take .

Eq (7)

We take λ = 5.2/day which is slightly larger than in [39]. We assume that the virus, although having decreased over time, is still effective in activating dendritic cells, so that λ V is comparable to λ when V ≈ 10−10 g/cm3. Accordingly, we take λ = 5.2 × 1010 cm3/g ⋅ day.

Sensitivity analysis

We performed sensitivity analysis with respect to the tumor volume at day 30 for two sets of parameters. The first set consists of parameters marked by “*” in Tables 2 and 3. These parameters were not derived, or not carefully estimated, from experimental or clinical data. These parameters are: δ, δ, δ, , , λ, μ, ε, d, , T10, T80, , . Following the method of [53], 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 (except the diffusion coefficients) from 1/2 to twice its value in Tables 2 or 3. In the simulations of the model we assumed that the diffusion coefficients of all the cell types are equal. What may cause a significant difference in the simulations is actually the differences between the diffusion coefficients of cell types, rather than their actual values. Since macrophages are highly mobile, we chose to include only δ in the sensitivity analysis, keeping all other diffusion coefficient equal, and randomly increasing δ by up to a factor of 4. The results are shown in Fig 8. We see that increasing the source of T cells () decreases the tumor volume, as does the depletion rate (μ) of PD-1 by the PD-1 inhibitor. On the other hand the production rate of PD-L1 by the cancer (ε) increases the tumor volume. An increase of the random mobility of macrophages, by a factor up to 4, only slightly increases the tumor volume.

The second set of parameters in the sensitivity analysis are some production parameters, namely , , λ, and , and the parameters β, , μ, η8, and which play important roles in the dynamics of C. Here again we sampled all the parameters by taking the range of each parameter for 1/2 to twice its value in Tables 2 and 3. The results are shown in Fig 9.

It is interesting to see from Fig 9 that the parameters that promote killing of infected cancer cells, such as μ, μ, and are positively correlated with the tumor volume, while the parameters that promote viral infection, such as β, and , are negatively correlated with the tumor volume. We also see that the production/activation rates that promote effector T cells, namely, λ, and , are negatively correlated to the tumor volume, while the killing rate of uninfected cancer cells cells by CD8+ T cells, η8, is negatively correlated with tumor volume.

Computational method

We employ moving mesh method [48] to numerically solve the free boundary problem for the tumor proliferation model. To illustrate this method, we take Eq (2) as example and rewrite it as the following form: where F represents the term in the right hand side of Eq (2). Let and denote numerical approximations of i-th grid point and , respectively, where τ is the size of time-step. The discretization of Eq (22) is derived by the fully implicit finite difference scheme: where , , , and . The mesh moves by , where is solved by the velocity equation.