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 anticancer treatment with chemotherapy and targeted therapy [1][2][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]  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][18][19][20].
• 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.
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.

The model equations
The list of variables is given in Table 1 in unit of g/cm 3 ; 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 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 @D @t þ r � ðuDÞ |ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl } velocity À d D r 2 D |ffl ffl ffl {zffl ffl ffl } difusion ¼ l DCD0 C K C þ C |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } activation by necrotic cells whereD 0 is the density of immature dendritic cells which we assume to be constant, δ D is the dispersion (or diffusion) coefficients of DCs, and d D is the death rate of DCs.

Equation for tumor cells (C)
We assume a logistic growth for cancer cells with carrying capacity (C M ) 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:

PLOS ONE
where η are the killing rates of cancer cells by T, and d C is the natural death rate of cancer cells (by apoptosis).

Equation for IL-12 (I 12 )
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 A T a . 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 (Q 1 )
PD-1 is expressed on the surface of activated CD8 + T cells [41,42]. We denote by ρ T 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 ρ L the mass of all the PD-L1 in one T cell divided by the mass of the cell, and by ρ L ε C 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 ρ P is constant when no anti-PD-1 drug is administered. And in this case, to a change in T, given by @T @t , there corresponds a change of P, given by r P @T @t . For the same reason, r � (uP) = ρ P r �(uP) and r 2 P = ρ P r 2 P when no anti-PD-1 drug is injected. Hence, P satisfies the equation

PLOS ONE
TNF-α inhibitor reduces drug-resistance to anti-PD-1: A mathematical model Recalling Eq (3) for T, we get @P @t þ r � ðuPÞ À d T r 2 P ¼ r P ½RHS of Eq: ð3Þ�: When anti-PD-1 drug (A) is applied, PD-1 is depleted by A, and the ratio P T 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 . Hence, where μ PA 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 (Q 1 ) on the T cells [42,43]. Denoting the association and disassociation rates of Q 1 by α PL and d Q 1 , respectively, we can write The half-life of Q 1 is less then 1 second (i.e. 1.16 × 10 −5 day) [44], so that d Q 1 is very large. Hence we may approximate the dynamical equation for Q 1 by the steady state equation, where

Equation for TIM-3 (T M ), Gal-9 (G) and TIM-3-Gal-9 (Q 2 )
TIM-3 is expressed on the surface of activated CD8 + T cells. If we denote by ρ M 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 ρ G 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, K 0 (11), we see that resistance to anti-PD-1 is due to upregulation of TIM-3 [24], as represented by the parameter α MA , 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 γ A at days t 1 , t 2 , . . ., and set g A ðtÞ ¼ g A P t i <t e À bðtÀ t i Þ 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 A .
Hence A satisfies the following equation:

Equation for anti-TNF-α (B)
Similarly to Eq (14), we take the dynamics of B to be:

Equation for cells velocity (u)
As estimated in the section on parameter estimation, the average steady state densities of the immune cells D, T 8 and C are taken to be (in units of g/cm 3 ) 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 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 inactive CTL cells that migrated from the lymph nodes into the tumor microenvironment have constant densityT 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 s T ðI 12 Þ ¼ s 0 I 12 K I 12 þI 12 . We also assume no-flux for D; I 12 ; T a ; A and B at r ¼ RðtÞ; 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 P ¼ r P T at t ¼ 0: The last condition ensures that in the control case (where A � 0), the function P = ρ 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 X 0 . In estimating some production parameters (in section of parameter estimation) we made the assumption that each "half-saturation" parameter K X is equal to X 0 , and X 0 is estimated from clinical/experimental data. By comparing the values of K X in Table 2 with the values X 0 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 (I 12 , 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.

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][47][48] we shall evaluate the efficacy of a treatment by the tumor volume reduction rate (TVRR) at the end-point (final day) of  Table 2, for the mouse model.  K T half-saturation of T cells 1 × 10 −3 g/cm 3 [50] K C half-saturation of tumor cells 0.4 g/cm 3 [49] K I 12 half-saturation of IL-12 8 × 10 −10 g/cm 3 [50,53] K T a half-saturation of TNF-α 3 × 10 −11 g/cm 3 [54] K TQ 1 inhibition of function of T cells by PD-1-PD-L1 and V(0) = initial tumor volume, V(t e ) = end-point tumor volume. We note that the doses γ A and γ B 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 γ B is several orders of magnitude larger than γ A does not mean that it is relatively larger than γ A 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 γ A , γ B 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 γ A , γ B 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%. In

11
−1 × 10 −11 g/cm 3 � day and 0.1 × 10 −6 −1 × 10 −6 g/cm 3 � 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 (γ B ) 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 γ A is injected in day 1 of the second  Table 2, for the mouse model.

PLOS ONE
TNF-α inhibitor reduces drug-resistance to anti-PD-1: A mathematical model week of a cycle will occur only 2 weeks later when γ B 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 (γ A ) most effectively, the anti-TNF-α drug (γ B ) should be injected not too soon but definitely not too late after the injection of γ A .
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.
where M p is the molecular weight of p and A is a constant. In particular where M V and D V are respectively the molecular weight and diffusion coefficient of VEGF: M V = 24kDa [56] and D V = 8.64 × 10 −2 cm −2 d −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) Death/Degradation rates. From [58] we take d T = 0.18/d, and from [49] we take d D = 0.1/ d and d C = 0.17/d. The half-life of I 12 in peripheral blood is 352 minutes and in chord blood 614 minutes [51]. We accordingly take the half-life of I 12 in tissue to be 5 hours, so that d I 12 ¼ 1:38=d. The half-life of TNF-α is 4.6 minutes [52]; hence, d T a ¼ 216=d.
Production parameters. In order to estimate production parameters, we use the "steadystate" of each equation, that is, we equate the right-hand side of each equation to zero, replace each species X in this equation by X 0 and each K X by X 0 .

Sensitivity analysis
We performed sensitivity analysis on the production parameters in Eqs (2)-(6), the depletion rates m BT a and μ PA 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 A T a , λ Tα T are positively correlated and μ Tα B is negatively correlated. On the other hand, as T increases, cancer decreases. Hence λ DC , l TI 12 , λ I12 D and μ PA (in Fig 6) are negatively correlated. Finally the killing rate η of C by T is negatively correlated while the growth rate λ C of C is highly positively correlated. was partially supported by the Fundamental Research Funds for the Central Universities (19XNLG14), and the Research Funds of Renmin University of China, and the National Natural Science Foundation of China (No.11571364).