Combination therapy for mCRPC with immune checkpoint inhibitors, ADT and vaccine: A mathematical model

Metastatic castration resistant prostate cancer (mCRPC) is commonly treated by androgen deprivation therapy (ADT) in combination with chemotherapy. Immune therapy by checkpoint inhibitors, has become a powerful new tool in the treatment of melanoma and lung cancer, and it is currently being used in clinical trials in other cancers, including mCRPC. However, so far, clinical trials with PD-1 and CTLA-4 inhibitors have been disappointing. In the present paper we develop a mathematical model to assess the efficacy of any combination of ADT with cancer vaccine, PD-1 inhibitor, and CTLA-4 inhibitor. The model is represented by a system of partial differential equations (PDEs) for cells, cytokines and drugs whose density/concentration evolves in time within the tumor. Efficacy of treatment is determined by the reduction in tumor volume at the endpoint of treatment. In mice experiments with ADT and various combinations of PD-1 and CTLA-4 inhibitors, tumor volume at day 30 was always larger than the initial tumor. Our model, however, shows that we can decrease tumor volume with large enough dose; for example, with 10 fold increase in the dose of anti-PD-1, initial tumor volume will decrease by 60%. Although the treatment with ADT in combination with PD-1 inhibitor or CTLA-4 inhibitor has been disappointing in clinical trials, our simulations suggest that, disregarding negative effects, combinations of ADT with checkpoint inhibitors can be effective in reducing tumor volume if larger doses are used. This points to the need for determining the optimal combination and amounts of dose for individual patients.


Introduction
Prostate cancer is a major public health concern in the United States, with 248,000 new cases annually, and 34,000 deaths [1]. In metastatic prostate cancer, 5-year survival is 35% [2]. Androgen is a group of sex hormones that give men their 'male' characteristics. A major sex hormone is testosterone which is produced mainly in the testes. Prostate cells need androgen for their growth and function [3,4]. Androgen affects the immune system by increasing the proliferation of T regulatory cells (Tregs) through secretion of IL-10 [3,5,6]. Testoterone, upon entering prostate cells, is enzymatically converted into a more potent androgen, dihydrotestoterone (DHT), which binds to androgen receptor with more affinity [7].
Cancer vaccines serve to enlarge the pool of tumor-specific T cells from the naive repertoire, and also to activate tumor specific T cells which are dormant [17]. GM-CSF can activate dendritic cells, and is commonly used as a cancer vaccine [18][19][20].
PD-1 is an immunoinhibitory receptor predominantly expressed on activated T cells [21,22]. Its ligand PD-L1 is upregulated on the same activated T cells, and in some human cancer cells [21,23]. The compex PD-1-PD-L1 is known to inhibit T cells function [22]. Immune checkpoints are regulatory pathways in the immune system that inhibit its active response against specific targets. In case of cancer, the complex PD-1-PD-L1 functions as an immune checkpoint for anti-tumor T cells. CTLA-4 is another immunoinhibitory receptor expressed on activated T cells, the complex CTLA-4-B7 acts as a checkpoint inhibitor for anti-tumor T cells [24,25]. There has been much progress in recent years in developing checkpoint inhibitors, primarily anti-PD-1 and anti-PD-L1 (e.g., Nivolumab) [26], and anti-CTLA-4 (e.g., Ipilimumab) [27,28].
The standard care of metastatic prostate cancer is androgen deprivation therapy (ADT), commonly referred to as medical castration. Under ADT, blood tests show that patients develop adaptive immunity [29], and the level of effective T cells (Th1 and CD8 + T cells) increases. Enzalutamide (ENZ) is anti-androgen drug (approved in 2018) that inhibits androgen binding to androgen receptor on prostate cells, and it also inhibits androgen receptor from entering into the nucleus [30]. Clinical trials show that ENZ has significantly longer progression-free and overall survival than 'standard care' of androgen suppression [31]. ENZ is administered orally, once daily, with tablets or capsules [32].
In this paper we consider metastatic castration resistant prostate cancer (mCRPC), that is, metastatic prostate cancer with androgen-independent cancer cells. Sipuleucel-T (Provenge) (Sip-T) is a cancer vaccine (approved in 2010) for treatment of men with symptomatic or minimally symptomatic mCRPC. The vaccine is made by drawing immune cells from patients and culturing them with combinant fusion protein containing prostatic acid phosphotase (PAP) and GM-CSF. It is administered intravenously to activate dendritic cells [33], which indirectly increases antigen-specific T cells [34,35].
Treatment of mCRPC with ADT and PD-1 inhibitor has been disappointing [42], conferring only modest benefits [43], even though PD-L1 is increased under ADT [44]. Treatment with ADT and CTLA-4 inhibitor was also disappointing, since it did not increase the overall survival time [42]. Challenges and rationales for immune checkpoint inhibitors in the treatment of mCRPC are discussed in [44].
Preclinical trials with androgen ablation (ADT) and cancer vaccine show increase in both CD8 + T cells and Tregs [45]. Such a combination therapy is most effective when vaccine is delivered after ADT [46,47].
Vaccine Sip-T activates dendritic cells, and hence indirectly activates T cells. When ligand B7 on the activated dendritic cells combines with CTLA-4 or effective T cells, it initiates a signaling cascade that blocks the activation and proliferation of the T cells. This suggests that a combination therapy with ADT, Sip-T and CTLA-4 or PD-1 inhibitors may be effective in treatment of mCRPC.
Ardiani et al. [48] treated prostate cancer in mice with a combination of ENZ and a vaccine that targets the Twist antigen (involved in the epithelial-to-mesenchymal transition and metastasis) and increases the functional Twist-specific CD8 + T cells. ENZ was found to be immune inert since no changes were seen in CD4 + T cell proliferation and Treg functional assays, and ENZ did not also diminish the Twist vaccine's ability to generate CD4 + and CD8 + Twist-specific T cells responses. However, the combination of ENZ with Twist vaccine resulted in significantly increased overall survival of the mice compared to treatments with Twist vaccine alone (27.5 weeks vs 10.3 weeks). This suggests that combination of ENZ and immunotherapy is a promising treatment strategy for mCRPC.
In other mice experiments, Shen et al. [49] found that combination of ADT with anti-PD-1 and/or anti-CTLA-4 significantly delayed the development of castration resistance, reduced tumor volume and prolonged survival of tumor-bearing mice in some cases. Immunotherapy alone did not improve survival, and was ineffective if not administered in the peri-castration period.
There are several mathematical models of combination therapy with checkpoint inhibitors, for either generic or specific cancers. Lai and Friedman [54] considered combination therapy for melanoma with BRAF and PD-1 inhibitors. They showed that the combination is effective, in terms of tumor volume reduction, in "small" doses, but not in "large" doses. In [55] they considered treatment of a generic tumor with cancer vaccine (GVAX) and anti-PD-1. The vaccine produces GM-CSF which promotes activation of anti-cancer T cells. They addressed the question of which dose amounts and proportions to inject in order to increase synergy and efficacy. In another paper [56] they considered combination of PD-1 inhibitor with oncolytic virus (OV); the virus infects only cancer cells and replicates in them. Since CD8 + T cells kill both infected and uninfected cancer cells, they may either promote or suppress the tumor. They showed that anti-PD-1 in dose γ P in combination with OV in dose γ O is anti-cancer for one set of pairs (γ P , γ O ), while in the complementary set the combination is pro-cancer.
In [57] they considered combination of PD-1 and VEGF inhibitors and addressed the question in which order to administer the drugs in cases where VEGF inhibitor is known to affect the perfusion of other drugs. They showed that non-overlapping schedule of injections of the two drugs is significantly more effective than simultaneous injections. In Lai et al. [57] they considered treatment of breast cancer with CTLA-4 inhibitor in combination with BET inhibitor. They noted that more effective combinations to reduce the tumor volume result in higher level of toxicity, as measured by overexpression of TNF-α.
Cancer resistance was considered in Lai et al. [58] and Siewe and Friedman [59]. In [58] it was shown that anti-TNF-α reduces cancer resistance to anti-PD-1, and it is more effective if injected after anti-PD-1 injection, rather than simultaneously. In [59] it was shown that initial resistance to anti-PD-1, which is quite common, can be overcome by combination with TGF-β inhibitor, but the efficacy of the combination depends on two specific biomarkers.
ENZ inhibits androgen (A). It also inhibits androgen receptor (AR) from entering into the nucleus, which we take, in the model, as inhibiting AR. For simplicity, we shall simplify these two different activities of ENZ by combining "androgen" with "androgen receptor", and referring to it as androgen/receptor (A/AR) or, briefly, as androgen A.
In the present paper we develop a mathematical model to explore the efficacy of different combination therapies. The model includes androgen-dependent prostate cancer cells (N) and androgen-independent (castration-resistant) cancer cells (M), dendritic cells (D), Th1 cells (T 1 ), CD8 + T cells (T 8 ), T regulatory cells (Tregs, or T r ), and cytokines IL-12 (I 12 ), IL-10 (I 10 ) and IL-2 (I 2 ); the model includes also checkpoints PD-1 and CTLA-4 and their ligands PD-L1 and B7, respectively, and drugs. The M cells are cancer cells that underwent changes (e.g., epigenetic) so that they are adapted to survive and proliferate with (or little) androgen; for simplicity we refer to them as mutated cancer cells.
Androgen blockade increases the death rate of N cells [60] and the mutation rate of N to M [61][62][63]. Dendritic cells (D) are activated by the high mobility group box 1 (HMGB-1) expressed on necrotic cancer cells [9,10,64]. The activated dendritic cells secrete pro-inflammatory cytokine I 12 which induces the differentiation of naive T cells into T 1 cells and T 8 cells [11,12,65,66], a process inhibited by I 10 [67] and T r cells [68]. I 10 is secreted by cancer cells [3,5,6] and by Tregs [69,70], and Tregs differentiate from naive T cells under activation by Fox3p + transcription factor, a process enhanced by I 10 [69,70]. T 1 cells secrete cytokine IL-2 (I 2 ) which enhances the proliferation of T 1 and T 8 cells.
PD-1 and PD-L1 are expressed on T cells, and PD-L1 is expressed also on cancer cells. The complex PD1/PD-L1 blocks the anti-cancer activity of T 1 and T 8 cells [71], but also increases the proliferation of T r by mediating a phenotype change from T 1 to T r [72,73]. CTLA-4 is expressed on T cells, and B7 is expressed dendritic cells. The complex CTLA-4/B7 blocks the anti-cancer activity of T 1 and T 8 cells [74], and at the same time it also increases the proliferation of T r [75]; we assume that this increase in T r is caused by a change from T 1 to T r phenotype, as in the case of PD-1/PD-L1.
In this paper, we develop for the first time a mathematical model for cancer therapy that combines checkpoint inhibitors, vaccine and chemical castration. The mathematical model is represented by a system of partial differential equations based on Fig 1, which is a network describing the interactions among the cells, cytokines and checkpoints. The list of variables used in the model is given in Table 1. The list includes the following drugs: anti-PD-1, A 1 (nivolumab); anti-CTLA-4, A 4 (ipilimumab); ENZ (E) and Sip-T (S).

Mathematical model
The mathematical model is based on the network shown in Fig 1, with variables listed in Table 1. The variables satisfy a system of partial differential equations in a domain O(t), the region occupied by the cancer cells, which varies with time t.
We assume that the combined densities of cells within the prostate tumor O(t) remains constant in space and time: for some constant θ > 0. We assume that the densities of immature dendritic cells and naive CD4 + and CD8 + T cells remain constant throughout the tumor tissue. Under Assumption (1), proliferation of cancer cells and immigration of immune cells into the tumor, give rise to internal pressure which results in cells movement with velocity, u; u depends on space and time and will be taken in units of cm/day. We assume that cytokines and anti-tumor drugs are diffusing within the tumor, and that also cells undergo diffusion (i.e., dispersion).
In what follows, we denote by Y X K X þX a quantity proportional to the rate of proliferation/ activation of species Y by species X, and by Y 1 1þX=K X the rate proportional to the inhibition of Y by X. If Y is activated by two species, X 1 and X 2 , then we separately add each of the

PLOS ONE
Therapy for mCRPC with ICI, ADT and vaccine activated terms, but if Y is inhibited by X 1 and X 2 , then its total inhibition is proportional to Y

Equation for androgen-dependent cancer cells (N)
We assume a logistic growth for the androgen-dependent cancer cells with carrying capacity K NM , to account for the competition for space and nutrients among cancer cells. Androgens are primary regulators of prostate cancer cell growth and proliferation [60]. We accordingly model cancer cell growth rate β as an increasing saturating function of A, taking where A 0 is a level that corresponds to physiologically normal androgen concentration [76]. The drug ENZ (E) inhibits androgen binding to androgen receptor [30]. We represent its effect by multiplying β(A) by a factor cancer cells N mutate into androgen-independent cells M, at a rate that increases with decreasing androgen level [61][62][63]; we take this mutation rate to be proportional to Additional mutation from N to M results from the blockade of androgen receptor by ENZ, which we take to be proportional to 1 À Cancer cells are killed by CD8 + T cells [77,78]. We write the equation for N in the following form: |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } Growth of N À l NM 1 þ A=K A N |ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl } |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } where δ N is the diffusion coefficient, m T 8 N is the killing rate of cancer cells by T 8 and μ N is the natural death rate of cancer cells.

Equation for mutated androgen-independent tumor cells (M)
The dynamics of mutated androgen-independent cancer cells is given by |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } where δ M is a diffusion coefficient, q is their growth rate as they are recruited from the mutation of N, m T 8 M is the killing rate of cancer cells by T 8 , and their death rate is independent of androgen. Note that independent proliferation of castration-resistant cancer cells are included in the term ql NM =ð1 þ A=K A Þ when A is small.

Equation for dendritic cells (D)
The binding of extracellular high mobility box 1 (HMGB-1) to toll-like receptor 4 (TLR4) convert the immature dendritic cells, D 0 , into the activated tumor-associated dendritic cells [9,10,64] at a rate proportional to HMGB-1/(H 0 +HMGB-1), where H 0 is constant. Assuming that the concentration of HMGB-1 is proportional to the density of cancer cells, this activation rate is proportional to a linear combination of , where K N and K M are constants. The vaccine Sip-T (S) augments the activation of dendritic cells [33] by a factor λ DS S/ (K S + S), for some constants λ DS , K S . The dynamics of dendritic cells is given by |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } activation by HMGBÀ 1 where δ D is a diffusion coefficient, μ D is the death rate of dendritic cells, and the activation rates λ DN and λ DC are constants.

Equation for Th1 cells (T 1 )
Naive CD4 + T cells, T 10 , differentiate into Th1 cells under IL-12 inducement [11,12,65], and this process is inhibited by IL-10 [67] and Tregs [68]. The proliferation of activated CD4 + T cells is enhanced by IL-2 [79]. Activation and proliferation of T 1 cells are inhibited by the complex PD-1/PD-L1 (Q 1 ), represented by a factor 1 1þQ 1 =K TQ 1 [71], and by the complex CTLA-4/B7 (Q 2 ) as a factor 1 1þQ 2 =K TQ 2 [74]. The complex Q 1 also mediates phenotype change from T 1 cells to Tregs [72,73], at a rate l T r Q 1 T 1 þQ 1 , and Q 2 enhances naive Th cells to become Tregs [75], at a rate l T r Q 2 T 1 |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }

Equation for Tregs (T r )
Naive CD4 + T cells differentiate into Tregs under activation by Fox3p+ transcription factor, a process enhanced by IL-10 [69,70]. We have the following equation for T r : @T r @t þ r � ðuT r Þ À d T r 2 T r ¼ l T r I 10 T 10 I 10 K I 10 þ I 10 |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } where the second term in the right-hand side is the same as in Eq (5).

PLOS ONE
Therapy for mCRPC with ICI, ADT and vaccine

Equations for androgen (A)
Androgen is consumed by prostate cancer cells N at a rate proportional to β(A)E [3,80]. Hence, A satisfies the following equation where λ A is the constant production rate and μ A is the degradation rate.

Equations for PD-1 (P D ), PD-L1 (P L ) and PD-1/PD-L1 (Q 1 )
PD-1 is expressed on the membrane of activated CD4 + T cells, activated CD8 + T cells. We assume that the number of PD-1 proteins per cell is the same for T 1 , T r and T 8 cells. If we denote by r P D the ratio between the mass of the PD-1 proteins in one T cell to the mass of the cell, so that The coefficient r P D is constant when no anti-PD-1 drug is administered. In this case, to a change in T = T 1 + T 8 + T r , given by @T/@t, there corresponds a change in P D , given by r P D @T=@t. For the same reason, r � ðuP D Þ ¼ r P D r � ðuTÞ and r 2 P D ¼ r P D r 2 T when no anti-PD-1 drug is injected. Hence, P D satisfies the equation: Recalling Eqs (5)-(7) for T 1 , T 8 and T r , we get When anti-PD-1 drug (A 1 ) 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 r ) may change. In order to include in the model both cases of with and without anti-PD-1, we replace r P D in the above equation by P D /(T 1 + T 8 + T r ). where m P D A 1 is the depletion rate of PD-1 by anti-PD-1.

Hence,
We assume that the number of PD-L1 proteins in one T 1 cell is the same as in one T r cell and one T 8 cell, and denote by r 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 on cancer cells is r P L ε C . Hence, PD-L1 from T cells or cancer cells combines with PD-1 on the plasma membrane of T cells, forming a complex PD-1/PD-L1 (Q 1 ) on the T cells [21,23]. Denoting the association and disassociation rates of Q 1 by a P D P L and m Q 1 , respectively, we write Since the half-life of Q 1 is less than 1 second (i.e., 1.16 × 10 −5 day) [81], we may approximate the dynamical equation for Q 1 by the steady state equation a P D P L P D P L ¼ m Q 1 Q 1 , or where s 1 ¼ a P D P L =m Q 1 .

Equation for CTLA-4 (P A ), B7 (B 7 ) and CTLA-4/B7 (Q 2 ).
CTLA-4 is a receptor expressed on activated T 1 and T 8 cells [82] and the complex CTLA-4/B7 blocks the activities of these cells [74,82]. CTLA-4 is constitutively expressed on T r cells, but its activity is not blocked by the complex CTLA-4/B7 [83]. We assume that the number of CTLA-4 proteins per cell is the same for T 1 and T 8 cells, but different for T r cells, by a factor κ T . We denote by r P A the ratio between the mass of all CTLA-4 proteins in one T cell to the mass of this cell, so that The coefficient r P A is constant when no anti-CTLA-4 drug is administered. In this case, to a change in T = T 1 + T 8 + T r , given by @T/@t, there corresponds a change of P A , given by

PLOS ONE
Therapy for mCRPC with ICI, ADT and vaccine r P A @T=@t. Similar changes in P A arises from the terms of diffusion and advection, so that When anti-CTLA-4 drug (A 4 ) is applied, CTLA-4 is depleted at a rate proportional to A 4 , and, in this case, the ratio P A /(T 1 + T 8 + κ T T r ) may change. In order to include in the model both cases, with and without anti-CTLA-4, we replace r P A in the above equation by P A /(T 1 + T 8 + κ T T r ). Hence, where m P A A 4 is the depletion rate of CTLA-4 by anti-CTLA-4.
The ligand B7 is expressed on dendritic cells, so that

PLOS ONE
Therapy for mCRPC with ICI, ADT and vaccine CTLA-4 and B7 from the complex CTLA-4/B7 (Q 2 ) with association and disassociation rates a P A B 7 and m Q 2 , respectively: We assume that the half-life of Q 2 is very short [81,84], so that we may approximate the dynamics Q 2 by the steady state, a P A B 7 P A B 7 ¼ m Q 2 Q 2 , or Equations for anti-PD-1 (A 1 ) and anti-CTLA-4 (A 4 ). If a drug X with dose γ X and halflife t 1/2 is injected at time t 0 , we assume that its effect at time t (t > t 0 ) continues to be effective We shall compare our simulations with experimental results in [49], where PD-1 inhibitor and CTLA-4 inhibitor were injected at fixed dose in days 0, 3 and 6. The half-life of PD-1 inhibitor (nivolumab) is 26.7 days [85], and A 1 is depleted in the process of blocking PD-1, hence

Equation for ENZ (E).
In [49] the ADT drug was degarelix (G) and it was injected once every 30 days. The half-life of degarelix is 53 days [87], so its effective level at time t is g G e À ln 2 53 t , where γ G is the initial dose, with average 0.7γ G . In our model we let ENZ (E) take the role of degarelix. The drug ENZ has similar effect as degarelix, but is somewhat different in its mechanisms, and its half-life is 5.8 days [32]. In mice experiment [48] it was given in a way that maintained the level of daily dose (γ E ) constant. Since E is depleted in the process of inhibiting androgen, we have:

Equation for Sip-T (S).
In mCRPC clinical trials [88] Sip-T was administered with three infusions, two weeks apart. We approximate the effective level of the dose by a constant γ S . The drug is depleted in the process of activating dendritic cells, so that

Equation for cells velocity (u)
We assume that all cells have approximately the same diffusion coefficient. Adding Eqs (2)-(7) and using Eq (1), we get ½RHS of Eq: ð2:jÞ�: ð20Þ To simplify the computations, we assume that the tumor is spherical, and that all the densities and concentrations are radially symmetric, that is, functions of (r, t), 0 � r � R(t), where r = R(t) is the boundary of the tumor, and that u = u(r, t)e r , where e r is the unit radial vector.
Equation for the 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 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: it is tacitly assumed here that the receptors PD-1 and CTLA-4, and ligands PD-L1 and B7 become active only after the T cells are already inside the tumor.

Numerical simulations
All the computations were done using Python 3.5.4. The parameter values of the model equations are estimated in S1 File Section 1 and are listed in S1 File Tables 1 and 2. Parameter sensitivity analysis was performed in S1 File Section 2, and the techniques used for the simulations are in described in S1 File Section 3.

Model calibration
We simulated the model Eqs (2)-(21) with boundary conditions (23) and initial conditions, in units of g/cm 3 , I 2 ¼ 1:2 � 10 À 12 ; I 10 ¼ 3 � 10 À 10 ; I 12 ¼ 7 � 10 À 10 ; A ¼ 2:2 � 10 À 10 ; R ¼ 0:5 cm: We let the program run for 5 days (t = −5 to t = 0) before we began therapy. Fig 2 shows the profiles of the average densities/concentrations of the variables of the model, and of the tumor volume, with/without ADT. Without ENZ, the density of mutated cells (M) remains small, and tumor volume grows exponentially. With ENZ, given daily from t = 0 to t = 30, the tumor volume is first increasing, then decreasing during days 2-21, and finally it is again increasing. These changes in monotonicity can be explained by the fact that there is sharp decrease in androgen-dependent density (N) and slow increase in androgen-independent density (M) during an intermediate period, as seen in    Tables 1 and 2. "Pre-C" represents the level of the species at the time when the tumor volume attains its first maximum before decline due to ENZ, with γ E = 10 −7 g/ cm 3 �d, "ENZ-Effective" is the level of the species at the time when the tumor volume attains its lowest value under ENZ, and "C-Resistant" represents the level of the species at day 30 of treatment with ENZ, when androgen-resistance cells density (M) is at highest level.

PLOS ONE
Therapy for mCRPC with ICI, ADT and vaccine Degarelix is an androgen-receptor antagonist, which can be viewed as somewhat similar to ENZ in our model. Shen et al. [49] conducted mice experiments with treatment of prostate cancer using degarelix. The levels of T cells and DCs in Fig 3 are in qualitative agreement with Fig 3A in [49], and the levels of PD-1 and CTLA-4 are in qualitative agreement with Fig 4B in [49]. More precisely: As in [49], the level of DCs is decreasing through days 0 (Pre-C), 7 (ENZ Effective), 30 (C-Resistant); T 1 is decreasing-increasing; T r , P D and P A are increasing-decreasing. The profile of T 8 is increasing-decreasing while in [49] the profile of T 8 is constant; however, in [49] they also include the profile of NK cells which is increasing-decreasing while in our model we did not include NK, and, instead, let T 8 be the only cells who kill cancer cells. Hence, the T 8 in our model functions as T 8 + NK in the experimental results of [49]; and since in [49] NK is increasing-decreasing while T 8 is flat, there is a fit of our profile of T 8 with [49].
On the other hand, the concentrations of cytokine IL-2 in the microenvironment (outside the tumor) in Fig 5B of [49], cannot be compared with the concentrations in Fig 2, which is taken within the tumor, because of the large diffusion of cytokines.
In Fig 4A, we see that various combinations without ENZ do not reduce tumor volume significantly. This is in agreement with clinical trials referenced in [49]. In Fig 4B, we see that the combinations with ENZ increase efficacy, from 89.08% to 96.52%; the largest benefits are with combination of all the four drugs, A 1 +A 4 +ENZ+SipT. In particular, the combination with A 1 +A 4 increases efficacy from 89.08% to 94.41%; this moderate increase is in agreement with Fig 5A of [49], where degarelix was combined with α-PD-1 and α-CLTA-4 (ND).
We also note that the increase-decrease-increase profiles of the tumor volumes in

Therapy predictions
The parameter q is the ratio of growth rate of M to growth rate of N. According to [50], q is slightly smaller than 1 if the concentrations of DHT-activated androgen receptors and of testosterone-activated receptors are both the same for N and M. In our model we view q as a "personalized" parameter (a parameter in personalized, or precision, medicine), and let it vary in the interval 0.6 < q < 1.2.
We consider the case where the ENZ level is constant for 30 days, and it is either delivered as single agent or in combination with A 1 , A 4 or Sip-T by the same protocol as in Fig 4. Fig 5A shows the profile of tumor volume as function of q and time, 0 < t < 30.
We introduce two definitions to measure the benefit of treatment. Defining V drug (t) and V(t) as the tumor volume at time t under treatment and without treatment, respectively, the first definition, in terms of efficacy, is the following: The second definition is in terms of tumor volume reduction (TVR): Efficacy tells us how much we can reduce the tumor volume by treatment compared to no treatment; increased efficacy means improved treatment. But even very high efficacy does not inform whether the initial tumor was actually decreased. To get this information we look at TVR. With TVR, the larger it is the more the tumor volume was reduced compared to the initial volume, and TVR negative means that the treatment did not decrease the initial tumor volume. Clearly, a drug that increases efficacy also increases TVR. Fig 5A is a map showing the benefit of treatment with ENZ, as γ E varies in the range 0.5-1.8 × 10 −7 g/cm 3 �d, and q varies in the range 0.6-1.2. Fig 5B shows a similar map of benefits in terms of TVR. We see that, as γ E is increased and q is decreased, both efficacy and TVR increase. The range in benefits for efficacy is 70-94%, while for TVR it is −138% to 41%; for

PLOS ONE
Therapy for mCRPC with ICI, ADT and vaccine in Fig 4), and γ E varies in the range 0.6 − 1.8 × 10 −7 g/cm 3 �d; the dose amount in Fig 4 was 10 −7 g/cm 3 �d.
We see that efficacy of 95% corresponds, approximately, to 50% of TVR. Keeping γ E at the level of 10 −7 g/cm 3 �d, as in Fig 4, we can reduce tumor volume by nearly 60% if we increase g A 1 by 10 fold of its amount in Fig 4. The situation in Fig 6B with ENZ+A 4 is similar. We can decrease tumor volume by 50% if we increase g A 4 15 fold of its value of 2 × 10 −8 g/cm 3 �d in Fig 4.  Fig 6C shows that we can achieve 50% tumor volume reduction with ENZ+γ S if we use half the dose amount that was taken in Fig 4.

Conclusion
Androgen deprivation therapy (ADT) in combination with chemotherapy significantly increased overall survival time in patients with metastatic prostate cancer [89]. More recently, immune therapy by checkpoint inhibitors, has become a powerful new tool in the treatment of melanoma and lung cancer, and is currently used in clinical trials in other cancers, including metastatic castration resistant prostate cancer (mCRPC). Clinical trials, in increasing number, consider ADT in combination with cancer vaccine and immune checkpoint inhibitors (ICI), particularly for checkpoints CTLA-4 and PD-1 [39]. In the present paper, we developed a mathematical model to assess the efficacy of such combinations, as we vary the dose amounts and proportions of each agent in a combination. The model includes CD4 + and CD8 + T cells, dendritic cells, and cytokines by which these cells interact, as well as cancer cells (androgenindependent (M) and androgen-dependent (N)), and drugs. The densities/concentrations of these species are evolving within the tumor, and their evolution is described by a system of partial differential equations (PDEs); the tumor region is also evolving in time, and its volume growth is used to assess the effectiveness of treatments.
In previous work on metastatic castration resistant prostate cancer (mCRPC), Jain et al. [50] introduced several parameters as personalized parameters. In the present paper, we introduce one such parameter, q, which is the ratio of the growth rate of M cells to the growth rate of N cells.
Simulations of the model for 30 days are shown to be in qualitative agreement with experimental results for mice [49], where we used the same protocol of treatment, and took doses γ E = 10 −7 of ENZ (for ADT), g A 1 ¼ 4 � 10 À 9 (for anti-PD-1), g A 4 ¼ 2 � 10 À 8 (for anti-CTLA-4) in units of g/cm 3 �d, and q = 0.8. We then proceeded to evaluate (in Fig 6) the effectiveness of various combinations of γ E with g A 1 ; g A 4 and γ S (vaccine).
The experimental results in [49] show a tumor volume reduction of only 5-10%. On the other hand, the simulations in Fig 6 show that, in the mice model protocol of [49], we can achieve a much better tumor reduction by increasing the values of g E ; g A 1 ; and g A 4 . In particular, with fixed γ E and q as above, if g A 1 is increased 10 fold, the treatment with ðg E ; g A 1 Þ reduces tumor volume by nearly 60% (at day 30). Similarly, if g A 4 is increased 15 fold, the treatment with ðg E ; g A 4 Þ reduces tumor volume by 50%.
The model has several limitations: 1. We made a simplification by combining androgen with androgen receptor into one variable, which we just referred to it as androgen. This however does not affect the interactions associated with ADT by ENZ.

2.
The assumption (1) is another simplification, since it implies that non-cancerous prostate cells within the tumor have constant density, as if they were in homeostasis.
We did not discuss the question whether the PDE system of the model has a solution. This is indeed the case, and be proved by the same method as in [90].
Clinical trials of ADT and immune checkpoint inhibitors have been disappointing [42][43][44]. The simulations in Fig 6, based on mice experiments, suggest that combination of ADT with PD-1 and CTLA-4 inhibitors would have much more benefits if we increase significantly the dose of these checkpoint inhibitors.
We note however that in terms of clinical applications, PD-1 inhibition is associated with adverse events such as thyroid dysfunction and pneumonitis, CTLA-4 inhibition is closely associated with colitis and hypophysitis, and both drugs are associated with rash and hepatitis [91], and ENZ adverse events includes seizure and ischemic heart disease. This raises the question of determining the maximum dosages, in combinations of ICI and ENZ, that will reduce significantly these side effects. Another question that needs to be addressed in clinical setting is drug resistance, which is primary obstacle to successful cancer treatment. These issues are beyond the scope of the present work. However, the present paper can be used as a first step in addressing these clinical issues.
Supporting information S1 File. Parameters estimates, sensitivity analysis, numerical methods and tables of parameters. (PDF)