TNF-α inhibitor reduces drug-resistance to anti-PD-1: A mathematical model.

Drug resistance is a primary obstacle in cancer treatment. In many patients who at first respond well to treatment, relapse occurs later on. Various mechanisms have been explored to explain drug resistance in specific cancers and for specific drugs. In this paper, we consider resistance to anti-PD-1, a drug that enhances the activity of anti-cancer T cells. Based on results in experimental melanoma, it is shown, by a mathematical model, that resistances to anti-PD-1 can be significantly reduced by combining it with anti-TNF-α. The model is used to simulate the efficacy of the combined therapy with different range of doses, different initial tumor volume, and different schedules. In particular, it is shown that under a course of treatment with 3-week cycles where each drug is injected in the first day of either week 1 or week 2, injecting anti-TNF-α one week after anti-PD-1 is the most effective schedule in reducing tumor volume.

Introduction

Drug resistance is one of the primary causes for suboptimal outcomes in cancer therapy [1]. Patients respond well at first, but relapse occurs for many cancer patients [2]. One common resistance is associated with the ATP-binding cassette (ABC) transporters. These protein pumps act to protect cells by ejecting variety of toxins, and they are overexpressed during anti-cancer treatment with chemotherapy and targeted therapy [1-3].

Other mechanisms of drug resistance are associated with mutations and epigenetic changes [4]. There are many mathematical models of drug resistance, particularly resistance to chemotherapy; e.g., a model of evolution of resistance [5], and a model of resistance associated with symmetric/asymmetric division of stem cells [6]. Article [7] reviews the mathematical models up to 2011, and very recent reviews are found in [8, 9]. Recent models considered multi-mutations in drug resistance for specific cases [10], and optimal therapy design to reduce drug resistance [11]. A list of the mathematical and computational methods used to simulate models of drug resistance are given in [12, 13]; in particular, articles [4, 14] use PDE models, as do the recent papers [15, 16] and the present one.

There are different mechanisms of resistance to different drugs. In this paper we consider resistance to the immunotherapy drug anti-PD-1, and show that injection of TNF-α will reduce the resistance.

PD-1 is a checkpoint on T cells. Its ligand PD-L1 is expressed on both T cells and cancer cells. The formation of the complex PD-1-PD-L1 initiates a signaling cascade that results in blocking the anti-cancer activity of T cells. PD-1 blockade by anti-PD-1 drugs, is currently increasing used in the treatment of cancers.

Internal mechanisms of drug resistance to anti-PD-1 include the following:

B2M mutation: Loss of b2-microglobulin (B2M) expression results in impaired cell surface expression of HLA class I (MHC class I), which in turn impairs antigen presentation to cytotoxic T cells, and thereby leads to anti-PD-1 resistance [17-20].

JAK-1/JAK-2 mutation: INF-γ released by T cells activates the signalling pathway. JAK1/JAK2-STAT1/STAT2/STAT3-IRF1, which leads to upregulation of PD-L1 on cancer cells [21]. Acquired resistance to PD-1 blockade immunotherapy in patients with melanoma was associated with defects in the pathways involved in interferon-receptor signaling, such as mutation of interferon-receptor-associated Janus kinase 1 (JAK1) or Janus kinase 2 (JAK2) [17, 20, 22]. Inability to respond to IFN-γ leads to reduced expression of PD-L1 on cancer and hence to reduced effectiveness of anti-PD-1.

Loss of neoantigen: Genetic mutations that make cells become cancer cells may result in production of proteins that immune system can recognize as antigen; these are called neoantigens. Anti-PD-1 treatment may cause mutations in cancer cells that result in loss of neoantigens, and, correspondingly, to a decrease in immune response, including decreases in the number of active anti-cancer T cells, and hence a decrease in the effectiveness of anti-PD-1 [23].

Another form of drug resistance is associated with a change undergoing in the anti-cancer CD8+ T cells (cytotoxic T cells, CTLs):

TIM-3 checkpoint: Under anti-PD-1 treatment, T-cell immunoglobulin mucin-3 (TIM-3) is upregulated on T cells [24]. Its ligand, Galectin-9 (Gal-9) is expressed on cancer cells, and the complex formed by the checkpoint TIM-3 and its ligand Gal-9 induces apoptosis of Th1 cells [25, 26], and, consequently, a reduction in CD8+ T cells whose proliferation depends on Th1-secreted IL-2.

Recent approaches to overcome anti-PD-1 drug resistance focus on combination strategies [27]. A prime example is the combination of anti-PD-1 and anti-CTL-4 [28-31].

TNF-α is a pleiotropic cytokine that is involved in diverse functions. In immune response to cancer it acts as immunosuppressive [32]. TNF-α has been shown to enhance the expression of PD-L1 on CD8+ T cells in cancer [33], including melanoma [34]. TNF-α elicits an increase in TIM-3+ CD8+ T cells, and anti-PD-1 triggers TIM-3 expression in TNF-α dependent manner [32]. Hence blockade of TNF-α overcomes resistance to anti-PD-1 by reducing TIM-3 and PD-L1 expression.

For simplicity, we shall combine CD4+ Th1 cells with CD8+ T cells, since they play similar roles in the response to anti-PD-1 drug resistance, and refer to the combined populaltion of T cells as cytotoxic lymphocytes (CTL).

In the present paper we develop a mathematical model of combination therapy with anti-PD-1 and anti-TNF-α. The model will be used to simulate the efficacy of anti-PD-1 under different amount of anti-TNF-α, and under different schedules of treatment.

The model includes the densities of cancer cells (C), dendritic cells (D) and CTL cells (T), the concentrations of cytokines TNF-α and IL-12, and the proteins PD-1, PD-L1, TIM-3 and Gal-9. Resistance to anti-PD-1 treatment is modeled by upregulation of TIM-3 [24]. The model is based on the network shown in Fig 1, and will be represented by a system of partial differential equations.

Fig 1

Interaction of immune cells with cancer cells.

Sharp arrows indicate proliferation/activation, blocked arrows indicate killing/blocking. C: cancer cells, D: dendritic cells, T: activated CTL cells. T cells express PD-1, PD-L1, TIM-3; cancer cells express PD-L1 and Gal-9.

The model equations

The list of variables is given in Table 1 in unit of g/cm3; the unit of time is taken to be 1 day. We assume that the total density of cells within three tumor remains constant in space and time:

Table 1

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

Notation	Description
C	density of cancer cells
D	density of dendritic cells
T	density of activated CTL cells
Tα	TNF-α concentration
I12	IL-12 concentration
P	PD-1 concentration
L	PD-L1 concentration
Q1	PD-1-PD-L1 concentration
Q2	TIM-3-Gal-9 concentration
TM	TIM-3 concentration
G	Gal-9 concentration
A	anti-PD-1 (e.g. nivolumab)
B	anti-TNF-α (e.g. infliximab)

Under the assumption (1), proliferation of cancer cells and immigration of immune cells into the tumor give rise to internal pressure among the cells, which results in cells movement. We assume that all the cells move with the same velocity, u; u depends on space and time and will be taken in units of cm/day. We also assume that all the cells undergo dispersion (i.e., diffusion), and that all the cytokines and anti-tumor drugs are also diffusing within the tumor.

Equations for DCs (D)

By necrotic cancer cells we mean cancer cells undergoing the process of necrosis. Necrotic cancer cells release HMGB-1 [35]. Dendritic cells are activated by HMGB-1 [36, 37], and we assume that the density of HMGB-1 is proportional to the density of cancer cells. Hence the dynamics of activated dendritic cells is given by where is the density of immature dendritic cells which we assume to be constant, δ is the dispersion (or diffusion) coefficients of DCs, and d is the death rate of DCs.

Equations for CTL cells (T)

Inactive T cells are activated by IL-12 [38, 39]. The processes of activation and proliferation of T cells are inhibited by the complexes PD-1-PD-L1 (Q1) and TIM-3-Gal-9 (Q2). We represents these inhibitions by factors and , respectively, and take (see Eq (10)) and, similarly, . Hence T satisfies the following equation: where is the density of the inactive T cells.

Equation for tumor cells (C)

We assume a logistic growth for cancer cells with carrying capacity (C) in order to account for competition for space among these cells. Cancer cells are killed by T cells. The equation for C takes the form: where η are the killing rates of cancer cells by T, and d is the natural death rate of cancer cells (by apoptosis).

Equation for IL-12 (I12)

The proinflammatory cytokine IL-12 is secreted by activated DCs [38, 39], so that

Since the diffusion coefficients of proteins are very large compared to those of cells, their advection velocity may be neglected.

Equation for TNF-α (T)

Macrophages are attracted to the tumor, and they produce TNF-α [40]. For simplicity, we do not include macrophages in the model, and represent their TNF-α contribution as a source . The cytokine TNF-α is also secreted by T cells, so that where B is the TNF-α inhibitor.

Equation for PD-1 (P), PD-L1 (L) and PD-1-PD-L1 (Q1)

PD-1 is expressed on the surface of activated CD8+ T cells [41, 42]. We denote by ρ the mass of all the PD-1 in one T cell divided by the mass of a T cell, so that

PD-L1 is expressed on the surface of activated CD8+ T cells [42] and cancer cells [42, 43]. We denote by ρ the mass of all the PD-L1 in one T cell divided by the mass of the cell, and by ρ ε the mass of all the PD-L1 in one cancer cell divided by the mass of the cell, so that

The expression of PD-L1 is upregulated by TNF-α [33], so that

The coefficient ρ is constant when no anti-PD-1 drug is administered. And in this case, to a change in T, given by , there corresponds a change of P, given by . For the same reason, ∇ · (uP) = ρ ∇ ·(uP) and ∇2 P = ρ∇2 P when no anti-PD-1 drug is injected. Hence, P satisfies the equation

Recalling Eq (3) for T, we get

When anti-PD-1 drug (A) is applied, PD-1 is depleted by A, and the ratio may change. In order to include in the model both cases of with and without anti-PD-1, we replace ρ in the previous equation by . Hence, where μ is the depletion rate of PD-1 by anti-PD-1.

PD-L1 from T cells or cancer cells combines with PD-1 on the plasma membrane of T cells to form a complex PD-1-PD-L1 (Q1) on the T cells [42, 43]. Denoting the association and disassociation rates of Q1 by α and , respectively, we can write

The half-life of Q1 is less then 1 second (i.e. 1.16 × 10−5 ) [44], so that is very large. Hence we may approximate the dynamical equation for Q1 by the steady state equation, , or where . Hence, .

Equation for TIM-3 (T), Gal-9 (G) and TIM-3-Gal-9 (Q2)

TIM-3 is expressed on the surface of activated CD8+ T cells. If we denote by ρ the ratio between the mass of all the TIM-3 proteins in one T cell and the mass of one T cell, then

The expression of TIM-3 is enhanced by TNF-α [32], and further upregulated by anti-PD-1 [24, 32], so that

Gal-9 is expressed on the surface of cancer cells. We denote by ρ the ratio between the mass of all the Gal-9 proteins in one cancer cell and the mass of a cancer cell, so that

Similarly to Eq (10), we represent the concentration of the complex TIM-3-Gal-9 in the form: Hence, .

From Eq (11), we see that resistance to anti-PD-1 is due to upregulation of TIM-3 [24], as represented by the parameter α, but the level of resistance depends on the concentration of TNF-α. This suggests that anti-TNF-α will reduce resistance to anti-PD-1.

Equation for anti-PD-1 (A)

We shall consider simple mice experiments where the anti-PD-1 drug is administered in equal amount γ at days t1, t2, …, and set for some β > 0. Some of the drug is depleted in the process of blocking PD-1 and some is degraded, or washed out at rate d. Hence A satisfies the following equation:

Equation for anti-TNF-α (B)

Similarly to Eq (14), we take the dynamics of B to be: where .

Equation for cells velocity (u)

As estimated in the section on parameter estimation, the average steady state densities of the immune cells D, T8 and C are taken to be (in units of g/cm3)

To be consistent with Eq (1), we take the constant on the RHS of Eq (1) to be 0.4014. We further assume that all cells have approximately the same diffusion coefficient. Adding Eqs (2)–(4), we get

To simplify the computations, we assume that the tumor is spherical and denote its radius by r = R(t). We also assume that all the densities and concentrations are radially symmetric, that is, they are functions of (r, t), where 0 ≤ r ≤ R(t). In particular, u = u(r, t) er, where er 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 inactive CTL cells that migrated from the lymph nodes into the tumor microenvironment have constant density at the tumor boundary, and, upon crossing the boundary, they are activated by IL-12. We then have the following flux conditions for T cells and P at the tumor boundary: where we take . We also assume the boundary condition for C is determined by Eq (1).

Initial conditions

The initial conditions must be consistent with Eqs (1) and (16), namely, and

The last condition ensures that in the control case (where A ≡ 0), the function P = ρ T is the solution of Eq (9) with the boundary condition (20). The specific choice of initial conditions for our simulations is

However, other choices give similar 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 [45] (see the section on computational method). Fig 2, in the control case (no drugs), shows the average densities of cells and concentrations of cytokines for the first 30 days, and the growth of tumor volume; We note that the average concentrations of each species, X, tends to steady state X0. In estimating some production parameters (in section of parameter estimation) we made the assumption that each “half-saturation” parameter K is equal to X0, and X0 is estimated from clinical/experimental data. By comparing the values of K in Table 2 with the values X0 from Fig 2 we see that this assumption is approximately satisfied. Fig 2 shows that, as tumor progresses, the average densities of the anti-cancer immune cells (D, T) increase, as do the average concentrations of the inflammatory cytokines (I12, T), while the tumor volume is increasing exponentially. We can also simulate the spatial variations of the variables, but we are interested, in this paper, only in the average profiles of the variables in order to determine the tumor volume growth (or shrinkage, under treatment). We note however that the diffusion coefficients of cells and cytokines differ by several orders of magnitude, and the boundary condition (19) brings new T cells into the tumor. For these reasons an ODE system cannot adequately represent the dynamics of the tumor volume.

Fig 2

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

All parameter values are the same as in Table 2, for the mouse model.

Table 2

Summary of parameter values.

Notation	Description	Value used	References
δD	diffusion coefficient of DCs	8.64 × 10−7 cm2 day−1	[49]
δT	diffusion coefficient of T cells	8.64 × 10−7 cm2 day−1	[49]
δC	diffusion coefficient of tumor cells	8.64 × 10−7 cm2 day−1	[49]
δ_I_12	diffusion coefficient of IL-12	7.17 × 10−2 cm2 day−1	estimated
δ_T_α	diffusion coefficient of TNF-α	8.46 × 10−2 cm2 day−1	estimated
δA	diffusion coefficient of anti-PD-1	4.76 × 10−2 cm2 day−1	estimated
δB	diffusion coefficient of anti-TNF-α	4.75 × 10−2 cm2 day−1	estimated
σ0	flux rate of T cells at the boundary	1 cm−1	[49]
η	killing rate of tumor cells by T cells	328.55 cm3/g · day	[50]
μPA	blocking rate of PD-1 by anti-PD-1	1.03 × 107 cm3/g · day	estimated
μTαB	blocking rate of TNF-α by anti-TNF-α	2.56 × 108 cm3/g · day	estimated
ρP	expression of PD-1 in T cells	2.49 × 10−7	[50]
ρL	expression of PD-L1 in T cells	3.25 × 10−7	[50]
ρM	expression of TIM in T cells	1.5 × 10−7	estimated
ρG	expression of Gal-9 in cancer cells	2 × 10−8	estimated
εC	the ratio of mass of PD-L1 in cancer cell and T cell	0.1	estimated
αL	upregulation rate of PD-L1 by TNF-α	1	estimated
αT	upregulation rate of TIM-3 by TNF-α	1	estimated
αMA	enhancement of upregulation of TIM-3 by anti-PD-1	2	estimated
λ_DCD^_0	activation rate of DCs by tumor cells	8 × 10−5 g/cm3 · day	[50]
λ_TI_12T^_0	activation rate of T cells by IL-12	6.48 × 10−4 g/cm3 · day	estimated
λC	growth rate of cancer cells in mice	1.295 day−1	estimated
λC	growth rate of cancer cells in humans	0.895 day−1	estimated
λI12D	production rate of IL-12 by DCs	2.76 × 10−6 day−1	estimated
λTαT	production rate of TNF-α by T cells	6.48 × 10−4 day−1	estimated
dD	death rate of DCs	0.1 day−1	estimated
dT	death rate of T cells	0.197 day−1	estimated
dC	death rate of tumor cells	0.17 day−1	[49]
d_I_12	degradation rate of IL-12	1.38 day−1	[51]
d_T_α	degradation rate of TNF-α	216 day−1	[52]
dA	degradation rate of anti-PD-1	0.046 day−1	estimated
dB	degradation rate of anti-TNF-α	0.069 day−1	estimated
KD	half-saturation of dendritic cells	4 × 10−4 g/cm3	[50]
KT	half-saturation of T cells	1 × 10−3 g/cm3	[50]
KC	half-saturation of tumor cells	0.4 g/cm3	[49]
K_I_12	half-saturation of IL-12	8 × 10−10 g/cm3	[50, 53]
K_T_α	half-saturation of TNF-α	3 × 10−11 g/cm3	[54]
K_TQ_1	inhibition of function of T cells by PD-1-PD-L1	1.36 × 10−18 g2/cm6	estimated
K_TQ_2	inhibition of function of T cells by Tim-3-Gal-9	1.365 × 10−18 g2/cm6	estimated
CM	carrying capacity of cancer cells	0.8 g/cm3	estimated
T^	density of CD8+ T cells from lymph node	2 × 10−3 g/cm3	estimated
β	drug control parameter	1.55 g/cm3	fitted

Bertrand et al. [32] demonstrated, in a mouse model, that TNF-α blockade overcomes resistance to anti-PD-1 in melanoma. In [32] mice were injected with 3 × 105 melanoma cells, and the tumor became detectible by day 6. Drugs (anti-PD-1 and anti-TNF-α) were injected at days 6,10,13 in (see Fig 6(d) of [32]) and at days 13, 16, 19 in (see Fig 6(e) of [32]); the tumor volume increase was in the range of 100-350 mm3.

Fig 3 shows the growth of tumor volume under treatment with anti-PD-1 (γ) and anti-TNF-α (γ). Taking γ = 10−10 g/cm3 · day−1, γ = 10−6 g/cm3 · day−1, and using the same schedules of treatment as in Bertrand et al. [32], Fig 3(a) is in qualitative agreement with Figure 6(d) of Bertrand et al. [32], and Fig 3(b) is in qualitative agreement with Figure 6(e) of [32].

Fig 3

Growth of tumor volume without treatment, or under treatment with γ, or combination (γ, γ).

(a) The treatment is at days 6, 10 and 13. (b) The treatment is at days 13, 16 and 19. γ = 10−10 g/cm3 · day, γ = 10−6 g/cm3 · day. All other parameter values are the same as in Table 2, for the mouse model.

Clinical trials in silico

We proceed with treatment of humans, given in cycle of 3 weeks. We compare the different scheduling options in a cycle:

(S1) Both anti-PD-1 (nivolumab) and anti-TNF-α (infliximab) are injected in the first day of week 1;

(S2) Anti-PD-1 is injected in the first day of week 1 and anti-TNF-α is injected in the first day of week 2;

(S3) Anti-TNF-α is injected in the first day of week 1 and anti-PD-1 is injected in the first day of week 2.

A course of cancer treatment may take 3 to 8 cycles. As in [46-48] we shall evaluate the efficacy of a treatment by the tumor volume reduction rate (TVRR) at the end-point (final day) of the treatment, where we define and V(0) = initial tumor volume, V(t) = end-point tumor volume.

We note that the doses γ and γ used in the model are correlated to, but not the same as, the dose amounts used in actual treatments of patients; in particular, the fact that in our model γ is several orders of magnitude larger than γ does not mean that it is relatively larger than γ in clinical treatment. But it is reasonable to expect that the reduction in the doses in the model will have the same effect on the tumor volume growth if the corresponding reduction was made in the treatment of patients. The doses of γ, γ in Fig 3 were chosen in order to fit the simulations with experimental results. For humans we took a smaller tumor cells growth, and smaller ranges of γ, γ in Figs 4 and 5 (below) in order to get the reduction of tumor volume to be in the same range as in clinical data; In [47], under two 3-week cycles in treatment of breast cancer, TVRR was >65%, and in [48], under four 2-week cycles treatment of rectal cancer, TVRR median was 52% but, for significant number of patients, it was >65%.

Fig 4

Efficacy map for three schedules: S1 (left), S2 (middle) and S3 (right).

The color columns indicate the tumor volume reduction rate (TVRR) after 5 cycles; tumor is initially with the radius R(0) = 0.5 cm (upper row) and R(0) = 1 cm (lower row).

Fig 5

Efficacy map for three schedules: S1 (left), S2 (middle) and S3 (right).

The color columns indicate the tumor volume reduction rate (TVRR) after 5 cycles (upper row) and 10 cycles (lower row); tumor is initially with radius R(0) = 1 cm.

In Fig 4 we simulated the TVRR for a tumor with initial radius R(0) = 0.5 cm (volume 0.5236 cm3), in the upper row, and R(0) = 1 cm (volume 3.1456 cm3) in the lower row, under 5 cycle treatment course, undergoing 5 cycles with doses of γ and γ in the ranges of 0.1 × 10−11−1 × 10−11 g/cm3 · day and 0.1 × 10−6−1 × 10−6 g/cm3 · day, respectively. We see that, in the treatment of both tumors, S2 is somewhat more effective than S1, while S1 is significantly more effective than S3. We can explain this as follows:

Under S1, resistance to anti-PD-1 begins at day 1 of a cycle, but its effect may be still small in the first few days. So it is more effective to apply the resistance-blocking drug (γ) in the first day of week 2 of a cycle rather than immediately in day 1 of week 1. On the other hand, under schedule S3, blocking anti-PD-1 resistance when γ is injected in day 1 of the second week of a cycle will occur only 2 weeks later when γ is injected in day 1 of week 1 of the following cycle. Hence the blockade of resistance comes “too late” and is not very effective. We conclude that, in order to overcome the resistance to anti-PD-1 drug (γ) most effectively, the anti-TNF-α drug (γ) should be injected not too soon but definitely not too late after the injection of γ.

In Fig 5 we repeated the simulations of Fig 4 with R(0) = 1 cm and a smaller range of doses; we did it for one treatment with 5 cycles and an extended treatment with 10 cycles. We first see that, once again, S2 is somewhat more effective than S1, and S1 is significantly more effective than S3. But we can also see from Fig 5 the benefits of extended treatment. For example, under schedule S2, with R(0) = 1 cm, the TVRR will increase from 75%, 85% and 90% to 96%, 97.8% and 99%, respectively, if the treatment is extended from 5 cycles to 10 cycles. By comparing the treatments with different ranges of dose, in the case of schedule S2 with R(0) = 1 cm, we see that if the range is decreased by a factor of 1/10 then the TVRR will decrease, for instances, from 98–99% to 82–85%. We also see from Fig 4 that under the sane treatment, the TVRR is larger for the tumor with the smaller initial volume. This can be explained by the fact that concentration of the same dose injected into a small tumor is larger than that injected into a larger tumor.

Conclusion

Drug resistance is a major obstacle to cancer treatment. Patients respond well for several months, or even years, but relapse eventually occurs. This is the case with both chemotherapy and immune therapy drugs. Anti-Cancer T cells have ‘brakes’ (called checkpoints) on their activities in the form of membranous proteins such as PD-1 and TIM-3; when PD-L1 (or Gal-9) combines with PD-1 (or TIM-3) the T cell activation is inhibited. A drug, anti-PD-1, is increasing used in the treatment of cancer, and in the present paper we considered the problem of how to deal with cancer resistance to anti-PD-1. We focused on the role of TNF-α in reducing this resistance.

TNF-α increases the production of PD-L1, which renders anti-PD-1 less effective. TNF-α also activates expression of TIM-3, which enables Gal-3, expressed on cancer cells, to block T cells by combing with TIM-3 on T cells membrane. These facts suggest that anti-TNF-α could be an effective drug in reducing resistance to anti-PD-1.

We developed a mathematical model where anti-PD-1 is combined with anti-TNF-α. The model is represented by a system of PDEs. The model is ‘minimal’ in the sense that it includes only cancer cells, T cells and dendritic cells, and the most relevant protein populations needed to address the effectiveness of the combined therapy. Simulations of the model were shown to be in agreement with mice experiments.

The model can be used to assess the efficacy of various protocols of treatment with the combination of the two drugs. As an example, we considered the following scheduling cycles.: Both drugs administered in day 1 of week 1 (S1), anti-PD-1 injected in day 1 of week 1 and anti-TNF-α injected in day 1 of week 2 (S2), or in the inverse order (S3). For each pair of (anti-PD-1, anti-TNF-α) in some range we marked on the color columns the tumor volume reduction rate (TVRR) at the end of the treatment. In Figs 4 and 5 it was shown that schedule S2 is somewhat more beneficial than schedule S1, while schedule S1 is significantly more beneficial than schedule S3; these results mean that in order for the blockade of the drug resistance to anti-PD-1, anti-TNF-α should not be injected too late after the injection of anti-PD-1, but also not too early. This result should be viewed as an hypothesis, to be validated in clinical trials.

Figs 4 and 5 also show how TVRR increases as the dose level increases when the range of doses increases by a factor of 10, or when the treatment, with the same dose range, is given to a tumor with a smaller volume.

The model developed in this paper can be used also be explore the efficacy of schedules with shorter cycles when such clinical trials become more prevalent.

Supplementary materials

Computational method

We consider the following convection-diffusion equation to illustrate our computational method: where F accounts for all the right-hand side terms. Since the model we consider is a free boundary problem, we employ the moving mesh method to compute it. We write Eq (25) can be written in the total derivative form

Let , denote numerical approximations of i-th grid point and , respectively, when t = nτ, where τ is the time stepsize. The discretization is derived by the explicit Euler finite difference scheme, i.e., where , , and , . Then the mesh is moving by , where is solved by the velocity equation. In order to make the Euler method stable, we take .

Parameter estimation

Half-saturation

In an expression which represents an activation of a species Y by a species X, the parameter K is called the “half-saturation” of X. We assume that the average density/concentration of each species, X, tends to steady state X0, and that the quotient X0/(K + X0) is not too close to 0 or to 1, and take

The values X0 will be estimated from clinical/experimental data, and then the K are determined by (26).

Diffusion coefficients

Young [55] established the following formula for the diffusion coefficient, D, of any protein p: where M is the molecular weight of p and A is a constant. In particular where M and D are respectively the molecular weight and diffusion coefficient of VEGF: M = 24kDa [56] and D = 8.64 × 10−2 cm−2d−1 [57]. The following table lists the molecular weight of the proteins in our model, taken from [56], and their corresponding diffusion coefficients computed from the formula (28): The molecular weights of nivalumab (A) and infliximab (B) are 143.6 kDa and 144.2 kDa, respectively. Hence, by formula (28), we get DA = 4.76 × 10−2 cm2d−1 and DB = 4.75 × 10−2 cm2d−1.

Diffusion coefficients of cells

For simplicity, we assume that all cell types have the same diffusion coefficient, and as in [57], we take

Death/Degradation rates

From [58] we take d = 0.18/d, and from [49] we take d = 0.1/d and d = 0.17/d. The half-life of I12 in peripheral blood is 352 minutes and in chord blood 614 minutes [51]. We accordingly take the half-life of I12 in tissue to be 5 hours, so that . The half-life of TNF-α is 4.6 minutes [52]; hence, .

Inhibition of T by Q1, Q2

From the notation , we deduce, as in [59], that . We assume a similar inhibition parameter for the checkpoint T, taking .

Production parameters

In order to estimate production parameters, we use the “steady-state” of each equation, that is, we equate the right-hand side of each equation to zero, replace each species X in this equation by X0 and each K by X0.

Eq (2). In steady state, we have

The number of dendritic cells in advanced melanoma can increase to 400 cells/mm3 [60]. Taking it to be 200 cells/mm3 and assuming that the mass of one DC is 5 × 10−10 g, we get D0 = 4 × 10−4. Since d = 0.1/day [49], we obtain

Eq (3). we assume that in steady state,

Hence, in steady state of Eq (3),

The number of CD8+ T cells in melanoma is in the range of 7 − 10 × 105 cells/cm3 [61]. Taking the mass of one cells to be 5 × 10−10, we estimate the density of CD8+ T cells at 5 × 10−4 g/cm3. Since T includes Th1 cells, we take K = T0 = 10−3 g/cm3. Since d = 0.18/day [58], we get g/cm3/day.

Eq (4). We take, in steady state, C = 0.4 g/cm3 [49] and C = 0.8 g/cm3, η = 328 0.8 g/cm3/g · day−1 [50], d = 0.17/day [58]. Hence, in steady state, so that λ = 0.996/day. Since, however, the tumor volume is increasing on account of the abnormal proliferation of cancer cells, we increase this value by a factor, taking λ = 1.2×0.996 = 1.295/day. This value will be used in mouse models, as in Figs 2 and 3. But for humans (Figs 4 and 5) it is assumed that the tumor cells proliferation is slower, taking λ = 0.885/day.

Eq (5). From the steady state equation we get

Since K = 4 × 10−4 g/cm3, [50, 53], and [51], we get λ = 2.76 × 10−6/day.

Eq (6). From the steady state equation in the control case we get

We take g/cm3day. Since [54] and [52], we get λ = 1.08 × 10−6/day.

Remaining parameters. We take ρ = 2.49 × 10−7 [50]. Then, in steady state with T = K = 10−3 g/cm3 and P = 2.49 × 10−10 g/cm3. We take ρ = 3.25 × 10−7 g/cm3 [50] and assume that the expression of PD-L1 in cancer cells is less than it is on T cells, taking ε = 0.1 in Eq (8). We also take α = 1 which means that TNF-α increases the expression of PD-L1 by 50%, in steady state.

We assume that TIM-3 is expressed on T cells at a lower level than PD-1, that is, ρ < ρ, and take

We also take (in Eq (11)) α = 1 and α = 2. We similarly assume that Gal-9 is expressed in cancer cells at a lower level than PD-L1 (in Eq (8)), that is, ρ < ρ ε, and take ρ = 2 × 10−8.

The half-life od nivolumab is 10-20 days. Taking it to be 15 days, we get

The half-life of infliximab is 7-12 days. Taking half-life of 10 days, we get

We assume that 10% of A is used in blocking PD-L1 while the remaining 90% degrades naturally or decreased by washout, that is

Recalling that, in steady state, with P = 2.49 × 10−10, g/cm3, we get

Under a similar assumption on the depletion of B, in steady state

Sensitivity analysis

We performed sensitivity analysis on the production parameters in Eqs (2)–(6), the depletion rates and μ in Eqs (6) and (9), and η in Eq (4). Following the method in [62], we performed Latin hypercube sampling and generated 1000 samples to calculate the partial rank correlation (PRCC) and the p-values with respect to the radius of the tumor at day 30. The results are shown in Fig 6 (The p-value was <0.01). Note that as T increases, cancer increases. Hence , λ are positively correlated and μ is negatively correlated. On the other hand, as T increases, cancer decreases. Hence λ, , λ and μ (in Fig 6) are negatively correlated. Finally the killing rate η of C by T is negatively correlated while the growth rate λ of C is highly positively correlated.

Fig 6

Sensitivity analysis of the PRCC values with respect to the radius of the tumor.

X	I12	Tα
Molecular Weight(kDa)	37.2	25.6
DX (10−2 cm2d−1)	7.17	8.46