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][10][11][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][16][17][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][30][31][32][33]. There are also several clinical trials with OV and checkpoint inhibitors [34][35][36][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 γ V the dose amount of the injected OV and by γ A the dose amount of the checkpoint inhibitor, and define the efficacy of the treatment by (γ V , γ A ) 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 (γ V , γ A ) plane where an increase in γ A results in actual decrease in the efficacy. We denote by l V i the replication rate of viruses within infected cancer cells. We then construct efficacy maps for (l V i , γ A ) and find regions where an increase in γ A 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.
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 (l V i , λ DV ) and The mathematical model is represented by a system of partial differential equations based on

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/cm 3 . The time unit is 1 day.
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 (C).
We assume a logistic growth for cancer cells, and that cancer cells are killed primarily by CD8 + T cells at a rate η 8 T 8 C where η 8 is a constant. Cancer cells become infected by V e at a rate proportional to CV e . Therefore C satisfies the following equation: where δ C is the dispersion coefficient and d C is the death rate by apoptosis.

Equation for infected cancer cells (C i ).
We assume that CD8 + T cells kill infected cancer cells at a rate Z 8C i T 8 C i , where Z 8C i is a constant larger than η 8 . We also assume that macrophages kill infected cancer cells by phagocytosis [40] at a rate proportional to C i M. The death rate of infected cancer cells is larger than the death rate of uninfected cancer cells by a factor m V i V i which represents the effect of viral burden. We take the rate by which cancer cells become infected to be β C CV e . We finally assume that the density of the V i cells is proportional to the density of the C i cells within which the V i reside, so that they have the same dispersion coefficient. Hence the equation for infected cancer cells is given by Equation for extracellular virus (V e ). We assume that virus at amount γ V is injected into the tumor at successive days t 1 , t 2 , …, t n . Thus at each day t j we have to increase V e by an amount γ V , that is, V e (t j + 0) − V e (t j − 0) = γ V . 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 T 8 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 e . Hence, the equation for V e 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 β V is related to the coefficient β C in Eq (2) by the equation where m VC is the ratio of the mass of one virus to one cancer cell.

Equation for intercellular virus (V i ).
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 l V i C i . The equation for V i is the following: The last two terms represent a loss of V i due to death of their host C i by the macrophages and CD8 + T cells. Note that V i moves with the same velocityũ as C i .

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

Equation for dentritic cells (D).
Oncolytic virus is often armed to elicit adaptive immune response [19,24]. In particular, we assume that inactive dendritic cells with density D 0 are activated by intracellular armed viruses at a rate proportional to D 0 V i . 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 D 0 C K C þC , 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 + T cells (T 1 ). Naive CD4 + T cells are activated by IL-12 while in direct contact with dendritic cells. IL-2 induces proliferation of activated T 1 cells [44,45]. Both processes are inhibited by the complex PD-1-PD-L1 (Q) [46], by a factor 1 1þQ=K TQ . Hence T 1 satisfies the following equation:

Equation for CD8 + T cells (T 8 ).
IL-12 activates CD8 + T cells and IL-2 induces proliferation of CD8 + T cells [44,45]. Hence, similarly to the equation for T 1 , T 8 satisfies the following equation: Equation for IL-12 (I 12 ). IL-12 is produced by activated DCs, so that The diffusion coefficient of I 12 is several orders of magnitude larger than the diffusion coefficient of cells. Hence the transport term r Á (uI 12 ) is negligible compared to the diffusion term d I 12 r 2 I 12 , and it was therefore omitted.
Equation for IL-2 (I 2 ). IL-2 is produced by activated CD4 + T cells. Hence, Here again the transport term was omitted.

Equation for PD-1 (P), PD-L1 (L) and PD-1-PD-L1 (Q).
PD-1 is expressed on the surface of activated CD4 + T cells and activated CD8 + T cells. Hence, P is given by P = ρ P (T 1 + T 8 ), where ρ P 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 r P ¼ P T 1 þT 8 . 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 PA, so that In the sequal we take the dimension of μ PA to be cm 3 /g Á day so that A is given in unit of g/cm 3 .
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 (T 1 + T 8 ) and C: where ρ L 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 α PL and d Q , 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 PL = d Q Q, or where σ = α PL /d Q .

Equation for anti-PD-1 (A).
We assume that anti-PD-1 is injected intraperitoneally in the amount γ A at the same days t 1 , t 2 , …, t n as in the injection of virus. The PK/PD effect of the drug is assumed to be The drug A is depleted in the process of blocking PD-1. Hence,

Equation for cells velocity (u):
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 y Â r Á u ¼ Right À hand side of Eqs ð2Þ; ð3Þ and ð6Þ-ð9Þ: 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 r Á u.
We assume that the average density of each cell type eventually stabilizes with the following values: For cancer cells, 0.4 g/cm 3 ; for dendritic cells 0.4 × 10 −4 g/cm 3 ; for macrophage, 0.2 g/ cm 3 ; for T 1 cells, 2 × 10 −3 g/cm 3 ; and for T 8 cells, 1 × 10 −3 g/cm 3 . 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 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 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 densitiesT 8 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 s T I 12 ð Þ ¼ a T It is implicity assumed that receptors P become expressed only after T 1 and T 8 cells were already inside the tumor.

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). 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 γ V and γ A which represent the effective doses for this subject. We also note that γ V and γ A 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 γ V and γ A , 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]. 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.
We can characterize the anticancer effectiveness of a virus by (i) its ability to replicate within cancer cells, as represented by the parameters l V i in Eq (5), and (ii) by its ability to stimulate the anticancer immune response, as represented by the activation rate λ DV in Eq (7). For any pair (l V i , λ DV ) 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 γ V , and the anti-PD-1 is given at the beginning of weeks 1,4,7,10, 13, 16 at an amount γ A . We denote by V 24 (γ V , γ A ) the  Tables  2 and 3. Initial values are as in (21). 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 V 24 (γ V , γ A ) is decreased. Fig 4 shows an efficacy map for the parameters l V i ¼ 5 Â 10 À 4 =day, λ DV = 5.2 × 10 10 cm 3 / g Á day. For clarity we marked tumor volumes V 24 (γ V , γ A ) on the equi-efficacy curves. We see that as γ V increases so does the efficacy. However, the same is not true of γ A : there are regions where the efficacy decreases as γ A increases. To understand what happens in such regions we take two points with the same γ V : (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 γ A is somewhat larger than the tumor volume for the smaller γ A . 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.
We find the same phenomenon in Fig 7, which is an efficacy map for (γ A , l V i ), for specific values of γ V = 2.5 × 10 −7 g/cm 3 /day and λ DV = 5.2 × 10 10 cm 3 /g Á day. The tumor volume decreases as l V i increases, but there are values of l V i for which the tumor volume increases when γ A is increased.
We note that λ DV and γ A 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 γ A 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

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][30][31][32][33]. There are recent clinical trials with OV and checkpoint inhibitors [34][35][36][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 (l V i ) and its immunogenicity potential (λ DV ). We first simulated a treatment corresponding to mice experiments, where the dose γ A of anti-PD-1 and the dose γ V 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 (γ V , γ A ) was given for 14 weeks, and we observed the results of the treatment for 24 weeks. We simulated the tumor volume V 24 (γ V , γ A ) at the end of 24 weeks for a range of (γ V , γ A ). We found that there are regions in the (γ V , γ A )-plane where an increase in γ A 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 V 24 (l V i , γ A ) for a fixed parameter γ V and variable (l V i , γ A ). We again found regions of antagonism where an increase in γ A 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. 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;  [39] l T 1 I 2 activation rate of CD4 + T cells by IL-2 0.25 day −1 [39] l T 8 I 12 activation rate of CD8 + T cells by IL-12 8.30 day −1 [39] l T 8 I 2 activation rate of CD8 + T cells by IL-2 0.25 day −1 [39] λ I 12 D production rate of IL-12 by DCs 2.76 × 10 −6 day −1 [39] l I 2 T 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 × 10 4 cm 3 /g Á day estimated  [50] Ã In this reference the value was estimated but not obtained directly from experimental results.
K C half-saturation of tumor cells 0.4 g/cm 3 [50] K D half-saturation of DCs 0.4 × 10 −4 g/cm 3 [39] K I 12 half-saturation of IL-12 1.5 × 10 −10 g/cm 3 [50] K I 2 half-saturation of IL-2 2.37 × 10 −11 g/cm 3 [50] K T 1 half-saturation of CD4 + T cells 2 × 10 −3 g/cm 3 [39] K T 8 half-saturation of CD8 + T cells 1 × 10 −3 g/cm 3 [39] K 0 TQ inhibition of function of T cells by PD-1-PD-L1 but the dispersion coefficient of macrophages was increased by up to a factor of 4 because they are highly mobile. 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 δ p and molecular weight M p . 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 Y X K X þX where Y is activated by X, the parameter K X is called the half-saturation of X. If X reaches a steady state X 0 , we expect that X 0 /(K X + X 0 ) will not be "too close" to 0 and not "too close" to 1. For definiteness we take X 0 /(K X + X 0 ) to be 1/2, so that the steady state X 0 is derived from experimental or clinical data. Combination therapy with oncolytic virus and checkpoint inhibitor Eq (2). We take λ C = 0.65/day, which is slightly smaller than in [38], and β C = 9 × 10 4 cm 3 / g Á day, which is slightly larger than in [52], and we take η 8 = 1.38 × 10 2 cm 3 /g Á day, which is slightly larger than in [39].
Eq (3). We assume that T 8 cells kill C i much more efficiently than they kill C, and take Z 8C i ¼ 55Z 8 ¼ 7:59 Â 10 3 cm 3 =g Á day. 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 m V i ¼ 2 Â 10 7 cm 3 =g so that m V i V i 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/cm 3 . Macrophages engulf infected cancer cells [40], but we assume that the rate is extremely small compared to the rate by which T 8 kill infected cancer cells. We accordingly take μ C i M = 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 VC = 10 −6 . Hence β V = β C m VC = 0.09 cm 3 /g Á day. We assume that the clearance rate of V e by macrophages is much larger than the rate by which V e invades uninfected cancer cells, and take μ V e M = 2 cm 3 /g Á day. We assume that the ratio of C i /V i in the first 12 days averages 3 × 10 6 , and that the replication of an intracellular virus occurs approximately every 22-23 hours, so that the growth rate per day is 1.5 × 2 9 V i .  3 . Hence λ M = 0.003/day. We note however that in estimating λ M , we ignored the contribution of r Á (uM), whose integral over the tumor {r < R(t)} is R r¼RðtÞ dRðtÞ dt Á M, which is a positive quantity. Hence, @M @t is actually decreased when we equate to zero the right-hand side (RHS) of Eq (6); we therefore need to increase λ M ; we take λ M = 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, l MC i MC i , is assumed to be 'relatively' small; we take l MC i ¼ 0:04 cm 3 =g. Eq (7). We take λ DC = 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 λ DV V i is comparable to λ DC when V i % 10 −10 g/cm 3 . Accordingly, we take λ DV = 5.2 × 10 10 cm 3 /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: TQ , T 10 , T 80 ,T 1 ,T 8 . 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 δ M in the sensitivity analysis, keeping all other diffusion coefficient equal, and randomly increasing δ M by up to a factor of 4. The results are shown in Fig 8. We see that increasing the source of T cells (T 10 ; T 80 ;T 10 ;T 80 ) decreases the tumor volume, as does the depletion rate (μ PA ) 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 l V i , l MC i , λ DV , l T 1 I 12 and l T 8 I 12 , and the parameters β C , m V i , μ C i M , η 8 , Z 8C i and d V e M 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 μ C i M , μ V e M , Z 8C i and l MC i are positively correlated with the tumor volume, while the parameters that promote viral infection, such as β C , m V i and l V i , are negatively correlated with the tumor volume. We also see that the production/activation rates that promote effector T cells, namely, λ DV , l T 1 I 12 and l T 8 I 12 , 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 r k i and C k i denote numerical approximations of i-th grid point and Cðr k i ; ntÞ, respectively, where τ is the size of time-step. The discretization of Eq (22) is derived by the fully implicit finite difference scheme: where C r ¼