Predicting Efficacy of 5-Fluorouracil Therapy via a Mathematical Model with Fuzzy Uncertain Parameters.

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 tumor-immune 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 5-fluorouracil (5-FU) chemotherapy. Result: FODE model investigates how parametric uncertainty affects the uncertainty band of cell dynamics in the presence and absence of 5-FU treatment. In silico experiments revealed that the frequent 5-FU 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 myeloid-derived 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. Copyright:

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.[234] There is an inherent noise in biological systems in various scales from subcellular networks to cell-cell 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.[567] In stochastic models, this band of uncertainty is limited to the heterogeneous environment, and behavioral and interactive diversity of system components.[89101112] On the other hand, deterministic models can't produce aleatory uncertainty band of model agents.[131415] Epistemic and aleatory uncertainties can be quantified by sensitivity analysis.[11617] Sensitivity analysis is an essential step in understanding system behavior and interpreting its findings.[1819202122] 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 tumor-immune 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.[242526] 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 L-arginine and 5-fluorouracil (5-FU) 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 anti-PD-1 therapy by a discrete-time mathematical model was evaluated and temporal dynamics of TIS agents were captured.[31] Furthermore, in another study, the combination of a vaccine (GVAX) and anti-PD-1 therapy by a set of partial differential equations was evaluated and spatio-temporal dynamics of TIS constituents in-silico environment were assessed.[32] Furthermore, in another study, the effect of anti-PD-1/PD-L1 therapy and anti-CTLA-4 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 5-FU 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 5-FU 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 5-FU therapy affect the uncertainty band of population of tumor cells/immune cells.

Cancer-related 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.[3536] Natural-killer 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 death-inducing Ligand (Fas), and tumor necrotic factor-related apoptosis-inducing ligand, production of proinflammatory factors such as interferon-gamma and nitrite oxide.[37] Regarding adaptive immunity, cytotoxic T-cells (CTL) are the most competent cells against tumor cells. Regulatory T-cells (Treg) and myeloid-derived suppressor cells (MDSC) are also recruited to the tumor microenvironment for modulating the immune responses in the tumor site.[3839] 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 (IL-10), transforming growth factor β, the production of reactive oxygen species, indoleamine oxidase, induction of Treg cells and inhibitory effect on the anti-tumor function of NK cells.[4041] Chemotherapy drugs such as gemcitabine, 5-FU, and paclitaxel suppress the activity and production of MDSCs, enhancing the protective immune responses against tumors.[4243] Although many studies have pointed to the positive effects of inhibiting MDSCs for tumor treatment, the efficacy of the 5-FU treatment has remained questionable and requires further investigation.[4445] The FODE model of the present study can predict 5-FU 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 tumor-immune 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]

Table 1

Crisp values of the model parameters

Parameter	Unit	Value	Biological description	Reference
a		1.45×10-1	El4-luc2 tumor growth rate in absence of treatment	[30]
Cmax	cell	1×1010	Tumor-carrying capacity	[30]
b		3.23×10-7	NK-mediated tumor killing rate	[30]
η		1.1×10-7	CTL-mediated tumor killing rate	[30]
d		0.7	EL4-luc2 tumor apoptosis rate by low dose 5-FU treatment	[30]
l		100	Depth of access of immune cells to tumor cells	[30]
σ	cell/day	1.4×104	Constant influx rate of NK cells	[30]
f		4.12×10-2	Exponential death rate of NK cells	[30]
g		2.5×10-2	Maximum recruitment rate of NK cells	[30]
h	cell 2	2.02×107	Steepness coefficient of the NK cell recruitment curve	[30]
p		10^-7	Inactivation rate of NK cells after encounter with the tumor	[30]
q	cell	1×1010	Steepness coefficient of the MDSCs production curve	[30]
m		2×10-2	Exponential death rate of CTLs	[30]
j		1×10-1	Maximum recruitment rate of CTLs	[30]
k	cell 2	2.02×107	Steepness coefficient of the CTL recruitment curve	[30]
r		1.1×10-7	Stimulation rate of CTLs as a result of tumor cell and NK cell interactions	[30]
u		1×10-10	Inactivation rate of CTLs after an encounter with the tumor	[30]
v		1×10-2	Differentiation rate of CTL cells to other T-cells	[46]
Smin	none	0.18	Minimum CTL proliferation factor due to inhibition by MDSCs	[30]
γ	Cell−2	6×10-3	Inhibition rate of CTL stimulation by MDSCs	[30]
Mmin	cell	2.5×106	Normal number of MDSCs in C57/BL6 mice	[30]
ρ	cell/day	0.25×Mmin	Normal production rate of MDSCs	[30]
β		0.25	Exponential death rate of MDSCs	[30]
Α		7×106	MDSC expansion rate in EL4-luc2 tumor-bearing mice	[30]

NK - Natural killer; CTL - Cytotoxic T-cells; 5-FU - 5-fluorouracil; MDSCs - Mediated by myeloid-derived suppressor cells

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 NK-mediated tumor cell and CTL-mediated tumor killing rate, respectively. The term represents the therapeutic effect of 5-FU, which reduces the tumor cell population by up to 3 days after treatment (d = 0.7), and then the effect of the 5-FU 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 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.[3046]

Myeloid-derived 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, βM 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 setA 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 level of a is i ∈{1,2…, k}. Also,ai of the fuzzy set A is a crisp subset of X, i.e., Aαi ={x| (μA) (x)≥, ai, x ϵ X, ai ϵ [0, 1]}. Then the values x of crisp subset Aα of the i α-cuts (ai) 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

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, ai,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 = max1,2...,I(MaxUncertaintyoutput)

13: Lower Band = min1,2...,I(MinUncertaintyoutput)

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.

Table 2

Fuzzy values of model parameters

Parameter	Triangular fuzzy membership function parameters
a	(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
a	(0.9,1,1.1)×7×106

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 5-FU treatment by in-silico 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 re-growth, 18–20 days after tumor injection. Therefore, 5-FU treatment from day 10 after tumor inoculation to inhibit the immune-suppressive effect of MDSCs in the inflammatory environment.[47]

Figure 2a and b depict the dynamics of tumor cells, NK cells, CTLs, and MDSCs in the crisp setting and by applying 5-FU chemotherapy on days 10,16, 22 after tumor inoculation on day 0 (under 5-FU treatment on days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76). Figure 2c and d show the inhibition percentage of the instantaneous tumor cell population by applying the 5-FU 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 5-FU chemotherapy injection timings revealed that by 3 injections of 5-FU (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 5-FU 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 5-FU injections, the tumor no longer grows and is almost eliminated while with fewer 5-FU injections the tumor regrowth and metastasize will be probable. To investigate the efficacy of 5-FU treatment for inhibiting the instantaneous tumor cell population, the inhibition percentage of tumor cells affected by different timings of 5-FU treatment (10, 16, 22 in Figure 2c and 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76 in Figure 2d) have been calculated. As it's depicted in Figure 2c, 3-time injection of 5-FU 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 5-FU remains until three days after induction), and after that, the effect of the drug disappeared and tumor regrowth. In according to the Figure 2d by increasing 5-FU 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 5-FU treatment on day 76) and the tumor cells are eliminated.

Figure 2

Dynamics of tumor cells, natural killer (NK), cytotoxic T cells (CTL), and myeloid-derived suppressor cell (MDSC) in response to different timing of 5-fluorouracil (5-FU) treatment over time along with efficacy assessment of 5-FU 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 5-FU treatment on days 10, 16, and 22 after tumor injection on day 0 and (b) by 5-FU 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 5-FU on days 10, 16, 22 and (d) by 5-FU treatment on days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76

It has been found that frequent 5-FU treatment would inhibit the tumor, so for further investigation, the frequency of 5-FU 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 3a] and MDSCs [Figure 3b] and the dynamics of relative population of NK cells [Figure 3c] and CTLs [Figure 3d] in treatment case to no treatment case. As the frequency of 5-FU 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 B16-F10 melanoma-bearing 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 T-cells.[48] In another study, the authors found that 5-FU at low doses can boost circulating NK cells.[49] Also, 5-FU has been used in acute pancreatitis to minimize the abnormal immune cytokines.[50] In this study, in silico experiments revealed that 5-FU can enhance the infiltration of NK cells [as depicted in [Figure 3c] and CTLs [Figure 3d] in the tumor microenvironment and subsequently suppress MDSCs [Figure 3b] and tumor cells [Figure 3a].

Figure 3

Prediction of the effect of the different timing of 5-fluorouracil (5-FU) treatment on TISTIS cells population. The treatment efficacy for different timing of 5-FU was plotted as the percentage of instantaneous tumor growth inhibition (a) and percentage of instantaneous myeloid-derived suppressor cell expansion inhibition (b). The ratio of the population of natural killer cells (c) and cytotoxic T cells (d) in different timing of 5-FU treatment to their population in the control case (no treatment)

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 5-FU treatment. To assess the effect of uncertainty of kinetic parameters on dynamics of cells and for efficacy assessment of different timings of 5-FU 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 5-FU injection. Figures 4a-4d 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 5-FU 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, 5-FU 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 5-FU 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 5-FU 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 three-dimensional simulation plot of cancer cells (a), natural killer cells (b), cytotoxic T cells (c), and myeloid-derived suppressor cells (d) for two different timing of 5-fluorouracil 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

Figure 5

A three-dimensional simulation plot of cancer cells (a), natural killer cells (b), myeloid-derived suppressor cells (c), and cytotoxic T cells (d) for two different timing of 5-FU 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 × 106

Figure 6

The membership function of the average of dynamics of cancer cells (a), Natural killer cells (b), Cytotoxic T cells (c), and myeloid-derived suppressor cells (d) in the time interval from day 10 to day 100 (after first 5-fluorouracil [5-FU] injection) for different timing of 5-FU 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

Figure 7

The membership function of average of dynamics of cancer cells (a), Natural killer cells (b), cytotoxic T cells (c) and myeloid-derived suppressor cells (d) in the time interval from day 10 to day 100 (after first 5-fluorouracil [5-FU] injection) for different timing of 5-FU 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 × 10-2 and α = (09,,1,1.1) × 7 × 106

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 8c (membership function of CTLs with regarding 5-FU injection on days 10, 16, 22), 6d (membership function of MDSCs with regarding 5-FU injection on days 10, 16, 22), Figure 9c (membership function of CTLs with regarding 5-FU injection on days 10, 16, 22, 28, 34, 40, 46, 52, 58, 64, 70, 76) and Figure 9d (membership function of MDSCs with regarding 5-FU 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 5-FU injection) in the control group (no treatment) and 5-FU 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 5-fluorouracil (5-FU) 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). (a-d) Depict the membership function of averaged dynamics of cancer cells, natural killer (NK) cells, cytotoxic T cells (CTLs), and myeloid-derived 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 5-FU treatment (three times) and two fuzzy parameters, the red dotted line shows for the control group and four fuzzy uncertain parameters, and green dash-dotted line depicts for 5-FU treatment (three times) and four uncertain parameters

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 5-fluorouracil [5-FU] injection) in the control group (no treatment) and 5-FU 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 5-FU 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). (a-d) Depict the membership function of averaged dynamics of cancer cells, natural killer (NK) cells, cytotoxic T cells (CTLs), and myeloid-derived suppressor cell (MDSCs), respectively. The blue solid line depicts the membership function of fuzzy uncertain outcome measures (cancer cells, NK cells, CTLs, and myeloid-derived 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 5-FU treatment (10 times) and two fuzzy parameters, the red dotted line shows for the control group and four fuzzy uncertain parameters, and green dash-dotted line depicts for 5-FU treatment (10 times) and four uncertain parameters

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 one-half 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 Figures 10-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 Figures 10-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 Figures 10-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 T-cells, and myeloid-derived 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)

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 T-cells, and myeloid-derived 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)

Statistically significant partial rank correlation coefficient (PRCC) values (P < 0.05) for tumor cells, natural killer cells, cytotoxic T-cells, and myeloid-derived 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)

Statistically significant partial rank correlation coefficient (PRCC) values (P < 0.05) for tumor cells, natural killer cells, cytotoxic T-cells, and myeloid-derived 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)

As illustrated in Figure 10a, 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 myeloid-derived suppressor cells (d) from day 0 to day 100 as a read-out. 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

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 read-out. 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 cell-mediated tumor cell killing rate, respectively, have been identified as being important for the tumor cell output [Figure 14a]. 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 ƞ (CTL-mediated 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 14b]. The parameter a is the most leading parameter on mean population of CTLs and has the most interaction and linear effects [Figure 14c]. 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 14d].

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,[2851] stochastic differential equations (SDEs),[5253] agent-based models with stochastic rules,[811125455] probabilistic Boolean networks,[5657] Markov model,[5859] etc. The SDEs with random parameters belonging to specific probability density functions (PDFs) create stochastic perturbation terms,[6061] 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 a-cuts and discretize each a-cut 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[2862] 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.

Conclusion

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 5-FU chemotherapy which leads to suppression of MDSCs and tumor cells. The FODE model of the present study through mechanistically modeling the different immune-tumor cell interactions predicts the uncertainty band of the cell populations. The model simulates the effect of 5-FU treatment for the improvement of the immune system performance in the inflammatory tumor microenvironment. The FODE is configurable for 5-FU chemotherapy injection timing and was used to investigate the efficacy of 5-FU chemotherapy in the fuzzy setting. Simulation results revealed that 5-FU 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 5-FU therapy limits tumor growth and induces anti-tumor immunity. Since fuzzy analysis takes advantage of parametric uncertainty of TIS, in silico assessment of 5-FU therapy, robust suggestions for the protocol of 5-FU treatment can be generated. The results of chemotherapy by 5-FU injection can provide a practical tool for the medical community to conduct experiments for validation in the laboratory environment. In silico design of 5-FU 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 5-FU injection are qualitatively consistent with the results of in vivo studies[424849] 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.

Financial support and sponsorship

None.

Conflicts of interest

There are no conflicts of interest.