Skip to main content

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

Abstract

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.

Peer Review reports

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 anti-tumor 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. 1.

    Tumor cells have the logistic growth in deprivation of the immune system [28].

  2. 2.

    Tumor cells stimulate the activation of cytotoxic T lymphocytes [19].

  3. 3.

    In the tumor microenvironment, NK cells also stimulate the activation of cytotoxic T lymphocytes [19].

  4. 4.

    NK cells and activated cytotoxic T lymphocytes induce tumor cell regression [19, 29].

  5. 5.

    The activity of NK cells and cytotoxic T lymphocytes declines over time after encountering tumor cells [19, 30].

  6. 6.

    Tumor cells also recruit MDSCs and stimulate their proliferation in the tumor microenvironment [19, 31].

  7. 7.

    MDSCs suppress cytotoxic T lymphocyte activation mediated by NK cells in the tumor microenvironment [19, 26].

  8. 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. 9.

    Low doses of 5-FU chemotherapy inhibit tumor progression through the deactivation of MDSCs [19, 27].

  10. 10.

    Anti-CD25 immunotherapy inhibits tumor progression through the depletion of Treg cells [34].

Fig. 1
figure 1

Conceptual model of tumor-immune system interactions. The arrows depict activation/induction and blocked arrows indicate blocking/inhibiting

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 et al. study in the control group (no treatment) and 5-FU therapy group are as follows. Control group: C57BL/6 mice were subcutaneously implanted with Panc02 cells (6 × 105 per mouse) in the right flank and tumor size is recorded every 5 days up to 25 days beginning 5 days after tumor inoculation. 5-FU 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 5-FU (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 5-FU therapy in pancreatic cancer. In the experiment of Pu et al., C57BL/6 mice were subcutaneously inoculated with Panc02 cells (2 × 106 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 anti-CD25 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 5-FU and anti-CD25 therapies on pancreatic cancer. The parameterized mathematical model using control, 5-FU, and anti-CD25 data sets, reproduced the published in vivo data of both studies to describe tumor-immune cell interactions through 5-FU and anti-CD25 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 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 factor-beta (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. 19) 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 time-dependent; hence, for example, for the sake of brevity, we have written C instead of C(t).

$$ \frac{dC}{dt}={a}_1C\ \log \left(\frac{C_{max}}{C}\right)-{b}_1N\frac{\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}-\frac{\frac{c_1T\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}}{\left(\frac{1+{e}_1R}{1+{d}_1{(antiCD25)}^2{R}^3}\right)\left(1+{f}_1S\right)}-{g}_1C-{h}_1\log \left(1+{k}_1\mathit{\operatorname{sign}}(antiCD25){\int}_0^t antiCD25\left(\tau \right)C\left(\tau \right) d\tau \right) $$
(1)

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 a1 is the growth rate of cancer cells in absence of treatment and Cmax 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 (\( \frac{\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}} \)) to NKs and CTLs, respectively, and by the depth of access l1 [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 NK-mediated and CTL-mediated tumor cell killing are parameters b1 and c1, respectively [19]. The terms \( \left(\frac{1+{e}_1R}{1+{d}_1{(antiCD25)}^2{R}^3}\right) \) and (1 + f1S) in the denominator of the third term describe the inhibitory effects of Tregs and TGF-β on CTL-mediated cancer cell killing. The proportional parameter d1 regulates the weight (intensity) of anti-CD25 monotherapy and parameters e1 and f1 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 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 g1 represents the apoptotic rate of tumor cells by 5-FU and the parameters h1 and k1 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 \( {\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 anti-CD25 and include this term in a logarithmic function as a non-monotonic 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].

$$ \frac{dN}{dt}={a}_2-{b}_2N+\frac{{c^{\bullet}}_2 IN}{d_2+I}+\frac{e_2 FN}{f_2+F}-\frac{g_2N\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}-{h}_2 RN $$
(2)

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 a2 and normal death rate of NKs with an exponential rate b2, 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 c2, \( {d}_2\frac{\tau_1}{\alpha_1} \), e2 and f2 [34]. The fifth term represents the inactivation of NK cells by interacting with accessible tumor cells with maximum inactivation rate g2 [19]. The last term describes Treg-mediated NK cell killing (inactivation) in a granzyme-B-dependent fashion with a constant rate h2 [34].

$$ \frac{dT}{dt}=-{a}_3T+\frac{b_3T{C}^2}{c_{3+}{C}^2}+\frac{d_3N\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}\left(\frac{1-{m}_3}{1+{n}_3{\left(M-{p}_3\right)}^2}+{m}_3\right)-\frac{e_3T\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}+\frac{{f^{\bullet}}_3 IT}{g_3+I}+\frac{h_3 FT}{k_3+F}-{l}_3 RT $$
(3)

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 a3 [19]. The second term describes the tumor recruitment of CTLs that has a Michaelis-Menten form with parameters b3 and c3 [19]. The third term shows the activation of CTLs as a result of interactions of NKs and accessible tumor cells with a constant rate d3 [19]. Also, MDSCs negatively regulate this interaction to prevent CTL activation. The parameter p3 is the normal population of MDSCs when there is no tumor and parameters m3 and n3 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 e3 [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 f3, g3, h3, k3 and with parameters f3, \( {g}_3\frac{\tau_1}{\alpha_1} \)h3, k3 (eq. 15) [34]. The last term represents Treg-mediated CTL inhibition with a constant rate l3 [34].

$$ \frac{dM}{dt}={a}_4-{b}_4M+\frac{c_4C}{d_4+C} $$
(4)

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 a4 of bone marrow produced MDSCs to the spleen [19]. The second term represents the death rate of MDSCs with a constant rate (b4) in normal conditions or during 5-FU treatment [19]. The third term describes the tumor-induced expansion of MDSCs with parameters c4 and d4 [19].

$$ \frac{dH}{dt}={a}_5-{b}_5H+\frac{{c^{\bullet}}_5 IH}{d_5+I}+\frac{e_5 FH}{f_5+F}-{g}_5 RH $$
(5)

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 (a5) of T helper cells in the thymus [34]. The second term describes the exponential degradation rate (b5) 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 c5, d5, e5, f5 and parameters c5, \( {d}_5\frac{\tau_1}{\alpha_1} \), e5 and f5 (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 g5 [34].

$$ \frac{dR}{dt}={a}_6-{b}_6R+{c}_6T+{d}_6H+\frac{{e^{\bullet}}_6 IR}{f_6+I}-{g}_6 NR-{h}_6(antiCD25)R $$
(6)

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 (a6) [34]. The second term shows the exponential death rate of Tregs (b6) [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 (c6T + d6H) that is proportional to CTLs population with a constant rate c6 and T helper population with a constant rate d6 [34]. The fifth term describes the IL-2-mediated growth rate of Tregs that is modeled by a Michaelis-Menten form with parameters e6, f6, and parameters e6 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 g6 and the last term models anti-CD25-mediated Treg depletion with a constant rate h6 [34].

$$ \frac{dI}{dt}={\alpha}_1H-{\tau}_1I $$
(7)

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].

$$ \frac{dF}{dt}={\beta}_1T+{\beta}_2N+{\beta}_3H-{\tau}_2F $$
(8)

Interferon-gamma (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 half-life [34].

$$ \frac{dS}{dt}={\lambda}_1C-{\tau}_3S $$
(9)

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. 1012) for cytokine concentrations to simplify our model.

$$ \frac{dI}{dt}=0;\kern0.75em I=\frac{\alpha_1}{\tau_1}H $$
(10)
$$ \frac{dF}{dt}=0;\kern1em F=\frac{\beta_1}{\tau_2}T+\frac{\beta_2}{\tau_2}N+\frac{\beta_3}{\tau_2}H $$
(11)
$$ \frac{dS}{dt}=0;\kern0.5em S=\frac{\lambda_1}{\tau_3}C $$
(12)

By substituting eqs. 1012 into eqs. 16, the simplified model, as described using eqs. 1318, 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.

$$ \frac{dC}{dt}={a}_1C\ \log \left(\frac{C_{max}}{C}\right)-{b}_1N\frac{\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}-\frac{\frac{c_1T\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}}{\left(\frac{1+{e}_1R}{1+{d}_1{(antiCD25)}^2{R}^3}\right)\left(1+{f}_1\frac{\lambda_1}{\tau_3}C\right)}-{g}_1C-{h}_1\log \left(1+{k}_1\mathit{\operatorname{sign}}(antiCD25){\int}_0^t antiCD25\left(\tau \right)C\left(\tau \right) d\tau \right) $$
(13)
$$ \frac{dN}{dt}={a}_2-{b}_2N+\frac{c_2 HN}{d_2\frac{\tau_1}{\alpha_1}+H}+\frac{e_2\left(\frac{\beta_1}{\tau_2}T+\frac{\beta_2}{\tau_2}N+\frac{\beta_3}{\tau_2}H\right)N}{f_2+\frac{\beta_1}{\tau_2}T+\frac{\beta_2}{\tau_2}N+\frac{\beta_3}{\tau_2}H}-{g}_2N\frac{\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}-{h}_2 RN $$
(14)
$$ \frac{dT}{dt}=-{a}_3T+\frac{b_3T{C}^2}{c_{3+}{C}^2}+\frac{d_3N\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}\left(\frac{1-{m}_3}{1+{n}_3{\left(M-{p}_3\right)}^2}+{m}_3\right)-{e}_3T\frac{\mathrm{C}}{1+\frac{1}{l_1}{C}^{\frac{1}{3}}}+\frac{f_3 HT}{g_3\frac{\tau_1}{\alpha_1}+H}+\frac{h_3\left(\frac{\beta_1}{\tau_2}T+\frac{\beta_2}{\tau_2}N+\frac{\beta_3}{\tau_2}H\right)T}{k_3+\left(\frac{\beta_1}{\tau_2}T+\frac{\beta_2}{\tau_2}N+\frac{\beta_3}{\tau_2}H\right)}-{l}_3 RT $$
(15)
$$ \frac{dM}{dt}={a}_4-{b}_4M+\frac{c_4C}{d_4+C} $$
(16)
$$ \frac{dH}{dt}={a}_5-{b}_5H+\frac{c_5{H}^2}{d_5\frac{\tau_1}{\alpha_1}+H}+\frac{e_5\left(\frac{\beta_1}{\tau_2}T+\frac{\beta_2}{\tau_2}N+\frac{\beta_3}{\tau_2}H\right)H}{f_5+\frac{\beta_1}{\tau_2}T+\frac{\beta_2}{\tau_2}N+\frac{\beta_3}{\tau_2}H}-{g}_5 RH $$
(17)
$$ \frac{dR}{dt}={a}_6-{b}_6R+{c}_6T+{d}_6H+\frac{e_6 HR}{f_6\frac{\tau_1}{\alpha_1}+H}-{g}_6 NR-{h}_6(antiCD25)R $$
(18)

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 (NRMSE), defined as follow:

T1, T2, …, T5 were the tumor volumes measuring at time points of t1, t2, …, t5. 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−6mm3 [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 (g1,b4) 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 (d1, h1, k1, h6) 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.

Table 1 Summary of parameter values
Fig. 2
figure 2

Data fitting. The blue line shows predicted Panc02 tumor volume dynamics in no treatment case and blue stars are records of Panc02 tumor volume on days 5, 10, 15, 20, and 25 in the control group (Panc02 tumor inoculation is carried out on day 0 and normalized root mean square error (NRMSE = 0.1107) is used as a measure of goodness of fit). The purple dashed line shows predicted dynamics of Panc02 tumor volume and purple ‘>‘are data points gathered from experimental data in the 5-FU treatment group (5-FU therapy is carried out on days 1, 2, 3, and 4 after Panc02 tumor injection on day 0, and data record is carried out on days 5, 10, 15, 20 and 25, and NRMSE of 0.1085 is computed to assess model fitting). The black dotted line shows predicted dynamics of Panc02 tumor cells and black ‘circles’ are data points gathered from in vivo experiments in the anti-CD25 treatment group (anti-CD25 therapy is carried out on days 3, 6, 10, 13, 17, and 20 after tumor inoculation on day 0 and data record is carried out on days 7, 14, 21, 28 and 35, and NRMSE of 0.0928 is computed to assess model fitting)

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 × 105, 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 × 104, 3 × 103, 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. 1012. 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.

Fig. 3
figure 3

Dynamics of all species in strategies: control, 5-FU, Anti-CD25, and combination therapy. The blue dashed lines show predicted dynamics of cells/cytokines in the control group, and the red lines show with 5-FU treatment (on days 1, 2, 3, and 4 after tumor inoculation), the yellow lines show with anti-CD25 therapy (on days 3, 6, 10, 13, 17 and 20 after tumor inoculation) and purple lines depict the dynamics of cells/cytokines under combination therapy. In each subplot (except for the second and third figures in the last panel), the y-axis represents the number/concentration of cell population/cytokine, and the x-axis represents the time in days after tumor inoculation. The initial population/concentration of tumor cells, NKs, CTLs, MDSCs, TCD4+, Treg, IL-2, IFN-γ, and TGF-β have all been set to 6 × 105, 105840, 21 × 104, 3 × 103, 564480, 42336, 2.5107 × 10−6, 3.9348 × 10−7, 5.363 × 10−7, respectively. The second and third figures in the fourth panel are parameters and terms related to 5-FU and anti-CD25 therapies, respectively

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:

$$ E\left({\phi}_a=i,{\phi}_5=\mathrm{j},\mathrm{t}\right)=\frac{C\left({\phi}_a=0,\kern0.5em {\phi}_5=0,\mathrm{t}\right)-C\left({\phi}_a=i,\kern0.75em {\phi}_5=j,\kern0.75em \mathrm{t}\right)}{C\left({\phi}_a=0,\kern0.75em {\phi}_5=0,\kern0.75em \mathrm{t}\right)},i,j=0,1 $$

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 tC(ϕ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, shows a strong synergistic effect of the combination of 5-FU and anti-CD25 (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 5-FU and anti-CD25 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.

Fig. 4
figure 4

Comparison treatment efficacies and Analysis of interactions among treatments. 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) and average tumor size (second panel) as outcomes. The blue line represents the inhibition percentage of tumor growth under 5-FU treatment on days 1, 2, 3, and 4 after tumor injection calculated based on the instantaneous tumor size during time and the red line under anti-CD25 treatment on days 3, 6, 10, 13, 17 and 20 after tumor inoculation and yellow line under combination therapy. (Third panel) The Bliss combination index (CI) for 5-FU and Anti-CD25treatments combination using instantaneous tumor size (CI inst) and average tumor size (CI ave) as an outcome. CI < 1 represents the synergistic effect of two treatments, CI = 1 additive, and CI > 1 antagonistic effect. The green dash-line indicates the threshold below which the 5-FU treatment and anti-CD25 treatment have a synergistic effect

Fig. 5
figure 5

Comparison treatment efficacies and Analysis of interactions among treatments. 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 Fig. 5. A, B, C, and D) and average tumor size (second panel of Figs. 5. A, B, C, and D) as outcomes. The blue lines represent the inhibition percentage of tumor growth under 5-FU treatment on a specified time setting in each subplot that is calculated based on the instantaneous tumor size during time and the red lines under anti-CD25 treatment and yellow lines under combination therapy. (Third panel of Fig. 5. A, B, C, and D) The Bliss combination index (CI) for 5-FU and Anti-CD25 treatments combination using instantaneous tumor size (CI inst) and average tumor size (CI ave) as outcome shows the synergistic interaction of treatments

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 t5F and tCD, 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 5-FU and anti-CD25 will be injected, respectively. GA estimate the optimal timings of 5-FU (t5F) and anti-CD25 (tCD) by minimizing the following cost function to minimize the population of cancer cells in treatment groups:

$$ cost=\frac{1}{3N}\sum \limits_{i=1}^N{\left(\frac{C_{5F}\left({t}_i\right)}{C_{ctrl}\left({t}_i\right)}\right)}^2+{\left(\frac{C_{CD}\left({t}_i\right)}{C_{ctrl}\left({t}_i\right)}\right)}^2+{\left(\frac{C_{comb}\left({t}_i\right)}{C_{ctrl}\left({t}_i\right)}\right)}^2 $$
$$ \left[{C_{5F}}^{\ast}\left({{\boldsymbol{t}}_{\mathbf{5}\boldsymbol{F}}}^{\ast}\right),{C_{CD}}^{\ast}\left({{\boldsymbol{t}}_{\boldsymbol{CD}}}^{\ast}\right),{C_{comb}}^{\ast}\left({{\boldsymbol{t}}_{\mathbf{5}\boldsymbol{F}}}^{\ast },{{\boldsymbol{t}}_{\boldsymbol{CD}}}^{\ast}\right)\kern0.5em \right]=\underset{{{\boldsymbol{t}}_{\mathbf{5}\boldsymbol{F}}}^{\ast },{{\boldsymbol{t}}_{\boldsymbol{CD}}}^{\ast }}{\min }(cost) $$

which Cctrl(ti), C5F(ti), CCD(ti), and Ccomb(ti) represent the population of cancer cells at time point ti, i = 1, 2, …, N in control group, 5-FU group, anti-CD25 group, and combination therapy group, respectively. Also, tN is the end time of simulation which is day 100. The quantities C5F(t5F), CCD(tCD), Ccomb(t5F, tCD) 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.

Fig. 6
figure 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 x-axis 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. A, C, and E) and average tumor size (second panel of figs. 6. A, C and E) as outcomes. The blue line represents the inhibition percentage of tumor growth under 5-FU treatment on specified time points after tumor injection calculated based on the instantaneous tumor size during time and the red line under anti-CD25 treatment on specified time points after tumor inoculation and the yellow line under combination therapy. (Third panel of 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

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.

Fig. 7
figure 7

In silico assessment of treatments in the fuzzy setting. The membership function of the average of dynamics of cancer cells (panel A of figs. 7. A, B, C, and D), NK cells (panel B of figs. 7. A, B, C, and D), CTLs (panel C of figs. 7. A, B, C, and D), MDSCs (panel D of figs. 7. A, B, C, and D), T helper (panel E of figs. 7. A, B, C, and D), Treg (panel F of figs. 7. A, B, C, and D), cytokine IL-2 (panel G of figs. 7. A, B, C, and D), IFN-γ (panel H of figs. 7. A, B, C, and D), and TGF-β (panel K of figs. 7. A, B, C, and D) in the time interval from the start of therapies to day 150 for injection of 5-FU (on days 1, 2, 3, and 4) and anti-CD25 (on days 3, 6, 10, 13, 17 and 20) in fig. 7A, the different timing of 5-FU injection (fig. 7B), and for different timing of anti-CD25 injection (fig. 7C), and different timings of the combination of 5-FU and anti-CD25 (fig. 7D), in the fuzzy setting of kinetic parameter a1 = (0.9, 1, 1.1) × 4.3992×102

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 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 anti-CD25 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 anti-CD25 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 anti-CD25 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 5-FU and anti-CD25 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 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 c, only the significant correlation values (p-value< 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 pro-tumor parameters a1 (tumor growth rate), Cmax (maximum tumor size), and anti-immune parameter b2 (death rate of NK cells) while it is negatively correlated to the anti-tumor parameters b1 (NK-mediated killing rate of tumor cells), l1 (depth of access of immune cell to tumor cells) and a2 (the constant source of NK cells). In this instance, the negative correlation between tumor population at day 20 and parameter l1 represents that if parameter l1 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, a1, Cmax, and b2 are negatively correlated to the NK cell population at day 20 while anti-tumor parameters such as b1 and a2 are positively correlated to the NK cell population at day 20. The NK cell population at day 20 is negatively correlated to the l1 (depth of access of immune cells to the tumor) and g2 (inactivation rate of NKs by tumor cells). Since there is a strong correlation between parameter b3 (Maximum tumor-mediated 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 a1 and Cmax. According to panel A of Fig. 8. A, there exist a negative correlation between b1 (maximum killing rate of tumor cells by NKs), b2 (death rate of NKs), g2 (maximum inactivation rate of NKs by NK-tumor 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.

Fig. 8
figure 8

GSA analysis. Statistically significant PRCC values (p-value< 0.05) for tumor cells, NK cells, CTLs, MDSCs, T helper cells, Tregs, IL-2, IFN-γ, and TGF-β at days 20 (A. A), 50 (A. B) and 100 (A. C) after tumor injection. The mean of PRCC values for five replications of PRCC analysis was depicted in each pixel. Black pixels (‘NaN’) show ‘not a number’ and represent no significant correlation between outcome measures (population/concentration of cells/cytokines, elements in the vertical axis) and kinetic parameters of the model (elements in the horizontal axis). The standard deviation of significant PRCC values (p-value< 0.05) for five replications of PRCC analysis for tumor cells, NK cells, CTLs, MDSCs, T helper cells, Tregs, IL-2, IFN-γ and TGF-β at days 20 (B. A), 50 (B. B) and 100 (B. C) after tumor injection. The standard deviation of significant PRCC values for five replications of PRCC analysis was depicted in each pixel. P-values of PRCC analysis for tumor cells, NK cells, CTLs, MDSCs, T helper cells, Tregs, IL-2, IFN-γ, and TGF-β at days 20 (C. A), 50 (C. B), and 100 (C. C) after tumor injection. The maximum of p-values for five replications of PRCC analysis was depicted in each pixel. (D) The absolute mean value and standard deviation of the elementary effects test. Figures 8. D presents the relative importance of kinetic parameters of the TIS model, considering the population/concentration of cells/ cytokines at day 100 as the read-out, including the population of tumor cells (D. A), NK cells (D. B), CTLs (D. C), MDSCs (D. D), T helper cells (D. E), Tregs (D. F), IL-2 (D. G), IFN-γ (D. 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

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 IL-2, IFN-γ, and TGF-β (in no treatment case) at day 100 of simulation as a read-out.

Our data revealed that the parameters a1Cmax, 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 a1 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 a1, Cmax, b1 (NK-mediated tumor cell killing rate), c1 (CTL-mediated tumor cell killing rate), e1 (Rate of the suppressive effect of Treg on CTL-mediated tumor cell killing) were predicted to play an important role in NK cells dynamics. The parameters a1, Cmax, b1 have both linear and interaction effects while the parameters c1, e1 have interaction effects on the dynamics of NK cells. Also parameter a1 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 a1, Cmax have linear effects while parameter b1 has an interaction effect on the dynamics of MDSCs (panel D of Fig. 8. D). Also, parameters Cmax, b1 have an interaction effect, and parameter a1 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 a1 have interaction and both linear and interaction effects, respectively. As shown in panels F, G, and H of Fig. 8. D, the parameter a1 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 a1, Cmax have linear effects and parameter b1 has an interaction effect on the dynamics of TGF-β.

Discussion

Mathematical modeling of complex networks of tumor-immune 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. Tumor-induced 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 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).

Fig. 9
figure 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

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 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.

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:

Myeloid-derived suppressor cell

5-FU:

5-fluorouracil

CTL:

Cytotoxic T lymphocyte

NK:

Natural killer

IL-2:

Interleukin-2

IFN-γ:

Interferon-γ

TGF-β:

Transforming growth factor-β

TIS:

Tumor-immune system

GUI:

Graphical user interface

ODE:

Ordinary differential equation

NRMSE:

Normalized root-mean 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.

    Article  Google Scholar 

  2. Torre LA, Bray F, Siegel RL, Ferlay J, Lortet-Tieulent J, Jemal A. Global cancer statistics, 2012. CA Cancer J Clin. 2015;65(2):87–108. https://doi.org/10.3322/caac.21262.

    Article  Google Scholar 

  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/0008-5472.CAN-14-0155.

    Article  CAS  PubMed  Google Scholar 

  4. Yang MX, Coates RF, Ambaye A, Gardner J-A, Zubarick R, Gao Y, et al. Investigation of HNF-1B as a diagnostic biomarker for pancreatic ductal adenocarcinoma. Biomark Res. 2018;6(1):25. https://doi.org/10.1186/s40364-018-0139-6.

    Article  PubMed  PubMed Central  Google Scholar 

  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.

    Article  PubMed  Google Scholar 

  6. Corbo V, Tortora G, Scarpa A. Molecular pathology of pancreatic cancer: from bench-to-bedside translation. Curr Drug Targets. 2012;13(6):744–52. https://doi.org/10.2174/138945012800564103.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Guillen-Ponce C, Blazquez J, Gonzalez I. de-Madaria 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/s12094-017-1681-7.

    Article  CAS  PubMed  Google Scholar 

  8. Mohammed S, George Van Buren II, Fisher WE. Pancreatic cancer: advances in treatment. World J Gastroenterol. 2014;20:9354.

    PubMed  PubMed Central  Google Scholar 

  9. Li X, Xu J-X. 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.

    Article  PubMed  Google Scholar 

  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.0529-5815.2017.05.018.

    Article  CAS  PubMed  Google Scholar 

  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.

    CAS  PubMed  Google Scholar 

  12. Herman JM, Swartz MJ, Hsu CC, Winter J, Pawlik TM, Sugar E, et al. Analysis of fluorouracil-based 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.

    Article  PubMed  Google Scholar 

  13. Aroldi F, Zaniboni A. Immunotherapy for pancreatic cancer: present and future. Immunotherapy. 2017;9(7):607–16. https://doi.org/10.2217/imt-2016-0142.

    Article  CAS  PubMed  Google Scholar 

  14. Tang Y, Xu X, Guo S, Zhang C, Tang Y, Tian Y, et al. An increased abundance of tumor-infiltrating 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, Yamazaki-Itoh 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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/1078-0432.CCR-06-0369.

    Article  CAS  PubMed  Google Scholar 

  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/2051-1426-1-10.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Vincent J, Mignot G, Chalmin F, Ladoire S, Bruchard M, Chevriaux A, et al. 5-Fluorouracil selectively kills Tumor-associated Myeloid-derived suppressor cells resulting in enhanced t cell-dependent antitumor immunity. Cancer Res. 2010 [cited 2 Apr 2020]. https://doi.org/10.1158/0008-5472.CAN-09-3690.

  19. Shariatpanahi SP, Shariatpanahi SP, Madjidzadeh K, Hassan M, Abedi-Valugerdi M. Mathematical modeling of tumor-induced immunosuppression by myeloid-derived suppressor cells: implications for therapeutic targeting strategies. J Theor Biol. 2018;442:1–10. https://doi.org/10.1016/j.jtbi.2018.01.006.

    Article  CAS  PubMed  Google Scholar 

  20. Loizides C, Iacovides D, Hadjiandreou MM, Rizki G, Achilleos A, Strati K, et al. Model-based 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Castillo-Montiel E, Chimal-Eguia JC, Tello JI, Piñon-Zaráte G, Herrera-Enríquez M, Castell-Rodrí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/s12976-015-0007-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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.

    Article  PubMed  PubMed Central  Google Scholar 

  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/s12967-018-1673-6.

    Article  CAS  Google Scholar 

  24. Haeno H, Gonen M, Davis MB, Herman JM, Iacobuzio-Donahue CA, Michor F. Computational modeling of pancreatic cancer reveals kinetics of metastasis suggesting optimum treatment strategies. Cell. 2012;148(1-2):362–75. https://doi.org/10.1016/j.cell.2011.11.060.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Condamine T, Gabrilovich DI. Molecular mechanisms regulating myeloid-derived suppressor cell differentiation and function. Trends Immunol. 2011;32(1):19–25. https://doi.org/10.1016/j.it.2010.10.002.

    Article  CAS  PubMed  Google Scholar 

  27. Abedi-Valugerdi M, Wolfsberger J, Pillai PR, Zheng W, Sadeghi B, Zhao Y, Hassan M. Suppressive effects of lowdose 5-fluorouracil, busulfan or treosulfan on the expansion of circulatory neutrophils and myeloid derived immunosuppressor cells in tumor-bearing 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. NK-and CD8+ T cell-mediated eradication of established tumors by peritumoral injection of CpG-containing oligodeoxynucleotides. J Immunol. 2001;167(9):5247–53. https://doi.org/10.4049/jimmunol.167.9.5247.

    Article  CAS  PubMed  Google Scholar 

  30. Kuznetsov VA. Basic models of tumor-immune system interactions identification, analysis and predictions. A survey of models for Tumor-Immune system dynamics. Springer; 1997:237–294, DOI: https://doi.org/10.1007/978-0-8176-8119-7_6.

  31. Youn J, Collazo M, Shalova IN, Biswas SK, Gabrilovich DI. Characterization of the nature of granulocytic myeloid-derived suppressor cells in tumor-bearing mice. J Leukoc Biol. 2012;91(1):167–81. https://doi.org/10.1189/jlb.0311177.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Kyewski B, Suri-Payer 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.

    Article  Google Scholar 

  34. He D-H, Xu J-X. 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.

    Article  CAS  Google Scholar 

  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.

    Article  CAS  PubMed  Google Scholar 

  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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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.

    Article  Google Scholar 

  39. Lo W-C, 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/s11538-013-9853-2.

    Article  PubMed  Google Scholar 

  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.

    Article  Google Scholar 

  41. Zaid M, Elganainy D, Dogra P, Dai A, Widmann L, Fernandes P, et al. Imaging-based subtypes of pancreatic ductal adenocarcinoma exhibit differential growth and metabolic patterns in the pre-diagnostic period: implications for early detection. Front Oncol. 2020;10:2629. https://doi.org/10.3389/fonc.2020.596931.

    Article  Google Scholar 

  42. Srivastava MK, Zhu L, Harris-White 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Shafiekhani S, Rahbar S, Akbarian F, Jafari AH. Fuzzy Stochastic Petri Net with Uncertain Kinetic Parameters for Modeling Tumor-Immune 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/s12918-018-0568-8.

    Article  Google Scholar 

  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 myeloid-derived suppressor cells and regulatory T cells in pancreatic ductal adenocarcinoma. Front Immunol. 2020;10:3070. https://doi.org/10.3389/fimmu.2019.03070.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Twyman-Saint 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):373-7. https://doi.org/10.1038/nature14292.

    Article  CAS  Google Scholar 

  50. Fujimura T, Ring S, Umansky V, Mahnke K, Enk AH. Regulatory T cells stimulate B7-H1 expression in myeloid-derived suppressor cells in ret melanomas. J Invest Dermatol. 2012;132(4):1239–46. https://doi.org/10.1038/jid.2011.416.

    Article  CAS  PubMed  Google Scholar 

  51. Schlecker E, Stojanovic A, Eisen C, Quack C, Falk CS, Umansky V, et al. Tumor-infiltrating monocytic myeloid-derived suppressor cells mediate CCR5-dependent recruitment of regulatory T cells favoring tumor growth. J Immunol. 2012;189(12):5602–11. https://doi.org/10.4049/jimmunol.1201018.

    Article  CAS  PubMed  Google Scholar 

  52. Serafini P, Mgebroff S, Noonan K, Borrello I. Myeloid-derived suppressor cells promote cross-tolerance in B-cell lymphoma by expanding regulatory T cells. Cancer Res. 2008;68(13):5439–49. https://doi.org/10.1158/0008-5472.CAN-07-6621.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Centuori SM, Trad M, LaCasse CJ, Alizadeh D, Larmonier CB, Hanke NT, et al. Myeloid-derived suppressor cells from tumor-bearing 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Holmgaard RB, Zamarin D, Li Y, Gasmi B, Munn DH, Allison JP, et al. Tumor-expressed IDO recruits and activates MDSCs in a Treg-dependent manner. Cell Rep. 2015;13(2):412–24. https://doi.org/10.1016/j.celrep.2015.08.077.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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 interleukin-2 use in solid tumors: A feasibility off-label experience. Cytokine. 2019;113:50–60.

    Article  Google Scholar 

Download references

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

Authors and Affiliations

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.

Corresponding author

Correspondence to Amir Homayoun Jafari.

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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Shafiekhani, S., Dehghanbanadaki, H., Fatemi, A.S. et al. Prediction of anti-CD25 and 5-FU treatments efficacy for pancreatic cancer using a mathematical model. BMC Cancer 21, 1226 (2021). https://doi.org/10.1186/s12885-021-08770-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-021-08770-z

Keywords