Prediction of anti-CD25 and 5-FU treatments efficacy for pancreatic cancer using a mathematical model

Background Pancreatic ductal adenocarcinoma (PDAC) is a highly lethal disease with rising incidence and with 5-years overall survival of less than 8%. PDAC creates an immune-suppressive tumor microenvironment to escape immune-mediated eradication. Regulatory T (Treg) cells and myeloid-derived suppressor cells (MDSC) are critical components of the immune-suppressive tumor microenvironment. Shifting from tumor escape or tolerance to elimination is the major challenge in the treatment of PDAC. Results In a mathematical model, we combine distinct treatment modalities for PDAC, including 5-FU chemotherapy and anti- CD25 immunotherapy to improve clinical outcome and therapeutic efficacy. To address and optimize 5-FU and anti- CD25 treatment (to suppress MDSCs and Tregs, respectively) schedule in-silico and simultaneously unravel the processes driving therapeutic responses, we designed an in vivo calibrated mathematical model of tumor-immune system (TIS) interactions. We designed a user-friendly graphical user interface (GUI) unit which is configurable for treatment timings to implement an in-silico clinical trial to test different timings of both 5-FU and anti- CD25 therapies. By optimizing combination regimens, we improved treatment efficacy. In-silico assessment of 5-FU and anti- CD25 combination therapy for PDAC significantly showed better treatment outcomes when compared to 5-FU and anti- CD25 therapies separately. Due to imprecise, missing, or incomplete experimental data, the kinetic parameters of the TIS model are uncertain that this can be captured by the fuzzy theorem. We have predicted the uncertainty band of cell/cytokines dynamics based on the parametric uncertainty, and we have shown the effect of the treatments on the displacement of the uncertainty band of the cells/cytokines. We performed global sensitivity analysis methods to identify the most influential kinetic parameters and simulate the effect of the perturbation on kinetic parameters on the dynamics of cells/cytokines. Conclusion Our findings outline a rational approach to therapy optimization with meaningful consequences for how we effectively design treatment schedules (timing) to maximize their success, and how we treat PDAC with combined 5-FU and anti- CD25 therapies. Our data revealed that a synergistic combinatorial regimen targeting the Tregs and MDSCs in both crisp and fuzzy settings of model parameters can lead to tumor eradication. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-021-08770-z.


Introduction
Pancreatic cancer is the seventh leading cause of mortality related to cancer around the world [1] and the fourth leading cause in the United States [2], and it is predicted to be the second leading cause by 2030 [3]. Pancreatic ductal adenocarcinoma (PDAC) constitutes 85% of histological diagnoses of pancreatic cancer, and this subtype mostly emerges from the exocrine glands of the neck and head of the pancreas [4]. PDAC has a very poor prognosis with the lowest surveillance among all cancers, a five-year overall relative survival rate of 8% [1]. This disappointing surveillance is due to a delay in diagnosis of this disease because of having no specific symptom; thus, patients are usually diagnosed with metastasis or an advanced unresectable mass [5]. Additionally, the only current curative therapy for this disease is surgery, but only 20% of patients have locally resectable mass when diagnosed [6,7]. So, various therapeutic modalities like immunotherapy, neoadjuvant chemotherapy, radiotherapy, and surgery are performed to improve surveillance and alleviate the patient's discomfort, but there is no curative approach for an advanced PDAC already [8]. Recent findings showed that modulating the effector cells in the tumor microenvironment by immunotherapy or chemotherapy would lead to impressive clinical and experimental therapeutic effects on PDAC [9][10][11][12]. Besides, some specific features of PDAC, like the existence of many immunosuppressive mediators in its microenvironment that are surrounded by a dense stroma cause the tumor to physically block the penetration of the drug; thus, this issue emphasizes the promising efficacy of immunotherapy specifically when combined with chemotherapy [13]. Till now, the experimental studies showed the beneficial therapeutic effect of Anti-CD25 immunotherapy targeting Treg cells [14][15][16] as well as 5-fluorouracil (5-FU) chemotherapy targeting MDSCs [12,17,18] on PDAC, but to our knowledge, there is no experimental study on the efficacy of Anti-CD25 and 5-FU combination therapy for PDAC. Also, several mathematical models have been conducted on the antitumor effect of Anti-CD25 and 5-FU. Shariatpanahi et al. designed an ODE model to simulate the effect of MDSC depletion by 5-FU on tumor-immune system dynamics and to evaluate the effect of replication of this treatment on tumor degradation. In their study, they designed a simulation framework to capture the dynamics of tumor cells, MDSC, CTL, and NK cells with and without 5-FU treatment. Their study using in silico assessment of 5-FU treatment proposed a testable hypothesis in vivo/in vitro environments [19]. In another study, Loizides et al. constructed an ODE model to capture the tumorigenesis process and tumor interactions with the immune system.
They also modeled 5-FU therapy using a two compartmental pharmacokinetic/pharmacodynamics model. Their Gompertz model simulated overall characteristics of the inherent variability of in vivo tumor growth rates and 5-FU killing effects that was observed in the de novo animal cancer model and predicted by mathematical modeling [20]. Montiel et al. developed a mathematical model to interrogate the effects of immunotherapy using dendritic cells (DCs) on tumor-immune system interactions. Their model consists of five delay differential equations that are calibrated by experimental data and are used to test different immunotherapy protocols. By in silico assessment of DC therapy, they suggest that changing the infusion time and using more doses of DCs causes more degradation of tumor cells [21]. Although various mathematical model analyses were conducted to determine the efficacy of immunotherapy or chemotherapy for pancreatic cancer separately [22][23][24][25], the effect of anti-CD25 immunotherapy in combination with 5-FU chemotherapy on pancreatic ductal adenocarcinoma has not assessed so far. Therefore, in the present study, we constructed a mathematical model based on ordinary differential equations (ODEs) to describe the dynamic interactions among dominant cells and their cytokines in the pancreatic tumor microenvironment during different phases of treatments and provided a quantitative prediction of anti-CD25 and 5-FU efficacy for pancreatic ductal adenocarcinoma.

Model development
Biological concept Figure 1 demonstrates a simplified biological concept of the immunosuppressive mechanisms and antitumor activities mediated by various effector cells in the tumor microenvironment as well as the influence of 5-FU chemotherapy and anti-CD25 immunotherapy on these interactions which are used for mathematical modeling. Besides, we considered the following biological assumptions in our model with regards to the previous studies [19,[26][27][28][29][30][31][32][33][34].
1. Tumor cells have the logistic growth in deprivation of the immune system [28]. 2. Tumor cells stimulate the activation of cytotoxic T lymphocytes [19]. 3. In the tumor microenvironment, NK cells also stimulate the activation of cytotoxic T lymphocytes [19]. 4. NK cells and activated cytotoxic T lymphocytes induce tumor cell regression [19,29]. 5. The activity of NK cells and cytotoxic T lymphocytes declines over time after encountering tumor cells [19,30]. 6. Tumor cells also recruit MDSCs and stimulate their proliferation in the tumor microenvironment [19,31]. 7. MDSCs suppress cytotoxic T lymphocyte activation mediated by NK cells in the tumor microenvironment [19,26]. 8. Treg cells that are increased in the tumor microenvironment suppress the proliferation of NK cells, cytotoxic T lymphocytes, and T helper cells [32]. 9. Low doses of 5-FU chemotherapy inhibit tumor progression through the deactivation of MDSCs [19,27]. 10. Anti-CD25 immunotherapy inhibits tumor progression through the depletion of Treg cells [34].

Methodology and experiment setting
We confirm that no animals were involved in the present study. In this study, we proposed a mathematical model that is calibrated by in vivo data of two studies. Wu et al. [35] and Pu et al. [23] presented in vivo data for change in tumor size during 5-FU chemotherapy and anti-CD25 immunotherapy, respectively. The experiments of Wu  understanding of complex dynamics and regulatory mechanisms in the tumor microenvironment.
General characteristics of the mathematical model used for anticipating the therapeutic effects of 5FU and anti-CD25 on pancreatic cancer The mathematical model of the tumor-immune system of this study is constructed using ODEs to describe the dynamics of interactions between pancreatic cancer cells and major innate and adaptive immune system components includes cytotoxic T lymphocytes (CTLs), natural killer (NK) cells, myeloid-derived suppressor cells (MDSCs), TCD4+ cells (T helper), T regulatory cells (Tregs) and cytokines e.g. interleukin-2 (IL-2), interferon-gamma (IFN-γ) and transforming growth factorbeta (TGF-β). The model considered anti-CD25 monotherapy as well as low-dose 5-fluorouracil (5-FU) chemotherapy to predict their efficacy for pancreatic cancer. The anti-CD25 and 5-FU treatments modeling techniques include the techniques that were used in the studies of Dan-Hua HE et al. [34] and SP Shariatpanahi et al. [19], respectively. The mathematical model is constructed using nine ODEs (eqs. 1-9) that each equation describes the change in population/concentration of cells/cytokines. All model variables including C, N, T, M, H, R, I, F, S, and antiCD25 are timedependent; hence, for example, for the sake of brevity, we have written C instead of C(t).
Pancreatic cancer cell (C): In eq. (1), dC dt describes the time derivative of cancer cells, and this equation represents the dynamics of pancreatic cancer cells. The first term on the right side of the equation describes the Gompertzian growth of cancer cells, where a 1 is the growth rate of cancer cells in absence of treatment and C max is the carrying capacity or the maximum size of the tumor in this model [19]. The second and third terms represent the NK-mediated and CTL-mediated killing of accessible tumor cells ( C ) to NKs and CTLs, respectively, and by the depth of access l 1 [19].
represents the surface layers of tumor cells in tumor mass [19]. The maximum rates of NK-mediated and CTL-mediated tumor cell killing are parameters b 1 and c 1 , respectively [19]. The terms ð 1þe 1 R 1þd 1 ðantiCD25Þ 2 R 3 Þ and (1 + f 1 S) in the denominator of the third term describe the inhibitory effects of Tregs and TGF-β on CTL-mediated cancer cell killing. The proportional parameter d 1 regulates the weight (intensity) of anti-CD25 monotherapy and parameters e 1 and f 1 reflect the measure of the inhibitory effect of Tregs and TGF-β on CTL cytotoxicity, respectively [34]. The denominator of term ð 1þe 1 R 1þd 1 ðantiCD25Þ 2 R 3 Þ represents the Treg depletion by anti-CD25 antibody clone that it improves CTL cytotoxicity by inhibition of Tregs. The powers 2 and 3 for antiCD25 treatment and Treg cells (R) were used to adjust model for better calibration and parametrization. We assumed that chemotherapy using low dose 5-FU and monotherapy using anti-CD25 directly impacts the pancreatic tumor volume that is reflected in the fourth and fifth terms of eq. (1), respectively. The parameter g 1 represents the apoptotic rate of tumor cells by 5-FU and the parameters h 1 and k 1 are the apoptotic rate of pancreatic cancer cells by anti-CD25 treatment. The function sign(antiCD25) is used to indicate whether or not the anti-CD25 monotherapy is involved; therefore, this function is one when anti-CD25 treatment is applied and is zeros otherwise [37,38]. Term antiCD25(τ) presents the administration pattern of anti-CD25 therapy [34]. Term R t 0 antiCD25ðτÞCðτÞdτ represents the total dose of antiCD25(τ), 0 ≤ τ ≤ t injected until time t. This integral term computes the total dose of anti-CD25 and include this term in a logarithmic function as a nonmonotonic function to simulate the decay of affected cancer cell by treatment. Actually, this term indicates that the population of cancer cells at time t is depended on total dose of injected antiCD25(τ) at time τ ≤ t. The use of integral terms and non-monotonic functions for simulation of kinetic-dynamics of biological phenomena is convenient, for example, in a study, integral terms were used to simulate the dynamics of cytokines [39]. Also, in other study, Islam et al. used sliding model approach to determine anticancer agent dosing and they used non-monotonic functions to simulate anticancer therapeutic effect on cancer cells [40].
Natural killer cell (N): Eq. (2) describes the dynamic of NKs (N). The first and second terms on the right side of the equation represent the constant production rate of NKs in bone marrow with parameter a 2 and normal death rate of NKs with an exponential rate b 2 , respectively [19]. The third and fourth terms describe the IL-2 mediated and IFN-γ mediated stimulation of NK cells, respectively, that are modeled in the simplified version of the model (eq. 14) by a Michaelis-Menten form with parameters c 2 , d 2 τ 1 α 1 , e 2 and f 2 [34]. The fifth term represents the inactivation of NK cells by interacting with accessible tumor cells with maximum inactivation rate g 2 [19]. The last term describes Treg-mediated NK cell killing (inactivation) in a granzyme-B-dependent fashion with a constant rate h 2 [34].
Cytotoxic T lymphocytes (T): Eq. (3) describes the dynamics of CTLs as the major component of the adaptive immune system. The first term on the right side of the equation stands for the death of CTLs with an exponential rate a 3 [19]. The second term describes the tumor recruitment of CTLs that has a Michaelis-Menten form with parameters b 3 and c 3 [19]. The third term shows the activation of CTLs as a result of interactions of NKs and accessible tumor cells with a constant rate d 3 [19]. Also, MDSCs negatively regulate this interaction to prevent CTL activation. The parameter p 3 is the normal population of MDSCs when there is no tumor and parameters m 3 and n 3 are effectiveness factors and proportional parameters related to MDSC-mediated inhibition of CTL activation [19]. The fourth term represents the inactivation rate of interacting CTLs with accessible tumor cells with a constant rate e 3 [19]. The fifth and sixth terms show the stimulatory effects of cytokines IL-2 and IFN-γ on CTL proliferation/activation that they are modeled by Michaelis-Menten form with parameters f • 3 , g 3 , h 3 , k 3 and with parameters f 3 , g 3 τ 1 α 1 , h 3 , k 3 (eq. 15) [34]. The last term represents Treg-mediated CTL inhibition with a constant rate l 3 [34].
Myeloid derives suppressor cell (M): Eq. (4) describes the dynamics of tumor-induced MDSCs. The first term on the right side of the equation shows the constant recruitment rate a 4 of bone marrow produced MDSCs to the spleen [19]. The second term represents the death rate of MDSCs with a constant rate (b 4 ) in normal conditions or during 5-FU treatment [19]. The third term describes the tumor-induced expansion of MDSCs with parameters c 4 and d 4 [19].
T helper cells (H): Eq. (5) represents the dynamics of T helper cells. The first term on the right side of eq. (5) describes the production rate (a 5 ) of T helper cells in the thymus [34]. The second term describes the exponential degradation rate (b 5 ) of T helper cells [34]. The third and fourth terms show the activation rate of T helper cells by IL-2 and IFN-γ, respectively, that are modeled by a Michaelis-Menten form with parameters c • 5 , d 5 , e 5 , f 5 and parameters c 5 , d 5 τ 1 α 1 , e 5 and f 5 (the third and fourth terms of eq. 17) [34]. The last term models the inactivation/degradation of T helper cells by Tregs with a constant rate g 5 [34].
Regulatory T cell (R): eq. (6) describes the dynamic of Tregs as an important component of the immune system for the induction of peripheral tolerance. The first term on the right side of the equation describes the constant production rate of Tregs from their origin in the thymus and peripheral (a 6 ) [34]. The second term shows the exponential death rate of Tregs (b 6 ) [34]. The third and fourth terms show that Tregs originate from both TCD8+ and T helper cells and the population of Tregs increases with the sum of these two terms (c 6 T + d 6 H) that is proportional to CTLs population with a constant rate c 6 and T helper population with a constant rate d 6 [34]. The fifth term describes the IL-2-mediated growth rate of Tregs that is modeled by a Michaelis-Menten form with parameters e • 6 , f 6 , and parameters e 6 and f 6 τ 1 α 1 (the fifth term in eq. 18) [34]. The sixth term describes that NKs degrade Treg cells with a constant rate g 6 and the last term models anti-CD25-mediated Treg depletion with a constant rate h 6 [34].
Interleukin-2 (I): Eq. (7) represents the dynamic of IL-2. The first term on the right side describes the IL-2 secretion from T helper with a constant rate α 1 , and the second term describes the degradation with a rate τ 1 [34].
Transforming growth factor-beta (S): Eq. (9) describes the dynamics of TGF-β which is produced by cancer cells with a constant rate λ 1 and is degraded with a rate τ 3 based on its half-life [34].
The changes in concentration of cytokines compared to the population of cells occur on a very short time scale; therefore, we can use quasi-steady-state approximations (eqs. [10][11][12] for cytokine concentrations to simplify our model.
By substituting eqs. 10-12 into eqs. 1-6, the simplified model, as described using eqs. 13-18, is achieved. Finally, the simplified model of the tumor-immune system includes 6 ODEs and 65 kinetic parameters. Some kinetic parameters are estimated based upon the experimental data of 5-FU chemotherapy and anti-CD25 monotherapy in a murine model of PDAC in the C57/ BL6 mouse using Panc02 cell line.
The goodness of fit assessment For model fitting, we used a Genetic algorithm (GA) method to estimate model parameters. The cost function of GA is a normalized root-mean-square error (NRMS E), defined as follow: T 1 , T 2 , …, T 5 were the tumor volumes measuring at time points of t 1 , t 2 , …, t 5 . Also,T 1 ;T 2 ; …;T 5 were the predicted tumor volumes by the mathematical model at the same time points. The formulation of NRMSE is de- . Actually our model predicts the populaion of tumor cells that we convert it into tumor volume by considering the volume of a single cancerous cell is 10 −6 mm 3 [41,42].

Uncertainty analysis
There are two types of uncertainty in biological network modeling including fuzzy uncertainty and random uncertainty [43]. Due to imprecise, missing, or incomplete experimental data, the kinetic parameters of computational models are uncertain, and assigning a fuzzy uncertain number instead of crisp value to kinetic parameters seems a better choice [36][37][38][39][40][41][42][43][44][45]. Also, the effect of random uncertainty on dynamics of model constituents and robustness of model against perturbations can be determined by sensitivity analysis. For this aim, we used two different methods including partial rank correlation coefficient (PRCC) [46] and elementary effect (EE) test [47]. In the results section, we assess the effect of the fuzzy uncertainty of parameters on the dynamics of tumor cells in the presence and absence of treatments.
Also, the results of model robustness and the relation between dynamics of cells/cytokines and kinetic parameters are computed by sensitivity analyses that are provided in the results section.

Results
Model calibration for the prediction of the dynamics of pancreatic tumor cells in control case, 5-FU and anti-CD25 therapies We used GA to estimate model parameters in no treatment case to predict the dynamics of tumor-immune system constituents in the control group. We recruited the experimental data of the control group of study [35] (the recorded Panc02 tumor size on days 5, 10, 15, 20, and 25 after tumor inoculation on day 0) for model fitting and the resulted NRMSE was 0.1107. To estimate model parameters (g 1 ,b 4 ) related to the inhibitory effect of 5-FU therapy on pancreatic tumor cells, we used the data of the 5-FU treatment group of study [35] (the recorded Panc02 tumor cell population on days 5, 10, 20, and 25 after tumor inoculation on day 0, with regarding 5-FU therapy is carried out on days 1, 2, 3 and 4 after tumor inoculation). Similarly, the model was fitted to in vivo data by NRMSE = 0.1085. Finally, we used the experimental data of the anti-CD25 treatment group of study [23] (the recorded Panc02 tumor cell population on days 7, 14, 21, 28, and 35 after tumor inoculation on day 0, with regarding that anti-CD25 therapy is carried out on days 3, 6, 10, 13, 17 and 20 after tumor inoculation) to estimate model parameters (d 1 , h 1 , k 1 , h 6 ) related to the effect of anti-CD25 therapy on pancreatic tumor cells. The model was fitted to in vivo data by NRMSE of 0.0928. The model parameters, their description, units, and references are provided in Table 1. As shown in Fig. 2, model fitting is carried out and predicted dynamics of tumor population by the parameterized model is fitted to tumor sizes measured in the control group, 5-FU therapy group, and anti-CD25 therapy group.
We demonstrated the cytotoxic effects of anti-CD25 and 5-FU therapies in a separate or in a combinatorial manner on the dynamics of pancreatic tumor cells. We considered the initial population of the tumor cells to be 6 × 10 5 , based on in vivo data of study [23], and the initial condition for NKs, CTLs, MDSCs, TCD4+ cells, and Treg cells to be 105840, 21 × 10 4 , 3 × 10 3 , 564480, 42336, respectively [34]. Also, the initial concentration of cytokines IL-2, IFN-γ, and TGF-β are 2.5107 × 10 −6 , 3.9348 × 10 −7 , 5.363 × 10 −7 , respectively, that are computed by eqs. 10-12. The dynamics of all cells/cytokines in different strategies including control, 5-FU, anti-CD25, and combination therapy are shown in different panels of Fig. 3. As shown in Fig. 3, injection of 5-FU on days 1, 2, 3, and 4 after tumor inoculation on day 0, and applying anti-CD25 therapy on days 3, 6, 10, 13, 17, and 20, affected dynamics of tumor cells and other factors of the model. To further investigate the effects of treatments on the dynamics of tumor cells, we computed the dynamics of inhibition percentage of instantaneous tumor size and average tumor size. Also, analysis of interaction among treatments was carried out to investigate the synergistic, additive, or antagonistic effects of treatments in a combinatorial manner. Finally, we optimized the timings of 5-FU and anti-CD25 injections by GA which results are provided in the next sections.

Analysis of interactions among treatments
In this section, we aim to analyze the interaction among treatments. To determine synergistic, additive, or antagonistic effects of the combinatorial regimen of anti-CD25 and 5-FU therapies, we used the Bliss combination index (CI) [37]. We considered the instantaneous tumor size and the average tumor size as outcomes and we computed the CI values by setting different administration timings for anti-CD25 and 5-FU therapies. The interpretation of this analysis was based on the curve distance from the reference line and also the direction of the curve based on the reference line. For example, CI < 1 represents two treatments have synergistic effect if combined, CI = 1 indicates the additive effect of two treatment combination, and CI > 1 infers the antagonistic effect. The green dashed line indicates the threshold below which the 5-FU treatment and anti-CD25 treatment have a synergistic effect [37].
We defined the efficacy of anti-CD25 treatment (ϕ a ) at time point t by E(ϕ a = 1, ϕ 5 = 0, t) and the efficacy of 5-FU treatment (ϕ 5 ) at time point t by E(ϕ a = 0, ϕ 5 = 1, t) and their combination efficacy at time point t by E(ϕ a = 1, ϕ 5 = 1, t) that the formula of efficacy is as follows: where C(ϕ a = i, ϕ 5 = j, t) represents the tumor cell population at day t in control case (i = 0, j = 0), in 5-FU treatment case (i = 0, j = 1), in anti-CD25 treatment case (i = 1, j = 0) and in combination therapy case (i = 1, j = 1). If the tumor population at day t, C(ϕ a = i, ϕ 5 = j, t), is smaller than the population in the control case, C(ϕ a = 0, ϕ 5 = 0, t), then the efficacy is a positive number and its value is between 0 and 100%. As depicted in the first and second panels of Figs. 4 and 5. A, B, and C, we found although 5-FU therapy could suppress the tumor progression, anti-CD25 therapy has more regression impact than 5-FU therapy on the pancreatic tumor cell dynamics. Also, in silico assessments revealed the synergistic combinatorial manner has the most killing effect on tumor cells. This finding is consistent with in vivo data set that is used for model calibration (Fig. 2). The third panel of Figs. 4 and 5. A,  we realized that the combined therapy has a more synergistic effect on the instantaneous tumor size rather than averaged tumor size.

Optimization of therapy
In this section, we aim to answer this question that what is the optimal combination and frequency of drug administration to present clinically significant conclusions? Dosing and timing of drug exposure determine the toxicity sustained in the tumor. Since the present mathematical model is configurable for 5-FU and anti-CD25 injection timing, we can optimize the timing of drug administration with GA to minimize the tumor burden.
Since it has been shown that accumulation of pro-tumor immune cells such as Tregs and MDSCs mediate tumor cell regrowth, 18 to 20 days after tumor injection [48], we apply 5-FU and anti-CD25 therapies before day 20 to suppress the MDSCs and Tregs, respectively. To prevent toxic side effects of drugs, we assume that a maximum of two doses of each drug is allowed in each injection. The infusion times of 5-FU and anti-CD25 are denoted by t 5F and t CD , respectively, that t 5 F ¼½t 5 F 1 ; t 5 F 2 ; …; t 5 F N1 , t CD ¼½t CD 1 ; t CD 2 ; …; t CD N2 and N1 (in the present study we assume 4), N2 (we assume 6) represent the preassigned number of times that 5-FU and anti-CD25 will be injected, respectively. GA estimate the optimal timings of 5-FU (t 5F * ) and anti-CD25 (t CD * ) by minimizing the following cost function to minimize the population of cancer cells in treatment groups: , and C comb (t i ) represent the population of cancer cells at time point t i , i = 1, 2, …, N in control group, 5-FU group, anti-CD25 group, and combination therapy group, respectively. Also, t N is the end time of simulation which is day 100. The quantities C 5F * (t 5F * ), C CD * (t CD * ), C comb * (t 5F * , t CD * ) represent the minimized tumor population by applying 5-FU, anti-CD25, and combination therapy in optimal times, respectively.
As shown in first and second panels of Figs. 6. A, 6.C and 6.E, optimization of 5-FU anti-CD25 therapies caused the tumor inhibition percentage in combination therapy regimen to reach its maximum value (100%). As shown in third panel of Figs. 6. A, B and C, optimization of timing of 5-FU and anti-CD25 treatments caused these treatments to have a strong synergistic effect.

Results of fuzzy analysis
Due to the parametric uncertainty in the model, it seems that the evaluation of treatments in the fuzzy setting is more appropriate than the crisp setting. In this section, we aim to predict uncertain dynamics of cells/cytokines of the TIS model affected by the uncertainty of parameters. We assign a triangular membership function to some kinetic parameters and investigate the effect of that parametric uncertainty on the dynamics of model constituents. The simulation of parametric uncertainty with a fuzzy theorem in an ODE model of this study is similar to that used in other models including stochastic Petri net and continuous Petri net [43,45]. We simulated the ODE model with fuzzy parameters and in the presence and absence of 5-FU and anti-CD25 therapies.
In silico assessment of model in fuzzy setting revealed that 5-FU therapy on days 1, 2, 3, and 4 and anti-CD25 therapy on days 3, 6, 10, 13, 17, and 20 and combinatorial manner caused the uncertainty band of cancer cells (panel A of Fig. 7. A), MDSCs (panel D of Fig. 7. A), and Tregs (panel F of Fig. 7. A), shift left, toward lower population (volume) of these cells and the uncertainty band of NK cells (panel B of Fig. 7. A), and T helper cells (panel E of Fig. 7. A), shift right, toward an upper population of these cells.
The results of increasing the frequency of 5-FU injections are depicted in Fig. 7. B. By increasing the frequency of 5-FU injections, the uncertainty band of cancer cells (panel A of Fig. 7. B), MDSCs (panel D of Fig. 7. B), and Tregs (panel F of Fig. 7. B) shift left  Fig. 7. D) to be shifted right toward the higher population of these cells.

Results of global sensitivity analysis
We carried out global sensitivity analysis (GSA) regarding the population/concentration of cells/cytokines at days 20, 50, and 100 in the no treatment case as well as all kinetic parameters of the system described in Table 1. By using the GSA method in study [46], we carried out Latin hypercube sampling (LHS) and produced 10,000 samples to compute the partial rank correlation coefficient (PRCC) and the p-values regarding the population/ concentration of cells/cytokines at days 20, 50 and 100 after tumor inoculation. In LHS, we took the range of parameters ½ to twice their values in Table 1. We performed the PRCC analysis five times and calculated the mean (Fig. 8a) and standard deviation (Fig. 8b) of the statistically significant correlation values (p-value< 0.05). The computed p-values (Fig. 8c) are the maximum values between five times replications. In Fig. 8a, b and  Fig. 8. A, we see that the population of tumor cells at day 20 after tumor injection is positively correlated to the pro-tumor parameters a 1 (tumor growth rate), C max (maximum tumor size), and anti-immune parameter b 2 (death rate of NK cells) while it is negatively correlated to the anti-tumor parameters b 1 (NK-mediated killing rate of tumor cells), l 1 (depth of access of immune cell to tumor cells) and a 2 (the constant source of NK cells). In this instance, the negative correlation between tumor population at day 20 and parameter l 1 represents that if parameter l 1 is increased, the depth of immune cells access to tumor cells will increase, resulting in the killing of more tumor cells and consequently a decrease in the tumor cell population. Also, we see that the pro-tumor parameters, namely, a 1 , C max , and b 2 are negatively correlated to the NK cell population at day 20 while anti-tumor parameters such as b 1 and a 2 are positively correlated to the NK cell population at day 20. The NK cell population at day 20 is negatively correlated to the l 1 (depth of access of immune cells to the tumor) and g 2 (inactivation rate of NKs by tumor cells). Since there is a strong correlation between parameter b 3 (Maximum tumor-mediated CTL recruitment rate) and the CTL population at day 20, and according to the second term of eq. (15) that  Morris GSA was used to identify which of 52 kinetic parameters of the model have a significant effect on cell dynamics. For 52 kinetic parameters (with 10% perturbations), EE test using Morris sampling strategy were taken into account by setting 10 levels in the sampling grid and 1000 trajectories to compute the mean μ * and standard deviation σ. The identified most influential kinetic parameters with respect to the mean μ * and interaction effect σ are depicted in Fig. 8d. The parameters with large μ * values indicate the linear or additive effects while the parameters with large σ values indicate interaction effects. The dashed line μ Ã ¼ 2σ sqrtðrÞ (r is the number of trajectories) which all parameters are below that, translates into 95% confidence interval. Morris analysis was performed by considering the population of tumor cells, NK cells, CTLs, MDSCs, T helper, Tregs, and cytokines IL-2, IFN-γ, and TGF-β (in no treatment case) at day 100 of simulation as a read-out.
Our data revealed that the parameters a 1 , Cmax, reflecting the tumor growth rate and carrying capacity of tumor cells, respectively, are the most influential parameters for the tumor cell output (panel A of Fig. 8. D). The parameter a 1 has the most linear effect while the parameter Cmax has the most interaction effect on the dynamics of tumor cells. As depicted in panel B of Fig.  8. D, the parameters a 1 , Cmax, b 1 (NK-mediated tumor cell killing rate), c 1 (CTL-mediated tumor cell killing rate), e 1 (Rate of the suppressive effect of Treg on CTLmediated tumor cell killing) were predicted to play an important role in NK cells dynamics. The parameters a 1 , Cmax, b 1 have both linear and interaction effects while the parameters c 1 , e 1 have interaction effects on the dynamics of NK cells. Also parameter a 1 has strong interaction and liner effects on CTL's dynamics (panel C of Fig. 8. D). the results of the EE test show that parameters a 1 , Cmax have linear effects while parameter b 1 has an interaction effect on the dynamics of MDSCs (panel D of Fig. 8. D). Also, parameters Cmax, b 1 have an interaction effect, and parameter a 1 has both linear and interaction effects on the dynamics of T helper cells (panel E of Fig. 8. D). As depicted in panel F of Fig. 8. D, the parameters Cmax and a 1 have interaction and both linear and interaction effects, respectively. As shown in panels F, G, and H of Fig. 8. D, the parameter a 1 has most both linear and interaction effects on dynamics of Tregs, IL-2, and IFN-γ while the parameter Cmax has both interaction and linear effects on dynamics of Tregs and IL-2, respectively. As shown in panel K of Fig. 8. D, the parameters a 1 , Cmax have linear effects and parameter b 1 has an interaction effect on the dynamics of TGF-β.

Discussion
Mathematical modeling of complex networks of tumorimmune system (TIS) interactions is constantly advancing. Mathematical oncology helps us to better understand cancer biology, to assess the efficacy and toxicity of different treatment planning, and to predict dynamics of cancer and immune system behaviors during treatment. Systematic analysis of TIS by sophisticated modeling approaches can be used to refine and optimize drug dosing and scheduling. Dedicated modeling of chemotherapy/immunotherapy drug regimens by pharmacokinetic/pharmacodynamics models in oncology should result in the improved clinical efficacy of therapy and decreased toxicity. Mathematical modeling provides a (See figure on previous page.) Fig. 6 Dynamics of all cells/cytokines in strategies: control, 5-FU, Anti-CD25, and combination therapy with an optimized treatment schedule and comparison treatment efficacies and analysis of interactions among treatments. The blue lines in Fig. 6. B, D, and F show predicted dynamics of cells/cytokines in the control group, and the red lines show with 5-FU treatment (on specified time points), the yellow lines show with anti-CD25 therapy (on specified time points) and purple lines depict the dynamics under combinatorial manner. In each subplot in Fig. 6. B, D, and F (except for the second and third figures in the last panel), the y-axis represents the number/concentration of cell population/cytokine, and the xaxis represents the time in days after tumor inoculation. The initial condition is the same as those given in Fig. 3. The second and third figures in the fourth panels of fig. 6. B, D, and F show the terms related to 5-FU and anti-CD25 therapies, respectively. The overall efficacy of 5-FU treatment, anti-CD25 treatment, and combination therapy was plotted as the percentage of tumor growth inhibition using instantaneous tumor size (first panel of figs. 6 fig. 6. A, C and E) The Bliss combination index (CI) for 5-FU and Anti-treatments combination using instantaneous tumor size (CI inst) and average tumor size (CI ave) as an outcome personalized medicine approach to design a better efficacy-toxicity balance of therapy.
PDAC creates a tumor microenvironment that enhances tumor progression and metastasis. Tumorinduced dysfunction of immune cells is a critical issue in this microenvironment. Tumor cells induce immune suppression by many different mechanisms, including accumulation of regulatory T cells (Treg) and myeloid-derived suppressor cells (MDSCs). Interactions of immune system constituents including pro/anti-tumor cells and cytokines with tumor cells create a complex network with unknown behaviors that can be predicted by mathematical modeling. Mathematical modeling of the mutual interactions of tumor cells with the immune system in the tumor microenvironment helps us to understand the various processes involved in tumor growth and metastasis. The advent of mathematical modeling can help us to address treatment scheduling while simultaneously helping to unravel the processes driving therapeutic responses. The understanding complex interactions of tumor-immune system agents will elucidate the mechanisms of action of chemotherapy and immunotherapy drugs and lead to modify treatment schedules. Theoretically, a combined anti-CD25 immunotherapy and 5-FU chemotherapy would elicit a greater immune response. However, regimen scheduling H) and TGF-β (D. K). Each kinetic parameter is specified by two Morris indices, μ *(horizontal axis) and σ (vertical axis), which describe the significance of the effects and the interaction effects, respectively of combination 5-FU and anti-CD25 therapies has yet to be established. Here, we calibrate a mathematical model of the tumor-immune system (TIS) to in vivo data from 5-FU and anti-CD25 therapies to simulate the complex interplay between the invading tumor cells and innate and adaptive immune system constituents. The model helps us to implement an in silico clinical trial to test combination 5-FU and anti-CD25 therapies and optimize combination regimens to improve treatment efficacy. For this aim, we designed a user-friendly graphical user interface (GUI) unit (Fig. 9) that is configurable for 5-FU and anti-CD25 treatment timing in both crisp/ fuzzy settings. The present GUI, as a rigorous simulation framework, help us to predict dynamics of TIS constituents (cell/ cytokine) in different schedules of 5-FU and anti-CD25 therapies or absence of treatment with (fuzzy setting) and without (crisp setting) regarding parametric uncertainty (MATLAB codes of GUI and the additional file are in supplementary file).
In a study, Peng et al. designed a mathematical model to evaluate the efficacy of different treatments for prostate cancer and measured the interactions between different treatments based on tumor size at a given time point [37]. In the present study, we evaluated the effect of 5-FU and anti-CD25 treatments for pancreatic cancer by an in-vivo parameterized mathematical model and calculated the dynamics of interactions between these treatments based on instantaneous and average tumor size over time. Using the present model, we can capture the dynamics of TIS constituents and investigate the effect of different treatment schedules on the dynamics of TIS agents. Also, Peng et al analyzed the effect of parameter perturbations on cell dynamics by local sensitivity analysis methods while in the present study, we performed PRCC and EE tests which are global sensitivity analysis methods. Also, in the present study, we used fuzzy theorem to assess the effect of fuzzy uncertainty of kinetic parameters on dynamics of cells/ cytokines in the presence and absence of treatment modalities. The present model can be used as a rigorous simulation framework to predict the dynamics of TIS interactions and identify different behaviors of TIS in response to Fig. 9 The Graphical user interface (GUI) of the model. The user-friendly GUI of TIS with regarding fuzzy/crisp kinetic parameters for in silico assessment of 5-FU and anti-CD25 therapies treatments. The kinetic parameters of the present model are estimated by GA and based on published data from 5-FU and anti-CD25 pancreatic cancer treatments. The model of the present study is formulated by ordinary differential equations that by deterministic rules simulate the complex behaviors and interactions of TIS agents. Although deterministic models cannot capture uncertain cell-to-cell interactions and stochastic behaviors of TIS agents and also intrinsic noise in signaling pathways that control different behaviors of TIS agents, they have been widely used to simulate complex interplay between the invading tumor cells, innate and adaptive immune system constituents to predict dynamics of TIS agents with or without different treatment modalities. In the present study, we simulated the effect of uncertainty of model kinetic parameters through fuzzy theorem and we assessed the effect of different 5-FU and anti-CD25 treatment schedules in both crisp and fuzzy settings of kinetic parameters. Fuzzification of model parameters can help us to capture fuzzy uncertainty of model parameters which is due to imprecise, missing, or incomplete data.
Although no experimental or clinical study was conducted on the combination therapy of 5-FU and anti-CD25 so far, we evaluate the crosstalk between the immune system components on which these therapies mediated their consequences in the tumor microenvironment. Recently, Siret C et al. unraveled the underlying interaction between immune system cells within the PDAC microenvironment with a focus on Treg cells and MDSCs. They demonstrated that in the PDAC microenvironment, MDSCs suppress CTL proliferation, induce CTL death and produce arginase 1 and reactive oxygen species (ROS) and simultaneously Treg cells inhibit the proliferation of T helper cells that all these consequences provided strong immunosuppression in the PDAC microenvironment. Also, they showed in-vivo depletion of MDSCs inhibits induction and recruitment of Treg cells in PDAC microenvironment and also ex-vivo co-culture assays of Treg cells and MDSCs revealed tumoral MDSCs induce the development and proliferation of Treg cells mediated by cell-to-cell crosstalk and conversely the presence of Treg cells leads to survival and increase of tumoral MDSCs [48]. Moreover, the positive direct interactions between MDSCs and Treg cells in the tumor microenvironment were shown in other cancers [49][50][51][52][53]. In melanoma, Treg cells induce tumoral MDSCs differentiation through the expression of B7H1 and also the expression of Indoleamine 2,3-dioxygenase (IDO) by tumor cells in a Treg cells dependent manner recruits and activate MDSCs in the tumor microenvironment [50][51][52][53][54]. Moreover, Re GL et al. surveyed the combination therapy of Cyclophosphamide, 5-FU, and IL-2 for solid tumors, and as we know IL-2 against anti-CD-25 therapy induces Treg cell proliferation.
They reported the response duration of this combination therapy for pancreatic cancer was over 18 months and during this period, intravenous IL-2 in compared to subcutaneous administration leads to more platelet decrease, less platelet/lymphocyte decrease, and less Treg cells increase. However, the total number of lymphocytes and Treg cells increased after therapy [55].

Conclusion
We integrate a mathematical model of 5-FU and anti-CD25 into a simulation framework to optimize their administration in combination therapy. Using this framework, we inferred a combination schedule for the treatment of PDAC that significantly improved treatment outcomes when compared to 5-FU and anti-CD25 separately and provided a standard combination regimen. Our findings outline a rational approach to therapy optimization with meaningful consequences for how we effectively design treatment schedules to maximize their success, and how we treat PDAC with combined 5-FU and anti-CD25 therapy. In silico assessment of the model reveals that the combination of 5-FU and anti-CD25 treatments has potentially improved therapeutic effects through preventing tumor-induced immune suppressive mechanisms within the PDAC microenvironment.

Acknowledgments
We thank Hamid Reza Mirzaei, Leila Jafarzadeh, and all our collagenous in departments of Immunology and Biomedical Engineering at Tehran University of Medical Sciences for assistance with this study.
Authors' contributions SS conceived the initial idea of the study, developed the design of the study, collected, analyzed, and visualized the data and carried out program coding and model simulations, contributed to the interpretation of the findings, drafted the manuscript. HD helped draft the manuscript and contributed to the interpretation of the findings. AHJ supervised this research, designed the study, and critically revised the manuscript. AF and SR found the experimental data and participate in the design of the study. JH participated in the design of the study and critically revised the manuscript. All authors have provided critiques and feedbacks during the manuscript preparation and have read and approved the final version of the manuscript and are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work can be appropriately investigated and resolved.