|Home | About | Journals | Submit | Contact Us | Français|
The B-raf gene is mutated in up to 66% of human malignant melanomas, and its protein product, BRAF kinase, is a key part of RAS-RAF-MEK-ERK (MAPK) pathway of cancer cell proliferation. BRAF-targeted therapy induces significant responses in the majority of patients, and the combination BRAF/MEK inhibitor enhances clinical efficacy, but the response to BRAF inhibitor and to BRAF/MEK inhibitor is short lived. On the other hand, treatment of melanoma with an immune checkpoint inhibitor, such as anti-PD-1, has lower response rate but the response is much more durable, lasting for years. For this reason, it was suggested that combination of BRAF/MEK and PD-1 inhibitors will significantly improve overall survival time.
This paper develops a mathematical model to address the question of the correlation between BRAF/MEK inhibitor and PD-1 inhibitor in melanoma therapy. The model includes dendritic and cancer cells, CD 4+ and CD 8+ T cells, MDSC cells, interleukins IL-12, IL-2, IL-6, IL-10 and TGF- β, PD-1 and PD-L1, and the two drugs: BRAF/MEK inhibitor (with concentration γ B) and PD-1 inhibitor (with concentration γ A). The model is represented by a system of partial differential equations, and is used to develop an efficacy map for the combined concentrations (γ B,γ A). It is shown that the two drugs are positively correlated if γ B and γ A are at low doses, that is, the growth of the tumor volume decreases if either γ B or γ A is increased. On the other hand, the two drugs are antagonistic at some high doses, that is, there are zones of (γ B,γ A) where an increase in one of the two drugs will increase the tumor volume growth, rather than decrease it.
It will be important to identify, by animal experiments or by early clinical trials, the zones of (γ B,γ A) where antagonism occurs, in order to avoid these zones in more advanced clinical trials.
PD-1 is an immunoinhibitory receptor predominantly expressed on activated T cells [1, 2]. Its ligand PD-L1 is upregulated on the same activated T cells, and is also expressed by myeloid-derived suppressor cells (MDSCs) [2–5] and in some human cancer cells, including melanoma, lung cancer, colon cancer, and leukemia [2, 3]. The complex PD-1-PD-L1 is known to inhibit T cell function . Immune checkpoints are regulatory pathways in the immune system that inhibit its active response against specific targets. In the case of cancer, the complex PD-1-PD-L1 functions as an immune checkpoint for anti-tumor T cells. There has been much progress in recent years in developing checkpoint inhibitors, primarily anti-PD-1 and anti-PD-L1 .
The B-raf gene is mutated in up to 66% of human malignant melanomas, and its protein product, BRAF kinase, is a key part of the RAS-RAF-MEK-ERK (MAPK) pathway of cancer cell proliferation . BRAF-targeted therapy induces significant response in the majority of patients but the response is short lived (about 6 months) [7–9]. Initial clinical trials indicate that concurrent inhibition of BRAF with MEK decreases MAPK-driven acquired resistance, resulting in enhanced clinical efficacy and decreased toxicity [10, 11]. This provides a rational for using combined BRAF/MEK inhibition instead of BRAF inhibition alone . Treatment of melanoma with immune checkpoint inhibitors has a lower response rate compared to treatment with BRAF/MEK inhibitors, but the response tends to be more durable, lasting for years [11–13]. It was therefore suggested that BRAF/MEK-targeted therapy may synergize with the PD-1 pathway blockade to enhance anti-tumor immunity [4, 11, 14, 15]. Meta-Analysis of randomized clinical trials show that compared with other treatments of advanced BRAF-mutated melanoma, combined BRAF/MEK and PD-1 inhibitions significantly improved overall survival time .
In this paper we develop a mathematical model to address the efficacy of the combination of BRAF/MEK inhibitor (BRAF/MEKi) and anti-PD-1 (e.g. nivolumab). The model includes several types of T cells, MDSCs, and dendritic cells, as well as signaling molecules involved in the crosstalks among these cells.
Melanoma-derived factors alter the maturation and activation of tissue-resident dendritic cells, thus favoring tumor immune escape . In BRAF mutant melanoma, BRAF inhibitor restores the compromised dendritic cells function, and, in particular, the production of IL-12 by dendritic cells . Although MEK inhibitor (e.g. trametinib), as single agent, negatively impacts DC function, when combined with BRAF inhibitor (e.g. vemurafenib or dabrafenib), the functionality of DCs is restored, as well as their production of IL-12 [18, 19].
Dendritic cell-derived IL-12 activates effector T cells (Th1 and CD 8+ T cells) [20, 21]. Th1 produces IL-2 which further promotes proliferation of effector T cells. CD 4+ T cells (Th1) can kill cancer cell directly, for example, through FAS- or TRAIL-dependent pathway [22–25], while CD 8+ T cell is more effective in killing cancer cells . Cancer cells suppress the functions of effector T cells by producing immunosuppressor cytokines TGF- β, IL-6, CCL2 and IL-10 . IL-10 inhibits the activation of Th1 and CD 8+ T cells . IL-6 and CCL2 recruit MDSCs into tumor [19, 28, 29], and MDSCs produce TGF- β and IL-10. IL-6 and CCL2 also recruit T regulatory T cells (Tregs) [15, 28, 29]. TGF- β is produced not only by cancer cells and MDSCs, but also by Tregs , and Tregs become activated by TGF- β [30, 31]. Tregs modulate Th1 and CD 8+ T cells , thus promoting tumor growth.
One of the checkpoints on T cells is the membrane protein PD-1. Its ligand PD-L1 is expressed on activated effector T cells, on MDSCs and on cancer cells [2–5]. The complex PD-1-PD-L1 inhibits the function of effector T cells , but enhances the activation of Tregs  and thus promoting cancer.
The above interactions between cancer cells and the immune cells are summarized in Fig. Fig.1.1. The mathematical model developed in the present paper is based on Fig. Fig.1,1, and it includes BRAF/MEK and PD-1 inhibitors. Simulations of the model show that at low doses the two drugs are positively correlated, in the sense that the tumor volume decreases as each of the drugs is increased. However, at high doses the two drugs may become antagonistic, that is, an increase in dose of one of the drugs may actually result in an increase in the tumor volume.
The mathematical model is based on the network shown in Fig. Fig.1.1. The list of variables is given in Table Table1.1. Since CCL2 and IL-6 are both produced by cancer cells and both recruit MDSCs and Tregs into tumor environment, we shall consider, for simplicity, only IL-6 in our model.
We assume that the total density of cells within the tumor remains constant in space and time:
We assume that the density of debris of dead cells from necrosis or apoptosis is constant. We also assume that the densities of immature dendritic cells and naive CD 4+ and CD 8+ T cells remain constant throughout the tumor tissue. Under the assumption (1), proliferation of cancer cells and immigration of immune cells into the tumor, give rise to internal pressure 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 unit 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 diffusing within the tumor.
By necrotic cancer cells (N C) we mean cancer cells undergoing the process of necrosis. Necrotic cancer cells release HMGB-1 (H) . We model the dynamics of N C and H by the following equations:
where λNCC is the rate at which cancer cells become necrotic, d N is the rate at which necrotic cells turn into debris, and λHNC 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) , and we assume that the process of necrosis is also fast. We may then approximate the two dynamical equations by the steady statenλNCCC − dNCNC = 0 and λHNCNC − dHH = 0, so that H is proportional to C.
Dendritic cells are activated by HMGB-1 [35, 36]. Hence, the activation rate of immature dendritic cells, with density D 0, is proportional to , or to , since H is proportional to C. Here, the Michaelis-Menten law is used to account for the limited rate of receptor recycling time which takes place in the process of DCs activation. Hence, the dynamics of DCs is given by
where δ D is the diffusion coefficient and d D is the death rate of DCs.
Naive CD 4+ T cells differentiate into Th1 cells (T 1) under IL-12 (I 12) environment [20, 21], while IL-10 and Tregs inhibit the differentiation of naive CD 4+ T cells into T 1 cells [27, 30]. The proliferation of activated T 1 cells is enhanced by IL-2. Both processes of activation and proliferation of T 1 are assumed to be inhibited by PD-1-PD-L1 (Q) by a factor . Hence T 1 satisfies the following equation:
where T 10 is the density of the naive CD 4+ T cells.
Inactive CD 8+ T cells are activated by IL-12 [20, 21], and this process is resisted by IL-10 and Tregs [27, 30]. IL-2 enhances the proliferation of activated CD 8+ T cells. Similarly to the equation for T 1, T 8 satisfies the following equation:
where T 80 is the density of the inactive CD 8+ T cells.
Naive CD 4+ T cells differentiate into Tregs (T r) under activation by Fox3+ transcription factor. The complex PD-1-PD-L1 enhances the expression of PTEN, which results in upregulation of Fox3+, and hence in increased production of Tregs . The production of T r is also enhanced by TGF- β (T β) [30, 31]. The activated Tregs are recruited into tumor by tumor-derived immunosuppressive cytokine IL-6 (and CCL2)[15, 28, 29]. Representing this chemoattraction by ·(χ T rI 6), we get the following equation for T r:
Tumor recruits macrophages and “educates” them to become tumor-associated-macrophages (TAMs), which behave like MDSCs [37, 38]. MDSCs are also chemotactically attracted to the tumor microenvironment by IL-6 (and CCL2) [15, 19, 28, 29, 39]. As in , the Eq. of MDSCs is taken to be the following form:
where M 0 is the source/influx of macrophages from the blood.
Cancer cells are killed by T 1 and T 8 cells. We assume a logistic growth with carrying capacity (C M) in order to account for competition for space among cancer cells. BRAF/MEK inhibitor (B), for example vemurafenib/dabrafenib, is used for treatment of advanced melanoma. Its mechanism of action involves selective inhibition of the mutated BRAF kinase that leads to reduced signaling through the aberrant RAS-RAF-MEK-ERK (MAPK) pathway. We assume that BRAF/MEK inhibitor suppresses the abnormal proliferation of tumor cells by a factor . Then, the equation for C takes the form:
where η 1 and η 8 are the killing rates of cancer cells by T 1 and T 8, and d C is the natural death rate of cancer cells.
The maturation and activation of dendritic cells is interrupted by melanoma cells, which means that the production rate coefficient λI12D is small. However, in BRAF mutant melanoma, BRAF inhibitor alone or in combination with MEK inhibitor, restores the compromised dendritic cells function, and in particular, the production of IL-12 by dendritic cells [18, 19], and the corresponding equation for I 12 then takes the form:
IL-2 is produced by activated CD 4+ T cells (T 1) . Hence,
IL-10 is produced by cancer cells and MDSCs . Hence it satisfies the following equation:
PD-1 is expressed on the surface of activated CD 4+ T cells, activated CD 8+ T cells and Tregs. We assume that the number of PD-1 per cell is the same for T 1 and T 8 cells, but is smaller, by a factor ε T, for T r cells. If we denote by ρ P the ratio between the mass of one PD-1 protein to the mass of one T cell, then
The coefficient ρ P is constant when no anti-PD-1 drug is administered. And in this case, to a change in T=T 1+T 8+ε T T r, given by , there corresponds a change of P, given by . For the same reason, ·(u P)=ρ P·(u T) and 2 P=ρ P2 T when no anti-PD-1 drug is injected. Hence, P satisfies the equation
When anti-PD-1 drug (A) is applied, PD-1 is depleted (or blocked) by A. In this case, the ratio 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 . Hence,
where μ PA is the depletion rate of PD-1 by anti-PD-1.
PD-L1 is expressed on the surface of activated CD 4+ T cells, activated CD 8+ T cells, MDSCs, and tumor cells. We assume that the number of PD-L1 per cell is the same for T 1, T 8 and M cells, and denote the ratio between the mass of one PD-L1 protein to the mass of one cell by ρ L. Then
where ε C depends on the specific type of tumor.
PD-L1 from T cells or cancer cells combines with PD-1 on the plasma membrane of T cells, thus forming a complex PD-1-PD-L1 (Q) on the T cells [2, 3]. Denoting the association and disassociation rates of Q by α PL and d Q, respectively, we can write
The half-life of Q is less then 1 second (i.e. 1.16×10−5 day) , so that d Q is very large. Hence we may approximate the dynamical equation for Q by the steady state equation, α PL P L=d Q Q, or
where σ=α PL/d Q.
We assume that anti-PD-1 is injected intradermally every three days for 30 days (as in mouse experiments ), providing a source Â(t) of anti-PD-1:
where γ A is the effective level of the drug; although the level of the drug varies between injections, for simplicity we take it to be constant. The drug A is depleted in the process of blocking PD-1. Hence,
We assume that the BRAF/MEK inhibitor is injected intradermally every days for 30 days, providing a source of BRAF/MKEi:
Assuming that BRAF/MEKi is absorbed by C at a rate , we get the following equation for B:
We assume that a part of the tumor consists of extracellular matrix, ECM (approximately, 0.4 g/ cm3), cancer cells (approximately, C=0.4 g/cm3) and MDSCs (approximately, M=0.2 g/cm3). We assume (in the section of parameter estimation) that the densities of the immune cells D, T 1, T 8 and T r are approximately 4×10−4, 2×10−3, 1×10−3 g/ cm3 and 5×10−4 g/ cm3, respectively, and, for consistency, take the constant in Eq. (1) to be 0.6039. We further assume that all cells have approximately the same volume and surface area, so that the diffusion coefficients of all the cells are the same. Adding Eqs. (2)-(7), we then 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, 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.
We assume that the free boundary r=R(t) moves with the velocity of cells, so that
We assume that the naive CD 4+ T cells and inactive CD 8+ T cells that migrated from the lymph nodes into the tumor microenvironment have constant densities and at the tumor boundary, and that T 1 and T 8 are activated by IL-12 upon entering the tumor. We then have the following flux conditions at the tumor boundary:
We impose a no-flux boundary condition for all the remaining variables:
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.
Later on we shall compare the simulations of the model with experimental results for mice, for 60 days. Accordingly, we take initial values whereby the average density of cancer cells has not yet increased to its steady state. Then, by Eq. (1), the total density of the immune cells is initially above its steady state. We take (in unit of g/ cm3):
Note that the initial conditions for the cells satisfy Eq. (1).
We assume that initially B=0 and A=0, and take the initial condition for I 12, I 2, T β, I 6, I 10 and P to be close to their steady state values, which are computed in the section on parameter estimation. One choice of initial conditions is given as follows (in unit of g/ cm3):
However, other choices of these initial conditions do not affect the simulations of the model after a few days.
The simulations of the model were performed by Matlab based on the moving mesh method for solving partial differential equations with free boundary  (see the section on computational method).
Figure Figure22 is a simulation of the model with no drugs (the control case) for the first 60 days. The average density or concentration of a species is computed as its total mass in the tumor divided by the tumor volume. The simulation shows consistency in the choice of the model parameters. Indeed, as can be quickly checked, the steady states of all the cytokines and cells are approximately equal to the half-saturation values that we assumed in estimating the parameters of the model. Furthermore, the volume of the tumor doubles approximately every 10 days, as was assumed in the choice of the parameter λ 0 (used in estimating some parameters of Eq. (7)). It is interesting to note that the initial increase in TGF- β more than compensates for the initial decrease in P and L, as evident by the initial increase in T r. This initial increase of T r results in initial decrease in the T 1 and T 8 cells. We also note that the initial increase in cancer cells results in an increase in the D cells.
Figure Figure33 shows the growth of the tumor radius during 60 days when drug is administered. With no drugs, the radius increases from 0.01 cm to 0.037 cm. Treatment with BRAF/MEK inhibitor alone decreased the radius growth more than anti-PD-1 alone, and the combined therapy did better than anti-PD-1 alone. These results agree with mouse experiments reported in .
We next consider combination therapy for a range of values of BRAF/MEK inhibitor and anti-PD-1. We define the efficacy of a combination therapy, with (γ B,γ A), by the formula:
where the tumor radius R 60=R 60(γ B,γ A) is computed at day 60; R 60(0,0) is the radius at day 60 in the control case (no drugs). The efficacy is a positive number, and its value lies between 0 and 1 (or between 0 and 100%). Figure Figure44 is the efficacy map of the combined therapy, with γ B in the range of 0−5×10−9 g/cm3·day and γ A in the range of 0−1.4×10−9 g/cm3·day. The color column shows the efficacy for any pair of (γ B,γ A); the maximum efficacy is 0.97 (97%).
As the number of cancer cells increases, the tumor radius increases. Hence, if T 1 and T 8 were monotone increasing functions of γ A(or of γ B), then we should see that R 60(γ B,γ A) is a decreasing function of γ A(or of γ B), and E(γ B,γ A) would then also be an increasing function of γ A(or of γ B). But Fig. Fig.44 shows that this is not generally the case; indeed there are small oscillations in “northeast” corner of the figure. This means that the functions T 1 and T 8 cannot be monotone increasing with respect to γ B for fixed γ A>0.5×10−9 g/cm3·day, and also cannot be monotone increasing in γ A for fixed γ B>1.5×10−9 g/cm3·day. Indeed, for example, Fig. Fig.5a5a shows that the average densities of T 1 and T 8 are decreasing functions of γ B, for fixed γ A=1.26×10−9 g/cm3·day; however, for smaller values of γ A, T 1 and T 8 may become monotone increasing, as seen, for example, in Fig. Fig.5b5b with γ A=0.14×10−9 g/cm3·day. Similarly, Fig. Fig.6a6a shows that, for fixed γ B=3×10−9 g/cm3·day, there is a γ A-interval where T 1 and T 8 are decreasing as γ A increases. The γ A-interval where T 1 and T 8 are decreasing may shrink as we take a smaller fixed γ B, as seen, for example, in Fig. Fig.6b6b with γ B=0.1×10−9 g/cm3·day.
A possible explanation for Fig. Fig.5a5a is based on the antagonistic pathway shown in Fig. Fig.7.7. When γ B increases, the population of cancer cells decreases, and then, by Eqs. (2)-(4) and (8), so does the signal to activate T cells by dendritic cells-derived IL-12 (since the number of activated dendritic cells decrease with decreased cancer cell density) and thus the densities of T 1 and T 8 decrease. As for Fig. Fig.6a,6a, when γ A begins to increase, T 1 and T 8 also begin to increase, which results in a decrease of cancer cells. Then, as explained in the case of Fig. Fig.5a,5a, this leads to a decrease in dendritic cells-derived IL-12 and, hence, the density of activated T 1 and T 8 cells will begin to decrease as γ A continues to increase for a while.
If we inject IL-12 directly into tumor (as an additional drug), the influence of dendritic cells-secreted IL-12 diminishes, and the antagonism between BRAF/MEKi and anti-PD-1 also diminishes and it disappears already at very small amount of injection, e.g., an injection of order of magnitude 10−14 gcm3·day.
We performed sensitivity analysis, with respect to the rumor radius R at day 60 in the control case, with respect to some of the production parameters of the system (2)-(16), namely, λ DC, λT1I12, λT8I12, λTrTβ, λTβC, λI6C, λI10C,and the parameters K TQ, η 1 and η 8 which play important role in the dynamics of C. Following the method of , we performed Latin hypercube sampling and generated 1000 samples to calculate the partial rank correlation coefficients (PRCC) and the p-values with respect to the tumor radius at day 60. In sampling all the parameters, we took the range of each from 1/2 to twice its values in Tables Tables22 and and3.3. The results are shown in Fig. Fig.88.
We see that the production/activation rates that promote effector T cells, namely, λ DC, λT1I12 and λT8I12, are negatively correlated to the tumor radius, while the production/activation rates of the effector T cell-suppressors, such as λTrTβ, λI10C, λTβC and λI6C, are positively correlated to the tumor radius. The killing rate of effector T cells, η 1 and η 8 are negatively correlated to the tumor radius, and the correlation with η 8 is higher than with η 1.
BRAF mutation occurs in up to 66% of human malignant melanomas and for this reason BRAF has been one of the primary targets in melanoma therapy. Treatment with BRAF inhibitors (such as vemurafenib or dabrafenib) encounters MAPK-driven resistance, but combining it with MEK inhibitor (e.g. trametinib) significantly reduces this resistance as well as toxicity. While the response to the combined BRAF/MEK inhibitor is significant, it is short lived. On the other hand, PD-1 antibody (nivolumab) has lower response rate but a far greater durability. It was therefore suggested that BRAF/MEK inhibitor should positively correlate with anti-PD-1.
In the present paper we developed a mathematical model to test this hypothesis, in silico, by computing the efficacy of the combined therapy. The model is represented by a system of partial differential equations within the tumor tissue. The model includes immune cells (Th1 and CD 8+ T cells, Tregs, MDSCs and dendritic cells), cytokines (IL-12, IL-2, IL-6, IL-10 and TGF- β), and PD-1, PD-L1 and the complex PD-1-PD-L1. We simulated the model with combination of drugs, BRAF/MEK inhibitor at the ‘level’ γ B and PD-1 antibody at the ‘level’ γ A, and computed the tumor radius R 60=R 60(γ A,γ B) at day 60, and the efficacy ; the efficacy is an expression that quantifies the reduction in tumor size compared to the control case (no drugs).
The efficacy map in Fig. Fig.44 shows that for low levels of γ B and γ A, the two drugs are positively correlated, in the sense that tumor volume decreases as each of the drugs is increased. However, in the ‘northeast’ corner of Fig. Fig.44 we see that for higher levels of γ B and γ A there are zones where the drugs are antagonistic in the sense that when γ B and γ A in these zones are increased, the efficacy actually decreases. The antagonism between the combined drugs can be explained by the pathway shown in Fig. Fig.7.7. An increase in the number of effector T cells (Th1 and CD 8+) results in decrease in cancer cells and necrotic cancer cells, hence in decreased signals to activate dendritic cells. This results in a decrease in IL-12 production by dendritic cells, and hence in a decrease in effector T cells.
The parameter λI12B may be viewed as the immune system response to BRAF/MEK inhibitor. When this parameter is increased, the antagonism in the combined therapy is reduced, but it does not completely disappear (not shown here).
The mathematical model presented in this paper has several limitations:
A general study of synergistic and antagonistic networks in drug combinations appeared in . Clinical records on combination therapy show that the number of drugs that are synergistic far exceeds the number of drugs that are antagonistic .
In our model, the combination (γ B,γ A) is antagonistic when the drugs are administered in high doses, but not in low doses. For this reason it will be important to identify more carefully the zones of antagonism, by animal experiments or by early clinical trials, in order to avoid those zones in more advanced clinical trials.
In an expression of the form where Y is activated by X, the half-saturation parameter K X is taken to be the approximate steady state concentration of species X.
By , we have the following relation for estimating the diffusion coefficients of a protein p:
where M V and δ V are respectively the molecular weight and diffusion coefficient of VEGF, M p is the molecular weight of p, and M V=24kDa  and δ V=8.64×10−2 cm2 day−1 . Since, MI2 = 17.6kDa , MI12 = 70kDa , MTβ = 25kDa , MI6 = 21kDa [55, 56], MI10 = 20.5kDa , M A=32kDa  and M B=489.93Da , we get δI2 = 9.58 × 10−2 cm2 day−1, δI12 = 6.05 × 10−2 cm2 day−1, δTβ = 8.52 × 10−2 cm2 day−1, δI6 = 9.03 × 10−2 cm2 day−1, δI10 = 9.11 × 10−2 cm2 day−1, δ A=7.85×10−2 cm2 day−1 and δ B=3.16×10−1 cm2 day−1.
The number of DCs in various organs (heart, kidney, pancreas and liver) in mouse varies from 1.1×106 cells/ cm3 to 6.6×106 cells/ cm3 . In the dermal tissue, the number of DCs is larger (600-1500 cells/ mm2) [61, 62], but we do not specify where the melanoma cancer is located; it may be at its initial dermal tissue or in another organ where it metastasized. Mature DCs are approximately 10 to 15 μm in diameter . Accordingly, we estimate the steady state of DCs to be K D=4×10−4 g/cm3. We assume that there are always immature dendritic cells, some coming from the blood as tumor infiltrating dendritic cells (TID) [20, 21, 64]. We also assume that the density of immature DCs to be smaller than the density of active DCs, and take . From the steady statenof Eq. (2), we get λ DC=2d D D/D 0=4/day, since d D=0.1/day . We take K C=0.4 g/cm3.
The number of lymphocytes is approximately twice the number of DCs . T cells are approximately 14 to 20 μm in diameter. Assuming that the number of Th1 cells is 1/4 the number of lymphocytes, we estimate steady state density of Th1 cells to be KT1 = 2 × 10−3 g/cm3. We assume that the density of naive CD 4+ T cells to be less than the density of Th1, and take T10 = ⅕KT = 4 × 10−4 g/cm3. As in , we choose KTTr to be half-saturation of T r, that is, KTTr = 5 × 10−4 g/cm3, and as in , we choose KTI10 to be half-saturation of I 10, namely, KTI10 = 2 × 10−7 g/cm3. We assume that in steady state, Q/K TQ=2 (the value of K TQ is derived in the estimates of Eqs. (13)-(15)). From the steady state of Eq. (3), we get
The CD4/CD8 ratio in the blood is 2:1. We assume a similar ratio in tissue, and take T80 = ½T10 = 2 × 10−4 g/cm3. We also take steady state of T 8 to be the half of steady state of T 1, i.e., KT8 = ½KT1 = 1 × 10−3 g/cm3. From the steady state of Eq. (4), we have
We assume that TGF- β activates Tregs more than PD-1-PD-L1 does, and take λTrTβ = 5λTrQ. From the steady state of Eq. (5), we get, , where T 10=1×10−3 g/cm3, Tr = KTr = 5 × 10−4 g/cm3 , and dTr = 0.2/day . Hence λTrQ = 0.083/day and λTrTβ = 0.415/day.
The density of tumor-associated macrophages in melanoma can be up to 30% of the tumor tissue density ; we take MDSC density to be 20% of the tumor tissue density, so that M=0.2 g/cm3 in steady state. From the steady state of Eq. (6), we get, ½λM(M0 − M) = dMM, where d M=0.015/day , λ M=20/19=1.05 , and M=K M=0.2 g/cm3. Hence, M 0=0.21 g/cm3.
We take d C=0.17 day−1 and C M=0.8 g/cm3 . In the control case (no anti-tumor drugs), the tumor grows according to
during the volume doubling time in the control case, we conclude from Eq. (23) that
where . We assume that without immune responses and BRAF/MEK inhibitor,
We further assume that with immune response and BRAF/MEK inhibitor, the density of cancer cell still grows,
We take λ 0=0.069/day, and assume that in steady state, C is approximately 0.4 g/cm3, so that from Eq. (25) we get ½λC − dC = 2λ0, or λ C=2(2λ 0+d C)=0.616/day. By comparing Eq. (24) to Eq. (25), we see that η 1 T 1+η 8 T 8=λ 0. Noting that T 8 cells kill cancer cells more effectively than T 1 cells, we take η 8=4η 1, so that (with T1 = KT1 = 2 × 10−3 g/cm3 and T8 = KT8 = 1 × 10−3 g/cm3) and η 8=46 cm3/g·day. From Eq. (26), we have . Since λ C=2(2λ 0+d C) and η 1 T 1+η 8 T 8=λ 0, we get , so that (with B=K B=6.69×10−10 g/cm3) .
The serum level of IL-12 in melanoma patients varies from 7.5×10−11−9.6×10−11 g/cm3 [71, 72]. We assume that the IL-12 level in tissue is higher, and take I12 = KI12 = 8 × 10−10 g/cm3. In the control case (no drugs), from the steady state of Eq. (8), we get λI12DD − dI12I12 = 0, where dI12 = 1.38/day  and D=K D=4×10−4 g/cm3. Hence, λI12D = 2.76 × 10−6/day. In the simulations we take λI12B = 1, but simulations do not change qualitatively with smaller or larger values of λI12B.
The half-life of TGF- β is about 2 min , that is, t 1/2=0.0014 day, so that dTβ = ln2/t1/2 = 499.07 day−1. The concentration of serum TGF- β in melanoma is 2.68×10−14 g/cm3 . We assume that the concentration of TGF- β in tissue is higher than in serum, and take T β=2.68×10−13 g/cm3. By , λTβTr = 5.57 × 10−9/day. According to [27, 42], melanoma cells secrete more TGF- β than MDSC, and we assume that λTβCC = 2λTβMM. Hence, from the steady state of Eq. (10) we have, λTβCC + λTβMM + λTβTrTr = dTβTβ, or 3λTβMM + λTβTrTr = dTβTβ. Thus λTβM = (dTβTβ − λTβTrTr)/(3M) = 2.18 × 10−10/day, and λTβC = 2λTβMM/C = 2.18 × 10−10/day.
The half-life of IL-6 is less than 6 hours , and we take it to be 4 hours, that is, t 1/2=0.17 day, so that dI6 = ln2/t1/2 = 4.16 day−1. The concentration of serum IL-6 in melanoma is 3.4×10−12 g/cm3 . We assume that the concentration of IL-6 in tissue is higher than in serum, and take I 6=3.4×10−11 g/cm3. From the steady state of Eq. (11), we get λI6C = dI6I6/C = 3.54 × 10−10/day.
The half-life of IL-10 ranges from 1.1 to 2.6 hours ; we take it to be 2 hours, that is, t 1/2=0.08 day, so that dI10 = 8.32 day−1. The concentration of serum IL-10 in melanoma is 8.75×10−12 g/cm3 . We assume that the concentration of IL-10 in tissue is higher than in serum, and take I 10=8.75×10−11 g/cm3. In melanoma, the tissue concentrations of IL-10 secreted by tumor cells and by macrophages are similar , and, accordingly, we assume that λI10CC = λI10MM in steady state. Hence, from the steady state of Eq. (12) we get, 2λI10CC − dI10I10 = 0, so that λI10C = dI10I10/2C = 9.10 × 10−10/day, and λI10M = λI10CC/M = 1.82 × 10−9/day.
In order to estimate the parameters K TQ (in Eqs. (3) and (4)) and K Q (in Eq. (5)), we need to determine the steady state concentrations of P and L in the control case (no drugs). To do that, we begin by estimating ρ P and ρ L.
By , the mass of one PD-1 is m P=8.3×10−8 pg= 8.3×10−20 g, and by  the mass of one PD-L1 is m L=5.8×10−8 pg= 5.8×10−20 g. We assume that the mass of one T cell is m T=10−9 g. By , there are 3000 PD-1 proteins and 9000 PD-L1 proteins on one T cell (T 1 or T 8). Since ρ P T is the density of PD-1 (without anti-PD-1 drug), we get , and .
In order to estimate steady state concentration of P, we assume that the average densities of T 1, T 8 and T r are approximately 2×10−3, 1×10−3 and 5×10−4 g/ cm3, respectively. PD-1 is expressed by Tregs at higher or lower level than in T 1 and T 8 cells depending on the type of the cancer ; we assume that ε T=0.8. Hence, in steady state,
The parameter ε C in Eq. (14) depends on the type of cancer. We take ε C=0.01 . MDSCs express PD-L1 at lower level than tumor cells , and accordingly, we assume that εMM = ¼εCC, so that ε M=ε C C/4M=ε C/2=0.005. Then, by Eq. (14), we get
In mice experiments [44, 86] different amounts of drugs were injected, and the amount of BRAF/MEK inhibitor was larger than the amount of anti-PD-1. It is difficult to compare the amounts injected into mice with the actual levels of the drugs which appear in Eqs. (16) and (17), since there is no information available on the PK/PD of the drugs. For the choice of γ A=0.3×10−9 g/cm3·day and γ B=0.5×10−9 g/cm3·day, we found that the simulations are in qualitative agreement with experiments reported in . We shall accordingly take γ A in the range of n.4×10−9 g/cm3·day and γ B in the range of 0−5×10−9 g/cm3·day.
By , the half-life of anti-PD-1 is 15 days, so that . We assume that 10% of A is used in blocking PD-1, while the remaining 90% degrades naturally. Hence, μ PA P A/10%=d A A/90%, so that
The half-life of BRAF inhibitor (dabrafenib) is 8 hours , and the half-life of MEK inhibitor (trametinib) is 33 h . In the combination of BRAF/MEKi, the proportion of MEKi is smaller than BRAFi , and accordingly we take the half-life of BRAF/MEKi to be 10 h, so that . We assume that 10% of B is absorbed by cancer cells, while the remaining 90% degrades naturally, so that . From Eq. (17), we get B≥γ B/d B, and we assume that
where d B=1.66/day. We take γ B to be order of magnitude 10−9 g/cm3·day in the simulations. Hence, B=K B=6.69×10−10 g/cm3 in steady state, and μ BC=2d B B/9C=6.17×10−10/day.
Eqs. (20): We assume that is larger than KT1 and take . Similarly, we also assume that is larger than KT8 and take .
We employ moving mesh method  to numerically solve the free boundary problem for the tumor proliferation model. To illustrate this method, we take Eq. (2) as example and rewrite it as the following form:
where F represents the term in the right hand side of Eq. (2). Let and denote numerical approximations of i-th grid point and , respectively, where τ is the size of time-step. The discretization of Eq. (27) is derived by the fully implicit finite difference scheme:
where , ,, and . The mesh moves by , where is solved by the velocity equation.
This work is supported by the Mathematical Biosciences Institute and the National Science Foundation (Grant DMS 0931642), and by the Renmin University of China and the National Natural Science Foundation of China (Grant No. 11501568), and the International Postdoctoral Exchange Fellowship Program 2016 by the Office of China Postdoctoral Council.
The dataset supporting the conclusions of this article is included within the article.
XL and AF developed and simulated the model, and wrote the final manuscript. Both authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Xiulan Lai, Email: nc.ude.cur@ialnaluix.
Avner Friedman, Email: ude.uso.htam@namdeirfa.