Combination therapy of cancer with cancer vaccine and immune checkpoint inhibitors: A mathematical model

In this paper we consider a combination therapy of cancer. One drug is a vaccine which activates dendritic cells so that they induce more T cells to infiltrate the tumor. The other drug is a checkpoint inhibitor, which enables the T cells to remain active against the cancer cells. The two drugs are positively correlated in the sense that an increase in the amount of each drug results in a reduction in the tumor volume. We consider the question whether a treatment with combination of the two drugs at certain levels is preferable to a treatment by one of the drugs alone at ‘roughly’ twice the dosage level; if that is the case, then we say that there is a positive ‘synergy’ for this combination of dosages. To address this question, we develop a mathematical model using a system of partial differential equations. The variables include dendritic and cancer cells, CD4+ and CD8+ T cells, IL-12 and IL-2, GM-CSF produced by the vaccine, and a T cell checkpoint inhibitor associated with PD-1. We use the model to explore the efficacy of the two drugs, separately and in combination, and compare the simulations with data from mouse experiments. We next introduce the concept of synergy between the drugs and develop a synergy map which suggests in what proportion to administer the drugs in order to achieve the maximum reduction of tumor volume under the constraint of maximum tolerated dose.


Introduction
When cancer cells undergo necrosis, they release high mobility group box-1 (HMGB-1) which activates dendritic cells [1][2][3]. Activated dendritic cells (DCs) mature as APC cells and play a critical role in the communication between the innate and adaptive immune responses. Once activated, dendritic cells produce IL-12, which activates effector T cells CD4 + Th1 and CD8 + T [4,5]. Th1 produces IL-2 which further promotes proliferation of the effector T cells. Both CD4 + Th1 and CD8 + T cells kill cancer cells [6][7][8]. CD8 + T cells are more effective in killing cancer cells, but the helper function of CD4 + Th1 cells improves the efficacy of tumor-reactive CD8 + T cells [9]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 a mathematical model which combines a cancer vaccine with a checkpoint inhibitor: the vaccine (GVAX) increases the pool of T cells and the checkpoint inhibitor (anti-PD-1) enables the T cells to remain fully active in killing cancer cells.

Mathematical model
The mathematical model is based on the network shown in Fig 1. The list of variables, with units, is given in Table 1.

Fig 1. Interaction of immune cells with tumor cells.
Sharp arrows indicate proliferation/activation, blocked arrows indicate killing/blocking, and dashed lines indicate proteins on T cells. GM-CSF activates dendritic cells; activated dendritic cells produce IL-12; IL-12 activates naive CD4 + and CD8 + T cells; activated CD4 + T cells (Th1) produce IL-2 which induces proliferation of activated CD4 + and CD8 + T cells. Activated CD4 + and CD8 + T cells kill cancer cells. Activated CD4 + and CD8 + T cells express PD-1 and PD-L1, and cancer cells express PD-L1. The complex PD-1-PD-L1 inhibits the function of active CD4 + and CD8 + T cells.
https://doi.org/10.1371/journal.pone.0178479.g001 A anti-PD-1 concentration g/cm 3 We assume that the total density of cancer cells (C), active dendritic cells (D), CD4 + T cells (T 1 ) and CD8 + T cells (T 8 ) within the tumor remains constant in space and time: It is tacitly assumed that the debris of dead cells, including cancer cells undergoing necrosis or apoptosis, is quickly cleared from the tumor tissue. It is also tacitly assumed that the densities of immature dendritic cells and naive CD4 + and CD8 + T cells are constant throughout the tumor tissue. Under the assumption Eq (1), cancer cell proliferation, and migration of immune cells into the tumor give rise to internal pressure which results in cell movement, and we assume that all the cells move with the same velocity, u; u depends on space and time. We also assume that all the cells undergo diffusion, and that all the cytokines and drugs are diffusing within the tumor.
Equation for DCs (D). When cancer cells undergo necrosis, they release HMGB-1 [1]. We can model the dynamics of the necrotic cells and of HMGB-1 by the following equations: where λ N C C is the rate at which cancer cells become necrotic and λ HN C is the rate at which necrotic cells produce HMGB-1. We note that although molecules like HMGB-1, or other proteins, may be affected by the velocity u, their diffusion coefficients are several orders of magnitude larger than the diffusion coefficients of cells; hence their velocity terms may be neglected. The degradation of HMGB-1 is fast (*0.01/day) [31], and we assume that the removal of N C is also fast. We can therefore approximate the two dynamical equations by the steady state equations λ N C C C − d N C N C = 0 and λ HN C N C − d H H = 0, so that H is proportional to C, i.e., H = constant × C. Dendritic cells are activated by HMGB-1 [2,3]. Hence, the activation rate of immature dendritic cell D 0 is proportional to D 0 C K C þC , where the Michaelis-Menten law is used to account for the limited rate of receptor recycling time which occurs in the process of DCs activation. In the same way, GM-CSF, produced by the cancer vaccine, activates DCs at rate proportional to D 0 . Hence, the dynamics of DCs is given by where δ D is the diffusion coefficient and d D is the death rate of DCs. Equation for CD4 + T cells (T 1 ). Naive CD4 + T cells are activated by IL-12, and IL-2 induces proliferation of activated T 1 cells [4,5]. Both processes are assumed to be inhibited by the complex PD-1-PD-L1 (Q) [14], which reduces the production of T 1 cells by a factor 1 1þQ=K TQ .
Hence T 1 satisfies the following equation: where T 10 is the density of naive CD4 + T cells.

Equation for activated CD8 + T cells (T 8 ).
IL-12 activates CD8 + T cells and IL-2 induces the proliferation of CD8 + T cells [4,5]. Hence, similarly to the equation for T 1 , T 8 satisfies the equation where T 80 is the density of naive CD8 + T cells.

Equation for tumor cells (C).
Cancer cells are killed by T 1 and T 8 [6][7][8]. We assume a logistic growth with carrying capacity (C M ) in order to account for competition for space among the cancer cells. Hence, where η 1 , η 8 are the killing rates of cancer cells by T 1 and T 8 , and d C is the natural death rate of cancer cells.

Equation for GM-CSF (G).
We assume that GVAX is injected intradermally every 3 days for 30 days (as in mouse experiments [21]) providing a sourceĜðtÞ of GM-CSF, which we represent byĜ where γ G is the effective level of the drug; although the level of the drug varies between injections, for simplicity we take it to be constant. The concentration of GM-CSF in tissue is very small [32], and accordingly, we assume a low rate of constant source λ G for GM-CSF. Hence G satisfies the following equation: where d G is the degradation rate of GM-CSF.

Equation for IL-12 (I 12 )
. IL-12 is produced by activated DCs [4,5], so that @I 12 @t À d I 12 r 2 I 12 ¼ l I 12 D D |fflffl{zfflffl} production by DCs Equation for IL-2 (I 2 ). IL-2 is produced by activated CD4 + T cells (T 1 ) [4,5]. Hence, 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 [14,15]. We assume that the expression level of PD-1 is the same for activated CD4 + and CD8 + T cells. Hence, P is given by where ρ P is the ratio between the mass of one PD-1 protein to the mass of one T cell. The coefficient ρ P is constant when no anti-PD-1 drug is injected. In that case, to a change in T = T 1 + T 8 , given by @T @t , there corresponds a change of P, given by r P @T @t . For the same reason, r Á ðuPÞ ¼ r P r Á ðuTÞ and r 2 P = ρ P r 2 T when no anti-PD-1 drug is injected. Hence, P satisfies the equation Recalling Eqs (3) and (4) for T 1 and T 8 , we get @P @t þ r Á ðuPÞ À d T r 2 P ¼ r P ðl T 1 I 12 T 10 þ l T 8 I 12 T 80 Þ I 12 When anti-PD-1 drug (A) is applied, PD-1 is depleted (or blocked) by A. In this case, the ratio P T 1 þT 8 may change. In order to include in the model both cases of with and without anti-PD-1, we replace ρ P in the previous equation by P T 1 þT 8 . Hence, @P @t þ r Á ðuPÞ À d T r 2 P ¼ P where μ PA represents the rate at which P is depleted/blocked by A. PD-L1 is expressed on the surface of activated CD4 + and CD8 + T cells [14,15] and on cancer cells [15,16]. Hence, the concentration of PD-L1 (L) is proportional to (T 1 + T 8 ) and C: where ε depends on the specific type of tumor.
PD-L1 from T cells or cancer cells ligands to PD-1 on the plasma membrane of T cells, thus forming a complex PD-1-PD-L1 (Q) on the T cells [15,16]. Denoting the association and disassociation rates of Q by α PL and d Q , respectively, we can write so that The half-life of Q is less then 1 second (i.e. 1.16 × 10 −5 day) [33], and hence d Q is very large, and we may approximate the dynamical equation by the steady state equation,

Equation for anti-PD-1 (A).
We assume that anti-PD-1 is injected intradermally every 3 days for 30 days (as in mouse experiments [21]) providing a sourceÂðtÞ of anti-PD-1: where γ A is the effective level of the drug; although the level of the drug varies between injections, for simplicity we take it to be constant. The drug A is depleted in the process of blocking PD-1. Hence, where μ PA is the rate at which A is degraded in the process of blocking PD-1.

Equation for cell velocity (u):
We assume that most of the tumor consists of the extracellular matrix, ECM (approximately, 0.6 g/cm 3 ), and cancer cells (approximately, C = 0.4 g/cm 3 ), and that the densities of D, T 1 and T 8 are approximately 4 × 10 −4 , 2 × 10 −3 and 1 × 10 −3 g/cm 3 , respectively (as explained in the section on parameter estimation). We further assume that all cells are approximately of the same volume and surface area, so that their diffusion coefficients are the same. For definiteness, we take the constant in Eq (1) to be 0.4034. Adding Eqs (2)-(5), we then get To simplify the computations, we shall assume that the tumor is spherical and denote its moving boundary (i.e. its radius) by r = R(t). We also assume 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, and each of the equations of the form

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 the naive CD4 + T cells and naive CD8 + T cells which migrated from the lymph nodes into the tumor microenvironment have constant densi-tiesT 1 andT 8 at the tumor boundary, and that they are activated by IL-12 upon entering the tumor. We represent this process by the flux conditions at the boundary r = R(t): where s T ðI 12 Þ ¼ s 0 We impose the no-flux boundary conditions for all the remaining variables: No À flux for Dðr; tÞ; Cðr; tÞ; Gðr; tÞ; I 12 ðr; tÞ; I 2 ðr; tÞ; Pðr; tÞ and Aðr; tÞ; at r ¼ RðtÞ: It is tacitly assumed here that the receptors PD-1 and ligands PD-L1 become active only after the T cells are inside the tumor.
Initial conditions. We take initial values of all the variables to be constant. Later on we shall compare the simulations of the model with mouse experimental results, for 60 days. Accordingly, we take initial values whereby the average density of cancer cells has not yet increased to its steady state, 0.4 g/cm 3 , and, in view of Eq (1), the total density of the immune cells is initially above its steady state. We take (in unit of g/cm 3 ): We assume that initially A = 0, and G ¼ 2:61 Â 10 À 10 g= cm 3 ; I 12 ¼ 1:8 Â 10 À 10 g= cm 3 ; I 2 ¼ 4:74 Â 10 À 11 g= cm 3 ; P ¼ 11:2 Â 10 À 10 g= cm 3 : These values are close to the steady state values which are computed in the section on parameter estimation.

Results
The simulations of the model were performed by Matlab based on the moving mesh method for solving partial differential equations with free boundary [34] (see the section on computational method). All the computations are done in dimensionless form, but displayed in dimensional form.
The average density or concentration of a species is defined as its total mass in the tumor divided by the tumor volume. Fig 2 shows the average concentrations of all the species of the model over a period of 60 days in the control case, that is, when no drugs are administered; the parameter values are given in Tables 2 and 3. The radius of the tumor is increasing, from 0.01 cm to 0.0313 cm. The average density of the cancer cells is initially increasing and later it stabilizes while the densities of the immune cells are first decreasing and later stabilize. Correspondingly, the concentrations of the cytokines produced by the immune system also first decrease and later stabilize. Some of the parameters in the model were estimated by assuming the immune cells and cytokines are in steady state. Their steady states in Fig 2 approximately agree with those which we assumed in estimating the parameters, thus establishing the consistency of our assumed steady-state values.
We can also simulate the spatial distribution of each of the variables. Fig 3 shows the distribution of cancer cells and T cells (T 1 + T 8 ) at times t = 15, 30, 60 days. We see that the density of T cells increases toward the boundary; this is a result of the influx of T cells from the lymph nodes. Correspondingly, the density of cancer cells decreases toward the boundary. In our model, we tacitly assumed avascular conditions, since we wanted to focus primarily on the difference between control and treatment. Since, however, the tumor radius reaches approximately 313 μm at day 60, hypoxic conditions may actually reduce the cancer cells' density at the core of the tumor (both in control and treatment).
In Figs 2 and 3, we have taken ε = 0.01. Some cancer cells may express more or less PD-L1 [17,[35][36][37]. By decreasing or increasing ε, the radius of the tumor will decrease or increase, respectively, but the profiles of the densities/concentrations remain qualitatively the same (not shown here).
We next proceed to explore the effect of treatment with GM-CSF-secreting vaccine (GVAX) and anti-PD-1 drug. We are unable to make a direct connection between the levels of drugs administered to the patient, and their 'effective strengths' γ G and γ A in the model, since these data are not available. Based on the estimate of the concentration of GM-CSF in normal healthy tissue (see Parameter Estimations for Eq (6)), we chose the order of magnitude of GVAX to be 10 −10 (g/cm 3 Á day). The order of magnitude for anti-PD-1 drug (γ A ) is chosen so as to get the best agreement with the mouse experiments. In Fig 4(a), the 'effective strength' of the vaccine is given by γ G = 0.87 × 10 −10 g/cm 3 Á day and the 'effective strength' of the anti-PD-1 is given by γ A = 2 × 10 −10 g/cm 3 Á day. We see that, as a single-agent, anti-PD-1 is more effective than GVAX, and with the combined therapy the tumor radius is still increasing. This is in agreement with the mouse experiments reported in Fig 1-(A) (with melanoma) in [21], and Figs 3-(b) (with colon carcinoma) and 3-(c) (with melanoma) in [22]. T 1 density of CD4 + T cells from lymph node 4 × 10 −3 g/cm 3 estimated T 8 density of CD8 + T cells from lymph node 2 × 10 −3 g/cm 3 estimated γ G source of GM-CSF from the vaccine 1 × 10 −10 g/cm 3 Á day* estimated γ A source of anti-PD-1 1 × 10 −10 g/cm 3 Á day* estimated In Fig 4(b), we increased γ A by a factor 3/2. As a result, the tumor radius begins to decrease around day 50, even when administering anti-PD-1 as single-agent. This is in agreement with mouse experiments (with colon carcinoma) reported in Fig 1-(D) of [23].
In Fig 4(a), GM-CSF-secreting vaccine alone did not reduce the tumor radius as much as anti-PD-1 alone. However, if we increase the strength of the vaccine and decrease the anti-PD-1, taking for example, γ G = 3.84 × 10 −10 g/cm 3 Á day and γ A = 1 × 10 −10 g/cm 3 Á day, we then find that the vaccine decreases the tumor radius more than anti-PD-1 does; this is shown in Fig 5, and it is in agreement with mouse experiments (with melanoma) reported in Fig 1-(B) of [23].  The results of Figs 4 and 5 show that a combination therapy of GVAX and anti-PD-1 in appropriate amounts could significantly slow the growth of a tumor.
There is some uncertainty in the estimates of some of the parameters of the model. With somewhat different choices of these parameters, the simulation results will change quantitatively, and sensitivity analysis indicates the direction and intensity of the change (see the section on sensitivity analysis). In particular, the choice of γ G and γ A will affect the relative reduction in the growth of the tumor radius. In Figs 4 and 5 we made specific choices of γ G and γ A , in order to get simulations that agree with experimental results.
We next consider combination therapy for any values of GVAX and anti-PD-1. We define the efficacy of the combination therapy with (γ G , γ A ) by the formula: where R 60 = R 60 (γ G , γ A ) represents the tumor radius computed at day 60. If the tumor radius at day 60, R 60 (γ G , γ A ), is smaller than the radius in the control case, R 60 (0, 0), then the efficacy is a positive number, and its value is between 0 and 1 (or between 0% and 100%); the efficacy increases to 1 (or to 100%) when the tumor radius R 60 (γ G , γ A ) decreases to 0 by day 60. Fig 6  (a) is the efficacy map of the combined therapy with γ G in the range of 0 − 4.8 × 10 −10 ) g/cm 3 Á day and γ A in the range of 0 − 4 × 10 −10 g/cm 3 Á day. The color column in Fig 6(a) shows the efficacy for any pair of (γ G , γ A ); the efficacy is positive, and its maximum is 0.95 (95%). We see that an increase in either γ G or γ A improves the efficacy of the treatment. Early stage clinical trials consider the safety of each of the drugs, GVAX and anti-PD-1, separately. In the GVAX treatment of patients with pancreatic cancer, no dose-limiting toxicity or minimal toxicity was observed [38][39][40]. On the other hand, in treatment with anti-PD-1, mild to moderate toxicity (grades 1 and 2) was observed in melanoma [41,42]. In non-small-cell lung cancer, 14% of the patients were observed to have more severe toxicity (grades 3 and 4) [43]. Based on these observations we conclude that anti-PD-1 causes more toxicity than Combination therapy of cancer vaccine and immune checkpoint inhibitors GVAX. However, clinical trials on the safety and efficacy of the combined drugs are limited [23,44,45].
The amount of drug in clinical trials is constrained by the maximum tolerated dose (MTD). In combination therapy this constraint may depend on the proportion between the amounts of the drugs. We note that there is a large literature on the trade-off between efficacy and toxicity [46][47][48][49]. Here we consider, as an example, two treatments, ðg Ã G ; g Ã A Þ and ðg ÃÃ where both satisfy the MTD requirement. The question is then which of the two treatments is more effective in reducing the tumor volume. We can use the efficacy map to address such a question. We illustrate this in one special case.
We compare treatment of combination (γ G , γ A ) with monotherapy GVAX and monotherapy anti-PD-1. For monotherapy with GVAX, we take (1 + θ G )γ G , and for monotherapy with anti-PD-1 we take (1 + θ A )γ A , with θ A < θ G , to reflect the higher toxicity associated with anti-PD-1. If E(γ G , γ A ) is larger than both E((1 + θ G )γ G , 0) and E(0, (1 + θ A )γ A ), then we say that the synergy for the combination (γ G , γ A ) is positive, and otherwise, we say that the synergy is negative. More generally, we define the synergy σ = σ(γ G , γ A ) by the formula: Thus σ(γ G , γ A ) > 0 (positive synergy) if the combination (γ G , γ A ) reduces tumor growth more than either one of the single agents with (1 + θ G )γ G or (1 + θ A )γ A . Negative synergy occurs in the reverse case where instead of a combination therapy with (γ G , γ A ) we achieve better reduction of the tumor radius R 60 if we apply only one drug, (1 + θ G )γ G or (1 + θ A )γ A . The above concept of synergy is somewhat different from the usual definitions of synergy. For definiteness we take θ G = 1 and γ A = 0.5, but this choice, which is somewhat arbitrary, could be made more precise as more clinical data become available. Fig 6(b) is the synergy map for (γ G , γ A ) in the same range as in Fig 6(a); the color column shows the synergy σ(γ G , γ A ), with values that vary from -0.38 to 0.28.
We first note that the synergy is negative if γ G < 0.2 × 10 −10 g/cm 3 Á day. The reason is the following: if γ G is small then the numbers of T cells is also small, so instead of introducing a drug γ A which blocks the relatively small number of PD-1, it is more effective to increase the number of T cells by increasing γ G , i.e. E(2γ G , 0) > E(γ G , γ A ).
Next, for any fixed γ G , as seen in Fig 6(b), the synergy first increases with γ A and then decreases. In order to explain this occurrence, we note that with γ G fixed, the number of T cells that arrive into the tumor microenvironment is limited, and so is the number of their PD-1. Thus, in order to block the PD-1 there is a need for only a limited amount of anti-PD-1 drug; i.e. it is 'wasteful' to administer too much of γ A . We conclude that the maximum synergy is achieved when the amount of γ A is appropriately dependent on the amount of γ G , as indicated by the solid curve shown in Fig 6(b). Fig 7 shows that as the amounts of γ G and γ A increase the average density of T cells (T 1 + T 8 ) increases, and correspondingly the average density of cancer cells decreases. We also see that the density of T cells shown in the color column increases approximately by a factor of 6, whereas the density of cancer cells decreases approximately by a factor of 1.075. The proportion of these changes (i.e. 5/0.075) is similar to the proportion of densities of cancer cells to T cells.

Conclusion
The introduction of immune checkpoint inhibitors has been a very promising approach to cancer treatment. Blockage of the programmed death PD-1/PD-L1 is increasingly explored in single-agent studies [16,18]. However, because of the lack of tumor-infiltrating effector T cells, many patients do not respond to checkpoint inhibitor treatment as single-agent [18]. On the other hand, cancer vaccines have been shown to induce effector T-cells infiltration into tumors [19], although, to be fully effective, cancer vaccines have to overcome immune evasion [10]. It was recently suggested that the combination of a cancer vaccine and an immune checkpoint inhibitor may function synergistically to induce more effective antitumor immune response [18,20]. Clinical trials to test such combinations are currently ongoing [18,20].
In the present paper we developed a mathematical model to study the efficacy of the combination of a GM-CSF-secreting cancer vaccine (GVAX) and an anti-PD-1 drug, and to address the following question: at what proportion should the two drugs be administered in order to achieve the best efficacy when the amount of drugs is limited by MTD. The mathematical model is represented by a system of partial differential equations, based on the interactions among cancer cells, dendritic cells, CD4 + and CD8 + T cells, cytokines IL-12 and IL-2, GM-CSF, PD-1, PD-L1, PD-1-PD-L1 complex and anti-PD-1. Simulations of the model are shown to be in qualitative agreement with mouse experiments [21][22][23].
The cancer vaccine and the anti-PD-1 work in collaboration: The vaccine increases the number of tumor-infiltrating effector T cells and the anti-PD-1 ensures that these cells remain active. Given a fixed amount γ G of the vaccine and γ A of the anti-PD-1, it is important to estimate the level of synergy between these amounts in order to administer them in the most effective proportion. We introduced a specific concept of synergy σ(γ G , γ A ) and developed accordingly a synergy map in Fig 6(b). The map shows that if γ G is small, single-agent treatment (with γ A = 0) is the best treatment. Fig 6(b) also shows that, for any γ G , there is a unique value γ AG , such that the synergy increases as γ A increases as long as γ A remains smaller than γ AG , but the synergy decreases as γ A increases above γ AG ; the points (γ G , γ AG ) form the solid curve shown in Fig 6(b). We suggest that for optimal efficacy under MTD constraint, the level of dosage of anti-PD-1 (γ A ) should be related to the level of dosage of GVAX (γ G ) by setting γ A = γ AG , as indicated by the solid curve in Fig 6(b).
The mathematical model presented in this paper has several limitations: 1. In order to focus on the combined therapy of a cancer vaccine and an anti-PD-1 drug, we did not include some other important cells and cytokines that are found in the tumor microenvironment, such as T regulatory cells, macrophages, endothelial cells, and IL-10, IL-6 and TGF-β. We also did not include blood vessels and oxygen, thus assuming that the tumor is avascular. We tacitly assumed that the effect of these omissions is not significant in comparing the results of therapy to no therapy.
2. We assumed that the densities of immature, or naive, immune cells remain constant throughout the progression of the cancer and that dead cells are quickly removed from the tumor.
3. In estimating parameters we made a steady state assumption in some of the differential equations.
4. In the definition of synergy, and in the synergy map, we included in a crude way the fact that anti-PD-1 causes more toxicity than GVAX. Our aim was to develop a concept that will take account not only of efficacy but also of toxicity. For this reason we compared the treatment benefits for combination (γ G , γ A ) with the single agents (2γ G , 0) and (0, 1.5γ A ).

5.
We did not make any direct connection between drugs administered to the patient, and their 'effective strengths' γ G and γ A , as they appear in the differential equations, since these data are not available. The order of magnitude of GVAX (10 −10 ) was based on the estimate of the concentration of GM-CSF in normal healthy tissue (see Parameter Estimations for Eq (6) Clinical data on efficacy and toxicity are quite limited at this time. Our model should be viewed as setting up a computational framework, to address the question of optimal efficacy in combination therapy with cancer vaccine and checkpoint inhibitor.

Parameter estimation
Half-saturation. In an expression of the form Y X K X þX where Y is activated by X, the half-saturation parameter K X is taken to be the approximate steady state concentration of species X.

Sensitivity analysis
We performed sensitivity analysis, with respect to the tumor radius R at day 60, with respect to some of the production parameters of the System (2)-(13), namely, λ DC , λ DG , λ T 1 I 12 , λ T 1 I 2 , λ T 8 I 12 , λ T 8 I 2 , and the important parameters K TQ , η 1 and η 8 . Following the method in [64], we performed Latin hypercube sampling and generated 1000 samples to calculate the partial rank correlation coefficients (PRCC) and the p-values with respect to the tumor radius at day 60. In sampling all the parameters, we took the range of each from 1/2 to twice its values in Tables 2  and 3. The results are shown in Fig 8. Not surprisingly all the parameters are negatively correlated with the tumor radius. We note that the highest negatively correlated parameters are the activation rate of dendritic cells by cancer cells (λ DC ) and the inhibition of T cells activation by PD-1-PD-L1 complex (K TQ ). However, with different values of γ G and γ A the parameter λ DG can exceed λ DC .

Computational method
We employ moving mesh method [34] 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: @Dðr; tÞ @t ¼ d D DDðr; tÞ À divðuDÞ þ F; where F represents the term in the right hand side of Eq (2). Let r k i and D k i denote numerical approximations of i-th grid point and Dðr k i ; ntÞ, respectively, where τ is the size of time-step.