In Silico Assessment of 5-FU Therapy via A Mathematical Model with Fuzzy Uncertain Parameters


 Background: Ordinary differential equation (ODE) models widely have been used in mathematical oncology to capture dynamics of tumor and immune cells and evaluate the efficacy of treatments. However, for dynamic models of tumor-immune system (TIS), some parameters are uncertain due to inaccurate, missing or incomplete data, which has hindered the application of ODEs that require accurate parameters. Methods: We extended an available ODE model of TIS interactions via fuzzy logic to illustrate the fuzzification procedure of an ODE model. Fuzzy ODE (FODE) models, in comparison with the stochastic differential equation (SDE) models, assigns a fuzzy number instead of a random number (from a specific probability density function) to the parameters, to capture parametric uncertainty. We used FODE model to predict tumor and immune cells dynamics and assess the efficacy of 5-FU. The present model is configurable for 5-FU chemotherapy injection timing and propose testable hypothesis in vitro/ in vivo experiments. Result: FODE model was used to explore the uncertainty of cells dynamics resulting from parametric uncertainty in presence and absence of 5-FU therapy. 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 NK cells, and decreasing that of myeloid-derived suppressor (MDSC) cells. We investigate the effect of perturbation on model parameters on dynamics of cells through global sensitivity analysis (GSA) and compute correlation between model parameters and cell dynamics. Conclusion: ODE models with fuzzy uncertain kinetic parameters cope with insufficient experimental data in the field of mathematical oncology and can predict cells dynamics uncertainty band. In silico assessment of treatments considering parameter uncertainty and investigating the effect of the drugs on movement of cells dynamics uncertainty band may be more appropriate than in crisp setting.

Cancer related death 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 [1]. 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 [2] [3].
Natural killer cells (NK cells) are cells of the innate immunity which act as first line of defense against cancer cells.
Natural killer 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 (TNF)-related apoptosis-inducing ligand (TRAIL)), production of proinflammatory factors such as Interferon-gamma (IFN-γ) and Nitrite oxide (NO) [4] . 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 [5] [6]. MDSCs abundant in tumor tissues and secondary lymph nodes. MDSCs have inhibitory effects on the immune response to tumor cells through production of inhibitory cytokines such as interleukin 10 (IL-10), transforming growth factor β (TGF-β), by production of reactive oxygen species (ROS), NO, Indoleamine oxidase (IDO), induction of Treg cells and inhibitory effect on anti-tumor function of NK cells [7][8] [9]. Chemotherapy drugs such as Gemcitabine, 5-Fluorouracil (5-FU), and paclitaxel suppress the activity and production of MDSCs, enhancing the protective immune responses against tumors [10][11] [12]. 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 [13] [14] [15].
In addition, due to the complex dynamic responses elicited by the immune system against tumor cells [16] [17], and also inherent noise and fluctuations in signaling pathways and regulatory networks that control different functions of cells, and uncertainty in the kinetic/dynamic rate of tumor-immune system interactions [18], there is a demanding need for new class of computational models to predict dynamics of TIS agents in uncertain environment.
In general, computational models for tumor-immune system interactions can be categorized into two groups including deterministic models and stochastic models [19] [20]. Taking into account all the interactions of the system in the form of pre-defined relationships with exact equations, deterministic models can simulate the dynamics of the model components without regarding uncertainty in behaviors of cells and in their interactions. The ordinary differential equation (ODE) models as deterministic models, widely have been used in systems biology [21] [22] [23][24] [25]. Due to the lack of comprehensive knowledge about how biological processes (from subcellular networks to cell-cell interactions) occurs, there is a structural uncertainty. Also, because of dynamic features of biological networks during time and from patient to patient, error in data acquisition, incomplete or missing data, etc., estimating the exact (crisp value) rate of different behaviors of cells (for example, the rate of apoptosis, expansion, recruitment, etc.) is inconceivable [26][27] [28]. Therefore, deterministic models that can simulate this amount of biological detail will have very complex relationships and many parameters that are constrained by the lack of enough precise experimental data. On the other hand, stochastic models with assigning specific probability density functions (pdf) for different behaviors of cells and by probabilistic rules for simulating cell-cell interactions can simulate the behavioral uncertainty and inherent noise in tumor-immune system [29] [30][31] [32]. Often, stochastic models versus deterministic models require fewer kinetic parameters to predict dynamics, but they are computationally cost. It seems that, by assigning fuzzy uncertain numbers instead of crisp values in kinetic/dynamic rate of cells in deterministic models, and without using stochastic rules or assigning specific pdf for the rate of different behaviors of cells, we can capture dynamics of tumor-immune system in uncertain environment. Therefore, we can use deterministic models with fuzzy uncertain kinetic parameters instead of stochastic models that often are computationally cost. In order to investigate the efficacy of 5-FU treatment in suppressing MDSCs and tumor cells in the inflammatory environment and dynamically analyzing the tumor-immune system interactions, a mathematical model with fuzzy uncertain parameters for the NK cells, MDSCs and CTLs interactions with tumor cells has been developed. By considering some of the kinetic parameters of the model as being fuzzy, current mathematical model studied the effect of uncertainty on kinetic parameters on the dynamics of cancer and immune cells and also analyzed the effect of 5-FU therapy in uncertain environment.
Fuzzy theorem describes 'possibility', different from 'probability' theorem which studies random processes [33]. Fuzzy sets have the capability to 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 the systems biology are based on fuzzy knowledge [34][35] [36]. In a study, a fuzzy inference system (FIS) was used to calculate the interaction rates of the continuous Petri net model [37]. It was shown that this model with the lowest kinetic parameters (often unavailable due to lack of experimental data) and using linguistic rules (qualitative description) about the interacting cells was able to capture the dynamics of the TIS agents and simulate its behaviors [37]. In two recent studies, fuzzy parameters were used to model the uncertainties in the kinetic parameters of stochastic Petri net [38]and continuous Petri net [39].
Current study aimed to use fuzzy uncertain kinetic parameters for the ODE model and create a FODE model. This model was used to simulate TIS behaviors assigning fuzzy and crisp values for kinetic parameters and evaluate the efficacy of 5-FU treatment in both crisp and fuzzy setting.
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-FU therapies for the treatment of lymphoma (El4-luc2 cell line), breast cancer (4T1 cell line) and lung carcinoma (3LL cell line) by a system of ODEs was evaluated [40]. For the same purpose, in another study, combination of radiotherapy and anti-PD-1 therapy by a discreet-time mathematical model was evaluated and temporal dynamics of TIS agents were captured [41]. Also in an another study, the combination of 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 was assessed [42]. Also, 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 [43]. All of these studies explored the efficacy of different treatment strategies regarding crisp values for kinetic parameters, whereas in biological networks including TIS, various sources of uncertainty exist that should be considered in the computational model [44] [45]. Current study aimed to evaluate the efficacy of 5-FU treatment in fuzzy uncertain environment (in different setting of fuzzy parameters) and explore how uncertainty in kinetic parameters of model affects the dynamics of TIS agents in different modalities of 5-FU treatment. For this purpose, present study developed a fuzzy ordinary differential equation (FODE) model to capture uncertain dynamical behavior of TIS agents and investigate how different regimens of 5-FU therapy causes the uncertainty band of tumor cells and immune cells to move.

Methods:
The first part of methods describes the tumor-immune system and explain the detail of ODE formulation of system.
In the next section, we describe how fuzzify the kinetic parameters of an ODE to capture uncertainty band of model's constituents in response to fuzzy uncertainty of kinetic parameters. After that, the results of this study are presented.

Structure of ODE model of TIS:
The mathematical model of tumor-immune system interactions of this study was adapted from the model developed by shariatpanhi et al [40]. The structure of model is based on ordinary differential equation that with deterministic rates simulates the biological and biophysical/biochemical behaviors of TIS agents. TIS of this study consist of tumor cells, NK cells, CTLs and MDSCs. 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 , the term * and term * describes NK-mediated tumor cell killing rate and CTL-mediated tumor killing rate, respectively.

(FODE) Fuzzy Ordinary Differential Equation
A fuzzy set of universal set χ is defined by its membership function: which for an element ∈ χ, determines the value µ ( ) as the membership degree of element in fuzzy set . The value µ ( ) = 0 means that the element is not a member of a fuzzy set and the value µ ( ) = 1 means that the element fully belongs to the fuzzy set . The values 0 < µ ( ) < 1 characterize fuzzy members, which partially belong to the fuzzy set .
In first step, the fuzzy number A of uncertain kinetic parameter is partitioned into α-cuts, α ϵ [0, 1] with k levels, which Therefore, with this simple method, we can apply uncertainty to model kinetic parameters and compute the dynamic of cells to find uncertainty bands of cells/cytokines dynamics.
In the following, we describe the algorithm of fuzzification of kinetic parameters of ODE model.  table 1 and table 2, respectively. All simulations in both crisp and fuzzy settings were run in MATLAB 2019a.

Results:
We begin with evaluation of model in crisp setting (Table 1)  Therefore, we start 5-FU treatment from day 10 after tumor inoculation to inhibit the immune-suppressive effect of MDSCs in inflammatory environment [47].  that of myeloid-derived suppressor cells 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  Figure 4. A, Figure 4. B, Figure 4. C and Figure 4. D shows the uncertainty region of cancer cells, NK cells, CTLs and MDSCs, respectively, with regarding two uncertain kinetic parameters ( , ). The membership function of uncertain parameters and are triangular (as described in Table 2). We assign three α levels

Conclusion:
In this study, we aimed to capture the uncertainty of kinetic parameters of an ODE model of tumor The results of the model simulations with 5-FU injection are qualitatively consistent with that results of in vivo studies [10][48] [49]. Besides, to help understand the mechanisms of tumor-immune system interactions, the model can also provide testable hypotheses in vitro/in vivo experiments. So that the model can be developed via in vitro/in vivo data and parameterized (learned) model can be evaluated in silico environment and model refinement can be done through in vivo/in silico (or in vitro/in silico) iterative process. Regarding the dynamic information obtained from in silico experiments, a more detailed study of the system can be conducted in vitro/in vivo experiments.

Discussion:
One of the major challenges in the field of mathematical oncology is the lack of sufficient precise empirical data (in vitro/ in vivo) for model parameterization and estimating crisp values for parameters. The high sensitivity of differential equation models to kinetic parameters caused the model calibration to be a difficult task. Models such as Petri nets, Agent-based models, Boolean networks, and the cellular automata models have less kinetic parameters compared to the differential equation models. Whereas in differential equation models, many parameters are required to model different behaviors of the system. This matter causes the models based on differential equations to encounter major limitations when they are used to model the systems with insufficient imprecise experimental data. While kinetic parameters' fuzzification in differential equation models eliminates such limitations.
The ODE model of TIS of this study is taken from [40]which is parameterized by imprecise in vivo data sets from Therefore, this study presents the procedure of fuzzification of kinetic parameters of ODE models and although it is illustrated for the ODE model of TIS but is not confined to this system and can be used for any ODE model.
Uncertainty as an inherent feature of biological systems should be considered in computational models. There is two types of uncertainty, randomness and fuzziness. Random uncertainty is simulated using stochastic models such as stochastic Petri net, stochastic differential equations, agent-based models with stochastic rules, probabilistic Boolean networks, Markov model, etc. The stochastic differential equations (SDEs) with random parameters belonging to specific probability density functions (PDFs), create stochastic perturbation terms. By sampling the PDFs of the parameters, and simulating the model with random parameters, the dynamic uncertainty band of the model components (cells) is obtained. In present study, we analyzed the effect of randomness of kinetic parameters through global sensitivity analysis (GSA). In GSA, we performed elementary effect (EE) test [52] and partial rank correlation coefficient (PRCC) test [51] that by Morris sampling and by Latin hypercube sampling strategies, respectively, from uniform distribution of parameters, execute ODE model of TIS (with sampled parameters) and quantify the dynamical behavior of cells. In present study, we wanted to quantify the effect of fuzziness of kinetic parameters, therefore we designed a fuzzy ODE model. The ODE model with fuzzy kinetic parameters creates a framework to include the quantities with imprecise values. In FODE, we decompose the membership function of kinetic parameters into itscuts and discretize each -cut to different levels and execute ODE model with that parameters to compose -cuts of cell's dynamics (to form membership function of cells' dynamics). This fuzzification method of kinetic parameters is similar with that used in stochastic petri net [38][53] and continuous petri net [39].
Using the present model, researchers can test different hypotheses by in silico experiments. Also, they can predict the dynamics of different behaviors/interactions of TIS through the nonlinear complex FODE model. The results of chemotherapy by 5-FU injection can provide a practical tool for the medical community to conduct experiments in the laboratory environment. Finally, we need to point out that even though the fuzzification method for ODE parameters described in this study has been used in a TIS model, but it is not confined to this system and as a powerful tool can be applied to ODE model of any biological network. In silico design of 5-FU treatment conducted by novel FODE model help us to test different schedules of this treatment by virtual clinical therapy and also significantly reduce the cost.

Data accessibility.
MATLAB code and description of GUI is available in the electronic supplementary material.

Ethics.
Since this study was the mathematical modeling using the findings of the other published articles, no ethical approval was necessary.

Competing interests.
We have no competing interests.

Funding.
This work partially was supported by a grant from the Tehran University of Medical Sciences (TUMS)