 Research
 Open Access
 Published:
Prediction of antiCD25 and 5FU treatments efficacy for pancreatic cancer using a mathematical model
BMC Cancer volume 21, Article number: 1226 (2021)
Abstract
Background
Pancreatic ductal adenocarcinoma (PDAC) is a highly lethal disease with rising incidence and with 5years overall survival of less than 8%. PDAC creates an immunesuppressive tumor microenvironment to escape immunemediated eradication. Regulatory T (Treg) cells and myeloidderived suppressor cells (MDSC) are critical components of the immunesuppressive 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 5FU chemotherapy and anti CD25 immunotherapy to improve clinical outcome and therapeutic efficacy. To address and optimize 5FU and anti CD25 treatment (to suppress MDSCs and Tregs, respectively) schedule insilico and simultaneously unravel the processes driving therapeutic responses, we designed an in vivo calibrated mathematical model of tumorimmune system (TIS) interactions. We designed a userfriendly graphical user interface (GUI) unit which is configurable for treatment timings to implement an insilico clinical trial to test different timings of both 5FU and anti CD25 therapies. By optimizing combination regimens, we improved treatment efficacy. Insilico assessment of 5FU and anti CD25 combination therapy for PDAC significantly showed better treatment outcomes when compared to 5FU 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 5FU 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.
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 fiveyear 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 AntiCD25 immunotherapy targeting Treg cells [14,15,16] as well as 5fluorouracil (5FU) chemotherapy targeting MDSCs [12, 17, 18] on PDAC, but to our knowledge, there is no experimental study on the efficacy of AntiCD25 and 5FU combination therapy for PDAC. Also, several mathematical models have been conducted on the antitumor effect of AntiCD25 and 5FU. Shariatpanahi et al. designed an ODE model to simulate the effect of MDSC depletion by 5FU on tumorimmune 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 5FU treatment. Their study using in silico assessment of 5FU 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 5FU 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 5FU 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 tumorimmune 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 antiCD25 immunotherapy in combination with 5FU 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 antiCD25 and 5FU 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 5FU chemotherapy and antiCD25 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 5FU chemotherapy inhibit tumor progression through the deactivation of MDSCs [19, 27].

10.
AntiCD25 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 5FU chemotherapy and antiCD25 immunotherapy, respectively. The experiments of Wu et al. study in the control group (no treatment) and 5FU therapy group are as follows. Control group: C57BL/6 mice were subcutaneously implanted with Panc02 cells (6 × 10^{5} per mouse) in the right flank and tumor size is recorded every 5 days up to 25 days beginning 5 days after tumor inoculation. 5FU therapy group: this group is similar to the control group, except that on the first to fourth days, at each day a single dose of 5FU (30 mg/kg) was injected. We used the trend of tumor size recorded over time (data records are extracted from figures by WebPlotDigitizer tool [36]) to estimate some parameters of the model which are describing the effect of 5FU therapy in pancreatic cancer. In the experiment of Pu et al., C57BL/6 mice were subcutaneously inoculated with Panc02 cells (2 × 10^{6} per mouse) and receiving CD25 antibody (50 μg per mouse) twice a week for 3 weeks beginning 3 days after tumor induction. The trend of tumor size measured in the Pu et al. study was used to estimate parameters reflecting the effect of antiCD25 therapy. Both studies of Wu et al. and Pu et al. were performed on a specific cell line (Panc02) and the same mice (C57BL/6 mice). We used the data of these experiments to evaluate the effect of the combinatorial manner of 5FU and antiCD25 therapies on pancreatic cancer. The parameterized mathematical model using control, 5FU, and antiCD25 data sets, reproduced the published in vivo data of both studies to describe tumorimmune cell interactions through 5FU and antiCD25 therapy. This mathematical model enables us to investigate the interaction effects of treatments (synergistic, additive, or antagonistic effect), optimize time schedules of treatments, and deepen our 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 antiCD25 on pancreatic cancer
The mathematical model of the tumorimmune 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, myeloidderived suppressor cells (MDSCs), TCD4+ cells (T helper), T regulatory cells (Tregs) and cytokines e.g. interleukin2 (IL2), interferongamma (IFNγ) and transforming growth factorbeta (TGFβ). The model considered antiCD25 monotherapy as well as lowdose 5fluorouracil (5FU) chemotherapy to predict their efficacy for pancreatic cancer. The antiCD25 and 5FU treatments modeling techniques include the techniques that were used in the studies of DanHua 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), \( \frac{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 NKmediated and CTLmediated killing of accessible tumor cells (\( \frac{\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}} \)) to NKs and CTLs, respectively, and by the depth of access l_{1} [19]. Term \( \frac{\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}} \) represents the surface layers of tumor cells in tumor mass [19]. The maximum rates of NKmediated and CTLmediated tumor cell killing are parameters b_{1} and c_{1}, respectively [19]. The terms \( \left(\frac{1+{e}_1R}{1+{d}_1{(antiCD25)}^2{R}^3}\right) \) and (1 + f_{1}S) in the denominator of the third term describe the inhibitory effects of Tregs and TGFβ on CTLmediated cancer cell killing. The proportional parameter d_{1} regulates the weight (intensity) of antiCD25 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 \( \left(\frac{1+{e}_1R}{1+{d}_1{(antiCD25)}^2{R}^3}\right) \) represents the Treg depletion by antiCD25 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 5FU and monotherapy using antiCD25 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 5FU and the parameters h_{1} and k_{1} are the apoptotic rate of pancreatic cancer cells by antiCD25 treatment. The function sign(antiCD25) is used to indicate whether or not the antiCD25 monotherapy is involved; therefore, this function is one when antiCD25 treatment is applied and is zeros otherwise [37, 38]. Term antiCD25(τ) presents the administration pattern of antiCD25 therapy [34]. Term \( {\int}_0^t antiCD25\left(\tau \right)C\left(\tau \right) d\tau \) represents the total dose of antiCD25(τ), 0 ≤ τ ≤ t injected until time t. This integral term computes the total dose of antiCD25 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 nonmonotonic functions for simulation of kineticdynamics 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 nonmonotonic 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 IL2 mediated and IFNγ mediated stimulation of NK cells, respectively, that are modeled in the simplified version of the model (eq. 14) by a MichaelisMenten form with parameters c_{2}, \( {d}_2\frac{\tau_1}{\alpha_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 Tregmediated NK cell killing (inactivation) in a granzymeBdependent 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 MichaelisMenten 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 MDSCmediated 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 IL2 and IFNγ on CTL proliferation/activation that they are modeled by MichaelisMenten form with parameters f^{●}_{3}, g_{3}, h_{3}, k_{3} and with parameters f_{3}, \( {g}_3\frac{\tau_1}{\alpha_1} \), h_{3}, k_{3} (eq. 15) [34]. The last term represents Tregmediated CTL inhibition with a constant rate l_{3} [34].
Myeloid derives suppressor cell (M): Eq. (4) describes the dynamics of tumorinduced 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 5FU treatment [19]. The third term describes the tumorinduced 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 IL2 and IFNγ, respectively, that are modeled by a MichaelisMenten form with parameters c^{●}_{5}, d_{5}, e_{5}, f_{5} and parameters c_{5}, \( {d}_5\frac{\tau_1}{\alpha_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 IL2mediated growth rate of Tregs that is modeled by a MichaelisMenten form with parameters e^{●}_{6}, f_{6}, and parameters e_{6} and \( {f}_6\frac{\tau_1}{\alpha_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 antiCD25mediated Treg depletion with a constant rate h_{6} [34].
Interleukin2 (I): Eq. (7) represents the dynamic of IL2. The first term on the right side describes the IL2 secretion from T helper with a constant rate α_{1}, and the second term describes the degradation with a rate τ_{1} [34].
Interferongamma (F): Eq. (8) shows the dynamic of IFNγ. According to terms 1–3 of eq. (8), IFNγ is produced by CD8 + T cells, NK cells and T helper cells with constant rates β_{1}, β_{2} and β_{3}, respectively [34]. The last term of the equation describes the natural death rate of IFNγ based on its halflife [34].
Transforming growth factorbeta (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 halflife [34].
The changes in concentration of cytokines compared to the population of cells occur on a very short time scale; therefore, we can use quasisteadystate approximations (eqs. 10–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 tumorimmune system includes 6 ODEs and 65 kinetic parameters. Some kinetic parameters are estimated based upon the experimental data of 5FU chemotherapy and antiCD25 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 rootmeansquare error (NRMSE), 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, \( {\hat{\mathrm{T}}}_1,{\hat{\mathrm{T}}}_2,\dots, {\hat{\ \mathrm{T}}}_5 \) were the predicted tumor volumes by the mathematical model at the same time points. The formulation of NRMSE is defined by \( NRMSE=\sqrt{\frac{\sum \limits_{i=1}^5{\left({T}_i{\hat{T}}_i\right)}^2}{{T_1}^2}} \). 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, 5FU and antiCD25 therapies
We used GA to estimate model parameters in no treatment case to predict the dynamics of tumorimmune 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 5FU therapy on pancreatic tumor cells, we used the data of the 5FU 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 5FU 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 antiCD25 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 antiCD25 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 antiCD25 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, 5FU therapy group, and antiCD25 therapy group.
We demonstrated the cytotoxic effects of antiCD25 and 5FU 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 IL2, 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, 5FU, antiCD25, and combination therapy are shown in different panels of Fig. 3. As shown in Fig. 3, injection of 5FU on days 1, 2, 3, and 4 after tumor inoculation on day 0, and applying antiCD25 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 5FU and antiCD25 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 antiCD25 and 5FU 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 antiCD25 and 5FU 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 5FU treatment and antiCD25 treatment have a synergistic effect [37].
We defined the efficacy of antiCD25 treatment (ϕ_{a}) at time point t by E(ϕ_{a} = 1, ϕ_{5} = 0, t) and the efficacy of 5FU 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 5FU treatment case (i = 0, j = 1), in antiCD25 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 5FU therapy could suppress the tumor progression, antiCD25 therapy has more regression impact than 5FU 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, shows a strong synergistic effect of the combination of 5FU and antiCD25 (in specified efficient design), while as depicted in the third panel of Figs. 5. B, C, and D, the ineffective administration of treatments (time settings) resulted in a poor synergy among treatments and consequently low efficacy in tumor degradation. Therefore, as shown in the third panel of Figs. 4 and 5. A, an appropriate schedule of 5FU and antiCD25 caused a strong synergistic effect of combined therapy on pancreatic tumor size and as time goes by, this effect becomes stronger (more distance from 1), while as shown in Figs. 5. B, C, and D, the specified timings of combination therapy led to an additive, an antagonistic, or poor synergy of treatments (lead to 1 or become > 1). Moreover, 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 5FU and antiCD25 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 protumor immune cells such as Tregs and MDSCs mediate tumor cell regrowth, 18 to 20 days after tumor injection [48], we apply 5FU and antiCD25 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 5FU and antiCD25 are denoted by t_{5F} and t_{CD}, respectively, that \( {\boldsymbol{t}}_{\mathbf{5}\boldsymbol{F}}=\left[{t}_{5{F}_1},{t}_{5{F}_2},\dots, {t}_{5{F}_{N1}}\right] \), \( {\boldsymbol{t}}_{\boldsymbol{CD}}=\left[{t}_{CD_1},{t}_{CD_2},\dots, {t}_{CD_{N2}}\right] \) and N1 (in the present study we assume 4), N2 (we assume 6) represent the preassigned number of times that 5FU and antiCD25 will be injected, respectively. GA estimate the optimal timings of 5FU (t_{5F}^{∗}) and antiCD25 (t_{CD}^{∗}) by minimizing the following cost function to minimize the population of cancer cells in treatment groups:
which C_{ctrl}(t_{i}), C_{5F}(t_{i}), C_{CD}(t_{i}), and C_{comb}(t_{i}) represent the population of cancer cells at time point t_{i}, i = 1, 2, …, N in control group, 5FU group, antiCD25 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 5FU, antiCD25, 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 5FU antiCD25 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 5FU and antiCD25 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 5FU and antiCD25 therapies. In silico assessment of model in fuzzy setting revealed that 5FU therapy on days 1, 2, 3, and 4 and antiCD25 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 5FU injections are depicted in Fig. 7. B. By increasing the frequency of 5FU 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 toward the lower population of these cells, respectively while the uncertainty band of NK cells (panel B of Fig. 7. B) and T helper cells (panel E of Fig. 7. B) shift right toward the upper population of these cells. As shown in Fig. 7. C, the antiCD25 treatment caused the uncertainty band of cancer cells (panel A of Fig. 7. C), the MDSCs (panel D of Fig. 7. C), and Tregs (panel F of Fig. 7. C) for no treatment case to be shifted left toward the lower population of these cells and the uncertainty band of NK cells to be shifted right toward the upper population of NK cells, but by increasing frequency of antiCD25 treatment, the uncertainty band of cancer cells, NK cells, and MDSCs cells does not move effectively. On the other hand, by increasing the frequency of antiCD25 treatments, the uncertainty band of T helper (panel F of Fig. 7. C) and Tregs (panel E of Fig. 7. C) shift right and left, respectively. Finally, as depicted in Fig. 7. D, increasing the regimens of both 5FU and antiCD25 in a combinatorial manner caused the uncertainty band of cancer cells (panel A of Fig. 7. D), MDSCs (panel D of Fig. 7. D), and Tregs (panel F of Fig. 7. D) to be shifted left toward the lower population of these cells while caused the uncertainty band of NK cells (panel B of Fig. 7. D) and T helper cells (panel E of 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 pvalues 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 (pvalue< 0.05). The computed pvalues (Fig. 8c) are the maximum values between five times replications. In Fig. 8a, b and c, only the significant correlation values (pvalue< 0.05) are reported. As shown in panel A of Fig. 8. A, we see that the population of tumor cells at day 20 after tumor injection is positively correlated to the protumor parameters a_{1} (tumor growth rate), C_{max} (maximum tumor size), and antiimmune parameter b_{2} (death rate of NK cells) while it is negatively correlated to the antitumor parameters b_{1} (NKmediated 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 protumor parameters, namely, a_{1}, C_{max}, and b_{2} are negatively correlated to the NK cell population at day 20 while antitumor 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 tumormediated CTL recruitment rate) and the CTL population at day 20, and according to the second term of eq. (15) that describe CTL recruitment, the CTL population is strongly affected by the tumor cell population (as the tumor cell population increases/decreases, the immune cell population increases/decreases). Since both tumor cells and CTLs have similar dynamics, their populations at day 20 are positively correlated with parameters a_{1} and C_{max}. According to panel A of Fig. 8. A, there exist a negative correlation between b_{1} (maximum killing rate of tumor cells by NKs), b_{2} (death rate of NKs), g_{2} (maximum inactivation rate of NKs by NKtumor interactions), and CTL population at day 20. The other correlation values are depicted in panels A, B, and C of Fig. 8. A, are similarly interpretable.
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 \( {\upmu}^{\ast }=\frac{2\sigma }{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 IL2, IFNγ, and TGFβ (in no treatment case) at day 100 of simulation as a readout.
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} (NKmediated tumor cell killing rate), c_{1} (CTLmediated 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, IL2, and IFNγ while the parameter Cmax has both interaction and linear effects on dynamics of Tregs and IL2, 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 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 myeloidderived suppressor cells (MDSCs). Interactions of immune system constituents including pro/antitumor 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 tumorimmune system agents will elucidate the mechanisms of action of chemotherapy and immunotherapy drugs and lead to modify treatment schedules. Theoretically, a combined antiCD25 immunotherapy and 5FU chemotherapy would elicit a greater immune response. However, regimen scheduling of combination 5FU and antiCD25 therapies has yet to be established. Here, we calibrate a mathematical model of the tumorimmune system (TIS) to in vivo data from 5FU and antiCD25 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 5FU and antiCD25 therapies and optimize combination regimens to improve treatment efficacy. For this aim, we designed a userfriendly graphical user interface (GUI) unit (Fig. 9) that is configurable for 5FU and antiCD25 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 5FU and antiCD25 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 antiCD25 treatments for pancreatic cancer by an invivo 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 treatments. The kinetic parameters of the present model are estimated by GA and based on published data from 5FU and antiCD25 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 celltocell 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 5FU and antiCD25 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 5FU and antiCD25 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 invivo depletion of MDSCs inhibits induction and recruitment of Treg cells in PDAC microenvironment and also exvivo coculture assays of Treg cells and MDSCs revealed tumoral MDSCs induce the development and proliferation of Treg cells mediated by celltocell 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,3dioxygenase (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, 5FU, and IL2 for solid tumors, and as we know IL2 against antiCD25 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 IL2 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 5FU and antiCD25 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 5FU and antiCD25 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 5FU and antiCD25 therapy. In silico assessment of the model reveals that the combination of 5FU and antiCD25 treatments has potentially improved therapeutic effects through preventing tumorinduced immune suppressive mechanisms within the PDAC microenvironment.
Availability of data and materials
All data generated or analysed during this study are included in this manuscript and its supplementary information files which includes MATLAB code (Please open MATLAB codes by MATLAB software or see the MATLAB codes in the supplementary file) and description of GUI.
Abbreviations
 PDAC:

Pancreatic ductal adenocarcinoma
 Treg:

Regulatory T cell
 MDSC:

Myeloidderived suppressor cell
 5FU:

5fluorouracil
 CTL:

Cytotoxic T lymphocyte
 NK:

Natural killer
 IL2:

Interleukin2
 IFNγ:

Interferonγ
 TGFβ:

Transforming growth factorβ
 TIS:

Tumorimmune system
 GUI:

Graphical user interface
 ODE:

Ordinary differential equation
 NRMSE:

Normalized rootmean square error
 PRCC:

Partial rank correlation coefficient
 EE:

Elementary effect
 CI:

Combination index
References
 1.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424. https://doi.org/10.3322/caac.21492.
 2.
Torre LA, Bray F, Siegel RL, Ferlay J, LortetTieulent J, Jemal A. Global cancer statistics, 2012. CA Cancer J Clin. 2015;65(2):87–108. https://doi.org/10.3322/caac.21262.
 3.
Rahib L, Smith BD, Aizenberg R, Rosenzweig AB, Fleshman JM, Matrisian LM. Projecting cancer incidence and deaths to 2030: the unexpected burden of thyroid, liver, and pancreas cancers in the United States. Cancer Res. 2014;74(11):2913–21. https://doi.org/10.1158/00085472.CAN140155.
 4.
Yang MX, Coates RF, Ambaye A, Gardner JA, Zubarick R, Gao Y, et al. Investigation of HNF1B as a diagnostic biomarker for pancreatic ductal adenocarcinoma. Biomark Res. 2018;6(1):25. https://doi.org/10.1186/s4036401801396.
 5.
Hackert T, Büchler MW. Pancreatic cancer: advances in treatment, results and limitations. Dig Dis. 2013;31(1):51–6. https://doi.org/10.1159/000347178.
 6.
Corbo V, Tortora G, Scarpa A. Molecular pathology of pancreatic cancer: from benchtobedside translation. Curr Drug Targets. 2012;13(6):744–52. https://doi.org/10.2174/138945012800564103.
 7.
GuillenPonce C, Blazquez J, Gonzalez I. deMadaria E, Montáns J, Carrato a. diagnosis and staging of pancreatic ductal adenocarcinoma. Clin Transl Oncol. 2017;19(10):1205–16. https://doi.org/10.1007/s1209401716817.
 8.
Mohammed S, George Van Buren II, Fisher WE. Pancreatic cancer: advances in treatment. World J Gastroenterol. 2014;20:9354.
 9.
Li X, Xu JX. A mathematical prognosis model for pancreatic cancer patients receiving immunotherapy. J Theor Biol. 2016;406:42–51. https://doi.org/10.1016/j.jtbi.2016.06.021.
 10.
Yang JQ, Wei T, Chen YW, Bai XL, Liang TB. Advances in immunotherapy of pancreatic ductal adenocarcinoma. Zhonghua Wai Ke Za Zhi. 2017;55(5):396–400. https://doi.org/10.3760/cma.j.issn.05295815.2017.05.018.
 11.
Corbett TH, Roberts BJ, Leopold WR, Peckham JC, Wilkoff LJ, Griswold DP, et al. Induction and chemotherapeutic response of two transplantable ductal adenocarcinomas of the pancreas in C57BL/6 mice. Cancer Res. 1984;44(2):717–26.
 12.
Herman JM, Swartz MJ, Hsu CC, Winter J, Pawlik TM, Sugar E, et al. Analysis of fluorouracilbased adjuvant chemotherapy and radiation after pancreaticoduodenectomy for ductal adenocarcinoma of the pancreas: results of a large, prospectively collected database at the Johns Hopkins Hospital. J Clin Oncol. 2008;26(21):3503–10. https://doi.org/10.1200/JCO.2007.15.8469.
 13.
Aroldi F, Zaniboni A. Immunotherapy for pancreatic cancer: present and future. Immunotherapy. 2017;9(7):607–16. https://doi.org/10.2217/imt20160142.
 14.
Tang Y, Xu X, Guo S, Zhang C, Tang Y, Tian Y, et al. An increased abundance of tumorinfiltrating regulatory T cells is correlated with the progression and prognosis of pancreatic ductal adenocarcinoma. PLoS One. 2014;9(3). https://doi.org/10.1371/journal.pone.0091551.
 15.
Ino Y, YamazakiItoh R, Shimada K, Iwasaki M, Kosuge T, Kanai Y, et al. Immune cell infiltration as an indicator of the immune microenvironment of pancreatic cancer. Br J Cancer. 2013;108(4):914–23. https://doi.org/10.1038/bjc.2013.32.
 16.
Hiraoka N, Onozato K, Kosuge T, Hirohashi S. Prevalence of FOXP3+ regulatory T cells increases during the progression of pancreatic ductal adenocarcinoma and its premalignant lesions. Clin Cancer Res. 2006;12(18):5423–34. https://doi.org/10.1158/10780432.CCR060369.
 17.
Wesolowski R, Markowitz J, Carson WE. Myeloid derived suppressor cells–a new therapeutic target in the treatment of cancer. J Immunother Cancer. 2013;1(1):10. https://doi.org/10.1186/20511426110.
 18.
Vincent J, Mignot G, Chalmin F, Ladoire S, Bruchard M, Chevriaux A, et al. 5Fluorouracil selectively kills Tumorassociated Myeloidderived suppressor cells resulting in enhanced t celldependent antitumor immunity. Cancer Res. 2010 [cited 2 Apr 2020]. https://doi.org/10.1158/00085472.CAN093690.
 19.
Shariatpanahi SP, Shariatpanahi SP, Madjidzadeh K, Hassan M, AbediValugerdi M. Mathematical modeling of tumorinduced immunosuppression by myeloidderived suppressor cells: implications for therapeutic targeting strategies. J Theor Biol. 2018;442:1–10. https://doi.org/10.1016/j.jtbi.2018.01.006.
 20.
Loizides C, Iacovides D, Hadjiandreou MM, Rizki G, Achilleos A, Strati K, et al. Modelbased tumor growth dynamics and therapy response in a mouse model of de novo carcinogenesis. PLoS One. 2015;10(12):e0143840. https://doi.org/10.1371/journal.pone.0143840.
 21.
CastilloMontiel E, ChimalEguia JC, Tello JI, PiñonZaráte G, HerreraEnríquez M, CastellRodríguez AE. Enhancing dendritic cell immunotherapy for melanoma using a simple mathematical model. Theor Biol Med Model. 2015;12(1):11. https://doi.org/10.1186/s1297601500070.
 22.
Louzoun Y, Xue C, Lesinski GB, Friedman A. A mathematical model for pancreatic cancer growth and treatments. J Theor Biol. 2014;351:74–82. https://doi.org/10.1016/j.jtbi.2014.02.028.
 23.
Pu N, Zhao G, Yin H, Li J, Nuerxiati A, Wang D, et al. CD25 and TGFβ blockade based on predictive integrated immune ratio inhibits tumor growth in pancreatic cancer. J Transl Med. 2018;16(1):1–13. https://doi.org/10.1186/s1296701816736.
 24.
Haeno H, Gonen M, Davis MB, Herman JM, IacobuzioDonahue CA, Michor F. Computational modeling of pancreatic cancer reveals kinetics of metastasis suggesting optimum treatment strategies. Cell. 2012;148(12):362–75. https://doi.org/10.1016/j.cell.2011.11.060.
 25.
Lee JJ, Huang J, England CG, McNally LR, Frieboes HB. Predictive modeling of in vivo response to gemcitabine in pancreatic cancer. PLoS Comput Biol. 2013;9(9):e1003231. https://doi.org/10.1371/journal.pcbi.1003231.
 26.
Condamine T, Gabrilovich DI. Molecular mechanisms regulating myeloidderived suppressor cell differentiation and function. Trends Immunol. 2011;32(1):19–25. https://doi.org/10.1016/j.it.2010.10.002.
 27.
AbediValugerdi M, Wolfsberger J, Pillai PR, Zheng W, Sadeghi B, Zhao Y, Hassan M. Suppressive effects of lowdose 5fluorouracil, busulfan or treosulfan on the expansion of circulatory neutrophils and myeloid derived immunosuppressor cells in tumorbearing mice. Int Immunopharmacol. 2016;40:41–9.
 28.
de Pillis LG, Radunskaya A. A mathematical model of immune response to tumor invasion. In: Computational Fluid and Solid Mechanics 2003: Elsevier; 2003. p. 1661–8.
 29.
Kawarada Y, Ganss R, Garbi N, Sacher T, Arnold B, Hämmerling GJ. NKand CD8+ T cellmediated eradication of established tumors by peritumoral injection of CpGcontaining oligodeoxynucleotides. J Immunol. 2001;167(9):5247–53. https://doi.org/10.4049/jimmunol.167.9.5247.
 30.
Kuznetsov VA. Basic models of tumorimmune system interactions identification, analysis and predictions. A survey of models for TumorImmune system dynamics. Springer; 1997:237–294, DOI: https://doi.org/10.1007/9780817681197_6.
 31.
Youn J, Collazo M, Shalova IN, Biswas SK, Gabrilovich DI. Characterization of the nature of granulocytic myeloidderived suppressor cells in tumorbearing mice. J Leukoc Biol. 2012;91(1):167–81. https://doi.org/10.1189/jlb.0311177.
 32.
Kyewski B, SuriPayer E, eds. CD4+ CD25+ regulatory T cells: origin, function and therapeutic potential. Vol. 293. Springer Science & Business Media; 2005.
 33.
Ghaffari A, Bahmaie B, Nazari M. A mixed radiotherapy and chemotherapy model for treatment of cancer with metastasis. Math Methods Appl Sci. 2016;39(15):4603–17. https://doi.org/10.1002/mma.3887.
 34.
He DH, Xu JX. A mathematical model of pancreatic cancer with two kinds of treatments. J Biol Syst. 2017;25(01):83–104. https://doi.org/10.1142/S021833901750005X.
 35.
Wu Y, Gan Y, Yuan H, Wang Q, Fan Y, Li G, et al. Enriched environment housing enhances the sensitivity of mouse pancreatic cancer to chemotherapeutic agents. Biochem Biophys Res Commun. 2016;473(2):593–9. https://doi.org/10.1016/j.bbrc.2016.03.128.
 36.
Rohatgi A. WebPlotDigitizer: Web based tool to extract data from plots, images, and maps. Austin; 2017. Available online at: https://arohatgi.info/WebPlotDigitizer.
 37.
Peng H, Zhao W, Tan H, Ji Z, Li J, Li K, et al. Prediction of treatment efficacy for prostate cancer using a mathematical model. Sci Rep. 2016;6(1):21599. https://doi.org/10.1038/srep21599.
 38.
Coletti R, Leonardelli L, Parolo S, Marchetti L. A QSP model of prostate cancer immunotherapy to identify effective combination therapies. Sci Rep. 2020;10:1–18.
 39.
Lo WC, Arsenescu RI, Friedman A. Mathematical model of the roles of T cells in inflammatory bowel disease. Bull Math Biol. 2013;75(9):1417–33. https://doi.org/10.1007/s1153801398532.
 40.
Islam Y, Ahmad I, Zubair M, Shahzad K. Double integral sliding mode control of leukemia therapy. Biomed Signal Process Contr. 2020;61:102046. https://doi.org/10.1016/j.bspc.2020.102046.
 41.
Zaid M, Elganainy D, Dogra P, Dai A, Widmann L, Fernandes P, et al. Imagingbased subtypes of pancreatic ductal adenocarcinoma exhibit differential growth and metabolic patterns in the prediagnostic period: implications for early detection. Front Oncol. 2020;10:2629. https://doi.org/10.3389/fonc.2020.596931.
 42.
Srivastava MK, Zhu L, HarrisWhite M, Kar U, Huang M, Johnson MF, et al. Myeloid suppressor cell depletion augments antitumor activity in lung cancer. PLoS One. 2012;7(7):e40677. https://doi.org/10.1371/journal.pone.0040677.
 43.
Liu F, Heiner M, Yang M. Fuzzy stochastic petri nets for modeling biological systems with uncertain kinetic parameters. PLoS One. 2016;11(2):e0149674. https://doi.org/10.1371/journal.pone.0149674.
 44.
Shafiekhani S, Rahbar S, Akbarian F, Jafari AH. Fuzzy Stochastic Petri Net with Uncertain Kinetic Parameters for Modeling TumorImmune System. In: 2018 25th Iranian conference on biomedical engineering and 2018 3rd international Iranian conference on biomedical engineering: ICBME 2018; 2019. https://doi.org/10.1109/ICBME.2018.8703573.
 45.
Liu F, Chen S, Heiner M, Song H. Modeling biological systems with uncertain kinetic data using fuzzy continuous petri nets. BMC Syst Biol. 2018;12(S4):63–74. https://doi.org/10.1186/s1291801805688.
 46.
Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theoretical Biol. 2008:178–96. Academic Press. https://doi.org/10.1016/j.jtbi.2008.04.011.
 47.
Pianosi F, Sarrazin F, Wagener T. A Matlab toolbox for global sensitivity analysis. Environ Model Softw. 2015;70:80–5.
 48.
Siret C, Collignon A, Silvy F, Robert S, Cheyrol T, André P, et al. Deciphering the crosstalk between myeloidderived suppressor cells and regulatory T cells in pancreatic ductal adenocarcinoma. Front Immunol. 2020;10:3070. https://doi.org/10.3389/fimmu.2019.03070.
 49.
TwymanSaint Victor C, Rech AJ, Maity A, Rengan R, Pauken KE, Stelekati E, et al. Radiation and dual checkpoint blockade activate nonredundant immune mechanisms in cancer. Nature. 2015;520 (7547):3737. https://doi.org/10.1038/nature14292.
 50.
Fujimura T, Ring S, Umansky V, Mahnke K, Enk AH. Regulatory T cells stimulate B7H1 expression in myeloidderived suppressor cells in ret melanomas. J Invest Dermatol. 2012;132(4):1239–46. https://doi.org/10.1038/jid.2011.416.
 51.
Schlecker E, Stojanovic A, Eisen C, Quack C, Falk CS, Umansky V, et al. Tumorinfiltrating monocytic myeloidderived suppressor cells mediate CCR5dependent recruitment of regulatory T cells favoring tumor growth. J Immunol. 2012;189(12):5602–11. https://doi.org/10.4049/jimmunol.1201018.
 52.
Serafini P, Mgebroff S, Noonan K, Borrello I. Myeloidderived suppressor cells promote crosstolerance in Bcell lymphoma by expanding regulatory T cells. Cancer Res. 2008;68(13):5439–49. https://doi.org/10.1158/00085472.CAN076621.
 53.
Centuori SM, Trad M, LaCasse CJ, Alizadeh D, Larmonier CB, Hanke NT, et al. Myeloidderived suppressor cells from tumorbearing mice impair TGFβinduced differentiation of CD4+ CD25+ FoxP3+ Tregs from CD4+ CD25− FoxP3− T cells. J Leukoc Biol. 2012;92(5):987–97. https://doi.org/10.1189/jlb.0911465.
 54.
Holmgaard RB, Zamarin D, Li Y, Gasmi B, Munn DH, Allison JP, et al. Tumorexpressed IDO recruits and activates MDSCs in a Tregdependent manner. Cell Rep. 2015;13(2):412–24. https://doi.org/10.1016/j.celrep.2015.08.077.
 55.
Re GL, Re FL, Doretto P, Del Conte A, Amadio M, Cozzi C, et al. Cyclophosphamide with or without fluorouracil followed by subcutaneous or intravenous interleukin2 use in solid tumors: A feasibility offlabel experience. Cytokine. 2019;113:50–60.
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.
Funding
No fund.
Author information
Affiliations
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.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
We confirm that no animals were involved in the present study. Since this study was the mathematical modeling using the experimental data of the other published articles, no ethical approval was necessary for mice models and experimental protocols. All methods were carried out in this study are in accordance with relevant guidelines provided by the Institutional Review Board of Tehran University of Medical Sciences. All experimental data used in this study are adapted from other published papers and we did not have any experimental studies.
Consent for publication
All authors have been read the submitted manuscript and agree with its submission.
Competing interests
All authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Shafiekhani, S., Dehghanbanadaki, H., Fatemi, A.S. et al. Prediction of antiCD25 and 5FU treatments efficacy for pancreatic cancer using a mathematical model. BMC Cancer 21, 1226 (2021). https://doi.org/10.1186/s1288502108770z
Received:
Accepted:
Published:
Keywords
 5FU
 AntiCD25
 GUI
 Fuzzy
 ODE