

ORIGINAL ARTICLE 

Year : 2022  Volume
: 12
 Issue : 3  Page : 202218 

Predicting efficacy of 5fluorouracil therapy via a mathematical model with fuzzy uncertain parameters
Sajad Shafiekhani^{1}, Amir Homayoun Jafari^{2}, Leila Jafarzadeh^{3}, Vahid Sadeghi^{4}, Nematollah Gheibi^{5}
^{1} Department of Biomedical Engineering, Research Center for Biomedical Technologies and Robotics, School of Medicine; Students' Scientific Research Center, Tehran University of Medical Sciences, Tehran, Iran ^{2} Departments of Biomedical Engineering, Research Center for Biomedical Technologies and Robotics, School of Medicine, Tehran University of Medical Sciences, Tehran, Iran ^{3} Department of Medical Immunology, School of Medicine, Tehran University of Medical Sciences, Tehran, Iran ^{4} Department of Biomedical Engineering, School of Medicine, Isfahan University of Medical Sciences, Isfahan, Iran ^{5} Department of Medical Biotechnology, School of Paramedical Sciences; Cellular and Molecular Research Center, Research Institute for Prevention of Noncommunicable Diseases, Qazvin University of Medical Sciences, Qazvin, Iran
Date of Submission  21Feb2021 
Date of Decision  12Nov2021 
Date of Acceptance  21Jan2022 
Date of Web Publication  26Jul2022 
Correspondence Address:
Source of Support: None, Conflict of Interest: None
DOI: 10.4103/jmss.jmss_92_21
Background: Due to imprecise/missing data used for parameterization of ordinary differential equations (ODEs), model parameters are uncertain. Uncertainty of parameters has hindered the application of ODEs that require accurate parameters. Methods: We extended an available ODE model of tumorimmune system interactions via fuzzy logic to illustrate the fuzzification procedure of an ODE model. The fuzzy ODE (FODE) model assigns a fuzzy number to the parameters, to capture parametric uncertainty. We used the FODE model to predict tumor and immune cell dynamics and to assess the efficacy of 5fluorouracil (5FU) chemotherapy. Result: FODE model investigates how parametric uncertainty affects the uncertainty band of cell dynamics in the presence and absence of 5FU treatment. In silico experiments revealed that the frequent 5FU injection created a beneficial tumor microenvironment that exerted detrimental effects on tumor cells by enhancing the infiltration of CD8+ T cells, and natural killer cells, and decreasing that of myeloidderived suppressor cells. The global sensitivity analysis was proved model robustness against random perturbation to parameters. Conclusion: ODE models with fuzzy uncertain kinetic parameters cope with insufficient/imprecise experimental data in the field of mathematical oncology and can predict cell dynamics uncertainty band.
Keywords: 5fluorouracil, fuzzy, ordinary differential equation, uncertain
How to cite this article: Shafiekhani S, Jafari AH, Jafarzadeh L, Sadeghi V, Gheibi N. Predicting efficacy of 5fluorouracil therapy via a mathematical model with fuzzy uncertain parameters. J Med Signals Sens 2022;12:20218 
How to cite this URL: Shafiekhani S, Jafari AH, Jafarzadeh L, Sadeghi V, Gheibi N. Predicting efficacy of 5fluorouracil therapy via a mathematical model with fuzzy uncertain parameters. J Med Signals Sens [serial online] 2022 [cited 2022 Dec 10];12:20218. Available from: https://www.jmssjournal.net/text.asp?2022/12/3/202/351889 
Introduction   
Epistemic and aleatory uncertainties are the two sources of uncertainty in computational models.^{[1]} Error in parameter estimation due to the lack of enough and accurate experimental data is the origin of epistemic uncertainty in dynamic models. The second source of uncertainty involves the uncertainty caused by random behaviors in the agents of stochastic/probabilistic models. The aleatory uncertainty can be considered as an intrinsic noise in systems biology.^{[2],[3],[4]} There is an inherent noise in biological systems in various scales from subcellular networks to cellcell interaction scale. In stochastic models, the behaviors and interactions of system components occur based on probabilities with specific functions, so with each execution of the stochastic model, we will find different but bounded dynamics.^{[5],[6],[7]} In stochastic models, this band of uncertainty is limited to the heterogeneous environment, and behavioral and interactive diversity of system components.^{[8],[9],[10],[11],[12]} On the other hand, deterministic models can't produce aleatory uncertainty band of model agents.^{[13],[14],[15]} Epistemic and aleatory uncertainties can be quantified by sensitivity analysis.^{[1],[16],[17]} Sensitivity analysis is an essential step in understanding system behavior and interpreting its findings.^{[18],[19],[20],[21],[22]} Stochastic models usually are computationally cost and they have not been developed as well as deterministic models in the field of immune/oncology mathematical modeling. Developing a method to capture epistemic uncertainty in deterministic models is a good step forward finding reliable outcomes in immune/oncological systems. In this study, we combined fuzzy theorem with an ordinary differential equation (ODE) model of tumorimmune system (TIS) to illustrate the process of fuzzification of a deterministic model. Fuzzy theorem can capture epistemic uncertainty of mathematical models. The fuzzy ODE (FODE) model can predict outcomes of treatments for tumor eradication by considering the effect of epistemic uncertainty. Quantifying the effect of epistemic uncertainty on interactions of agents help us to predict the reliable outcomes of treatments.
The fuzzy theorem describes “possibility,” different from the “probability” theorem which studies random processes.^{[23]} Fuzzy sets can deal with uncertain information. Fuzzy sets describe uncertainty caused by ambiguity, lack of knowledge, incomplete or missing data, imprecision, and errors of measurements. Since fuzzy uncertainty is an inherent feature of biological networks, many models in systems biology are based on fuzzy knowledge.^{[24],[25],[26]} In a study, a fuzzy inference system was used to calculate the interaction rates of the continuous Petri net model.^{[27]} It was shown that this model with the lowest kinetic parameters and using linguistic rules (qualitative description) about the interacting cells was able to capture the dynamics of the agents and simulate their behaviors.^{[27]} In two studies, fuzzy parameters were used to model the uncertainties in the kinetic parameters of stochastic Petri net^{[28]} and continuous Petri net.^{[29]} The current study aimed to use fuzzy uncertain kinetic parameters for an ODE model and create a FODE model. This model was used to simulate the behaviors of a TIS in the fuzzy and crisp setting of kinetic parameters.
Mathematical modeling widely has been used to investigate the efficacy of different treatment strategies for various cancers. For instance, in a recent study, the efficacy of Larginine and 5fluorouracil (5FU) therapies for the treatment of cancer was evaluated by a system of ODEs.^{[30]} For the same purpose, in another study, a combination of radiotherapy and antiPD1 therapy by a discretetime mathematical model was evaluated and temporal dynamics of TIS agents were captured.^{[31]} Furthermore, in another study, the combination of a vaccine (GVAX) and antiPD1 therapy by a set of partial differential equations was evaluated and spatiotemporal dynamics of TIS constituents insilico environment were assessed.^{[32]} Furthermore, in another study, the effect of antiPD1/PDL1 therapy and antiCTLA4 using a set of ODEs and pharmacokinetic/pharmacodynamics equations was investigated.^{[33]} All of these studies explored the efficacy of different treatment strategies regarding crisp values for kinetic parameters, whereas there exist various sources of uncertainty in biological networks including epistemic uncertainty. Predicting treatments outcomes in the presence of epistemic uncertainty help us to generate reliable outcomes consistent with what seen in experimental studies or in clinical trials. The current study aimed to evaluate the efficacy of 5FU treatment in a fuzzy uncertain environment (in the different setting of fuzzy parameters) and explore how uncertainty of kinetic parameters affects the dynamics of agents in different time schedules of 5FU treatment. For this purpose, the present study developed FODE model to capture the uncertain dynamical behavior of TIS agents and investigate how different regimens of 5FU therapy affect the uncertainty band of population of tumor cells/immune cells.
Cancerrelated mortality is regarded as one of the leading causes of death around the world. According to global statistics, 27.5 million people will be diagnosed with cancer by 2040.^{[34]} The immune system can organize immune responses and eliminate tumor cells by identifying tumor antigens. However, tumor cells have evolved into different pathways to escape immune surveillance and metastasize to other tissues. The immune system consists of two general types: innate immunity and adaptive immunity.^{[35],[36]} Naturalkiller cells (NK cells) are cells of innate immunity that act as the first line of defense against cancer cells. NK cells prevent tumor growth with various mechanisms such as direct cell destruction, induction of programmed death through the expression of deathinducing Ligand (Fas), and tumor necrotic factorrelated apoptosisinducing ligand, production of proinflammatory factors such as interferongamma and nitrite oxide.^{[37]} Regarding adaptive immunity, cytotoxic Tcells (CTL) are the most competent cells against tumor cells. Regulatory Tcells (Treg) and myeloidderived suppressor cells (MDSC) are also recruited to the tumor microenvironment for modulating the immune responses in the tumor site.^{[38],[39]} MDSCs abundant in tumor tissues and secondary lymph nodes. MDSCs have inhibitory effects on the immune response through the production of inhibitory cytokines such as interleukin 10 (IL10), transforming growth factor β, the production of reactive oxygen species, indoleamine oxidase, induction of Treg cells and inhibitory effect on the antitumor function of NK cells.^{[40],[41]} Chemotherapy drugs such as gemcitabine, 5FU, and paclitaxel suppress the activity and production of MDSCs, enhancing the protective immune responses against tumors.^{[42],[43]} Although many studies have pointed to the positive effects of inhibiting MDSCs for tumor treatment, the efficacy of the 5FU treatment has remained questionable and requires further investigation.^{[44],[45]} The FODE model of the present study can predict 5FU treatment outcomes and generate testable hypothesis for experimental studies.
Methods   
In the first part of the methods the TIS and details of the ODE formulation of the system have been explained. In the following section, the fuzzification of the kinetic parameters of an ODE model for capturing the uncertainty band of the model's constituents in response to the fuzzy uncertainty of kinetic parameters has been described. Finally, the results of this manuscript will be presented.
Structure of the ordinary differential equation model of tumorimmune system
The mathematical model for TIS interplays of this study has been derived from the model developed by Shariatpanahi et al.^{[30]} The model is an ODE with deterministic rates in which the biological behaviors of TIS agents have been simulated. TIS of this study consists of tumor and NK cells, CTLs, and MDSCs. In the upcoming, the model equations will be expressed and the details description of the model parameters is provided in [Table 1], but for more details the readers are referred to the.^{[30]}
Cancer cells (C)
Equation (1) describes the dynamics of tumor cells, which consists of four terms. The term describes tumor cell growth rate in absence of treatment with carrying capacity (Cmax), the term bNC^{*} and ηTC^{*} shows NKmediated tumor cell and CTLmediated tumor killing rate, respectively. The term represents the therapeutic effect of 5FU, which reduces the tumor cell population by up to 3 days after treatment (d = 0.7), and then the effect of the 5FU disappears (d = 0).^{[30]}
Natural killer cells (N)
Equation (2) portrays the dynamics of NK cells, which consists of four terms. σ is the representative of the constant influx rate of NK cells in the tumor microenvironment, fN can be interpreted as the exponential apoptosis rate of NK cells, and pNC^{*} suggest the recruitment rate of NK cells into the tumor microenvironment and inactivation rate of encountered NK cells with tumor cells, respectively.^{[30]}
Cytotoxic T cells (T)
In equation (3) the dynamics of CTLs has been modelled, which constitute of five distinct parts. The exponential apoptosis rate of CTLs has been described by mT, is the expressive term of the recruitment rate of CTLs into the tumor microenvironment, rNC^{*}S explains the activation rate of CTLs as a result of interactions of NK cells and accessible tumor cells that this stimulation rate is inhibited by MDSCs.^{[30]} The uTC^{*} and vT describe the inactivation rate of CTLs after encountering accessible tumor cells and the differentiation rate of CTLs into other phenotypes of T cells such as regulatory T cells (Tregs), respectively.^{[30],[46]}
Myeloidderived suppressor cells (M)
The suppressive effect of MDSCs on the stimulation rate of CTLs has been described in equation.^{[4]} The last formula explains the dynamics of MDSCs in which is the production rate of splenic MDSCs, and are the MDSC's exponential apoptosis rate and MDSC's expansion rate in an inflammatory environment, respectively.^{[30]} [Table 1] contains more details about the kinetic parameters of the mentioned equations and their references.
Fuzzy ordinary differential equation
A fuzzy set of universal set χ is defined by its membership function:
For an element determines the membership degree of element χ in fuzzy set A. The term μ_{A}(x = 0 means that the element x is not a member of a fuzzy set A and in contrast μ_{A}(x = 1 can be interpreted that the element x fully belongs to the fuzzy set A. The values 0<μ_{A}(x)<1characterize fuzzy members, which partially belong to the fuzzy set A.
As shown in [Figure 1], in the first step, the fuzzy number A of an uncertain kinetic parameter is partitioned into αcuts, α ϵ [0, 1] with k levels, which the i^{th} level of a is i ∈{1, 2, …, k}. Also, a_{i }of the fuzzy set A is a crisp subset of X, i.e., Aα_{i} ={x (μ_{A}) (x)≥, a_{i}, x ϵ X, a_{i} ϵ [0, 1]}. Then the values x of crisp subset Aα_{i} of the i^{th} αcuts (a_{i}) has been discretized to J points, therefore a subset has been generated (in this study the number of quantization points for each αcut has been set to J = 5). ODE model has been executed with parameters Minimum and maximum values of population/concentration for each cells/cytokines have been found to capture lower and upper band of the uncertainty of outcome measures in this αcut. By increasing the α value, the uncertainty band in the input parameter of the model declines, and if, α = 1, there is no uncertainty for the kinetic parameter (similar to the ODE model with crisp kinetic parameters). Therefore, with such a simple method, uncertainty to model kinetic parameters can be applied and the dynamic of cells may be calculated to find uncertainty bands of cells/cytokines dynamics.  Figure 1: (Left) Decomposition of a fuzzy uncertain kinetic parameter to its αcuts and (right) composition of a set of αcuts to create a membership function for fuzzy uncertain outcome measure
Click here to view 
In the following, fuzzification algorithm of kinetic parameters for the ODE model has been accurately described.
Input: A FODE model.
Output: An uncertain band of each outcome measure.
1: For each α level, a_{i},i = 1,2..., k do
2: For each fuzzy number from uncertain kinetic parameters, denoted by do
3: Compute αcuts of the fuzzy parameter, represented as
4: Discretize each αcut, and obtain J crisp values for each fuzzy uncertain number.
5: End for
6: For each combination of crisp values for all fuzzy uncertain parameters do
7: Run ODE model for each combination to obtain dynamics
8: Obtain the minimum (MinUncertaintyoutputi= min) and maximum (MaxUncertaintyoutputi= max).
9: End for
10: Compute the uncertainty band (membership function) for each outcome measure
11: For each record of dynamics of model do
12: Upper Band = max_{i=}_{1,2...,I}(MaxUncertaintyoutput^{i})
13: Lower Band = min_{i=}_{1,2...,I}(MinUncertaintyoutput^{i})
14: end for
All parameters of the model are uncertain and fuzzy numbers should be assigned to all of them, but by increasing the number of fuzzy parameters, the computational load of the model simulation increases exponentially. Therefore, for simplicity for each of the equations 1, 2, 3, and 5, one parameter has been chosen and the simulated model has been executed with the assigned fuzzy value to the selected parameter. As shown in [Table 2], the fuzzy numbers with the triangular membership function are set to take into account ten percent more and less (as an uncertainty band) than that estimated in [Table 1]. All simulations in both crisp and fuzzy settings have been executed in MATLAB 2019a.
Results   
The evaluation of the model will be begun in a crisp setting [Table 1] in which no uncertainty exists in kinetic parameters. So as to predict the dynamics of cells and evaluate the efficacy of different regimens (timing) of 5FU treatment by insilico experiments the TIS interactions have been simulated. It can be derived that slow accumulation of immune suppressive cells such as MDSCs and Treg mediates tumor cell regrowth, 18–20 days after tumor injection. Therefore, 5FU treatment from day 10 after tumor inoculation to inhibit the immunesuppressive effect of MDSCs in the inflammatory environment.^{[47]}
[Figure 2]a and [Figure 2]b depict the dynamics of tumor cells, NK cells, CTLs, and MDSCs in the crisp setting and by applying 5FU chemotherapy on days 10, 16, 22 after tumor inoculation on day 0 (under 5FU treatment on days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76). [Figure 2]c and d show the inhibition percentage of the instantaneous tumor cell population by applying the 5FU treatment on days 10, 16, 22 and 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76, respectively. Assessment the model for different 5FU chemotherapy injection timings revealed that by 3 injections of 5FU (induction on days 10, 16, 22 after tumor inoculation on day 0), the tumor volume on day 100 had been 600 times smaller than of its initial value, whereas with ten 5FU injections (in days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76 after tumor inoculation on day 0), 100 days after tumor induction, the tumor volume was reduced to 1% of its initial inoculated value. Therefore, the simulations showed that with increasing 5FU injections, the tumor no longer grows and is almost eliminated while with fewer 5FU injections the tumor regrowth and metastasize will be probable. To investigate the efficacy of 5FU treatment for inhibiting the instantaneous tumor cell population, the inhibition percentage of tumor cells affected by different timings of 5FU treatment (10, 16, 22 in [Figure 2]c and 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76 in [Figure 2]d) have been calculated. As it's depicted in [Figure 2]c, 3time injection of 5FU causes the tumor inhibition percentage on days 25 to reach its maximum level (because the last injection was on day 22 and the effect of 5FU remains until three days after induction), and after that, the effect of the drug disappeared and tumor regrowth. In according to the [Figure 2]d by increasing 5FU injections (10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) the instantaneous tumor inhibition percentage remains close to 98% (even after discontinuation of 5FU treatment on day 76) and the tumor cells are eliminated.  Figure 2: Dynamics of tumor cells, natural killer (NK), cytotoxic T cells (CTL), and myeloidderived suppressor cell (MDSC) in response to different timing of 5fluorouracil (5FU) treatment over time along with efficacy assessment of 5FU treatment regarding crisp values for kinetic parameters of the model. (a) The dynamics of cancer cells, NK cells, CTLs, and MDSCs (normalized to initial population) affected by 5FU treatment on days 10, 16, and 22 after tumor injection on day 0 and (b) by 5FU treatment on days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70 and 76 after tumor injection on day 0. (c) The instantaneous inhibition percentage of tumor cells affected by 5FU on days 10, 16, 22 and (d) by 5FU treatment on days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76
Click here to view 
It has been found that frequent 5FU treatment would inhibit the tumor, so for further investigation, the frequency of 5FU treatment has been increased from three times (days 10, 16, 22) to twelve times (days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) to capture the dynamics of inhibition percentage of tumor cells [Figure 3]a and MDSCs [Figure 3]b and the dynamics of relative population of NK cells [Figure 3]c and CTLs [Figure 3]d in treatment case to no treatment case. As the frequency of 5FU treatment increases, the percentage of tumor cell, MDSC suppression, and population of NK cells and CTLs increase, respectively. This that antitumor treatment can improve the tumor microenvironment to exert detrimental effects on tumor cells, has been investigated in previous experimental studies. For example in a B16F10 melanomabearing mice model, it was shown that the combination of radiotherapy and hyperthermia enhanced the infiltration of CD8+ T cells, NK cells, and CD11c +/MHCII +/CD86+ dendritic cells, and decreased that of MDSCs and regulatory Tcells.^{[48]} In another study, the authors found that 5FU at low doses can boost circulating NK cells.^{[49]} Also, 5FU has been used in acute pancreatitis to minimize the abnormal immune cytokines.^{[50]} In this study, in silico experiments revealed that 5FU can enhance the infiltration of NK cells [as depicted in [Figure 3]c and CTLs [Figure 3]d in the tumor microenvironment and subsequently suppress MDSCs [Figure 3]b and tumor cells [Figure 3]a.  Figure 3: Prediction of the effect of the different timing of 5fluorouracil (5FU) treatment on TISTIS cells population. The treatment efficacy for different timing of 5FU was plotted as the percentage of instantaneous tumor growth inhibition (a) and percentage of instantaneous myeloidderived suppressor cell expansion inhibition (b). The ratio of the population of natural killer cells (c) and cytotoxic T cells (d) in different timing of 5FU treatment to their population in the control case (no treatment)
Click here to view 
The TIS model has been simulated in the fuzzy uncertain setting to capture the uncertainty band of tumor cells, NK cells, CTLs, and MDSCs in the presence or absence of 5FU treatment. To assess the effect of uncertainty of kinetic parameters on dynamics of cells and for efficacy assessment of different timings of 5FU treatment in the fuzzy uncertain environment, TIS has been simulated by setting the values of parameters (a,f) and (a,f,m,a) to fuzzy numbers, and three times (10, 16, 22) and twelve times (10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) of 5FU injection. [Figure 4]a, [Figure 4]b, [Figure 4]c, [Figure 4]d portrays the uncertainty region of cancer cells, NK cells, CTLs, and MDSCs, respectively, regarding two uncertain kinetic parameters (a,f). The membership function of uncertain a parameters f and are triangular (as described in [Table 2]). Three α levels have been dedicated for fuzzy numbers (α = 0, α = 0.5 and α = 1) and it has been concluded that there is negative correlation between the α and the uncertainty band of kinetic parameter. Actually, α = 0 corresponds to maximum uncertainty for parameters and α = 1 corresponds to the crisp setting of the model (no uncertainty exists). It has been proven that with increasing the frequency of 5FU treatment from 3 to 12 times, the population of cancer cells and MDSCs will be lessen (their uncertainty band shift left toward the lower population of cells), and the population of NK cells and CTLs increase (their uncertainty band shift right toward the higher population of cells). Also, with increasing the uncertainty level α from 1 to 0, the uncertainty band of all cells increases. So, 5FU efficacy was demonstrated in both crisp and fuzzy settings. [Figure 5] depicts the results of model simulation with four fuzzy uncertain numbers as mentioned in [Table 2]. To explore the effect of different timings of 5FU injection (first timing: 10, 16, 22, and second timing: 10, 16, 22, 28, 34, and third timing: 10, 16, 22, 28, 34, 40, 46, and fourth timing: 10, 16, 22, 28, 34, 40, 46, 52, 58, and fifth timing: 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) on the uncertainty band of the tumor and immune cells, TIS was simulated with two and four uncertain parameters. in accordance with [Figure 6] (two uncertain parameters) and [Figure 7] (four uncertain parameters), by boosting the times of 5FU injection shift the uncertainty band of tumor cells and MDSCs to left and shift the uncertainty band of NK cells and MDSCs to right.  Figure 4: A threedimensional simulation plot of cancer cells (a), natural killer cells (b), cytotoxic T cells (c), and myeloidderived suppressor cells (d) for two different timing of 5fluorouracil treatment (same as those given in Figure 2) in the fuzzy setting of kinetic parameters. The two fuzzy uncertain numbers are given as follow: α = (0.9, 1, 1.1) × 1.45× 10^{1} and f = (0.9, 1, 1.1) × 4.12 × 10^{2}
Click here to view 
 Figure 5: A threedimensional simulation plot of cancer cells (a), natural killer cells (b), myeloidderived suppressor cells (c), and cytotoxic T cells (d) for two different timing of 5FU treatment (same as those given in Figure 2) in the fuzzy setting of kinetic parameters. The four fuzzy uncertain numbers are given as follow: α = (0.9, 1, 1.1) × 1.45 × 10^{1} , f = (0.9, 1, 1.1) × 4.12 ×10^{2} , m = (0.9, 1, 1.1) × 2 × 10^{2} and α = (09, 1, 1.1) × 7 × 10^{6}
Click here to view 
 Figure 6: The membership function of the average of dynamics of cancer cells (a), Natural killer cells (b), Cytotoxic T cells (c), and myeloidderived suppressor cells (d) in the time interval from day 10 to day 100 (after first 5fluorouracil [5FU] injection) for different timing of 5FU injection (first timing: 10, 16, 22, and second timing: 10, 16, 22, 28, 34, and third timing: 10, 16, 22, 28, 34, 40, 46, and fourth timing: 10, 16, 22, 28, 34, 40, 46, 52, 5 8, and fifth timing: 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) in the fuzzy setting of kinetic parameters. The two fuzzy uncertain numbers are given as follow: α = (0.9, 1, 1.1) × 1.45 × 10^{1} and f = (0.9, 1, 1.1) × 4.12 × 10^{2}
Click here to view 
 Figure 7: The membership function of average of dynamics of cancer cells (a), Natural killer cells (b), cytotoxic T cells (c) and myeloidderived suppressor cells (d) in the time interval from day 10 to day 100 (after first 5fluorouracil [5FU] injection) for different timing of 5FU injection (first timing: 10, 16, 22, and second timing: 10, 16, 22, 28, 34, and third timing: 10, 16, 22, 28, 34, 40, 46, and fourth timing: 10, 16, 22, 28, 34, 40, 46, 52, 58, and fifth timing: 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) in fuzzy setting of kinetic parameters. The four fuzzy uncertain numbers are given as follow: α = (0.9,1,1.1) × 1.45 × 10^{1} ,f = (0.9,1,1.1) × 4.12 × 10^{2} ,m = (0.9,1,1.1) × 2 × 102 and α = (09,1,1.1) × 7 × 10^{6}
Click here to view 
In the upcoming survey, the effect of increasing the number of fuzzy uncertain parameters (from two parameters to four) has been explored. It has been expected that the number of fuzzy parameters has direct impact on the cell uncertainty band. The result has been illustrated in [Figure 8]c (membership function of CTLs with regarding 5FU injection on days 10, 16, 22), 6d (membership function of MDSCs with regarding 5FU injection on days 10, 16, 22), [Figure 9]c (membership function of CTLs with regarding 5FU injection on days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) and [Figure 9]d (membership function of MDSCs with regarding 5FU injection on days 10, 16, 22). With the increase of fuzzy uncertain numbers from two to four, the uncertainty band of cancer cells and NK cells do not change (due to insensitivity of cancer cells and NK cells to parameters ), while the uncertainty band of CTLs and MDSCs expand.  Figure 8: The membership function of the average of dynamics of TIS' cells in the time interval from day 10 to day 100 (after first 5FU injection) in the control group (no treatment) and 5FU treatment group (on day 10:6:22 after tumor inoculation) and in two different fuzzy settings (two or four fuzzy uncertain kinetic parameters). (a) The membership function of averaged dynamics of cancer cells in control compared with 5fluorouracil (5FU) treatment group and with regarding two fuzzy parameters (same as those given in Figure 4) or four fuzzy parameters (same as those given in Figure 5). (ad) Depict the membership function of averaged dynamics of cancer cells, natural killer (NK) cells, cytotoxic T cells (CTLs), and myeloidderived suppressor cell (MDSCs), respectively. The blue solid line depicts the membership function of fuzzy uncertain outcome measures (cancer cells, NK cells, CTLs, and MDSCs) in the control group (no treatment) with regarding two fuzzy uncertain kinetic parameters (same as Figure 4), the brown dashed lines shows for 5FU treatment (three times) and two fuzzy parameters, the red dotted line shows for the control group and four fuzzy uncertain parameters, and green dashdotted line depicts for 5FU treatment (three times) and four uncertain parameters
Click here to view 
 Figure 9: The membership function of the average of dynamics of TIS' cells in the time interval from day 10 to day 100 (after first 5fluorouracil [5FU] injection) in the control group (no treatment) and 5FU treatment group (on day 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76 after tumor inoculation) and in two different fuzzy settings (two or four fuzzy uncertain kinetic parameters). (a) The membership function of averaged dynamics of cancer cells in control compared with 5FU treatment group and with regarding two fuzzy parameters (same as those given in Figure 4) or four fuzzy parameters (same as those given in Figure 5). (ad) Depict the membership function of averaged dynamics of cancer cells, natural killer (NK) cells, cytotoxic T cells (CTLs), and myeloidderived suppressor cell (MDSCs), respectively. The blue solid line depicts the membership function of fuzzy uncertain outcome measures (cancer cells, NK cells, CTLs, and myeloidderived suppressor cells) in the control group (no treatment) with regarding two fuzzy uncertain kinetic parameters (same as Figure 4), the brown dashed lines shows for 5FU treatment (10 times) and two fuzzy parameters, the red dotted line shows for the control group and four fuzzy uncertain parameters, and green dashdotted line depicts for 5FU treatment (10 times) and four uncertain parameters
Click here to view 
Global sensitivity analysis
Global sensitivity analysis (GSA) identifies a few most influential kinetic parameters from a model with a large number of parameters, which is critical for optimization and structural design. In this section, GSA has been performed to investigate the impact of variation in the models' kinetic parameters on the dynamics of tumor cells, NK cells, CTLs, and MDSCs. To this end, the partial rank correlation coefficient (PRCC) method has been used to compute the correlation between outcome measures (population of cells) concerning all kinetic parameters of the model that have been listed in [Table 1]. Following sensitivity analysis developed in,^{[17]} uniform distribution has been dedicated for all parameters of the model [in the range of onehalf to twice its value in [Table 1]] and generate 1000 samples from these distributions using Latin hypercube sampling (LHS). Then the model has been evaluated by these samples, subsequently the PRCC values and the corresponding P values (significance level) concerning the dynamic of all cells at days 50, 100, and 200 of the model simulation and cells' average dynamics from day 0 to day 200 have been calculated.
The heatmaps related to the GSA include the mean and standard deviation of the PRCC values (five replication) and the corresponding P values (maximum of P values for five replication). The first panel of [Figure 10], [Figure 11], [Figure 12], [Figure 13] illustrate the 5 replicated average of PRCC values for cells (cancer cells, NK cells, CTLs, and MDSCs) populations' record on days 50, 100, 200 after tumor inoculation and for averaged dynamics of cells in the time interval from day 0 to day 200, respectively. The second and third panel of [Figure 10], [Figure 11], [Figure 12], [Figure 13] depicted the standard deviation of PRCC values for 5 replication and corresponding P values (maximum of P values for 5 replication) for record of cell populations on days 50, 100, 200 after tumor inoculation and for averaged dynamics of cells in the time interval from day 0 to day 200, respectively. Each pixel shows the correlation between the population of cells (on the vertical axis) and kinetic parameters (on the horizontal axis). Correlation values range from −1 to +1 and measure linear trend between two variables (population of cell and kinetic parameter). In [Figure 10], [Figure 11], [Figure 12], [Figure 13], only meaningful correlation values (P < 0.05) have been presented.  Figure 10: Statistically significant partial rank correlation coefficient (PRCC) values (P < 0.05) for tumor cells, natural killer cells, cytotoxic Tcells, and myeloidderived suppressor cells at day 50 after tumor injection. The first (a) and second (b) panels show the mean and standard deviation of PRCC values forfive replications of PRCC analysis and the third (c) panel depicts the maximum of their corresponding P values. Black pixels (NaN) show “not a number” and represent no significant correlation between outcome measures (population of cells, elements in the vertical axis) and kinetic parameters of the model (elements in the horizontal axis)
Click here to view 
 Figure 11: Statistically significant partial rank correlation coefficient (PRCC) values (P < 0.05) for tumor cells, natural killer cells, cytotoxic Tcells, and myeloidderived suppressor cells at day 100 after tumor injection. The first (a) and second (b) panels show the mean and standard deviation of PRCC values forfive replications of PRCC analysis and the third (c) panel depicts the maximum of PRCC corresponding P values. Black pixels (NaN) show “not a number” and represents no significant correlation between outcome measures (population of cells, its element in the vertical axis) and kinetic parameters of the model (element in the horizontal axis)
Click here to view 
 Figure 12: Statistically significant partial rank correlation coefficient (PRCC) values (P < 0.05) for tumor cells, natural killer cells, cytotoxic Tcells, and myeloidderived suppressor cells at day 200 after tumor injection. The first (a) and second (b) panels show the mean and standard deviation of PRCC values forfive replications of PRCC analysis and the third (c) panel depicts the maximum of PRCC corresponding P values. Black pixels (NaN) show “not a number” and represents no significant correlation between outcome measures (population of cells, elements in the vertical axis) and kinetic parameters of the model (element in the horizontal axis)
Click here to view 
 Figure 13: Statistically significant partial rank correlation coefficient (PRCC) values (P < 0.05) for an average of dynamics of tumor cells, natural killer cells, cytotoxic Tcells, and myeloidderived suppressor cells in the time interval from day 0 to day 200 after tumor injection. The first (a) and second (b) panels show the mean and standard deviation of PRCC values forfive replications of PRCC analysis and the third (c) panel depicts the maximum of PRCC corresponding P values. Black pixels (NaN) show “not a number” and represents no significant correlation between outcome measures (population of cells, elements in the vertical axis) and kinetic parameters of the model (element in the horizontal axis)
Click here to view 
As illustrated in [Figure 10]a, there is a strong correlation between the population of cancer cells at day 50 and parameters a (growth rate of cancer cells) and Cmax (carrying capacity of tumor). Also, the population of MDSCs at day 50 and parameters a,Cmax have the same correlation. As it has been shown, there is a negative correlation between the population of NK cells (and also CTLs) and parameters a,Cmax and l (depth of access of immune cells to tumor cells). Also, there is a positive correlation between the population of NK cells at day 50 and parameter σ (constant influx rate of NK cells following an encounter with tumor cells). a negative correlation between the population of NK cells at day 50 and parameter p (inactivation rate of NK cells) and also between the population of CTLs at the same time and parameter m (death rate of CTLs) can be observed. Results of PRCC analysis show that there is an inverse correlation between the population of cancer cells (and also MDSCs) at day 50 and parameter j (maximum recruitment rate of CTLs), while there is a positive correlation between the population of CTLs (and also NK cells) with this parameter. The population of MDSCs at day 50 has a positive correlation with parameters a (MDSC's expansion rate), ρ (MDSC's production rate) and also has an inverse correlation with β(MDSC's apoptosis rate), q (Steepness coefficient of the MDSCs production curve). The results of PRCC analysis revealed that the population of CTLs at day 50 has an inverse correlation with parameter v (differentiation rate of CTLs into other types of T cells such as Treg).
In this part, the elementary effects test have been applied to identify which of the kinetic parameters have nonlinear or linear effects.^{[16]} Morris GSA is used to screen and identify which of the 22 kinetic parameters of the model are most influential and have a significant effect on outcome measures (cell dynamics). For 22 kinetic parameters, elementary effects tests using the Morris sampling strategy were taken into account by setting 6 levels in the sampling grid and 1000 trajectories to compute the mean μ^{*} and standard deviation σ. The identified most influential kinetic parameters concerning the mean μ^{*} and interaction effect σ have been depicted in [Figure 14]. The parameters with large σ values indicate nonlinear and interaction effects, while the parameters with large μ^{*}values represent the linear or additive effects. The dashed line (r is the number of trajectories) which all parameters are below than that, translates into a 95% confidence interval.  Figure 14: The absolute mean value and standard deviation of Morris GSA (elementary effects analysis). Figures present the relative importance of kinetic parameters of the TIS model, considering the mean population of tumor cells (a), mean population of natural killer cells (b), mean population of cytotoxic T cells (c), and mean population of myeloidderived suppressor cells (d) from day 0 to day 100 as a readout. Each kinetic parameter is specified by two Morris indices, σ (vertical axis) and * μ* (horizontal axis), which describe the interaction or nonlinear effects and the significance of the effects, respectively
Click here to view 
Morris analysis has been performed by considering the mean population of tumor cells, NK cells, CTLs, and MDSCs (in no treatment case) from day 0 to day 100 of simulation as a readout. The results of the sensitivity analysis with 10% perturbations are showed in [Figure 14]. Findings revealed that the parameters a, Cmax, and, b reflecting the tumor growth rate, carrying capacity or maximal population of tumor cells, and NK cellmediated tumor cell killing rate, respectively, have been identified as being important for the tumor cell output [Figure 14]a. The parameter has the most interaction effect, while parameter is the most linear affecting factor on the dynamics of tumor cells.
The parameters a, Cmax, b and ƞ (CTLmediated tumor killing rate), have been predicted to play a crucial role in the regulation of the NK cells population and all of them have the most interaction and linear effects [Figure 14]b. The parameter a is the most leading parameter on mean population of CTLs and has the most interaction and linear effects [Figure 14]c. The parameters a, Cmax, b ƞ and l, (depth of access of immune cells to tumor cells) have been identified as the most influential kinetic parameters for the mean of MDSCs from day 0 to day 100 of simulation and all of these parameters have interaction and linear effect for MDSCs [Figure 14]d.
Discussion   
One of the major challenges in the mathematical oncology domain is the lack of precise experimental data (in vitro/in vivo) for model parameterization and estimating crisp values for parameters. The high sensitivity of differential equation models to the kinetic parameters caused the model calibration to be a challenging task. This issue ends up with the models based on differential equations that have been widely used in mathematical oncology to overcome major limitations when they are used to model biological systems with insufficient and inaccurate experimental data. The fuzzification of kinetic parameters of differential equation models and assigning a fuzzy uncertain number instead of crisp values for parameters left behind such limitations. In the present study, fuzzy theorem has been applied to illustrate the fuzzification procedure of parameters for an ODE model of TIS interactions.
Uncertainty is an inherent feature of biological systems that should be considered in computational models. There are two types of uncertainty, randomness, and fuzziness. Random uncertainty is simulated using stochastic models such as stochastic Petri net,^{[28],[51]} stochastic differential equations (SDEs),^{[52],[53]} agentbased models with stochastic rules,^{[8],[11],[12],[54],[55]} probabilistic Boolean networks,^{[56],[57]} Markov model,^{[58],[59]} etc. The SDEs with random parameters belonging to specific probability density functions (PDFs) create stochastic perturbation terms,^{[60],[61]} By sampling the PDFs of the parameters, and simulating the model with random parameters, the dynamic uncertainty band of the model components (agents) will be obtained. In this study, the effect of the randomness of kinetic parameters through GSA has been investigated. In GSA, an elementary effect (EE)^{[16]} and PRCC test^{[17]} was implemented. The Morris and LHS strategies, respectively, was carried out to create a sample size from uniform distributions allotted to model parameters. The ODE model was executed by these samples and relation between parameters and cell dynamics was identified. In the present study, to quantify the effect of fuzziness of kinetic parameters, a FODE model has been designed. The ODE model with fuzzy kinetic parameters creates a framework to include the quantities with imprecise values. In FODE, the membership function of kinetic parameters has been decomposed into its acuts and discretize each acut to different levels and execute the ODE model with that level to compose a cuts of cell dynamics (to form membership function of cells' dynamics). This fuzzification method of kinetic parameters is similar to stochastic Petri net^{[28],[62]} and continuous Petri net.^{[29]}
The ODE model of TIS of this article has been taken from^{[30]} which is parameterized by imprecise in vivo data sets. Final aim is the reconstruction of ODE model with the fuzzy theorem to capture the fuzzy uncertainty of kinetic parameters. Due to error in data gathering, inaccurate, incomplete or missing data, natural variability between patients and variable environmental factors, etc., the kinetic parameters of the TIS model are uncertain, and assigning fuzzy uncertain numbers instead of crisp values, can help to capture uncertainty band of the tumor and immune cells (compose the membership function of cells as a result of the uncertainty of kinetic parameters). This study presents the procedure of fuzzification of the kinetic parameters of an ODE model and illustrates this procedure for an ODE model of the TIS model. Even though the described fuzzification method in this study has been used for a TIS model, but it is not confined to this system and this method as a powerful tool can be applied to any ODE model of any system/network.
In this study, capturing the uncertainty of the kinetic parameters of an ODE model of TIS was the ultimate goal. The ODE model mechanistically simulates the interactions of tumor cells, CTLs, NKs, and MDSCs. CTLs and NK cells are the most prominent components of the adaptive and innate immune system play a vital role in encountering with tumor cells, while MDSCs as immature immune cells suppress the immune responses in the inflammatory tumor microenvironment. To consider the parametric uncertainty of the ODE model, the fuzzy theorem has been applied and fuzzy numbers with triangular membership functions instead of deterministic values have been used. Fuzzy parameters are applied as the input source of uncertainty to the TIS and cause uncertainty in tumor and immune cell dynamics. Hence, the uncertainty source of TIS agents' dynamics is the uncertainty of the ODE model's kinetic parameters. The of uncertainty of kinetic parameters originated from the errors in experimental data acquisition, missing or incomplete data values. Furthermore, some specific features of TIS cause the kinetic/dynamic rate of different actions, behaviors, and interaction of system to be uncertain, including, variability and dynamics of TIS from patient to patient and during time and treatment, dynamical features of TIS including dynamic cell size, cell density, various and unpredictable patterns of extracellular ligands and receptors, evolutionary mutation types during time and treatment, diversity in phenotypic patterns, chaotic and complex patterns of vasculature status and so on. All of the mentioned features make the TIS a complex system which requires very sophisticated mathematical functions along with many kinetic parameters that must be estimated by in vitro/in vivo data. For the sake of brevity and simplicity of the mathematical model, ODE model with fuzzy uncertain kinetic parameters has been proposed to construct a FODE model. FODE model can be used for the dynamical analysis of the TIS interactions and in silico assessment of 5FU chemotherapy which leads to suppression of MDSCs and tumor cells. The FODE model of the present study through mechanistically modeling the different immunetumor cell interactions predicts the uncertainty band of the cell populations. The model simulates the effect of 5FU treatment for the improvement of the immune system performance in the inflammatory tumor microenvironment. The FODE is configurable for 5FU chemotherapy injection timing and was used to investigate the efficacy of 5FU chemotherapy in the fuzzy setting. Simulation results revealed that 5FU therapy caused the uncertainty band of the population of cancer cells and MDSCs to shift to the smaller populations while the uncertainty band of the population of NK cells and CTLs shifted to the larger populations. Our data reveals that increasing/decreasing the uncertainty band of the model's fuzzy parameters increases/decreases the uncertainty region of the dynamics of species. It can be concluded that 5FU therapy limits tumor growth and induces antitumor immunity. Since fuzzy analysis takes advantage of parametric uncertainty of TIS, in silico assessment of 5FU therapy, robust suggestions for the protocol of 5FU treatment can be generated. The results of chemotherapy by 5FU injection can provide a practical tool for the medical community to conduct experiments for validation in the laboratory environment. In silico design of 5FU treatment conducted by the novel FODE model help us to test different schedules of this treatment by virtual experiments and significantly reduce the cost and time of real experiments. The results of the model simulations with 5FU injection are qualitatively consistent with the results of in vivo studies^{[42],[48],[49]} Besides, for better understanding the mechanisms of TIS interactions, the model can also provide testable hypotheses in vitro/in vivo experiments. So that the model can be extended via in vitro/in vivo data and parameterized (learned) model may be evaluated in silico environment and model has the ability to be refined through in vivo/in silico (or in vitro/in silico) iterative process. Taking into consideration the dynamical information obtained from in silico experiments, a more detailed study of the system can be conducted in vitro/in vivo experiments.
Acknowledgments
The authors thank Philip Maini at Oxford University and Roobina Boghozian for their critical comments on the manuscript.
Financial support and sponsorship
None.
Conflicts of interest
There are no conflicts of interest.
References   
1.  Renardy M, Hult C, Evans S, Linderman JJ, Kirschner DE. Global sensitivity analysis of biological multiscale models. Curr Opin Biomed Eng 2019;11:10916. 
2.  Eling N, Morgan MD, Marioni JC. Challenges in measuring and understanding biological noise. Nat Rev Genet 2019;20:53648. 
3.  Tsimring LS. Noise in biology. Rep Prog Phys 2014;77:26601. 
4.  Casanova MP. Noise and synthetic biology: How to deal with stochasticity? Nanoethics 2020;14:11322. 
5.  Moore N, Doty D, Zielstorff M, Kariv I, Moy LY, Gimbel A, et al. A multiplexed microfluidic system for evaluation of dynamics of immunetumor interactions. Lab Chip 2018;18:184458. 
6.  Rihan FA, Hashish A, AlMaskari F, Hussein MS, Ahmed E, Riaz MB, et al. “Dynamics of tumorimmune system with fractionalorder.” Journal of Tumor Research 2, no. 1 (2016): 109115. 
7.  Hara A, Iwasa Y. Coupled dynamics of intestinal microbiome and immune system – A mathematical study. J Theory Biol 2019;464:920. 
8.  Allahverdy A, Moghaddam AK, Rahbar S, Shafiekhani S, Mirzaie HR, Amanpour S, et al. An agentbased model for investigating the effect of myeloidderived suppressor cells and its depletion on tumor immune surveillance. J Med Signals Sens 2019;9:1523. [ PUBMED] [Full text] 
9.  Pennisi M, Pappalardo F, Motta S. Agent Based Modeling of Lung MetastasisImmune System Competition. In: International Conference on Artificial Immune Systems; 2009. p. 13. 
10.  Baldazzi V, Castiglione F, Bernaschi M. An enhanced agent based model of the immune system response. Cell Immunol 2006;244:779. 
11.  Gong C, Milberg O, Wang B, Vicini P, Narwal R, Roskos L, et al. A computational multiscale agentbased model for simulating spatiotemporal tumour immune response to PD1 and PDL1 inhibition. J R Soc Interface 2017;14:20170320. 
12.  Cosgrove J, Butler J, Alden K, Read M, Kumar V, CucurullSanchez L, et al. Agentbased modeling in systems pharmacology. CPT Pharmacometrics Syst Pharmacol 2015;4:61529. 
13.  da Silva JG, de Morais RM, da Silva IC, de Arruda Mancera PF. Mathematical models applied to thyroid cancer. Biophys Rev 2019;11:1839. 
14.  Mahasa KJ, Ouifki R, Eladdadi A, Pillis L. Mathematical model of tumorimmune surveillance. J Theor Biol 2016;404:31230. 
15.  de Pillis LG, Radunskaya AE, Wiseman CL. A validated mathematical model of cellmediated immune response to tumor growth. Cancer Res 2005;65:79508. 
16.  Pianosi, Francesca, Fanny Sarrazin, and Thorsten Wagener. “A Matlab toolbox for global sensitivity analysis.” Environmental Modelling & Software 70 (2015):8085. 
17.  Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol 2008;254:17896. 
18.  Lebedeva G, Sorokin A, Faratian D, Mullen P, Goltsov A, Langdon SP, et al. Modelbased global sensitivity analysis as applied to identification of anticancer drug targets and biomarkers of drug resistance in the ErbB2/3 network. Eur J Pharm Sci 2012;46:24458. 
19.  Cândea D, Halanay A, Rădulescu R, Tălmaci R. Parameter estimation and sensitivity analysis for a mathematical model with time delays of leukemia. AIP Conf Proc 2017;1798:20034. 
20.  Poleszczuk J, Hahnfeldt P, Enderling H. Therapeutic implications from sensitivity analysis of tumor angiogenesis models. PLoS One 2015;10:e0120007. 
21.  Alam M, Deng X, Philipson C, BassaganyaRiera J, Bisset K, Carbo A, et al. Sensitivity analysis of an ENteric immunity SImulator (ENISI)based model of immune responses to Helicobacter pylori infection. PLoS One 2015;10:e0136139. 
22.  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:5939. 
23.  Puri, Madan L, Ralescu DA, Zadeh L. “Fuzzy random variables.” In Readings in fuzzy sets for intelligent systems, Morgan Kaufmann, 1993. pp. 265271. 
24.  Shafiekhani S, Poursheykhani A, Rahbar S, Jafari AH. Simulating ATO mechanism and EGFR signaling with fuzzy logic and petri net. J Biomed Phys Eng 2021;11:32536. 
25.  Liu, Fei, Monika Heiner, and David Gilbert. “Fuzzy Petri nets for modelling of uncertain biological systems.” Briefings in bioinformatics 21, no. 1 (2020): 198210. 
26.  Liu F, Sun W, Heiner M, Gilbert D. Hybrid modelling of biological systems using fuzzy continuous Petri nets. Brief Bioinform 2021;22:43850. 
27.  Park, Inho, Dokyun Na, Doheon Lee, and Kwang H. Lee. “Fuzzy continuous Petri Netbased approach for modeling immune systems.” In Neural Nets, Springer, Berlin, Heidelberg, 2005. pp. 278285. 
28.  Liu F, Heiner M, Yang M. Fuzzy stochastic petri nets for modeling biological systems with uncertain kinetic parameters. PLoS One 2016;11:e0149674. 
29.  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:42. 
30.  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:110. 
31.  Serre R, Benzekry S, Padovani L, Meille C, André N, Ciccolini J, et al. Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy. Cancer Res 2016;76:493140. 
32.  Lai X, Friedman A. Combination therapy of cancer with cancer vaccine and immune checkpoint inhibitors: A mathematical model. PLoS One 2017;12:e0178479. 
33.  Milberg O, Gong C, Jafarnejad M, Bartelink IH, Wang B, Vicini P, et al. A QSP model for predicting clinical responses to monotherapy, combination and sequential therapy following CTLA4, PD1, and PDL1 checkpoint blockade. Sci Rep 2019;9:11286. 
34.  
35.  Abbas, Abul K., Andrew H. Lichtman, and Shiv Pillai. Cellular and molecular immunology Ebook. Elsevier Health Sciences, 2014. 
36.  Gajewski TF, Schreiber H, Fu YX. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol 2013;14:101422. 
37.  Guillerey C, Huntington ND, Smyth MJ. Targeting natural killer cells in cancer immunotherapy. Nat Immunol 2016;17:102536. 
38.  Parker KH, Beury DW, OstrandRosenberg S. Myeloidderived suppressor cells: Critical cells driving immune suppression in the tumor microenvironment. Adv Cancer Res 2015;128:95139. 
39.  Liu C, Workman CJ, Vignali DA. Targeting regulatory T cells in tumors. FEBS J 2016;283:273148. 
40.  Groth C, Hu X, Weber R, Fleming V, Altevogt P, Utikal J, et al. Immunosuppression mediated by myeloidderived suppressor cells (MDSCs) during tumour progression. Br J Cancer 2019;120:1625. 
41.  Munder M, Schneider H, Luckner C, Giese T, Langhans CD, Fuentes JM, et al. Suppression of Tcell functions by human granulocyte arginase. Blood 2006;108:162734. 
42.  AbediValugerdi M, Wolfsberger J, Pillai PR, Zheng W, Sadeghi B, Zhao Y, et al. 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:419. 
43.  Umansky V, Blattner C, Gebhardt C, Utikal J. The role of myeloidderived suppressor cells (MDSC) in cancer progression. Vaccines (Basel) 2016;4:E36. 
44.  Srivastava MK, Zhu L, HarrisWhite M, Kar UK, Huang M, Johnson MF, et al. Myeloid suppressor cell depletion augments antitumor activity in lung cancer. PLoS One 2012;7:e40677. 
45.  Si Y, Merz SF, Jansen P, Wang B, Bruderek K, Altenhoff P, Mattheis S, et al. Multidimensional imaging provides evidence for downregulation of T cell effector function by MDSC in human cancer tissue. Sci Immunol 2019;4:eaaw9159. 
46.  Wilson S, Levy D. A mathematical model of the enhancement of tumor vaccine efficacy by immunotherapy. Bull Math Biol 2012;74:1485500. 
47.  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:3737. 
48.  Werthmöller N, Frey B, Rückert M, Lotter M, Fietkau R, Gaipl US. Combination of ionising radiation with hyperthermia increases the immunogenic potential of B16F10 melanoma cells in vitro and in vivo. Int J Hyperthermia 2016;32:2330. 
49.  Orecchioni S, Talarico G, Labanca V, Calleri A, Mancuso P, Bertolini F. Vinorelbine, cyclophosphamide, and 5FU effects on the circulating and intratumoral landscape of immune cells improve antiPDL1 efficacy in preclinical models of breast cancer and lymphoma. Br J Cancer 2018;118:132936. 
50.  Chen XL, Ciren SZ, Zhang H, Duan LG, Wesley AJ. Effect of 5FU on modulation of disarrangement of immuneassociated cytokines in experimental acute pancreatitis. World J Gastroenterol 2009;15:20327. 
51.  Tüysüz F, Kahraman C. Modeling a flexible manufacturing cell using stochastic Petri nets with fuzzy parameters. Expert Syst Appl 2010;37:391020. 
52.  Chen KC, Wang TY, Tseng HH, Huang CY, Kao CY. A stochastic differential equation model for quantifying transcriptional regulatory network in Saccharomyces cerevisiae. Bioinformatics 2005;21:288390. 
53.  Manninen T, Linne ML, Ruohonen K. Developing Itô stochastic differential equation models for neuronal signal transduction pathways. Comput Biol Chem 2006;30:28091. 
54.  Bogle G, Dunbar PR. Agentbased simulation of Tcell activation and proliferation within a lymph node. Immunol Cell Biol 2010;88:1729. 
55.  Castro C, Flores DL, CervantesVásquez D, VargasViveros E, GutiérrezLópez E, MuñozMuñoz F. An agentbased model of the fission yeast cell cycle. Curr Genet 2019;65:193200. 
56.  Zhang SQ, Ching WK, Ng MK, Akutsu T. Simulation study in Probabilistic Boolean Network models for genetic regulatory networks. Int J Data Min Bioinform 2007;1:21740. 
57.  Trairatphisan P, Wiesinger M, Bahlawane C, Haan S, Sauter T. A Probabilistic Boolean Network approach for the analysis of cancerspecific signalling: A case study of deregulated PDGF signalling in GIST. PLoS One 2016;11:e0156223. 
58.  Shafiekhani S, Kraikivski P, Gheibi N, Ahmadian M, Jafari AH. Dynamical analysis of the fission yeast cell cycle via Markov chain. Curr Genet 2021;67:78597. 
59.  Shafiekhani S, Shafiekhani M, Rahbar S, Jafari AH. Extended robust Boolean network of budding yeast cell cycle. J Med Signals Sens 2020;10:95105. 
60.  Tajmirriahi M, Amini Z, Hamidi A, Zam A, Rabbani H. Modeling of retinal optical coherence tomography based on stochastic differential equations: Application to Denoising. IEEE Trans Med Imaging 2021;40:212941. 
61.  Arnold, Ludwig. Stochastic differential equations. New York 1974. 
62.  Shafiekhani, Sajad, Rahbar S, Akbarian F, Jafari AH. “Fuzzy stochastic petri net with uncertain kinetic parameters for modeling tumorimmune system.” In 2018 25th National and 3rd International Iranian Conference on Biomedical Engineering (ICBME), pp. 15. IEEE, 2018. 
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8], [Figure 9], [Figure 10], [Figure 11], [Figure 12], [Figure 13], [Figure 14]
[Table 1], [Table 2]
