ANTAGONISM AND NEGATIVE SIDE-EFFECTS IN COMBINATION THERAPY FOR CANCER

. Most clinical trials with combination therapy fail. One of the reasons is that not enough forethought is given to the interaction between the diﬀerent agents, as well as the potential negative side-eﬀects that may arise in the combined therapy. In the present paper we consider a generic cancer model with combination therapy consisting of chemotherapy agent X and checkpoint inhibitor A . We use a mathematical model to investigate the results of injecting diﬀerent amounts γ X of X and γ A of A . We show that there are some regions in the ( γ A ,γ X )-plane where as increase in γ X or γ A actually decreases the tumor volume; such ‘regions of antagonism’ should be avoided in clinical trials. We also show how to achieve the same level of tumor volume reduction with least negative-side eﬀects, where the side-eﬀects are represented by the level of inﬂammation of the tumor microenvironment.


1.
Introduction. Most cancer clinical trials with combination therapy that are successful in phase II fail in phase III [19]. One of the reasons for this failure is the fact that not sufficient forethought is given to the interaction between the different agents [28]. Interaction between two different anti-cancer agents in combination therapy may not be completely positive in the sense that increasing one of the drugs may actually result in an increase of the tumor volume, thus decreasing the efficacy of the treatment [7,15]. Clinical trials may fail if they show no significant efficacy. They may also fail because of unacceptable negative side-effects. In [16] it was shown how one can reduce negative side-effects without reducing the efficacy.
In the present paper we introduce a model of combination therapy and show how to avoid 'zones of the antagonism' between the drugs while, at the same time, reducing the inflammation without sacrificing efficacy.
The second agent is a checkpoint inhibitor. Anti-cancer T cells have programmed death protein 1 (PD-1) receptors. Activated T cells, as well as some cancer cells, express ligand PD-L1 which, when combined with PD-1, it initiates a signaling cascade that leads to inactivation of the T cells. The second agent in our model is PD-1 antibody, for instance, pembrolitumab, which was a recently approved by FDA; we shall refer to it as anti-PD-1 or PD-1 inhibitor. This drug blocks the formation of the complex PD-1-PD-L1 and thus keeps the T cells active.
We take the negative side-effect to be the increase in inflammation within the tumor microenvironment, and we represent the inflammation by the concentration of TNF-α.
We designate by γ A the injected amount of anti-PD-1 per unit time and by γ X the injected amount of chemotherapy drugs per unit time. We denote by V (γ A , γ X ) the tumor volume at some prescribed time, e.g., 60 days from the beginning of the treatment.
We address the following question: (i) Is V (γ A , γ X ) monotone decreasing in both variables for all pairs γ A , γ X in the (γ A , γ X )-plane, and if not, what are the "zones of antagonism", that is, the regions where V (γ A , γ X ) increases as γ A or γ X is increased? (ii) Given any equi-efficacy curve in the (γ A , γ X )-plane (i.e., a curve on which V (γ A , γ X ) is constant), how to choose a pair (γ A , γ X ) on this curve so that TNF-α is at its minimum? 2. Mathematical model. The variables of the model are listed in Table 1 in unit of g/cm 3 . We note that CD8 + T cells (also called cytotoxic T cells) kill cancer cells; they are helped by a subclass of CD4 + T cells, called Th1 cells. For simplicity we combine both cells, and refer to them as effective T cells, or briefly, T cells. The mathematical model is based on the network shown in Fig. 1. The model will be represented by a system of partial differential equations with constant parameters. The description and values of these parameters are given in Tables 2 and  3.

Notation Description Notation Description
We assume that the total density of all cells is constant within the tumor, Sharp arrows indicate proliferation/activation, blocked arrow indicates killing/blocking, inverted sharp arrow indicates recruitment/chemoattraction. C: cancer cells, T : effector T cells.
As the tumor evolves and cancer cells proliferate, cells are pushing against each other. The resulting pressure among the cells give rise to velocity u, measured in unit of cm/day.

Equation for
DCs (D). By necrotic cancer cells we mean cancer cells undergoing the process of necrosis. Necrotic cancer cells release HMGB-1 [31]. We model the dynamics of the necrotic cells (N C ) and HMGB-1 (H) 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 order of magnitude larger than the diffusion coefficients of cells, hence their velocity term may be neglected. The degradation of HMGB-1 is fast (∼0.01/day) [37], and we assume that the process of necrosis is also fast. We may then approximate the two dynamical equations by the steady state equations Dendritic cells are activated by HMGB-1 [25,27]. Hence, the activation rate of immature dendritic cells, with density D 0 , is proportional to D 0 H K H +H , or to D 0 C K C +C since H is proportional to C. Here, the Michaelis-Menten law was used to account for the limited rate of receptor recycling time which takes place in the process of DCs activation. Chemotherapy treatment with some drugs, for example, Dolastatin 10, and Ansamitocin [22], induce immunogenic cell death that make dying cells present more antigen to DCs. Thus the DC-activation term should include the contribution of cancer cells dying from the chemotherapy drug. We assume that the density of these dying cells is proportional to the concentration X of the chemotherapy drug, and represent the corresponding enhancement in the DC-activation by a factor 1 + ε DX X K X +X . The equation for DC is then given by where δ D is the diffusion coefficient and d D is the death rate of DCs. Equation for effector T cells (T ). Naive T cells differentiate into active effector T cells (T ) under IL-12 (I 12 ) environment [12,18], while IL-10 (I 10 ) inhibits the differentiation of naive T cells into effector T cells [26,35]. The activation and proliferation of T are assumed to be inhibited by PD-1-PD-L1 (Q), by a factor 1 1+Q/K T Q . Hence T satisfies the following equation: activation by IL-12 · 1 1 + I 10 /K T I 10 inhibition by IL-10 where T 0 is the density of the naive T cells.
Equations for M1 macrophages (M 1 ) and M2 macrophages (M 2 ). The equation for M1 macrophages takes the following form: where we used the notation X + = X if X > 0, and X + = 0 if X ≤ 0. The first term on the RHS represents a source of macrophages differentiated from monocytes which are activated by MPC-1 (M P ) and the second term represents the chemoattraction of macrophages to MCP-1 [24]. The third and forth terms on the RHS are phenotype changes from M1 to from M2 under the influence of TGF-β, and M2 to M1 under the influence of TNF-α (T α ) [2,24]. Similarly,

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. The proliferation rate depends on the density of oxygen (W ) [2]. Docetacel (X) reduces proliferation of tumor cells. Cancer cells are killed by T cells.

ANTAGONISM AND NEGATIVE SIDE-EFFECTS 2241
The equation for C takes the form: where η is the killing rates of cancer cells by T , d C is the natural death rate of cancer cells, and 1 1+X/K CX represents the inhibition in cancer cells growth under the treatment with the chemotherapy agent.
Equation for IL-12 (I 12 ). The proinflammatory cytokine IL-12 is secreted by activated DCs [12,18] and by M1 macrophages [24], so that production by DCs and M1 − dI 12 I12. degradation (7) Here, we included increased production of I 12 by D under the influence of X [22,34]. Equation for TGF-β (T β ). TGF-β is produced by tumor cells [26] and M2 macrophages [1,6]: production by cancer and M2 Equation for IL-10 (I 10 ). IL-10 is produced by cancer cells and by M2 macrophages [26]. Hence I 10 satisfies the following equation: production by cancer and M2 − dI 10 I10. degradation (9) Equation for TNF-α (T α ). TNF-α is produced by primarily M1 macrophages and T cells [3,38]. Hence Equation for MCP-1 (M P ). MCP-1 is produced by cancer cells and by M2 macrophages [20,33], so that production by cancer and M2 − dM P MP . degradation (11) Equation for PD-1 (P ), PD-L1 (L) and PD-1-PD-L1 (Q). PD-1 is expressed on the surface of activated T cells [5,29]. If we denote by ρ P the ratio between the mass of one PD-1 protein to the mass of one T cell, then PD-L1 is expressed on the surface of activated T cells [29], M2 macrophages [23,29], and cancer cells [23,29]. We denote the ratio between the mass of one PD-L1 protein to the mass of one T cell by ρ L . Then where ε M and ε C depends on the specific type of tumor.
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 ρ P ∂T ∂t . For the same reason, ∇ · (uP ) = ρ P ∇ · (uT ) and ∇ 2 P = ρ P ∇ 2 T when no anti-PD-1 drug is injected. Hence, P satisfies the equation Anti-PD-1 drug (A) may change the ratio P T . 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, depletion by anti-PD-1 (13) where µ P A is the depletion rate of PD-1 by anti-PD-1. PD-L1 from T cells, M2 macrophages or cancer cells combines with PD-1 on the plasma membrane of T cells, to form a complex PD-1-PD-L1 (Q) on the T cells [23,29]. Denoting the association and disassociation rates of Q by α P L and d Q , respectively, we can write The half-life of Q is less then 1 second (i.e. 1.16 × 10 −5 day) [21], so that d Q is very large. Hence we may approximate the dynamical equation for Q by the steady state equation, α P L P L = d Q Q, or where σ = α P L /d Q .

Equation for cells velocity (u).
We assume that the densities of cells in the growing tumor tend to steady state, and take the density of cancer cells in steady state to be 0.4g/cm 3 . We also assume that most of the macrophages are "tumor associated macrophages" (TAM) [32], which we identify as M2 macrophages.
We take, as in the parameter estimation section, the steady state densities of the immune cells to be in unit of g/cm 3 , and C = 0.4 g/cm 3 . We assume that all cells have approximately the same diffusion coefficient. Adding the equations of all the cells, we get where the constant 0.4067 follows from Eqs.(1) and (14). Equation for anti-PD-1 (A). For simplicity, we assume that the level of the injected drug is constant, and denote by γ A the "effective level" of the injection The drug A is depleted in the process of blocking PD-1 (P ). Hence, Equation for chemotherapy agent (X). We assume a constant level of injection and denote by γ B its effective level. We represent the rate of absorption by cancer cells by µ XC C X K X +X , and hence obtain the following equation for X: (18) 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, functions of (r, t), where 0 ≤ r ≤ R(t). In particular, u = u(r, t)e r , where e r is the unit radial vector.
Equation for free boundary. We assume that the free boundary r = R(t) moves with the velocity of cells, so that where u is determined indirectly by Eq. (16). Boundary conditions. We assume that the inactive T cells that migrated from the lymph nodes into the tumor microenvironment have constant densitiesT at the tumor boundary, and that they activated by IL-12 upon entering the tumor. We then have the following flux condition at the tumor boundary: where we take σ T (I 12 ) = σ 0 I12 K I 12 +I12 . We impose a no-flux boundary condition for all the remaining variables: No-flux for D, M1, M2, C, I12, T β , I10, Tα, MP , P, A, and X 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 already inside the tumor.
3. Numerical results. The simulations of the model were performed by Matlab based on the moving mesh method for solving partial differential equations with free boundary as in [16]. Fig. 2 shows the profiles of the average densities/concentrations of all the variables of the model in the control case (no drugs), and the tumor volume over 30 days; the parameter values are as in Tables 2 and 3. Some of these parameters were estimated by assuming the convergence of the averages to steady state. The nearly steady states in Fig. 2 are consistent with these assumed in the parameters estimation. We see that the average density of cancer cells is continuously increasing, while the average density of the anticancer T cells is decreasing. The average density of dendritic cells is also decreasing after a few days, and so is the average density of the cytokine I 12 which activates the T cells and is produced by dendritic cells and M 1 macrophages. The initial increase or decrease in some of the species is a result of the choice of the initial conditions which was somewhat arbitrary, but this does not affect the dynamics of the average densities after a few days. It is interesting to note that the tumor volume increases first slowly, but later on exponentially. Fig. 3 shows two examples where the anti-PD-1 is more effective than the chemotherapy in reducing cancer volume (Fig. 3(a)), or less effective (Fig. 3(b)).
We next consider the efficacy of the combination therapy for a range of values of chemotherapy agent and anti-PD-1. We define the efficacy of a combination  Tables 2 and 3, for a mouse model.  Figure 3. Growth of tumor volume under treatment with γ X or γ A , or combination (γ X , γ A ). The chemotherapy or/and anti-PD-1 treatments. (a) γ X = 5 × 10 −13 g/cm 3 · day, γ A = 2 × 10 −11 g/cm 3 · day; (b) γ X = 3 × 10 −13 g/cm 3 · day, γ A = 5 × 10 −11 g/cm 3 · day. All other parameter values are the same as in Tables  2 and 3, for a mouse model. therapy, with (γ X , γ A ), by the formula: where the tumor volume V 30 = V 30 (γ X , γ A ) is computed at day 30; V 30 (0, 0) is the tumor volume at day 30 in the control case. The efficacy values lie between 0 and 1 (or between 0% and 100%). Fig. 4(a) is an efficacy map for γ X in the range of 0 − 2 × 10 −12 g/cm 3 · day and γ A in the range of 0 − 1.8 × 10 −10 g/cm 3 · day. We see, in the southwest corner that the two drugs are positively correlated, that is, the efficacy increases if γ X or γ A is increased. Hence, every set of points (doses) (γ X , γ A ) of equal efficacy E can be represented by a graph γ A = f E (γ X ) where f E (s) is a monotone decreasing function, and f E1 (s) > f E2 (s) if E 1 > E 2 . This is no longer the case in the northeast corner of the efficacy map, as seen in Fig. 4(b), in an area around a northwest/southeast diagonal. We refer to this area as a "zone of antagonism." In this area colors representing different levels of efficacy are seen protruding into each other. This means that, in this area, if one drug is increased then the tumor volume may actually increase. This antagonism between the two drugs is also reflected in the density of the T cells, as seen in Fig. 5, where the peninsula-like region of red color is surrounded on both sides by yellow-green-blue colors.
Treatment with (γ X , γ A ) increases the inflammation, which is represented by the proinflammatory cytokine TNF-α, and this may increase the risk of side effects, such as inflammatory bowel disease, psoraisis and rheumutoid arthritis. Fig. 6 shows the level of TNF-α for each pair (γ X , γ A ) in the same range as in Fig. 4(a). By comparing the two figures we conclude that, roughly, in order to achieve the same efficacy E with the smallest level of TNF-α we should take the smallest γ X from among all the pairs (γ X , γ A ) with equal efficacy. We note, however, that taking the smallest γ X , or the largest γ A , may cause other negative side-effects which are not considered here. Figure 4. Efficacy of combination therapy at day 30 for different pair of (γ X , γ A ). Here (a) γ X = 0 − 2 × 10 −12 g/cm 3 and γ A = 0 − 1.8 × 10 −10 g/cm 3 . (b) γ X = 0.9 × 10 −12 − 2 × 10 −12 g/cm 3 and γ A = 1 × 10 −10 − 1.8 × 10 −10 g/cm 3 . All other parameter values are the same as in Tables 2 and 3. 4. Conclusion. Before designing a cancer clinical trial with combination therapy, we need to gain a good understanding how the agents to be administered will interact with each other, and what will be their impact on negative side effects. In the present paper this was illustrated in the case where one drug is anti-PD-1, an immunotherapy drug, and the other one is a chemotherapy drug, which activates dendritic cells. Denoting by γ X and γ A the amount, per time, of injected doses of the chemotherapy and anti-PD-1, we developed a mathematical model by a system of PDEs, and used it to simulate an efficacy map for a range of γ X and γ A . We found that the two drugs are "mostly" positively correlated in the sense that, when each Figure 5. Average T cell density at day 30 for different pair of (γ X , γ A ). Here γ X = 0 − 2 × 10 −12 g/cm 3 and γ A = 0 − 1.8 × 10 −10 g/cm 3 . All other parameter values are the same as in Tables 2 and  3. Figure 6. Average concentration of TNF-α at day 30 for different pair of (γ X , γ A ). Here γ X = 0 − 2 × 10 −12 g/cm 3 and γ A = 0 − 1.8 × 10 −10 g/cm 3 . All other parameter values are the same as in Tables 2 and 3. of them is increased, the tumor volume decreases. However, there are (γ X , γ A )regions where the two agents are antagonistic, that is, an increase in one of the agents actually increases tumor volume. We also found that in order to decrease the negative side effects of inflammation, γ X should generally be taken as small as possible on the equi-efficacy curve.
The results of this paper suggest that the (γ X , γ A ) to be chosen in clinical trials, for our particular cancer model, should avoid the areas where antagonism between the drugs exists and, at the same time, should have the smallest feasible γ X on equi-efficacy curve in order to minimize the inflammation.
The method and concepts developed in this paper should be applicable to other clinical trial with combination therapy.

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.
Diffusion coefficients. By [36], we have the following relation for estimating the diffusion coefficients of a protein p: where M G and δ G are respectively the molecular weight and diffusion coefficient of VEGF, M p is the molecular weight of p, and M G = 24kDa [30] and δ G = 8.64×10 −2 cm 2 day −1 [17]. Since, M A = 32kDa [15], we get δ A = 4.73 × 10 −2 cm 2 day −1 . We take δ X = 0.27 cm 2 day −1 .
Eq. (6). We take d C = 0.17 day −1 , C M = 0.8 g/cm 3 [8] and λ C = 1.6/day [4], and in the steady state of the control case (no anti-tumor drugs), C = 0.4 g/cm 3 , and T = 3 × 10 −3 g/cm 3 . Hence by the steady state of Eq. (6), so that η = (λ C /2 − d C )/T = 210 cm 3 /g · day. In using the steady state equation we ignored the effect of the advection term and the fact that the tumor grows. We compensate for this omission by increasing the growth rate of cancer cells, taking λ C = 1.92/day.