Mathematical Modeling of Interleukin-35 Promoting Tumor Growth and Angiogenesis

Interleukin-35 (IL-35), a cytokine from the Interleukin-12 cytokine family, has been considered as an anti-inflammatory cytokine which promotes tumor progression and tumor immune evasion. It has also been demonstrated that IL-35 is secreted by regulatory T cells. Recent mouse experiments have shown that IL-35 produced by cancer cells promotes tumor growth via enhancing myeloid cell accumulation and angiogenesis, and reducing the infiltration of activated CD8 T cells into tumor microenvironment. In the present paper we develop a mathematical model based on these experimental results. We include in the model an anti-IL-35 drug as treatment. The extended model (with drug) is used to design protocols of anti-IL-35 injections for treatment of cancer. We find that with a fixed total amount of drug, continuous injection has better efficacy than intermittent injections in reducing the tumor load while the treatment is ongoing. We also find that the percentage of tumor reduction under anti-IL-35 treatment improves when the production of IL-35 by cancer is increased.


Introduction
Interleukin-35 (IL-35) is a member of the IL-12 cytokine family. It is produced in human cancer tissues such as in melanoma, B cell lymphoma [1], lung cancer, colon cancer, esophageal carcinoma, hepatocellular carcinoma, cervical carcinoma, and colorectal cancer [2,3], and it plays important roles in tumor progression and tumor immune evasion [1]. Fox3 z regulatory T cells (T reg ) are common in tumor microenvironment [4,5], where they induce immune-suppression. They do so by producing various cytokines, including TGF-b, IL-10 [6], and IL-9 [7], thereby promoting tumor growth. It was also shown that T reg secrete IL-35 [8][9][10][11][12][13][14]. IL-35 functions through IL-35R on various cell types, and is a potent immune-suppressor. Indeed, T reg -derived IL-35 was shown to inhibit antitumor T cell response [15], whereas IL-35-deficient T reg have significantly reduced activity in vitro and in vivo [8]. Stable expression of EBI3, a gene that codes for IL-35 subunit, confers growth-promoting activity in lung cancer, whereas small interfering RNA silencing of EBI3 inhibits proliferation of lung cancer [16].
Recently Wang et al. [1] generated IL-35 producing plasmacytoma cancer cells and showed that the expression of IL-35 in tumor microenvironment increased the number of myeloid derived suppressor cells (MDSCs), and promoted tumor angiogenesis; furthermore, IL-35 inhibited the infiltration of cytotoxic T lymphocytes into the tumor microenvironment and rendered the cancer cells less susceptible to CTL destruction.
These experimental results suggest that blocking IL-35 may be an effective therapeutic approach to human cancer. To explore this possibility we develop in the present paper a mathematical model and then conduct in silica experiments to evaluate to what extend blocking IL-35 reduces tumor growth.
The model consists of a system of partial differential equations (PDEs) that involve interactions among cells (tumor cells, MDSCs, T cells, T reg s, endothelial cells) and cytokines (M-CSF, TGF-b, VEGF, IL-35). We first consider the situation which corresponds to the experiments in Wang et al. [1]. In these experiments two kinds of plasmacytoma cells were injected into wild type mice: tumor cells that have been transfected with IL-35 (J558-IL-35) so that tumor secretes high amount of IL-35 into the microenvironment, and ''normal'' plasmacytoma cells (J558-Ctrl) that secrete very small amount of IL-35. There is also a small amount of IL-35 produced by MDSC [17,18] as well as IL-35 produced by T reg [8][9][10][11][12][13][14]. We show that the model simulations agree with the experimental data in [1]. We also introduce, in this model, the effect of a drug which inhibits production of IL-35, and simulate various protocols for administering the drug. We find, that administering the drug frequently in small amounts yields better results than administering it infrequently in larger amounts. We also find that the percentage of tumor reduction under anti-IL-35 drug improves when the production of IL-35 by cancer is increased.

Mathematical model
The mathematical model is based on the network schematically shown in Figure 1. Cancer cells secrete M-CSF which attracts MDSCs; cancer cells and MDSCs secrete VEGF which triggers angiogenesis by attracting endothelial cells and enhancing their proliferation. The additional roles of MDSC are described in the caption of Figure 1. In particular, MDSC, inhibits the activation CD8 z T cells via IL-10 and a variety of other mechanisms.
As mentioned in the Introduction, Wang et al. [1] considered two kinds of tumor cells injected into mice: J558-IL-35 and J558-Ctrl. In the case of J558-IL-35, IL-35 is produced mostly by tumor cells, less by T reg , and little by MDSC. In the case of J558-Ctrl, cancer cells produce very small amount of IL-35 so that IL-35 mainly comes from T reg and MDSC. MDSC secretes TGF-b and IL-10 which promote T reg [19,20], and there is a positive feedback loop where the last activation is activated by TGF-b and IL-10.
We use the network described in Figure 1 to construct a system of partial differential equations. In order to simplify the computations we assume that the tumor and all the variables are radially symmetric. The variables of the model and their dimension are listed below. We proceed to write down the differential equation of each of the variables. Most of the parameters are taken from the literatures, as indicated; in Methods we explain how the remaining parameters were estimated.
Tumor cell (c). The density c(r,t) of tumor cells satisfies the following equation: MDSCs promote T reg s, but also secrete MCP-1 to attract macrophages into the tumor microenvironment. Macrophages secrete IL-12 to activate CD4 z T cells, and CD4 z T cells secrete IL-2 which activates CD8 z T cells. MDSCs also produce large amount of IL-10, which inhibits the chemotaxis and activation of CD4 z T cells. doi:10.1371/journal.pone.0110126.g001 Lc Lt~D where w 0 is the oxygen level in heathy tissue, and the levels of oxygen for necrotic, extremely hypoxic, and intermediate hypoxic states vary in the intervals ½0,w n , (w n ,w h and (w h ,w 0 , respectively. The first term on the right-hand side of Equation (1) represents the dispersion (or diffusion) of tumor cells with diffusion coefficient D c . The second term accounts for the tumor proliferation, which depends on the concentration of oxygen w(r,t) and tissue carrying capacity c Ã . The third and fourth terms represent the death of tumor cells by necrosis and apoptosis, respectively. The last term accounts for the killing of tumor cells by activated CD8 z T cells [21]. The parameters in Equation (1) are listed in Table 1.
M-CSF (q). The concentration of M-CSF is given by the equation: The first term on the right-hand side is the diffusion of M-CSF with coefficient D q . The second term represents the M-CSF secreted by tumor cells [19,22], and the last term is the decay of M-CSF. The parameters in Equation (2) are listed in Table 2.
The first and last terms on the right-hand side account for the source and death of MDSCs. MDSCs undergo dispersion as well as chemotaxis driven by M-CSF (the third and fourth terms) [23][24][25]. It was reported in [1], that MDSCs do not undergo chemotaxis by IL-35 in vitro experiments. However, it has been observed that differentiation of MDSCs from myeloid precursor cells is enhanced by IL-35, although the mechanism is currently unknown [1]. We assume that this mechanism results in the second term on the right-hand side of Equation (3). The fifth term accounts for the differentiation of MDSCs from myeloid cells promoted by M-CSF [26]. The parameters in Equation (3) are listed in Table 3.

Regulatory T cell (R). The equation for the density of regulatory T cells is given by
T reg is activated by TGF-b (the third term on the right-hand side) and by IL-10. IL-10 is secreted by MDSC [19,20] and, for simplicity, we do not introduce IL-10 explicitly, and represent the activation of T reg by IL-10 by the term d M M=(Mzs R ). The parameters in Equation (5) are listed in Table 5.
Activated CD8 z T cell (T). Cytotoxic T cells (CTL), or CD8 z T cells, satisfy the equation: MDSC secretes MCP-1 which exerts chemotactic force on macrophages [39,40], while macrophages secrete IL-12 which activates CD4 z T cells [41] and CD4 z T cells produce IL-2 [42,43] which activates CD8 z T cells. The activation of CD8 z T cells is inhibited by TGF-b [44][45][46]. For simplicity we combine all these process by attributing the chemotactic force or CD8 z T cells and activation source of CD8 z T cells to MDSC (the terms in square brackets in Equation (7)). The factor s M =(s M za 1 M) represents the fact that MDSC suppresses CD8 z T cells proliferation by amino acid metabolism. The parameters in Equation (7) are listed in Table 7. VEGF (h). The concentration of VEGF evolves according to the equation where l 5 (w)~l 5 w(w) and l 6 (w)~l 6 w(w) depend on the oxygen concentration w, as follows: and w Ã [ (w h ,w 0 ) is the threshold at which the hypoxic effect on VEGF production by tumor cells and MDSCs is maximal. The function w(w) is chosen such that tumor cells and MDSCs can secrete VEGF under mild hypoxic conditions. The second term on the right-hand side of Equation (8) represents the VEGF produced by tumor cells and enhanced by I 35 [1], and the third term accounts for VEGF produced by MDSCs and enhanced by M-CSF [47]; accordingly, the ratios k 1 =s h and k 2 =q 0 should be small. The parameters in Equation (8) are listed in Table 8.
Endothelial cell (EC) (e). The equation of the density of EC includes dispersion, chemotaxis by VEGF, and proliferation by VEGF: Le Lt~D Here e 1 is the maximal density of EC inside the tumor, and H( : ) is defined by The last term, taken from [22], reflects the fact that VEGF induces proliferation of EC when the concentration of VEGF is higher than the threshold h 1 . The parameters in Equation (9) are given in Table 9.
Oxygen (w). We model the concentration of oxygen by the equation: Oxygen is delivered by EC (the first term) and is taken up by CD8 z T cells (the third term), MDSCs (the fourth term), T reg s (the fifth term), and tumor cells (the last term). The parameters in Equation (10) are listed in Table 10.
We assume that the tumor is radially symmetric and is contained in a sphere 0ƒrƒL, where L~1:5 cm.
We next introduce the initial and boundary conditions for each of the variables.
Initial conditions. We assume that the tumor cells are concentrated initially near r~0, and take with a positive parameter E, 0vEƒ1, and scaling parameters c 0~7 :2|10 8 cell=cm 3 and L 0~0 :5 cm. Since M-CSF is secreted by tumor cells, we take the initial concentration of M-CSF to be similar to the density of tumor cells, where the constant a q =m q comes from the steady state equation for q.
Since tumor cells are concentrated at the center r~0, we assume that the MDSC is higher at the center and negligible near the boundary r~L, where the constant s 0 =m M comes from the steady state equation of Equation (3). We assume that initially there are no activated CD8 z T cells, and take The activation of T reg s and the productions of I 35 and VEGF are triggered by tumor cells and MDSCs; accordingly, we take

<
: and I 0 35~1 0 2 pg=cm 3 , and h 0~1 0 3 pg=cm 3 . Similarly, I b is produced by tumor cells and T reg s, so accordingly we take where I 0 b~2 :4|10 3 pg=cm 3 . Endothelial cells migrate into the tumor from the surrounding normal healthy tissue, so we take where e 0 is the density of endothelial cell in normal healthy tissue. Finally, since endothelial cells represent capillaries through which oxygen is delivered, we prescribe where w 0 is the oxygen concentration in normal healthy tissue. Boundary conditions. Since we assume radial symmetry, the first r-derivative of each variable vanishes at r~0. We assume no-flux condition at r~L for all the variables except for the oxygen and endothelial cells, and we take where m is the flux rate of EC from healthy normal tissue into the tumor microenvironment. Parameters nondimensionalization.. We nondimensionalizate the Equations (1)-(10) by the following scaling:   fl l 1 (ŵ w),l l 2 (ŵ w),l l 5 (ŵ w),l l 6 (ŵ w)gt fl 1 (w),l 2 (w),c 0 l 5 (w)=h 0 ,M 0 l 6 (w)=h 0 g, fl l 7 ,l l 8 ,l l 9 ,l l 10 ,l l 11 ,l l 12 g~tfe 0 l 7 =w 0 ,T 0 l 8 ,M 0 l 9 ,R 0 l 10 ,c 0 l 11 ,l 12 g, The dimensional and nondimensional values of all the parameters of Tables 1-10 are summarized in Tables 11 and 12. After dropping the symbol ''^'', the model equations in the nondimensional form are as follows:

Numerical simulation
In accordance with the experiments in Wang et al. [1], we consider two types of mice plasmacytoma J558 cells in wild type mice: (i) J558-Ctrl tumor cells that secrete a very small amount of I 35 .
(ii) J558-IL-35 tumor cells that secrete a large amount of I 35 .
We use matlab with dr~1=40 and dt~7=216000 in nondimensional variables (i.e., dr~1=80 cm and dt~7=72000 day in dimensional variables). Figure 2 displays the spatial distributions of tumor cell density in cases (i)-(ii) at different times. We note that, in Figure 2, as time goes on, tumor cells migrate toward the boundary r~1:5 cm, where oxygen is rich while tumor cell density is lower near the center r~0 cm, where oxygen is sparse. The migration speeds of these two cases (i)-(ii) are similar to each other, but tumor cells with larger I 35 production (i.e., J558-IL-35 case) have higher peak during migration. The results of Wang et al. [1] were reported 2 weeks after injection of tumor cells into mice. Hence, we compare our simulations at the end of the second week with the results in [1]. In Figure 3(C), the ratio for MDSC of J558-IL-35 to J558-Ctrl is 2, which is the same as Figure five A in [1]. In Figure 3(H), the ratio for VEGF of J558-IL-35 to J558-Ctrl is 17, which is the approximately same as Figure four D in [1]. Next, we compare the ratio for T reg /CD8 z T cells of J558-IL-35 to J558-Ctrl with the result in [1]. But, in [1], they only showed the percentages of CD8 z /CD45 z , of CD4 z /CD45 z , and of Foxp3 z /CD4 z . By combining these results (Figures seven B, seven D, and seven E in [1]), we find that this ratio (for T reg /CD8 z T cells) is 0:54.
n c c |{z} production by tumor z n R R |{z} production by Treg   arrival of MDSCs to the tumor microenvironment is somewhat delayed and therefore the number of CD8 z T cells in the control case is significantly less than in the J558-IL-35 case, while (for simplicity) our model does not include such a time delay. The subunits of IL-35, EBI3 and IL-12p35, are highly expressed in cancers such as lung cancer, colorectal cancer, and esophageal carcinoma [2,3]. Anti-IL-35 drug blocks the expression of IL-35 and could be an agent in treating these cancers [48]. To determine the effect of anti-IL-35 drug on cancer growth, we proceed to introduce it, as a drug, into our model. If we denote its concentration by f (r,t) then all we need to do is to modify Equation (4) by    We make the pharmacokinetic assumption that f (r,t) decreases in r from the outer boundary of the tumor (r~1:5 cm) towards the center of the tumor (r~0), and take where a~L 2 (~2:25 cm 2 ) and F~10. We shall compare several dosing schedules: (i) no dosing of anti-IL-35, i.e., f (r,t)~1, for all t and 0ƒrƒL; (ii) continuous dosing with anti-IL-35 at fixed level F for 2 months, f (r,t)~F | r 2 za L 2 za , for 0ƒrƒL and 0ƒtƒ2 months ; ð16Þ (iii) intermittent dosing for 2 months, at double level 2F , one week at a time with one week spacing between dosing, f (r,t)~2 F | r 2 za L 2 za , for 0ƒrƒL and t 2i ƒtvt 2iz1 , 0, for 0ƒrƒL and t 2iz1 ƒtvt 2(iz1) , for i~0, 1, 2, 3, where t 0~0 and the length of each interval ½t j ,t jz1 is one week.
We use matlab with dr~1=80 cm and dt~7=24000 day in dimensional variables.  Figure 4 also shows that the reduction rate by anti-IL-35 is larger when tumor cells secrete higher amount of IL-35 as in Lung cancer and colorectal cancer [2,3] than lower amount of IL-35 as in plasmacytoma [1]. Accordingly, as a 35 increases, the reduction in total tumor population becomes increasingly significant.

Sensitivity analysis
In this section we perform sensitivity analysis on the parameters (in dimensional form) including those that were only roughly estimated and those that play important role in the model. We list these parameters with their ranges, baselines, and units in Table 13. We use the method described in Marino et al. [49], using the Latin hypercube sampling to generated 500 samples with dr~1=40 cm and dt~7=12000 day. Since we focus on how anti-IL-35 drug inhibits tumor growth, we calculate the partial rank correlation coefficients (PRCC) and p-value, corresponding to the ratio C :~Ð  Table 14 lists the PRCC and their p-values. Figure 5 plots the PRCC of the parameters with p-values smaller than 0:01. A negative PRCC (i.e. negative correlation) with pvalue smaller than 0:01 means that increasing this parameter value will decrease the value of C and hence increase the (relative) efficacy of the drug. A positive PRCC with p-value smaller than 0:01 has the opposite meaning, that is, it will decrease the efficacy of the drug. In Table 14, only g c , e 1 , l 5 , s M , s b , a 35 , and b 35 have negative PRCC with p-value smaller than 0:01. The most significant negatively correlated parameter is g c . Larger l 5 increases the production of VEGF and larger a 35 increases the production of I 35 and both increase tumor load. The negative correlation of these parameters shows that the drug is more effective for tumor with higher rate of production of VEGF and IL-35. On the other hand, the negative correlation of g c shows that the efficacy of the drug improves when the CD8 z T cells are more affective in killing tumor cells. However, it is not true to conclude that, in general, the drug efficacy increases with larger tumor load, since larger g c and s b shrink the tumor load but yield better drug efficacy. Similar results hold for the parameters with positive PRCC. For example, larger l 1 and s 0 lead to higher tumor cell population while the tumor efficacy is decreased.

Discussion
IL-35 is the most anti-inflammatory cytokine within the IL-12 cytokine family. In this paper we addressed the questions to what extend IL-35 is involved in tumor microenvironment and how effective is anti-IL-35 drug in reducing tumor growth. It is well known that T reg s are presented in the tumor microenvironment We next extended the model to include anti-IL-35 as an anticancer drug. We compared the efficacy of the drug under two schedules: continuous versus intermittent injections of the same total amount of the drug. We found that continuous injection has better efficacy while the treatment is ongoing. Since it is well known that some cancers including lung and colorectal cancers most likely secrete large amounts of IL-35, we also investigated the efficacy of the drug for such cancers. We found that the percentage of tumor reduction under anti-IL-35 drug improves when the production of IL-35 by cancer is increased.
There are currently only few experimental results by which our model can be tested. In recent experiments by Nicholl et al. [50] it was demonstrated that IL-35 promotes pancreatic cancer cells proliferation while anti-IL-35 reduces this promotion. More specifically, in Figure three of Nicholl et al. [50] it is shown that IL-35 (50 ng=ml) increases, on the average, by 100% the proliferation of colonies of several pancreatic cancer cell lines, while in the presence of anti-IL-35 (200 ng=ml) this increase is reduced to 50%. These in vitro results are in qualitative agreement with our results in Figure three (at week 8). Another example is taken from colorectal cancer in patients. As reported in Zeng et al. [2]. Foxp3 z T reg increases linearly with IL-35, and this is in qualitative agreement with Figures 3D and 3E of our simulations. As more experimental and clinical data become available, we should be able to test our model in more quantitative way, so that the model can further be refined.
In this paper we focused on the role of IL-35, although T reg secrete besides IL-35 also other cytokines that promote tumor, such as IL-10 and IL-9 [7,[51][52][53][54]; these were not included directly in the present model, since we wanted to base the model on the recent experimental data by Wang et al. [1]. When data for other cytokines become available to the same precision as, for instance, in [1], our model could then be extended to include these cytokines, and to obtain a more comprehensive evaluation of anti-IL-35 efficacy in combination with other drugs. We assume that, due to the secretion of IL-35, the production of MDSC in the present model is larger than the production assumed in [56], so we have taken s 0 and a M to be larger than in [56].
Estimate D I 35 and m 35 in Equation (4) Since IL-35 belongs to the IL-12 family, we assume that its diffusion coefficient and its degradation rate are the same as for IL-12 [60][61][62][63]: In the in vivo experiments of Wang et al. [1] the initial number of cancer cells that were injected was 5|10 6 and we assume that they occupy a volume of 50 mm 3 , so that There is no data in [1] on the density of the tumor cells in day 15, but the tumor cells were observed to grow rapidly in the first 15 days. We assume that the average of the density of tumor cells in the first 15 days is very close to the maximal capacity 10 9 cell=cm 3 and take, in (23), c~10 9 cell=cm 3 for J558-IL-35 tumor cells. Recalling Equation (18), we get, with m 35~2 =day (Table 4) Table 14 [11,13,14,27], little by MDSCs, and very little by tumor cells. Hence, in the J558-Ctrl case, we take the production rate of I 35 by tumor cells to be a 35~1 0 {7 pg=cell=day.
The production rate of I 35 by T reg is estimated to be b 35~1 :67|10 {3 pg=cell=day [34] and we take the production rate of I 35 by MDSCs to be small enough, i.e., c 35~1 0 {4 pg=cell=day, so that the production of I 35 in the J558-IL-35 case satisfies:  (5) In [38], the cytokine signalling by TGF-b on T reg is modeled bỹ whered d b~3 3:27=day which has dimension per day ands s b~1 which is nondimension. In our Equation (5), the dimension of d b is cell=cm 3 =day and the dimension of s b is pg=cm 3 . Correspondingly, we take d b~d d b |R 0~3 3:27=day|10 5 cell=cm 33 :327|10 6 cell=cm 3 =day, s b~s s b |I 0 b~1 |2:4|10 3 pg=cm 3~2 :4|10 3 pg=cm 3 , where I 0 b &2:4|10 3 pg=cm 3 [64]. MDSC also activates T reg population. We assume that the activation of T reg by MDSC is weaker than the activation of T reg by TGF-b, and hence take it to be d M~3 8 d b &1:25|10 6 cell=cm 3 =day: We also take s R~1 0 7 cell=cm 3 .
Estimate n c and n R in Equation (6) We assume as before that the initial tumor occupies a volume of 50 mm 3 and, accordingly, also T reg occupies the same volume. In [34], the production of I b by tumor cells and T reg s are Estimate s M , b 1 , b 2 , a 1 , a 2 , a 3 , c 5 in Equation (7) Since IL-35 enhances the population of MDSC, the concentration of IL-10, which we represent by a 1 M, is larger than the one in [56]. Hence, we chose s M to be larger than the corresponding value of s M in [56]. Moreover, since IL-35 promotes tumor growth, we expect a stronger immune response by T cells than in [56] and hence we take b 1 and b 2 larger than the corresponding value in [56]. The parameter c 5 is taken from [56]. Since the chemotaxis and activation of CD8 z T cells are indirect, we take a 2 and a 3 to be smaller than a 1 : a 1~2 pg=cell and a 2~a3~0 :01 pg=cell.