TGF-β inhibition can overcome cancer primary resistance to PD-1 blockade: A mathematical model

Immune checkpoint inhibitors have demonstrated, over the recent years, impressive clinical response in cancer patients, but some patients do not respond at all to checkpoint blockade, exhibiting primary resistance. Primary resistance to PD-1 blockade is reported to occur under conditions of immunosuppressive tumor environment, a condition caused by myeloid derived suppressor cells (MDSCs), and by T cells exclusion, due to increased level of T regulatory cells (Tregs). Since TGF-β activates Tregs, TGF-β inhibitor may overcome primary resistance to anti-PD-1. Indeed, recent mice experiments show that combining anti-PD-1 with anti-TGF-β yields significant therapeutic improvements compared to anti-TGF-β alone. The present paper introduces two cancer-specific parameters and, correspondingly, develops a mathematical model which explains how primary resistance to PD-1 blockade occurs, in terms of the two cancer-specific parameters, and how, in combination with anti-TGF-β, anti-PD-1 provides significant benefits. The model is represented by a system of partial differential equations and the simulations are in agreement with the recent mice experiments. In some cancer patients, treatment with anti-PD-1 results in rapid progression of the disease, known as hyperprogression disease (HPD). The mathematical model can also explain how this situation arises, and it predicts that HPD may be reversed by combining anti-TGF-β to anti-PD-1. The model is used to demonstrate how the two cancer-specific parameters may serve as biomarkers in predicting the efficacy of combination therapy with PD-1 and TGF-β inhibitors.

1 Model Equations 1 We assume that the combined densities of cells within the tumor remains constant in space and time: for some constant θ > 0. We assume that the densities of immature dendritic cells and 2 naive CD4 + and CD8 + T cells remain constant throughout the tumor tissue. Under 3 4 tumor, give rise to internal pressure which results in cells movement. We assume that 5 all the cells move with the same velocity, u; u depends on space and time and will be 6 taken in units of cm/day. We assume that cytokines and anti-tumor drugs are diffusing 7 within the tumor, and that also cells undergo diffusion (i.e., dispersion). 8

Equation for tumor cells (C)
We write the equation for C in the following form: Equation for dendritic cells (D) Dendritic cells are activated by the high mobility group box 1 (HMGB-1) expressed on necrotic cancer cells [1,2]. We assume that the density of HMGB-1 is proportional to the density of C. Hence, the activation rate of immature dendritic cells, with density D 0 , is proportional C K C +C , for some constant K C . The dynamics of dendritic cells is given by where δ D is the diffusion coefficient and µ D is the death rate of dendritic cells, and the 9 activation rate λ DC is a constant. 10 April 27, 2021 1/17 Equation for M1 macrophages (M 1 ) The equation for M1 macrophages has the following form:

14
T cells, but smaller, by a factor ε T , for Treg cells. If we denote by ρ P D the ratio 15 between the mass of the PD-1 proteins in one T cell to the mass of this cell, then The coefficient ρ P D is constant when no anti-PD-1 drug is administered. In this case, to 17 a change in T = T 1 + T 8 + ε T T r , given by ∂T /∂t, there corresponds a change of P D , 18 given by ρ P D ∂T /∂t. For the same reason, ∇ · (uP D ) = ρ P D ∇ · (uT ) and 19 ∇ 2 P D = ρ P D ∇ 2 T when no anti-PD-1 drug is injected. Hence, P D satisfies the equation: 20 Recalling Eqs. (6)-(8) for T 1 , T 8 and T r , we get When anti-PD-1 drug is applied, PD-1 is depleted at a rate proportional to A 1 , and, in this case, the ratio P D /(T 1 + T 8 + ε T T r ) may change. In order to include in the model both cases of with and without anti-PD-1, we replace ρ P D in Eq. (14) by P D /(T 1 + T 8 + ε T T r ). Hence, where µ P DA1 is the depletion rate of PD-1 by anti-PD-1. 21 We assume that the number of PD-L1 proteins in one T 1 cell is the same as for one T 8 cell, and denote by ρ P L the ratio of the mass of all the PD-L1 proteins in one T 1 cell to the mass of one cell. We assume that this ratio is smaller for the PD-L1 proteins on MDSC and cancer cells, by factors ε M and ε C , respectively, so that PD-L1 from T cells or cancer cells combines with PD-1 on the plasma membrane of 22 T cells, forming a complex PD-1/PD-L1 (Q) on the T cells [13,14]. Denoting the 23 association and disassociation rates of Q by α P D P L and µ Q , respectively, we write Since the half-life of Q is less than 1 second (i.e., 1.16 × 10 −5 day) [15], we may approximate the dynamical equation for Q by the steady state equation where σ = α P D P L /µ Q .

Equation for cells velocity (u).
We assume that all cells have approximately the same diffusion coefficient. Adding Eqs. (2)-(8) and using Eq, (1), we get

26
To simplify the computations, we assume that the tumor is spherical and that all the 27 densities and concentrations are radially symmetric, that is, functions of (r, t), is the boundary of the tumor, and that u = u(r, t)e r , 29 where e r is the unit radial vector. 30 Boundary conditions. We assume that the inactive CD4 + and CD8 + T cells that migrated from the lymph nodes into the tumor microenvironment have constant densitiesT 1 andT 8 , respectively, at the tumor boundary, and that they are activated by IL-12 upon entering the tumor. We then have the following conditions at the tumor boundary: We impose no-flux boundary condition on all the remaining variables: 31 No flux for M 1 , M 2 , D, T r , C, P, I 10 , I 12 , T β , P D , A 1 and A β at r = R(t).
It is tacitly assumed here that the receptors PD-1 and ligands PD-L1 become active 32 only after the T cells are already inside the tumor.

33
Initial conditions (at t = t 0 ). We take the initial values of cells (in units of g/cm 3 ) as follows: April 27, 2021 5/17 Then, by Eq. (1), θ = 0.632 in Eq. (17). 34 We assume that initially A 1 = 0 and A β = 0, and P = 1.3 × 10 −7 , I 10 = 2.5 × 10 −9 , I 12 = 1.4 × 10 −10 , Other nearby choices of initial conditions do not affect the simulations of the model 35 after a few days. 36 2 Parameter Estimates 37 Half-saturation. In an expression of the form the parameter K X is called the half-saturation of X. We denote by X 0 the average 39 concentration of X, and assume that to be not too close to 0 or to 1, and, for simplicity, take it to be 1/2, so that Estimate for K D . Vescovi et al. [16] measured the density of plasmacytoid Estimates for K I10 andK T I10 . The expression of IL-10 in human malignant 52 melanoma was investigated by Kruger-Krasagakes et al. [19]. They found that the level 53 of IL-10 ranged between 0.6 − 6 × 10 −9 g/cm 3 . We take 54 K I10 = 5 × 10 −9 g/cm 3 andK I10 = 5 × 10 −9 g/cm 3 .
April 27, 2021 6/17 Estimate ofK T Q and K Q . Denoting the association and disassociation rates of 63 the complex Q=PD-1/PD-L1 by α P D P L and d Q , respectively, we can write The half-life of Q is less than 1 second (i.e., 1.16 × 10 −5 day) [15], so that µ Q is very large. Hence, we may approximate the dynamical equation for Q by the steady-state equation where σ = α P D P L /µ Q . We can then write the inhibition of Th1 and CD8 + T cells by Q 65 in the form and take, as in [22], We assume that 68 K Q = 4.86 × 10 −20 g 2 /cm 6 .
Estimates for M 0 , D 0 and T 0 . We estimate the levels of macrophages, dendritic 69 cells and T cells in the mouse liver, and assume that these levels are similar in other 70 mouse organs. 71 We note that in mouse liver there are 1.5 × 10 7 macrophages [23], 10 7 T cells [24] 72 and 7 × 10 6 DCs [25]. The weight of a mouse liver is 2 grams [26,27], which corresponds, 73 approximately, to a volume of 2 cm 3 . Hence the density of macrophages in mouse liver 74 is 7.5 × 10 6 cells/cm 3 , the density of T cells is 5 × 10 6 cells/cm 3 , and the density of DC 75 is 3.5 × 10 6 cells/cm 3 .

83
We assume that the formula (26)  Estimate for λ C . The growth rate of human melanoma is λ h C = 0.616 d −1 [22]. We 91 assume that melanoma grows faster in mice than in humans, and take 92 λ C = 1.2λ h C = 0.739 d −1 .

138
We performed sensitivity analysis with a group of parameters, mostly production rates, 139 which were roughly estimated. The computations were done using Latin Hypercube

140
Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) with a Matlab package 141 by [50,51]. The range for the parameters in the sensitivity analysis were between ±50% 142 of their baseline values in Tables 1-2, except for λ C which was chosen between ±10% of 143 the baseline in Table 1.

144
The parameter λ C is clearly positively correlated, and so is the cancer-specific